stim_pixe_tomography/scripts/BinToStd_gamma_position.C

0001 //***********************************************************************************************************
0002 // BinToStd_gamma_position.C
0003 // Root command file
0004 // Type: root BinToStd_gamma_position.C
0005 //
0006 // Read the X-ray output file that is generated by Geant4 tomography
0007 // simulation. It reads gamma information, either at creation, or at exit, and rewrite the events
0008 // in a binary file PixeEvent_std.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 bool IsEqual(double a, double b, double eps, double releps)
0023 {
0024   if (a == b) {
0025     return true;
0026   }
0027
0028   if (fabs(a - b) <= releps * fabs(b)) {
0029     return true;
0030   }
0031
0032   if (fabs(a - b) < eps) {
0033     return true;
0034   }
0035
0036   return false;
0037 }
0038 double eps = 1e-20;  // absolut difference
0039 double releps = 1e-10;  // relative difference
0040
0041 // Define a structure to read and write each event in the required binary format
0042 struct PixeEvent
0043 {
0044   uint16_t energy_10eV;
0045   uint16_t pixelIndex;
0046   uint16_t sliceIndex;
0047   uint8_t projectionIndex;
0048 };
0049 struct ParticleInfo
0050 {
0051   float energy_keV;
0052   float mx;
0053   float my;
0054   float mz;
0055   float x;
0056   float y;
0057   float z;
0058 };
0059 struct RunInfo
0060 {
0061   // uint_16t
0062   uint8_t projectionIndex;  // 1 byte
0063   uint16_t sliceIndex;  //
0064   uint16_t pixelIndex;
0065   uint32_t nbParticle;  // 4 bytes int
0066 };
0067 struct Point
0068 {
0069   double m_x;
0070   double m_y;
0071   double m_z;
0072 };
0073
0074 bool IsDetected(Point poi1, Point poi2, double theta)
0075 {
0076   double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
0077              / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
0078              / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
0079   if (a > 1.0) a = 1;
0080   if (a < -1.0) a = -1;
0081   double r = acos(a);
0082   if (r > theta)
0083     return false;
0084   else {
0085     // printf("      acos: %f, radius: %f\n", r, theta);
0086     return true;
0087   }
0088 }
0089
0090 bool IsDetected_position(Point poi1, Point poi2, double r)
0091 {
0092   double a = sqrt((poi1.m_x - poi2.m_x) * (poi1.m_x - poi2.m_x)
0093                   + (poi1.m_y - poi2.m_y) * (poi1.m_y - poi2.m_y)
0094                   + (poi1.m_z - poi2.m_z) * (poi1.m_z - poi2.m_z));
0095
0096   // if(a <= r)    return true;
0097   if (a > r)
0098     return false;
0099
0100   else {
0101     // printf("      distance of two points: %f, radius: %f\n", a, r);
0102     return true;
0103   }
0104 }
0105
0106 void BinToStd_gamma_position()
0107 {
0108   //***********************************************************************
0109   //**************************Detection parameters (begin)*****************
0110   //***********************************************************************
0111
0112   const int nbProjection = 1;
0113   const int nbSlice = 1;
0114   const int nbPixel = 1;
0115   double totalAngleSpan = 180.;  // in degree
0116
0117   double angleOfDetector = 135.;  // angle of detector relative to the incident
0118
0119   double distanceObjectDetector = 22000.;  // um
0120
0121   // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex
0122   // angle of the right circular cone in radian
0123   double theta = 14.726 * TMath::DegToRad();  // in radian
0124   double radiusOfDetector = distanceObjectDetector * tan(theta);
0125   bool usePosition = true;
0126
0127   //***********************************************************************
0128   //**************************Detection parameters (end)*******************
0129   //***********************************************************************
0130
0131   FILE* input = fopen("../build/GammaAtExit.dat", "rb");
0132   FILE* out = fopen("../build/PixeEvent_std_AtExit.DAT", "wb");
0133
0134   if (input == NULL) {
0135     printf("error for opening the input file\n");
0136     return;
0137   }
0138
0139   RunInfo runInfo;
0140   PixeEvent pixeEvent;
0141   Point centerOfDetector;
0142   Point gammaMomentum;
0143   Point gammaPosition;
0144   Point intersectionPoint;
0145   long long count = 0;
0146   int runID = -1;  // index of simulations, namely runID, starting from 0
0147
0148   // while(!feof(input)) //if not the end, read
0149   while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
0150     runID++;
0151     int nbParticle = runInfo.nbParticle;
0152
0153     //(begin)****************************************************************
0154     // the following codes are used only when in the simulation
0155     // the index of projection, slice and pixel is not
0156     // correctly configured
0157     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0158     int remain = runID % (nbSlice * nbPixel);
0159     runInfo.sliceIndex = remain / nbPixel;
0160     runInfo.pixelIndex = remain % nbPixel;
0161
0162     //(end)******************************************************************
0163
0164     //***********************************************************************
0165     //**************************Print information (begin)********************
0166     //***********************************************************************
0167
0168     printf(
0169       "---------RunID=%d:\nProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,"
0170       "nbParticle = %d\n",
0171       runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0172
0173     //***********************************************************************
0174     //**************************Print information (end)**********************
0175     //***********************************************************************
0176
0177     if (!nbParticle) continue;
0178     std::vector<ParticleInfo> gammaAtExit(nbParticle);
0179     fread(&gammaAtExit[0], sizeof(ParticleInfo), nbParticle, input);
0180
0181     // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means
0182     // the angle between source direction and detector, which should be constant
0183     // when source is rotating
0184     double ra = TMath::DegToRad()
0185                 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
0186     centerOfDetector.m_x = distanceObjectDetector * cos(ra);
0187     centerOfDetector.m_y = distanceObjectDetector * sin(ra);
0188     centerOfDetector.m_z = 0;
0189
0190     for (int i = 0; i < nbParticle; ++i) {
0191       // gamma selection: energy should be lower than 4095*10eV = 49.45 keV
0192       if (gammaAtExit[i].energy_keV >= 40.95 || gammaAtExit[i].energy_keV <= 0.9) continue;
0193
0194       gammaMomentum.m_x = gammaAtExit[i].mx;
0195       gammaMomentum.m_y = gammaAtExit[i].my;
0196       gammaMomentum.m_z = gammaAtExit[i].mz;
0197
0198       if (!usePosition) {
0199         if (!IsDetected(centerOfDetector, gammaMomentum, theta)) continue;
0200       }
0201       else {
0202         double c =
0203           distanceObjectDetector * (gammaMomentum.m_x * cos(ra) + gammaMomentum.m_y * sin(ra));
0204         if (IsEqual(0, c, eps, releps)) continue;  // parallel
0205
0206         gammaPosition.m_x = gammaAtExit[i].x;
0207         gammaPosition.m_y = gammaAtExit[i].y;
0208         gammaPosition.m_z = gammaAtExit[i].z;
0209
0210         double t = (distanceObjectDetector * distanceObjectDetector
0211                     - gammaPosition.m_x * distanceObjectDetector * cos(ra)
0212                     - gammaPosition.m_y * distanceObjectDetector * sin(ra))
0213                    / c;
0214
0215         intersectionPoint.m_x = gammaPosition.m_x + gammaMomentum.m_x * t;
0216         intersectionPoint.m_y = gammaPosition.m_y + gammaMomentum.m_y * t;
0217         intersectionPoint.m_z = gammaPosition.m_z + gammaMomentum.m_z * t;
0218
0219         if (!IsDetected_position(centerOfDetector, intersectionPoint, radiusOfDetector)) continue;
0220
0221         // printf("      t = %f, intersection point: (%f, %f, %f) centor of detector: (%f, %f, %f)
0222         // 111=%f, 222=%f \n", t, intersectionPoint.m_x,intersectionPoint.m_y,intersectionPoint.m_z,
0223         // centerOfDetector.m_x,centerOfDetector.m_y,centerOfDetector.m_z,
0224         // (distanceObjectDetector*distanceObjectDetector-gammaPosition.m_x*distanceObjectDetector*cos(ra)
0225         // -gammaPosition.m_y*distanceObjectDetector*sin(ra)), c);
0226
0227         // printf("      distanceObjectDetector = %f, gammaPosition.m_x=%f,
0228         // distanceObjectDetector*cos(ra)=%f, gammaPosition.m_y=%f,
0229         // distanceObjectDetector*sin(ra)=%f\n", distanceObjectDetector, gammaPosition.m_x,
0230         // distanceObjectDetector*cos(ra),
0231         // gammaPosition.m_y,
0232         // distanceObjectDetector*sin(ra));
0233
0234         double tt = (intersectionPoint.m_x - gammaPosition.m_x) * gammaMomentum.m_x
0235                     + (intersectionPoint.m_y - gammaPosition.m_y) * gammaMomentum.m_y
0236                     + (intersectionPoint.m_z - gammaPosition.m_z) * gammaMomentum.m_z;
0237         if (tt < 0) continue;
0238       }
0239
0240       pixeEvent.energy_10eV = floor(100 * gammaAtExit[i].energy_keV + 0.5);
0241       pixeEvent.projectionIndex = runInfo.projectionIndex;
0242       pixeEvent.sliceIndex = runInfo.sliceIndex;
0243       pixeEvent.pixelIndex = runInfo.pixelIndex;
0244       fwrite(&pixeEvent, 7, 1, out);
0245       count++;
0246
0247       //***********************************************************************
0248       //**************************Print information (begin)********************
0249       //***********************************************************************
0250
0251       if (!usePosition) {
0252         printf(
0253           "---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, momentum: "
0254           "(%f, %f, %f), energy: %f keV\n",
0255           i, runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex,
0256           gammaAtExit[i].mx, gammaAtExit[i].my, gammaAtExit[i].mz, gammaAtExit[i].energy_keV);
0257       }
0258       else {
0259         // printf("---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,
0260         // momentum: (%f, %f, %f), energy: %f keV, position: (%f, %f, %f)\n", i, runID,
0261         // runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, gammaAtExit[i].mx,
0262         // gammaAtExit[i].my, gammaAtExit[i].mz, gammaAtExit[i].energy_keV, gammaAtExit[i].x,
0263         // gammaAtExit[i].y, gammaAtExit[i].z);
0264
0265         printf(
0266           "---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, momentum: "
0267           "(%f, %f, %f), energy: %f keV\n",
0268           i, runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex,
0269           gammaAtExit[i].mx, gammaAtExit[i].my, gammaAtExit[i].mz, gammaAtExit[i].energy_keV);
0270       }
0271
0272       //***********************************************************************
0273       //**************************Print information (end)**********************
0274       //***********************************************************************
0275     }
0276   }
0277   printf(
0278     "\n---------------Number of PixeEvent in total: "
0279     "%lld------------------------\n",
0280     count);
0281   fclose(input);
0282   fclose(out);
0283
0284   // Recheck the output file in case
0285   // FILE* input2;
0286   // input2 = fopen("PixeEvent_std_AtExit.DAT","rb");
0287   // PixeEvent p;
0288   // while(fread(&p, 7, 1, input2))
0289   // {
0290   // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,
0291   // Energy_10eV=%d\n", p.projectionIndex, p.sliceIndex, p.pixelIndex,
0292   // p.energy_10eV);
0293
0294   // }
0295   // fclose(input2);
0296 }