stim_pixe_tomography/scripts/BinToStd_ProtonAtExit.C

0001 //***********************************************************************************************************
0002 // BinToStd_ProtonAtExit.C
0003 // Root command file
0004 // Type: root BinToStd_ProtonAtExit.C
0005 //
0006 // Read the output file ProtonAtExit.dat that is generated by Geant4 tomography simulation
0007 // It reads proton at exit information, and rewrite the events in a binary file StimEvent_std.DAT
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 a structure to read and write each event in the required binary format
0022 struct StimEvent
0023 {
0024   uint16_t energy_keV;  // different from Pixe Event, it is in keV
0025   uint16_t pixelIndex;
0026   uint16_t sliceIndex;
0027   uint8_t projectionIndex;
0028 };
0029
0030 struct ParticleInfo
0031 {
0032   float energy_keV;
0033   float mx;
0034   float my;
0035   float mz;
0036 };
0037 struct RunInfo
0038 {
0039   // uint_16t
0040   uint8_t projectionIndex;  // 1 byte
0041   uint16_t sliceIndex;  //
0042   uint16_t pixelIndex;
0043   uint32_t nbParticle;  // 4 bytes int
0044 };
0045 struct Point
0046 {
0047   double m_x;
0048   double m_y;
0049   double m_z;
0050 };
0051
0052 bool IsDetected(Point poi1, Point poi2, double theta)
0053 {
0054   double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
0055              / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
0056              / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
0057   if (a > 1.0) a = 1;
0058   if (a < -1.0) a = -1;
0059   double r = acos(a);
0060   if (r > theta)
0061     return false;
0062   else
0063     return true;
0064 }
0065 void BinToStd_ProtonAtExit()
0066 {
0067   //***********************************************************************
0068   //**************************Detection parameters (begin)****************
0069   //***********************************************************************
0070
0071   const int nbProjection = 10;
0072   const int nbSlice = 1;
0073   const int nbPixel = 20;
0074   double totalAngleSpan = 180.;  // in degree
0075
0076   // angle of detector relative to the incident direction of the primary protons at first projection
0077   // for proton, it is fixed to 0 degree, namely opposite to the source
0078   double angleOfDetector = 0.;
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 = 10.2 * TMath::DegToRad();  // in radian
0084
0085   //***********************************************************************
0086   //**************************Detection parameters (end)*******************
0087   //***********************************************************************
0088
0089   FILE* input = fopen("../build/ProtonAtExit.dat", "rb");
0090   FILE* out = fopen("../build/StimEvent_std.DAT", "wb");
0091
0092   if (input == NULL) {
0093     printf("error for opening the input ProtonAtExit.dat file\n");
0094     return;
0095   }
0096
0097   RunInfo runInfo;
0098   StimEvent stimEvent;
0099   Point centerOfDetector;
0100   Point protonMomentum;
0101   long long count = 0;
0102   int runID = -1;
0103
0104   // while(!feof(input)) //if not the end, read
0105   while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
0106     runID++;
0107     int nbParticle = runInfo.nbParticle;
0108
0109     //(begin)***************************************************************
0110     // the following codes are used only when in the simulation
0111     // the index of projection, slice and pixel is not
0112     // correctly configured
0113     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0114     int remain = runID % (nbSlice * nbPixel);
0115     runInfo.sliceIndex = remain / nbPixel;
0116     runInfo.pixelIndex = remain % nbPixel;
0117     //(end)******************************************************************
0118
0119     //***********************************************************************
0120     //**************************Print information (begin)********************
0121     //***********************************************************************
0122
0123     printf("---------RunID=%d: ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, nbParticle = %d\n",
0124            runID, runInfo.projectionIndex, runInfo.sliceIndex, runInfo.pixelIndex, nbParticle);
0125
0126     //***********************************************************************
0127     //**************************Print information (end)**********************
0128     //***********************************************************************
0129
0130     if (!nbParticle) continue;
0131     std::vector<ParticleInfo> protonAtExit(nbParticle);
0132     fread(&protonAtExit[0], sizeof(ParticleInfo), nbParticle, input);
0133
0134     // angleOfDetector+totalAngleSpan/nbProjection*runInfo.projectionIndex means the angle between
0135     // source direction and detector, which should be constant when source is rotating
0136     double ra = TMath::DegToRad()
0137                 * (angleOfDetector + totalAngleSpan / nbProjection * runInfo.projectionIndex);
0138     centerOfDetector.m_x = distanceObjectDetector * cos(ra);
0139     centerOfDetector.m_y = distanceObjectDetector * sin(ra);
0140     centerOfDetector.m_z = 0;
0141
0142     for (int i = 0; i < nbParticle; ++i) {
0143       // proton selection: energy should be lower than 4095 keV
0144       if (protonAtExit[i].energy_keV >= 4095) continue;  // proton selection
0145
0146       protonMomentum.m_x = protonAtExit[i].mx;
0147       protonMomentum.m_y = protonAtExit[i].my;
0148       protonMomentum.m_z = protonAtExit[i].mz;
0149
0150       if (!IsDetected(centerOfDetector, protonMomentum, theta))
0151         continue;
0152       else {
0153         stimEvent.energy_keV = floor(protonAtExit[i].energy_keV + 0.5);
0154         stimEvent.projectionIndex = runInfo.projectionIndex;
0155         stimEvent.sliceIndex = runInfo.sliceIndex;
0156         stimEvent.pixelIndex = runInfo.pixelIndex;
0157         fwrite(&stimEvent, 7, 1, out);
0158         count++;
0159         // printf("energy=%f keV\n",protonAtExit[i].energy_keV);
0160       }
0161     }
0162   }
0163   printf("---------------Number of StimEvent in total: %lld------------------------\n", count);
0164   fclose(input);
0165   fclose(out);
0166
0167   // FILE* input2;
0168   // input2 = fopen("StimEvent_std.DAT","rb");
0169   // StimEvent p;
0170   // double eventId = -1;
0171   // while(fread(&p, 7, 1, input2))
0172   // {
0173
0174   // if(p.projectionIndex == 8 &&p.sliceIndex ==64 && p.pixelIndex==64)
0175   // {
0176   // eventId++;
0177   // printf("StimEvent_%.0f ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_keV=%d keV\n",
0178   // eventId, p.projectionIndex, p.sliceIndex, p.pixelIndex, p.energy_keV);
0179
0180   // }
0181
0182   // }
0183   // fclose(input2);
0184 }