Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/doiPET/src/doiPETAnalysis.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/doiPET/src/doiPETAnalysis.cc (Version 11.3.0) and /examples/advanced/doiPET/src/doiPETAnalysis.cc (Version 11.1.3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25                                                    25 
 26 //GEANT4 - Depth-of-Interaction enabled Positr     26 //GEANT4 - Depth-of-Interaction enabled Positron emission tomography (PET) advanced example 
 27                                                    27 
 28 //Authors and contributors                         28 //Authors and contributors
 29                                                    29 
 30 // Author list to be updated, with names of co     30 // Author list to be updated, with names of co-authors and contributors from National Institute of Radiological Sciences (NIRS)
 31                                                    31 
 32 // Abdella M. Ahmed (1, 2), Andrew Chacon (1,      32 // Abdella M. Ahmed (1, 2), Andrew Chacon (1, 2), Harley Rutherford (1, 2),
 33 // Hideaki Tashima (3), Go Akamatsu (3), Akram     33 // Hideaki Tashima (3), Go Akamatsu (3), Akram Mohammadi (3), Eiji Yoshida (3), Taiga Yamaya (3)
 34 // Susanna Guatelli (2), and Mitra Safavi-Naei     34 // Susanna Guatelli (2), and Mitra Safavi-Naeini (1, 2)
 35                                                    35 
 36 // (1) Australian Nuclear Science and Technolo     36 // (1) Australian Nuclear Science and Technology Organisation, Australia
 37 // (2) University of Wollongong, Australia         37 // (2) University of Wollongong, Australia
 38 // (3) National Institute of Radiological Scie     38 // (3) National Institute of Radiological Sciences, Japan
 39                                                    39 
 40 //Implemetation of the doiPETAnalysis.cc class     40 //Implemetation of the doiPETAnalysis.cc class
 41 //This implementation mimics readout (or digit     41 //This implementation mimics readout (or digitizer) of the PET scanner. To mimic realistic PET detector, the signals are blurred. Blurring 
 42 //parameters are given in inputParameters.txt      42 //parameters are given in inputParameters.txt file. Crystal based energy resolution and quantum efficiency has been applied. Deadtime (on 
 43 //each detector block and axially multiplexed      43 //each detector block and axially multiplexed detector) is also applied before the event is rejected by the energy window. The units for 
 44 //blurring parameters are in keV (for energy)      44 //blurring parameters are in keV (for energy) and ns (nano sec) for dead time. If the units are different, exception will be thrown and the
 45 //program quits. In this class, an ideal PMT i     45 //program quits. In this class, an ideal PMT is assumed to be placed at the corners of the crystal block. First, the ideal interaction position
 46 //(obtained by G4) is used to determine the di     46 //(obtained by G4) is used to determine the distance of the PMT from the interaction point. The signal on each PMT depends on the lateral (2D)
 47 //distance from the PMTs. Light sharing method     47 //distance from the PMTs. Light sharing method (reflector based) DOI identification method has been used. If the crystal ID is out of bound, 
 48 //error message will be displayed and the even     48 //error message will be displayed and the event will be rejected. The output file is single based list-mode ASCII file and can be then be 
 49 //processed into coinsidence list-mode data. A     49 //processed into coinsidence list-mode data. As, an option, binary output method is also given.  
 50 //Explanation is given for the methods provide     50 //Explanation is given for the methods provided. 
 51                                                    51 
 52                                                    52 
 53 #include "doiPETAnalysis.hh"                       53 #include "doiPETAnalysis.hh"
 54 #include "G4SystemOfUnits.hh"                      54 #include "G4SystemOfUnits.hh"
 55 #include "G4PhysicalConstants.hh"                  55 #include "G4PhysicalConstants.hh"
 56 #include <iomanip>                                 56 #include <iomanip>
 57 #include "Randomize.hh"                            57 #include "Randomize.hh"
 58 #include "G4SPSRandomGenerator.hh"                 58 #include "G4SPSRandomGenerator.hh"
 59 #include "doiPETAnalysisMessenger.hh"              59 #include "doiPETAnalysisMessenger.hh"
 60 #include "G4UnitsTable.hh"                         60 #include "G4UnitsTable.hh"
 61 #include "globals.hh"                              61 #include "globals.hh"
 62                                                    62 
 63 doiPETAnalysis* doiPETAnalysis::instance=0;        63 doiPETAnalysis* doiPETAnalysis::instance=0;
 64                                                    64 
 65 /////////// Constructor //////////////////////     65 /////////// Constructor /////////////////////////////////////////////
 66 doiPETAnalysis::doiPETAnalysis()                   66 doiPETAnalysis::doiPETAnalysis()
 67 {                                                  67 {
 68   factoryOn = false; //G4ROOT                      68   factoryOn = false; //G4ROOT
 69   // Initialization ntuple                         69   // Initialization ntuple
 70     for (G4int k=0; k<MaxNtCol; k++) fNtColId[     70     for (G4int k=0; k<MaxNtCol; k++) fNtColId[k] = 0;
 71                                                    71 
 72   fAnalysisMessenger = new doiPETAnalysisMesse     72   fAnalysisMessenger = new doiPETAnalysisMessenger(this);
 73                                                    73 
 74   //Set energy window                              74   //Set energy window
 75   lowerThreshold = 400*keV;                        75   lowerThreshold = 400*keV;
 76   upperThreshold = 600*keV;                        76   upperThreshold = 600*keV;
 77   triggerEnergy = 50*eV;                           77   triggerEnergy = 50*eV;
 78                                                    78   
 79   //give default initial activity. Activity st     79   //give default initial activity. Activity strength is changed in the .mac file
 80   InitialActivity = 1000000*becquerel;             80   InitialActivity = 1000000*becquerel;
 81                                                    81 
 82   //In NEMA NU2, all test is done with F-18        82   //In NEMA NU2, all test is done with F-18
 83   halfLife = 109.771*60 * s;//The halfLife of      83   halfLife = 109.771*60 * s;//The halfLife of a given isotope can be changed via the run.mac file
 84                                                    84 
 85   //                                               85   //
 86   totalTime = 0 * s;                               86   totalTime = 0 * s;
 87   prev_totalTime = 0 * s;                          87   prev_totalTime = 0 * s;
 88   prev_eventID = 0;                                88   prev_eventID = 0;
 89                                                    89 
 90   //                                               90   //
 91   //Initialize crystal ID                          91   //Initialize crystal ID
 92   crystalIDNew_tan = -1;                           92   crystalIDNew_tan = -1;
 93   crystalIDNew_axial = -1;                         93   crystalIDNew_axial = -1;
 94   crystalIDNew_DOI = -1;                           94   crystalIDNew_DOI = -1;
 95                                                    95 
 96   //                                               96   //
 97   scatterIndex = 0;                                97   scatterIndex = 0;
 98                                                    98 
 99   //                                               99   //
100   numberOfPixel_tan = 2*numberOfCrystal_tangen    100   numberOfPixel_tan = 2*numberOfCrystal_tangential; //32;
101   numberOfPixel_axial = 2*numberOfCrystal_axia    101   numberOfPixel_axial = 2*numberOfCrystal_axial; //32;
102                                                   102 
103   //Default value for deadtime.                   103   //Default value for deadtime.
104   block_DeadTime = 256*ns;                        104   block_DeadTime = 256*ns;
105   module_DeadTime = 0*ns;                         105   module_DeadTime = 0*ns;
106   //                                              106   //
107                                                   107 
108   //Crystal blurring parameters. One detector     108   //Crystal blurring parameters. One detector has 1024 crystals. All the crystals have different energy resolution. 
109   //So, a range of energy resolution is applie    109   //So, a range of energy resolution is applied between minumun and maximum values. 
110   //The energy resolution can be set in the in    110   //The energy resolution can be set in the inputParameter.txt file
111   crystalResolutionMin = 0.13;//13%               111   crystalResolutionMin = 0.13;//13%
112   crystalResolutionMax = 0.17;//17%               112   crystalResolutionMax = 0.17;//17%
113                                                   113 
114   crystalEnergyRef = 511 * keV;//Energy of ref    114   crystalEnergyRef = 511 * keV;//Energy of reference in which the energy resolution of the crystal is computed
115                                                   115 
116   //The quantum efficiency models the probabil    116   //The quantum efficiency models the probability for the event to be detected by the photo-detector.
117   //The quantum efficiency can be set in the i    117   //The quantum efficiency can be set in the inputParameter.txt file
118   crystalQuantumEfficiency = 1;//100%             118   crystalQuantumEfficiency = 1;//100% 
119   //                                              119   //
120                                                   120 
121   //intialize deadtime for blocks and modules     121   //intialize deadtime for blocks and modules
122   numberOfBlocks_total = numberOfRings * numbe    122   numberOfBlocks_total = numberOfRings * numberOfDetector_perRing; 
123   blockTime = new double[numberOfBlocks_total]    123   blockTime = new double[numberOfBlocks_total];//for each individual block.
124   moduleTime = new double[numberOfBlocks_total    124   moduleTime = new double[numberOfBlocks_total];//for axially multiplexed detectors.
125                                                   125 
126   //Initialize the deadtime for each detector     126   //Initialize the deadtime for each detector and axially multiplexed detector (also called modules)
127   for(G4int i = 0; i<numberOfBlocks_total; i++    127   for(G4int i = 0; i<numberOfBlocks_total; i++){
128     blockTime [i] = 0.0;                          128     blockTime [i] = 0.0;
129     moduleTime [i] = 0.0;                         129     moduleTime [i] = 0.0;
130   }                                               130   }
131                                                   131 
132   //Initialize type of output. The default out    132   //Initialize type of output. The default output is single events
133   getSinglesData  = false; //default value        133   getSinglesData  = false; //default value 
134   getCoincidenceData = false;                     134   getCoincidenceData = false;
135                                                   135 
136   ApplyAngerLogic = true;                         136   ApplyAngerLogic = true;
137   isDOIlookUpTablePrepared = false;               137   isDOIlookUpTablePrepared = false;
138                                                   138 
139   numberOfHit = 0;                                139   numberOfHit = 0;
140                                                   140 
141   //This value is based on the assumption that    141   //This value is based on the assumption that the shift of the response due to the reflector is half distance from the interaction position to the air gap.
142   shiftCoeff = 0.5;                               142   shiftCoeff = 0.5;
143                                                   143 
144   //                                              144   //
145   PMTblurring_tan = 0.0;                          145   PMTblurring_tan = 0.0;
146   PMTblurring_axial = 0.0;                        146   PMTblurring_axial = 0.0;
147 }                                                 147 }
148 ////////// Destructor ////////////////////////    148 ////////// Destructor ///////////////////////////////////////////////
149 doiPETAnalysis::~doiPETAnalysis()                 149 doiPETAnalysis::~doiPETAnalysis()
150 {                                                 150 {
151   delete fAnalysisMessenger;                      151   delete fAnalysisMessenger;
152   delete [] blockTime;                            152   delete [] blockTime;
153   delete [] moduleTime;                           153   delete [] moduleTime;
154 }                                                 154 }
155                                                   155 
156 ////////// GetInstance ///////////////////////    156 ////////// GetInstance /////////////////////////////////////////////
157 doiPETAnalysis* doiPETAnalysis::GetInstance()     157 doiPETAnalysis* doiPETAnalysis::GetInstance()
158 {                                                 158 {
159   if(instance==0) instance = new doiPETAnalysi    159   if(instance==0) instance = new doiPETAnalysis();
160   return instance;                                160   return instance;
161 }                                                 161 }
162 void doiPETAnalysis::Delete()                     162 void doiPETAnalysis::Delete()
163 {                                                 163 {
164   delete instance;                                164   delete instance;
165 }                                                 165 }
166                                                   166 
167 //If there is energy deposition in the phantom    167 //If there is energy deposition in the phantom by the photon, the scatter index is 1, otherwise it is 0
168 //Use this for checking                           168 //Use this for checking
169 void doiPETAnalysis::SetScatterIndexInPhantom(    169 void doiPETAnalysis::SetScatterIndexInPhantom(G4int sci){
170   scatterIndex = sci;                             170   scatterIndex = sci;
171 }                                                 171 }
172                                                   172 
173 //Get the source position if the process is an    173 //Get the source position if the process is annihilation.
174 //Use this for checking                           174 //Use this for checking
175 void doiPETAnalysis::SetSourcePosition(G4Three    175 void doiPETAnalysis::SetSourcePosition(G4ThreeVector spos){
176   spositionX = spos.x();                          176   spositionX = spos.x();
177   spositionY = spos.y();                          177   spositionY = spos.y();
178   spositionZ = spos.z();                          178   spositionZ = spos.z();
179 }                                                 179 }
180                                                   180 
181                                                   181 
182 //Set the event ID                                182 //Set the event ID
183 //Use this for checking. eventID should not be    183 //Use this for checking. eventID should not be used to sort coincidence events if realistic simulation is done
184 void doiPETAnalysis::SetEventID(G4int evID){      184 void doiPETAnalysis::SetEventID(G4int evID){
185   eventID = evID;                                 185   eventID = evID;
186 }                                                 186 }
187                                                   187 
188 //                                                188 //
189                                                   189 
190 void doiPETAnalysis::GetSizeOfDetector(G4doubl    190 void doiPETAnalysis::GetSizeOfDetector(G4double detSizeDoi, G4double detSizeTan, G4double detSizeAxial){
191   sizeOfDetector_DOI = detSizeDoi;                191   sizeOfDetector_DOI = detSizeDoi;
192   sizeOfDetector_axial = detSizeTan;              192   sizeOfDetector_axial = detSizeTan;
193   sizeOfDetector_tangential = detSizeAxial;       193   sizeOfDetector_tangential = detSizeAxial;
194 }                                                 194 }
195                                                   195 
196 //                                                196 //
197 void doiPETAnalysis::SetActivity(G4double newA    197 void doiPETAnalysis::SetActivity(G4double newActivity){
198   InitialActivity = newActivity;                  198   InitialActivity = newActivity;
199   G4cout<<"Initial activity: "<<InitialActivit    199   G4cout<<"Initial activity: "<<InitialActivity/becquerel<<" Bq."<<G4endl;
200 }                                                 200 }
201 void doiPETAnalysis::SetIsotopeHalfLife(G4doub    201 void doiPETAnalysis::SetIsotopeHalfLife(G4double newHalfLife){
202   halfLife = newHalfLife;                         202   halfLife = newHalfLife;
203   G4cout<<"Half life of the isotope "<<halfLif    203   G4cout<<"Half life of the isotope "<<halfLife/s<<" sec."<<G4endl;
204 }                                                 204 }
205                                                   205 
206 //The time is based on random time intervals b    206 //The time is based on random time intervals between events. See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3267383/
207 //This time mimics acquisition time of a PET s    207 //This time mimics acquisition time of a PET scanner for a given number of particles. 
208 void doiPETAnalysis::CalulateAcquisitionTime()    208 void doiPETAnalysis::CalulateAcquisitionTime(){
209   //Calculate the strength of activity at t=to    209   //Calculate the strength of activity at t=totaltime using decay equation 
210   activityNow = InitialActivity * std::exp(-((    210   activityNow = InitialActivity * std::exp(-((0.693147/halfLife)*totalTime)); //ln(2) = 0.693147181
211                                                   211 
212   //Activity based time interval.                 212   //Activity based time interval. 
213   timeInterval = -std::log(G4UniformRand())*(1    213   timeInterval = -std::log(G4UniformRand())*(1./activityNow);
214   totalTime = timeInterval+prev_totalTime;        214   totalTime = timeInterval+prev_totalTime;
215   prev_totalTime = totalTime;                     215   prev_totalTime = totalTime; 
216 }                                                 216 }
217                                                   217 
218 //Apply energy blurring on the crystals. The v    218 //Apply energy blurring on the crystals. The value of the energy blurring with respect to a reference energy is given in the inputParameter.txt file
219 G4double doiPETAnalysis::QuantumEffifciency(G4    219 G4double doiPETAnalysis::QuantumEffifciency(G4double edep, G4int blkID, G4int cysID)
220 {                                                 220 {
221   if(fixedResolution){                            221   if(fixedResolution){
222     crystalResolution = energyResolution_fixed    222     crystalResolution = energyResolution_fixed;
223   }                                               223   }
224   else{                                           224   else{
225     crystalResolution = energyResolution_cryDe    225     crystalResolution = energyResolution_cryDependent[blkID][cysID];
226   }                                               226   }
227   crystalCoeff = crystalResolution * std::sqrt    227   crystalCoeff = crystalResolution * std::sqrt(crystalEnergyRef);
228                                                   228 
229   G4double QE = G4UniformRand();                  229   G4double QE = G4UniformRand();
230                                                   230 
231   //The quantum efficiency models the probabil    231   //The quantum efficiency models the probability for the event to be detected by the photo-detector. It can be changed in the inputParameter.txt file
232   if(QE <= crystalQuantumEfficiency)              232   if(QE <= crystalQuantumEfficiency)
233   {                                               233   {
234     edep_AfterCrystalBlurring = G4RandGauss::s    234     edep_AfterCrystalBlurring = G4RandGauss::shoot(edep,crystalCoeff*std::sqrt(edep)/2.35);
235   }                                               235   }
236   else {                                          236   else {
237     //not detected by the photodetector, event    237     //not detected by the photodetector, eventhough there was an interaction
238     edep_AfterCrystalBlurring = 0 *keV;           238     edep_AfterCrystalBlurring = 0 *keV;
239   }                                               239   }
240   return edep_AfterCrystalBlurring;               240   return edep_AfterCrystalBlurring;
241 }                                                 241 }
242                                                   242 
243 ///////// ReadOut ////////////////////////////    243 ///////// ReadOut ///////////////////////////////////
244                                                   244 
245 void doiPETAnalysis::ReadOut(G4int blkID, G4in    245 void doiPETAnalysis::ReadOut(G4int blkID, G4int cryID, G4double interTime, G4double timeAnnih, G4ThreeVector interPos, G4double edep)
246 {                                                 246 {
247   blockID = blkID;                                247   blockID = blkID;
248   crystalID = cryID;                              248   crystalID = cryID;
249   interactionTime = interTime;                    249   interactionTime = interTime;
250   time_annihil = timeAnnih;                       250   time_annihil = timeAnnih;
251   interactionPos = interPos;                      251   interactionPos = interPos;
252   totalEdep = edep;                               252   totalEdep = edep;
253                                                   253 
254   //Get the time of flight. This is the durati    254   //Get the time of flight. This is the duration from the annihilation process to the detection of the photons by the scintillator. 
255   time_tof = interactionTime - time_annihil;      255   time_tof = interactionTime - time_annihil;
256                                                   256 
257   //time of the event when detected (timerTag)    257   //time of the event when detected (timerTag)
258   timeStamp = totalTime + time_tof;               258   timeStamp = totalTime + time_tof;
259                                                   259 
260   //triggerEnergy is the energy deposited in t    260   //triggerEnergy is the energy deposited in the detector below which the detector is insensitive to any interaction.
261   if(totalEdep<triggerEnergy)return;              261   if(totalEdep<triggerEnergy)return;
262                                                   262 
263   //************************************** App    263   //************************************** Apply dead-time ********************************************//
264   //Apply paralizable dead-time in the block b    264   //Apply paralizable dead-time in the block beofore events are rejected by the energy window
265   if(std::fabs(timeStamp - blockTime[blockID])    265   if(std::fabs(timeStamp - blockTime[blockID]) >=  block_DeadTime){ //If true, the event is accepted
266     blockTime[blockID] = timeStamp;               266     blockTime[blockID] = timeStamp;
267   }                                               267   }
268   else {                                          268   else {
269     //If the time difference is less than the     269     //If the time difference is less than the processing time of the detector (dead time), then the dead time (blockTime) of the block is extended.
270     blockTime[blockID] = timeStamp;               270     blockTime[blockID] = timeStamp;
271                                                   271 
272     //the event is then rejected                  272     //the event is then rejected
273     //continue;                                   273     //continue;
274     return;                                       274     return;
275   }                                               275   }
276                                                   276 
277   //Apply Non-paralyzable dead-time on axially    277   //Apply Non-paralyzable dead-time on axially multiplexed detectors (4 detectors are arranged axailly)
278   //If the time difference is less than the pr    278   //If the time difference is less than the processing time of the module,  the event is rejected without extending the dead time of the module
279   if(std::fabs(timeStamp - moduleTime[blockID]    279   if(std::fabs(timeStamp - moduleTime[blockID]) > module_DeadTime){
280                                                   280 
281     //The following finds the block id's of fo    281     //The following finds the block id's of four blocks which are arranged axially
282     for (G4int r_ring = 0; r_ring < numberOfRi    282     for (G4int r_ring = 0; r_ring < numberOfRings; r_ring++){
283       if (blockID >= r_ring*numberOfDetector_p    283       if (blockID >= r_ring*numberOfDetector_perRing && blockID <(r_ring + 1)*numberOfDetector_perRing){
284         for (G4int m_module = 0; m_module < nu    284         for (G4int m_module = 0; m_module < numberOfRings; m_module++){
285                                                   285 
286           //Set the time of the module (four b    286           //Set the time of the module (four blocks) the same
287           moduleTime[blockID + (m_module - r_r    287           moduleTime[blockID + (m_module - r_ring)*numberOfDetector_perRing] = timeStamp;
288         }                                         288         }
289       }                                           289       }
290     }                                             290     }
291   }                                               291   }
292   else return;                                    292   else return;
293                                                   293 
294                                                   294 
295   /////////////////////////   Write qualified     295   /////////////////////////   Write qualified single events based the energy deposition in the detector   ///////////
296                                                   296 
297   if(totalEdep>lowerThreshold && totalEdep<upp    297   if(totalEdep>lowerThreshold && totalEdep<upperThreshold ){
298                                                   298 
299     //identifiy the layer                         299     //identifiy the layer
300     DOI_ID =  G4int(crystalID/(numberOfCrystal    300     DOI_ID =  G4int(crystalID/(numberOfCrystal_tangential * numberOfCrystal_axial));
301                                                   301 
302     //identify the crystal id for each Layer.     302     //identify the crystal id for each Layer. Now, crystalID_2D can take  0,1, ... numberOfCrystal_tangential x numberOfCrystal_axial 
303     crystalID_2D = crystalID - (DOI_ID*numberO    303     crystalID_2D = crystalID - (DOI_ID*numberOfCrystal_tangential * numberOfCrystal_axial);
304                                                   304 
305     //identify the crystal ID in the tangentia    305     //identify the crystal ID in the tangential and axial direction
306     crystalID_axial = crystalID_2D/numberOfCry    306     crystalID_axial = crystalID_2D/numberOfCrystal_axial;
307     crystalID_tangential = crystalID_2D%number    307     crystalID_tangential = crystalID_2D%numberOfCrystal_tangential;
308                                                   308 
309     intPosX = interactionPos.x();                 309     intPosX = interactionPos.x();
310     intPosY = interactionPos.y();                 310     intPosY = interactionPos.y();
311     intPosZ = interactionPos.z();                 311     intPosZ = interactionPos.z();
312                                                   312 
313                                                   313 
314     if(ApplyAngerLogic){                          314     if(ApplyAngerLogic){
315       //shiftCoeff = 0.5 is used. This value i    315       //shiftCoeff = 0.5 is used. This value is based on the assumption that the shift due to the reflector is half distance from the interaction position to the air gap. 
316       AngerLogic(intPosX, intPosY, intPosZ, to    316       AngerLogic(intPosX, intPosY, intPosZ, totalEdep, shiftCoeff, isDOIlookUpTablePrepared);//
317     }                                             317     }
318                                                   318 
319     //Single event output. Coincidence events     319     //Single event output. Coincidence events can then be made using the single events.
320     if(getSinglesData) WriteOutput();             320     if(getSinglesData) WriteOutput();
321                                                   321 
322     //Coincidence output                          322     //Coincidence output
323     if(getCoincidenceData){                       323     if(getCoincidenceData){
324       eventID_coin.push_back(eventID);            324       eventID_coin.push_back(eventID);
325       blockID_coin.push_back(blockID);            325       blockID_coin.push_back(blockID);
326       cryID_axial_coin.push_back(crystalID_axi    326       cryID_axial_coin.push_back(crystalID_axial);
327       cryID_tan_coin.push_back(crystalID_tange    327       cryID_tan_coin.push_back(crystalID_tangential);
328       edep_coin.push_back(totalEdep);             328       edep_coin.push_back(totalEdep);
329       cryDOI_coin.push_back(DOI_ID);              329       cryDOI_coin.push_back(DOI_ID);
330       time_coin.push_back(timeStamp);             330       time_coin.push_back(timeStamp);
331                                                   331 
332       numberOfHit++;                              332       numberOfHit++;
333                                                   333 
334       if(numberOfHit == 2){ //two events withi    334       if(numberOfHit == 2){ //two events within the energy window are qualified.
335         WriteOutput();                            335         WriteOutput();
336         ResetNumberOfHits();                      336         ResetNumberOfHits();
337       }                                           337       }
338     }                                             338     }
339   }                                               339   }
340 }                                                 340 }
341                                                   341 
342                                                   342 
343                                                   343 
344 ////////// Clear /////////////////////////////    344 ////////// Clear ///////////////////////////////////////////////////
345 void doiPETAnalysis::ResetNumberOfHits()          345 void doiPETAnalysis::ResetNumberOfHits()
346 {                                                 346 {
347   numberOfHit = 0;                                347   numberOfHit = 0;
348   eventID_coin.clear();                           348   eventID_coin.clear();
349   blockID_coin.clear();                           349   blockID_coin.clear();
350   cryID_axial_coin.clear();                       350   cryID_axial_coin.clear();
351   cryID_tan_coin.clear();                         351   cryID_tan_coin.clear();
352   edep_coin.clear();                              352   edep_coin.clear();
353   cryDOI_coin.clear();                            353   cryDOI_coin.clear();
354   time_coin.clear();                              354   time_coin.clear();
355                                                   355 
356 }                                                 356 }
357                                                   357 
358 //                                                358 //
359 void doiPETAnalysis::Open(G4String fileName)      359 void doiPETAnalysis::Open(G4String fileName)
360 {                                                 360 {
361   if(getSinglesData){                             361   if(getSinglesData){
362     asciiFileName = fileName + "Singles.data";    362     asciiFileName = fileName + "Singles.data";
363     rootFileName = fileName+"Singles.root";       363     rootFileName = fileName+"Singles.root";
364   }                                               364   }
365   if(getCoincidenceData){                         365   if(getCoincidenceData){
366     asciiFileName = fileName + "Coincidence.da    366     asciiFileName = fileName + "Coincidence.data";
367     rootFileName = fileName+"Coincidence.root"    367     rootFileName = fileName+"Coincidence.root";
368   }                                               368   }
369                                                   369 
370   ofs.open(asciiFileName.c_str());                370   ofs.open(asciiFileName.c_str());
371   if(!ofs.is_open()){                             371   if(!ofs.is_open()){
372     G4cerr<<"=== \n File opening Error to writ    372     G4cerr<<"=== \n File opening Error to write the output ===="<<G4endl;
373     exit(0);                                      373     exit(0);
374   }                                               374   }
375                                                   375 
376 }                                                 376 }
377                                                   377 
378 void doiPETAnalysis::WriteOutput(){               378 void doiPETAnalysis::WriteOutput(){
379   if(getSinglesData){                             379   if(getSinglesData){
380     ofs<<eventID<<" "<<blockID<<" "<<std::setp    380     ofs<<eventID<<" "<<blockID<<" "<<std::setprecision(17)<<timeStamp/s<<" "<<std::setprecision(7)<<totalEdep/keV<<" "<<intPosX<<" "<<intPosY<<" "<<intPosZ<<" "<<spositionX<<" "<<spositionY<<" "<<spositionZ<<G4endl;
381     //ofs<<eventID<<" "<<blockID<<" "<<crystal    381     //ofs<<eventID<<" "<<blockID<<" "<<crystalID_axial<<" "<<crystalID_tangential<<" "<<DOI_ID<<" "<<std::setprecision(17)<<timeStamp/s<<" "<<std::setprecision(7)<<totalEdep/keV<<G4endl;
382                                                   382 
383   }                                               383   }
384   if(getCoincidenceData){                         384   if(getCoincidenceData){
385     //2 singles will qualify to be in coincide    385     //2 singles will qualify to be in coincidence within the energy window.
386     for(G4int i=0; i<2; i++){                     386     for(G4int i=0; i<2; i++){
387                                                   387 
388       //First Single                              388       //First Single
389       if(i==0){                                   389       if(i==0){
390         eventID0        = eventID_coin[0];        390         eventID0        = eventID_coin[0];
391         blockID0        = blockID_coin[0];        391         blockID0        = blockID_coin[0];
392         crystalID_axial0    = cryID_axial_coin    392         crystalID_axial0    = cryID_axial_coin[0];
393         crystalID_tangential0 = cryID_tan_coin    393         crystalID_tangential0 = cryID_tan_coin[0];
394         DOI_ID0         = cryDOI_coin[0];         394         DOI_ID0         = cryDOI_coin[0];
395         timeStamp0        = time_coin[0];         395         timeStamp0        = time_coin[0];
396         totalEdep0        = edep_coin[0];         396         totalEdep0        = edep_coin[0];
397       }                                           397       }
398       if(i==1){                                   398       if(i==1){
399         //Second Single                           399         //Second Single
400         eventID1        = eventID_coin[1];        400         eventID1        = eventID_coin[1];
401         blockID1        = blockID_coin[1];        401         blockID1        = blockID_coin[1];
402         crystalID_axial1    = cryID_axial_coin    402         crystalID_axial1    = cryID_axial_coin[1];
403         crystalID_tangential1 = cryID_tan_coin    403         crystalID_tangential1 = cryID_tan_coin[1];
404         DOI_ID1         = cryDOI_coin[1];         404         DOI_ID1         = cryDOI_coin[1];
405         timeStamp1        = time_coin[1];         405         timeStamp1        = time_coin[1];
406         totalEdep1        = edep_coin[1];         406         totalEdep1        = edep_coin[1];
407       }                                           407       }
408     }                                             408     }
409                                                   409 
410     ofs<<eventID0<<" "<<blockID0<<" "<<crystal    410     ofs<<eventID0<<" "<<blockID0<<" "<<crystalID_axial0<<" "<<crystalID_tangential0<<" "<<DOI_ID0<<" "<<std::setprecision(17)<<timeStamp0/s<<" "<<std::setprecision(7)<<totalEdep0/keV<<" "
411       <<eventID1<<" "<<blockID1<<" "<<crystalI    411       <<eventID1<<" "<<blockID1<<" "<<crystalID_axial1<<" "<<crystalID_tangential1<<" "<<DOI_ID1<<" "<<std::setprecision(17)<<timeStamp1/s<<" "<<std::setprecision(7)<<totalEdep1/keV<<" "
412       <<spositionX<<" "<<spositionY<<" "<<spos    412       <<spositionX<<" "<<spositionY<<" "<<spositionZ<<G4endl;
413                                                   413 
414   }                                               414   }
415 #ifdef ANALYSIS_USE                               415 #ifdef ANALYSIS_USE
416   FillListModeEvent();                            416   FillListModeEvent();
417 #endif                                            417 #endif
418                                                   418 
419 }                                                 419 }
420                                                   420 
421 //                                                421 //
422 ///////// Close //////////////////////////////    422 ///////// Close /////////////////////////////////////////////////////
423 void doiPETAnalysis::Close()                      423 void doiPETAnalysis::Close()
424 {                                                 424 {
425   //close ascii file                              425   //close ascii file
426   ofs.close();                                    426   ofs.close();
427                                                   427 
428 }                                                 428 }
429                                                   429 
430 //Place the photomultiplier tube (PMT) at each    430 //Place the photomultiplier tube (PMT) at each corner of the detector. The positions of the PMT is with respect to the axis of the detector block
431 //All the PMTs are placed at the same doi (x)     431 //All the PMTs are placed at the same doi (x) position (at +sizeOfDetector_DOI/2 which is at the top of the detector). 
432                                                   432 
433 //The PMT is placed at each corner of the crys    433 //The PMT is placed at each corner of the crystal block and is assumed to be an ideal PMT.
434 //The signal (energy deposition) of each PMT d    434 //The signal (energy deposition) of each PMT depends on  the distance of the respective PMT from the interaction point
435 void doiPETAnalysis::PMTPosition(){               435 void doiPETAnalysis::PMTPosition(){
436                                                   436 
437   sizeOfDetector_DOI = (numberOfCrystal_DOI *     437   sizeOfDetector_DOI = (numberOfCrystal_DOI * sizeOfCrystal_DOI) + (numberOfCrystal_DOI - 1)*crystalGap_DOI;
438   sizeOfDetector_axial = (numberOfCrystal_axia    438   sizeOfDetector_axial = (numberOfCrystal_axial * sizeOfCrystal_axial) + (numberOfCrystal_axial - 1)*crystalGap_axial;
439   sizeOfDetector_tangential = (numberOfCrystal    439   sizeOfDetector_tangential = (numberOfCrystal_tangential * sizeOfCrystal_tangential) + (numberOfCrystal_tangential - 1)*crystalGap_tangential;
440                                                   440 
441   //Position of PMT1.                             441   //Position of PMT1. 
442   posPMT1x = sizeOfDetector_DOI/2;//mm            442   posPMT1x = sizeOfDetector_DOI/2;//mm
443   posPMT1y = -sizeOfDetector_tangential/2;        443   posPMT1y = -sizeOfDetector_tangential/2;
444   posPMT1z = -sizeOfDetector_axial/2;             444   posPMT1z = -sizeOfDetector_axial/2;
445                                                   445 
446   //Position of PMT2                              446   //Position of PMT2
447   posPMT2x = sizeOfDetector_DOI/2;                447   posPMT2x = sizeOfDetector_DOI/2;
448   posPMT2y = sizeOfDetector_tangential/2;         448   posPMT2y = sizeOfDetector_tangential/2;
449   posPMT2z = -sizeOfDetector_axial/2;             449   posPMT2z = -sizeOfDetector_axial/2;
450                                                   450 
451   //Position of PMT3                              451   //Position of PMT3 
452   posPMT3x = sizeOfDetector_DOI/2;                452   posPMT3x = sizeOfDetector_DOI/2;
453   posPMT3y = -sizeOfDetector_tangential/2;        453   posPMT3y = -sizeOfDetector_tangential/2;
454   posPMT3z = sizeOfDetector_axial/2;              454   posPMT3z = sizeOfDetector_axial/2;
455                                                   455 
456   //Position of PMT4                              456   //Position of PMT4 
457   posPMT4x = sizeOfDetector_DOI/2;                457   posPMT4x = sizeOfDetector_DOI/2;
458   posPMT4y = sizeOfDetector_tangential/2;         458   posPMT4y = sizeOfDetector_tangential/2;
459   posPMT4z = sizeOfDetector_axial/2;              459   posPMT4z = sizeOfDetector_axial/2;
460                                                   460 
461   G4cout<<"PMT positions: "<<G4endl;              461   G4cout<<"PMT positions: "<<G4endl;
462   G4cout<<"PMT1 (mm) ("<<posPMT1x<<", "<<posPM    462   G4cout<<"PMT1 (mm) ("<<posPMT1x<<", "<<posPMT1y<<", "<<posPMT1z<<")"<<G4endl;
463   G4cout<<"PMT2 (mm) ("<<posPMT2x<<", "<<posPM    463   G4cout<<"PMT2 (mm) ("<<posPMT2x<<", "<<posPMT2y<<", "<<posPMT2z<<")"<<G4endl;
464   G4cout<<"PMT3 (mm) ("<<posPMT3x<<", "<<posPM    464   G4cout<<"PMT3 (mm) ("<<posPMT3x<<", "<<posPMT3y<<", "<<posPMT3z<<")"<<G4endl;
465   G4cout<<"PMT4 (mm) ("<<posPMT4x<<", "<<posPM    465   G4cout<<"PMT4 (mm) ("<<posPMT4x<<", "<<posPMT4y<<", "<<posPMT4z<<")"<<G4endl;
466                                                   466 
467 }                                                 467 }
468                                                   468 
469 //The blurring parameters are given and can be    469 //The blurring parameters are given and can be changed in the inputParameter.txt file
470 void doiPETAnalysis::BlurringParameters(){        470 void doiPETAnalysis::BlurringParameters(){
471   char inputChar[256];                            471   char inputChar[256];
472   std::string inputLine;                          472   std::string inputLine;
473   G4String value[7];                              473   G4String value[7];
474   std::string  filename = "inputParameter.txt"    474   std::string  filename = "inputParameter.txt";
475   ifs.open(filename.c_str());                     475   ifs.open(filename.c_str());
476   if(!ifs.good()){                                476   if(!ifs.good()){
477     G4cerr<<"File opening Error: Could not ope    477     G4cerr<<"File opening Error: Could not open "<<filename<<G4endl;
478     exit(0);                                      478     exit(0);
479   }                                               479   }
480   while(!ifs.eof()){                              480   while(!ifs.eof()){
481     ifs.getline(inputChar,256);                   481     ifs.getline(inputChar,256);
482     inputLine = inputChar;                        482     inputLine = inputChar;
483     if(inputChar[0]!='#' && inputLine.length()    483     if(inputChar[0]!='#' && inputLine.length()!=0 ){
484       if( (std::string::size_type)inputLine.fi    484       if( (std::string::size_type)inputLine.find("block_DeadTime:")!=std::string::npos){
485         std::istringstream tmpStream(inputLine    485         std::istringstream tmpStream(inputLine);
486         tmpStream >> value[0] >> value[1] >> v    486         tmpStream >> value[0] >> value[1] >> value[2];
487         block_DeadTime = atof(value[1].c_str()    487         block_DeadTime = atof(value[1].c_str());
488         if(value[2] != "ns"){                     488         if(value[2] != "ns"){
489           G4cerr<<" Dead time unit is not in n    489           G4cerr<<" Dead time unit is not in nano seconds (ns), Make it in 'ns' "<<G4endl; 
490           exit(0);                                490           exit(0);
491         }                                         491         }
492         block_DeadTime = block_DeadTime*ns;       492         block_DeadTime = block_DeadTime*ns;
493         G4cout<<"Dead time of the detector: "<    493         G4cout<<"Dead time of the detector: "<<block_DeadTime <<" ns."<<G4endl;
494       }                                           494       }
495       if( (std::string::size_type)inputLine.fi    495       if( (std::string::size_type)inputLine.find("module_DeadTime:")!=std::string::npos){
496         std::istringstream tmpStream(inputLine    496         std::istringstream tmpStream(inputLine);
497         tmpStream >> value[0] >> value[1] >> v    497         tmpStream >> value[0] >> value[1] >> value[2];
498         module_DeadTime = atof(value[1].c_str(    498         module_DeadTime = atof(value[1].c_str());
499         if(value[2] != "ns"){                     499         if(value[2] != "ns"){
500           G4cerr<<" Dead time unit is not in n    500           G4cerr<<" Dead time unit is not in nano seconds (ns), Make it in 'ns' "<<G4endl; 
501           exit(0);                                501           exit(0);
502         }                                         502         }
503         module_DeadTime = module_DeadTime*ns;     503         module_DeadTime = module_DeadTime*ns;
504         G4cout<<"Dead time of the module (axia    504         G4cout<<"Dead time of the module (axially multiplexed detectors): "<<module_DeadTime <<" ns."<<G4endl;
505       }                                           505       }
506       //                                          506       //
507       if( (std::string::size_type)inputLine.fi    507       if( (std::string::size_type)inputLine.find("crystalResolutionMin:")!=std::string::npos){
508         std::istringstream tmpStream(inputLine    508         std::istringstream tmpStream(inputLine);
509         tmpStream >> value[0] >> value[1];        509         tmpStream >> value[0] >> value[1];
510         crystalResolutionMin = atof(value[1].c    510         crystalResolutionMin = atof(value[1].c_str());
511         G4cout<<"crystal Resolution (Min.): "<    511         G4cout<<"crystal Resolution (Min.): "<<crystalResolutionMin*100<< " %." <<G4endl;
512       }                                           512       }
513       if( (std::string::size_type)inputLine.fi    513       if( (std::string::size_type)inputLine.find("crystalResolutionMax:")!=std::string::npos){
514         std::istringstream tmpStream(inputLine    514         std::istringstream tmpStream(inputLine);
515         tmpStream >> value[0] >> value[1];        515         tmpStream >> value[0] >> value[1];
516         crystalResolutionMax = atof(value[1].c    516         crystalResolutionMax = atof(value[1].c_str());
517         G4cout<<"crystal Resolution (Max.): "<    517         G4cout<<"crystal Resolution (Max.): "<<crystalResolutionMax*100<<" %"<<G4endl;
518       }                                           518       }
519                                                   519 
520       //                                          520       //
521       if( (std::string::size_type)inputLine.fi    521       if( (std::string::size_type)inputLine.find("fixedResolution:")!=std::string::npos){
522         std::istringstream tmpStream(inputLine    522         std::istringstream tmpStream(inputLine);
523         tmpStream >> value[0] >> value[1];        523         tmpStream >> value[0] >> value[1];
524         if(value[1]=="true"){                     524         if(value[1]=="true"){
525           fixedResolution = true;                 525           fixedResolution = true;
526           energyResolution_fixed = (crystalRes    526           energyResolution_fixed = (crystalResolutionMin + crystalResolutionMax)*0.5;
527           G4cout<<"Fixed crystal resolution is    527           G4cout<<"Fixed crystal resolution is used. "<<G4endl;
528         }                                         528         }
529         else {                                    529         else {
530           fixedResolution = false;                530           fixedResolution = false;
531           //Store into a file if needed.          531           //Store into a file if needed.
532           std::string fname = "crystalDependen    532           std::string fname = "crystalDependentResolution.txt";
533           std::ofstream outFname(fname.c_str()    533           std::ofstream outFname(fname.c_str());
534                                                   534 
535           G4cout<<" \n Crystal dependent resol    535           G4cout<<" \n Crystal dependent resolution is used. preparing look-up table .... "<<G4endl;
536           energyResolution_cryDependent.resize    536           energyResolution_cryDependent.resize(numberOfBlocks_total,std::vector<G4double>(numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI,0)); 
537           for(G4int i_blk = 0; i_blk < numberO    537           for(G4int i_blk = 0; i_blk < numberOfBlocks_total; i_blk++){
538             for(G4int i_cry = 0; i_cry < numbe    538             for(G4int i_cry = 0; i_cry < numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI; i_cry++){
539               energyResolution_cryDependent[i_    539               energyResolution_cryDependent[i_blk][i_cry] = crystalResolutionMin + (crystalResolutionMax - crystalResolutionMin)*G4UniformRand();
540               //store into a file                 540               //store into a file
541               outFname<<i_blk<<" "<<i_cry<<" "    541               outFname<<i_blk<<" "<<i_cry<<" "<<energyResolution_cryDependent[i_blk][i_cry]<<G4endl;
542             }                                     542             }
543           }                                       543           }
544           G4cout<<"Done. \n"<<G4endl;             544           G4cout<<"Done. \n"<<G4endl;
545           outFname.close();                       545           outFname.close();
546         }                                         546         }
547                                                   547 
548       }                                           548       }
549       //                                          549       //
550                                                   550 
551       if( (std::string::size_type)inputLine.fi    551       if( (std::string::size_type)inputLine.find("crystalEnergyRef:")!=std::string::npos){
552         std::istringstream tmpStream(inputLine    552         std::istringstream tmpStream(inputLine);
553         tmpStream >> value[0] >> value[1] >> v    553         tmpStream >> value[0] >> value[1] >> value[2];
554         crystalEnergyRef = atof(value[1].c_str    554         crystalEnergyRef = atof(value[1].c_str());
555         if(value[2] != "keV"){                    555         if(value[2] != "keV"){
556           G4cerr<<" The unit of reference ener    556           G4cerr<<" The unit of reference energy is not in keV, Make it in 'keV' "<<G4endl; 
557           exit(0);                                557           exit(0);
558         }                                         558         }
559         crystalEnergyRef = crystalEnergyRef*ke    559         crystalEnergyRef = crystalEnergyRef*keV;
560         G4cout<<"Energy of refernce: "<<crysta    560         G4cout<<"Energy of refernce: "<<crystalEnergyRef/keV<<" keV."<<G4endl;
561       }                                           561       }
562       if( (std::string::size_type)inputLine.fi    562       if( (std::string::size_type)inputLine.find("crystalQuantumEfficiency:")!=std::string::npos){
563         std::istringstream tmpStream(inputLine    563         std::istringstream tmpStream(inputLine);
564         tmpStream >> value[0] >> value[1];        564         tmpStream >> value[0] >> value[1];
565         crystalQuantumEfficiency = atof(value[    565         crystalQuantumEfficiency = atof(value[1].c_str());
566         G4cout<<"Quantum Efficiency "<<crystal    566         G4cout<<"Quantum Efficiency "<<crystalQuantumEfficiency*100<< " % "<<G4endl;
567       }                                           567       }
568       if( (std::string::size_type)inputLine.fi    568       if( (std::string::size_type)inputLine.find("lowerThreshold:")!=std::string::npos){
569         std::istringstream tmpStream(inputLine    569         std::istringstream tmpStream(inputLine);
570         tmpStream >> value[0] >> value[1] >> v    570         tmpStream >> value[0] >> value[1] >> value[2];
571         lowerThreshold = atof(value[1].c_str()    571         lowerThreshold = atof(value[1].c_str());
572         if(value[2] != "keV"){                    572         if(value[2] != "keV"){
573           G4cerr<<" The unit of Lower energy t    573           G4cerr<<" The unit of Lower energy threshold is not in keV, Make it in 'keV' "<<G4endl; 
574           exit(0);                                574           exit(0);
575         }                                         575         }
576         lowerThreshold = lowerThreshold*keV;      576         lowerThreshold = lowerThreshold*keV;
577         G4cout<<"Lower energy threshold: "<<lo    577         G4cout<<"Lower energy threshold: "<<lowerThreshold/keV<<" keV."<<G4endl;
578                                                   578 
579       }                                           579       }
580       if( (std::string::size_type)inputLine.fi    580       if( (std::string::size_type)inputLine.find("upperThreshold:")!=std::string::npos){
581         std::istringstream tmpStream(inputLine    581         std::istringstream tmpStream(inputLine);
582         tmpStream >> value[0] >> value[1] >> v    582         tmpStream >> value[0] >> value[1] >> value[2];
583         upperThreshold = atof(value[1].c_str()    583         upperThreshold = atof(value[1].c_str());
584         if(value[2] != "keV"){                    584         if(value[2] != "keV"){
585           G4cerr<<" The unit of Upper energy t    585           G4cerr<<" The unit of Upper energy threshold is not in keV, Make it in 'keV' "<<G4endl; 
586           exit(0);                                586           exit(0);
587         }                                         587         }
588         upperThreshold = upperThreshold*keV;      588         upperThreshold = upperThreshold*keV;
589         G4cout<<"Upper energy threshold: "<<up    589         G4cout<<"Upper energy threshold: "<<upperThreshold/keV<<" keV."<<G4endl;
590                                                   590 
591       }                                           591       }
592                                                   592 
593       //                                          593       //
594       if( (std::string::size_type)inputLine.fi    594       if( (std::string::size_type)inputLine.find("triggerEnergy:")!=std::string::npos){
595         std::istringstream tmpStream(inputLine    595         std::istringstream tmpStream(inputLine);
596         tmpStream >> value[0] >> value[1] >> v    596         tmpStream >> value[0] >> value[1] >> value[2];
597         triggerEnergy = atof(value[1].c_str())    597         triggerEnergy = atof(value[1].c_str());
598         if(value[2] != "keV"){                    598         if(value[2] != "keV"){
599           G4cerr<<" The unit of Trigger energy    599           G4cerr<<" The unit of Trigger energy threshold is not in keV, Make it in 'keV' "<<G4endl; 
600           exit(0);                                600           exit(0);
601         }                                         601         }
602         triggerEnergy = triggerEnergy*keV;        602         triggerEnergy = triggerEnergy*keV;
603         G4cout<<"Trigger energy threshold: "<<    603         G4cout<<"Trigger energy threshold: "<<triggerEnergy/keV<<" keV."<<G4endl;
604                                                   604 
605       }                                           605       }
606                                                   606 
607       //Option to apply AngerLogic                607       //Option to apply AngerLogic
608       if( (std::string::size_type)inputLine.fi    608       if( (std::string::size_type)inputLine.find("ApplyAngerLogic:")!=std::string::npos){
609         std::istringstream tmpStream(inputLine    609         std::istringstream tmpStream(inputLine);
610         tmpStream >> value[0] >> value[1];        610         tmpStream >> value[0] >> value[1];
611         if(value[1]=="true"){                     611         if(value[1]=="true"){
612           ApplyAngerLogic = true;                 612           ApplyAngerLogic = true;
613           G4cout<<"Angler Logic calculation is    613           G4cout<<"Angler Logic calculation is applied. "<<G4endl;
614         }                                         614         }
615         else if(value[1]=="false") {              615         else if(value[1]=="false") {
616           ApplyAngerLogic = false;                616           ApplyAngerLogic = false;
617           G4cout<<"Angler Logic calculation is    617           G4cout<<"Angler Logic calculation is NOT applied. "<<G4endl;
618         }                                         618         }
619         else {                                    619         else {
620           ApplyAngerLogic = true;                 620           ApplyAngerLogic = true;
621           G4cout<<"Angler Logic calculation is    621           G4cout<<"Angler Logic calculation is applied (by defualt). "<<G4endl;
622         }                                         622         }
623       }                                           623       }
624                                                   624 
625       //PMT position calculation blurring at F    625       //PMT position calculation blurring at FWHM
626       if( (std::string::size_type)inputLine.fi    626       if( (std::string::size_type)inputLine.find("PMTblurring:")!=std::string::npos){
627         std::istringstream tmpStream(inputLine    627         std::istringstream tmpStream(inputLine);
628         tmpStream >> value[0] >> value[1]>>val    628         tmpStream >> value[0] >> value[1]>>value[2];
629         PMTblurring_tan = atof(value[1].c_str(    629         PMTblurring_tan = atof(value[1].c_str());
630         PMTblurring_axial = atof(value[2].c_st    630         PMTblurring_axial = atof(value[2].c_str());
631         G4cout<<"PMTblurring position response    631         G4cout<<"PMTblurring position response blurring at FWHM (tan x axial) "<<PMTblurring_tan<<" x " <<PMTblurring_axial<<" mm2"<<G4endl;
632       }                                           632       }
633                                                   633 
634       //                                          634       //
635       if( (std::string::size_type)inputLine.fi    635       if( (std::string::size_type)inputLine.find("TypeOfOutput:")!=std::string::npos){
636         std::istringstream tmpStream(inputLine    636         std::istringstream tmpStream(inputLine);
637         tmpStream >> value[0] >> value[1];        637         tmpStream >> value[0] >> value[1];
638         if(value[1]=="singlesOutput"){            638         if(value[1]=="singlesOutput"){
639           getSinglesData = true;                  639           getSinglesData = true;
640           G4cout<<"Single mode output enabled.    640           G4cout<<"Single mode output enabled. "<<G4endl;
641         }                                         641         }
642         else if(value[1]=="coincidenceOutput")    642         else if(value[1]=="coincidenceOutput") {
643           getCoincidenceData = true;              643           getCoincidenceData = true;
644           G4cout<<"Coicidence mode output enab    644           G4cout<<"Coicidence mode output enabled. "<<G4endl;
645         }                                         645         }
646                                                   646 
647       }                                           647       }
648                                                   648 
649     }                                             649     }
650   }                                               650   }
651   ifs.close();                                    651   ifs.close();
652 }                                                 652 }
653                                                   653 
654 //The following function reads the reflector p    654 //The following function reads the reflector pattern for each layer. Each layer has different patterns along the tangetial and axial positions.
655 //For defualt reflector pattern, see https://l    655 //For defualt reflector pattern, see https://link.springer.com/article/10.1007/s12194-013-0231-4
656 //The patter of the reflectors can be changed     656 //The patter of the reflectors can be changed in the inputParameter.txt file
657 //The pattern is given as 0 and 1. If there is    657 //The pattern is given as 0 and 1. If there is reflector the value is 1 and if there is no reflector, the value is 0.
658                                                   658 
659 void doiPETAnalysis::ReadReflectorPattern(){      659 void doiPETAnalysis::ReadReflectorPattern(){
660   G4cout<<" Reflector pattern is being read "<    660   G4cout<<" Reflector pattern is being read "<<G4endl;
661   //                                              661   //
662   std::vector<std::string> stringReflectorValu    662   std::vector<std::string> stringReflectorValue;
663   //                                              663   //
664   char inputChar[256];                            664   char inputChar[256];
665   std::string inputLine;                          665   std::string inputLine;
666                                                   666 
667   //open inputParameter.txt to read reflector     667   //open inputParameter.txt to read reflector pattern. 
668   std::string  filename = "inputParameter.txt"    668   std::string  filename = "inputParameter.txt";
669                                                   669 
670   G4String refValue;                              670   G4String refValue;
671                                                   671 
672   ifs.open(filename.c_str());                     672   ifs.open(filename.c_str());
673   if(!ifs.good()){                                673   if(!ifs.good()){
674     G4cerr<<"File opening Error: Could not ope    674     G4cerr<<"File opening Error: Could not open "<<filename<<G4endl;
675     exit(0);                                      675     exit(0);
676   }                                               676   }
677   while(!ifs.eof()){                              677   while(!ifs.eof()){
678     ifs.getline(inputChar,256);                   678     ifs.getline(inputChar,256);
679     inputLine = inputChar;                        679     inputLine = inputChar;
680                                                   680 
681     //The reflector patter in read from the in    681     //The reflector patter in read from the inputparamter.txt file      
682     if(inputChar[0]!='#' && inputLine.length()    682     if(inputChar[0]!='#' && inputLine.length()!=0 ){
683                                                   683 
684       //Reflector patter for Layer1 in the tan    684       //Reflector patter for Layer1 in the tangential direction
685       if( (std::string::size_type)inputLine.fi    685       if( (std::string::size_type)inputLine.find("reflectorLayer1_Tangential:")!=std::string::npos){
686         std::istringstream tmpStream(inputLine    686         std::istringstream tmpStream(inputLine);
687         while(tmpStream >> refValue){             687         while(tmpStream >> refValue){
688           stringReflectorValue.push_back(refVa    688           stringReflectorValue.push_back(refValue);
689           if(stringReflectorValue.size()>1){      689           if(stringReflectorValue.size()>1){
690             G4int tmp_value = atoi(stringRefle    690             G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
691             ireflectorLayer1_Tangential.push_b    691             ireflectorLayer1_Tangential.push_back(tmp_value);
692           }                                       692           }
693         }                                         693         }
694         stringReflectorValue.clear();             694         stringReflectorValue.clear();
695       }                                           695       }
696                                                   696       
697                                                   697 
698       //Reflector patter for Layer1 in the axi    698       //Reflector patter for Layer1 in the axial direction
699       if( (std::string::size_type)inputLine.fi    699       if( (std::string::size_type)inputLine.find("reflectorLayer1_Axial:")!=std::string::npos){
700         std::istringstream tmpStream(inputLine    700         std::istringstream tmpStream(inputLine);
701         while(tmpStream >> refValue){             701         while(tmpStream >> refValue){
702           stringReflectorValue.push_back(refVa    702           stringReflectorValue.push_back(refValue);
703           if(stringReflectorValue.size()>1){      703           if(stringReflectorValue.size()>1){
704             G4int tmp_value = atoi(stringRefle    704             G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
705             ireflectorLayer1_Axial.push_back(t    705             ireflectorLayer1_Axial.push_back(tmp_value);
706           }                                       706           }
707         }                                         707         }
708         stringReflectorValue.clear();             708         stringReflectorValue.clear();
709       }                                           709       }
710                                                   710 
711       //Reflector patter for Layer2 in the tan    711       //Reflector patter for Layer2 in the tangential direction
712       if( (std::string::size_type)inputLine.fi    712       if( (std::string::size_type)inputLine.find("reflectorLayer2_Tangential:")!=std::string::npos){
713         std::istringstream tmpStream(inputLine    713         std::istringstream tmpStream(inputLine);
714         while(tmpStream >> refValue){             714         while(tmpStream >> refValue){
715           stringReflectorValue.push_back(refVa    715           stringReflectorValue.push_back(refValue);
716           if(stringReflectorValue.size()>1){      716           if(stringReflectorValue.size()>1){
717             G4int tmp_value = atoi(stringRefle    717             G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
718             ireflectorLayer2_Tangential.push_b    718             ireflectorLayer2_Tangential.push_back(tmp_value);
719           }                                       719           }
720         }                                         720         }
721         stringReflectorValue.clear();             721         stringReflectorValue.clear();
722       }                                           722       }
723                                                   723 
724       //Reflector patter for Layer2 in the axi    724       //Reflector patter for Layer2 in the axial direction
725       if( (std::string::size_type)inputLine.fi    725       if( (std::string::size_type)inputLine.find("reflectorLayer2_Axial:")!=std::string::npos){
726         std::istringstream tmpStream(inputLine    726         std::istringstream tmpStream(inputLine);
727         while(tmpStream >> refValue){             727         while(tmpStream >> refValue){
728           stringReflectorValue.push_back(refVa    728           stringReflectorValue.push_back(refValue);
729           if(stringReflectorValue.size()>1){      729           if(stringReflectorValue.size()>1){
730             G4int tmp_value = atoi(stringRefle    730             G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
731             ireflectorLayer2_Axial.push_back(t    731             ireflectorLayer2_Axial.push_back(tmp_value);
732           }                                       732           }
733         }                                         733         }
734         stringReflectorValue.clear();             734         stringReflectorValue.clear();
735       }                                           735       }
736                                                   736 
737       //Reflector patter for Layer3 in the tan    737       //Reflector patter for Layer3 in the tangential direction
738       if( (std::string::size_type)inputLine.fi    738       if( (std::string::size_type)inputLine.find("reflectorLayer3_Tangential:")!=std::string::npos){
739         std::istringstream tmpStream(inputLine    739         std::istringstream tmpStream(inputLine);
740         while(tmpStream >> refValue){             740         while(tmpStream >> refValue){
741           stringReflectorValue.push_back(refVa    741           stringReflectorValue.push_back(refValue);
742           if(stringReflectorValue.size()>1){      742           if(stringReflectorValue.size()>1){
743             G4int tmp_value = atoi(stringRefle    743             G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
744             ireflectorLayer3_Tangential.push_b    744             ireflectorLayer3_Tangential.push_back(tmp_value);
745           }                                       745           }
746         }                                         746         }
747         stringReflectorValue.clear();             747         stringReflectorValue.clear();
748       }                                           748       }
749                                                   749 
750       //Reflector patter for Layer3 in the axi    750       //Reflector patter for Layer3 in the axial direction
751       if( (std::string::size_type)inputLine.fi    751       if( (std::string::size_type)inputLine.find("reflectorLayer3_Axial:")!=std::string::npos){
752         std::istringstream tmpStream(inputLine    752         std::istringstream tmpStream(inputLine);
753         while(tmpStream >> refValue){             753         while(tmpStream >> refValue){
754           stringReflectorValue.push_back(refVa    754           stringReflectorValue.push_back(refValue);
755           if(stringReflectorValue.size()>1){      755           if(stringReflectorValue.size()>1){
756             G4int tmp_value = atoi(stringRefle    756             G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
757             ireflectorLayer3_Axial.push_back(t    757             ireflectorLayer3_Axial.push_back(tmp_value);
758           }                                       758           }
759         }                                         759         }
760         stringReflectorValue.clear();             760         stringReflectorValue.clear();
761       }                                           761       }
762                                                   762 
763       //Reflector patter for Layer4 in the tan    763       //Reflector patter for Layer4 in the tangential direction
764       if( (std::string::size_type)inputLine.fi    764       if( (std::string::size_type)inputLine.find("reflectorLayer4_Tangential:")!=std::string::npos){
765         std::istringstream tmpStream(inputLine    765         std::istringstream tmpStream(inputLine);
766         while(tmpStream >> refValue){             766         while(tmpStream >> refValue){
767           stringReflectorValue.push_back(refVa    767           stringReflectorValue.push_back(refValue);
768           if(stringReflectorValue.size()>1){      768           if(stringReflectorValue.size()>1){
769             G4int tmp_value = atoi(stringRefle    769             G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
770             ireflectorLayer4_Tangential.push_b    770             ireflectorLayer4_Tangential.push_back(tmp_value);
771           }                                       771           }
772         }                                         772         }
773         stringReflectorValue.clear();             773         stringReflectorValue.clear();
774       }                                           774       }
775                                                   775 
776       //Reflector patter for Layer4 in the axi    776       //Reflector patter for Layer4 in the axial direction
777       if( (std::string::size_type)inputLine.fi    777       if( (std::string::size_type)inputLine.find("reflectorLayer4_Axial:")!=std::string::npos){
778         std::istringstream tmpStream(inputLine    778         std::istringstream tmpStream(inputLine);
779         while(tmpStream >> refValue){             779         while(tmpStream >> refValue){
780           stringReflectorValue.push_back(refVa    780           stringReflectorValue.push_back(refValue);
781           if(stringReflectorValue.size()>1){      781           if(stringReflectorValue.size()>1){
782             G4int tmp_value = atoi(stringRefle    782             G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
783             ireflectorLayer4_Axial.push_back(t    783             ireflectorLayer4_Axial.push_back(tmp_value);
784           }                                       784           }
785         }                                         785         }
786         stringReflectorValue.clear();             786         stringReflectorValue.clear();
787       }                                           787       }
788     }//#                                          788     }//#
789   }//while(eof)                                   789   }//while(eof)
790                                                   790 
791   //                                              791   //
792   //for debug                                     792   //for debug
793   G4cout<<"\n========= Reflector Pattern =====    793   G4cout<<"\n========= Reflector Pattern ==========="<<G4endl;
794   G4cout<<"Layer 1"<<G4endl;                      794   G4cout<<"Layer 1"<<G4endl;
795   for(unsigned int i = 0; i<ireflectorLayer1_T    795   for(unsigned int i = 0; i<ireflectorLayer1_Tangential.size();i++){
796     G4cout<<ireflectorLayer1_Tangential[i]<<"     796     G4cout<<ireflectorLayer1_Tangential[i]<<" ";
797   }G4cout<<G4endl;                                797   }G4cout<<G4endl;
798   for(unsigned int i = 0; i<ireflectorLayer1_A    798   for(unsigned int i = 0; i<ireflectorLayer1_Axial.size();i++){
799     G4cout<<ireflectorLayer1_Axial[i]<<" ";       799     G4cout<<ireflectorLayer1_Axial[i]<<" ";
800   }G4cout<<G4endl;                                800   }G4cout<<G4endl;
801   G4cout<<"Layer 2"<<G4endl;                      801   G4cout<<"Layer 2"<<G4endl;
802   for(unsigned int i = 0; i<ireflectorLayer2_T    802   for(unsigned int i = 0; i<ireflectorLayer2_Tangential.size();i++){
803     G4cout<<ireflectorLayer2_Tangential[i]<<"     803     G4cout<<ireflectorLayer2_Tangential[i]<<" ";
804   }G4cout<<G4endl;                                804   }G4cout<<G4endl;
805   for(unsigned int i = 0; i<ireflectorLayer2_A    805   for(unsigned int i = 0; i<ireflectorLayer2_Axial.size();i++){
806     G4cout<<ireflectorLayer2_Axial[i]<<" ";       806     G4cout<<ireflectorLayer2_Axial[i]<<" ";
807   }G4cout<<G4endl;                                807   }G4cout<<G4endl;
808   G4cout<<"Layer 3"<<G4endl;                      808   G4cout<<"Layer 3"<<G4endl;
809   for(unsigned int i = 0; i<ireflectorLayer3_T    809   for(unsigned int i = 0; i<ireflectorLayer3_Tangential.size();i++){
810     G4cout<<ireflectorLayer3_Tangential[i]<<"     810     G4cout<<ireflectorLayer3_Tangential[i]<<" ";
811   }G4cout<<G4endl;                                811   }G4cout<<G4endl;
812   for(unsigned int i = 0; i<ireflectorLayer3_A    812   for(unsigned int i = 0; i<ireflectorLayer3_Axial.size();i++){
813     G4cout<<ireflectorLayer3_Axial[i]<<" ";       813     G4cout<<ireflectorLayer3_Axial[i]<<" ";
814   }G4cout<<G4endl;                                814   }G4cout<<G4endl;
815   G4cout<<"Layer 4"<<G4endl;                      815   G4cout<<"Layer 4"<<G4endl;
816   for(unsigned int i = 0; i<ireflectorLayer4_T    816   for(unsigned int i = 0; i<ireflectorLayer4_Tangential.size();i++){
817     G4cout<<ireflectorLayer4_Tangential[i]<<"     817     G4cout<<ireflectorLayer4_Tangential[i]<<" ";
818   }G4cout<<G4endl;                                818   }G4cout<<G4endl;
819   for(unsigned int i = 0; i<ireflectorLayer4_A    819   for(unsigned int i = 0; i<ireflectorLayer4_Axial.size();i++){
820     G4cout<<ireflectorLayer4_Axial[i]<<" ";       820     G4cout<<ireflectorLayer4_Axial[i]<<" ";
821   }G4cout<<G4endl;                                821   }G4cout<<G4endl;
822   G4cout<<"========= Reflector Pattern Ended =    822   G4cout<<"========= Reflector Pattern Ended ===========\n"<<G4endl;
823                                                   823   
824   //Read DOI look-up-table. This look-up-table    824   //Read DOI look-up-table. This look-up-table is prepared based on the assumption that the interaction is occured at the center of the crystal.
825   G4int index_doi = 0, doiID;                     825   G4int index_doi = 0, doiID;
826   doi_table.resize(numberOfCrystal_tangential*    826   doi_table.resize(numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI,0);
827                                                   827 
828   std::string LUT_FileName = "look_up_table_DO    828   std::string LUT_FileName = "look_up_table_DOI.txt";
829   std::ifstream ifs_doiLUT;                       829   std::ifstream ifs_doiLUT;
830   std::ofstream ofs_doiLUT;                       830   std::ofstream ofs_doiLUT;
831   ifs_doiLUT.open(LUT_FileName.c_str());          831   ifs_doiLUT.open(LUT_FileName.c_str());
832   if(ifs_doiLUT.is_open()){                       832   if(ifs_doiLUT.is_open()){
833     G4cout<<" DOI Look-up table found and used    833     G4cout<<" DOI Look-up table found and used: File name: "<<LUT_FileName<<G4endl;
834     //Read from file                              834     //Read from file
835     while(ifs_doiLUT>>index_doi>>doiID && inde    835     while(ifs_doiLUT>>index_doi>>doiID && index_doi < int(doi_table.size())){
836       doi_table[index_doi] = doiID;               836       doi_table[index_doi] = doiID;
837     }                                             837     }
838     if(index_doi==int(doi_table.size())){         838     if(index_doi==int(doi_table.size())){
839       G4cout<<"!!!Warning: The DOI table index    839       G4cout<<"!!!Warning: The DOI table index is greater than the total number of crystals."<<G4endl;
840       PrepareDOILookUpTable(LUT_FileName);        840       PrepareDOILookUpTable(LUT_FileName);
841     }                                             841     }
842     isDOIlookUpTablePrepared = true; //           842     isDOIlookUpTablePrepared = true; //
843   }                                               843   }
844   else                                            844   else
845   {                                               845   {
846     PrepareDOILookUpTable(LUT_FileName);          846     PrepareDOILookUpTable(LUT_FileName);
847   }                                               847   }
848   //Write into a file.                            848   //Write into a file.
849   ofs_doiLUT.open(LUT_FileName.c_str());          849   ofs_doiLUT.open(LUT_FileName.c_str());
850   if(!ofs_doiLUT.is_open()){                      850   if(!ofs_doiLUT.is_open()){
851     G4cerr<<"Unable to open file to write doi_    851     G4cerr<<"Unable to open file to write doi_LUT"<<G4endl;
852     exit(0);                                      852     exit(0);
853   }                                               853   }
854   for(G4int i=0;i<int(doi_table.size()); i++){    854   for(G4int i=0;i<int(doi_table.size()); i++){
855     ofs_doiLUT<<i<<"\t"<<doi_table[i]<<G4endl;    855     ofs_doiLUT<<i<<"\t"<<doi_table[i]<<G4endl;
856   }                                               856   }
857   ifs_doiLUT.close();                             857   ifs_doiLUT.close();
858                                                   858 
859                                                   859 
860   ifs.close();                                    860   ifs.close();
861 }                                                 861 }
862 void doiPETAnalysis::PrepareDOILookUpTable(G4S    862 void doiPETAnalysis::PrepareDOILookUpTable(G4String){
863   isDOIlookUpTablePrepared = false;               863   isDOIlookUpTablePrepared = false; 
864     G4cout<<"Preparing DOI look-up table... "<    864     G4cout<<"Preparing DOI look-up table... "<<G4endl;
865     std::string outputFileName = "_check_2Dpos    865     std::string outputFileName = "_check_2Dposition.txt";// excuted only once
866     std::ofstream outFile(outputFileName.c_str    866     std::ofstream outFile(outputFileName.c_str());
867                                                   867 
868     G4double crystalPositionX;                    868     G4double crystalPositionX;
869     G4double crystalPositionY;                    869     G4double crystalPositionY;
870     G4double crystalPositionZ;                    870     G4double crystalPositionZ;
871     //doi_table.resize(numberOfCrystal_tangent    871     //doi_table.resize(numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI,0);
872                                                   872 
873     for(G4int i_DOI = 0; i_DOI<numberOfCrystal    873     for(G4int i_DOI = 0; i_DOI<numberOfCrystal_DOI; i_DOI++){
874       crystalPositionX=(i_DOI-((float)numberOf    874       crystalPositionX=(i_DOI-((float)numberOfCrystal_DOI)/2 + 0.5)*(sizeOfCrystal_DOI + crystalGap_DOI); //Becuase only lateral distances are used
875       for(G4int i_axial=0; i_axial< numberOfCr    875       for(G4int i_axial=0; i_axial< numberOfCrystal_axial;i_axial++){
876         crystalPositionZ = (i_axial-((float)nu    876         crystalPositionZ = (i_axial-((float)numberOfCrystal_axial)/2 + 0.5)*(sizeOfCrystal_axial + crystalGap_axial);
877         for(G4int i_tan=0; i_tan<numberOfCryst    877         for(G4int i_tan=0; i_tan<numberOfCrystal_tangential;i_tan++){
878           crystalPositionY=(i_tan-((float)numb    878           crystalPositionY=(i_tan-((float)numberOfCrystal_tangential)/2 + 0.5)*(sizeOfCrystal_tangential + crystalGap_tangential);
879           AngerLogic(crystalPositionX, crystal    879           AngerLogic(crystalPositionX, crystalPositionY, crystalPositionZ, 1, 0.5, isDOIlookUpTablePrepared);
880           outFile<<PositionAngerZ<<" "<<Positi    880           outFile<<PositionAngerZ<<" "<<PositionAngerY<<G4endl;
881           doi_table[crystalID_in2D_posHist]=i_    881           doi_table[crystalID_in2D_posHist]=i_DOI;
882                                                   882 
883         }                                         883         }
884       }                                           884       }
885     }                                             885     }
886     G4cout<<"done."<<G4endl;                      886     G4cout<<"done."<<G4endl;
887     isDOIlookUpTablePrepared = true;              887     isDOIlookUpTablePrepared = true;
888 }                                                 888 }
889                                                   889 
890                                                   890 
891 //Based on ideal photomultiplier tube (PMT) pl    891 //Based on ideal photomultiplier tube (PMT) placement, the interaction position of the photon with the detector is calculated using Anger Logic method. 
892 //The reflectors shifts the response by some d    892 //The reflectors shifts the response by some distance so that the response can be projected into 2D position histogram. 
893 //From this 2D position histogram, the new cry    893 //From this 2D position histogram, the new crystal ID (in 3D along the tangential (y), axial (z) and DOI (x)) (after Anger Logic method is applied) can be obtained.
894 //If the crystal ID after Anger method apllied    894 //If the crystal ID after Anger method apllied is out of the give number of crystals (in 3D), then an error message is displayed and the event will be rejected.
895 void doiPETAnalysis::AngerLogic(G4double posCr    895 void doiPETAnalysis::AngerLogic(G4double posCrystalX, G4double posCrystalY, G4double posCrystalZ, G4double Edep, G4double shiftDis, G4bool isDOI_LUT)
896 {                                                 896 { 
897   G4double crystalPitch_DOI = sizeOfCrystal_DO    897   G4double crystalPitch_DOI = sizeOfCrystal_DOI + crystalGap_DOI;
898   G4double crystalPitch_tan = sizeOfCrystal_ta    898   G4double crystalPitch_tan = sizeOfCrystal_tangential + crystalGap_tangential;
899   G4double crystalPitch_axial = sizeOfCrystal_    899   G4double crystalPitch_axial = sizeOfCrystal_axial + crystalGap_axial;
900                                                   900 
901   //The crystal ID are calculated based on the    901   //The crystal ID are calculated based on the center of mass 
902   G4int i_doi = posCrystalX/crystalPitch_DOI +    902   G4int i_doi = posCrystalX/crystalPitch_DOI + (float)numberOfCrystal_DOI*0.5;
903   G4int i_tan = posCrystalY/crystalPitch_tan +    903   G4int i_tan = posCrystalY/crystalPitch_tan + (float)numberOfCrystal_tangential*0.5;
904   G4int i_axial = posCrystalZ/crystalPitch_axi    904   G4int i_axial = posCrystalZ/crystalPitch_axial + (float)numberOfCrystal_axial*0.5;
905                                                   905 
906   //position of interaction is shifted  the ce    906   //position of interaction is shifted  the centre of the crystal. This is to use DOI-look up tables as the real scanner
907   posCrystalX = (i_doi-((float)numberOfCrystal    907   posCrystalX = (i_doi-((float)numberOfCrystal_DOI)/2 + 0.5)*crystalPitch_DOI;
908   posCrystalY = (i_tan-((float)numberOfCrystal    908   posCrystalY = (i_tan-((float)numberOfCrystal_tangential)/2 + 0.5)*crystalPitch_tan;
909   posCrystalZ = (i_axial-((float)numberOfCryst    909   posCrystalZ = (i_axial-((float)numberOfCrystal_axial)/2 + 0.5)*crystalPitch_axial;
910                                                   910 
911   //1z and 2z are at the same z distance; 3z a    911   //1z and 2z are at the same z distance; 3z and 4z are at the same z distance
912   //The signal (the energy deposition) is devi    912   //The signal (the energy deposition) is devided into the four PMTs depending on their lateral distances (in the axial and tangential directions) from the interaction position
913                                                   913 
914   //The following is based on symetrical placm    914   //The following is based on symetrical placment of the PMTs 
915                                                   915 
916   //Calculate the axial (z) distance from the     916   //Calculate the axial (z) distance from the position of interaction to each PMT
917   dist1z = std::fabs(posPMT1z - posCrystalZ);     917   dist1z = std::fabs(posPMT1z - posCrystalZ);
918   dist2z = std::fabs(posPMT2z - posCrystalZ);     918   dist2z = std::fabs(posPMT2z - posCrystalZ);
919   dist3z = std::fabs(posPMT3z - posCrystalZ);     919   dist3z = std::fabs(posPMT3z - posCrystalZ);
920   dist4z = std::fabs(posPMT4z - posCrystalZ);     920   dist4z = std::fabs(posPMT4z - posCrystalZ);
921                                                   921 
922   //Calculate the resultant distance              922   //Calculate the resultant distance
923   //dist1z = dist2z, and dist3z = dist4z, so o    923   //dist1z = dist2z, and dist3z = dist4z, so only take two of them or take the average
924   //distz = ((dist1z + dist2z) + (dist3z + dis    924   //distz = ((dist1z + dist2z) + (dist3z + dist4z))/2; 
925   distz = dist1z + dist3z;                        925   distz = dist1z + dist3z; 
926                                                   926 
927   //1y and 3y are at the same y distance; and     927   //1y and 3y are at the same y distance; and 2y and 4y are at the same y distance
928   //Calculate the tangential (y)  distance fro    928   //Calculate the tangential (y)  distance from the position of interaction to each PMT
929   dist1y = std::fabs(posPMT1y - posCrystalY);     929   dist1y = std::fabs(posPMT1y - posCrystalY);
930   dist2y = std::fabs(posPMT2y - posCrystalY);     930   dist2y = std::fabs(posPMT2y - posCrystalY);
931   dist3y = std::fabs(posPMT3y - posCrystalY);     931   dist3y = std::fabs(posPMT3y - posCrystalY);
932   dist4y = std::fabs(posPMT4y - posCrystalY);     932   dist4y = std::fabs(posPMT4y - posCrystalY);
933                                                   933 
934   //Calculate the resultant distance              934   //Calculate the resultant distance
935   //dist1y = dist3y, and dist2y = dist4y, so o    935   //dist1y = dist3y, and dist2y = dist4y, so only take two of them or take the average
936   //disty = ((dist1y + dist3y) + (dist2y+dist4    936   //disty = ((dist1y + dist3y) + (dist2y+dist4y))/2;
937   disty = dist1y + dist2y;                        937   disty = dist1y + dist2y;
938                                                   938 
939   //signalPMT1z = signalPMT2z, and signalPMT3z    939   //signalPMT1z = signalPMT2z, and signalPMT3z = signalPMT4z
940   signalPMT1z = Edep * dist3z/(dist1z + dist3z    940   signalPMT1z = Edep * dist3z/(dist1z + dist3z);
941   signalPMT3z = Edep * dist1z/(dist1z + dist3z    941   signalPMT3z = Edep * dist1z/(dist1z + dist3z);
942                                                   942 
943   signalPMT2z = Edep * dist4z/(dist2z + dist4z    943   signalPMT2z = Edep * dist4z/(dist2z + dist4z);
944   signalPMT4z = Edep * dist2z/(dist2z + dist4z    944   signalPMT4z = Edep * dist2z/(dist2z + dist4z);
945                                                   945 
946                                                   946 
947   //signalPMT1y = signalPMT3y, and signalPMT2y    947   //signalPMT1y = signalPMT3y, and signalPMT2y = signalPMT4y
948   signalPMT1y = Edep * dist2y/(dist1y + dist2y    948   signalPMT1y = Edep * dist2y/(dist1y + dist2y);
949   signalPMT2y = Edep * dist1y/(dist1y + dist2y    949   signalPMT2y = Edep * dist1y/(dist1y + dist2y);
950                                                   950 
951   signalPMT3y = Edep * dist4y/(dist3y + dist4y    951   signalPMT3y = Edep * dist4y/(dist3y + dist4y);
952   signalPMT4y = Edep * dist3y/(dist3y + dist4y    952   signalPMT4y = Edep * dist3y/(dist3y + dist4y);
953                                                   953 
954   //Calculate the signal on each PMT from the     954   //Calculate the signal on each PMT from the 'component' signal
955   signalPMT1 = (signalPMT1z +  signalPMT1y)*0.    955   signalPMT1 = (signalPMT1z +  signalPMT1y)*0.5; 
956   signalPMT2 = (signalPMT2z +  signalPMT2y)*0.    956   signalPMT2 = (signalPMT2z +  signalPMT2y)*0.5;
957   signalPMT3 = (signalPMT3z +  signalPMT3y)*0.    957   signalPMT3 = (signalPMT3z +  signalPMT3y)*0.5;
958   signalPMT4 = (signalPMT4z +  signalPMT4y)*0.    958   signalPMT4 = (signalPMT4z +  signalPMT4y)*0.5;
959                                                   959 
960                                                   960 
961   signalZplus = (signalPMT3 + signalPMT4);        961   signalZplus = (signalPMT3 + signalPMT4);
962   signalZminus = (signalPMT1 + signalPMT2);       962   signalZminus = (signalPMT1 + signalPMT2);
963   signalYplus = (signalPMT2 + signalPMT4);        963   signalYplus = (signalPMT2 + signalPMT4);
964   signalYminus = (signalPMT1 + signalPMT3);       964   signalYminus = (signalPMT1 + signalPMT3);
965                                                   965 
966                                                   966 
967   //Position of interaction is calculated base    967   //Position of interaction is calculated based on Anger logic method. 
968   //To get the position by Anger calculation,     968   //To get the position by Anger calculation, the result should be multiplied by the dimenion of the total distance.
969   PositionAngerZ = (signalZplus - signalZminus    969   PositionAngerZ = (signalZplus - signalZminus)/(signalZplus + signalZminus)*distz; 
970   PositionAngerY = (signalYplus - signalYminus    970   PositionAngerY = (signalYplus - signalYminus)/(signalYplus + signalYminus)*disty;
971                                                   971 
972                                                   972 
973   //For detectors with reflector insertion (li    973   //For detectors with reflector insertion (light sharing), the response is shifted depending on the reflector patter.
974   //Here, it is assumed that the shift of the     974   //Here, it is assumed that the shift of the response is equal to half of the distance from the interaction position to the airgap in the lateral (transversal direction)
975                                                   975 
976   //If reflector is only in the left side of t    976   //If reflector is only in the left side of the crystal, then response shift is to the right side (away from the reflector).
977   //If reflector is only in the right side of     977   //If reflector is only in the right side of the crystal, then response shift is to the left side (away from the reflector).
978                                                   978 
979   //Response shift for 1st Layer                  979   //Response shift for 1st Layer
980   if(i_doi == 0){                                 980   if(i_doi == 0){
981     //If reflector is only in one (left) side     981     //If reflector is only in one (left) side of the crystal, then response shifts to the right side (away from the reflector)
982     if(ireflectorLayer1_Tangential[i_tan] == 1    982     if(ireflectorLayer1_Tangential[i_tan] == 1 && ireflectorLayer1_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2)*shiftDis;
983                                                   983 
984     //If reflector is only in one (right) side    984     //If reflector is only in one (right) side of the crystal, then response shifts to the left side (away from the reflector)
985     if(ireflectorLayer1_Tangential[i_tan] == 0    985     if(ireflectorLayer1_Tangential[i_tan] == 0 && ireflectorLayer1_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
986                                                   986 
987     if(ireflectorLayer1_Axial[i_axial] == 1 &&    987     if(ireflectorLayer1_Axial[i_axial] == 1 && ireflectorLayer1_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2)*shiftDis;
988     if(ireflectorLayer1_Axial[i_axial] == 0 &&    988     if(ireflectorLayer1_Axial[i_axial] == 0 && ireflectorLayer1_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
989   }                                               989   }
990   if(i_doi == 1){ //Response shift for 2nd Lay    990   if(i_doi == 1){ //Response shift for 2nd Layer
991     if(ireflectorLayer2_Tangential[i_tan] == 1    991     if(ireflectorLayer2_Tangential[i_tan] == 1 && ireflectorLayer2_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2)*shiftDis;
992     if(ireflectorLayer2_Tangential[i_tan] == 0    992     if(ireflectorLayer2_Tangential[i_tan] == 0 && ireflectorLayer2_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
993                                                   993 
994     if(ireflectorLayer2_Axial[i_axial] == 1 &&    994     if(ireflectorLayer2_Axial[i_axial] == 1 && ireflectorLayer2_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2)*shiftDis;
995     if(ireflectorLayer2_Axial[i_axial] == 0 &&    995     if(ireflectorLayer2_Axial[i_axial] == 0 && ireflectorLayer2_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
996   }                                               996   }
997   if(i_doi == 2){ //Response shift for 3rd Lay    997   if(i_doi == 2){ //Response shift for 3rd Layer
998     if(ireflectorLayer3_Tangential[i_tan] == 1    998     if(ireflectorLayer3_Tangential[i_tan] == 1 && ireflectorLayer3_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2)*shiftDis;
999     if(ireflectorLayer3_Tangential[i_tan] == 0    999     if(ireflectorLayer3_Tangential[i_tan] == 0 && ireflectorLayer3_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
1000                                                  1000 
1001     if(ireflectorLayer3_Axial[i_axial] == 1 &    1001     if(ireflectorLayer3_Axial[i_axial] == 1 && ireflectorLayer3_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2)*shiftDis;
1002     if(ireflectorLayer3_Axial[i_axial] == 0 &    1002     if(ireflectorLayer3_Axial[i_axial] == 0 && ireflectorLayer3_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
1003   }                                              1003   }
1004   if(i_doi == 3){ //Response shift for 4th La    1004   if(i_doi == 3){ //Response shift for 4th Layer
1005     if(ireflectorLayer4_Tangential[i_tan] ==     1005     if(ireflectorLayer4_Tangential[i_tan] == 1 && ireflectorLayer4_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2)*shiftDis;
1006     if(ireflectorLayer4_Tangential[i_tan] ==     1006     if(ireflectorLayer4_Tangential[i_tan] == 0 && ireflectorLayer4_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
1007                                                  1007 
1008     if(ireflectorLayer4_Axial[i_axial] == 1 &    1008     if(ireflectorLayer4_Axial[i_axial] == 1 && ireflectorLayer4_Axial [i_axial + 1] == 0) PositionAngerZ += (crystalPitch_axial/2)*shiftDis;
1009     if(ireflectorLayer4_Axial[i_axial] == 0 &    1009     if(ireflectorLayer4_Axial[i_axial] == 0 && ireflectorLayer4_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
1010   }                                              1010   }
1011                                                  1011    
1012    //Blur the 2D position (obtained by ANger     1012    //Blur the 2D position (obtained by ANger Logic method) to include uncertainity of the PMT position response.
1013   if(isDOI_LUT){                                 1013   if(isDOI_LUT){
1014     PositionAngerZ = G4RandGauss::shoot(Posit    1014     PositionAngerZ = G4RandGauss::shoot(PositionAngerZ,PMTblurring_axial/2.35);
1015     PositionAngerY = G4RandGauss::shoot(Posit    1015     PositionAngerY = G4RandGauss::shoot(PositionAngerY,PMTblurring_tan/2.35);
1016   }                                              1016   }
1017   //The main purpose of shifting the response    1017   //The main purpose of shifting the response is to be able to project the response of all the crytal elements into a 2D position histogram so that we can identify the DOI layer 
1018   //by comparing with a look-up-table which i    1018   //by comparing with a look-up-table which is prepared based on the reflector insertion. 
1019                                                  1019 
1020   //The crystal ID in 2D position histogram a    1020   //The crystal ID in 2D position histogram along the axial (z) direction. It can have values of: 0, 1, .. , 31, in 32x32 pixel position histogram
1021   crystalID_in2D_posHist_axial = (G4int)(Posi    1021   crystalID_in2D_posHist_axial = (G4int)(PositionAngerZ/(crystalPitch_axial*0.5) + (G4double)numberOfPixel_axial*0.5);//Note! crystalPitch_axial*0.5 is the pitch for the 32x32 2D pixel space, and 0.5 is added for round off
1022                                                  1022 
1023   //The crystal ID in 2D position histogram a    1023   //The crystal ID in 2D position histogram along the tangential (y) direction. It can have values of: 0, 1, .. , 31, in 32x32 pixel position histogram
1024   crystalID_in2D_posHist_tan =   (G4int)(Posi    1024   crystalID_in2D_posHist_tan =   (G4int)(PositionAngerY/(crystalPitch_tan*0.5) + (G4double)numberOfPixel_tan * 0.5);
1025                                                  1025   
1026   //continuous crystal ID in the 2D position     1026   //continuous crystal ID in the 2D position histogram. It will be from 0 to 1023 (in the case of 16x16x4 crystal array). 
1027   crystalID_in2D_posHist = crystalID_in2D_pos    1027   crystalID_in2D_posHist = crystalID_in2D_posHist_axial + crystalID_in2D_posHist_tan * numberOfPixel_tan;//32;
1028                                                  1028 
1029                                                  1029 
1030   //Now, lets find the crystal ID in 3D after    1030   //Now, lets find the crystal ID in 3D after applying Anger Logic calculation. NOTE that its value can be the same as the original crystal ID or not.
1031                                                  1031 
1032   //Crystal ID along the tangential diraction    1032   //Crystal ID along the tangential diraction after Anger Logic calculation
1033   crystalIDNew_tan = (G4int)(crystalID_in2D_p    1033   crystalIDNew_tan = (G4int)(crystalID_in2D_posHist_tan/2);
1034                                                  1034 
1035   //Crystal ID along the axial diraction afte    1035   //Crystal ID along the axial diraction after Anger Logic calculation
1036   crystalIDNew_axial = (G4int)(crystalID_in2D    1036   crystalIDNew_axial = (G4int)(crystalID_in2D_posHist_axial/2);
1037                                                  1037 
1038   ////Crystal ID along the DOI diraction afte    1038   ////Crystal ID along the DOI diraction after Anger Logic calculation
1039   if(crystalID_in2D_posHist>numberOfCrystal_t    1039   if(crystalID_in2D_posHist>numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI) return;
1040   crystalIDNew_DOI = doi_table[crystalID_in2D    1040   crystalIDNew_DOI = doi_table[crystalID_in2D_posHist];
1041                                                  1041 
1042   //If the crsytal ID is beyond the given the    1042   //If the crsytal ID is beyond the given the number of crystal in the detector, the following is excecuted and the event will be rejected
1043   if(crystalIDNew_tan < 0 || crystalIDNew_axi    1043   if(crystalIDNew_tan < 0 || crystalIDNew_axial < 0 || crystalIDNew_DOI < 0 ||
1044     crystalIDNew_tan >= numberOfCrystal_tange    1044     crystalIDNew_tan >= numberOfCrystal_tangential || crystalIDNew_axial >= numberOfCrystal_axial || crystalIDNew_DOI >= numberOfCrystal_DOI){
1045       return;                                    1045       return;
1046   }                                              1046   }
1047                                                  1047 
1048   CrystalIDAfterAngerLogic(crystalIDNew_tan,c    1048   CrystalIDAfterAngerLogic(crystalIDNew_tan,crystalIDNew_axial,crystalIDNew_DOI); 
1049 }                                                1049 }
1050                                                  1050 
1051 /////                                            1051 /////
1052 void doiPETAnalysis::CrystalIDAfterAngerLogic    1052 void doiPETAnalysis::CrystalIDAfterAngerLogic(G4int i_tan, G4int i_axial, G4int i_doi){
1053   crystalID_tangential = i_tan;                  1053   crystalID_tangential = i_tan;
1054   crystalID_axial = i_axial;                     1054   crystalID_axial = i_axial;
1055   DOI_ID = i_doi;                                1055   DOI_ID = i_doi;
1056 }                                                1056 }
1057                                                  1057 
1058 void doiPETAnalysis::book()                      1058 void doiPETAnalysis::book() 
1059 {                                                1059 { 
1060   auto manager = G4AnalysisManager::Instance(    1060   auto manager = G4AnalysisManager::Instance();
1061                                                  1061   
1062   //manager->SetVerboseLevel(2);                 1062   //manager->SetVerboseLevel(2);
1063                                                  1063  
1064   G4bool fileOpen = manager->OpenFile(rootFil    1064   G4bool fileOpen = manager->OpenFile(rootFileName);
1065   if (!fileOpen) {                               1065   if (!fileOpen) {
1066     G4cout << "\n---> HistoManager::book(): c    1066     G4cout << "\n---> HistoManager::book(): cannot open " 
1067            << rootFileName << G4endl;            1067            << rootFileName << G4endl;
1068     return;                                      1068     return;
1069   }                                              1069   }
1070   // Create directories                          1070   // Create directories  
1071   //manager->SetNtupleDirectoryName("ListMode    1071   //manager->SetNtupleDirectoryName("ListModeData");
1072                                                  1072 
1073   manager->SetFirstNtupleId(1);                  1073   manager->SetFirstNtupleId(1);
1074                                                  1074 
1075   if(getSinglesData){                            1075   if(getSinglesData){
1076   manager -> CreateNtuple("Singles", "Singles    1076   manager -> CreateNtuple("Singles", "Singles");
1077   fNtColId[0] = manager -> CreateNtupleIColum    1077   fNtColId[0] = manager -> CreateNtupleIColumn("eventID");
1078   fNtColId[1] = manager -> CreateNtupleIColum    1078   fNtColId[1] = manager -> CreateNtupleIColumn("blockID");
1079   //fNtColId[2] = manager -> CreateNtupleDCol    1079   //fNtColId[2] = manager -> CreateNtupleDColumn("crystalID_axial");
1080   //fNtColId[3] = manager -> CreateNtupleDCol    1080   //fNtColId[3] = manager -> CreateNtupleDColumn("crystalID_tangential");
1081   //fNtColId[4] = manager -> CreateNtupleDCol    1081   //fNtColId[4] = manager -> CreateNtupleDColumn("DOI_ID");
1082   fNtColId[2] = manager -> CreateNtupleDColum    1082   fNtColId[2] = manager -> CreateNtupleDColumn("timeStamp");
1083   fNtColId[3] = manager -> CreateNtupleDColum    1083   fNtColId[3] = manager -> CreateNtupleDColumn("totalEdep");
1084                                                  1084 
1085   //Interaction position of the photon with t    1085   //Interaction position of the photon with the detector
1086   fNtColId[4] = manager -> CreateNtupleDColum    1086   fNtColId[4] = manager -> CreateNtupleDColumn("intPosX");
1087   fNtColId[5] = manager -> CreateNtupleDColum    1087   fNtColId[5] = manager -> CreateNtupleDColumn("intPosY");
1088   fNtColId[6] = manager -> CreateNtupleDColum    1088   fNtColId[6] = manager -> CreateNtupleDColumn("intPosZ");
1089                                                  1089 
1090   ////source position (annihilation position)    1090   ////source position (annihilation position)
1091   fNtColId[7] = manager -> CreateNtupleDColum    1091   fNtColId[7] = manager -> CreateNtupleDColumn("spositionX");
1092   fNtColId[8] = manager -> CreateNtupleDColum    1092   fNtColId[8] = manager -> CreateNtupleDColumn("spositionY");
1093   fNtColId[9] = manager -> CreateNtupleDColum    1093   fNtColId[9] = manager -> CreateNtupleDColumn("spositionZ");
1094                                                  1094 
1095   manager -> FinishNtuple();                     1095   manager -> FinishNtuple();    
1096   }                                              1096   }
1097   if(getCoincidenceData){                        1097   if(getCoincidenceData){
1098     manager -> CreateNtuple("Coincidence", "C    1098     manager -> CreateNtuple("Coincidence", "Coincidence");
1099   fNtColId[0] = manager -> CreateNtupleIColum    1099   fNtColId[0] = manager -> CreateNtupleIColumn("eventID0");
1100   fNtColId[1] = manager -> CreateNtupleIColum    1100   fNtColId[1] = manager -> CreateNtupleIColumn("blockID0");
1101   fNtColId[2] = manager -> CreateNtupleIColum    1101   fNtColId[2] = manager -> CreateNtupleIColumn("crystalID_axial0");
1102   fNtColId[3] = manager -> CreateNtupleIColum    1102   fNtColId[3] = manager -> CreateNtupleIColumn("crystalID_tangential0");
1103   fNtColId[4] = manager -> CreateNtupleIColum    1103   fNtColId[4] = manager -> CreateNtupleIColumn("DOI_ID0");
1104   fNtColId[5] = manager -> CreateNtupleDColum    1104   fNtColId[5] = manager -> CreateNtupleDColumn("timeStamp0");
1105   fNtColId[6] = manager -> CreateNtupleDColum    1105   fNtColId[6] = manager -> CreateNtupleDColumn("totalEdep0");
1106                                                  1106 
1107   fNtColId[7] = manager -> CreateNtupleIColum    1107   fNtColId[7] = manager -> CreateNtupleIColumn("eventID1");
1108   fNtColId[8] = manager -> CreateNtupleIColum    1108   fNtColId[8] = manager -> CreateNtupleIColumn("blockID1");
1109   fNtColId[9] = manager -> CreateNtupleIColum    1109   fNtColId[9] = manager -> CreateNtupleIColumn("crystalID_axial1");
1110   fNtColId[10] = manager -> CreateNtupleIColu    1110   fNtColId[10] = manager -> CreateNtupleIColumn("crystalID_tangential1");
1111   fNtColId[11] = manager -> CreateNtupleIColu    1111   fNtColId[11] = manager -> CreateNtupleIColumn("DOI_ID1");
1112   fNtColId[12] = manager -> CreateNtupleDColu    1112   fNtColId[12] = manager -> CreateNtupleDColumn("timeStamp1");
1113   fNtColId[13] = manager -> CreateNtupleDColu    1113   fNtColId[13] = manager -> CreateNtupleDColumn("totalEdep1");
1114                                                  1114 
1115   //source position                              1115   //source position
1116   fNtColId[14] = manager -> CreateNtupleDColu    1116   fNtColId[14] = manager -> CreateNtupleDColumn("spositionX");
1117   fNtColId[15] = manager -> CreateNtupleDColu    1117   fNtColId[15] = manager -> CreateNtupleDColumn("spositionY");
1118   fNtColId[16] = manager -> CreateNtupleDColu    1118   fNtColId[16] = manager -> CreateNtupleDColumn("spositionZ");
1119                                                  1119 
1120   manager -> FinishNtuple();                     1120   manager -> FinishNtuple();
1121                                                  1121 
1122   }                                              1122   }
1123                                                  1123   
1124                                                  1124   
1125   factoryOn = true;                              1125   factoryOn = true;    
1126 }                                                1126 }
1127 void doiPETAnalysis::FillListModeEvent()         1127 void doiPETAnalysis::FillListModeEvent()
1128 {                                                1128 {
1129                                                  1129 
1130   auto manager = G4AnalysisManager::Instance(    1130   auto manager = G4AnalysisManager::Instance();
1131   if(getSinglesData){                            1131   if(getSinglesData){
1132     manager -> FillNtupleIColumn(1, fNtColId[    1132     manager -> FillNtupleIColumn(1, fNtColId[0], G4int(eventID));
1133     manager -> FillNtupleIColumn(1, fNtColId[    1133     manager -> FillNtupleIColumn(1, fNtColId[1], G4int(blockID));
1134     //manager -> FillNtupleDColumn(1, fNtColI    1134     //manager -> FillNtupleDColumn(1, fNtColId[2], crystalID_axial);
1135     //manager -> FillNtupleDColumn(1, fNtColI    1135     //manager -> FillNtupleDColumn(1, fNtColId[3], crystalID_tangential);
1136     //manager -> FillNtupleDColumn(1, fNtColI    1136     //manager -> FillNtupleDColumn(1, fNtColId[4], DOI_ID);
1137     manager -> FillNtupleDColumn(1, fNtColId[    1137     manager -> FillNtupleDColumn(1, fNtColId[2], timeStamp/s);// in second
1138     manager -> FillNtupleDColumn(1, fNtColId[    1138     manager -> FillNtupleDColumn(1, fNtColId[3], totalEdep/keV); //in keV
1139                                                  1139     
1140     //Interaction position of the photon in t    1140     //Interaction position of the photon in the detector
1141     manager -> FillNtupleDColumn(1, fNtColId[    1141     manager -> FillNtupleDColumn(1, fNtColId[4], intPosX); //mm
1142     manager -> FillNtupleDColumn(1, fNtColId[    1142     manager -> FillNtupleDColumn(1, fNtColId[5], intPosY); //mm
1143     manager -> FillNtupleDColumn(1, fNtColId[    1143     manager -> FillNtupleDColumn(1, fNtColId[6], intPosZ); //mm
1144                                                  1144 
1145     //                                           1145     //
1146     //Add source position                        1146     //Add source position
1147     manager -> FillNtupleDColumn(1, fNtColId[    1147     manager -> FillNtupleDColumn(1, fNtColId[7], spositionX);
1148     manager -> FillNtupleDColumn(1, fNtColId[    1148     manager -> FillNtupleDColumn(1, fNtColId[8], spositionY);
1149     manager -> FillNtupleDColumn(1, fNtColId[    1149     manager -> FillNtupleDColumn(1, fNtColId[9], spositionZ);
1150                                                  1150 
1151     manager -> AddNtupleRow(1);                  1151     manager -> AddNtupleRow(1);
1152   }                                              1152   }
1153                                                  1153 
1154   if(getCoincidenceData){                        1154   if(getCoincidenceData){
1155     //First Single                               1155     //First Single
1156     manager -> FillNtupleIColumn(1, fNtColId[    1156     manager -> FillNtupleIColumn(1, fNtColId[0], eventID0);
1157     manager -> FillNtupleIColumn(1, fNtColId[    1157     manager -> FillNtupleIColumn(1, fNtColId[1], blockID0);
1158     manager -> FillNtupleIColumn(1, fNtColId[    1158     manager -> FillNtupleIColumn(1, fNtColId[2], crystalID_axial0);
1159     manager -> FillNtupleIColumn(1, fNtColId[    1159     manager -> FillNtupleIColumn(1, fNtColId[3], crystalID_tangential0);
1160     manager -> FillNtupleIColumn(1, fNtColId[    1160     manager -> FillNtupleIColumn(1, fNtColId[4], DOI_ID0);
1161     manager -> FillNtupleDColumn(1, fNtColId[    1161     manager -> FillNtupleDColumn(1, fNtColId[5], timeStamp0/s);
1162     manager -> FillNtupleDColumn(1, fNtColId[    1162     manager -> FillNtupleDColumn(1, fNtColId[6], totalEdep0/keV);
1163                                                  1163   
1164   //Second Single                                1164   //Second Single
1165     manager -> FillNtupleIColumn(1, fNtColId[    1165     manager -> FillNtupleIColumn(1, fNtColId[7], eventID1);
1166     manager -> FillNtupleIColumn(1, fNtColId[    1166     manager -> FillNtupleIColumn(1, fNtColId[8], blockID1);
1167     manager -> FillNtupleIColumn(1, fNtColId[    1167     manager -> FillNtupleIColumn(1, fNtColId[9], crystalID_axial1);
1168     manager -> FillNtupleIColumn(1, fNtColId[    1168     manager -> FillNtupleIColumn(1, fNtColId[10], crystalID_tangential1);
1169     manager -> FillNtupleIColumn(1, fNtColId[    1169     manager -> FillNtupleIColumn(1, fNtColId[11], DOI_ID1);
1170     manager -> FillNtupleDColumn(1, fNtColId[    1170     manager -> FillNtupleDColumn(1, fNtColId[12], timeStamp1/s);
1171     manager -> FillNtupleDColumn(1, fNtColId[    1171     manager -> FillNtupleDColumn(1, fNtColId[13], totalEdep1/keV);
1172                                                  1172   
1173   //Add source position                          1173   //Add source position
1174     manager -> FillNtupleDColumn(1, fNtColId[    1174     manager -> FillNtupleDColumn(1, fNtColId[14], spositionX);
1175     manager -> FillNtupleDColumn(1, fNtColId[    1175     manager -> FillNtupleDColumn(1, fNtColId[15], spositionY);
1176     manager -> FillNtupleDColumn(1, fNtColId[    1176     manager -> FillNtupleDColumn(1, fNtColId[16], spositionZ);
1177                                                  1177 
1178     manager -> AddNtupleRow(1);                  1178     manager -> AddNtupleRow(1);
1179   }                                              1179   }
1180                                                  1180     
1181 }                                                1181 }
1182 void doiPETAnalysis::finish()                    1182 void doiPETAnalysis::finish() 
1183 {                                                1183 {   
1184  if (factoryOn)                                  1184  if (factoryOn) 
1185    {                                             1185    {
1186     auto manager = G4AnalysisManager::Instanc    1186     auto manager = G4AnalysisManager::Instance();    
1187     manager -> Write();                          1187     manager -> Write();
1188     manager -> CloseFile();                      1188     manager -> CloseFile();  
1189                                                  1189       
1190    // delete G4AnalysisManager::Instance();      1190    // delete G4AnalysisManager::Instance();
1191     factoryOn = false;                           1191     factoryOn = false;
1192    }                                             1192    }
1193 }                                                1193 }