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