Nuclear Medicine Detector Applications

There are three detector applications in GAMOS, two related to Nuclear Medicine, i.e. PET, SPECT and another one that, while also used in Nuclear Medicine, it has also an extensive use in other fields, Compton camera. We describe in this chapter the utilities that are common to all detector applications, while we leave those that are specific to each of the in the corresponding chapter.

Identifying Compton interactions

It is frequent that a gamma that enters a detector suffers one or more Compton interactions before the photoelectric one, and then it may leave several hits in different detector elements. GAMOS offers several algorithms to identify which of the hits corresponds to the first gamma interaction, so that the line to the origin of the positron annihilation can be properly reconstructed. This is specially useful for detectors that have 3D identification capabilities, like solid state or liquid xenon detectors. For those without these capabilities the only algorithm of choice is likely the one that is based on the energy of the hits.

The following algorithms are currently available

  • Det1stHitByEnergy: Classifies the reconstructed hits in order of decreasing energy and selects the first one as the one corresponding to the 1st Compton interaction. A parameter can be used to select instead the second, third, …, highest energy hit: /gamos/setParam Det1stHitByEnergy:Order ORDER
  • Det1stHitByXYPos: If the detector has cylindrical symmetry and the gamma travels preferentially increasing the position in the XY plane (i.e. the cylinder radius) you may use this algorithm, which classifies the reconstructed hits in order of increasing position in the XY plane (i.e. sqrt(x*x+y*y) and selects the first one as the one corresponding to the 1st Compton interaction. A parameter can be used to select instead the second, third, …, highest energy hit: /gamos/setParam Det1stHitByXYPos:Order ORDER
  • Det1stHitByXPos: If the detector is a block and the gamma travels preferentiallly increasing the position along the positive or negative X axis you may use this algorithm, which classifies the reconstructed hits in order of increasing absolute value of X position and selects the first one as the one corresponding to the 1st Compton interaction. A parameter can be used to select instead the second, third, …, highest energy hit: /gamos/setParam Det1stHitByXPos:Order ORDER
  • Det1stHitByXYZPos: Even if the detector has some symmetry, it may happen that the best variable is the 3D position (i.e. the spherical radius). You may then use this algorithm, which classifies the reconstructed hits in order of increasing 3D position (i.e. sqrt(x*x+y*y+z*z) and selects the first one as the one corresponding to the 1st Compton interaction. A parameter can be used to select instead the second, third, …, highest energy hit: /gamos/setParam Det1stHitByXYZPos:Order ORDER
  • Det1stHitByDistanceToOther: This algorithm is a generalization of the two above: it calculate the minimum distance from each reconstructed hit to any of the reconstructed hit that belongs to the other sets (a PET event has two sets, corresponding to the two gammas from positron annihilation) and classifies the reconstructed hits in order of increasing value of this minimum distance. It then the first one as the one corresponding to the 1st Compton interaction. A parameter can be used to select instead the second, third, …, highest energy hit: /gamos/setParam Det1stHitByDistanceToOther:Order ORDER This algorithm cannot be used for SPECT detectors
  • Det1stHitByComptonCone: The idea of this algorithm is to profit from the fixed relationship between the initial and final energy and the deviation angle in Compton interactions. It selects a pair of reconstructed hits and between the two chooses one as corresponding to the first interaction. Assuming that the gamma suffered the first interaction with energy equal to the electron mass; for SPECT detectors you may change this value with the parameter: /gamos/setParam DetRecHitCone:InitialHitEnergy ENERGY It then calculates the deviation angle using the initial energy and final energy (= initial - hit energy). Using the line joining the two reconstructed hits and this angle it can build the cone of possible directions of the gamma before the first interaction. If the hit selected as first is indeed the one corresponding to the first gamma interaction, this cone points ideally to the origin of the gamma; as the origin of the gamma cannot be determined in a real detector, it is assumed that there is another gamma from the positron annihilation with the same direction and opposite sense and that if has left a hit in the opposite side of the detector before any other interaction. As it cannot be known which of the hits in the other side of the detector corresponds to the first interaction, the algorithm computes the distance from the cone surface to each of these reconstructed hits (only those that do not belong to the hit set being analysed are considered); the hit with smallest distance is selected, and this minimal distance is stored. The whole algorithm is repeated selecting as first the other hit of the pair of reconstructed hits, and then making other pairs of hits with the rest of hits in the set. By looking at all the minimal distances stored, it selects as good combination the one with the smallest distance. For PET detectors when the algorithms is applied to the second set of hits, the position of the first set is used, instead of looking at the position of all hits in the other side of the detector.

These algorithms act only on those event that have been accepted by the PET/SPECT/ComptonCamera classifier. To use them the first step is to merge the reconstructed hits that are supposed to correspond to the interaction of the same gamma into a single one. This is done by setting the distance between hits to be merged with the parameter (substitute PET by SPECT or ComptonCamera):

/gamos/setParam PET:EvtClass:ComptonRecHitDist DISTANCE

which by default takes a value of 0. By default the position of the set of hits is the one of the hit with highest energy (i.e. the Det1stHitByEnergy algorithm is used). Other algorithms may be selected instead with the parameter:

/gamos/setParam PET:EvtClass:1stHitAlgorithm ALGORITHM

where ALGORITHM is one of those enumerated above.

For the case of PET and Compton camera events it is possible to use a different algorithm for the first and second set of hits, what may be specially useful if the Compton cone algorithm is used, but also in other cases. To do it two different parameters should be set:

/gamos/setParam PET:EvtClass:1stHitAlgorithmFirst ALGORITHM

/gamos/setParam PET:EvtClass:1stHitAlgorithmSecond ALGORITHM

Compton studies histograms

Several histograms can be filled to help you in determining which is the most efficient algorithm to identify the hit corresponding to the first gamma interaction. To produce them you have to set the parameter (substitute PET by SPECT or CC):

/gamos/setParam PET:EvtClass:ComptonStudyHistos 1

which by default takes a value of 0.

The histograms count the proportion of times that the algorithm used selected correctly the hit closest to the first gamma interaction and how are the hits are from the gamma interactions (what can serve you to have an estimate of the precision you may reach).

The logic of this class is the following: the interactions of the original gammas are stored (we understand by an original gamma the one that is created as primary particle, comes from the annihilation of a positron that is a primary particle or, if a radioactive ion was the primary particle, the gamma originated by the ion decay or by the annihilation of the positron that was created by the ion decay). Each gamma interaction is associated to the closest reconstructed hit. If the first interaction identification algorithm chooses as first hit the one that is associated to the first gamma interaction, it is considered that the choice is correct.

The following histograms are filled:

  • N gamma interactions: Number of original gamma interactions in one event
  • N rec hits: Number of reconstructed hits in one event
  • N gamma interactions - N rec hits: Number of original gamma interactions minus number of reconstructed hits in one event

Histograms of data about the interactions and the reconstructed hits

Another set of histograms may help to evaluate if the energy is a good criteria to identify the reconstructed hit that corresponds to the first gamma interaction. The energy lost at each first interaction of an “original gamma” is plotted, and also the energy of the second one, the third one, …. To avoid reserving space for an unlimited number of interactions, the third and later interactions are grouped in a unique histogram. You may change the interaction number to start this grouping with the parameter:

/gamos/setParam DetCAlgoEnergy:HistoGroupingNumber NUMBER

Also histograms of the difference between the first interaction and each of the other ones, the second and each of the other ones, etc. are done. And two-dimensional histograms of the energy of the first vs the rest, the second vs the rest, etc. The name of the histograms are (assuming the grouping number is three):

  • DetCompton:Algo:Interaction: Energy Energy of gamma interactions
  • DetCompton:Algo:Interaction: Energy: 1st Energy of first gamma interaction
  • DetCompton:Algo:Interaction: Energy: 2nd Energy of second gamma interaction
  • DetCompton:Algo:Interaction: Energy: 3rd+ Energy of third and higher gamma interactions
  • DetCompton:Algo:Interaction: Energy: 1st - others Energy of first gamma interaction minus energy of other interaction
  • DetCompton:Algo:Interaction: Energy: 2nd - others Energy of second gamma interaction minus energy of other interaction
  • DetCompton:Algo:Interaction: Energy: 3rd+ - others Energy of third and higher gamma interactions minus energy of other interaction
  • DetCompton:Algo:Interaction: Energy: 1st - others Energy of first gamma interaction vs energy of other interaction
  • DetCompton:Algo:Interaction: Energy: 2nd - others Energy of second gamma interaction vs energy of other interaction
  • DetCompton:Algo:Interaction: Energy: 3rd+ - others Energy of third and higher gamma interactions vs energy of other interaction

By default all these histograms have 100 bins between 0. and 1.

The same set of histograms can done for the position instead of energy. Several position variables can be selected: X, Y, Z, XY ( sqrt(X*X+Y*Y) ), XZ, YZ, XYZ.

To fill these histograms the following parameter should be used:

/gamos/setParam DetComptonStudyHistosUA:AlgorithmVariables VAR1 VAR2 … VARN

where the variables can be Energy, , XPos, YPos, ZPos, XYPos, XZPos, YZPos or XYZPos.

At the same time that the histograms about gamma interaction data, they are filled also histograms about the reconstructed hits that are associated to each interaction are filled (so the first hit is the one closest to the first interaction, etc.) The only difference is that instead of the word Interaction in their name, they have the word RecHit.

The classes that manage this histograms are plug-ins, so that other variables could easily be created.

Classification of good and bad identification histograms as a function of variable

It may happen that different algorithms to identify the first gamma interaction may be optimal for different types of reconstructed hit set. For example for sets with only two hits the Compton cone algorithm may be optimal, while for sets with more hits it may give worse results; or the minimal position algorithm may be optimal when hits have small energy, but if there is one with high energy (and therefore high deviation angle) it may give worse results. To help you in doing this optimization several histograms about the number of good and bad first gamma interaction identification can be done for different variables, e.g. for different number of reconstructed hits in the set, for different maximum hit energy bins, for different position variables. The histograms are about the distance between the reconstructed hit sets (using the position of the hit identified as corresponding to the first interaction) and the interactions. For example, for the “number of reconstructed hits” variable the following histograms are filled:

  • DetCompton:Classif: NRecHit: INTERVAL_LOWER_EDGE-INTERVAL_UPPER_EDGE: Dist RecHitSet - Associated Interaction Distance between a reconstructed hit set and the interaction to which it is associated
  • DetCompton:Classif: NRecHit: INTERVAL_LOWER_EDGE-INTERVAL_UPPER_EDGE: Dist RecHitSet - 1st Interaction Distance between a reconstructed hit set and the first interaction
  • DetCompton:Classif: NRecHit: INTERVAL_LOWER_EDGE-INTERVAL_UPPER_EDGE: Dist RecHitSet - 1st Interaction, Wrong assoc. Distance between a reconstructed hit set and the first interaction, only for those sets with wrong association (other hit is closest to the first interaction than the one identifed by the algorithm)
  • DetCompton:Classif: NRecHit: INTERVAL_LOWER_EDGE-INTERVAL_UPPER_EDGE: Dist RecHitSet - 1st Interaction, Good assoc. Distance between a reconstructed hit set and the first interaction, only for those sets with good association

These four histograms are filled for each interval of the variable, being INTERVAL_LOWER_EDGE the lower edge of each interval and INTERVAL_UPPER_EDGE the upper edge of each interval. The intervals are set by the parameters:

/gamos/setParam DetCClassifEnergy:Min MIN_VALUE

/gamos/setParam DetCClassifEnergy:Max MAX_VALUE

/gamos/setParam DetCClassifEnergy:Step STEP_VALUE

The intervals will be built between MIN_VALUE and MAX_VALUE each STEP_VALUE. Alternatively you can define the N+1 limits of the N intervals with the parameter:

/gamos/setParam DetCClassifEnergy:IntervalLimits VALUE_1 VALUE_2 … VALUE_N+1

The same set of histograms can done for the maximum, minimum or average reconstructed hit energy, and minimum position X, Y, Z, XY ( sqrt(X*X+Y*Y) ), XZ, YZ, XYZ.

To fill these histograms the following parameter should be used:

/gamos/setParam DetComptonStudyHistosUA:ClassificationVariables VAR1 VAR2 … VARN

where the variables can be EnergyMin, EnergyMax, EnergyAverage, XPosMin, YPosMin, ZPosMin, XYPosMin, XZPosMin, YZPosMin or XYZPosMin.

Another variable, but in this case it is always filled is the variable ALL, which indeed is not a variable and it makes no classification, but fill all reconstructed hit in a unique set of histograms.

The classes that manage this histograms are plug-ins, so that other variables could easily be created.

Histograms of gammas at sensitive detectors

We describe here the GmHistosGammaAtSD, an utility that prints information about the interaction of ‘original’ gammas in the sensitive detector. It is a user action and therefore it can be activated with the command

/gamos/userAction GmHistosGammaAtSD

At the end of run a table like the following one is printed:

$$$$$$$$$$ Classification of Gamma Interactions in SD $$$$$$$$
$$GC: nEvents       : 1000000
$$GC: n gamma in SD : 516170 : 25.8085 %
$$GC:  n PE           : 321588 : 62.30273 %
$$GC:    PE 0 COMP      : 157614 : 49.011157 %
$$GC:    PE 1 COMP      : 111304 : 34.610744 %
$$GC:    PE >1 COMP     : 52670 : 16.378099 %
$$GC:  n COMP         : 222590 : 43.12339 %
$$GC:    1 COMP         : 159193 : 71.518487 %
$$GC:    >1 COMP        : 63397 : 28.481513 %
$$GC: nGamma_COMP/event: 0.585931
  • nEvents : total number of events.
  • n gamma in SD : number of ‘original’ gammas reaching one sensitive detector.
  • n PE : number of ‘original’ gammas with photoelectric interaction in SD
  • PE 0 COMP : number of ‘original’ gammas with photoelectric interaction and no Compton interactions in SD (percentage relates to ``n PE’‘).
  • PE 1 COMP : number of ‘original’ gammas with photoelectric interaction and one Compton interaction in SD (percentage relates to ``n PE’‘).
  • PE >1 COMP : number of ‘original’ gammas with photoelectric interaction and more than one Compton interaction in SD (percentage relates to ``n PE’‘).
  • n COMP : number of ‘original’ gammas with Compton interactions in SD.
  • 1 COMP : number of ‘original’ gammas with only one Compton interaction in SD (percentage relates to ``n COMP’‘).
  • >1 COMP : number of ‘original’ gammas with more than one Compton interactions in SD (percentage relates to ``n COMP’‘).
  • nCOMP/event : number of Compton interactions in SD per ‘original’ gamma.

This class also makes several histograms, all of them have the words Gamma At SD: included in their names. And all are written to the file gammaSD.root/csv. The following histograms are defined:

  • Event Type / 100, i.e. no Rayleigh counting (“Event Type”)
  • Number of photoelectric interactions per ‘original’ gamma (“N PhotoElec”)
  • Number of Compton interactions per ‘original’ gamma (“N Compton”)
  • Number of Rayleigh interactions per ‘original’ gamma (“N Rayleigh”)
  • Number of photoelectric interactions vs Number of Compton+Rayleigh interactions per ‘original’ gamma (“N PhotoElec vs Compton+Rayleigh”)
  • Energy of gamma upon entering SD(“Energy at entering SD (keV)”)
  • Energy lost in the photoelectric interactions (“Energy lost PhotoElec (keV)”);
  • Difference in position from the gamma entering SD to the point where a photoelectric interaction happened (“Diff pos when PhotoElec (mm)”)
  • Difference in direction from the gamma entering SD to the point where a photoelectric interaction happened (“Diff dir when PhotoElec (mm)”)
  • Difference in kinetic energy from the gamma entering SD to the point where a photoelectric interaction happened (“Diff energy when PhotoElec (mm)”)
  • Energy lost in the Compton interactions (“Energy lost Compton (keV)”);
  • Angle variation in the Compton interactions (“Angle variation Compton (mrad)”)
  • Energy lost in the Rayleigh interactions (must be 0.)(“Energy lost Rayleigh (eV)”)
  • Angle variation in the Rayleigh interactions (“Angle variation Rayleigh (mrad)”)

All these histograms are repeated for the different types of events (a prefix in the name marks the event type the histogram refers to):

  • “ALL: “: every event
  • “No PE: ” events where no photoelectric interaction occurred
  • “Only PE: ” events where only photoelectric interaction occurred
  • “PE + 1 Compt: ” events where one photoelectric interaction plus only one Compton interaction occurred
  • “PE + >1 Compt: ” events where one photoelectric interaction plus more than one Compton interaction occurred

Automatic determination of production cuts for a detector

The method used in GAMOS to determine the best production cuts is based on what we can call an ‘inverse reasoning’. We count each particle that reaches the sensitive detector and we calculate first the range of the particle in the region where it is created. Then we can know that if we put a range cut in that region smaller than the calculated range, that particle would not reach our target plane. We also compute the range of the mother particle in the region where it was created and the same consecutively for all the ancestors. We know then that if we set in any of the regions where each of the ancestor particles are created a cut smaller than the corresponding range, we would stop the chain of particles and therefore we would have no particle in the target plane. After running a good number of tracks we can know for each particle type and for each region which is the biggest range we can put if we do not want to lose any particle. Indeed we may allow to lose a small amount of particles if this speeds up our simulation. To know easily which is the biggest cut you can use to lose less than a given percentage of particles, GAMOS provides a set of plots (one per each particle type and per each region) and a simple script to get automatically the cut values.

One warning is due here: as mentioned above when a track reaches the target, its range fills a histogram, but also the range of all the ancestors of this track. It may happen then that when you set a certain cut and the above-mentioned script gives you how many tracks would be killed, more than one killed track correspond to the same track reaching the target (i.e., with a cut you kill the track that reaches the target and the parent track). Therefore you might have an over-counting of the number of tracks killed by a cut. To avoid this the total number of tracks (the last lines of output) is not computed as the sum of tracks in the region. This number uses a histogram that contains only one entry per track reaching the target, the one corresponding to the track with the smallest range. If you want to set a different cut for each region and are worried for this double counting, you may have a look at the histogram named “trackInfos per Track in target”, that plots per each track reaching the target how many track informations are kept in the histograms. Another useful histogram for this case may be the 2D histogram “trackInfo Region vs trackInfo Region”, that plots all the region number of all the pairs of track informations that correspond to the same track reaching the target (you can get a list of which region number corresponds to which region at the end of the standard output file).

To use this utility in GAMOS it is only needed to add this command in your script:

/gamos/userAction GmProdCutsStudyUA DetCutsStudyFilter

that will use as target condition that a track enters a sensitive detector

This command will produce at the end of run a table and a histogram file with the needed information. The table will contain the minimum range that can be applied for each region/particle/process not to lose any track reaching the target, and it will look like this

%%%%% PRODUCTION CUTS STUDY RESULTS
GmProdCutsStudyUA: REGION= DefaultRegionForTheWorld PARTICLE=
gamma PROCESS= ALL MIN RANGE= 353161.38
GmProdCutsStudyUA: REGION= DefaultRegionForTheWorld PARTICLE=
gamma PROCESS= eBrem MIN RANGE= 353161.38

To get the cuts values for not losing a given percentage of particles in the target plane you can execute the ROOT script that can be found at GamosCore/GamosPhysics/Cuts/getProdCutsEffect.C :

*root -b -p -q .x getProdCutsEffect.C++\\(\\"prodcuts.root\\",percentage\\)*

and look at the last lines of output, those that contain the word ‘FINAL’, like the following ones

PARTICLE= e+ FINAL= 17 / 19
PARTICLE= e- FINAL= 72 / 34185
PARTICLE= gamma FINAL= 72 / 34184

Automatic determination of user limits for a detector

The method used in GAMOS to determine the minimum range user limits is similar to the one used to determine the best production cuts. The main difference is that when a track reaches the target we do not have to look at the range it had when created, but at the range it had in every step. This is because even if we want the minimum step, the track may have crossed several regions and the smallest range may not correspond to the last step. What we do nevertheless is to consider only the last step when there are a set of contiguous steps in the same region. Also for the ancestor tracks we have to store the information of each step, starting of course with the one when the track that reached the target (or its n-th ancestor if we are looking at the (n+1)-th ancestor) was created.

The same warning as for the production cuts should be mentioned here, but in this case it is more than a mere warning: when a track reaches the target, we accumulate one-track information of the last step in each region, for each of the ancestor tracks. Therefore it is very likely that there are more than one track information per track reaching the target, and therefore there will be overcounting of the number of tracks killed by a cut. As for the production cuts you should keep an eye on this.

To use this utility in GAMOS it is only needed to add the command in your script:

/gamos/userAction GmMinRangeLimitsStudyUA DetCutsStudyFilter

what will produce at the end of run a table and a histogram file with the needed information.

root -b -p -q .x getMinRangeCutsEffect.C++\\(\\"prodcuts.root\\",percentage\\) \|& tee out

and look at the last lines of output, those that contain the word ‘FINAL’, like the following ones

PARTICLE= e+ FINAL= 17 / 19
PARTICLE= e- FINAL= 72 / 34185
PARTICLE= gamma FINAL= 72 / 34184