Analysis utilities

This is a set of utilities that may serve in the analysis of phase space and dose files, for example to sum phase space or dose files from different jobs, to get basic information of the file contents, to fill histograms out of the files or to compare files from two jobs.

All these utilities are under the directory analysis. They are compiled by default with the rest of GAMOS code and after that they are available as executables as mentioned in the following subsections.

Summing phase space files

This utility serves to sum phase space files corresponding to different jobs with the same setup. To use it you have to write a file containing the list of header phase space files, one file per line, for example


	ps.20000.IAEAheader
	ps.20001.IAEAheader
	ps.20002.IAEAheader

Then you just have to run the executable

sumPS INPUT_FILE_LIST_NAME OUTPUT_FILE_NAME

where INPUT_FILE_LIST_NAME is the name of the file containing the list of files to add and OUTPUT_FILE_NAME is the name of the output file that will contain the sum of all the files (two files indeed as usual for IAEA phase space files: OUTPUT_FILE_NAME.IAEAheader and OUTPUT_FILE_NAME.IAEAphsp ).

When running you will see on the screen something similar to this:

	
	Opening phase space contained in     ps.20000.IAEAheader of type IAEA
	PARTICLES 225437 NPART_TOT 225437 NPARTORIG_TOT 5000000 
     RATIO 0.0450874 +- 0.000101661 RATIO_TOT 0.0450874
	Opening phase space contained in     ps.20001.IAEAheader of type IAEA
 	PARTICLES 224635 NPART_TOT 450072 NPARTORIG_TOT 10000000 
     RATIO 0.044927 +- 0.000101455 RATIO_TOT 0.0450072
	Opening phase space contained in     ps.20002.IAEAheader of type IAEA
	PARTICLES 224813 NPART_TOT 674885 NPARTORIG_TOT 15000000 
     RATIO 0.0449626 +- 0.000101501 RATIO_TOT 0.0449923
	* * * * N Particles   = 674885
	* * * * N Photons   = 673804
	* * * * N Electrons = 1057
	* * * * N Positrons = 24
	* * * * N Original Histories =  1.5e+07
      

For each phase space files after the name of the file comes a line with the file statistics: number of particles, accumulated number of particles of all files, accumulated number of original histories, ratio particles/original histories +- ratio error, accumulated ratio particles/original histories. And at the end the statistics on the total number of particles in total, photons, electron and positrons, and the total number of original histories in all summed files.

Making histograms out of a phase space file

You can make histograms out of the particles contained in a phase space file by running gamos with the input script that you can file in RadioTherapy/analysis/phaseSpace/analysePS/rt.analysePS.in. If you edit it and change the name of the input phase space file, at /gamos/setParam RTGeneratorPhaseSpace:FileName, and run

gamos rt.analysePS.in

You will get a file named phaseSpace.root, that contains the same histograms that you get when you run the job to write the phase space file (see section on Phase space histograms).

Alternatively you can directly analyse a phase space file with the executable analysePS. You just have to type in your shell:

analysePS FILE_NAME

where FILE_NAME is the name of the phase space header file (the one that ends with .IAEAheader; you may indeed omit the suffix .IAEAheader if you want).

When running you will see on the screen something similar to this:

	
	READING FILE ps.ubs.10x10.6
	
	ORIGINAL HISTORIES= 2000000
	N PARTICLES= 5.13303e+07 PER ORIG_HIST= 25.6651 +- 0.0184982
	N GAMMAS= 5.11882e+07 PER ORIG_HIST= 25.5941 +- 0.0184479
	N ELECTRONS= 137257 PER ORIG_HIST 0.0686285 +- 0.000191492
	N POSITRONS= 4839 PER ORIG_HIST= 0.0024195 +- 3.48235e-05
	Reading PHSP track 0
        ... 
	Reading PHSP track 51000000
	=== saving histograms in file === phaseSpace.root
      

The output contains the total number of original histories and after the number of particles, gammas, electrons and positrons in the phase space file, and these numbers divided by the number of original histories, with the error. Also an histogram file equal to the file that was created when producing the phase space file will be created.

Several arguments can be supplied to the executable in the standard Unix format:

Merging 'sqdose' files

This utility serves to merge dose files in the format sqdose corresponding to different jobs with the same setup, and getting the average dose of them. To use it you have to write a file containing the list of header phase space files, one file per line, for example


	sqddose.water.20000.out
	sqdose.water.20001.out
	sqdose.water.20002.out
      

Then you just have to run the executable

mergeSqdose INPUT_FILE_LIST_NAME OUTPUT_FILE_NAME

where INPUT_FILE_LIST_NAME is the name of the file containing the list of files to merge and OUTPUT_FILE_NAME is the name of the output file that will contain the merging of all the files.

When merging files it will be checked that they correspond to the same phantom by checking the number of voxels and voxel limits.

Merging '3ddose' files

This utility serves to merge dose files in the format 3ddose corresponding to different jobs with the same setup, and getting the average dose of them. To use it you have to write a file containing the list of header phase space files, one file per line, for example

3ddose.water.20000.out
	3ddose.water.20001.out
	3ddose.water.20002.out
      

Then you just have to run the executable

merge3ddose INPUT_FILE_LIST_NAME OUTPUT_FILE_NAME

where INPUT_FILE_LIST_NAME is the name of the file containing the list of files to merge and OUTPUT_FILE_NAME is the name of the output file that will contain the merging of all the files.

When merging files it will be checked that they correspond to the same phantom by checking the number of voxels and voxel limits.

