| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25
26 //GEANT4 - Depth-of-Interaction enabled Positron emission tomography (PET) advanced example
27
28 //Authors and contributors
29
30 // Author list to be updated, with names of co-authors and contributors from National Institute of Radiological Sciences (NIRS)
31
32 // Abdella M. Ahmed (1, 2), Andrew Chacon (1, 2), Harley Rutherford (1, 2),
33 // Hideaki Tashima (3), Go Akamatsu (3), Akram Mohammadi (3), Eiji Yoshida (3), Taiga Yamaya (3)
34 // Susanna Guatelli (2), and Mitra Safavi-Naeini (1, 2)
35
36 // (1) Australian Nuclear Science and Technology Organisation, Australia
37 // (2) University of Wollongong, Australia
38 // (3) National Institute of Radiological Sciences, Japan
39
40 //Implemetation of the doiPETAnalysis.cc class
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 file. Crystal based energy resolution and quantum efficiency has been applied. Deadtime (on
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) 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 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 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 (reflector based) DOI identification method has been used. If the crystal ID is out of bound,
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. As, an option, binary output method is also given.
50 //Explanation is given for the methods provided.
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[k] = 0;
71
72 fAnalysisMessenger = new doiPETAnalysisMessenger(this);
73
74 //Set energy window
75 lowerThreshold = 400*keV;
76 upperThreshold = 600*keV;
77 triggerEnergy = 50*eV;
78
79 //give default initial activity. Activity strength is changed in the .mac file
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 a given isotope can be changed via the run.mac file
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_tangential; //32;
101 numberOfPixel_axial = 2*numberOfCrystal_axial; //32;
102
103 //Default value for deadtime.
104 block_DeadTime = 256*ns;
105 module_DeadTime = 0*ns;
106 //
107
108 //Crystal blurring parameters. One detector has 1024 crystals. All the crystals have different energy resolution.
109 //So, a range of energy resolution is applied between minumun and maximum values.
110 //The energy resolution can be set in the inputParameter.txt file
111 crystalResolutionMin = 0.13;//13%
112 crystalResolutionMax = 0.17;//17%
113
114 crystalEnergyRef = 511 * keV;//Energy of reference in which the energy resolution of the crystal is computed
115
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 inputParameter.txt file
118 crystalQuantumEfficiency = 1;//100%
119 //
120
121 //intialize deadtime for blocks and modules
122 numberOfBlocks_total = numberOfRings * numberOfDetector_perRing;
123 blockTime = new double[numberOfBlocks_total];//for each individual block.
124 moduleTime = new double[numberOfBlocks_total];//for axially multiplexed detectors.
125
126 //Initialize the deadtime for each detector and axially multiplexed detector (also called modules)
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 output is single events
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 the shift of the response due to the reflector is half distance from the interaction position to the air gap.
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 doiPETAnalysis();
160 return instance;
161 }
162 void doiPETAnalysis::Delete()
163 {
164 delete instance;
165 }
166
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
169 void doiPETAnalysis::SetScatterIndexInPhantom(G4int sci){
170 scatterIndex = sci;
171 }
172
173 //Get the source position if the process is annihilation.
174 //Use this for checking
175 void doiPETAnalysis::SetSourcePosition(G4ThreeVector spos){
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 used to sort coincidence events if realistic simulation is done
184 void doiPETAnalysis::SetEventID(G4int evID){
185 eventID = evID;
186 }
187
188 //
189
190 void doiPETAnalysis::GetSizeOfDetector(G4double detSizeDoi, G4double detSizeTan, G4double detSizeAxial){
191 sizeOfDetector_DOI = detSizeDoi;
192 sizeOfDetector_axial = detSizeTan;
193 sizeOfDetector_tangential = detSizeAxial;
194 }
195
196 //
197 void doiPETAnalysis::SetActivity(G4double newActivity){
198 InitialActivity = newActivity;
199 G4cout<<"Initial activity: "<<InitialActivity/becquerel<<" Bq."<<G4endl;
200 }
201 void doiPETAnalysis::SetIsotopeHalfLife(G4double newHalfLife){
202 halfLife = newHalfLife;
203 G4cout<<"Half life of the isotope "<<halfLife/s<<" sec."<<G4endl;
204 }
205
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 scanner for a given number of particles.
208 void doiPETAnalysis::CalulateAcquisitionTime(){
209 //Calculate the strength of activity at t=totaltime using decay equation
210 activityNow = InitialActivity * std::exp(-((0.693147/halfLife)*totalTime)); //ln(2) = 0.693147181
211
212 //Activity based time interval.
213 timeInterval = -std::log(G4UniformRand())*(1./activityNow);
214 totalTime = timeInterval+prev_totalTime;
215 prev_totalTime = totalTime;
216 }
217
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(G4double edep, G4int blkID, G4int cysID)
220 {
221 if(fixedResolution){
222 crystalResolution = energyResolution_fixed;
223 }
224 else{
225 crystalResolution = energyResolution_cryDependent[blkID][cysID];
226 }
227 crystalCoeff = crystalResolution * std::sqrt(crystalEnergyRef);
228
229 G4double QE = G4UniformRand();
230
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)
233 {
234 edep_AfterCrystalBlurring = G4RandGauss::shoot(edep,crystalCoeff*std::sqrt(edep)/2.35);
235 }
236 else {
237 //not detected by the photodetector, eventhough there was an interaction
238 edep_AfterCrystalBlurring = 0 *keV;
239 }
240 return edep_AfterCrystalBlurring;
241 }
242
243 ///////// ReadOut ///////////////////////////////////
244
245 void doiPETAnalysis::ReadOut(G4int blkID, G4int cryID, G4double interTime, G4double timeAnnih, G4ThreeVector interPos, G4double edep)
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 duration from the annihilation process to the detection of the photons by the scintillator.
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 the detector below which the detector is insensitive to any interaction.
261 if(totalEdep<triggerEnergy)return;
262
263 //************************************** Apply dead-time ********************************************//
264 //Apply paralizable dead-time in the block beofore events are rejected by the energy window
265 if(std::fabs(timeStamp - blockTime[blockID]) >= block_DeadTime){ //If true, the event is accepted
266 blockTime[blockID] = timeStamp;
267 }
268 else {
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;
271
272 //the event is then rejected
273 //continue;
274 return;
275 }
276
277 //Apply Non-paralyzable dead-time on axially multiplexed detectors (4 detectors are arranged axailly)
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]) > module_DeadTime){
280
281 //The following finds the block id's of four blocks which are arranged axially
282 for (G4int r_ring = 0; r_ring < numberOfRings; r_ring++){
283 if (blockID >= r_ring*numberOfDetector_perRing && blockID <(r_ring + 1)*numberOfDetector_perRing){
284 for (G4int m_module = 0; m_module < numberOfRings; m_module++){
285
286 //Set the time of the module (four blocks) the same
287 moduleTime[blockID + (m_module - r_ring)*numberOfDetector_perRing] = timeStamp;
288 }
289 }
290 }
291 }
292 else return;
293
294
295 ///////////////////////// Write qualified single events based the energy deposition in the detector ///////////
296
297 if(totalEdep>lowerThreshold && totalEdep<upperThreshold ){
298
299 //identifiy the layer
300 DOI_ID = G4int(crystalID/(numberOfCrystal_tangential * numberOfCrystal_axial));
301
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*numberOfCrystal_tangential * numberOfCrystal_axial);
304
305 //identify the crystal ID in the tangential and axial direction
306 crystalID_axial = crystalID_2D/numberOfCrystal_axial;
307 crystalID_tangential = crystalID_2D%numberOfCrystal_tangential;
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 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, totalEdep, shiftCoeff, isDOIlookUpTablePrepared);//
317 }
318
319 //Single event output. Coincidence events can then be made using the single 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_axial);
327 cryID_tan_coin.push_back(crystalID_tangential);
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 within the energy window are qualified.
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.data";
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 write the output ===="<<G4endl;
373 exit(0);
374 }
375
376 }
377
378 void doiPETAnalysis::WriteOutput(){
379 if(getSinglesData){
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<<" "<<crystalID_axial<<" "<<crystalID_tangential<<" "<<DOI_ID<<" "<<std::setprecision(17)<<timeStamp/s<<" "<<std::setprecision(7)<<totalEdep/keV<<G4endl;
382
383 }
384 if(getCoincidenceData){
385 //2 singles will qualify to be in coincidence within the energy window.
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[0];
393 crystalID_tangential0 = cryID_tan_coin[0];
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[1];
403 crystalID_tangential1 = cryID_tan_coin[1];
404 DOI_ID1 = cryDOI_coin[1];
405 timeStamp1 = time_coin[1];
406 totalEdep1 = edep_coin[1];
407 }
408 }
409
410 ofs<<eventID0<<" "<<blockID0<<" "<<crystalID_axial0<<" "<<crystalID_tangential0<<" "<<DOI_ID0<<" "<<std::setprecision(17)<<timeStamp0/s<<" "<<std::setprecision(7)<<totalEdep0/keV<<" "
411 <<eventID1<<" "<<blockID1<<" "<<crystalID_axial1<<" "<<crystalID_tangential1<<" "<<DOI_ID1<<" "<<std::setprecision(17)<<timeStamp1/s<<" "<<std::setprecision(7)<<totalEdep1/keV<<" "
412 <<spositionX<<" "<<spositionY<<" "<<spositionZ<<G4endl;
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 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) position (at +sizeOfDetector_DOI/2 which is at the top of the detector).
432
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 depends on the distance of the respective PMT from the interaction point
435 void doiPETAnalysis::PMTPosition(){
436
437 sizeOfDetector_DOI = (numberOfCrystal_DOI * sizeOfCrystal_DOI) + (numberOfCrystal_DOI - 1)*crystalGap_DOI;
438 sizeOfDetector_axial = (numberOfCrystal_axial * sizeOfCrystal_axial) + (numberOfCrystal_axial - 1)*crystalGap_axial;
439 sizeOfDetector_tangential = (numberOfCrystal_tangential * sizeOfCrystal_tangential) + (numberOfCrystal_tangential - 1)*crystalGap_tangential;
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<<", "<<posPMT1y<<", "<<posPMT1z<<")"<<G4endl;
463 G4cout<<"PMT2 (mm) ("<<posPMT2x<<", "<<posPMT2y<<", "<<posPMT2z<<")"<<G4endl;
464 G4cout<<"PMT3 (mm) ("<<posPMT3x<<", "<<posPMT3y<<", "<<posPMT3z<<")"<<G4endl;
465 G4cout<<"PMT4 (mm) ("<<posPMT4x<<", "<<posPMT4y<<", "<<posPMT4z<<")"<<G4endl;
466
467 }
468
469 //The blurring parameters are given and can be changed in the inputParameter.txt file
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 open "<<filename<<G4endl;
478 exit(0);
479 }
480 while(!ifs.eof()){
481 ifs.getline(inputChar,256);
482 inputLine = inputChar;
483 if(inputChar[0]!='#' && inputLine.length()!=0 ){
484 if( (std::string::size_type)inputLine.find("block_DeadTime:")!=std::string::npos){
485 std::istringstream tmpStream(inputLine);
486 tmpStream >> value[0] >> value[1] >> value[2];
487 block_DeadTime = atof(value[1].c_str());
488 if(value[2] != "ns"){
489 G4cerr<<" Dead time unit is not in nano seconds (ns), Make it in 'ns' "<<G4endl;
490 exit(0);
491 }
492 block_DeadTime = block_DeadTime*ns;
493 G4cout<<"Dead time of the detector: "<<block_DeadTime <<" ns."<<G4endl;
494 }
495 if( (std::string::size_type)inputLine.find("module_DeadTime:")!=std::string::npos){
496 std::istringstream tmpStream(inputLine);
497 tmpStream >> value[0] >> value[1] >> value[2];
498 module_DeadTime = atof(value[1].c_str());
499 if(value[2] != "ns"){
500 G4cerr<<" Dead time unit is not in nano seconds (ns), Make it in 'ns' "<<G4endl;
501 exit(0);
502 }
503 module_DeadTime = module_DeadTime*ns;
504 G4cout<<"Dead time of the module (axially multiplexed detectors): "<<module_DeadTime <<" ns."<<G4endl;
505 }
506 //
507 if( (std::string::size_type)inputLine.find("crystalResolutionMin:")!=std::string::npos){
508 std::istringstream tmpStream(inputLine);
509 tmpStream >> value[0] >> value[1];
510 crystalResolutionMin = atof(value[1].c_str());
511 G4cout<<"crystal Resolution (Min.): "<<crystalResolutionMin*100<< " %." <<G4endl;
512 }
513 if( (std::string::size_type)inputLine.find("crystalResolutionMax:")!=std::string::npos){
514 std::istringstream tmpStream(inputLine);
515 tmpStream >> value[0] >> value[1];
516 crystalResolutionMax = atof(value[1].c_str());
517 G4cout<<"crystal Resolution (Max.): "<<crystalResolutionMax*100<<" %"<<G4endl;
518 }
519
520 //
521 if( (std::string::size_type)inputLine.find("fixedResolution:")!=std::string::npos){
522 std::istringstream tmpStream(inputLine);
523 tmpStream >> value[0] >> value[1];
524 if(value[1]=="true"){
525 fixedResolution = true;
526 energyResolution_fixed = (crystalResolutionMin + crystalResolutionMax)*0.5;
527 G4cout<<"Fixed crystal resolution is used. "<<G4endl;
528 }
529 else {
530 fixedResolution = false;
531 //Store into a file if needed.
532 std::string fname = "crystalDependentResolution.txt";
533 std::ofstream outFname(fname.c_str());
534
535 G4cout<<" \n Crystal dependent resolution is used. preparing look-up table .... "<<G4endl;
536 energyResolution_cryDependent.resize(numberOfBlocks_total,std::vector<G4double>(numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI,0));
537 for(G4int i_blk = 0; i_blk < numberOfBlocks_total; i_blk++){
538 for(G4int i_cry = 0; i_cry < numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI; i_cry++){
539 energyResolution_cryDependent[i_blk][i_cry] = crystalResolutionMin + (crystalResolutionMax - crystalResolutionMin)*G4UniformRand();
540 //store into a file
541 outFname<<i_blk<<" "<<i_cry<<" "<<energyResolution_cryDependent[i_blk][i_cry]<<G4endl;
542 }
543 }
544 G4cout<<"Done. \n"<<G4endl;
545 outFname.close();
546 }
547
548 }
549 //
550
551 if( (std::string::size_type)inputLine.find("crystalEnergyRef:")!=std::string::npos){
552 std::istringstream tmpStream(inputLine);
553 tmpStream >> value[0] >> value[1] >> value[2];
554 crystalEnergyRef = atof(value[1].c_str());
555 if(value[2] != "keV"){
556 G4cerr<<" The unit of reference energy is not in keV, Make it in 'keV' "<<G4endl;
557 exit(0);
558 }
559 crystalEnergyRef = crystalEnergyRef*keV;
560 G4cout<<"Energy of refernce: "<<crystalEnergyRef/keV<<" keV."<<G4endl;
561 }
562 if( (std::string::size_type)inputLine.find("crystalQuantumEfficiency:")!=std::string::npos){
563 std::istringstream tmpStream(inputLine);
564 tmpStream >> value[0] >> value[1];
565 crystalQuantumEfficiency = atof(value[1].c_str());
566 G4cout<<"Quantum Efficiency "<<crystalQuantumEfficiency*100<< " % "<<G4endl;
567 }
568 if( (std::string::size_type)inputLine.find("lowerThreshold:")!=std::string::npos){
569 std::istringstream tmpStream(inputLine);
570 tmpStream >> value[0] >> value[1] >> value[2];
571 lowerThreshold = atof(value[1].c_str());
572 if(value[2] != "keV"){
573 G4cerr<<" The unit of Lower energy threshold is not in keV, Make it in 'keV' "<<G4endl;
574 exit(0);
575 }
576 lowerThreshold = lowerThreshold*keV;
577 G4cout<<"Lower energy threshold: "<<lowerThreshold/keV<<" keV."<<G4endl;
578
579 }
580 if( (std::string::size_type)inputLine.find("upperThreshold:")!=std::string::npos){
581 std::istringstream tmpStream(inputLine);
582 tmpStream >> value[0] >> value[1] >> value[2];
583 upperThreshold = atof(value[1].c_str());
584 if(value[2] != "keV"){
585 G4cerr<<" The unit of Upper energy threshold is not in keV, Make it in 'keV' "<<G4endl;
586 exit(0);
587 }
588 upperThreshold = upperThreshold*keV;
589 G4cout<<"Upper energy threshold: "<<upperThreshold/keV<<" keV."<<G4endl;
590
591 }
592
593 //
594 if( (std::string::size_type)inputLine.find("triggerEnergy:")!=std::string::npos){
595 std::istringstream tmpStream(inputLine);
596 tmpStream >> value[0] >> value[1] >> value[2];
597 triggerEnergy = atof(value[1].c_str());
598 if(value[2] != "keV"){
599 G4cerr<<" The unit of Trigger energy threshold is not in keV, Make it in 'keV' "<<G4endl;
600 exit(0);
601 }
602 triggerEnergy = triggerEnergy*keV;
603 G4cout<<"Trigger energy threshold: "<<triggerEnergy/keV<<" keV."<<G4endl;
604
605 }
606
607 //Option to apply AngerLogic
608 if( (std::string::size_type)inputLine.find("ApplyAngerLogic:")!=std::string::npos){
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 applied. "<<G4endl;
614 }
615 else if(value[1]=="false") {
616 ApplyAngerLogic = false;
617 G4cout<<"Angler Logic calculation is NOT applied. "<<G4endl;
618 }
619 else {
620 ApplyAngerLogic = true;
621 G4cout<<"Angler Logic calculation is applied (by defualt). "<<G4endl;
622 }
623 }
624
625 //PMT position calculation blurring at FWHM
626 if( (std::string::size_type)inputLine.find("PMTblurring:")!=std::string::npos){
627 std::istringstream tmpStream(inputLine);
628 tmpStream >> value[0] >> value[1]>>value[2];
629 PMTblurring_tan = atof(value[1].c_str());
630 PMTblurring_axial = atof(value[2].c_str());
631 G4cout<<"PMTblurring position response blurring at FWHM (tan x axial) "<<PMTblurring_tan<<" x " <<PMTblurring_axial<<" mm2"<<G4endl;
632 }
633
634 //
635 if( (std::string::size_type)inputLine.find("TypeOfOutput:")!=std::string::npos){
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. "<<G4endl;
641 }
642 else if(value[1]=="coincidenceOutput") {
643 getCoincidenceData = true;
644 G4cout<<"Coicidence mode output enabled. "<<G4endl;
645 }
646
647 }
648
649 }
650 }
651 ifs.close();
652 }
653
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://link.springer.com/article/10.1007/s12194-013-0231-4
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 reflector the value is 1 and if there is no reflector, the value is 0.
658
659 void doiPETAnalysis::ReadReflectorPattern(){
660 G4cout<<" Reflector pattern is being read "<<G4endl;
661 //
662 std::vector<std::string> stringReflectorValue;
663 //
664 char inputChar[256];
665 std::string inputLine;
666
667 //open inputParameter.txt to read reflector pattern.
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 open "<<filename<<G4endl;
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 inputparamter.txt file
682 if(inputChar[0]!='#' && inputLine.length()!=0 ){
683
684 //Reflector patter for Layer1 in the tangential direction
685 if( (std::string::size_type)inputLine.find("reflectorLayer1_Tangential:")!=std::string::npos){
686 std::istringstream tmpStream(inputLine);
687 while(tmpStream >> refValue){
688 stringReflectorValue.push_back(refValue);
689 if(stringReflectorValue.size()>1){
690 G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
691 ireflectorLayer1_Tangential.push_back(tmp_value);
692 }
693 }
694 stringReflectorValue.clear();
695 }
696
697
698 //Reflector patter for Layer1 in the axial direction
699 if( (std::string::size_type)inputLine.find("reflectorLayer1_Axial:")!=std::string::npos){
700 std::istringstream tmpStream(inputLine);
701 while(tmpStream >> refValue){
702 stringReflectorValue.push_back(refValue);
703 if(stringReflectorValue.size()>1){
704 G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
705 ireflectorLayer1_Axial.push_back(tmp_value);
706 }
707 }
708 stringReflectorValue.clear();
709 }
710
711 //Reflector patter for Layer2 in the tangential direction
712 if( (std::string::size_type)inputLine.find("reflectorLayer2_Tangential:")!=std::string::npos){
713 std::istringstream tmpStream(inputLine);
714 while(tmpStream >> refValue){
715 stringReflectorValue.push_back(refValue);
716 if(stringReflectorValue.size()>1){
717 G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
718 ireflectorLayer2_Tangential.push_back(tmp_value);
719 }
720 }
721 stringReflectorValue.clear();
722 }
723
724 //Reflector patter for Layer2 in the axial direction
725 if( (std::string::size_type)inputLine.find("reflectorLayer2_Axial:")!=std::string::npos){
726 std::istringstream tmpStream(inputLine);
727 while(tmpStream >> refValue){
728 stringReflectorValue.push_back(refValue);
729 if(stringReflectorValue.size()>1){
730 G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
731 ireflectorLayer2_Axial.push_back(tmp_value);
732 }
733 }
734 stringReflectorValue.clear();
735 }
736
737 //Reflector patter for Layer3 in the tangential direction
738 if( (std::string::size_type)inputLine.find("reflectorLayer3_Tangential:")!=std::string::npos){
739 std::istringstream tmpStream(inputLine);
740 while(tmpStream >> refValue){
741 stringReflectorValue.push_back(refValue);
742 if(stringReflectorValue.size()>1){
743 G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
744 ireflectorLayer3_Tangential.push_back(tmp_value);
745 }
746 }
747 stringReflectorValue.clear();
748 }
749
750 //Reflector patter for Layer3 in the axial direction
751 if( (std::string::size_type)inputLine.find("reflectorLayer3_Axial:")!=std::string::npos){
752 std::istringstream tmpStream(inputLine);
753 while(tmpStream >> refValue){
754 stringReflectorValue.push_back(refValue);
755 if(stringReflectorValue.size()>1){
756 G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
757 ireflectorLayer3_Axial.push_back(tmp_value);
758 }
759 }
760 stringReflectorValue.clear();
761 }
762
763 //Reflector patter for Layer4 in the tangential direction
764 if( (std::string::size_type)inputLine.find("reflectorLayer4_Tangential:")!=std::string::npos){
765 std::istringstream tmpStream(inputLine);
766 while(tmpStream >> refValue){
767 stringReflectorValue.push_back(refValue);
768 if(stringReflectorValue.size()>1){
769 G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
770 ireflectorLayer4_Tangential.push_back(tmp_value);
771 }
772 }
773 stringReflectorValue.clear();
774 }
775
776 //Reflector patter for Layer4 in the axial direction
777 if( (std::string::size_type)inputLine.find("reflectorLayer4_Axial:")!=std::string::npos){
778 std::istringstream tmpStream(inputLine);
779 while(tmpStream >> refValue){
780 stringReflectorValue.push_back(refValue);
781 if(stringReflectorValue.size()>1){
782 G4int tmp_value = atoi(stringReflectorValue[stringReflectorValue.size()-1].c_str());
783 ireflectorLayer4_Axial.push_back(tmp_value);
784 }
785 }
786 stringReflectorValue.clear();
787 }
788 }//#
789 }//while(eof)
790
791 //
792 //for debug
793 G4cout<<"\n========= Reflector Pattern ==========="<<G4endl;
794 G4cout<<"Layer 1"<<G4endl;
795 for(unsigned int i = 0; i<ireflectorLayer1_Tangential.size();i++){
796 G4cout<<ireflectorLayer1_Tangential[i]<<" ";
797 }G4cout<<G4endl;
798 for(unsigned int i = 0; i<ireflectorLayer1_Axial.size();i++){
799 G4cout<<ireflectorLayer1_Axial[i]<<" ";
800 }G4cout<<G4endl;
801 G4cout<<"Layer 2"<<G4endl;
802 for(unsigned int i = 0; i<ireflectorLayer2_Tangential.size();i++){
803 G4cout<<ireflectorLayer2_Tangential[i]<<" ";
804 }G4cout<<G4endl;
805 for(unsigned int i = 0; i<ireflectorLayer2_Axial.size();i++){
806 G4cout<<ireflectorLayer2_Axial[i]<<" ";
807 }G4cout<<G4endl;
808 G4cout<<"Layer 3"<<G4endl;
809 for(unsigned int i = 0; i<ireflectorLayer3_Tangential.size();i++){
810 G4cout<<ireflectorLayer3_Tangential[i]<<" ";
811 }G4cout<<G4endl;
812 for(unsigned int i = 0; i<ireflectorLayer3_Axial.size();i++){
813 G4cout<<ireflectorLayer3_Axial[i]<<" ";
814 }G4cout<<G4endl;
815 G4cout<<"Layer 4"<<G4endl;
816 for(unsigned int i = 0; i<ireflectorLayer4_Tangential.size();i++){
817 G4cout<<ireflectorLayer4_Tangential[i]<<" ";
818 }G4cout<<G4endl;
819 for(unsigned int i = 0; i<ireflectorLayer4_Axial.size();i++){
820 G4cout<<ireflectorLayer4_Axial[i]<<" ";
821 }G4cout<<G4endl;
822 G4cout<<"========= Reflector Pattern Ended ===========\n"<<G4endl;
823
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;
826 doi_table.resize(numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI,0);
827
828 std::string LUT_FileName = "look_up_table_DOI.txt";
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: File name: "<<LUT_FileName<<G4endl;
834 //Read from file
835 while(ifs_doiLUT>>index_doi>>doiID && index_doi < int(doi_table.size())){
836 doi_table[index_doi] = doiID;
837 }
838 if(index_doi==int(doi_table.size())){
839 G4cout<<"!!!Warning: The DOI table index is greater than the total number of crystals."<<G4endl;
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_LUT"<<G4endl;
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(G4String){
863 isDOIlookUpTablePrepared = false;
864 G4cout<<"Preparing DOI look-up table... "<<G4endl;
865 std::string outputFileName = "_check_2Dposition.txt";// excuted only once
866 std::ofstream outFile(outputFileName.c_str());
867
868 G4double crystalPositionX;
869 G4double crystalPositionY;
870 G4double crystalPositionZ;
871 //doi_table.resize(numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI,0);
872
873 for(G4int i_DOI = 0; i_DOI<numberOfCrystal_DOI; i_DOI++){
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< numberOfCrystal_axial;i_axial++){
876 crystalPositionZ = (i_axial-((float)numberOfCrystal_axial)/2 + 0.5)*(sizeOfCrystal_axial + crystalGap_axial);
877 for(G4int i_tan=0; i_tan<numberOfCrystal_tangential;i_tan++){
878 crystalPositionY=(i_tan-((float)numberOfCrystal_tangential)/2 + 0.5)*(sizeOfCrystal_tangential + crystalGap_tangential);
879 AngerLogic(crystalPositionX, crystalPositionY, crystalPositionZ, 1, 0.5, isDOIlookUpTablePrepared);
880 outFile<<PositionAngerZ<<" "<<PositionAngerY<<G4endl;
881 doi_table[crystalID_in2D_posHist]=i_DOI;
882
883 }
884 }
885 }
886 G4cout<<"done."<<G4endl;
887 isDOIlookUpTablePrepared = true;
888 }
889
890
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 distance so that the response can be projected into 2D position histogram.
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 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 posCrystalX, G4double posCrystalY, G4double posCrystalZ, G4double Edep, G4double shiftDis, G4bool isDOI_LUT)
896 {
897 G4double crystalPitch_DOI = sizeOfCrystal_DOI + crystalGap_DOI;
898 G4double crystalPitch_tan = sizeOfCrystal_tangential + crystalGap_tangential;
899 G4double crystalPitch_axial = sizeOfCrystal_axial + crystalGap_axial;
900
901 //The crystal ID are calculated based on the center of mass
902 G4int i_doi = posCrystalX/crystalPitch_DOI + (float)numberOfCrystal_DOI*0.5;
903 G4int i_tan = posCrystalY/crystalPitch_tan + (float)numberOfCrystal_tangential*0.5;
904 G4int i_axial = posCrystalZ/crystalPitch_axial + (float)numberOfCrystal_axial*0.5;
905
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_DOI)/2 + 0.5)*crystalPitch_DOI;
908 posCrystalY = (i_tan-((float)numberOfCrystal_tangential)/2 + 0.5)*crystalPitch_tan;
909 posCrystalZ = (i_axial-((float)numberOfCrystal_axial)/2 + 0.5)*crystalPitch_axial;
910
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 devided into the four PMTs depending on their lateral distances (in the axial and tangential directions) from the interaction position
913
914 //The following is based on symetrical placment of the PMTs
915
916 //Calculate the axial (z) distance from the position of interaction to each PMT
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 only take two of them or take the average
924 //distz = ((dist1z + dist2z) + (dist3z + dist4z))/2;
925 distz = dist1z + dist3z;
926
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 from the position of interaction to each PMT
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 only take two of them or take the average
936 //disty = ((dist1y + dist3y) + (dist2y+dist4y))/2;
937 disty = dist1y + dist2y;
938
939 //signalPMT1z = signalPMT2z, and signalPMT3z = signalPMT4z
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 = signalPMT4y
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 'component' signal
955 signalPMT1 = (signalPMT1z + signalPMT1y)*0.5;
956 signalPMT2 = (signalPMT2z + signalPMT2y)*0.5;
957 signalPMT3 = (signalPMT3z + signalPMT3y)*0.5;
958 signalPMT4 = (signalPMT4z + signalPMT4y)*0.5;
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 based on Anger logic method.
968 //To get the position by Anger calculation, the result should be multiplied by the dimenion of the total distance.
969 PositionAngerZ = (signalZplus - signalZminus)/(signalZplus + signalZminus)*distz;
970 PositionAngerY = (signalYplus - signalYminus)/(signalYplus + signalYminus)*disty;
971
972
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 response is equal to half of the distance from the interaction position to the airgap in the lateral (transversal direction)
975
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 the crystal, then response shift is to the left side (away from the reflector).
978
979 //Response shift for 1st Layer
980 if(i_doi == 0){
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 && ireflectorLayer1_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2)*shiftDis;
983
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 && ireflectorLayer1_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
986
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 && ireflectorLayer1_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
989 }
990 if(i_doi == 1){ //Response shift for 2nd Layer
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 && ireflectorLayer2_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
993
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 && ireflectorLayer2_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
996 }
997 if(i_doi == 2){ //Response shift for 3rd Layer
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 && ireflectorLayer3_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
1000
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 && ireflectorLayer3_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
1003 }
1004 if(i_doi == 3){ //Response shift for 4th Layer
1005 if(ireflectorLayer4_Tangential[i_tan] == 1 && ireflectorLayer4_Tangential[i_tan + 1] == 0) PositionAngerY += (crystalPitch_tan/2)*shiftDis;
1006 if(ireflectorLayer4_Tangential[i_tan] == 0 && ireflectorLayer4_Tangential[i_tan + 1] == 1) PositionAngerY -= (crystalPitch_tan/2)*shiftDis;
1007
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 && ireflectorLayer4_Axial [i_axial + 1] == 1) PositionAngerZ -= (crystalPitch_axial/2)*shiftDis;
1010 }
1011
1012 //Blur the 2D position (obtained by ANger Logic method) to include uncertainity of the PMT position response.
1013 if(isDOI_LUT){
1014 PositionAngerZ = G4RandGauss::shoot(PositionAngerZ,PMTblurring_axial/2.35);
1015 PositionAngerY = G4RandGauss::shoot(PositionAngerY,PMTblurring_tan/2.35);
1016 }
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 is prepared based on the reflector insertion.
1019
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)(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
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)(PositionAngerY/(crystalPitch_tan*0.5) + (G4double)numberOfPixel_tan * 0.5);
1025
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_posHist_axial + crystalID_in2D_posHist_tan * numberOfPixel_tan;//32;
1028
1029
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
1032 //Crystal ID along the tangential diraction after Anger Logic calculation
1033 crystalIDNew_tan = (G4int)(crystalID_in2D_posHist_tan/2);
1034
1035 //Crystal ID along the axial diraction after Anger Logic calculation
1036 crystalIDNew_axial = (G4int)(crystalID_in2D_posHist_axial/2);
1037
1038 ////Crystal ID along the DOI diraction after Anger Logic calculation
1039 if(crystalID_in2D_posHist>numberOfCrystal_tangential*numberOfCrystal_axial*numberOfCrystal_DOI) return;
1040 crystalIDNew_DOI = doi_table[crystalID_in2D_posHist];
1041
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_axial < 0 || crystalIDNew_DOI < 0 ||
1044 crystalIDNew_tan >= numberOfCrystal_tangential || crystalIDNew_axial >= numberOfCrystal_axial || crystalIDNew_DOI >= numberOfCrystal_DOI){
1045 return;
1046 }
1047
1048 CrystalIDAfterAngerLogic(crystalIDNew_tan,crystalIDNew_axial,crystalIDNew_DOI);
1049 }
1050
1051 /////
1052 void doiPETAnalysis::CrystalIDAfterAngerLogic(G4int i_tan, G4int i_axial, G4int i_doi){
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(rootFileName);
1065 if (!fileOpen) {
1066 G4cout << "\n---> HistoManager::book(): cannot open "
1067 << rootFileName << G4endl;
1068 return;
1069 }
1070 // Create directories
1071 //manager->SetNtupleDirectoryName("ListModeData");
1072
1073 manager->SetFirstNtupleId(1);
1074
1075 if(getSinglesData){
1076 manager -> CreateNtuple("Singles", "Singles");
1077 fNtColId[0] = manager -> CreateNtupleIColumn("eventID");
1078 fNtColId[1] = manager -> CreateNtupleIColumn("blockID");
1079 //fNtColId[2] = manager -> CreateNtupleDColumn("crystalID_axial");
1080 //fNtColId[3] = manager -> CreateNtupleDColumn("crystalID_tangential");
1081 //fNtColId[4] = manager -> CreateNtupleDColumn("DOI_ID");
1082 fNtColId[2] = manager -> CreateNtupleDColumn("timeStamp");
1083 fNtColId[3] = manager -> CreateNtupleDColumn("totalEdep");
1084
1085 //Interaction position of the photon with the detector
1086 fNtColId[4] = manager -> CreateNtupleDColumn("intPosX");
1087 fNtColId[5] = manager -> CreateNtupleDColumn("intPosY");
1088 fNtColId[6] = manager -> CreateNtupleDColumn("intPosZ");
1089
1090 ////source position (annihilation position)
1091 fNtColId[7] = manager -> CreateNtupleDColumn("spositionX");
1092 fNtColId[8] = manager -> CreateNtupleDColumn("spositionY");
1093 fNtColId[9] = manager -> CreateNtupleDColumn("spositionZ");
1094
1095 manager -> FinishNtuple();
1096 }
1097 if(getCoincidenceData){
1098 manager -> CreateNtuple("Coincidence", "Coincidence");
1099 fNtColId[0] = manager -> CreateNtupleIColumn("eventID0");
1100 fNtColId[1] = manager -> CreateNtupleIColumn("blockID0");
1101 fNtColId[2] = manager -> CreateNtupleIColumn("crystalID_axial0");
1102 fNtColId[3] = manager -> CreateNtupleIColumn("crystalID_tangential0");
1103 fNtColId[4] = manager -> CreateNtupleIColumn("DOI_ID0");
1104 fNtColId[5] = manager -> CreateNtupleDColumn("timeStamp0");
1105 fNtColId[6] = manager -> CreateNtupleDColumn("totalEdep0");
1106
1107 fNtColId[7] = manager -> CreateNtupleIColumn("eventID1");
1108 fNtColId[8] = manager -> CreateNtupleIColumn("blockID1");
1109 fNtColId[9] = manager -> CreateNtupleIColumn("crystalID_axial1");
1110 fNtColId[10] = manager -> CreateNtupleIColumn("crystalID_tangential1");
1111 fNtColId[11] = manager -> CreateNtupleIColumn("DOI_ID1");
1112 fNtColId[12] = manager -> CreateNtupleDColumn("timeStamp1");
1113 fNtColId[13] = manager -> CreateNtupleDColumn("totalEdep1");
1114
1115 //source position
1116 fNtColId[14] = manager -> CreateNtupleDColumn("spositionX");
1117 fNtColId[15] = manager -> CreateNtupleDColumn("spositionY");
1118 fNtColId[16] = manager -> CreateNtupleDColumn("spositionZ");
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[0], G4int(eventID));
1133 manager -> FillNtupleIColumn(1, fNtColId[1], G4int(blockID));
1134 //manager -> FillNtupleDColumn(1, fNtColId[2], crystalID_axial);
1135 //manager -> FillNtupleDColumn(1, fNtColId[3], crystalID_tangential);
1136 //manager -> FillNtupleDColumn(1, fNtColId[4], DOI_ID);
1137 manager -> FillNtupleDColumn(1, fNtColId[2], timeStamp/s);// in second
1138 manager -> FillNtupleDColumn(1, fNtColId[3], totalEdep/keV); //in keV
1139
1140 //Interaction position of the photon in the detector
1141 manager -> FillNtupleDColumn(1, fNtColId[4], intPosX); //mm
1142 manager -> FillNtupleDColumn(1, fNtColId[5], intPosY); //mm
1143 manager -> FillNtupleDColumn(1, fNtColId[6], intPosZ); //mm
1144
1145 //
1146 //Add source position
1147 manager -> FillNtupleDColumn(1, fNtColId[7], spositionX);
1148 manager -> FillNtupleDColumn(1, fNtColId[8], spositionY);
1149 manager -> FillNtupleDColumn(1, fNtColId[9], spositionZ);
1150
1151 manager -> AddNtupleRow(1);
1152 }
1153
1154 if(getCoincidenceData){
1155 //First Single
1156 manager -> FillNtupleIColumn(1, fNtColId[0], eventID0);
1157 manager -> FillNtupleIColumn(1, fNtColId[1], blockID0);
1158 manager -> FillNtupleIColumn(1, fNtColId[2], crystalID_axial0);
1159 manager -> FillNtupleIColumn(1, fNtColId[3], crystalID_tangential0);
1160 manager -> FillNtupleIColumn(1, fNtColId[4], DOI_ID0);
1161 manager -> FillNtupleDColumn(1, fNtColId[5], timeStamp0/s);
1162 manager -> FillNtupleDColumn(1, fNtColId[6], totalEdep0/keV);
1163
1164 //Second Single
1165 manager -> FillNtupleIColumn(1, fNtColId[7], eventID1);
1166 manager -> FillNtupleIColumn(1, fNtColId[8], blockID1);
1167 manager -> FillNtupleIColumn(1, fNtColId[9], crystalID_axial1);
1168 manager -> FillNtupleIColumn(1, fNtColId[10], crystalID_tangential1);
1169 manager -> FillNtupleIColumn(1, fNtColId[11], DOI_ID1);
1170 manager -> FillNtupleDColumn(1, fNtColId[12], timeStamp1/s);
1171 manager -> FillNtupleDColumn(1, fNtColId[13], totalEdep1/keV);
1172
1173 //Add source position
1174 manager -> FillNtupleDColumn(1, fNtColId[14], spositionX);
1175 manager -> FillNtupleDColumn(1, fNtColId[15], spositionY);
1176 manager -> FillNtupleDColumn(1, fNtColId[16], spositionZ);
1177
1178 manager -> AddNtupleRow(1);
1179 }
1180
1181 }
1182 void doiPETAnalysis::finish()
1183 {
1184 if (factoryOn)
1185 {
1186 auto manager = G4AnalysisManager::Instance();
1187 manager -> Write();
1188 manager -> CloseFile();
1189
1190 // delete G4AnalysisManager::Instance();
1191 factoryOn = false;
1192 }
1193 }