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