stim_pixe_tomography/scripts/BinToStd_proton_position.C

0001 //***********************************************************************************************************
0002 // BinToStd_proton_position.C
0003 // Root command file
0004 // Type: root BinToStd_proton_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 StimEvent_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 StimEvent
0043 {
0044   uint16_t energy_keV;  // different from Pixe Event, it is in keV
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_proton_position()
0107 {
0108   // printf("%f %f %f\n", acos(1), acos(-1), acos(0));
0109   // return;
0110
0111   //***********************************************************************
0112   //**************************Detection parameters (begin)*****************
0113   //***********************************************************************
0114
0115   const int nbProjection = 1;
0116   const int nbSlice = 1;
0117   const int nbPixel = 1;
0118   double totalAngleSpan = 180.;  // in degree
0119
0120   double angleOfDetector = 0.;  // angle of detector relative to the incident
0121
0122   double distanceObjectDetector = 22000.;  // um
0123
0124   // double theta = atan(radiusOfDetector/distanceObjectDetector); //half apex
0125   // angle of the right circular cone in radian
0126   double theta = 10.2 * TMath::DegToRad();  // in radian
0127   double radiusOfDetector = distanceObjectDetector * tan(theta);
0128   bool usePosition = true;
0129
0130   //***********************************************************************
0131   //**************************Detection parameters (end)*******************
0132   //***********************************************************************
0133
0134   FILE* input = fopen("../build/ProtonAtExit.dat", "rb");
0135   FILE* out = fopen("../build/StimEvent_std", "wb");
0136
0137   if (input == NULL) {
0138     printf("error for opening the input file\n");
0139     return;
0140   }
0141
0142   RunInfo runInfo;
0143   StimEvent stimEvent;
0144   Point centerOfDetector;
0145   Point protonMomentum;
0146   Point protonPosition;
0147   Point intersectionPoint;
0148   long long count = 0;
0149   int runID = -1;  // index of simulations, namely runID, starting from 0
0150
0151   // while(!feof(input)) //if not the end, read
0152   while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
0153     runID++;
0154     int nbParticle = runInfo.nbParticle;
0155
0156     //(begin)****************************************************************
0157     // the following codes are used only when in the simulation
0158     // the index of projection, slice and pixel is not
0159     // correctly configured
0160     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0161     int remain = runID % (nbSlice * nbPixel);
0162     runInfo.sliceIndex = remain / nbPixel;
0163     runInfo.pixelIndex = remain % nbPixel;
0164
0165     //(end)*******************************************************************
0166
0167     //***********************************************************************
0168     //**************************Print information (begin)********************
0169     //***********************************************************************
0170
0171     printf(
0172       "---------RunID=%d:\nProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,"
0173       "nbParticle = %d\n",
0174       runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0175
0176     //***********************************************************************
0177     //**************************Print information (end)**********************
0178     //***********************************************************************
0179
0180     if (!nbParticle) continue;
0181     std::vector<ParticleInfo> protonAtExit(nbParticle);
0182     fread(&protonAtExit[0], sizeof(ParticleInfo), nbParticle, input);
0183
0184     // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means
0185     // the angle between source direction and detector, which should be constant
0186     // when source is rotating
0187     double ra = TMath::DegToRad()
0188                 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
0189     centerOfDetector.m_x = distanceObjectDetector * cos(ra);
0190     centerOfDetector.m_y = distanceObjectDetector * sin(ra);
0191     centerOfDetector.m_z = 0;
0192
0193     for (int i = 0; i < nbParticle; ++i) {
0194       // proton selection: energy should be lower than 4095 keV
0195       if (protonAtExit[i].energy_keV >= 4095) continue;
0196
0197       protonMomentum.m_x = protonAtExit[i].mx;
0198       protonMomentum.m_y = protonAtExit[i].my;
0199       protonMomentum.m_z = protonAtExit[i].mz;
0200
0201       if (!usePosition) {
0202         if (!IsDetected(centerOfDetector, protonMomentum, theta)) continue;
0203       }
0204       else {
0205         double c =
0206           distanceObjectDetector * (protonMomentum.m_x * cos(ra) + protonMomentum.m_y * sin(ra));
0207         if (IsEqual(0, c, eps, releps)) continue;  // parallel
0208
0209         protonPosition.m_x = protonAtExit[i].x;
0210         protonPosition.m_y = protonAtExit[i].y;
0211         protonPosition.m_z = protonAtExit[i].z;
0212
0213         double t = (distanceObjectDetector * distanceObjectDetector
0214                     - protonPosition.m_x * distanceObjectDetector * cos(ra)
0215                     - protonPosition.m_y * distanceObjectDetector * sin(ra))
0216                    / c;
0217
0218         intersectionPoint.m_x = protonPosition.m_x + protonMomentum.m_x * t;
0219         intersectionPoint.m_y = protonPosition.m_y + protonMomentum.m_y * t;
0220         intersectionPoint.m_z = protonPosition.m_z + protonMomentum.m_z * t;
0221
0222         if (!IsDetected_position(centerOfDetector, intersectionPoint, radiusOfDetector)) continue;
0223
0224         // printf("      t = %f, intersection point: (%f, %f, %f) centor of detector: (%f, %f, %f)
0225         // 111=%f, 222=%f \n", t, intersectionPoint.m_x,intersectionPoint.m_y,intersectionPoint.m_z,
0226         // centerOfDetector.m_x,centerOfDetector.m_y,centerOfDetector.m_z,
0227         // (distanceObjectDetector*distanceObjectDetector-protonPosition.m_x*distanceObjectDetector*cos(ra)
0228         // -protonPosition.m_y*distanceObjectDetector*sin(ra)), c);
0229
0230         // printf("      distanceObjectDetector = %f, protonPosition.m_x=%f,
0231         // distanceObjectDetector*cos(ra)=%f, protonPosition.m_y=%f,
0232         // distanceObjectDetector*sin(ra)=%f\n", distanceObjectDetector, protonPosition.m_x,
0233         // distanceObjectDetector*cos(ra),
0234         // protonPosition.m_y,
0235         // distanceObjectDetector*sin(ra));
0236
0237         double tt = (intersectionPoint.m_x - protonPosition.m_x) * protonMomentum.m_x
0238                     + (intersectionPoint.m_y - protonPosition.m_y) * protonMomentum.m_y
0239                     + (intersectionPoint.m_z - protonPosition.m_z) * protonMomentum.m_z;
0240         if (tt < 0) continue;
0241       }
0242
0243       stimEvent.energy_10eV = floor(100 * protonAtExit[i].energy_keV + 0.5);
0244       stimEvent.projectionIndex = runInfo.projectionIndex;
0245       stimEvent.sliceIndex = runInfo.sliceIndex;
0246       stimEvent.pixelIndex = runInfo.pixelIndex;
0247       fwrite(&stimEvent, 7, 1, out);
0248       count++;
0249
0250       //***********************************************************************
0251       //**************************Print information
0252       //(begin)********************
0253       //***********************************************************************
0254
0255       if (!usePosition) {
0256         printf(
0257           "---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, momentum: "
0258           "(%f, %f, %f), energy: %f keV\n",
0259           i, runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex,
0260           protonAtExit[i].mx, protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV);
0261       }
0262       else {
0263         // printf("---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d,
0264         // momentum: (%f, %f, %f), energy: %f keV, position: (%f, %f, %f)\n", i, runID,
0265         // runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, protonAtExit[i].mx,
0266         // protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV, protonAtExit[i].x,
0267         // protonAtExit[i].y, protonAtExit[i].z);
0268
0269         printf(
0270           "---------id = %d, RunID=%d ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, momentum: "
0271           "(%f, %f, %f), energy: %f keV\n",
0272           i, runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex,
0273           protonAtExit[i].mx, protonAtExit[i].my, protonAtExit[i].mz, protonAtExit[i].energy_keV);
0274       }
0275
0276       //***********************************************************************
0277       //**************************Print information (end)**********************
0278       //***********************************************************************
0279     }
0280   }
0281   printf(
0282     "\n---------------Number of StimEvent in total: "
0283     "%lld------------------------\n",
0284     count);
0285   fclose(input);
0286   fclose(out);
0287
0288   // FILE* input2;
0289   // input2 = fopen("StimEvent_std.DAT","rb");
0290   // StimEvent p;
0291   // double eventId = -1;
0292   // while(fread(&p, 7, 1, input2))
0293   // {
0294
0295   // if(p.projectionIndex == 8 &&p.sliceIndex ==64 && p.pixelIndex==64)
0296   // {
0297   // eventId++;
0298   // printf("StimEvent_%.0f ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_keV=%d keV\n",
0299   // eventId, p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_keV);
0300
0301   // }
0302
0303   // }
0304   // fclose(input2);
0305 }