Hits XY position different from readout

[lobis]

REST version : v2.2.11
REST commit : a085569b

I am doing some analysis on restG4 simulations of IAXOD0 which has a squared readout on the XY plane rotated 45º.

Using the restViewReadout macro I can see the correct shape and rotation but when I plot the X and Y hits mean position I get a square but with a different rotation, not 45º. I don’t know where the problem can be, the readout seems to be correct but I have little experience with readouts. Maybe there is some issue in the analysis process? Something to do with the prism definition? I use the observable xMeanInPrism (or similar name).

This is the RMl file for the readout:

<?xml version="1.0" encoding="UTF-8" standalone="no" ?>

<readouts>

<globals>
    <variable name="PITCH" value="0.5" overwrite="false" /> <!-- The pitch value in mm -->
    <variable name="MODULES" value="0" overwrite="false" /> <!-- The number of modules per axis MxM -->
    <variable name="CHANNELS" value="10" overwrite="false" /> <!-- The number of pixels per side PxP -->
</globals>

<!-- Updated definition of IAXOD0 readout for final setup - Origin centered-->

<TRestReadout name="IAXOD0_Readout" title="IAXOD0 readout final setup" >
    <myParameter name="nChannels" value="120" />
    <myParameter name="nModules" value="1" />
    <myParameter name="pitch" value="0.5" />
    <myParameter name="pixelSize" value="0.5/sqrt(2.)" />
    <myParameter name="offset" value="0" />
    <myParameter name="extraside" value="0.125" />

    <parameter name="mappingNodes" value="0" />

    <!-- <readoutModule id="[mod]" plane="0" origin="(0.,-pitch*nChannels/sqrt(2.)-pixelSize/2.)" rotation="45" size="(pitch*nChannels+pitch,pitch*nChannels+pitch/2.)" > -->
    <!-- <readoutModule name="module" origin="(0,0)" size="(pitch*(nChannels-1)+pitch/2+2*extraside,pitch*(nChannels-1)+pitch/2+2*extraside)" tolerance="1.e-4" >-->

    <readoutModule name="module" size="(pitch*(nChannels-1)+pitch/2,pitch*(nChannels-1)+pitch/2)" tolerance="1.e-4" > <!-- size=(59.75,59.75)-->
        
        <!-- Y-strips -->
        <for variable="nCh" from="0" to="nChannels-1" step="1" >
            <readoutChannel id="[nCh]" >
                <for variable="nPix" from="0" to="nChannels-1" step="1" >
                <addPixel id="[nPix]" origin="([nCh]*pitch+pitch/2.,[nPix]*pitch-pitch/2.)" size="(pixelSize,pixelSize)" rotation="45" />
                </for>
            </readoutChannel>
        </for>
	    
        <!-- X-strips -->
        <for variable="nCh" from="0" to="nChannels-1" step="1" >
            <readoutChannel id="nChannels+[nCh]" >
            <for variable="nPix" from="0" to="nChannels-1" step="1" >
                <addPixel id="[nPix]" origin="([nPix]*pitch,[nCh]*pitch)" size="(pixelSize,pixelSize)" rotation="45" />
                </for>
            </readoutChannel>
        </for>

    </readoutModule>

    <readoutPlane position="(0,0,-15)mm" planeVector="(0,0,1)" chargeCollection="1" cathodePosition="(0,0,15)mm">
        <addReadoutModule id="0" name="module"
            origin="(0,-42.25)mm"
            rotation="45"
            decodingFile="IAXOD0.dec"/>

    </readoutPlane>

</TRestReadout>

<!-- Updated definition of IAXOD0 readout for final setup -->

