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