| Geant4 Cross Reference |
1 //******************************************** 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 pl
7 // simulation result GammaAtCreation.dat or Ga
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 e
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, d
51 {
52 auto mycanvas = new TCanvas("canvas", "canva
53 gPad->SetLeftMargin(0.15);
54
55 // unit is in keV
56 auto hist = new TH1D("hist (keV)", "Spectrum
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.da
74
75 if (input == NULL) {
76 printf("error for opening the input file\n
77 return;
78 }
79
80 //******************************************
81 //**************************Selection parame
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
89 int projection_index_end = 0; // end of the
90
91 int slice_index_begin = 0; // starter of th
92 int slice_index_end = 0; // end of the slic
93
94 //********************Parameters for spectru
95 int bin = 100;
96 double eMin = 0; // keV
97 double eMax = 0; // keV
98
99 //******************************************
100 //**************************Selection parame
101 //******************************************
102
103 RunInfo runInfo;
104 vector<double> energies;
105 int runID = -1; // index of simulations, na
106 // while(!feof(input)) //if not the end, rea
107 while (fread(&runInfo, sizeof(RunInfo), 1, i
108 runID++;
109 int nbParticle = runInfo.nbParticle;
110
111 // ***********the following codes are used
112 // if*************************************
113 // ***********the index of projection, sli
114 // configured in the simulation
115 runInfo.projectionIndex = runID / (nbSlice
116 int remain = runID % (nbSlice * nbPixel);
117 runInfo.sliceIndex = remain / nbPixel;
118 runInfo.pixelIndex = remain % nbPixel;
119 //****************************************
120
121 if (!nbParticle) continue;
122 std::vector<ParticleInfo> particles(nbPart
123 fread(&particles[0], sizeof(ParticleInfo),
124
125 if (runInfo.projectionIndex >= projection_
126 && runInfo.projectionIndex <= projecti
127 {
128 if (runInfo.sliceIndex >= slice_index_be
129 for (int i = 0; i < nbParticle; ++i) {
130 // printf("--%d, %.9e\n", i, particl
131
132 energies.push_back(particles[i].ener
133 if (particles[i].energy_keV > eMax)
134 }
135 }
136 }
137 else
138 break;
139 }
140
141 fclose(input);
142 Plot(energies, bin, eMin, eMax + 10);
143 }
144