| Geant4 Cross Reference |
1 //***********************************************************************************************************
2 // Concatenate_BinToStd_GammaAtExit.C
3 // Root command file
4 // Type: root Concatenate_BinToStd_GammaAtExit.C
5 //
6 // It is used in case of one interruption
7 // Read 2 output files GammaAtExit_1.dat and GammaAtExit_2.dat that are generated by Geant4
8 // tomography simulation It reads gamma at exit information, and rewrite the events in a binary file
9 // PixeEvent_std_AtExit.DAT
10 //
11 // More information is available in UserGuide
12 // Created by Z.LI LP2i Bordeaux 2022
13 //***********************************************************************************************************
14
15 #include <math.h>
16 #include <stdint.h>
17 #include <stdio.h>
18 #include <string.h>
19
20 #include <vector>
21 // using namespace std;
22
23 // Define a structure to read and write each event in the required binary format
24 struct PixeEvent
25 {
26 uint16_t energy_10eV;
27 uint16_t pixelIndex;
28 uint16_t sliceIndex;
29 uint8_t projectionIndex;
30 };
31 struct ParticleInfo
32 {
33 float energy_keV;
34 float mx;
35 float my;
36 float mz;
37 };
38 struct RunInfo
39 {
40 // uint_16t
41 uint8_t projectionIndex; // 1 byte
42 uint16_t sliceIndex; //
43 uint16_t pixelIndex;
44 uint32_t nbParticle; // 4 bytes int
45 };
46
47 struct Point
48 {
49 double m_x;
50 double m_y;
51 double m_z;
52 };
53 bool IsDetected(Point poi1, Point poi2, double theta)
54 {
55 double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
56 / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
57 / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
58 if (a > 1.0) a = 1;
59 if (a < -1.0) a = -1;
60 double r = acos(a);
61 if (r > theta)
62 return false;
63 else
64 return true;
65 }
66 void Concatenate_BinToStd_GammaAtExit()
67 {
68 //***********************************************************************
69 //**************************Detection parameters (begin)*****************
70 //***********************************************************************
71
72 const int nbProjection = 10;
73 const int nbSlice = 1;
74 const int nbPixel = 20;
75 double totalAngleSpan = 180.; // in degree
76
77 double angleOfDetector =
78 135.; // angle of detector relative to the incident direction of the primary protons //
79 double distanceObjectDetector = 22.; // 22 mm
80 double radiusOfDetector = 5.; // 5 mm
81 // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex angle of the right
82 // circular cone in radian
83 double theta = 70 * TMath::DegToRad(); // in radian
84
85 int P_interrupt = 6; // Projection of interruption
86
87 //***********************************************************************
88 //**************************Detection parameters (end)*******************
89 //***********************************************************************
90
91 // assuming there is one interruption
92 FILE* input1 = fopen("../build/GammaAtExit_1.dat", "rb");
93 FILE* input2 = fopen("../build/GammaAtExit_2.dat", "rb");
94 FILE* out = fopen("../build/PixeEvent_std_AtExit.DAT", "wb");
95
96 if (input1 == NULL) {
97 printf("error for opening the input GammaAtExit_1.dat file\n");
98 return;
99 }
100 if (input2 == NULL) {
101 printf("error for opening the input GammaAtExit_2.dat file\n");
102 return;
103 }
104
105 RunInfo runInfo;
106 PixeEvent pixeEvent;
107 Point centerOfDetector;
108 Point gammaMomentum;
109 long long count1 = 0;
110 long long count2 = 0;
111 int runID = -1; // index of simulations, namely runID, starting from 0
112
113 // ************************************************************
114 // **********************READ FIRST FILE (begin)***************
115 // ************************************************************
116 while (fread(&runInfo, sizeof(RunInfo), 1, input1)) {
117 runID++;
118
119 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
120 int remain = runID % (nbSlice * nbPixel);
121 runInfo.sliceIndex = remain / nbPixel;
122 runInfo.pixelIndex = remain % nbPixel;
123
124 if (runInfo.projectionIndex == P_interrupt) {
125 runID--;
126 break;
127 }
128 int nbParticle = runInfo.nbParticle;
129
130 //***********************************************************************
131 //**************************Print information (begin)********************
132 //***********************************************************************
133
134 printf("-1--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
135 runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
136
137 //***********************************************************************
138 //**************************Print information (end)**********************
139 //***********************************************************************
140
141 if (!nbParticle) continue;
142
143 std::vector<ParticleInfo> gammaAtExit(nbParticle);
144 fread(&gammaAtExit[0], sizeof(ParticleInfo), nbParticle, input1);
145
146 // if(runInfo.sliceIndex!=1) continue;
147 // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
148 // if(runInfo.sliceIndex!=31) continue;
149
150 // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
151 // source direction and detector, which should be constant when source is rotating
152 double ra = TMath::DegToRad()
153 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
154 centerOfDetector.m_x = distanceObjectDetector * cos(ra);
155 centerOfDetector.m_y = distanceObjectDetector * sin(ra);
156 centerOfDetector.m_z = 0;
157
158 for (int i = 0; i < nbParticle; ++i) {
159 // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
160 if (gammaAtExit[i].energy_keV >= 40.95 || gammaAtExit[i].energy_keV <= 0.9)
161 continue; // gamma selection
162
163 gammaMomentum.m_x = gammaAtExit[i].mx;
164 gammaMomentum.m_y = gammaAtExit[i].my;
165 gammaMomentum.m_z = gammaAtExit[i].mz;
166
167 if (!IsDetected(centerOfDetector, gammaMomentum, theta))
168 continue;
169 else {
170 pixeEvent.energy_10eV = floor(100 * gammaAtExit[i].energy_keV + 0.5);
171 pixeEvent.projectionIndex = runInfo.projectionIndex;
172 pixeEvent.sliceIndex = runInfo.sliceIndex;
173 pixeEvent.pixelIndex = runInfo.pixelIndex;
174 fwrite(&pixeEvent, 7, 1, out);
175 count1++;
176 }
177 }
178 }
179 printf("---------------Number of PixeEvent in the first file: %lld------------------------\n",
180 count1);
181 fclose(input1);
182
183 // ************************************************************
184 // **********************READ FIRST FILE (end)*****************
185 // ************************************************************
186
187 // ************************************************************
188 // **********************READ SECOND FILE (begin)**************
189 // ************************************************************
190 while (fread(&runInfo, sizeof(RunInfo), 1, input2)) {
191 runID++;
192
193 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
194 int remain = runID % (nbSlice * nbPixel);
195 runInfo.sliceIndex = remain / nbPixel;
196 runInfo.pixelIndex = remain % nbPixel;
197
198 int nbParticle = runInfo.nbParticle;
199
200 //***********************************************************************
201 //**************************Print information (begin)********************
202 //***********************************************************************
203
204 printf("-2--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
205 runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
206
207 //***********************************************************************
208 //**************************Print information (end)**********************
209 //***********************************************************************
210
211 if (!nbParticle) continue;
212 std::vector<ParticleInfo> gammaAtExit(nbParticle);
213 fread(&gammaAtExit[0], sizeof(ParticleInfo), nbParticle, input2);
214
215 // if(runInfo.sliceIndex!=1) continue;
216 // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
217 // if(runInfo.sliceIndex!=31) continue;
218
219 // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
220 // source direction and detector, which should be constant when source is rotating
221 double ra = TMath::DegToRad()
222 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
223 centerOfDetector.m_x = distanceObjectDetector * cos(ra);
224 centerOfDetector.m_y = distanceObjectDetector * sin(ra);
225 centerOfDetector.m_z = 0;
226
227 for (int i = 0; i < nbParticle; ++i) {
228 // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
229 if (gammaAtExit[i].energy_keV >= 40.95 || gammaAtExit[i].energy_keV <= 0.9)
230 continue; // gamma selection
231
232 gammaMomentum.m_x = gammaAtExit[i].mx;
233 gammaMomentum.m_y = gammaAtExit[i].my;
234 gammaMomentum.m_z = gammaAtExit[i].mz;
235
236 if (!IsDetected(centerOfDetector, gammaMomentum, theta))
237 continue;
238 else {
239 pixeEvent.energy_10eV = floor(100 * gammaAtExit[i].energy_keV + 0.5);
240 pixeEvent.projectionIndex = runInfo.projectionIndex;
241 pixeEvent.sliceIndex = runInfo.sliceIndex;
242 pixeEvent.pixelIndex = runInfo.pixelIndex;
243 fwrite(&pixeEvent, 7, 1, out);
244 count2++;
245 }
246 }
247 }
248 printf("---------------Number of PixeEvent in in the second file: %lld------------------------\n",
249 count2);
250
251 // ************************************************************
252 // **********************READ SECOND FILE (end)****************
253 // ************************************************************
254
255 printf("---------------Number of PixeEvent in total: %lld------------------------\n",
256 count1 + count2);
257 fclose(input2);
258 fclose(out);
259
260 // Recheck the output file in case
261 // FILE* input2 = fopen("PixeEvent_std_AtExit.DAT","rb");
262 // PixeEvent p;
263 // while(fread(&p, 7, 1, input2))
264 // {
265 // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_10eV=%d\n",
266 // p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_10eV);
267
268 // }
269 // fclose(input2);
270 }
271