| Geant4 Cross Reference |
1 //***********************************************************************************************************
2 // Concatenate_BinToStd_GammaAtCreation.C
3 // Root command file
4 // Type: root Concatenate_BinToStd_GammaAtCreation.C
5 //
6 // It is used in case of interruption
7 // Read 2 output files GammaAtCreation_1.dat and GammaAtCreation_2.dat that are generated by Geant4
8 // tomography simulation. It reads all the gamma at creation information, and rewrite the events in
9 // a binary file PixeEvent_std_AtCreation.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_GammaAtCreation()
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/GammaAtCreation_1.dat", "rb");
93 FILE* input2 = fopen("../build/GammaAtCreation_2.dat", "rb");
94 FILE* out = fopen("../build/PixeEvent_std_AtCreation.DAT", "wb");
95
96 if (input1 == NULL) {
97 printf("error for opening the input GammaAtCreation_1.dat file\n");
98 return;
99 }
100 if (input2 == NULL) {
101 printf("error for opening the input GammaAtCreation_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 // ************************************************************(begin)
114 // **********************READ FIRST FILE***********************
115 // ************************************************************
116 while (fread(&runInfo, sizeof(RunInfo), 1, input1)) {
117 runID++;
118
119 //(begin)***************************************************************
120 // the following codes are used only when in the simulation
121 // the index of projection, slice and pixel is not
122 // correctly configured
123 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
124 int remain = runID % (nbSlice * nbPixel);
125 runInfo.sliceIndex = remain / nbPixel;
126 runInfo.pixelIndex = remain % nbPixel;
127 //(end)******************************************************************
128
129 if (runInfo.projectionIndex == P_interrupt) {
130 runID--;
131 break;
132 }
133
134 int nbParticle = runInfo.nbParticle;
135
136 std::vector<ParticleInfo> gammaAtCreation(nbParticle);
137 fread(&gammaAtCreation[0], sizeof(ParticleInfo), nbParticle, input1);
138
139 // if(runInfo.sliceIndex!=1) continue;
140 // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
141 // if(runInfo.sliceIndex!=31) continue;
142
143 //***********************************************************************
144 //**************************Print information (begin)********************
145 //***********************************************************************
146
147 printf("-1--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
148 runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
149
150 //***********************************************************************
151 //**************************Print information (end)**********************
152 //***********************************************************************
153
154 // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
155 // source direction and detector, which should be constant when source is rotating
156 double ra = TMath::DegToRad()
157 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
158 centerOfDetector.m_x = distanceObjectDetector * cos(ra);
159 centerOfDetector.m_y = distanceObjectDetector * sin(ra);
160 centerOfDetector.m_z = 0;
161
162 for (int i = 0; i < nbParticle; ++i) {
163 // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
164 if (gammaAtCreation[i].energy_keV >= 40.95 || gammaAtCreation[i].energy_keV <= 0.9)
165 continue; // gamma selection
166
167 gammaMomentum.m_x = gammaAtCreation[i].mx;
168 gammaMomentum.m_y = gammaAtCreation[i].my;
169 gammaMomentum.m_z = gammaAtCreation[i].mz;
170
171 if (!IsDetected(centerOfDetector, gammaMomentum, theta))
172 continue;
173 else {
174 pixeEvent.energy_10eV = floor(100 * gammaAtCreation[i].energy_keV + 0.5);
175 pixeEvent.projectionIndex = runInfo.projectionIndex;
176 pixeEvent.sliceIndex = runInfo.sliceIndex;
177 pixeEvent.pixelIndex = runInfo.pixelIndex;
178 fwrite(&pixeEvent, 7, 1, out);
179 count1++;
180 }
181 }
182 }
183 printf("---------------Number of PixeEvent in the first file: %lld------------------------\n",
184 count1);
185 fclose(input1);
186
187 // ************************************************************
188 // **********************READ FIRST FILE (end)*****************
189 // ************************************************************
190
191 // ************************************************************
192 // **********************READ SECOND FILE (begin)**************
193 // ************************************************************
194 while (fread(&runInfo, sizeof(RunInfo), 1, input2)) {
195 runID++;
196
197 //(begin)***************************************************************
198 // the following codes are used only when in the simulation
199 // the index of projection, slice and pixel is not
200 // correctly configured
201 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
202 int remain = runID % (nbSlice * nbPixel);
203 runInfo.sliceIndex = remain / nbPixel;
204 runInfo.pixelIndex = remain % nbPixel;
205 //(end)******************************************************************
206
207 int nbParticle = runInfo.nbParticle;
208
209 //***********************************************************************
210 //**************************Print information (begin)********************
211 //***********************************************************************
212
213 printf("-2--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
214 runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
215
216 //***********************************************************************
217 //**************************Print information (end)**********************
218 //***********************************************************************
219
220 if (!nbParticle) continue;
221 std::vector<ParticleInfo> gammaAtCreation(nbParticle);
222 fread(&gammaAtCreation[0], sizeof(ParticleInfo), nbParticle, input2);
223
224 // if(runInfo.sliceIndex!=1) continue;
225 // if(runInfo.sliceIndex!=31) continue;
226 // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
227
228 // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
229 // source direction and detector, which should be constant when source is rotating
230 double ra = TMath::DegToRad()
231 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
232 centerOfDetector.m_x = distanceObjectDetector * cos(ra);
233 centerOfDetector.m_y = distanceObjectDetector * sin(ra);
234 centerOfDetector.m_z = 0;
235
236 for (int i = 0; i < nbParticle; ++i) {
237 // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
238 if (gammaAtCreation[i].energy_keV >= 40.95 || gammaAtCreation[i].energy_keV <= 0.9)
239 continue; // gamma selection
240
241 gammaMomentum.m_x = gammaAtCreation[i].mx;
242 gammaMomentum.m_y = gammaAtCreation[i].my;
243 gammaMomentum.m_z = gammaAtCreation[i].mz;
244
245 if (!IsDetected(centerOfDetector, gammaMomentum, theta))
246 continue;
247 else {
248 pixeEvent.energy_10eV = floor(100 * gammaAtCreation[i].energy_keV + 0.5);
249 pixeEvent.projectionIndex = runInfo.projectionIndex;
250 pixeEvent.sliceIndex = runInfo.sliceIndex;
251 pixeEvent.pixelIndex = runInfo.pixelIndex;
252 fwrite(&pixeEvent, 7, 1, out);
253 count2++;
254 }
255 }
256 }
257 printf("---------------Number of PixeEvent in in the second file: %lld------------------------\n",
258 count2);
259
260 // ************************************************************
261 // **********************READ SECOND FILE (end)****************
262 // ************************************************************
263
264 printf("---------------Number of PixeEvent in total: %lld------------------------\n",
265 count1 + count2);
266 fclose(input2);
267 fclose(out);
268
269 // Recheck the output file in case
270 // FILE* input2 = fopen("PixeEvent_std_AtCreation.DAT","rb");
271 // PixeEvent p;
272 // while(fread(&p, 7, 1, input2))
273 // {
274 // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_10eV=%d\n",
275 // p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_10eV);
276
277 // }
278 // fclose(input2);
279 }
280