| Geant4 Cross Reference |
1 //***********************************************************************************************************
2 // BinToStd_proton_position.C
3 // Root command file
4 // Type: root BinToStd_proton_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 StimEvent_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 StimEvent
43 {
44 uint16_t energy_keV; // different from Pixe Event, it is in keV
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_proton_position()
107 {
108 // printf("%f %f %f\n", acos(1), acos(-1), acos(0));
109 // return;
110
111 //***********************************************************************
112 //**************************Detection parameters (begin)*****************
113 //***********************************************************************
114
115 const int nbProjection = 1;
116 const int nbSlice = 1;
117 const int nbPixel = 1;
118 double totalAngleSpan = 180.; // in degree
119
120 double angleOfDetector = 0.; // angle of detector relative to the incident
121
122 double distanceObjectDetector = 22000.; // um
123
124 // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex
125 // angle of the right circular cone in radian
126 double theta = 10.2 * TMath::DegToRad(); // in radian
127 double radiusOfDetector = distanceObjectDetector * tan(theta);
128 bool usePosition = true;
129
130 //***********************************************************************
131 //**************************Detection parameters (end)*******************
132 //***********************************************************************
133
134 FILE* input = fopen("../build/ProtonAtExit.dat", "rb");
135 FILE* out = fopen("../build/StimEvent_std", "wb");
136
137 if (input == NULL) {
138 printf("error for opening the input file\n");
139 return;
140 }
141
142 RunInfo runInfo;
143 StimEvent stimEvent;
144 Point centerOfDetector;
145 Point protonMomentum;
146 Point protonPosition;
147 Point intersectionPoint;
148 long long count = 0;
149 int runID = -1; // index of simulations, namely runID, starting from 0
150
151 // while(!feof(input)) //if not the end, read
152 while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
153 runID++;
154 int nbParticle = runInfo.nbParticle;
155
156 //(begin)****************************************************************
157 // the following codes are used only when in the simulation
158 // the index of projection, slice and pixel is not
159 // correctly configured
160 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
161 int remain = runID % (nbSlice * nbPixel);
162 runInfo.sliceIndex = remain / nbPixel;
163 runInfo.pixelIndex = remain % nbPixel;
164
165 //(end)*******************************************************************
166
167 //***********************************************************************
168 //**************************Print information (begin)********************
169 //***********************************************************************
170
171 printf(
172 "---------RunID=%d:\nProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,"
173 "nbParticle = %d\n",
174 runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
175
176 //***********************************************************************
177 //**************************Print information (end)**********************
178 //***********************************************************************
179
180 if (!nbParticle) continue;
181 std::vector<ParticleInfo> protonAtExit(nbParticle);
182 fread(&protonAtExit[0], sizeof(ParticleInfo), nbParticle, input);
183
184 // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means
185 // the angle between source direction and detector, which should be constant
186 // when source is rotating
187 double ra = TMath::DegToRad()
188 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
189 centerOfDetector.m_x = distanceObjectDetector * cos(ra);
190 centerOfDetector.m_y = distanceObjectDetector * sin(ra);
191 centerOfDetector.m_z = 0;
192
193 for (int i = 0; i < nbParticle; ++i) {
194 // proton selection: energy should be lower than 4095 keV
195 if (protonAtExit[i].energy_keV >= 4095) continue;
196
197 protonMomentum.m_x = protonAtExit[i].mx;
198 protonMomentum.m_y = protonAtExit[i].my;
199 protonMomentum.m_z = protonAtExit[i].mz;
200
201 if (!usePosition) {
202 if (!IsDetected(centerOfDetector, protonMomentum, theta)) continue;
203 }
204 else {
205 double c =
206 distanceObjectDetector * (protonMomentum.m_x * cos(ra) + protonMomentum.m_y * sin(ra));
207 if (IsEqual(0, c, eps, releps)) continue; // parallel
208
209 protonPosition.m_x = protonAtExit[i].x;
210 protonPosition.m_y = protonAtExit[i].y;
211 protonPosition.m_z = protonAtExit[i].z;
212
213 double t = (distanceObjectDetector * distanceObjectDetector
214 - protonPosition.m_x * distanceObjectDetector * cos(ra)
215 - protonPosition.m_y * distanceObjectDetector * sin(ra))
216 / c;
217
218 intersectionPoint.m_x = protonPosition.m_x + protonMomentum.m_x * t;
219 intersectionPoint.m_y = protonPosition.m_y + protonMomentum.m_y * t;
220 intersectionPoint.m_z = protonPosition.m_z + protonMomentum.m_z * t;
221
222 if (!IsDetected_position(centerOfDetector, intersectionPoint, radiusOfDetector)) continue;
223
224 // printf(" t = %f, intersection point: (%f, %f, %f) centor of detector: (%f, %f, %f)
225 // 111=%f, 222=%f \n", t, intersectionPoint.m_x,intersectionPoint.m_y,intersectionPoint.m_z,
226 // centerOfDetector.m_x,centerOfDetector.m_y,centerOfDetector.m_z,
227 // (distanceObjectDetector*distanceObjectDetector-protonPosition.m_x*distanceObjectDetector*cos(ra)
228 // -protonPosition.m_y*distanceObjectDetector*sin(ra)), c);
229
230 // printf(" distanceObjectDetector = %f, protonPosition.m_x=%f,
231 // distanceObjectDetector*cos(ra)=%f, protonPosition.m_y=%f,
232 // distanceObjectDetector*sin(ra)=%f\n", distanceObjectDetector, protonPosition.m_x,
233 // distanceObjectDetector*cos(ra),
234 // protonPosition.m_y,
235 // distanceObjectDetector*sin(ra));
236
237 double tt = (intersectionPoint.m_x - protonPosition.m_x) * protonMomentum.m_x
238 + (intersectionPoint.m_y - protonPosition.m_y) * protonMomentum.m_y
239 + (intersectionPoint.m_z - protonPosition.m_z) * protonMomentum.m_z;
240 if (tt < 0) continue;
241 }
242
243 stimEvent.energy_10eV = floor(100 * protonAtExit[i].energy_keV + 0.5);
244 stimEvent.projectionIndex = runInfo.projectionIndex;
245 stimEvent.sliceIndex = runInfo.sliceIndex;
246 stimEvent.pixelIndex = runInfo.pixelIndex;
247 fwrite(&stimEvent, 7, 1, out);
248 count++;
249
250 //***********************************************************************
251 //**************************Print information
252 //(begin)********************
253 //***********************************************************************
254
255 if (!usePosition) {
256 printf(
257 "---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, momentum: "
258 "(%f, %f, %f), energy: %f keV\n",
259 i, runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex,
260 protonAtExit[i].mx, protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV);
261 }
262 else {
263 // printf("---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,
264 // momentum: (%f, %f, %f), energy: %f keV, position: (%f, %f, %f)\n", i, runID,
265 // runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, protonAtExit[i].mx,
266 // protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV, protonAtExit[i].x,
267 // protonAtExit[i].y, protonAtExit[i].z);
268
269 printf(
270 "---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, momentum: "
271 "(%f, %f, %f), energy: %f keV\n",
272 i, runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex,
273 protonAtExit[i].mx, protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV);
274 }
275
276 //***********************************************************************
277 //**************************Print information (end)**********************
278 //***********************************************************************
279 }
280 }
281 printf(
282 "\n---------------Number of StimEvent in total: "
283 "%lld------------------------\n",
284 count);
285 fclose(input);
286 fclose(out);
287
288 // FILE* input2;
289 // input2 = fopen("StimEvent_std.DAT","rb");
290 // StimEvent p;
291 // double eventId = -1;
292 // while(fread(&p, 7, 1, input2))
293 // {
294
295 // if(p.projectionIndex == 8 &&p.sliceIndex ==64 && p.pixelIndex==64)
296 // {
297 // eventId++;
298 // printf("StimEvent_%.0f ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_keV=%d keV\n",
299 // eventId, p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_keV);
300
301 // }
302
303 // }
304 // fclose(input2);
305 }
306