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