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