stim_pixe_tomography/scripts/Concatenate_BinToStd_ProtonAtExit.C

0001 //***********************************************************************************************************
0002 // Concatenate_BinToStd_ProtonAtExit.C
0003 // Root command file
0004 // Type: root Concatenate_BinToStd_ProtonAtExit.C
0005 //
0006 // It is used in case of interruption
0007 // Read 2 output files ProtonAtExit_1.dat and ProtonAtExit_2.dat that are generated by Geant4
0008 // tomography simulation It reads protons at exit information, and rewrite the events in a binary
0009 // file StimEvent_std.DAT
0010 //
0011 // More information is available in UserGuide
0012 // Created by Z.LI LP2i Bordeaux 2022
0013 //***********************************************************************************************************
0014
0015 #include <math.h>
0016 #include <stdint.h>
0017 #include <stdio.h>
0018 #include <string.h>
0019
0020 #include <vector>
0021 // using namespace std;
0022
0023 // Define a structure to read and write each event in the required binary format
0024 struct StimEvent
0025 {
0026   uint16_t energy_keV;  // different from Pixe Event, it is in keV
0027   uint16_t pixelIndex;
0028   uint16_t sliceIndex;
0029   uint8_t projectionIndex;
0030 };
0031 struct ParticleInfo
0032 {
0033   float energy_keV;
0034   float mx;
0035   float my;
0036   float mz;
0037 };
0038 struct RunInfo
0039 {
0040   // uint_16t
0041   uint8_t projectionIndex;  // 1 byte
0042   uint16_t sliceIndex;  //
0043   uint16_t pixelIndex;
0044   uint32_t nbParticle;  // 4 bytes int
0045 };
0046
0047 struct Point
0048 {
0049   double m_x;
0050   double m_y;
0051   double m_z;
0052 };
0053 bool IsDetected(Point poi1, Point poi2, double theta)
0054 {
0055   double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
0056              / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
0057              / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
0058   if (a > 1.0) a = 1;
0059   if (a < -1.0) a = -1;
0060   double r = acos(a);
0061   if (r > theta)
0062     return false;
0063   else
0064     return true;
0065 }
0066
0067 void Recheck()
0068 {
0069   // Recheck the output file in case
0070   FILE* input3 = fopen("../build/StimEvent_std_Detector0_Aperture10.2.DAT", "rb");
0071   StimEvent p;
0072   double eventId = -1;
0073   while (fread(&p, 7, 1, input3)) {
0074     if (p.projectionIndex == 8 && p.sliceIndex == 64 && p.pixelIndex == 10) {
0075       eventId++;
0076       printf("StimEvent_%.0f ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_keV=%d keV\n",
0077              eventId, p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_keV);
0078     }
0079   }
0080   fclose(input3);
0081 }
0082 void Concatenate_BinToStd_ProtonAtExit()
0083 {
0084   // Recheck();
0085   // return;
0086
0087   //***********************************************************************
0088   //**************************Detection parameters (begin)*****************
0089   //***********************************************************************
0090
0091   const int nbProjection = 10;
0092   const int nbSlice = 128;
0093   const int nbPixel = 20;
0094   double totalAngleSpan = 180.;  // in degree
0095
0096   // angle of detector relative to the incident direction of the primary protons at first projection
0097   // for proton, it is fixed to 0 degree, namely opposite to the source
0098   double angleOfDetector = 0.;
0099   double distanceObjectDetector = 22.;  // 22 mm
0100   double radiusOfDetector = 5.;  // 5 mm
0101   // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex angle of the right
0102   // circular cone in radian
0103   double theta = 10.2 * TMath::DegToRad();  // in radian
0104
0105   int P_interrupt = 2;  // Projection of interruption
0106
0107   //***********************************************************************
0108   //**************************Detection parameters (end)*******************
0109   //***********************************************************************
0110
0111   // assuming there is one interruption
0112   FILE* input1 = fopen("../build/ProtonAtExit_1.dat", "rb");
0113   FILE* input2 = fopen("../build/ProtonAtExit_2.dat", "rb");
0114   FILE* out = fopen("../build/StimEvent_std.DAT", "wb");
0115
0116   if (input1 == NULL) {
0117     printf("error for opening the input ProtonAtExit_1.dat file\n");
0118     return;
0119   }
0120   if (input2 == NULL) {
0121     printf("error for opening the input ProtonAtExit_2.dat file\n");
0122     return;
0123   }
0124
0125   RunInfo runInfo;
0126   StimEvent stimEvent;
0127   Point centerOfDetector;
0128   Point protonMomentum;
0129
0130   long long count1 = 0;
0131   long long count2 = 0;
0132   int runID = -1;  // index of simulations, namely runID, starting from 0
0133
0134   // ************************************************************(begin)
0135   // **********************READ FIRST FILE***********************
0136   // ************************************************************
0137   while (fread(&runInfo, sizeof(RunInfo), 1, input1)) {
0138     runID++;
0139     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0140     int remain = runID % (nbSlice * nbPixel);
0141     runInfo.sliceIndex = remain / nbPixel;
0142     runInfo.pixelIndex = remain % nbPixel;
0143     if (runInfo.projectionIndex == P_interrupt) {
0144       runID--;
0145       break;
0146     }
0147
0148     int nbParticle = runInfo.nbParticle;
0149
0150     //***********************************************************************
0151     //**************************Print information (begin)********************
0152     //***********************************************************************
0153
0154     printf("-1--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
0155            runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0156
0157     //***********************************************************************
0158     //**************************Print information (end)**********************
0159     //***********************************************************************
0160
0161     if (!nbParticle) continue;
0162     std::vector<ParticleInfo> protonAtExit(nbParticle);
0163     fread(&protonAtExit[0], sizeof(ParticleInfo), nbParticle, input1);
0164
0165     // if(runInfo.sliceIndex!=1) continue;
0166     // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
0167     // if(runInfo.sliceIndex!=31) continue;
0168
0169     // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
0170     // source direction and detector, which should be constant when source is rotating
0171     double ra = TMath::DegToRad()
0172                 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
0173     centerOfDetector.m_x = distanceObjectDetector * cos(ra);
0174     centerOfDetector.m_y = distanceObjectDetector * sin(ra);
0175     centerOfDetector.m_z = 0;
0176
0177     for (int i = 0; i < nbParticle; ++i) {
0178       // proton selection: energy should be lower than 4095 keV
0179       if (protonAtExit[i].energy_keV >= 4095) continue;  // proton selection
0180
0181       protonMomentum.m_x = protonAtExit[i].mx;
0182       protonMomentum.m_y = protonAtExit[i].my;
0183       protonMomentum.m_z = protonAtExit[i].mz;
0184
0185       if (!IsDetected(centerOfDetector, protonMomentum, theta))
0186         continue;
0187       else {
0188         stimEvent.energy_keV = floor(protonAtExit[i].energy_keV + 0.5);
0189         stimEvent.projectionIndex = runInfo.projectionIndex;
0190         stimEvent.sliceIndex = runInfo.sliceIndex;
0191         stimEvent.pixelIndex = runInfo.pixelIndex;
0192         fwrite(&stimEvent, 7, 1, out);
0193         count1++;
0194       }
0195     }
0196   }
0197   printf("---------------Number of StimEvent in the first file: %lld------------------------\n",
0198          count1);
0199   fclose(input1);
0200
0201   // ************************************************************
0202   // **********************READ FIRST FILE (end)*****************
0203   // ************************************************************
0204
0205   // ************************************************************
0206   // **********************READ SECOND FILE (begin)**************
0207   // ************************************************************
0208   while (fread(&runInfo, sizeof(RunInfo), 1, input2)) {
0209     runID++;
0210     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0211     int remain = runID % (nbSlice * nbPixel);
0212     runInfo.sliceIndex = remain / nbPixel;
0213     runInfo.pixelIndex = remain % nbPixel;
0214
0215     int nbParticle = runInfo.nbParticle;
0216
0217     //***********************************************************************
0218     //**************************Print information (begin)********************
0219     //***********************************************************************
0220
0221     printf("-2--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
0222            runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0223
0224     //***********************************************************************
0225     //**************************Print information (end)**********************
0226     //***********************************************************************
0227
0228     if (!nbParticle) continue;
0229
0230     std::vector<ParticleInfo> protonAtExit(nbParticle);
0231     fread(&protonAtExit[0], sizeof(ParticleInfo), nbParticle, input2);
0232
0233     // if(runInfo.sliceIndex!=1) continue;
0234     // if(runInfo.sliceIndex!=31) continue;
0235     // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
0236
0237     // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
0238     // source direction and detector, which should be constant when source is rotating
0239     double ra = TMath::DegToRad()
0240                 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
0241     centerOfDetector.m_x = distanceObjectDetector * cos(ra);
0242     centerOfDetector.m_y = distanceObjectDetector * sin(ra);
0243     centerOfDetector.m_z = 0;
0244
0245     for (int i = 0; i < nbParticle; ++i) {
0246       // proton selection: energy should be lower than 4095 keV
0247       if (protonAtExit[i].energy_keV >= 4095) continue;  // proton selection
0248
0249       protonMomentum.m_x = protonAtExit[i].mx;
0250       protonMomentum.m_y = protonAtExit[i].my;
0251       protonMomentum.m_z = protonAtExit[i].mz;
0252
0253       if (!IsDetected(centerOfDetector, protonMomentum, theta))
0254         continue;
0255       else {
0256         stimEvent.energy_keV = floor(protonAtExit[i].energy_keV + 0.5);
0257         stimEvent.projectionIndex = runInfo.projectionIndex;
0258         stimEvent.sliceIndex = runInfo.sliceIndex;
0259         stimEvent.pixelIndex = runInfo.pixelIndex;
0260         fwrite(&stimEvent, 7, 1, out);
0261         count2++;
0262       }
0263     }
0264   }
0265   printf("---------------Number of StimEvent in in the second file: %lld------------------------\n",
0266          count2);
0267
0268   // ************************************************************
0269   // **********************READ SECOND FILE (end)****************
0270   // ************************************************************
0271
0272   printf("---------------Number of StimEvent in total: %lld------------------------\n",
0273          count1 + count2);
0274   fclose(input2);
0275   fclose(out);
0276 }