stim_pixe_tomography/scripts/BinToStd_GammaAtCreation.C

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