<TRestReadout name="IAXOD0_Readout_notcen" title="IAXOD0 readout final setup" >
    <myParameter name="nChannels" value="120" />
    <myParameter name="nModules" value="1" />
    <myParameter name="pitch" value="0.5" />
    <myParameter name="pixelSize" value="0.5/sqrt(2.)" />
    <myParameter name="offset" value="0" />
    <myParameter name="extraside" value="0.125" />

    <parameter name="mappingNodes" value="0" />

    <!-- <readoutModule id="[mod]" plane="0" origin="(0.,-pitch*nChannels/sqrt(2.)-pixelSize/2.)" rotation="45" size="(pitch*nChannels+pitch,pitch*nChannels+pitch/2.)" > -->
    <!-- <readoutModule name="module" origin="(0,0)" size="(pitch*(nChannels-1)+pitch/2+2*extraside,pitch*(nChannels-1)+pitch/2+2*extraside)" tolerance="1.e-4" >-->

    <readoutModule name="module" size="(pitch*(nChannels-1)+pitch/2,pitch*(nChannels-1)+pitch/2)" tolerance="1.e-4" > <!-- size=(59.75,59.75)-->
        
        <!-- Y-strips -->
        <for variable="nCh" from="0" to="nChannels-1" step="1" >
            <readoutChannel id="[nCh]" >
                <for variable="nPix" from="0" to="nChannels-1" step="1" >
                <addPixel id="[nPix]" origin="([nCh]*pitch+pitch/2.,[nPix]*pitch-pitch/2.)" size="(pixelSize,pixelSize)" rotation="45" />
                </for>
            </readoutChannel>
        </for>
	    
        <!-- X-strips -->
        <for variable="nCh" from="0" to="nChannels-1" step="1" >
            <readoutChannel id="nChannels+[nCh]" >
            <for variable="nPix" from="0" to="nChannels-1" step="1" >
                <addPixel id="[nPix]" origin="([nPix]*pitch,[nCh]*pitch)" size="(pixelSize,pixelSize)" rotation="45" />
                </for>
            </readoutChannel>
        </for>

    </readoutModule>

    <readoutPlane position="(0,0,-15.6)mm" planeVector="(0,0,1)"
        chargeCollection="1" 
        cathodePosition="(0,0,22.5)mm" >

        <addReadoutModule id="0" name="module"
            origin="(0,0)mm"
            rotation="45"
            decodingFile="IAXOD0.dec"/>

    </readoutPlane>

</TRestReadout>

<!-- Updated definition of IAXOD0 readout (with 0123 decoding) -->

<TRestReadout name="IAXOD0_Readout_0123" title="IAXOD0 readout with 0123 decoding" >
    <myParameter name="nChannels" value="120" />
    <myParameter name="nModules" value="1" />
    <myParameter name="pitch" value="0.5" />
    <myParameter name="pixelSize" value="0.5/sqrt(2.)" />
    <myParameter name="offset" value="0" />
    <myParameter name="extraside" value="0.125" />

    <parameter name="mappingNodes" value="0" />

    <!-- <readoutModule id="[mod]" plane="0" origin="(0.,-pitch*nChannels/sqrt(2.)-pixelSize/2.)" rotation="45" size="(pitch*nChannels+pitch,pitch*nChannels+pitch/2.)" > -->
    <!-- <readoutModule name="module" origin="(0,0)" size="(pitch*(nChannels-1)+pitch/2+2*extraside,pitch*(nChannels-1)+pitch/2+2*extraside)" tolerance="1.e-4" >-->

    <readoutModule name="module" origin="(0,0)" size="(pitch*(nChannels-1)+pitch/2,pitch*(nChannels-1)+pitch/2)" tolerance="1.e-4" >
        
        <!-- Y-strips -->
        <for variable="nCh" from="0" to="nChannels-1" step="1" >
            <readoutChannel id="[nCh]" >
                <for variable="nPix" from="0" to="nChannels-1" step="1" >
                <addPixel id="[nPix]" origin="([nCh]*pitch+pitch/2.,[nPix]*pitch-pitch/2.)" size="(pixelSize,pixelSize)" rotation="45" />
                </for>
            </readoutChannel>
        </for>
	    
        <!-- X-strips -->
        <for variable="nCh" from="0" to="nChannels-1" step="1" >
            <readoutChannel id="nChannels+[nCh]" >
            <for variable="nPix" from="0" to="nChannels-1" step="1" >
                <addPixel id="[nPix]" origin="([nPix]*pitch,[nCh]*pitch)" size="(pixelSize,pixelSize)" rotation="45" />
                </for>
            </readoutChannel>
        </for>

    </readoutModule>

    <readoutPlane position="(0,0,-15.6)mm" planeVector="(0,0,1)"
        chargeCollection="1" 
        cathodePosition="(0,0,22.5)mm" >

        <addReadoutModule id="0" name="module"
            origin="(0,0)" 
            rotation="45"
            decodingFile="IAXOD0_0123.dec"/>

    </readoutPlane>

