| Geant4 Cross Reference |
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 }