| Geant4 Cross Reference |
1 //******************************************** 1
2 // BinToStd_GammaAtExit.C
3 // Root command file
4 // Type: root BinToStd_GammaAtExit.C
5 //
6 // Read the output file ProtonAtExit.dat that
7 // simulation It read all the gamma at exit in
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 e
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
55 {
56 double a = (poi1.m_x * poi2.m_x + poi1.m_y *
57 / sqrt(poi1.m_x * poi1.m_x + poi1
58 / sqrt(poi2.m_x * poi2.m_x + poi2
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 parame
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
79 // direction
80 double distanceObjectDetector = 22.; // 22
81 double radiusOfDetector = 5.; // 5 mm
82 // double theta = atan(radiusOfDetector/dist
83 // angle of the right circular cone in radia
84 double theta = 70 * TMath::DegToRad(); // i
85
86 //******************************************
87 //**************************Detection parame
88 //******************************************
89
90 FILE* input = fopen("../build/GammaAtExit.da
91 FILE* out = fopen("../build/PixeEvent_std_At
92
93 if (input == NULL) {
94 printf("error for opening the input GammaA
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, na
104
105 // while(!feof(input)) //if not the end, rea
106 while (fread(&runInfo, sizeof(RunInfo), 1, i
107 runID++;
108 int nbParticle = runInfo.nbParticle;
109
110 // the following codes are used only when
111 // ************(begin) the index of projec
112 // correctly configured
113 runInfo.projectionIndex = runID / (nbSlice
114 int remain = runID % (nbSlice * nbPixel);
115 runInfo.sliceIndex = remain / nbPixel;
116 runInfo.pixelIndex = remain % nbPixel;
117 //****************************************
118
119 //****************************************
120 //**************************Print informat
121 //****************************************
122
123 printf(
124 "---------RunID=%d:\nProjectionIndex=%d,
125 "nbParticle = %d\n",
126 runID, runInfo.projectionIndex, runInfo.
127
128 //****************************************
129 //**************************Print informat
130 //****************************************
131
132 if (!nbParticle) continue;
133 std::vector<ParticleInfo> gammaAtExit(nbPa
134 fread(&gammaAtExit[0], sizeof(ParticleInfo
135
136 // angleOfDetector+totalAngleSpan/nbProjec
137 // the angle between source direction and
138 // when source is rotating
139 double ra = TMath::DegToRad()
140 * (angleOfDetector + totalAngl
141 centerOfDetector.m_x = distanceObjectDetec
142 centerOfDetector.m_y = distanceObjectDetec
143 centerOfDetector.m_z = 0;
144
145 for (int i = 0; i < nbParticle; ++i) {
146 // gamma selection: energy should be low
147 if (gammaAtExit[i].energy_keV >= 40.95 |
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, gammaM
154 continue;
155 else {
156 pixeEvent.energy_10eV = floor(100 * ga
157 pixeEvent.projectionIndex = runInfo.pr
158 pixeEvent.sliceIndex = runInfo.sliceIn
159 pixeEvent.pixelIndex = runInfo.pixelIn
160 fwrite(&pixeEvent, 7, 1, out);
161 count++;
162
163 //************************************
164 //**************************Print info
165 //************************************
166
167 // printf("momentum: (%f, %f, %f), ene
168 // gammaAtExit[i].mx, gammaAtExit[i].m
169 // gammaAtExit[i].energy_keV, pixeEven
170
171 //************************************
172 //**************************Print info
173 //************************************
174 }
175 }
176 }
177 printf(
178 "\n---------------Number of PixeEvent in t
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"
187 // PixeEvent p;
188 // while(fread(&p, 7, 1, input2))
189 // {
190 // printf("__ProjectionIndex=%d, SliceIndex=
191 // Energy_10eV=%d\n", p.projectionIndex, p.s
192 // p.energy_10eV);
193
194 // }
195 // fclose(input2);
196 }
197