| Geant4 Cross Reference |
1 //***********************************************************************************************************
2 // Spectrum_gamma.C
3 // Root command file
4 // Type: root Spectrum_gamma.C
5 //
6 // It visualizes the spectrum of X-rays and plots a histogram by reading
7 // simulation result GammaAtCreation.dat or GammaAtExit.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 // struct ParticleInfo
40 //{
41 // float energy_keV;
42 // float mx;
43 // float my;
44 // float mz;
45 // float x;
46 // float y;
47 // float z;
48 //};
49
50 void Plot(vector<double>& energies, int bin, double eMin, double eMax)
51 {
52 auto mycanvas = new TCanvas("canvas", "canvas", 800, 50, 600, 600);
53 gPad->SetLeftMargin(0.15);
54
55 // unit is in keV
56 auto hist = new TH1D("hist (keV)", "Spectrum of photons", bin, eMin, eMax);
57
58 for (int i = 0; i < energies.size(); ++i) {
59 hist->Fill(energies[i]);
60 }
61
62 hist->Draw();
63 hist->GetXaxis()->SetTitle("Energy (keV)");
64 hist->GetYaxis()->SetTitle("Counts");
65 hist->GetXaxis()->CenterTitle();
66 hist->GetYaxis()->CenterTitle();
67
68 mycanvas->Print("spectrum_gamma.png");
69 }
70
71 void Spectrum_gamma()
72 {
73 FILE* input = fopen("../build/GammaAtExit.dat", "rb");
74
75 if (input == NULL) {
76 printf("error for opening the input file\n");
77 return;
78 }
79
80 //***********************************************************************
81 //**************************Selection parameters (begin)*****************
82 //***********************************************************************
83
84 const int nbProjection = 10;
85 const int nbSlice = 1;
86 const int nbPixel = 20;
87
88 int projection_index_begin = 0; // starter of the projection selected
89 int projection_index_end = 0; // end of the projection selected
90
91 int slice_index_begin = 0; // starter of the slice selected
92 int slice_index_end = 0; // end of the slice selected
93
94 //********************Parameters for spectrum***************************
95 int bin = 100;
96 double eMin = 0; // keV
97 double eMax = 0; // keV
98
99 //***********************************************************************
100 //**************************Selection parameters (end)*******************
101 //***********************************************************************
102
103 RunInfo runInfo;
104 vector<double> energies;
105 int runID = -1; // index of simulations, namely runID, starting from 0
106 // while(!feof(input)) //if not the end, read
107 while (fread(&runInfo, sizeof(RunInfo), 1, input)) {
108 runID++;
109 int nbParticle = runInfo.nbParticle;
110
111 // ***********the following codes are used
112 // if**************************************(begin)
113 // ***********the index of projection, slice and pixel is not correctly
114 // configured in the simulation
115 runInfo.projectionIndex = runID / (nbSlice * nbPixel);
116 int remain = runID % (nbSlice * nbPixel);
117 runInfo.sliceIndex = remain / nbPixel;
118 runInfo.pixelIndex = remain % nbPixel;
119 //************************************************************************(end)
120
121 if (!nbParticle) continue;
122 std::vector<ParticleInfo> particles(nbParticle);
123 fread(&particles[0], sizeof(ParticleInfo), nbParticle, input);
124
125 if (runInfo.projectionIndex >= projection_index_begin
126 && runInfo.projectionIndex <= projection_index_end)
127 {
128 if (runInfo.sliceIndex >= slice_index_begin && runInfo.sliceIndex <= slice_index_end) {
129 for (int i = 0; i < nbParticle; ++i) {
130 // printf("--%d, %.9e\n", i, particles[i].energy_keV);
131
132 energies.push_back(particles[i].energy_keV);
133 if (particles[i].energy_keV > eMax) eMax = particles[i].energy_keV;
134 }
135 }
136 }
137 else
138 break;
139 }
140
141 fclose(input);
142 Plot(energies, bin, eMin, eMax + 10);
143 }
144