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