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 7.0.p1)


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