| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25
26 //GEANT4 - Depth-of-Interaction enabled Positron emission tomography (PET) advanced example
27
28 //Authors and contributors
29
30 // Author list to be updated, with names of co-authors and contributors from National Institute of Radiological Sciences (NIRS)
31
32 // Abdella M. Ahmed (1, 2), Andrew Chacon (1, 2), Harley Rutherford (1, 2),
33 // Hideaki Tashima (3), Go Akamatsu (3), Akram Mohammadi (3), Eiji Yoshida (3), Taiga Yamaya (3)
34 // Susanna Guatelli (2), and Mitra Safavi-Naeini (1, 2)
35
36 // (1) Australian Nuclear Science and Technology Organisation, Australia
37 // (2) University of Wollongong, Australia
38 // (3) National Institute of Radiological Sciences, Japan
39
40
41 #ifndef doiPETAnalysis_h
42 #define doiPETAnalysis_h 1
43
44 #include "doiPETGlobalParameters.hh"
45 #include "globals.hh"
46 #include <vector>
47 #include <time.h>
48 #include <map>
49 #include <set>
50 #include "G4ThreeVector.hh"
51 #include <iostream>
52 #include <fstream>
53 #include <sstream>
54 #include <iterator>
55 #include <vector>
56 #include <algorithm>
57 #include "G4AnalysisManager.hh"
58
59 // Define the total number of columns in the ntuple
60 const G4int MaxNtCol = 17;
61
62 class doiPETAnalysisMessenger;
63
64 //class InteractionInformation;
65
66 class doiPETAnalysis
67 {
68 private:
69 doiPETAnalysis();
70
71 public:
72 ~doiPETAnalysis();
73 static doiPETAnalysis* GetInstance();
74 void FindInteractingCrystal();
75 void Open(G4String);
76 void Close();
77 void Delete();
78 void ResetNumberOfHits();
79 void Write(/*G4int, G4int, G4int, G4double*/);
80 void WriteOutput();
81
82 //void GetIntractionInfomation(InteractionInformation*);
83
84 void GetParentParticleName(G4String);
85 void GetSizeOfDetector (G4double, G4double, G4double);
86 void SetScatterIndexInPhantom(G4int);
87
88 void SetSourcePosition(G4ThreeVector);//
89 void SetEventID(G4int);
90
91 void BlurringParameters();
92 void GetTimeOfAnnihilation(G4double);
93
94
95 void PMTPosition();
96 void AngerLogic(G4double, G4double, G4double, G4double, G4double, G4bool);
97 void ReadReflectorPattern();
98 void PrepareDOILookUpTable(G4String);
99
100 void SetActivity(G4double);
101 void SetIsotopeHalfLife(G4double);
102 void CrystalIDAfterAngerLogic(G4int, G4int, G4int);
103 void TypeOfOutput(G4String);//Single or coincidence list-mode data
104 void CalulateAcquisitionTime();
105 //G4double QuantumEffifciency(G4double);
106 G4double QuantumEffifciency(G4double, G4int, G4int);
107 void ReadOut(G4int, G4int, G4double, G4double, G4ThreeVector, G4double);
108
109 //G4ROOT
110 void book(); // booking the ROOT file
111
112
113 void FillListModeEvent(); //Single or Coinsidence
114 void finish();
115 // Close the ROOT file with all the results stored in nutples
116
117
118 private:
119 static doiPETAnalysis* instance;
120 doiPETAnalysisMessenger* fAnalysisMessenger;
121 //std::multimap< G4int, InteractionInformation* > mapBlockInteraction;
122 std::set<G4int> setBlockInteraction;
123
124 G4double upperThreshold, lowerThreshold;
125 G4double triggerEnergy;
126
127 //G4ROOT
128 G4bool factoryOn;
129 G4int fNtColId[MaxNtCol];
130 //
131
132
133 //G4ThreeVector sourcePosition;
134
135
136
137 //
138 G4int scatterIndex;
139
140 G4String parentParticleName;//
141
142 //
143 G4int numberofInteractions;
144 G4int countCoincidence;
145
146 G4int numberOfBlocks_total;
147
148 G4double sizeOfDetector_DOI,sizeOfDetector_axial,sizeOfDetector_tangential;
149
150 //Virtual position of the PMT
151 G4double signalPMT1, signalPMT2, signalPMT3, signalPMT4;
152
153 G4double posPMT1x, posPMT2x, posPMT3x, posPMT4x;
154 G4double posPMT1y, posPMT2y, posPMT3y, posPMT4y;
155 G4double posPMT1z, posPMT2z, posPMT3z, posPMT4z;
156
157 //
158 G4double signalPMT1z, signalPMT2z, signalPMT3z, signalPMT4z;
159 G4double signalPMT1y, signalPMT2y, signalPMT3y, signalPMT4y;
160
161 //
162 G4double signalZplus, signalZminus;
163 G4double signalYplus, signalYminus;
164 //
165
166 G4double dist1z, dist2z, dist3z, dist4z, distz;
167 G4double dist1y, dist2y, dist3y, dist4y, disty;
168
169 G4double shiftCoeff;
170
171 G4double PositionAngerZ, PositionAngerY;
172
173 //reflector pattern
174 std::vector<G4int> ireflectorLayer1_Tangential;
175 std::vector<G4int> ireflectorLayer1_Axial;
176 std::vector<G4int> ireflectorLayer2_Tangential;
177 std::vector<G4int> ireflectorLayer2_Axial;
178 std::vector<G4int> ireflectorLayer3_Tangential;
179 std::vector<G4int> ireflectorLayer3_Axial;
180 std::vector<G4int> ireflectorLayer4_Tangential;
181 std::vector<G4int> ireflectorLayer4_Axial;
182 std::vector<G4int> doi_table;
183 //
184
185 //The number of pixes for the 2D position histogram after Anger Logic calculation
186 G4int numberOfPixel_axial;
187 G4int numberOfPixel_tan;
188
189 //source position
190 G4double spositionX;
191 G4double spositionY;
192 G4double spositionZ;
193
194 //interaction position with respect to the crystal axis
195 G4ThreeVector interactionPos;
196
197 //interaction position
198 G4double intPosX;
199 G4double intPosY;
200 G4double intPosZ;
201
202
203 G4double interactionTime;
204
205 G4int crystalID;//contineous crystal ID in 3D
206 G4int crystalID_2D;
207
208 G4int prev_eventID;
209
210 //Single output
211 G4int eventID;
212 G4int blockID;
213 G4int crystalID_axial;
214 G4int crystalID_tangential;
215 G4int DOI_ID;
216 G4double timeStamp;
217 G4double totalEdep;
218
219 //coincidence output
220 G4int eventID0, eventID1;
221 G4int blockID0, blockID1;
222 G4int crystalID_axial0, crystalID_axial1;
223 G4int crystalID_tangential0, crystalID_tangential1;
224 G4int DOI_ID0, DOI_ID1;
225 G4double timeStamp0, timeStamp1;
226 G4double totalEdep0, totalEdep1;
227 G4double sposX, sposY, sposZ;
228 //choice for the user
229 G4bool getSinglesData;
230 G4bool getCoincidenceData;
231
232 //
233 G4bool ApplyAngerLogic;
234
235 G4double PMTblurring_tan;
236 G4double PMTblurring_axial;
237
238 G4String outputData;
239 G4int numberOfHit;
240 std::vector<G4int> eventID_coin;
241 std::vector<G4double> edep_coin;
242 std::vector<G4int>blockID_coin;
243 std::vector<G4int> cryID_axial_coin;
244 std::vector<G4int> cryID_tan_coin;
245 std::vector<G4int> cryDOI_coin;
246 std::vector<G4double> time_coin;
247
248 //Crystal IDs after Anger Logic calculation
249 G4int crystalIDNew_DOI, crystalIDNew_tan, crystalIDNew_axial;
250
251 //Crystal ID in the 2D position histogram along the axial and tangetial direction
252 G4int crystalID_in2D_posHist_axial, crystalID_in2D_posHist_tan;
253
254 //continous crystal ID after after Anger Logic.
255 G4int crystalID_in2D_posHist;
256
257
258 //Crystal blurring
259 G4double crystalResolution;
260 G4double crystalResolutionMin;//
261 G4double crystalResolutionMax;//
262
263 //G4bool variableResolution;
264 G4bool fixedResolution;
265 G4bool isDOIlookUpTablePrepared;
266
267 G4double energyResolution_fixed;
268 std::vector<std::vector<G4double>> energyResolution_cryDependent;
269
270 G4double crystalEnergyRef;//This 511 keV
271 G4double crystalQuantumEfficiency;//
272 G4double edep_AfterCrystalBlurring;
273 G4double crystalCoeff;
274 G4double sigma_energyResolution;
275
276 G4double totalTime;
277 G4double prev_totalTime;
278 G4double timeInterval;
279 G4double time_annihil;
280 G4double time_tof;
281 G4double block_DeadTime;
282 G4double module_DeadTime;
283
284 G4double *blockTime;
285 G4double *moduleTime;
286
287 //
288 G4double activityNow;
289 G4double InitialActivity;
290 G4double halfLife;
291
292 G4String simulationType;
293
294 //
295 //for output file to write results
296 std::ofstream ofs;
297 G4String asciiFileName;
298 G4String rootFileName;
299
300 // #ifdef USEROOT
301 // TFile* file;
302 // TTree* tSingles;
303 // TTree* tCoincidence;
304 // //TH1F*hb;
305 // #endif
306
307 //input file to read reflector pattern
308 std::ifstream ifs;
309 };
310
311 #endif
312
313