stim_pixe_tomography/scripts/Spectrum_gamma.C

0001 //***********************************************************************************************************
0002 // Spectrum_gamma.C
0003 // Root command file
0004 // Type: root Spectrum_gamma.C
0005 //
0006 // It visualizes the spectrum of X-rays and plots a histogram by reading
0007 // simulation result GammaAtCreation.dat or GammaAtExit.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 RunInfo
0023 {
0024   // uint_16t
0025   uint8_t projectionIndex;  // 1 byte
0026   uint16_t sliceIndex;  //
0027   uint16_t pixelIndex;
0028   uint32_t nbParticle;  // 4 bytes int
0029 };
0030
0031 struct ParticleInfo
0032 {
0033   float energy_keV;
0034   float mx;
0035   float my;
0036   float mz;
0037 };
0038
0039 // struct ParticleInfo
0040 //{
0041 //  float energy_keV;
0042 //  float mx;
0043 //  float my;
0044 //  float mz;
0045 //  float x;
0046 //  float y;
0047 //  float z;
0048 //};
0049
0050 void Plot(vector<double>& energies, int bin, double eMin, double eMax)
0051 {
0052   auto mycanvas = new TCanvas("canvas", "canvas", 800, 50, 600, 600);
0053   gPad->SetLeftMargin(0.15);
0054
0055   // unit is in keV
0056   auto hist = new TH1D("hist (keV)", "Spectrum of photons", bin, eMin, eMax);
0057
0058   for (int i = 0; i < energies.size(); ++i) {
0059     hist->Fill(energies[i]);
0060   }
0061
0062   hist->Draw();
0063   hist->GetXaxis()->SetTitle("Energy (keV)");
0064   hist->GetYaxis()->SetTitle("Counts");
0065   hist->GetXaxis()->CenterTitle();
0066   hist->GetYaxis()->CenterTitle();
0067
0068   mycanvas->Print("spectrum_gamma.png");
0069 }
0070
0071 void Spectrum_gamma()
0072 {
0073   FILE* input = fopen("../build/GammaAtExit.dat", "rb");
0074
0075   if (input == NULL) {
0076     printf("error for opening the input file\n");
0077     return;
0078   }
0079
0080   //***********************************************************************
0081   //**************************Selection parameters (begin)*****************
0082   //***********************************************************************
0083
0084   const int nbProjection = 10;
0085   const int nbSlice = 1;
0086   const int nbPixel = 20;
0087
0088   int projection_index_begin = 0;  // starter of the projection selected
0089   int projection_index_end = 0;  // end of the projection selected
0090
0091   int slice_index_begin = 0;  // starter of the slice selected
0092   int slice_index_end = 0;  // end of the slice selected
0093
0094   //********************Parameters for spectrum***************************
0095   int bin = 100;
0096   double eMin = 0;  // keV
0097   double eMax = 0;  // keV
0098
0099   //***********************************************************************
0100   //**************************Selection parameters (end)*******************
0101   //***********************************************************************
0102
0103   RunInfo runInfo;
0104   vector<double> energies;
0105   int runID = -1;  // index of simulations, namely runID, starting from 0
0106   // while(!feof(input)) //if not the end, read
0107   while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
0108     runID++;
0109     int nbParticle = runInfo.nbParticle;
0110
0111     // ***********the following codes are used
0112     // if**************************************(begin)
0113     // ***********the index of projection, slice and pixel is not correctly
0114     // configured in the simulation
0115     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0116     int remain = runID % (nbSlice * nbPixel);
0117     runInfo.sliceIndex = remain / nbPixel;
0118     runInfo.pixelIndex = remain % nbPixel;
0119     //************************************************************************(end)
0120
0121     if (!nbParticle) continue;
0122     std::vector<ParticleInfo> particles(nbParticle);
0123     fread(&particles[0], sizeof(ParticleInfo), nbParticle, input);
0124
0125     if (runInfo.projectionIndex >= projection_index_begin
0126         && runInfo.projectionIndex <= projection_index_end)
0127     {
0128       if (runInfo.sliceIndex >= slice_index_begin && runInfo.sliceIndex <= slice_index_end) {
0129         for (int i = 0; i < nbParticle; ++i) {
0130           // printf("--%d, %.9e\n", i, particles[i].energy_keV);
0131
0132           energies.push_back(particles[i].energy_keV);
0133           if (particles[i].energy_keV > eMax) eMax = particles[i].energy_keV;
0134         }
0135       }
0136     }
0137     else
0138       break;
0139   }
0140
0141   fclose(input);
0142   Plot(energies, bin, eMin, eMax + 10);
0143 }