Geant4 Cross Reference |
1 //******************************************** 1 2 // BinToStd_ProtonAtExit.C 3 // Root command file 4 // Type: root BinToStd_ProtonAtExit.C 5 // 6 // Read the output file ProtonAtExit.dat that 7 // It reads proton at exit information, and re 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 e 22 struct StimEvent 23 { 24 uint16_t energy_keV; // different from Pixe 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 53 { 54 double a = (poi1.m_x * poi2.m_x + poi1.m_y * 55 / sqrt(poi1.m_x * poi1.m_x + poi1 56 / sqrt(poi2.m_x * poi2.m_x + poi2 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 parame 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 inciden 77 // for proton, it is fixed to 0 degree, name 78 double angleOfDetector = 0.; 79 double distanceObjectDetector = 22.; // 22 80 double radiusOfDetector = 5.; // 5 mm 81 // double theta = atan(radiusOfDetector/dist 82 // circular cone in radian 83 double theta = 10.2 * TMath::DegToRad(); // 84 85 //****************************************** 86 //**************************Detection parame 87 //****************************************** 88 89 FILE* input = fopen("../build/ProtonAtExit.d 90 FILE* out = fopen("../build/StimEvent_std.DA 91 92 if (input == NULL) { 93 printf("error for opening the input Proton 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, rea 105 while (fread(&runInfo, sizeof(RunInfo), 1, i 106 runID++; 107 int nbParticle = runInfo.nbParticle; 108 109 //(begin)********************************* 110 // the following codes are used only when 111 // the index of projection, slice and pixe 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 //(end)*********************************** 118 119 //**************************************** 120 //**************************Print informat 121 //**************************************** 122 123 printf("---------RunID=%d: ProjectionIndex 124 runID, runInfo.projectionIndex, run 125 126 //**************************************** 127 //**************************Print informat 128 //**************************************** 129 130 if (!nbParticle) continue; 131 std::vector<ParticleInfo> protonAtExit(nbP 132 fread(&protonAtExit[0], sizeof(ParticleInf 133 134 // angleOfDetector+totalAngleSpan/nbProjec 135 // source direction and detector, which sh 136 double ra = TMath::DegToRad() 137 * (angleOfDetector + totalAngl 138 centerOfDetector.m_x = distanceObjectDetec 139 centerOfDetector.m_y = distanceObjectDetec 140 centerOfDetector.m_z = 0; 141 142 for (int i = 0; i < nbParticle; ++i) { 143 // proton selection: energy should be lo 144 if (protonAtExit[i].energy_keV >= 4095) 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, proton 151 continue; 152 else { 153 stimEvent.energy_keV = floor(protonAtE 154 stimEvent.projectionIndex = runInfo.pr 155 stimEvent.sliceIndex = runInfo.sliceIn 156 stimEvent.pixelIndex = runInfo.pixelIn 157 fwrite(&stimEvent, 7, 1, out); 158 count++; 159 // printf("energy=%f keV\n",protonAtEx 160 } 161 } 162 } 163 printf("---------------Number of StimEvent i 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 175 // { 176 // eventId++; 177 // printf("StimEvent_%.0f ProjectionIndex=%d 178 // eventId, p.projectionIndex, p.sliceIndex, 179 180 // } 181 182 // } 183 // fclose(input2); 184 } 185