stim_pixe_tomography/scripts/TomoSpectrum_HIST_proton.C

0001 //***********************************************************************************************************
0002 // TomoSpectrum_HIST_proton.C
0003 // Root command file
0004 // Type: root TomoSpectrum_HIST_proton.C
0005 //
0006 // It visualizes the spectrum of protons and plots a histogram by reading StimEvent data
0007 //
0008 // More information is available in UserGuide
0009 // Created by Z.LI LP2i Bordeaux 2022
0010 //***********************************************************************************************************
0011
0012 #include <math.h>
0013 #include <stdint.h>
0014 #include <stdio.h>
0015 #include <string.h>
0016
0017 #include <vector>
0018 // using namespace std;
0019
0020 struct StimEvent
0021 {
0022   uint16_t energy_keV;  // different from Pixe Event, it is in keV
0023   uint16_t pixelIndex;
0024   uint16_t sliceIndex;
0025   uint8_t projectionIndex;
0026 };
0027
0028 void Plot(vector<double>& energies, int bin, double eMin, double eMax)
0029 {
0030   gROOT->Reset();
0031
0032   auto mycanvas = new TCanvas("canvas", "canvas", 800, 50, 600, 600);
0033   mycanvas->ToggleEventStatus();
0034
0035   gPad->SetLeftMargin(0.15);
0036
0037   auto hist = new TH1D("HIST", "Spectrum", bin, eMin, eMax);
0038
0039   for (int i = 0; i < energies.size(); ++i) {
0040     hist->Fill(energies[i]);
0041   }
0042
0043   hist->Draw();
0044   hist->SetTitle("TOMO Energy Spectrum");
0045   hist->GetXaxis()->SetTitle("ADC channels");
0046   hist->GetYaxis()->SetTitle("Nb events");
0047
0048   hist->GetXaxis()->CenterTitle();
0049   hist->GetYaxis()->CenterTitle();
0050
0051   // hist->GetYaxis()->SetTitleOffset(2);
0052
0053   mycanvas->Print("TomoSpectrum_hist_proton.png");
0054 }
0055
0056 void TomoSpectrum_HIST_proton()
0057 {
0058   FILE* input = fopen("../build/StimEvent_std_Detector0_Aperture10.2.DAT", "rb");
0059
0060   if (input == NULL) {
0061     printf("----------error for opening the input file--------------\n");
0062     return;
0063   }
0064
0065   //***********************************************************************
0066   //**************************Selection parameters (begin)*****************
0067   //***********************************************************************
0068   const int nbProjection = 10;
0069   const int nbSlice = 128;
0070   const int nbPixel = 20;
0071
0072   int projection_index_begin = 0;  // starter of the projection selected
0073   int projection_index_end = 0;  // end of the projection selected
0074
0075   int slice_index_begin = 64;  // starter of the slice selected
0076   int slice_index_end = 64;  // end of the slice selected
0077
0078   //********************Parameters for spectrum***************************
0079   int nbChannels = 4096;
0080   double eMin = 0;  // initialization
0081   double eMax = 0;  //
0082
0083   //***********************************************************************
0084   //**************************Selection parameters (end)*******************
0085   //***********************************************************************
0086
0087   vector<double> energies;
0088   StimEvent s;
0089   while (fread(&s, 7, 1, input)) {
0090     if (s.projectionIndex >= projection_index_begin && s.projectionIndex <= projection_index_end) {
0091       if (s.sliceIndex >= slice_index_begin && s.sliceIndex <= slice_index_end) {
0092         energies.push_back(s.energy_keV);
0093         if (s.energy_keV > eMax) eMax = s.energy_keV;
0094       }
0095     }
0096   }
0097
0098   fclose(input);
0099
0100   if (eMax > 4096) printf("---error in data----\n");
0101   long int size = energies.size();
0102   vector<int> X(nbChannels);  // save channels 1-4096
0103   vector<int> Y(nbChannels);
0104
0105   for (long int i = 0; i < size; ++i) {
0106     int energy = energies[i];
0107     Y[energy - 1] = Y[energy - 1] + 1;
0108   }
0109   for (int i = 0; i < nbChannels; ++i) {
0110     X[i] = 1 + i;
0111   }
0112
0113   FILE* out = fopen("Spectrum_hist_proton.txt", "wb");
0114
0115   for (int i = 0; i < nbChannels; ++i) {
0116     fprintf(out, "%d\t%d\n", X[i], Y[i]);
0117   }
0118
0119   fclose(out);
0120
0121   Plot(energies, nbChannels, 0, nbChannels);
0122 }