| Geant4 Cross Reference |
1 //***********************************************************************************************************
2 // TomoSpectrum_HIST_proton.C
3 // Root command file
4 // Type: root TomoSpectrum_HIST_proton.C
5 //
6 // It visualizes the spectrum of protons and plots a histogram by reading StimEvent 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 StimEvent
21 {
22 uint16_t energy_keV; // different from Pixe Event, it is in keV
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_proton.png");
54 }
55
56 void TomoSpectrum_HIST_proton()
57 {
58 FILE* input = fopen("../build/StimEvent_std_Detector0_Aperture10.2.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 = 128;
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 = 64; // starter of the slice selected
76 int slice_index_end = 64; // 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 StimEvent s;
89 while (fread(&s, 7, 1, input)) {
90 if (s.projectionIndex >= projection_index_begin && s.projectionIndex <= projection_index_end) {
91 if (s.sliceIndex >= slice_index_begin && s.sliceIndex <= slice_index_end) {
92 energies.push_back(s.energy_keV);
93 if (s.energy_keV > eMax) eMax = s.energy_keV;
94 }
95 }
96 }
97
98 fclose(input);
99
100 if (eMax > 4096) printf("---error in data----\n");
101 long int size = energies.size();
102 vector<int> X(nbChannels); // save channels 1-4096
103 vector<int> Y(nbChannels);
104
105 for (long int i = 0; i < size; ++i) {
106 int energy = energies[i];
107 Y[energy - 1] = Y[energy - 1] + 1;
108 }
109 for (int i = 0; i < nbChannels; ++i) {
110 X[i] = 1 + i;
111 }
112
113 FILE* out = fopen("Spectrum_hist_proton.txt", "wb");
114
115 for (int i = 0; i < nbChannels; ++i) {
116 fprintf(out, "%d\t%d\n", X[i], Y[i]);
117 }
118
119 fclose(out);
120
121 Plot(energies, nbChannels, 0, nbChannels);
122 }
123