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