| Geant4 Cross Reference |
1 //***********************************************************************************************************
2 // Concatenate_BinToStd_ProtonAtExit.C
3 // Root command file
4 // Type: root Concatenate_BinToStd_ProtonAtExit.C
5 //
6 // It is used in case of interruption
7 // Read 2 output files ProtonAtExit_1.dat and ProtonAtExit_2.dat that are generated by Geant4
8 // tomography simulation It reads protons at exit information, and rewrite the events in a binary
9 // file StimEvent_std.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 StimEvent
25 {
26 uint16_t energy_keV; // different from Pixe Event, it is in keV
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
67 void Recheck()
68 {
69 // Recheck the output file in case
70 FILE* input3 = fopen("../build/StimEvent_std_Detector0_Aperture10.2.DAT", "rb");
71 StimEvent p;
72 double eventId = -1;
73 while (fread(&p, 7, 1, input3)) {
74 if (p.projectionIndex == 8 && p.sliceIndex == 64 && p.pixelIndex == 10) {
75 eventId++;
76 printf("StimEvent_%.0f ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_keV=%d keV\n",
77 eventId, p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_keV);
78 }
79 }
80 fclose(input3);
81 }
82 void Concatenate_BinToStd_ProtonAtExit()
83 {
84 // Recheck();
85 // return;
86
87 //***********************************************************************
88 //**************************Detection parameters (begin)*****************
89 //***********************************************************************
90
91 const int nbProjection = 10;
92 const int nbSlice = 128;
93 const int nbPixel = 20;
94 double totalAngleSpan = 180.; // in degree
95
96 // angle of detector relative to the incident direction of the primary protons at first projection
97 // for proton, it is fixed to 0 degree, namely opposite to the source
98 double angleOfDetector = 0.;
99 double distanceObjectDetector = 22.; // 22 mm
100 double radiusOfDetector = 5.; // 5 mm
101 // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex angle of the right
102 // circular cone in radian
103 double theta = 10.2 * TMath::DegToRad(); // in radian
104
105 int P_interrupt = 2; // Projection of interruption
106
107 //***********************************************************************
108 //**************************Detection parameters (end)*******************
109 //***********************************************************************
110
111 // assuming there is one interruption
112 FILE* input1 = fopen("../build/ProtonAtExit_1.dat", "rb");
113 FILE* input2 = fopen("../build/ProtonAtExit_2.dat", "rb");
114 FILE* out = fopen("../build/StimEvent_std.DAT", "wb");
115
116 if (input1 == NULL) {
117 printf("error for opening the input ProtonAtExit_1.dat file\n");
118 return;
119 }
120 if (input2 == NULL) {
121 printf("error for opening the input ProtonAtExit_2.dat file\n");
122 return;
123 }
124
125 RunInfo runInfo;
126 StimEvent stimEvent;
127 Point centerOfDetector;
128 Point protonMomentum;
129
130 long long count1 = 0;
131 long long count2 = 0;
132 int runID = -1; // index of simulations, namely runID, starting from 0
133
134 // ************************************************************(begin)
135 // **********************READ FIRST FILE***********************
136 // ************************************************************
137 while (fread(&runInfo, sizeof(RunInfo), 1, input1)) {
138 runID++;
139 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
140 int remain = runID % (nbSlice * nbPixel);
141 runInfo.sliceIndex = remain / nbPixel;
142 runInfo.pixelIndex = remain % nbPixel;
143 if (runInfo.projectionIndex == P_interrupt) {
144 runID--;
145 break;
146 }
147
148 int nbParticle = runInfo.nbParticle;
149
150 //***********************************************************************
151 //**************************Print information (begin)********************
152 //***********************************************************************
153
154 printf("-1--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
155 runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
156
157 //***********************************************************************
158 //**************************Print information (end)**********************
159 //***********************************************************************
160
161 if (!nbParticle) continue;
162 std::vector<ParticleInfo> protonAtExit(nbParticle);
163 fread(&protonAtExit[0], sizeof(ParticleInfo), nbParticle, input1);
164
165 // if(runInfo.sliceIndex!=1) continue;
166 // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
167 // if(runInfo.sliceIndex!=31) continue;
168
169 // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
170 // source direction and detector, which should be constant when source is rotating
171 double ra = TMath::DegToRad()
172 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
173 centerOfDetector.m_x = distanceObjectDetector * cos(ra);
174 centerOfDetector.m_y = distanceObjectDetector * sin(ra);
175 centerOfDetector.m_z = 0;
176
177 for (int i = 0; i < nbParticle; ++i) {
178 // proton selection: energy should be lower than 4095 keV
179 if (protonAtExit[i].energy_keV >= 4095) continue; // proton selection
180
181 protonMomentum.m_x = protonAtExit[i].mx;
182 protonMomentum.m_y = protonAtExit[i].my;
183 protonMomentum.m_z = protonAtExit[i].mz;
184
185 if (!IsDetected(centerOfDetector, protonMomentum, theta))
186 continue;
187 else {
188 stimEvent.energy_keV = floor(protonAtExit[i].energy_keV + 0.5);
189 stimEvent.projectionIndex = runInfo.projectionIndex;
190 stimEvent.sliceIndex = runInfo.sliceIndex;
191 stimEvent.pixelIndex = runInfo.pixelIndex;
192 fwrite(&stimEvent, 7, 1, out);
193 count1++;
194 }
195 }
196 }
197 printf("---------------Number of StimEvent in the first file: %lld------------------------\n",
198 count1);
199 fclose(input1);
200
201 // ************************************************************
202 // **********************READ FIRST FILE (end)*****************
203 // ************************************************************
204
205 // ************************************************************
206 // **********************READ SECOND FILE (begin)**************
207 // ************************************************************
208 while (fread(&runInfo, sizeof(RunInfo), 1, input2)) {
209 runID++;
210 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
211 int remain = runID % (nbSlice * nbPixel);
212 runInfo.sliceIndex = remain / nbPixel;
213 runInfo.pixelIndex = remain % nbPixel;
214
215 int nbParticle = runInfo.nbParticle;
216
217 //***********************************************************************
218 //**************************Print information (begin)********************
219 //***********************************************************************
220
221 printf("-2--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
222 runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
223
224 //***********************************************************************
225 //**************************Print information (end)**********************
226 //***********************************************************************
227
228 if (!nbParticle) continue;
229
230 std::vector<ParticleInfo> protonAtExit(nbParticle);
231 fread(&protonAtExit[0], sizeof(ParticleInfo), nbParticle, input2);
232
233 // if(runInfo.sliceIndex!=1) continue;
234 // if(runInfo.sliceIndex!=31) continue;
235 // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
236
237 // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
238 // source direction and detector, which should be constant when source is rotating
239 double ra = TMath::DegToRad()
240 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
241 centerOfDetector.m_x = distanceObjectDetector * cos(ra);
242 centerOfDetector.m_y = distanceObjectDetector * sin(ra);
243 centerOfDetector.m_z = 0;
244
245 for (int i = 0; i < nbParticle; ++i) {
246 // proton selection: energy should be lower than 4095 keV
247 if (protonAtExit[i].energy_keV >= 4095) continue; // proton selection
248
249 protonMomentum.m_x = protonAtExit[i].mx;
250 protonMomentum.m_y = protonAtExit[i].my;
251 protonMomentum.m_z = protonAtExit[i].mz;
252
253 if (!IsDetected(centerOfDetector, protonMomentum, theta))
254 continue;
255 else {
256 stimEvent.energy_keV = floor(protonAtExit[i].energy_keV + 0.5);
257 stimEvent.projectionIndex = runInfo.projectionIndex;
258 stimEvent.sliceIndex = runInfo.sliceIndex;
259 stimEvent.pixelIndex = runInfo.pixelIndex;
260 fwrite(&stimEvent, 7, 1, out);
261 count2++;
262 }
263 }
264 }
265 printf("---------------Number of StimEvent in in the second file: %lld------------------------\n",
266 count2);
267
268 // ************************************************************
269 // **********************READ SECOND FILE (end)****************
270 // ************************************************************
271
272 printf("---------------Number of StimEvent in total: %lld------------------------\n",
273 count1 + count2);
274 fclose(input2);
275 fclose(out);
276 }
277