</TRestReadout>

<!-- Updated definition of IAXOD0 readout (with 1032 decoding) -->

<TRestReadout name="IAXOD0_Readout_1032" title="IAXOD0 readout with 1032 decoding" >
    <myParameter name="nChannels" value="120" />
    <myParameter name="nModules" value="1" />
    <myParameter name="pitch" value="0.5" />
    <myParameter name="pixelSize" value="0.5/sqrt(2.)" />
    <myParameter name="offset" value="0" />
    <myParameter name="extraside" value="0.125" />

    <parameter name="mappingNodes" value="0" />

    <!-- <readoutModule id="[mod]" plane="0" origin="(0.,-pitch*nChannels/sqrt(2.)-pixelSize/2.)" rotation="45" size="(pitch*nChannels+pitch,pitch*nChannels+pitch/2.)" > -->
    <!-- <readoutModule name="module" origin="(0,0)" size="(pitch*(nChannels-1)+pitch/2+2*extraside,pitch*(nChannels-1)+pitch/2+2*extraside)" tolerance="1.e-4" >-->

    <readoutModule name="module" origin="(0,0)" size="(pitch*(nChannels-1)+pitch/2,pitch*(nChannels-1)+pitch/2)" tolerance="1.e-4" >
        
        <!-- Y-strips -->
        <for variable="nCh" from="0" to="nChannels-1" step="1" >
            <readoutChannel id="[nCh]" >
                <for variable="nPix" from="0" to="nChannels-1" step="1" >
                <addPixel id="[nPix]" origin="([nCh]*pitch+pitch/2.,[nPix]*pitch-pitch/2.)" size="(pixelSize,pixelSize)" rotation="45" />
                </for>
            </readoutChannel>
        </for>
	    
        <!-- X-strips -->
        <for variable="nCh" from="0" to="nChannels-1" step="1" >
            <readoutChannel id="nChannels+[nCh]" >
            <for variable="nPix" from="0" to="nChannels-1" step="1" >
                <addPixel id="[nPix]" origin="([nPix]*pitch,[nCh]*pitch)" size="(pixelSize,pixelSize)" rotation="45" />
                </for>
            </readoutChannel>
        </for>

    </readoutModule>

    <readoutPlane position="(0,0,-15.6)mm" planeVector="(0,0,1)"
        chargeCollection="1" 
        cathodePosition="(0,0,22.5)mm" >

        <addReadoutModule id="0" name="module"
            origin="(0,0)" 
            rotation="45"
            decodingFile="IAXOD0_1032.dec"/>

    </readoutPlane>

</TRestReadout>

<!-- Old definition of IAXOD0 readout -->

