stim_pixe_tomography/scripts/TomoSpectrum_HIST.C

0001 //***********************************************************************************************************
0002 // TomoSpectrum_HIST.C
0003 // Root command file
0004 // Type: root TomoSpectrum_HIST.C
0005 //
0006 // It visualizes the spectrum of X-rays and plots a histogram by reading PixeEvent 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 PixeEvent
0021 {
0022   uint16_t energy_10eV;
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.png");
0054 }
0055
0056 void TomoSpectrum_HIST()
0057 {
0058   FILE* input = fopen("../build/PixeEvent_std_AtExit_Detector135_Aperture70.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 = 1;
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 = 0;  // starter of the slice selected
0076   int slice_index_end = 0;  // 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   PixeEvent p;
0089   while (fread(&p, 7, 1, input)) {
0090     if (p.projectionIndex >= projection_index_begin && p.projectionIndex <= projection_index_end) {
0091       if (p.sliceIndex >= slice_index_begin && p.sliceIndex <= slice_index_end) {
0092         energies.push_back(p.energy_10eV);
0093         if (p.energy_10eV > eMax) eMax = p.energy_10eV;
0094       }
0095     }
0096   }
0097
0098   fclose(input);
0099
0100   long int size = energies.size();
0101   vector<int> X(nbChannels);  // save channels 1-4096
0102   vector<int> Y(nbChannels);
0103   for (long int i = 0; i < size; ++i) {
0104     int energy = energies[i];
0105     Y[energy - 1] = Y[energy - 1] + 1;
0106   }
0107   for (int i = 0; i < nbChannels; ++i) {
0108     X[i] = 1 + i;
0109   }
0110
0111   FILE* out = fopen("Spectrum_hist.txt", "wb");
0112
0113   for (int i = 0; i < nbChannels; ++i) {
0114     fprintf(out, "%d\t%d\n", X[i], Y[i]);
0115   }
0116
0117   fclose(out);
0118
0119   Plot(energies, nbChannels, 0, nbChannels);
0120 }