| 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 40
41 #ifndef doiPETAnalysis_h 41 #ifndef doiPETAnalysis_h
42 #define doiPETAnalysis_h 1 42 #define doiPETAnalysis_h 1
43 43
44 #include "doiPETGlobalParameters.hh" 44 #include "doiPETGlobalParameters.hh"
45 #include "globals.hh" 45 #include "globals.hh"
46 #include <vector> 46 #include <vector>
47 #include <time.h> 47 #include <time.h>
48 #include <map> 48 #include <map>
49 #include <set> 49 #include <set>
50 #include "G4ThreeVector.hh" 50 #include "G4ThreeVector.hh"
51 #include <iostream> 51 #include <iostream>
52 #include <fstream> 52 #include <fstream>
53 #include <sstream> 53 #include <sstream>
54 #include <iterator> 54 #include <iterator>
55 #include <vector> 55 #include <vector>
56 #include <algorithm> 56 #include <algorithm>
57 #include "G4AnalysisManager.hh" 57 #include "G4AnalysisManager.hh"
58 58
59 // Define the total number of columns in the n 59 // Define the total number of columns in the ntuple
60 const G4int MaxNtCol = 17; 60 const G4int MaxNtCol = 17;
61 61
62 class doiPETAnalysisMessenger; 62 class doiPETAnalysisMessenger;
63 63
64 //class InteractionInformation; 64 //class InteractionInformation;
65 65
66 class doiPETAnalysis 66 class doiPETAnalysis
67 { 67 {
68 private: 68 private:
69 doiPETAnalysis(); 69 doiPETAnalysis();
70 70
71 public: 71 public:
72 ~doiPETAnalysis(); 72 ~doiPETAnalysis();
73 static doiPETAnalysis* GetInstance(); 73 static doiPETAnalysis* GetInstance();
74 void FindInteractingCrystal(); 74 void FindInteractingCrystal();
75 void Open(G4String); 75 void Open(G4String);
76 void Close(); 76 void Close();
77 void Delete(); 77 void Delete();
78 void ResetNumberOfHits(); 78 void ResetNumberOfHits();
79 void Write(/*G4int, G4int, G4int, G4double*/ 79 void Write(/*G4int, G4int, G4int, G4double*/);
80 void WriteOutput(); 80 void WriteOutput();
81 81
82 //void GetIntractionInfomation(InteractionIn 82 //void GetIntractionInfomation(InteractionInformation*);
83 83
84 void GetParentParticleName(G4String); 84 void GetParentParticleName(G4String);
85 void GetSizeOfDetector (G4double, G4double, 85 void GetSizeOfDetector (G4double, G4double, G4double);
86 void SetScatterIndexInPhantom(G4int); 86 void SetScatterIndexInPhantom(G4int);
87 87
88 void SetSourcePosition(G4ThreeVector);// 88 void SetSourcePosition(G4ThreeVector);//
89 void SetEventID(G4int); 89 void SetEventID(G4int);
90 90
91 void BlurringParameters(); 91 void BlurringParameters();
92 void GetTimeOfAnnihilation(G4double); 92 void GetTimeOfAnnihilation(G4double);
93 93
94 94
95 void PMTPosition(); 95 void PMTPosition();
96 void AngerLogic(G4double, G4double, G4double 96 void AngerLogic(G4double, G4double, G4double, G4double, G4double, G4bool);
97 void ReadReflectorPattern(); 97 void ReadReflectorPattern();
98 void PrepareDOILookUpTable(G4String); 98 void PrepareDOILookUpTable(G4String);
99 99
100 void SetActivity(G4double); 100 void SetActivity(G4double);
101 void SetIsotopeHalfLife(G4double); 101 void SetIsotopeHalfLife(G4double);
102 void CrystalIDAfterAngerLogic(G4int, G4int, 102 void CrystalIDAfterAngerLogic(G4int, G4int, G4int);
103 void TypeOfOutput(G4String);//Single or coin 103 void TypeOfOutput(G4String);//Single or coincidence list-mode data
104 void CalulateAcquisitionTime(); 104 void CalulateAcquisitionTime();
105 //G4double QuantumEffifciency(G4double); 105 //G4double QuantumEffifciency(G4double);
106 G4double QuantumEffifciency(G4double, G4int, 106 G4double QuantumEffifciency(G4double, G4int, G4int);
107 void ReadOut(G4int, G4int, G4double, G4doubl 107 void ReadOut(G4int, G4int, G4double, G4double, G4ThreeVector, G4double);
108 108
109 //G4ROOT 109 //G4ROOT
110 void book(); // booking the ROOT file 110 void book(); // booking the ROOT file
111 111
112 112
113 void FillListModeEvent(); //Single or Coinsi 113 void FillListModeEvent(); //Single or Coinsidence
114 void finish(); 114 void finish();
115 // Close the ROOT file with all the results 115 // Close the ROOT file with all the results stored in nutples
116 116
117 117
118 private: 118 private:
119 static doiPETAnalysis* instance; 119 static doiPETAnalysis* instance;
120 doiPETAnalysisMessenger* fAnalysisMessenger; 120 doiPETAnalysisMessenger* fAnalysisMessenger;
121 //std::multimap< G4int, InteractionInformati 121 //std::multimap< G4int, InteractionInformation* > mapBlockInteraction;
122 std::set<G4int> setBlockInteraction; 122 std::set<G4int> setBlockInteraction;
123 123
124 G4double upperThreshold, lowerThreshold; 124 G4double upperThreshold, lowerThreshold;
125 G4double triggerEnergy; 125 G4double triggerEnergy;
126 126
127 //G4ROOT 127 //G4ROOT
128 G4bool factoryOn; 128 G4bool factoryOn;
129 G4int fNtColId[MaxNtCol]; 129 G4int fNtColId[MaxNtCol];
130 // 130 //
131 131
132 132
133 //G4ThreeVector sourcePosition; 133 //G4ThreeVector sourcePosition;
134 134
135 135
136 136
137 // 137 //
138 G4int scatterIndex; 138 G4int scatterIndex;
139 139
140 G4String parentParticleName;// 140 G4String parentParticleName;//
141 141
142 // 142 //
143 G4int numberofInteractions; 143 G4int numberofInteractions;
144 G4int countCoincidence; 144 G4int countCoincidence;
145 145
146 G4int numberOfBlocks_total; 146 G4int numberOfBlocks_total;
147 147
148 G4double sizeOfDetector_DOI,sizeOfDetector_a 148 G4double sizeOfDetector_DOI,sizeOfDetector_axial,sizeOfDetector_tangential;
149 149
150 //Virtual position of the PMT 150 //Virtual position of the PMT
151 G4double signalPMT1, signalPMT2, signalPMT3, 151 G4double signalPMT1, signalPMT2, signalPMT3, signalPMT4;
152 152
153 G4double posPMT1x, posPMT2x, posPMT3x, posPM 153 G4double posPMT1x, posPMT2x, posPMT3x, posPMT4x;
154 G4double posPMT1y, posPMT2y, posPMT3y, posPM 154 G4double posPMT1y, posPMT2y, posPMT3y, posPMT4y;
155 G4double posPMT1z, posPMT2z, posPMT3z, posPM 155 G4double posPMT1z, posPMT2z, posPMT3z, posPMT4z;
156 156
157 // 157 //
158 G4double signalPMT1z, signalPMT2z, signalPMT 158 G4double signalPMT1z, signalPMT2z, signalPMT3z, signalPMT4z;
159 G4double signalPMT1y, signalPMT2y, signalPMT 159 G4double signalPMT1y, signalPMT2y, signalPMT3y, signalPMT4y;
160 160
161 // 161 //
162 G4double signalZplus, signalZminus; 162 G4double signalZplus, signalZminus;
163 G4double signalYplus, signalYminus; 163 G4double signalYplus, signalYminus;
164 // 164 //
165 165
166 G4double dist1z, dist2z, dist3z, dist4z, dis 166 G4double dist1z, dist2z, dist3z, dist4z, distz;
167 G4double dist1y, dist2y, dist3y, dist4y, dis 167 G4double dist1y, dist2y, dist3y, dist4y, disty;
168 168
169 G4double shiftCoeff; 169 G4double shiftCoeff;
170 170
171 G4double PositionAngerZ, PositionAngerY; 171 G4double PositionAngerZ, PositionAngerY;
172 172
173 //reflector pattern 173 //reflector pattern
174 std::vector<G4int> ireflectorLayer1_Tangenti 174 std::vector<G4int> ireflectorLayer1_Tangential;
175 std::vector<G4int> ireflectorLayer1_Axial; 175 std::vector<G4int> ireflectorLayer1_Axial;
176 std::vector<G4int> ireflectorLayer2_Tangenti 176 std::vector<G4int> ireflectorLayer2_Tangential;
177 std::vector<G4int> ireflectorLayer2_Axial; 177 std::vector<G4int> ireflectorLayer2_Axial;
178 std::vector<G4int> ireflectorLayer3_Tangenti 178 std::vector<G4int> ireflectorLayer3_Tangential;
179 std::vector<G4int> ireflectorLayer3_Axial; 179 std::vector<G4int> ireflectorLayer3_Axial;
180 std::vector<G4int> ireflectorLayer4_Tangenti 180 std::vector<G4int> ireflectorLayer4_Tangential;
181 std::vector<G4int> ireflectorLayer4_Axial; 181 std::vector<G4int> ireflectorLayer4_Axial;
182 std::vector<G4int> doi_table; 182 std::vector<G4int> doi_table;
183 // 183 //
184 184
185 //The number of pixes for the 2D position hi 185 //The number of pixes for the 2D position histogram after Anger Logic calculation
186 G4int numberOfPixel_axial; 186 G4int numberOfPixel_axial;
187 G4int numberOfPixel_tan; 187 G4int numberOfPixel_tan;
188 188
189 //source position 189 //source position
190 G4double spositionX; 190 G4double spositionX;
191 G4double spositionY; 191 G4double spositionY;
192 G4double spositionZ; 192 G4double spositionZ;
193 193
194 //interaction position with respect to the c 194 //interaction position with respect to the crystal axis
195 G4ThreeVector interactionPos; 195 G4ThreeVector interactionPos;
196 196
197 //interaction position 197 //interaction position
198 G4double intPosX; 198 G4double intPosX;
199 G4double intPosY; 199 G4double intPosY;
200 G4double intPosZ; 200 G4double intPosZ;
201 201
202 202
203 G4double interactionTime; 203 G4double interactionTime;
204 204
205 G4int crystalID;//contineous crystal ID in 3 205 G4int crystalID;//contineous crystal ID in 3D
206 G4int crystalID_2D; 206 G4int crystalID_2D;
207 207
208 G4int prev_eventID; 208 G4int prev_eventID;
209 209
210 //Single output 210 //Single output
211 G4int eventID; 211 G4int eventID;
212 G4int blockID; 212 G4int blockID;
213 G4int crystalID_axial; 213 G4int crystalID_axial;
214 G4int crystalID_tangential; 214 G4int crystalID_tangential;
215 G4int DOI_ID; 215 G4int DOI_ID;
216 G4double timeStamp; 216 G4double timeStamp;
217 G4double totalEdep; 217 G4double totalEdep;
218 218
219 //coincidence output 219 //coincidence output
220 G4int eventID0, eventID1; 220 G4int eventID0, eventID1;
221 G4int blockID0, blockID1; 221 G4int blockID0, blockID1;
222 G4int crystalID_axial0, crystalID_axial1 222 G4int crystalID_axial0, crystalID_axial1;
223 G4int crystalID_tangential0, crystalID_tang 223 G4int crystalID_tangential0, crystalID_tangential1;
224 G4int DOI_ID0, DOI_ID1; 224 G4int DOI_ID0, DOI_ID1;
225 G4double timeStamp0, timeStamp1; 225 G4double timeStamp0, timeStamp1;
226 G4double totalEdep0, totalEdep1; 226 G4double totalEdep0, totalEdep1;
227 G4double sposX, sposY, sposZ; 227 G4double sposX, sposY, sposZ;
228 //choice for the user 228 //choice for the user
229 G4bool getSinglesData; 229 G4bool getSinglesData;
230 G4bool getCoincidenceData; 230 G4bool getCoincidenceData;
231 231
232 // 232 //
233 G4bool ApplyAngerLogic; 233 G4bool ApplyAngerLogic;
234 234
235 G4double PMTblurring_tan; 235 G4double PMTblurring_tan;
236 G4double PMTblurring_axial; 236 G4double PMTblurring_axial;
237 237
238 G4String outputData; 238 G4String outputData;
239 G4int numberOfHit; 239 G4int numberOfHit;
240 std::vector<G4int> eventID_coin; 240 std::vector<G4int> eventID_coin;
241 std::vector<G4double> edep_coin; 241 std::vector<G4double> edep_coin;
242 std::vector<G4int>blockID_coin; 242 std::vector<G4int>blockID_coin;
243 std::vector<G4int> cryID_axial_coin; 243 std::vector<G4int> cryID_axial_coin;
244 std::vector<G4int> cryID_tan_coin; 244 std::vector<G4int> cryID_tan_coin;
245 std::vector<G4int> cryDOI_coin; 245 std::vector<G4int> cryDOI_coin;
246 std::vector<G4double> time_coin; 246 std::vector<G4double> time_coin;
247 247
248 //Crystal IDs after Anger Logic calculation 248 //Crystal IDs after Anger Logic calculation
249 G4int crystalIDNew_DOI, crystalIDNew_tan, cr 249 G4int crystalIDNew_DOI, crystalIDNew_tan, crystalIDNew_axial;
250 250
251 //Crystal ID in the 2D position histogram al 251 //Crystal ID in the 2D position histogram along the axial and tangetial direction
252 G4int crystalID_in2D_posHist_axial, crystalI 252 G4int crystalID_in2D_posHist_axial, crystalID_in2D_posHist_tan;
253 253
254 //continous crystal ID after after Anger Log 254 //continous crystal ID after after Anger Logic.
255 G4int crystalID_in2D_posHist; 255 G4int crystalID_in2D_posHist;
256 256
257 257
258 //Crystal blurring 258 //Crystal blurring
259 G4double crystalResolution; 259 G4double crystalResolution;
260 G4double crystalResolutionMin;// 260 G4double crystalResolutionMin;//
261 G4double crystalResolutionMax;// 261 G4double crystalResolutionMax;//
262 262
263 //G4bool variableResolution; 263 //G4bool variableResolution;
264 G4bool fixedResolution; 264 G4bool fixedResolution;
265 G4bool isDOIlookUpTablePrepared; 265 G4bool isDOIlookUpTablePrepared;
266 266
267 G4double energyResolution_fixed; 267 G4double energyResolution_fixed;
268 std::vector<std::vector<G4double>> energyRes 268 std::vector<std::vector<G4double>> energyResolution_cryDependent;
269 269
270 G4double crystalEnergyRef;//This 511 keV 270 G4double crystalEnergyRef;//This 511 keV
271 G4double crystalQuantumEfficiency;// 271 G4double crystalQuantumEfficiency;//
272 G4double edep_AfterCrystalBlurring; 272 G4double edep_AfterCrystalBlurring;
273 G4double crystalCoeff; 273 G4double crystalCoeff;
274 G4double sigma_energyResolution; 274 G4double sigma_energyResolution;
275 275
276 G4double totalTime; 276 G4double totalTime;
277 G4double prev_totalTime; 277 G4double prev_totalTime;
278 G4double timeInterval; 278 G4double timeInterval;
279 G4double time_annihil; 279 G4double time_annihil;
280 G4double time_tof; 280 G4double time_tof;
281 G4double block_DeadTime; 281 G4double block_DeadTime;
282 G4double module_DeadTime; 282 G4double module_DeadTime;
283 283
284 G4double *blockTime; 284 G4double *blockTime;
285 G4double *moduleTime; 285 G4double *moduleTime;
286 286
287 // 287 //
288 G4double activityNow; 288 G4double activityNow;
289 G4double InitialActivity; 289 G4double InitialActivity;
290 G4double halfLife; 290 G4double halfLife;
291 291
292 G4String simulationType; 292 G4String simulationType;
293 293
294 // 294 //
295 //for output file to write results 295 //for output file to write results
296 std::ofstream ofs; 296 std::ofstream ofs;
297 G4String asciiFileName; 297 G4String asciiFileName;
298 G4String rootFileName; 298 G4String rootFileName;
299 299
300 // #ifdef USEROOT 300 // #ifdef USEROOT
301 // TFile* file; 301 // TFile* file;
302 // TTree* tSingles; 302 // TTree* tSingles;
303 // TTree* tCoincidence; 303 // TTree* tCoincidence;
304 // //TH1F*hb; 304 // //TH1F*hb;
305 // #endif 305 // #endif
306 306
307 //input file to read reflector pattern 307 //input file to read reflector pattern
308 std::ifstream ifs; 308 std::ifstream ifs;
309 }; 309 };
310 310
311 #endif 311 #endif
312 312
313 313