stim_pixe_tomography/scripts/Concatenate_BinToStd_GammaAtCreation.C

0001 //***********************************************************************************************************
0002 // Concatenate_BinToStd_GammaAtCreation.C
0003 // Root command file
0004 // Type: root Concatenate_BinToStd_GammaAtCreation.C
0005 //
0006 // It is used in case of interruption
0007 // Read 2 output files GammaAtCreation_1.dat and GammaAtCreation_2.dat that are generated by Geant4
0008 // tomography simulation. It reads all the gamma at creation information, and rewrite the events in
0009 // a binary file PixeEvent_std_AtCreation.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 PixeEvent
0025 {
0026   uint16_t energy_10eV;
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 void Concatenate_BinToStd_GammaAtCreation()
0067 {
0068   //***********************************************************************
0069   //**************************Detection parameters (begin)*****************
0070   //***********************************************************************
0071
0072   const int nbProjection = 10;
0073   const int nbSlice = 1;
0074   const int nbPixel = 20;
0075   double totalAngleSpan = 180.;  // in degree
0076
0077   double angleOfDetector =
0078     135.;  // angle of detector relative to the incident direction of the primary protons //
0079   double distanceObjectDetector = 22.;  // 22 mm
0080   double radiusOfDetector = 5.;  // 5 mm
0081   // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex angle of the right
0082   // circular cone in radian
0083   double theta = 70 * TMath::DegToRad();  // in radian
0084
0085   int P_interrupt = 6;  // Projection of interruption
0086
0087   //***********************************************************************
0088   //**************************Detection parameters (end)*******************
0089   //***********************************************************************
0090
0091   // assuming there is one interruption
0092   FILE* input1 = fopen("../build/GammaAtCreation_1.dat", "rb");
0093   FILE* input2 = fopen("../build/GammaAtCreation_2.dat", "rb");
0094   FILE* out = fopen("../build/PixeEvent_std_AtCreation.DAT", "wb");
0095
0096   if (input1 == NULL) {
0097     printf("error for opening the input GammaAtCreation_1.dat file\n");
0098     return;
0099   }
0100   if (input2 == NULL) {
0101     printf("error for opening the input GammaAtCreation_2.dat file\n");
0102     return;
0103   }
0104
0105   RunInfo runInfo;
0106   PixeEvent pixeEvent;
0107   Point centerOfDetector;
0108   Point gammaMomentum;
0109   long long count1 = 0;
0110   long long count2 = 0;
0111   int runID = -1;  // index of simulations, namely runID, starting from 0
0112
0113   // ************************************************************(begin)
0114   // **********************READ FIRST FILE***********************
0115   // ************************************************************
0116   while (fread(&runInfo, sizeof(RunInfo), 1, input1)) {
0117     runID++;
0118
0119     //(begin)***************************************************************
0120     // the following codes are used only when in the simulation
0121     // the index of projection, slice and pixel is not
0122     // correctly configured
0123     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0124     int remain = runID % (nbSlice * nbPixel);
0125     runInfo.sliceIndex = remain / nbPixel;
0126     runInfo.pixelIndex = remain % nbPixel;
0127     //(end)******************************************************************
0128
0129     if (runInfo.projectionIndex == P_interrupt) {
0130       runID--;
0131       break;
0132     }
0133
0134     int nbParticle = runInfo.nbParticle;
0135
0136     std::vector<ParticleInfo> gammaAtCreation(nbParticle);
0137     fread(&gammaAtCreation[0], sizeof(ParticleInfo), nbParticle, input1);
0138
0139     // if(runInfo.sliceIndex!=1) continue;
0140     // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
0141     // if(runInfo.sliceIndex!=31) continue;
0142
0143     //***********************************************************************
0144     //**************************Print information (begin)********************
0145     //***********************************************************************
0146
0147     printf("-1--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
0148            runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0149
0150     //***********************************************************************
0151     //**************************Print information (end)**********************
0152     //***********************************************************************
0153
0154     // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
0155     // source direction and detector, which should be constant when source is rotating
0156     double ra = TMath::DegToRad()
0157                 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
0158     centerOfDetector.m_x = distanceObjectDetector * cos(ra);
0159     centerOfDetector.m_y = distanceObjectDetector * sin(ra);
0160     centerOfDetector.m_z = 0;
0161
0162     for (int i = 0; i < nbParticle; ++i) {
0163       // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
0164       if (gammaAtCreation[i].energy_keV >= 40.95 || gammaAtCreation[i].energy_keV <= 0.9)
0165         continue;  // gamma selection
0166
0167       gammaMomentum.m_x = gammaAtCreation[i].mx;
0168       gammaMomentum.m_y = gammaAtCreation[i].my;
0169       gammaMomentum.m_z = gammaAtCreation[i].mz;
0170
0171       if (!IsDetected(centerOfDetector, gammaMomentum, theta))
0172         continue;
0173       else {
0174         pixeEvent.energy_10eV = floor(100 * gammaAtCreation[i].energy_keV + 0.5);
0175         pixeEvent.projectionIndex = runInfo.projectionIndex;
0176         pixeEvent.sliceIndex = runInfo.sliceIndex;
0177         pixeEvent.pixelIndex = runInfo.pixelIndex;
0178         fwrite(&pixeEvent, 7, 1, out);
0179         count1++;
0180       }
0181     }
0182   }
0183   printf("---------------Number of PixeEvent in the first file: %lld------------------------\n",
0184          count1);
0185   fclose(input1);
0186
0187   // ************************************************************
0188   // **********************READ FIRST FILE (end)*****************
0189   // ************************************************************
0190
0191   // ************************************************************
0192   // **********************READ SECOND FILE (begin)**************
0193   // ************************************************************
0194   while (fread(&runInfo, sizeof(RunInfo), 1, input2)) {
0195     runID++;
0196
0197     //(begin)***************************************************************
0198     // the following codes are used only when in the simulation
0199     // the index of projection, slice and pixel is not
0200     // correctly configured
0201     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0202     int remain = runID % (nbSlice * nbPixel);
0203     runInfo.sliceIndex = remain / nbPixel;
0204     runInfo.pixelIndex = remain % nbPixel;
0205     //(end)******************************************************************
0206
0207     int nbParticle = runInfo.nbParticle;
0208
0209     //***********************************************************************
0210     //**************************Print information (begin)********************
0211     //***********************************************************************
0212
0213     printf("-2--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
0214            runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0215
0216     //***********************************************************************
0217     //**************************Print information (end)**********************
0218     //***********************************************************************
0219
0220     if (!nbParticle) continue;
0221     std::vector<ParticleInfo> gammaAtCreation(nbParticle);
0222     fread(&gammaAtCreation[0], sizeof(ParticleInfo), nbParticle, input2);
0223
0224     // if(runInfo.sliceIndex!=1) continue;
0225     // if(runInfo.sliceIndex!=31) continue;
0226     // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
0227
0228     // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
0229     // source direction and detector, which should be constant when source is rotating
0230     double ra = TMath::DegToRad()
0231                 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
0232     centerOfDetector.m_x = distanceObjectDetector * cos(ra);
0233     centerOfDetector.m_y = distanceObjectDetector * sin(ra);
0234     centerOfDetector.m_z = 0;
0235
0236     for (int i = 0; i < nbParticle; ++i) {
0237       // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
0238       if (gammaAtCreation[i].energy_keV >= 40.95 || gammaAtCreation[i].energy_keV <= 0.9)
0239         continue;  // gamma selection
0240
0241       gammaMomentum.m_x = gammaAtCreation[i].mx;
0242       gammaMomentum.m_y = gammaAtCreation[i].my;
0243       gammaMomentum.m_z = gammaAtCreation[i].mz;
0244
0245       if (!IsDetected(centerOfDetector, gammaMomentum, theta))
0246         continue;
0247       else {
0248         pixeEvent.energy_10eV = floor(100 * gammaAtCreation[i].energy_keV + 0.5);
0249         pixeEvent.projectionIndex = runInfo.projectionIndex;
0250         pixeEvent.sliceIndex = runInfo.sliceIndex;
0251         pixeEvent.pixelIndex = runInfo.pixelIndex;
0252         fwrite(&pixeEvent, 7, 1, out);
0253         count2++;
0254       }
0255     }
0256   }
0257   printf("---------------Number of PixeEvent in in the second file: %lld------------------------\n",
0258          count2);
0259
0260   // ************************************************************
0261   // **********************READ SECOND FILE (end)****************
0262   // ************************************************************
0263
0264   printf("---------------Number of PixeEvent in total: %lld------------------------\n",
0265          count1 + count2);
0266   fclose(input2);
0267   fclose(out);
0268
0269   // Recheck the output file in case
0270   // FILE* input2 = fopen("PixeEvent_std_AtCreation.DAT","rb");
0271   // PixeEvent p;
0272   // while(fread(&p, 7, 1, input2))
0273   // {
0274   // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_10eV=%d\n",
0275   // p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_10eV);
0276
0277   // }
0278   // fclose(input2);
0279 }