stim_pixe_tomography/scripts/BinToStd_GammaAtExit.C

0001 //***********************************************************************************************************
0002 // BinToStd_GammaAtExit.C
0003 // Root command file
0004 // Type: root BinToStd_GammaAtExit.C
0005 //
0006 // Read the output file ProtonAtExit.dat that is generated by Geant4 tomography
0007 // simulation It read all the gamma at exit information, and rewrite the events
0008 // in a binary file PixeEvent_std_AtExit.DAT
0009 //
0010 // More information is available in UserGuide
0011 // Created by Z.LI LP2i Bordeaux 2022
0012 //***********************************************************************************************************
0013
0014 #include <math.h>
0015 #include <stdint.h>
0016 #include <stdio.h>
0017 #include <string.h>
0018
0019 #include <vector>
0020 // using namespace std;
0021
0022 // Define a structure to read and write each event in the required binary format
0023 struct PixeEvent
0024 {
0025   uint16_t energy_10eV;
0026   uint16_t pixelIndex;
0027   uint16_t sliceIndex;
0028   uint8_t projectionIndex;
0029 };
0030
0031 struct ParticleInfo
0032 {
0033   float energy_keV;
0034   float mx;
0035   float my;
0036   float mz;
0037 };
0038
0039 struct RunInfo
0040 {
0041   // uint_16t
0042   uint8_t projectionIndex;  // 1 byte
0043   uint16_t sliceIndex;  //
0044   uint16_t pixelIndex;
0045   uint32_t nbParticle;  // 4 bytes int
0046 };
0047 struct Point
0048 {
0049   double m_x;
0050   double m_y;
0051   double m_z;
0052 };
0053
0054 bool IsDetected(Point poi1, Point poi2, double theta)
0055 {
0056   double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
0057              / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
0058              / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
0059   if (a > 1.0) a = 1;
0060   if (a < -1.0) a = -1;
0061   double r = acos(a);
0062   if (r > theta)
0063     return false;
0064   else
0065     return true;
0066 }
0067 void BinToStd_GammaAtExit()
0068 {
0069   //***********************************************************************
0070   //**************************Detection parameters (begin)*****************
0071   //***********************************************************************
0072
0073   const int nbProjection = 10;
0074   const int nbSlice = 1;
0075   const int nbPixel = 20;
0076   double totalAngleSpan = 180.;  // in degree
0077
0078   double angleOfDetector = 135.;  // angle of detector relative to the incident
0079                                   // direction of the primary protons //
0080   double distanceObjectDetector = 22.;  // 22 mm
0081   double radiusOfDetector = 5.;  // 5 mm
0082   // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex
0083   // angle of the right circular cone in radian
0084   double theta = 70 * TMath::DegToRad();  // in radian
0085
0086   //***********************************************************************
0087   //**************************Detection parameters (end)*******************
0088   //***********************************************************************
0089
0090   FILE* input = fopen("../build/GammaAtExit.dat", "rb");
0091   FILE* out = fopen("../build/PixeEvent_std_AtExit.DAT", "wb");
0092
0093   if (input == NULL) {
0094     printf("error for opening the input GammaAtExit.dat file\n");
0095     return;
0096   }
0097
0098   RunInfo runInfo;
0099   PixeEvent pixeEvent;
0100   Point centerOfDetector;
0101   Point gammaMomentum;
0102   long long count = 0;
0103   int runID = -1;  // index of simulations, namely runID, starting from 0
0104
0105   // while(!feof(input)) //if not the end, read
0106   while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
0107     runID++;
0108     int nbParticle = runInfo.nbParticle;
0109
0110     // the following codes are used only when in the simulation
0111     // ************(begin) the index of projection, slice and pixel is not
0112     // correctly configured
0113     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0114     int remain = runID % (nbSlice * nbPixel);
0115     runInfo.sliceIndex = remain / nbPixel;
0116     runInfo.pixelIndex = remain % nbPixel;
0117     //************************************************************************(end)
0118
0119     //***********************************************************************
0120     //**************************Print information (begin)********************
0121     //***********************************************************************
0122
0123     printf(
0124       "---------RunID=%d:\nProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,"
0125       "nbParticle = %d\n",
0126       runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0127
0128     //***********************************************************************
0129     //**************************Print information (end)**********************
0130     //***********************************************************************
0131
0132     if (!nbParticle) continue;
0133     std::vector<ParticleInfo> gammaAtExit(nbParticle);
0134     fread(&gammaAtExit[0], sizeof(ParticleInfo), nbParticle, input);
0135
0136     // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means
0137     // the angle between source direction and detector, which should be constant
0138     // when source is rotating
0139     double ra = TMath::DegToRad()
0140                 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
0141     centerOfDetector.m_x = distanceObjectDetector * cos(ra);
0142     centerOfDetector.m_y = distanceObjectDetector * sin(ra);
0143     centerOfDetector.m_z = 0;
0144
0145     for (int i = 0; i < nbParticle; ++i) {
0146       // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
0147       if (gammaAtExit[i].energy_keV >= 40.95 || gammaAtExit[i].energy_keV <= 0.9) continue;
0148
0149       gammaMomentum.m_x = gammaAtExit[i].mx;
0150       gammaMomentum.m_y = gammaAtExit[i].my;
0151       gammaMomentum.m_z = gammaAtExit[i].mz;
0152
0153       if (!IsDetected(centerOfDetector, gammaMomentum, theta))
0154         continue;
0155       else {
0156         pixeEvent.energy_10eV = floor(100 * gammaAtExit[i].energy_keV + 0.5);
0157         pixeEvent.projectionIndex = runInfo.projectionIndex;
0158         pixeEvent.sliceIndex = runInfo.sliceIndex;
0159         pixeEvent.pixelIndex = runInfo.pixelIndex;
0160         fwrite(&pixeEvent, 7, 1, out);
0161         count++;
0162
0163         //***********************************************************************
0164         //**************************Print information (begin)********************
0165         //***********************************************************************
0166
0167         // printf("momentum: (%f, %f, %f), energy: %f keV %d 10eV\n",
0168         // gammaAtExit[i].mx, gammaAtExit[i].my, gammaAtExit[i].mz,
0169         // gammaAtExit[i].energy_keV, pixeEvent.energy_10eV);
0170
0171         //***********************************************************************
0172         //**************************Print information (end)**********************
0173         //***********************************************************************
0174       }
0175     }
0176   }
0177   printf(
0178     "\n---------------Number of PixeEvent in total: "
0179     "%lld------------------------\n",
0180     count);
0181   fclose(input);
0182   fclose(out);
0183
0184   // Recheck the output file in case
0185   // FILE* input2;
0186   // input2 = fopen("PixeEvent_std_AtExit.DAT","rb");
0187   // PixeEvent p;
0188   // while(fread(&p, 7, 1, input2))
0189   // {
0190   // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,
0191   // Energy_10eV=%d\n", p.projectionIndex, p.sliceIndex, p.pixelIndex,
0192   // p.energy_10eV);
0193
0194   // }
0195   // fclose(input2);
0196 }