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