| Geant4 Cross Reference |
1 //The following code is for anylizing axial sensitivity from coincidence list-mode data
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
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 of the scanner
35
36 //Change this number based of the number of primary particles simulated in the "run.mac" file
37 #define NumberOfPrimaries 3000000 //Number of particles simulated
38
39 using namespace std;
40 int main(){
41
42 cout<<"============Axial Sensitivity Analysis =========================="<<endl;
43
44 int eventID0, blockID0, crystalID_axial0, crystalID_tangential0, DOI_ID0;
45 double timeStamp0, totalEdep0;
46 int eventID1, blockID1, crystalID_axial1, crystalID_tangential1, DOI_ID1;
47 double timeStamp1, totalEdep1;
48 double spositionX, spositionY, spositionZ; //source position
49
50 double z_offset = 0.0;//Axial offset position where the plane is located
51 double planeWidth = 3;// (mm)
52 int planeNumber;
53
54 float total_sensitivity = 0.0;
55
56 int numberOfSlices = int (AxialLength/planeWidth);
57 cout<<"Number of axial planes (slices) are: " <<numberOfSlices<<endl;
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 if it is stored in a different location.
63 ifstream InFile(filepath+filename);
64
65 if (!InFile.is_open())
66 {
67 cout << "Unable to open input file to read .... " << endl;
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)" << "," << "Sensitivity(%)" << endl;
77
78 for (int i_plane = 0; i_plane < numberOfSlices; i_plane++){
79 Counts_per_plane[i_plane] = 0;
80 }
81 #ifdef ASCII
82 cout<<"\nASCII coincidence list-mode data is being analised..."<<endl;
83 while (InFile >> eventID0 >> blockID0 >> crystalID_axial0 >> crystalID_tangential0 >> DOI_ID0 >> timeStamp0 >> totalEdep0 >>
84 eventID1 >> blockID1 >> crystalID_axial1 >> crystalID_tangential1 >> DOI_ID1 >> timeStamp1 >> totalEdep1 >>
85 spositionX >> spositionY >> spositionZ){
86
87 planeNumber = int(spositionZ/planeWidth + numberOfSlices/2 - 0.5 + 0.5);
88 Counts_per_plane[planeNumber]++;
89 total_sensitivity++;
90 }
91 #endif
92
93
94 #ifdef UseROOT
95 cout<<"\nROOT coincidence list-mode data is being analised..."<<endl;
96 TFile *f = new TFile("resultCoincidence.root","READ");
97 TTree *Singles = (TTree*)gDirectory->Get("Coincidence");
98
99 Singles->SetBranchAddress("eventID0",&eventID0);
100 Singles->SetBranchAddress("blockID0",&blockID0);
101 Singles->SetBranchAddress("crystalID_axial0",&crystalID_axial0);
102 Singles->SetBranchAddress("crystalID_tangential0",&crystalID_tangential0);
103 Singles->SetBranchAddress("DOI_ID0",&DOI_ID0);
104 Singles->SetBranchAddress("timeStamp0",&timeStamp0);
105 Singles->SetBranchAddress("totalEdep0",&totalEdep0);
106
107 Singles->SetBranchAddress("eventID1",&eventID1);
108 Singles->SetBranchAddress("blockID1",&blockID1);
109 Singles->SetBranchAddress("crystalID_axial1",&crystalID_axial1);
110 Singles->SetBranchAddress("crystalID_tangential1",&crystalID_tangential1);
111 Singles->SetBranchAddress("DOI_ID1",&DOI_ID1);
112 Singles->SetBranchAddress("timeStamp1",&timeStamp1);
113 Singles->SetBranchAddress("totalEdep1",&totalEdep1);
114
115 Singles->SetBranchAddress("spositionX",&spositionX);
116 Singles->SetBranchAddress("spositionY",&spositionY);
117 Singles->SetBranchAddress("spositionZ",&spositionZ);
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 + numberOfSlices/2 - 0.5 + 0.5);
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 < numberOfSlices; i_plane++){
133 //Axial mid-position of the plane or slice
134 z_offset = (i_plane - numberOfSlices/2 + 0.5)*planeWidth;
135
136 OutFile << i_plane << "," << z_offset << "," << (Counts_per_plane[i_plane]/NumberOfPrimaries)*100 << endl;
137 }
138 cout<<"\nSensitivity evaluation has completed."<<endl;
139 cout<<"\nTotal Sensitivity (Number of coinsidence per total number of pair of photons): "<<(total_sensitivity/NumberOfPrimaries)*100 << "%"<<endl;
140 InFile.close();
141 OutFile.close();
142
143 return 0;
144 }
145