stim_pixe_tomography/scripts/Concatenate_BinToStd_GammaAtExit.C

0001 //***********************************************************************************************************
0002 // Concatenate_BinToStd_GammaAtExit.C
0003 // Root command file
0004 // Type: root Concatenate_BinToStd_GammaAtExit.C
0005 //
0006 // It is used in case of one interruption
0007 // Read 2 output files GammaAtExit_1.dat and GammaAtExit_2.dat that are generated by Geant4
0008 // tomography simulation It reads gamma at exit information, and rewrite the events in a binary file
0009 // PixeEvent_std_AtExit.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_GammaAtExit()
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/GammaAtExit_1.dat", "rb");
0093   FILE* input2 = fopen("../build/GammaAtExit_2.dat", "rb");
0094   FILE* out = fopen("../build/PixeEvent_std_AtExit.DAT", "wb");
0095
0096   if (input1 == NULL) {
0097     printf("error for opening the input GammaAtExit_1.dat file\n");
0098     return;
0099   }
0100   if (input2 == NULL) {
0101     printf("error for opening the input GammaAtExit_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   // ************************************************************
0114   // **********************READ FIRST FILE (begin)***************
0115   // ************************************************************
0116   while (fread(&runInfo, sizeof(RunInfo), 1, input1)) {
0117     runID++;
0118
0119     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0120     int remain = runID % (nbSlice * nbPixel);
0121     runInfo.sliceIndex = remain / nbPixel;
0122     runInfo.pixelIndex = remain % nbPixel;
0123
0124     if (runInfo.projectionIndex == P_interrupt) {
0125       runID--;
0126       break;
0127     }
0128     int nbParticle = runInfo.nbParticle;
0129
0130     //***********************************************************************
0131     //**************************Print information (begin)********************
0132     //***********************************************************************
0133
0134     printf("-1--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
0135            runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0136
0137     //***********************************************************************
0138     //**************************Print information (end)**********************
0139     //***********************************************************************
0140
0141     if (!nbParticle) continue;
0142
0143     std::vector<ParticleInfo> gammaAtExit(nbParticle);
0144     fread(&gammaAtExit[0], sizeof(ParticleInfo), nbParticle, input1);
0145
0146     // if(runInfo.sliceIndex!=1) continue;
0147     // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
0148     // if(runInfo.sliceIndex!=31) continue;
0149
0150     // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
0151     // source direction and detector, which should be constant when source is rotating
0152     double ra = TMath::DegToRad()
0153                 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
0154     centerOfDetector.m_x = distanceObjectDetector * cos(ra);
0155     centerOfDetector.m_y = distanceObjectDetector * sin(ra);
0156     centerOfDetector.m_z = 0;
0157
0158     for (int i = 0; i < nbParticle; ++i) {
0159       // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
0160       if (gammaAtExit[i].energy_keV >= 40.95 || gammaAtExit[i].energy_keV <= 0.9)
0161         continue;  // gamma selection
0162
0163       gammaMomentum.m_x = gammaAtExit[i].mx;
0164       gammaMomentum.m_y = gammaAtExit[i].my;
0165       gammaMomentum.m_z = gammaAtExit[i].mz;
0166
0167       if (!IsDetected(centerOfDetector, gammaMomentum, theta))
0168         continue;
0169       else {
0170         pixeEvent.energy_10eV = floor(100 * gammaAtExit[i].energy_keV + 0.5);
0171         pixeEvent.projectionIndex = runInfo.projectionIndex;
0172         pixeEvent.sliceIndex = runInfo.sliceIndex;
0173         pixeEvent.pixelIndex = runInfo.pixelIndex;
0174         fwrite(&pixeEvent, 7, 1, out);
0175         count1++;
0176       }
0177     }
0178   }
0179   printf("---------------Number of PixeEvent in the first file: %lld------------------------\n",
0180          count1);
0181   fclose(input1);
0182
0183   // ************************************************************
0184   // **********************READ FIRST FILE (end)*****************
0185   // ************************************************************
0186
0187   // ************************************************************
0188   // **********************READ SECOND FILE (begin)**************
0189   // ************************************************************
0190   while (fread(&runInfo, sizeof(RunInfo), 1, input2)) {
0191     runID++;
0192
0193     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0194     int remain = runID % (nbSlice * nbPixel);
0195     runInfo.sliceIndex = remain / nbPixel;
0196     runInfo.pixelIndex = remain % nbPixel;
0197
0198     int nbParticle = runInfo.nbParticle;
0199
0200     //***********************************************************************
0201     //**************************Print information (begin)********************
0202     //***********************************************************************
0203
0204     printf("-2--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
0205            runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0206
0207     //***********************************************************************
0208     //**************************Print information (end)**********************
0209     //***********************************************************************
0210
0211     if (!nbParticle) continue;
0212     std::vector<ParticleInfo> gammaAtExit(nbParticle);
0213     fread(&gammaAtExit[0], sizeof(ParticleInfo), nbParticle, input2);
0214
0215     // if(runInfo.sliceIndex!=1) continue;
0216     // if(runInfo.sliceIndex!=31&&runInfo.sliceIndex!=32) continue;
0217     // if(runInfo.sliceIndex!=31) continue;
0218
0219     // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
0220     // source direction and detector, which should be constant when source is rotating
0221     double ra = TMath::DegToRad()
0222                 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
0223     centerOfDetector.m_x = distanceObjectDetector * cos(ra);
0224     centerOfDetector.m_y = distanceObjectDetector * sin(ra);
0225     centerOfDetector.m_z = 0;
0226
0227     for (int i = 0; i < nbParticle; ++i) {
0228       // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
0229       if (gammaAtExit[i].energy_keV >= 40.95 || gammaAtExit[i].energy_keV <= 0.9)
0230         continue;  // gamma selection
0231
0232       gammaMomentum.m_x = gammaAtExit[i].mx;
0233       gammaMomentum.m_y = gammaAtExit[i].my;
0234       gammaMomentum.m_z = gammaAtExit[i].mz;
0235
0236       if (!IsDetected(centerOfDetector, gammaMomentum, theta))
0237         continue;
0238       else {
0239         pixeEvent.energy_10eV = floor(100 * gammaAtExit[i].energy_keV + 0.5);
0240         pixeEvent.projectionIndex = runInfo.projectionIndex;
0241         pixeEvent.sliceIndex = runInfo.sliceIndex;
0242         pixeEvent.pixelIndex = runInfo.pixelIndex;
0243         fwrite(&pixeEvent, 7, 1, out);
0244         count2++;
0245       }
0246     }
0247   }
0248   printf("---------------Number of PixeEvent in in the second file: %lld------------------------\n",
0249          count2);
0250
0251   // ************************************************************
0252   // **********************READ SECOND FILE (end)****************
0253   // ************************************************************
0254
0255   printf("---------------Number of PixeEvent in total: %lld------------------------\n",
0256          count1 + count2);
0257   fclose(input2);
0258   fclose(out);
0259
0260   // Recheck the output file in case
0261   // FILE* input2 = fopen("PixeEvent_std_AtExit.DAT","rb");
0262   // PixeEvent p;
0263   // while(fread(&p, 7, 1, input2))
0264   // {
0265   // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_10eV=%d\n",
0266   // p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_10eV);
0267
0268   // }
0269   // fclose(input2);
0270 }