| Geant4 Cross Reference |
1 //******************************************************************************************
2 // BinToStd_GammaAtCreation.C
3 // Root command file
4 // Type: root BinToStd_GammaAtCreation.C
5 //
6 // Read the output file GammaAtCreation.dat that is generated by Geant4
7 // tomography simulation It read all the gamma at creation information, and
8 // rewrite the events in a binary file PixeEvent_std_AtCreation.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 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) { return (PI * degree / 180.); }
54
55 bool IsDetected(Point poi1, Point poi2, double theta)
56 {
57 double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
58 / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
59 / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
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 parameters (begin)*****************
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 detector relative to the incident
80 // direction of the primary protons //
81 double distanceObjectDetector = 22.; // 22 mm
82 double radiusOfDetector = 5.; // 5 mm
83 // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex
84 // angle of the right circular cone in radian
85 double theta = 70 * TMath::DegToRad(); // in radian
86
87 //***********************************************************************
88 //**************************Detection parameters (end)*******************
89 //***********************************************************************
90
91 FILE* input = fopen("../build/GammaAtCreation.dat", "rb");
92 FILE* out = fopen("../build/PixeEvent_std_AtCreation.DAT", "wb");
93
94 if (input == NULL) {
95 printf("error for opening the input GammaAtCreation.dat file\n");
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, namely runID, starting from 0
105
106 // while(!feof(input)) //if not the end, read
107 while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
108 runID++;
109 // if(runID==5) continue;
110 int nbParticle = runInfo.nbParticle;
111
112 //(begin)*****************************************************************
113 // the following codes are used only when in the simulation
114 // the index of projection, slice and pixel is not
115 // correctly configured
116 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
117 int remain = runID % (nbSlice * nbPixel);
118 runInfo.sliceIndex = remain / nbPixel;
119 runInfo.pixelIndex = remain % nbPixel;
120 //(end)******************************************************************
121
122 //***********************************************************************
123 //**************************Print information (begin)********************
124 //***********************************************************************
125
126 printf(
127 "---------RunID=%d:\nProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,"
128 "nbParticle = %d\n",
129 runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
130
131 //***********************************************************************
132 //**************************Print information (end)**********************
133 //***********************************************************************
134
135 if (!nbParticle) continue;
136 std::vector<ParticleInfo> gammaAtCreation(nbParticle);
137 fread(&gammaAtCreation[0], sizeof(ParticleInfo), nbParticle, input);
138
139 // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means
140 // the angle between source direction and detector, which should be constant
141 // when source is rotating
142 double ra = TMath::DegToRad()
143 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
144 centerOfDetector.m_x = distanceObjectDetector * cos(ra);
145 centerOfDetector.m_y = distanceObjectDetector * sin(ra);
146 centerOfDetector.m_z = 0;
147
148 for (int i = 0; i < nbParticle; ++i) {
149 // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
150 if (gammaAtCreation[i].energy_keV >= 40.95 || gammaAtCreation[i].energy_keV <= 0.9)
151 continue; // gamma selection
152
153 gammaMomentum.m_x = gammaAtCreation[i].mx;
154 gammaMomentum.m_y = gammaAtCreation[i].my;
155 gammaMomentum.m_z = gammaAtCreation[i].mz;
156
157 if (!IsDetected(centerOfDetector, gammaMomentum, theta))
158 continue;
159 else {
160 pixeEvent.energy_10eV = floor(100 * gammaAtCreation[i].energy_keV + 0.5);
161 pixeEvent.projectionIndex = runInfo.projectionIndex;
162 pixeEvent.sliceIndex = runInfo.sliceIndex;
163 pixeEvent.pixelIndex = runInfo.pixelIndex;
164 fwrite(&pixeEvent, 7, 1, out);
165 count++;
166
167 //***********************************************************************
168 //**************************Print information (begin)********************
169 //***********************************************************************
170
171 // printf("momentum: (%f, %f, %f), energy: %f keV %d 10eV\n",
172 // gammaAtCreation[i].mx, gammaAtCreation[i].my, gammaAtCreation[i].mz,
173 // gammaAtCreation[i].energy_keV, pixeEvent.energy_10eV);
174
175 //***********************************************************************
176 //**************************Print information (end)**********************
177 //***********************************************************************
178 }
179 }
180 }
181 printf(
182 "---------------Number of PixeEvent in total: "
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_AtCreation.DAT","rb");
190 // PixeEvent p;
191 // while(fread(&p, 7, 1, input2))
192 // {
193 // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,
194 // Energy_10eV=%d\n", p.projectionIndex, p.sliceIndex, p.pixelIndex,
195 // p.energy_10eV);
196
197 // }
198 // fclose(input2);
199 }
200