stim_pixe_tomography/scripts/Concatenate_BinToStd_GammaAtExit_fabricate.C

0001 //***********************************************************************************************************
0002 //     Concatenate_BinToStd_GammaAtExit_fabricate.C
0003 // Root command file
0004 // Use it by typing in the command line of Root terminal: root
0005 // Concatenate_BinToStd_GammaAtExit_fabricate.C
0006 //
0007 //
0008 //
0009 // More information is available in UserGuide
0010 // Created by Z.LI LP2i Bordeaux 2022
0011 //***********************************************************************************************************
0012
0013 #include <math.h>
0014 #include <stdint.h>
0015 #include <stdio.h>
0016 #include <string.h>
0017
0018 #include <vector>
0019 // using namespace std;
0020
0021 #define PI 3.14159265f
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 double DegreeToRadian(double degree)
0048 {
0049   return (PI * degree / 180.);
0050 }
0051
0052 struct Point
0053 {
0054   double m_x;
0055   double m_y;
0056   double m_z;
0057 };
0058 bool IsDetected(Point poi1, Point poi2, double theta)
0059 {
0060   double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
0061              / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
0062              / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
0063   if (a > 1.0) a = 1;
0064   if (a < -1.0) a = -1;
0065   double r = acos(a);
0066   if (r > theta)
0067     return false;
0068   else
0069     return true;
0070 }
0071
0072 void Concatenate_BinToStd_GammaAtExit_fabricate()
0073 {
0074   //***********************************************************************
0075   //**************************Detection parameters (begin)*****************
0076   //***********************************************************************
0077
0078   const int nbProjection = 100;
0079   const int nbSlice = 1;
0080   const int nbPixel = 128;
0081   double totalAngleSpan = 180.;  // in degree
0082
0083   double angleOfDetector =
0084     135.;  // angle of detector relative to the incident direction of the primary protons //
0085   double distanceObjectDetector = 22.;  // 22 mm
0086   double radiusOfDetector = 5.;  // 5 mm
0087   // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex angle of the right
0088   // circular cone in radian double theta = 14.726*TMath::DegToRad();    // in radian
0089   double theta = 70 * TMath::DegToRad();  // in radian
0090   // double theta = 70*TMath::DegToRad();    // in radian
0091   // double theta = DegreeToRadian(70);
0092
0093   int P_interrupt = 1;  // Projection of interruption
0094
0095   //***********************************************************************
0096   //**************************Detection parameters (end)*******************
0097   //***********************************************************************
0098
0099   // assuming there is one interruption
0100   FILE* input1 = fopen("../RT7_GDP_1Projs_1Slice_128Pixels_2000000_4MeV/GammaAtExit.dat", "rb");
0101   FILE* out =
0102     fopen("../RT7_GDP_1Projs_1Slice_128Pixels_2000000_4MeV/PixeEvent_std_AtExit.DAT", "wb");
0103   // FILE* temp;
0104   // temp =fopen("temp.DAT","wb");
0105
0106   if (input1 == NULL) {
0107     printf("error for opening the input GammaAtExit.dat file\n");
0108     return;
0109   }
0110
0111   RunInfo runInfo;
0112   PixeEvent pixeEvent;
0113   Point centerOfDetector;
0114   Point gammaMomentum;
0115   long long count1 = 0;
0116   long long count2 = 0;
0117   int runID = -1;  // index of simulations, namely runID, starting from 0
0118   std::vector<PixeEvent> eventVec;
0119
0120   // ************************************************************(begin)
0121   // **********************READ FIRST FILE***********************
0122   // ************************************************************
0123   while (fread(&runInfo, sizeof(RunInfo), 1, input1)) {
0124     runID++;
0125     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0126     int remain = runID % (nbSlice * nbPixel);
0127     runInfo.sliceIndex = remain / nbPixel;
0128     runInfo.pixelIndex = remain % nbPixel;
0129     if (runInfo.projectionIndex == P_interrupt) {
0130       runID--;
0131       break;
0132     }
0133
0134     int nbParticle = runInfo.nbParticle;
0135
0136     std::vector<ParticleInfo> gammaAtExit(nbParticle);
0137     fread(&gammaAtExit[0], sizeof(ParticleInfo), nbParticle, input1);
0138
0139     //***********************************************************************
0140     //**************************Print information (begin)********************
0141     //***********************************************************************
0142
0143     printf("-1--runId %d, ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
0144            runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0145
0146     //***********************************************************************
0147     //**************************Print information (end)**********************
0148     //***********************************************************************
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 =
0153       DegreeToRadian(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
0175         eventVec.push_back(pixeEvent);
0176         count1++;
0177       }
0178     }
0179   }
0180   printf("---------------Number of PixeEvent in the first file: %lld------------------------\n",
0181          count1);
0182   fclose(input1);
0183   // fclose(temp);
0184
0185   // ************************************************************(end)
0186   // **********************READ FIRST FILE***********************
0187   // ************************************************************
0188
0189   // ************************************************************(begin)
0190   // **********************READ SECOND FILE**********************
0191   // ************************************************************
0192   // temp =fopen("temp.DAT","rb");
0193
0194   PixeEvent pp;
0195   PixeEvent p;
0196
0197   for (int i = 0; i < nbProjection; ++i) {
0198     int size = eventVec.size();
0199     for (int j = 0; j < size; ++j) {
0200       p = eventVec[j];
0201       pp.energy_10eV = p.energy_10eV;
0202       pp.projectionIndex = p.projectionIndex + i;
0203       pp.sliceIndex = p.sliceIndex;  // index of slices should be reset, starting from 0
0204       pp.pixelIndex = p.pixelIndex;
0205       // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_10eV=%d\n",
0206       // pp.projectionIndex, pp.sliceIndex, pp.pixelIndex, pp.energy_10eV);
0207       fwrite(&pp, 7, 1, out);
0208     }
0209   }
0210
0211   // ************************************************************(end)
0212   // **********************READ SECOND FILE**********************
0213   // ************************************************************
0214
0215   // fclose(temp);
0216   fclose(out);
0217
0218   // Recheck the output file in case
0219   FILE* input2 =
0220     fopen("../RT7_GDP_1Projs_1Slice_128Pixels_2000000_4MeV/PixeEvent_std_AtExit.DAT", "rb");
0221   PixeEvent ppp;
0222   int proj = -1;
0223   while (fread(&ppp, 7, 1, input2)) {
0224     if (ppp.projectionIndex != proj) {
0225       printf("__ProjectionIndex=%d\n", ppp.projectionIndex);
0226       proj = ppp.projectionIndex;
0227     }
0228     // if(proj<20) printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_10eV=%d\n",
0229     // ppp.projectionIndex, ppp.sliceIndex, ppp.pixelIndex, ppp.energy_10eV);
0230   }
0231   fclose(input2);
0232 }