| Geant4 Cross Reference |
1 //***********************************************************************************************************
2 // BinToStd_GammaAtExit.C
3 // Root command file
4 // Type: root BinToStd_GammaAtExit.C
5 //
6 // Read the output file ProtonAtExit.dat that is generated by Geant4 tomography
7 // simulation It read all the gamma at exit information, and rewrite the events
8 // in a binary file PixeEvent_std_AtExit.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 // Define a structure to read and write each event in the required binary format
23 struct PixeEvent
24 {
25 uint16_t energy_10eV;
26 uint16_t pixelIndex;
27 uint16_t sliceIndex;
28 uint8_t projectionIndex;
29 };
30
31 struct ParticleInfo
32 {
33 float energy_keV;
34 float mx;
35 float my;
36 float mz;
37 };
38
39 struct RunInfo
40 {
41 // uint_16t
42 uint8_t projectionIndex; // 1 byte
43 uint16_t sliceIndex; //
44 uint16_t pixelIndex;
45 uint32_t nbParticle; // 4 bytes int
46 };
47 struct Point
48 {
49 double m_x;
50 double m_y;
51 double m_z;
52 };
53
54 bool IsDetected(Point poi1, Point poi2, double theta)
55 {
56 double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
57 / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
58 / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
59 if (a > 1.0) a = 1;
60 if (a < -1.0) a = -1;
61 double r = acos(a);
62 if (r > theta)
63 return false;
64 else
65 return true;
66 }
67 void BinToStd_GammaAtExit()
68 {
69 //***********************************************************************
70 //**************************Detection parameters (begin)*****************
71 //***********************************************************************
72
73 const int nbProjection = 10;
74 const int nbSlice = 1;
75 const int nbPixel = 20;
76 double totalAngleSpan = 180.; // in degree
77
78 double angleOfDetector = 135.; // angle of detector relative to the incident
79 // direction of the primary protons //
80 double distanceObjectDetector = 22.; // 22 mm
81 double radiusOfDetector = 5.; // 5 mm
82 // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex
83 // angle of the right circular cone in radian
84 double theta = 70 * TMath::DegToRad(); // in radian
85
86 //***********************************************************************
87 //**************************Detection parameters (end)*******************
88 //***********************************************************************
89
90 FILE* input = fopen("../build/GammaAtExit.dat", "rb");
91 FILE* out = fopen("../build/PixeEvent_std_AtExit.DAT", "wb");
92
93 if (input == NULL) {
94 printf("error for opening the input GammaAtExit.dat file\n");
95 return;
96 }
97
98 RunInfo runInfo;
99 PixeEvent pixeEvent;
100 Point centerOfDetector;
101 Point gammaMomentum;
102 long long count = 0;
103 int runID = -1; // index of simulations, namely runID, starting from 0
104
105 // while(!feof(input)) //if not the end, read
106 while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
107 runID++;
108 int nbParticle = runInfo.nbParticle;
109
110 // the following codes are used only when in the simulation
111 // ************(begin) the index of projection, slice and pixel is not
112 // correctly configured
113 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
114 int remain = runID % (nbSlice * nbPixel);
115 runInfo.sliceIndex = remain / nbPixel;
116 runInfo.pixelIndex = remain % nbPixel;
117 //************************************************************************(end)
118
119 //***********************************************************************
120 //**************************Print information (begin)********************
121 //***********************************************************************
122
123 printf(
124 "---------RunID=%d:\nProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,"
125 "nbParticle = %d\n",
126 runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
127
128 //***********************************************************************
129 //**************************Print information (end)**********************
130 //***********************************************************************
131
132 if (!nbParticle) continue;
133 std::vector<ParticleInfo> gammaAtExit(nbParticle);
134 fread(&gammaAtExit[0], sizeof(ParticleInfo), nbParticle, input);
135
136 // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means
137 // the angle between source direction and detector, which should be constant
138 // when source is rotating
139 double ra = TMath::DegToRad()
140 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
141 centerOfDetector.m_x = distanceObjectDetector * cos(ra);
142 centerOfDetector.m_y = distanceObjectDetector * sin(ra);
143 centerOfDetector.m_z = 0;
144
145 for (int i = 0; i < nbParticle; ++i) {
146 // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
147 if (gammaAtExit[i].energy_keV >= 40.95 || gammaAtExit[i].energy_keV <= 0.9) continue;
148
149 gammaMomentum.m_x = gammaAtExit[i].mx;
150 gammaMomentum.m_y = gammaAtExit[i].my;
151 gammaMomentum.m_z = gammaAtExit[i].mz;
152
153 if (!IsDetected(centerOfDetector, gammaMomentum, theta))
154 continue;
155 else {
156 pixeEvent.energy_10eV = floor(100 * gammaAtExit[i].energy_keV + 0.5);
157 pixeEvent.projectionIndex = runInfo.projectionIndex;
158 pixeEvent.sliceIndex = runInfo.sliceIndex;
159 pixeEvent.pixelIndex = runInfo.pixelIndex;
160 fwrite(&pixeEvent, 7, 1, out);
161 count++;
162
163 //***********************************************************************
164 //**************************Print information (begin)********************
165 //***********************************************************************
166
167 // printf("momentum: (%f, %f, %f), energy: %f keV %d 10eV\n",
168 // gammaAtExit[i].mx, gammaAtExit[i].my, gammaAtExit[i].mz,
169 // gammaAtExit[i].energy_keV, pixeEvent.energy_10eV);
170
171 //***********************************************************************
172 //**************************Print information (end)**********************
173 //***********************************************************************
174 }
175 }
176 }
177 printf(
178 "\n---------------Number of PixeEvent in total: "
179 "%lld------------------------\n",
180 count);
181 fclose(input);
182 fclose(out);
183
184 // Recheck the output file in case
185 // FILE* input2;
186 // input2 = fopen("PixeEvent_std_AtExit.DAT","rb");
187 // PixeEvent p;
188 // while(fread(&p, 7, 1, input2))
189 // {
190 // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,
191 // Energy_10eV=%d\n", p.projectionIndex, p.sliceIndex, p.pixelIndex,
192 // p.energy_10eV);
193
194 // }
195 // fclose(input2);
196 }
197