| Geant4 Cross Reference |
1 //The following code is for anylizing axial se 1
2 //It takes source position of the event (of th
3 //by Abdella M. Ahmed, 2020
4
5 #define _USE_MATH_DEFINES
6 #include <iostream>
7 #include <cfloat>
8 #include <cmath>
9 #include <fstream>
10 #include <string>
11 #include <limits>
12 #include <stdio.h>
13 #include <random>
14 #include <string>
15
16 #define UseROOT //ASCII or UseROOT
17
18 #ifdef UseROOT
19 //for root
20 #include "TRandom3.h"
21 #include "TFile.h"
22 #include "TNtuple.h"
23 #include "TROOT.h"
24 #include "TH1.h"
25 #include "TH2.h"
26 #include "TTree.h"
27 #include "TLeaf.h"
28 #include "TSystem.h"
29 #endif
30
31
32
33
34 #define AxialLength 216 //(mm), Axial length o
35
36 //Change this number based of the number of pr
37 #define NumberOfPrimaries 3000000 //Number of
38
39 using namespace std;
40 int main(){
41
42 cout<<"============Axial Sensitivity Analysi
43
44 int eventID0, blockID0, crystalID_axial0, cr
45 double timeStamp0, totalEdep0;
46 int eventID1, blockID1, crystalID_axial1,
47 double timeStamp1, totalEdep1;
48 double spositionX, spositionY, spositionZ;
49
50 double z_offset = 0.0;//Axial offset positio
51 double planeWidth = 3;// (mm)
52 int planeNumber;
53
54 float total_sensitivity = 0.0;
55
56 int numberOfSlices = int (AxialLength/planeW
57 cout<<"Number of axial planes (slices) are:
58 double Counts_per_plane[numberOfSlices];
59
60 ofstream OutFile("axial_sensitivity.csv");
61 string filename = "resultCoincidence.data";
62 string filepath = "";//provide the file path
63 ifstream InFile(filepath+filename);
64
65 if (!InFile.is_open())
66 {
67 cout << "Unable to open input file to read
68 return 1;
69 }
70 if (!OutFile.is_open())
71 {
72 cout << "Unable to open out file to write
73 return 1;
74 }
75
76 OutFile << "PlaneNmber" << "," << "Z(mm)" <<
77
78 for (int i_plane = 0; i_plane < numberOfSlic
79 Counts_per_plane[i_plane] = 0;
80 }
81 #ifdef ASCII
82 cout<<"\nASCII coincidence list-mode data
83 while (InFile >> eventID0 >> blockID0 >> cry
84 eventID1 >> blockID1 >> cryst
85 spositionX >> spositionY >> spositi
86
87 planeNumber = int(spositionZ/planeWidth +
88 Counts_per_plane[planeNumber]++;
89 total_sensitivity++;
90 }
91 #endif
92
93
94 #ifdef UseROOT
95 cout<<"\nROOT coincidence list-mode data is
96 TFile *f = new TFile("resultCoincidence.root
97 TTree *Singles = (TTree*)gDirectory->Get("Co
98
99 Singles->SetBranchAddress("eventID0",&eventI
100 Singles->SetBranchAddress("blockID0",&blockI
101 Singles->SetBranchAddress("crystalID_axial0"
102 Singles->SetBranchAddress("crystalID_tangent
103 Singles->SetBranchAddress("DOI_ID0",&DOI_ID0
104 Singles->SetBranchAddress("timeStamp0",&time
105 Singles->SetBranchAddress("totalEdep0",&tota
106
107 Singles->SetBranchAddress("eventID1",&eventI
108 Singles->SetBranchAddress("blockID1",&blockI
109 Singles->SetBranchAddress("crystalID_axial1"
110 Singles->SetBranchAddress("crystalID_tangent
111 Singles->SetBranchAddress("DOI_ID1",&DOI_ID1
112 Singles->SetBranchAddress("timeStamp1",&time
113 Singles->SetBranchAddress("totalEdep1",&tota
114
115 Singles->SetBranchAddress("spositionX",&spos
116 Singles->SetBranchAddress("spositionY",&spos
117 Singles->SetBranchAddress("spositionZ",&spos
118
119 //
120 int nentries = 0;
121 nentries = Singles->GetEntries();
122 for(int entry=0; entry<nentries; entry++){
123 Singles->GetEntry(entry);
124 planeNumber = int(spositionZ/planeWidth +
125 Counts_per_plane[planeNumber]++;
126 total_sensitivity++;
127 }
128 #endif
129
130
131 //Save it into CSV file
132 for (int i_plane = 0; i_plane < numberOfSlic
133 //Axial mid-position of the plane or slice
134 z_offset = (i_plane - numberOfSlices/2 + 0
135
136 OutFile << i_plane << "," << z_offset << "
137 }
138 cout<<"\nSensitivity evaluation has complete
139 cout<<"\nTotal Sensitivity (Number of coinsi
140 InFile.close();
141 OutFile.close();
142
143 return 0;
144 }
145