| Geant4 Cross Reference |
1 //***********************************************************************************************************
2 // Concatenate_BinToStd_GammaAtExit_fabricate.C
3 // Root command file
4 // Use it by typing in the command line of Root terminal: root
5 // Concatenate_BinToStd_GammaAtExit_fabricate.C
6 //
7 //
8 //
9 // More information is available in UserGuide
10 // Created by Z.LI LP2i Bordeaux 2022
11 //***********************************************************************************************************
12
13 #include <math.h>
14 #include <stdint.h>
15 #include <stdio.h>
16 #include <string.h>
17
18 #include <vector>
19 // using namespace std;
20
21 #define PI 3.14159265f
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 double DegreeToRadian(double degree)
48 {
49 return (PI * degree / 180.);
50 }
51
52 struct Point
53 {
54 double m_x;
55 double m_y;
56 double m_z;
57 };
58 bool IsDetected(Point poi1, Point poi2, double theta)
59 {
60 double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
61 / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
62 / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
63 if (a > 1.0) a = 1;
64 if (a < -1.0) a = -1;
65 double r = acos(a);
66 if (r > theta)
67 return false;
68 else
69 return true;
70 }
71
72 void Concatenate_BinToStd_GammaAtExit_fabricate()
73 {
74 //***********************************************************************
75 //**************************Detection parameters (begin)*****************
76 //***********************************************************************
77
78 const int nbProjection = 100;
79 const int nbSlice = 1;
80 const int nbPixel = 128;
81 double totalAngleSpan = 180.; // in degree
82
83 double angleOfDetector =
84 135.; // angle of detector relative to the incident direction of the primary protons //
85 double distanceObjectDetector = 22.; // 22 mm
86 double radiusOfDetector = 5.; // 5 mm
87 // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex angle of the right
88 // circular cone in radian double theta = 14.726*TMath::DegToRad(); // in radian
89 double theta = 70 * TMath::DegToRad(); // in radian
90 // double theta = 70*TMath::DegToRad(); // in radian
91 // double theta = DegreeToRadian(70);
92
93 int P_interrupt = 1; // Projection of interruption
94
95 //***********************************************************************
96 //**************************Detection parameters (end)*******************
97 //***********************************************************************
98
99 // assuming there is one interruption
100 FILE* input1 = fopen("../RT7_GDP_1Projs_1Slice_128Pixels_2000000_4MeV/GammaAtExit.dat", "rb");
101 FILE* out =
102 fopen("../RT7_GDP_1Projs_1Slice_128Pixels_2000000_4MeV/PixeEvent_std_AtExit.DAT", "wb");
103 // FILE* temp;
104 // temp =fopen("temp.DAT","wb");
105
106 if (input1 == NULL) {
107 printf("error for opening the input GammaAtExit.dat file\n");
108 return;
109 }
110
111 RunInfo runInfo;
112 PixeEvent pixeEvent;
113 Point centerOfDetector;
114 Point gammaMomentum;
115 long long count1 = 0;
116 long long count2 = 0;
117 int runID = -1; // index of simulations, namely runID, starting from 0
118 std::vector<PixeEvent> eventVec;
119
120 // ************************************************************(begin)
121 // **********************READ FIRST FILE***********************
122 // ************************************************************
123 while (fread(&runInfo, sizeof(RunInfo), 1, input1)) {
124 runID++;
125 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
126 int remain = runID % (nbSlice * nbPixel);
127 runInfo.sliceIndex = remain / nbPixel;
128 runInfo.pixelIndex = remain % nbPixel;
129 if (runInfo.projectionIndex == P_interrupt) {
130 runID--;
131 break;
132 }
133
134 int nbParticle = runInfo.nbParticle;
135
136 std::vector<ParticleInfo> gammaAtExit(nbParticle);
137 fread(&gammaAtExit[0], sizeof(ParticleInfo), nbParticle, input1);
138
139 //***********************************************************************
140 //**************************Print information (begin)********************
141 //***********************************************************************
142
143 printf("-1--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
144 runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
145
146 //***********************************************************************
147 //**************************Print information (end)**********************
148 //***********************************************************************
149
150 // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
151 // source direction and detector, which should be constant when source is rotating
152 double ra =
153 DegreeToRadian(angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
154 centerOfDetector.m_x = distanceObjectDetector * cos(ra);
155 centerOfDetector.m_y = distanceObjectDetector * sin(ra);
156 centerOfDetector.m_z = 0;
157
158 for (int i = 0; i < nbParticle; ++i) {
159 // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
160 if (gammaAtExit[i].energy_keV >= 40.95 || gammaAtExit[i].energy_keV <= 0.9)
161 continue; // gamma selection
162
163 gammaMomentum.m_x = gammaAtExit[i].mx;
164 gammaMomentum.m_y = gammaAtExit[i].my;
165 gammaMomentum.m_z = gammaAtExit[i].mz;
166
167 if (!IsDetected(centerOfDetector, gammaMomentum, theta))
168 continue;
169 else {
170 pixeEvent.energy_10eV = floor(100 * gammaAtExit[i].energy_keV + 0.5);
171 pixeEvent.projectionIndex = runInfo.projectionIndex;
172 pixeEvent.sliceIndex = runInfo.sliceIndex;
173 pixeEvent.pixelIndex = runInfo.pixelIndex;
174
175 eventVec.push_back(pixeEvent);
176 count1++;
177 }
178 }
179 }
180 printf("---------------Number of PixeEvent in the first file: %lld------------------------\n",
181 count1);
182 fclose(input1);
183 // fclose(temp);
184
185 // ************************************************************(end)
186 // **********************READ FIRST FILE***********************
187 // ************************************************************
188
189 // ************************************************************(begin)
190 // **********************READ SECOND FILE**********************
191 // ************************************************************
192 // temp =fopen("temp.DAT","rb");
193
194 PixeEvent pp;
195 PixeEvent p;
196
197 for (int i = 0; i < nbProjection; ++i) {
198 int size = eventVec.size();
199 for (int j = 0; j < size; ++j) {
200 p = eventVec[j];
201 pp.energy_10eV = p.energy_10eV;
202 pp.projectionIndex = p.projectionIndex + i;
203 pp.sliceIndex = p.sliceIndex; // index of slices should be reset, starting from 0
204 pp.pixelIndex = p.pixelIndex;
205 // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_10eV=%d\n",
206 // pp.projectionIndex, pp.sliceIndex, pp.pixelIndex, pp.energy_10eV);
207 fwrite(&pp, 7, 1, out);
208 }
209 }
210
211 // ************************************************************(end)
212 // **********************READ SECOND FILE**********************
213 // ************************************************************
214
215 // fclose(temp);
216 fclose(out);
217
218 // Recheck the output file in case
219 FILE* input2 =
220 fopen("../RT7_GDP_1Projs_1Slice_128Pixels_2000000_4MeV/PixeEvent_std_AtExit.DAT", "rb");
221 PixeEvent ppp;
222 int proj = -1;
223 while (fread(&ppp, 7, 1, input2)) {
224 if (ppp.projectionIndex != proj) {
225 printf("__ProjectionIndex=%d\n", ppp.projectionIndex);
226 proj = ppp.projectionIndex;
227 }
228 // if(proj<20) printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_10eV=%d\n",
229 // ppp.projectionIndex, ppp.sliceIndex, ppp.pixelIndex, ppp.energy_10eV);
230 }
231 fclose(input2);
232 }
233