stim_pixe_tomography/scripts/Check_PixeEventFile.C

0001 //***********************************************************************************************************
0002 //     Check_PixeEventFile.C
0003 // Root command file
0004 // Use it by typing in the command line of Root terminal: root Check_PixeEventFile.C
0005 //
0006 //
0007 // More information is available in UserGuide
0008 // Created by Z.LI LP2i Bordeaux 2022
0009 //***********************************************************************************************************
0010
0011 #include <math.h>
0012 #include <stdint.h>
0013 #include <stdio.h>
0014 #include <string.h>
0015
0016 #include <vector>
0017 // using namespace std;
0018
0019 #define PI 3.14159265f
0020
0021 // Define a structure to read and write each event in the required binary format
0022 struct PixeEvent
0023 {
0024   uint16_t energy_10eV;
0025   uint16_t pixelIndex;
0026   uint16_t sliceIndex;
0027   uint8_t projectionIndex;
0028 };
0029 struct ParticleInfo
0030 {
0031   float energy_keV;
0032   float mx;
0033   float my;
0034   float mz;
0035 };
0036 struct RunInfo
0037 {
0038   // uint_16t
0039   uint8_t projectionIndex;  // 1 byte
0040   uint16_t sliceIndex;  //
0041   uint16_t pixelIndex;
0042   uint32_t nbParticle;  // 4 bytes int
0043 };
0044
0045 double DegreeToRadian(double degree)
0046 {
0047   return (PI * degree / 180.);
0048 }
0049
0050 struct Point
0051 {
0052   double m_x;
0053   double m_y;
0054   double m_z;
0055 };
0056 bool IsDetected(Point poi1, Point poi2, double theta)
0057 {
0058   double a = (poi1.m_x * poi2.m_x + poi1.m_y * poi2.m_y + poi1.m_z * poi2.m_z)
0059              / sqrt(poi1.m_x * poi1.m_x + poi1.m_y * poi1.m_y + poi1.m_z * poi1.m_z)
0060              / sqrt(poi2.m_x * poi2.m_x + poi2.m_y * poi2.m_y + poi2.m_z * poi2.m_z);
0061   if (a > 1.0) a = 1;
0062   if (a < -1.0) a = -1;
0063   double r = acos(a);
0064   if (r > theta)
0065     return false;
0066   else
0067     return true;
0068 }
0069
0070 void Check_PixeEventFile()
0071 {
0072   FILE* input2 =
0073     fopen("../build/PixeEvent_std_AtExit_Detector135_Aperture70_50Projections.DAT", "rb");
0074   PixeEvent ppp;
0075   int proj = -1;
0076   while (fread(&ppp, 7, 1, input2)) {
0077     if (ppp.projectionIndex != proj) {
0078       printf("__ProjectionIndex=%d\n", ppp.projectionIndex);
0079       proj = ppp.projectionIndex;
0080     }
0081     // printf("__ProjectionIndex=%d, SliceIndex=%d, PixelIndex=%d, Energy_10eV=%d\n",
0082     // ppp.projectionIndex, ppp.sliceIndex, ppp.pixelIndex, ppp.energy_10eV);
0083   }
0084   fclose(input2);
0085 }