| Geant4 Cross Reference |
1 //******************************************** 1
2 // BinToStd_GammaAtCreation.C
3 // Root command file
4 // Type: root BinToStd_GammaAtCreation.C
5 //
6 // Read the output file GammaAtCreation.dat th
7 // tomography simulation It read all the gamma
8 // rewrite the events in a binary file PixeEve
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 e
23 struct PixeEvent
24 {
25 uint16_t energy_10eV;
26 uint16_t pixelIndex;
27 uint16_t sliceIndex;
28 uint8_t projectionIndex;
29 };
30 struct ParticleInfo
31 {
32 float energy_keV;
33 float mx;
34 float my;
35 float mz;
36 };
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 struct Point
47 {
48 double m_x;
49 double m_y;
50 double m_z;
51 };
52
53 // double DegreeToRadian(double degree) { retu
54
55 bool IsDetected(Point poi1, Point poi2, double
56 {
57 double a = (poi1.m_x * poi2.m_x + poi1.m_y *
58 / sqrt(poi1.m_x * poi1.m_x + poi1
59 / sqrt(poi2.m_x * poi2.m_x + poi2
60 if (a > 1.0) a = 1;
61 if (a < -1.0) a = -1;
62 double r = acos(a);
63 if (r > theta)
64 return false;
65 else
66 return true;
67 }
68 void BinToStd_GammaAtCreation()
69 {
70 //******************************************
71 //**************************Detection parame
72 //******************************************
73
74 const int nbProjection = 10;
75 const int nbSlice = 1;
76 const int nbPixel = 20;
77 double totalAngleSpan = 180.; // in degree
78
79 double angleOfDetector = 135.; // angle of
80 // direction
81 double distanceObjectDetector = 22.; // 22
82 double radiusOfDetector = 5.; // 5 mm
83 // double theta = atan(radiusOfDetector/dist
84 // angle of the right circular cone in radia
85 double theta = 70 * TMath::DegToRad(); // i
86
87 //******************************************
88 //**************************Detection parame
89 //******************************************
90
91 FILE* input = fopen("../build/GammaAtCreatio
92 FILE* out = fopen("../build/PixeEvent_std_At
93
94 if (input == NULL) {
95 printf("error for opening the input GammaA
96 return;
97 }
98
99 RunInfo runInfo;
100 PixeEvent pixeEvent;
101 Point centerOfDetector;
102 Point gammaMomentum;
103 long long count = 0;
104 int runID = -1; // index of simulations, na
105
106 // while(!feof(input)) //if not the end, rea
107 while (fread(&runInfo, sizeof(RunInfo), 1, i
108 runID++;
109 // if(runID==5) continue;
110 int nbParticle = runInfo.nbParticle;
111
112 //(begin)*********************************
113 // the following codes are used only when
114 // the index of projection, slice and pixe
115 // correctly configured
116 runInfo.projectionIndex = runID / (nbSlice
117 int remain = runID % (nbSlice * nbPixel);
118 runInfo.sliceIndex = remain / nbPixel;
119 runInfo.pixelIndex = remain % nbPixel;
120 //(end)***********************************
121
122 //****************************************
123 //**************************Print informat
124 //****************************************
125
126 printf(
127 "---------RunID=%d:\nProjectionIndex=%d,
128 "nbParticle = %d\n",
129 runID, runInfo.projectionIndex, runInfo.
130
131 //****************************************
132 //**************************Print informat
133 //****************************************
134
135 if (!nbParticle) continue;
136 std::vector<ParticleInfo> gammaAtCreation(
137 fread(&gammaAtCreation[0], sizeof(Particle
138
139 // angleOfDetector+totalAngleSpan/nbProjec
140 // the angle between source direction and
141 // when source is rotating
142 double ra = TMath::DegToRad()
143 * (angleOfDetector + totalAngl
144 centerOfDetector.m_x = distanceObjectDetec
145 centerOfDetector.m_y = distanceObjectDetec
146 centerOfDetector.m_z = 0;
147
148 for (int i = 0; i < nbParticle; ++i) {
149 // gamma selection: energy should be low
150 if (gammaAtCreation[i].energy_keV >= 40.
151 continue; // gamma selection
152
153 gammaMomentum.m_x = gammaAtCreation[i].m
154 gammaMomentum.m_y = gammaAtCreation[i].m
155 gammaMomentum.m_z = gammaAtCreation[i].m
156
157 if (!IsDetected(centerOfDetector, gammaM
158 continue;
159 else {
160 pixeEvent.energy_10eV = floor(100 * ga
161 pixeEvent.projectionIndex = runInfo.pr
162 pixeEvent.sliceIndex = runInfo.sliceIn
163 pixeEvent.pixelIndex = runInfo.pixelIn
164 fwrite(&pixeEvent, 7, 1, out);
165 count++;
166
167 //************************************
168 //**************************Print info
169 //************************************
170
171 // printf("momentum: (%f, %f, %f), ene
172 // gammaAtCreation[i].mx, gammaAtCreat
173 // gammaAtCreation[i].energy_keV, pixe
174
175 //************************************
176 //**************************Print info
177 //************************************
178 }
179 }
180 }
181 printf(
182 "---------------Number of PixeEvent in tot
183 "%lld------------------------\n",
184 count);
185 fclose(input);
186 fclose(out);
187
188 // Recheck the output file in case
189 // FILE* input2 = fopen("PixeEvent_std_AtCre
190 // PixeEvent p;
191 // while(fread(&p, 7, 1, input2))
192 // {
193 // printf("__ProjectionIndex=%d, SliceIndex=
194 // Energy_10eV=%d\n", p.projectionIndex, p.s
195 // p.energy_10eV);
196
197 // }
198 // fclose(input2);
199 }
200