stim_pixe_tomography/scripts/TomoSpectrum.C

0001 //***********************************************************************************************************
0002 // TomoSpectrum.C
0003 // Root command file
0004 // Type: root TomoSpectrum.C
0005 //
0006 // It visualizes the spectrum of X-rays and plots a graph 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(int nbChannels, vector<int>& X, vector<int>& Y)
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 graph = new TGraph(nbChannels, X.data(), Y.data());
0038   graph->SetLineColor(8);
0039
0040   graph->Draw("AL");
0041
0042   graph->SetLineColor(8);
0043   graph->SetTitle("TOMO Energy Spectrum");
0044   graph->GetXaxis()->SetTitle("ADC channels");
0045   graph->GetYaxis()->SetTitle("Nb events");
0046   graph->GetXaxis()->CenterTitle();
0047   graph->GetYaxis()->CenterTitle();
0048
0049   mycanvas->Print("TomoSpectrum.png");
0050 }
0051
0052 void TomoSpectrum()
0053 {
0054   FILE* input = fopen("../build/PixeEvent_std_AtExit_Detector135_Aperture70.DAT", "rb");
0055
0056   if (input == NULL) {
0057     printf("----------error for opening the input file--------------\n");
0058     return;
0059   }
0060
0061   //***********************************************************************
0062   //**************************Selection parameters (begin)*****************
0063   //***********************************************************************
0064   const int nbProjection = 10;
0065   const int nbSlice = 1;
0066   const int nbPixel = 20;
0067
0068   int projection_index_begin = 0;  // starter of the projection selected
0069   int projection_index_end = 0;  // end of the projection selected
0070
0071   int slice_index_begin = 0;  // starter of the slice selected
0072   int slice_index_end = 0;  // end of the slice selected
0073
0074   //***********************************************************************
0075   //**************************Selection parameters (end)*******************
0076   //***********************************************************************
0077
0078   int nbChannels = 4096;
0079   vector<int> X(nbChannels);  // save channels 1-4096, index X: 0-4095
0080   vector<int> Y(nbChannels);  // save event counts for channel 1-4096, index Y: 0-4095
0081   PixeEvent p;
0082
0083   while (fread(&p, 7, 1, input)) {
0084     if (p.projectionIndex >= projection_index_begin && p.projectionIndex <= projection_index_end) {
0085       if (p.sliceIndex >= slice_index_begin && p.sliceIndex <= slice_index_end) {
0086         // printf("%d %d %d\n",p.projectionIndex, p.sliceIndex, p.energy_10eV);
0087         Y[p.energy_10eV - 1] = Y[p.energy_10eV - 1] + 1;
0088       }
0089     }
0090   }
0091   fclose(input);
0092   for (int i = 0; i < nbChannels; ++i) {
0093     X[i] = 1 + i;
0094   }
0095
0096   FILE* out = fopen("Spectrum.txt", "wb");
0097   for (int i = 0; i < nbChannels; ++i) {
0098     fprintf(out, "%d\t%d\n", X[i], Y[i]);
0099   }
0100
0101   fclose(out);
0102   Plot(nbChannels, X, Y);
0103 }