stim_pixe_tomography/scripts/Spectrum_proton.C

0001 //***********************************************************************************************************
0002 // Spectrum_proton.C
0003 // Root command file
0004 // Type: root Spectrum_proton.C
0005 //
0006 // It visualizes the spectrum of protons and plots a histogram by reading
0007 // simulation result ProtonAtExit.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 void Plot(vector<double>& energies, int bin, double eMin, double eMax)
0040 {
0041   auto mycanvas = new TCanvas("canvas", "canvas", 800, 50, 600, 600);
0042   gPad->SetLeftMargin(0.15);
0043
0044   // unit is in keV
0045   auto hist = new TH1D("hist (keV)", "Spectrum of protons", bin, eMin, eMax);
0046
0047   for (int i = 0; i < energies.size(); ++i) {
0048     hist->Fill(energies[i]);
0049   }
0050
0051   hist->Draw();
0052   hist->GetXaxis()->SetTitle("Energy (keV)");
0053   hist->GetYaxis()->SetTitle("Counts");
0054   hist->GetXaxis()->CenterTitle();
0055   hist->GetYaxis()->CenterTitle();
0056
0057   mycanvas->Print("spectrum_proton.png");
0058 }
0059
0060 void Spectrum_proton()
0061 {
0062   FILE* input = fopen("../build/ProtonAtExit.dat", "rb");
0063   if (input == NULL) {
0064     printf("error for opening the input file\n");
0065     return;
0066   }
0067
0068   //***********************************************************************
0069   //**************************Selection parameters (begin)*****************
0070   //***********************************************************************
0071
0072   const int nbProjection = 10;
0073   const int nbSlice = 128;
0074   const int nbPixel = 20;
0075
0076   int projection_index_begin = 0;  // starter of the projection selected
0077   int projection_index_end = 0;  // end of the projection selected
0078
0079   int slice_index_begin = 64;  // starter of the slice selected
0080   int slice_index_end = 64;  // end of the slice selected
0081
0082   //********************Parameters for spectrum***************************
0083   int bin = 100;
0084   double eMin = 0;  // keV
0085   double eMax = 0;  // keV
0086
0087   //***********************************************************************
0088   //**************************Selection parameters (end)*******************
0089   //***********************************************************************
0090
0091   RunInfo runInfo;
0092   vector<double> energies;
0093   int runID = -1;  // index of simulations, namely runID, starting from 0
0094   // while(!feof(input)) //if not the end, read
0095   while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
0096     runID++;
0097     int nbParticle = runInfo.nbParticle;
0098
0099     // ***********the following codes are used
0100     // if**************************************(begin)
0101     // ***********the index of projection, slice and pixel is not correctly
0102     // configured in the simulation
0103     runInfo.projectionIndex = runID / (nbSlice * nbPixel);
0104     int remain = runID % (nbSlice * nbPixel);
0105     runInfo.sliceIndex = remain / nbPixel;
0106     runInfo.pixelIndex = remain % nbPixel;
0107     //******************************************************************************************(end)
0108
0109     if (!nbParticle) continue;
0110     std::vector<ParticleInfo> proton(nbParticle);
0111     fread(&proton[0], sizeof(ParticleInfo), nbParticle, input);
0112
0113     if (runInfo.projectionIndex >= projection_index_begin
0114         && runInfo.projectionIndex <= projection_index_end)
0115     {
0116       if (runInfo.sliceIndex >= slice_index_begin && runInfo.sliceIndex <= slice_index_end) {
0117         for (int i = 0; i < nbParticle; ++i) {
0118           energies.push_back(proton[i].energy_keV);
0119           if (proton[i].energy_keV > eMax) eMax = proton[i].energy_keV;
0120         }
0121       }
0122     }
0123     else
0124       break;
0125   }
0126
0127   fclose(input);
0128   Plot(energies, bin, eMin, eMax + 10);
0129 }