<TRestReadout name="Readout_IAXOD0_v1" title="IAXOD0 readout old" >
    <myParameter name="nChannels" value="120" />
    <myParameter name="nModules" value="1" />
    <myParameter name="pitch" value="0.5" />
    <myParameter name="pixelSize" value="0.5/sqrt(2.)" />
    <myParameter name="offset" value="0" />

    <parameter name="mappingNodes" value="0" />

    <!-- <readoutModule id="[mod]" plane="0" origin="(0.,-pitch*nChannels/sqrt(2.)-pixelSize/2.)" rotation="45" size="(pitch*nChannels+pitch,pitch*nChannels+pitch/2.)" > -->

    <readoutModule name="module" size="(pitch*nChannels+pitch/2,pitch*nChannels+pitch/2.)" tolerance="1.e-4" >

        <!-- Y-strips -->
        <for variable="nCh" from="0" to="nChannels-1" step="1" >
            <readoutChannel id="[nCh]" >
                <for variable="nPix" from="0" to="nChannels-1" step="1" >
                <addPixel id="[nPix]" origin="((1+[nCh])*pitch,[nPix]*pitch)" size="(pixelSize,pixelSize)" rotation="45" />
                </for>
            </readoutChannel>
        </for>
	    
        <!-- X-strips -->
        <for variable="nCh" from="0" to="nChannels-1" step="1" >
            <readoutChannel id="nChannels+[nCh]" >
            <for variable="nPix" from="0" to="nChannels-1" step="1" >
                <addPixel id="[nPix]+1" origin="([nPix]*pitch+pitch/2.,[nCh]*pitch+pitch/2.)" size="(pixelSize,pixelSize)" rotation="45" />
                </for>
            </readoutChannel>
        </for>

    </readoutModule>

    <readoutPlane position="(0,0,-15.6)mm" planeVector="(0,0,1)"
        chargeCollection="1" 
        cathodePosition="(0,0,22.5)mm" >

        <for variable="nModX" from="1" to="nModules" step="1" >
            <for variable="nModY" from="1" to="nModules" step="1" >
            <addReadoutModule id="[nModX]+nModules*([nModY]-1)-1" name="module"
                origin="(0.,-pitch*nChannels/sqrt(2.)-pixelSize/2.)" 
                rotation="45"/>
            </for>
        </for>
    </readoutPlane>

</TRestReadout>

        <!--<for variable="nModX" from="1" to="nModules" step="1" />
            <for variable="nModY" from="1" to="nModules" step="1" />
            <addReadoutModule id="0" name="module"
                size="(pitch*(nChannels-1)+pitch/2+2*extraside,pitch*(nChannels-1)+pitch/2+2*extraside)"
                origin="(-extraside,-extraside)" 
                rotation="0"  >
            </for>
        </for>-->
</readouts>

And this is the readout visualization.

This are my hits position in prism:

