| Geant4 Cross Reference |
1 //***********************************************************************************************************
2 // TomoSpectrum_HIST.C
3 // Root command file
4 // Type: root TomoSpectrum_HIST.C
5 //
6 // It visualizes the spectrum of X-rays and plots a histogram by reading PixeEvent data
7 //
8 // More information is available in UserGuide
9 // Created by Z.LI LP2i Bordeaux 2022
10 //***********************************************************************************************************
11
12 #include <math.h>
13 #include <stdint.h>
14 #include <stdio.h>
15 #include <string.h>
16
17 #include <vector>
18 // using namespace std;
19
20 struct PixeEvent
21 {
22 uint16_t energy_10eV;
23 uint16_t pixelIndex;
24 uint16_t sliceIndex;
25 uint8_t projectionIndex;
26 };
27
28 void Plot(vector<double>& energies, int bin, double eMin, double eMax)
29 {
30 gROOT->Reset();
31
32 auto mycanvas = new TCanvas("canvas", "canvas", 800, 50, 600, 600);
33 mycanvas->ToggleEventStatus();
34
35 gPad->SetLeftMargin(0.15);
36
37 auto hist = new TH1D("HIST", "Spectrum", bin, eMin, eMax);
38
39 for (int i = 0; i < energies.size(); ++i) {
40 hist->Fill(energies[i]);
41 }
42
43 hist->Draw();
44 hist->SetTitle("TOMO Energy Spectrum");
45 hist->GetXaxis()->SetTitle("ADC channels");
46 hist->GetYaxis()->SetTitle("Nb events");
47
48 hist->GetXaxis()->CenterTitle();
49 hist->GetYaxis()->CenterTitle();
50
51 // hist->GetYaxis()->SetTitleOffset(2);
52
53 mycanvas->Print("TomoSpectrum_hist.png");
54 }
55
56 void TomoSpectrum_HIST()
57 {
58 FILE* input = fopen("../build/PixeEvent_std_AtExit_Detector135_Aperture70.DAT", "rb");
59
60 if (input == NULL) {
61 printf("----------error for opening the input file--------------\n");
62 return;
63 }
64
65 //***********************************************************************
66 //**************************Selection parameters (begin)*****************
67 //***********************************************************************
68 const int nbProjection = 10;
69 const int nbSlice = 1;
70 const int nbPixel = 20;
71
72 int projection_index_begin = 0; // starter of the projection selected
73 int projection_index_end = 0; // end of the projection selected
74
75 int slice_index_begin = 0; // starter of the slice selected
76 int slice_index_end = 0; // end of the slice selected
77
78 //********************Parameters for spectrum***************************
79 int nbChannels = 4096;
80 double eMin = 0; // initialization
81 double eMax = 0; //
82
83 //***********************************************************************
84 //**************************Selection parameters (end)*******************
85 //***********************************************************************
86
87 vector<double> energies;
88 PixeEvent p;
89 while (fread(&p, 7, 1, input)) {
90 if (p.projectionIndex >= projection_index_begin && p.projectionIndex <= projection_index_end) {
91 if (p.sliceIndex >= slice_index_begin && p.sliceIndex <= slice_index_end) {
92 energies.push_back(p.energy_10eV);
93 if (p.energy_10eV > eMax) eMax = p.energy_10eV;
94 }
95 }
96 }
97
98 fclose(input);
99
100 long int size = energies.size();
101 vector<int> X(nbChannels); // save channels 1-4096
102 vector<int> Y(nbChannels);
103 for (long int i = 0; i < size; ++i) {
104 int energy = energies[i];
105 Y[energy - 1] = Y[energy - 1] + 1;
106 }
107 for (int i = 0; i < nbChannels; ++i) {
108 X[i] = 1 + i;
109 }
110
111 FILE* out = fopen("Spectrum_hist.txt", "wb");
112
113 for (int i = 0; i < nbChannels; ++i) {
114 fprintf(out, "%d\t%d\n", X[i], Y[i]);
115 }
116
117 fclose(out);
118
119 Plot(energies, nbChannels, 0, nbChannels);
120 }
121