Making histograms out of a 'sqdose' file

You can make histograms out of dose information contained in a sqdose file by running

analyseSqdose SQDOSE_FILE_NAME

When running you will see on the screen something similar to this:


	READING FILE sqdose.ubs.10x10.6.out
	GmSqdoseHeader::Read NEvent 2.5e+08
	GmSqdoseHeader::Read NVoxels 21 112 150
	GmSqdose::Read type 1
	USING std::map to store doses
	Number of voxels= 352800
	RTPSPDoseHistos nvoxel 21 112 150
	RTPSPDoseHistos dim 1 1 1
	RTPSPDoseHistos translation (0,0,-4)
	RTPSPDoseHistos rotation 
	[ (           1             0             0)
	(           0             1             0)
	(           0             0             1) ]
	MINIMUM DOSE 2.352e-15
	MAXIMUM DOSE 4.76001e-13
	RTPSPDoseHistos AVERAGE ERROR 20% = 0.0139839
	RTPSPDoseHistos AVERAGE ERROR 50% = 0.0111726
	RTPSPDoseHistos AVERAGE ERROR 90% = 0.00996564
	=== saving histograms in file === dose_analyseSqdose.root
	N EVENTS IN SOURCE 2.5e+08
      

First it is printer the sqdose file header information: number of events (original events run) and number of voxels in X, Y and Z. Then the type of the sqdose (see above). Then the STL container that will be used to store the doses and dose errors (see argument descriptions below). Then the total number of voxels to be read and the information from the class RTPSPDoseHistos: number of voxels in X, Y and Z, voxel dimensions in X, Y and Z, translation and rotation applied to phantom when dose file was written and the statistics, i.e. minimum and maximum voxel dose, and average dose error in 20%, 50% and 90% highest dose voxels.

You will also get a file named dose_analyseSqdose.root, that contains the same histograms that you get when you run the job to write the dose file using the scorer printer RTPSPDoseHistos.

Several arguments can be supplied to the executable in the standard Unix format:

Automatic determination of production cuts for an accelerator simulation

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 a given Z plane (corresponding to the phantom surface) 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 abovementioned 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 overcounting 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 RTCutsStudyFilter

or that will use as target condition that a track reaches a plane perpendicualr to the Z axis defined with the parameters

/gamos/setParam RTCutsStudyFilter:PlaneZ ZPOS

/gamos/setParam RTCutsStudyFilter:PlaneXDim XDIM

/gamos/setParam RTCutsStudyFilter:PlaneYDim YDIM

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/GamosCuts/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 production cuts for a dose in a phantom simulation

The use of production cuts in a dose computation may introduce a bias when a particle is killed (and then its energy is deposited locally) and it has enough energy to reach the next phantom voxel, or enough energy to create a particle that reaches the next phantom voxel (this happens mainly for electrons creating gammas, which have a much higher range).

To calculate automatically the best production cut, that is the one that gives the smallest CPU while biasing the dose computation a minimal amount, GAMOS uses an inverse reasoning. For a given set of cuts for electron and gamma it does not apply them but tags the particles that would have been killed by them. It also tags the voxel in which the particle is produced and then it computes all the dose deposited by the tagged particle or any of its children in a voxel that is not the same as the tagged voxel.

To use this utility in GAMOS you just have to associate to your dose scorer a filter of type GmProdCutOutsideVoxelFilter, passing to it as arguments the gamma cut and the electron cut, like in the following example

/gamos/scoring/addFilter2Scorer ProdCutFilter GmProdCutOutsideVoxelFilter PDDscorerPC10.1. 10.*mm 1.*mm

You should add another scorer without filter to get the total dose. After running your job with as many scorer-filter combinations as you like, you can look at the total dose deposited by each scorer and compare it with the total dose. It may happen that the dose lost with certain cuts, despite being a small proportion of the total dose, is distributed in a different manner than the total dose, introducing some bias in some region that you consider not acceptable. To check in detail the dose produced with a certain filter you can add a scorer printer of type RTPSPDoseHistos, that will produce several histograms of the dose (PDD, X & Y profiles, dose, dose-volume):

/gamos/scoring/addPrinter2Scorer RTPSPDoseHistos DoseScorerPC10.1.

The name of the printer will be passed to the name of the file containing the histograms.

Automatic determination of user limits for an accelerator simulation

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 only consider 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 RTCutsStudyFilter

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

Automatic determination of user limits for a dose in phantom simulation

The method used in GAMOS to determine the minimum range user limits is similar to the one used to determine the best production cut.

To use this utility in GAMOS you just have to associate to your dose scorer a filter of type GmProdCutOutsideVoxelFilter, passing to it as arguments the gamma cut and the electron cut, like in the following example

/gamos/filter ProdCutFilter GmMinRangeCutOutsideVoxelFilter 10.*mm 1.*mm

/gamos/scoring/addFilter2Scorer ProdCutFilter DoseScorerPC10.1.

After running your job with as many scorer-filter combinations as desired, you can look at the total dose deposited and compare it with the dose with very small cuts. It may happen that the dose with certain cuts, despite being a small number, is distributed in a different manner than with very small cuts, introducing some bias that you consider not acceptable. To check in detail the dose produced with a certain filter you can add a scorer printer of type RTPSPDoseHistos, that will produce several histograms of the dose (PDD, X & Y profiles, dose, dose-volume):

/gamos/scoring/addPrinter2Scorer RTPSPDoseHistos DoseScorerPC10.1.

The name of the printer will be passed to the name of the file containing the histograms.