| Geant4 Cross Reference |
1 //***********************************************************************************************************
2 // Spectrum_proton.C
3 // Root command file
4 // Type: root Spectrum_proton.C
5 //
6 // It visualizes the spectrum of protons and plots a histogram by reading
7 // simulation result ProtonAtExit.dat
8 //
9 // More information is available in UserGuide
10 // Created by Z.LI LP2i Bordeaux 2022
11 //***********************************************************************************************************
12
13 #include <math.h>
14 #include <stdint.h>
15 #include <stdio.h>
16 #include <string.h>
17
18 #include <vector>
19 // using namespace std;
20
21 // Define a structure to read and write each event in the required binary format
22 struct RunInfo
23 {
24 // uint_16t
25 uint8_t projectionIndex; // 1 byte
26 uint16_t sliceIndex; //
27 uint16_t pixelIndex;
28 uint32_t nbParticle; // 4 bytes int
29 };
30
31 struct ParticleInfo
32 {
33 float energy_keV;
34 float mx;
35 float my;
36 float mz;
37 };
38
39 void Plot(vector<double>& energies, int bin, double eMin, double eMax)
40 {
41 auto mycanvas = new TCanvas("canvas", "canvas", 800, 50, 600, 600);
42 gPad->SetLeftMargin(0.15);
43
44 // unit is in keV
45 auto hist = new TH1D("hist (keV)", "Spectrum of protons", bin, eMin, eMax);
46
47 for (int i = 0; i < energies.size(); ++i) {
48 hist->Fill(energies[i]);
49 }
50
51 hist->Draw();
52 hist->GetXaxis()->SetTitle("Energy (keV)");
53 hist->GetYaxis()->SetTitle("Counts");
54 hist->GetXaxis()->CenterTitle();
55 hist->GetYaxis()->CenterTitle();
56
57 mycanvas->Print("spectrum_proton.png");
58 }
59
60 void Spectrum_proton()
61 {
62 FILE* input = fopen("../build/ProtonAtExit.dat", "rb");
63 if (input == NULL) {
64 printf("error for opening the input file\n");
65 return;
66 }
67
68 //***********************************************************************
69 //**************************Selection parameters (begin)*****************
70 //***********************************************************************
71
72 const int nbProjection = 10;
73 const int nbSlice = 128;
74 const int nbPixel = 20;
75
76 int projection_index_begin = 0; // starter of the projection selected
77 int projection_index_end = 0; // end of the projection selected
78
79 int slice_index_begin = 64; // starter of the slice selected
80 int slice_index_end = 64; // end of the slice selected
81
82 //********************Parameters for spectrum***************************
83 int bin = 100;
84 double eMin = 0; // keV
85 double eMax = 0; // keV
86
87 //***********************************************************************
88 //**************************Selection parameters (end)*******************
89 //***********************************************************************
90
91 RunInfo runInfo;
92 vector<double> energies;
93 int runID = -1; // index of simulations, namely runID, starting from 0
94 // while(!feof(input)) //if not the end, read
95 while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
96 runID++;
97 int nbParticle = runInfo.nbParticle;
98
99 // ***********the following codes are used
100 // if**************************************(begin)
101 // ***********the index of projection, slice and pixel is not correctly
102 // configured in the simulation
103 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
104 int remain = runID % (nbSlice * nbPixel);
105 runInfo.sliceIndex = remain / nbPixel;
106 runInfo.pixelIndex = remain % nbPixel;
107 //******************************************************************************************(end)
108
109 if (!nbParticle) continue;
110 std::vector<ParticleInfo> proton(nbParticle);
111 fread(&proton[0], sizeof(ParticleInfo), nbParticle, input);
112
113 if (runInfo.projectionIndex >= projection_index_begin
114 && runInfo.projectionIndex <= projection_index_end)
115 {
116 if (runInfo.sliceIndex >= slice_index_begin && runInfo.sliceIndex <= slice_index_end) {
117 for (int i = 0; i < nbParticle; ++i) {
118 energies.push_back(proton[i].energy_keV);
119 if (proton[i].energy_keV > eMax) eMax = proton[i].energy_keV;
120 }
121 }
122 }
123 else
124 break;
125 }
126
127 fclose(input);
128 Plot(energies, bin, eMin, eMax + 10);
129 }
130