This is my hits analysis process, where the angle specified (theta) is also 45º

    <TRestHitsAnalysisProcess name="hitsAna" title="Hits analysis template" verboseLevel="info">
        <observable name="energy" value="ON"/>
        <observable name="energyX" value="ON"/>
        <observable name="energyY" value="ON"/>
        <observable name="balanceXYenergy" value="ON"/>
        <observable name="maxHitEnergy" value="ON"/>
        <observable name="minHitEnergy" value="ON"/>
        <observable name="meanHitEnergy" value="ON"/>
        <observable name="meanHitEnergyBalance" value="ON"/>
        <observable name="xMean" value="ON"/>
        <observable name="yMean" value="ON"/>
        <observable name="zMean" value="ON"/>
        <observable name="nHitsX" value="ON"/>
        <observable name="nHitsY" value="ON"/>
        <observable name="balanceXYnHits" value="ON"/>
        <observable name="nHitsSizeXY" value="ON"/>
        <observable name="xy2Sigma" value="ON"/>
        <observable name="xySigmaBalance" value="ON"/>
        <observable name="z2Sigma" value="ON"/>
        <observable name="xySkew" value="ON"/>
        <observable name="zSkew" value="ON"/>
        <!-- x0 and y0 are common coordinates for prism and cylinder. -->
        <parameter name="fiducial_x0" value="(0,0,-15)mm"/>
        <parameter name="fiducial_x1" value="(0,0,15)mm"/>
        <!-- This parameter allows to enable/disable the cylinder fiducialization. -->
        <parameter name="cylinderFiducialization" value="false"/>
        <!-- <parameter name="fiducial_R" value="49" units="mm" /> -->
        <parameter name="fiducial_R" value="48mm"/>
        <observable name="isInsideCylindricalVolume" value="OFF"/>
        <observable name="nInsideCylindricalVolume" value="OFF"/>
        <observable name="energyInsideCylindricalVolume" value="OFF"/>
        <observable name="distanceToCylinderWall" value="OFF"/>
        <observable name="distanceToCylinderTop" value="OFF"/>
        <observable name="distanceToCylinderBottom" value="OFF"/>
        <observable name="xMeanInCylinder" value="OFF"/>
        <observable name="yMeanInCylinder" value="OFF"/>
        <observable name="zMeanInCylinder" value="OFF"/>
        <!-- This parameter allows to enable/disable the prism fiducialization. -->
        <parameter name="prismFiducialization" value="true"/>
        <parameter name="fiducial_sX" value="60mm"/>
        <parameter name="fiducial_sY" value="60mm"/>
        <parameter name="fiducial_theta" value="45"/>
        <observable name="isInsidePrismVolume" value="ON"/>
        <observable name="nInsidePrismVolume" value="ON"/>
        <observable name="energyInsidePrismVolume" value="ON"/>
        <observable name="distanceToPrismWall" value="ON"/>
        <observable name="distanceToPrismTop" value="ON"/>
        <observable name="distanceToPrismBottom" value="ON"/>
        <observable name="xMeanInPrism" value="ON"/>
        <observable name="yMeanInPrism" value="ON"/>
        <observable name="zMeanInPrism" value="ON"/>
    </TRestHitsAnalysisProcess>

[jgalan]

@lobis

when I plot the X and Y hits mean position I get a square but with a different rotation

after which process those positions are obtained? Are they related to which observable from TRestHitsAnalysisProcess?

It is this connected to GetMeanPositionInPrism? An angle not related to the readout is provided there. Perhaps another angle is provided at some other step?

TRestHits::isHitNInsidePrism theta rotation along the axis?

Hi, having a look to the TRestHits::isHitNInsidePrism it seems to me that the theta rotation is intended to rotate the prism along the axis, and not along z, which is actually the case right now. Since I could place the prism in any direction, i.e. axis=x1-x0, the rotation I imagine should be along the axis? @luzon I imagine this worked up to now because we only used prisms defined on z axis. 134 Bool_t TRestHits::isHitNInsidePrism(Int_t n, TVector3 x0, TVector3 x1, Double_t sizeX, Double_t s…

[lobis]

@jgalan

after which process those positions are obtained? Are they related to which observable from TRestHitsAnalysisProcess ?

I just updated the main post with the process.

[jgalan]

@lobis

<parameter name=“fiducial_theta” value=“45”/>

I believe the problem is that the angle must be given in radians.

However, this mean position in prism is related to a fiducialization inside that prism. You could just plot the mean position without fiducialization. I.e. using xMean instead of xMeanInPrism.

[lobis]

@jgalan

I believe the problem is that the angle must be given in radians.

I will try this. Shouldn’t we modify this process so all angles are in degrees? its confusing if some angles are in degrees and others are in radians.

[jgalan]

@lobis

its confusing if some angles are in degrees and others are in radians.

Yes, probably we should adopt the convention that all values in RML should be in degrees, and then be transformed to radians internally? Because usually is easier for everybody to write in degrees, but all methods are usually written in radians. We should store the member directly in radians then.

[jgalan]

Also, I believe the TRestHitsAnalysisProcess should not become so complex. It would be better to have those observables in a separate TRestHitsFiducialAnalysisProcess. @luzon

Also, the methods defined in TRestHitsEvent

https://sultan.unizar.es/rest/classTRestHitsEvent.html

related to prism, cylinder and distances could probably be directly defined in that new TRestHitsFiducialAnalysisProcess.