advanced/doiPET/analysis.cpp

0001 //The following code is for anylizing axial sensitivity from coincidence list-mode data

0002 //It takes source position of the event (of those which are detected by the scanner) and analises axial sensitivty.

0003 //by Abdella M. Ahmed, 2020

0004
0005 #define _USE_MATH_DEFINES
0006 #include <iostream>
0007 #include <cfloat>
0008 #include <cmath>
0009 #include <fstream>
0010 #include <string>
0011 #include <limits>
0012 #include <stdio.h>
0013 #include <random>
0014 #include <string>
0015
0016 #define UseROOT //ASCII or UseROOT

0017
0018 #ifdef UseROOT
0019 //for root

0020 #include "TRandom3.h"
0021 #include "TFile.h"
0022 #include "TNtuple.h"
0023 #include "TROOT.h"
0024 #include "TH1.h"
0025 #include "TH2.h"
0026 #include "TTree.h"
0027 #include "TLeaf.h"
0028 #include "TSystem.h"
0029 #endif
0030
0031
0032
0033
0034 #define AxialLength 216 //(mm), Axial length of the scanner

0035
0036 //Change this number based of the number of primary particles simulated in the "run.mac" file

0037 #define NumberOfPrimaries 3000000 //Number of particles simulated

0038
0039 using namespace std;
0040 int main(){
0041
0042     cout<<"============Axial Sensitivity Analysis =========================="<<endl;
0043
0044     int eventID0, blockID0, crystalID_axial0, crystalID_tangential0, DOI_ID0;
0045     double timeStamp0, totalEdep0;
0046     int eventID1, blockID1, crystalID_axial1, crystalID_tangential1, DOI_ID1;
0047     double timeStamp1, totalEdep1;
0048     double  spositionX, spositionY, spositionZ; //source position

0049
0050     double z_offset = 0.0;//Axial offset position where the plane is located

0051     double planeWidth = 3;// (mm)

0052     int planeNumber;
0053
0054     float total_sensitivity = 0.0;
0055
0056     int numberOfSlices = int (AxialLength/planeWidth);
0057     cout<<"Number of axial planes (slices) are: " <<numberOfSlices<<endl;
0058     double Counts_per_plane[numberOfSlices];
0059
0060     ofstream OutFile("axial_sensitivity.csv");
0061     string filename = "resultCoincidence.data";
0062     string filepath = "";//provide the file path if it is stored in a different location.

0063     ifstream InFile(filepath+filename);
0064
0065     if (!InFile.is_open())
0066     {
0067         cout << "Unable to open input file to read .... " << endl;
0068         return 1;
0069     }
0070     if (!OutFile.is_open())
0071     {
0072         cout << "Unable to open out file to write ....";
0073         return 1;
0074     }
0075
0076     OutFile << "PlaneNmber" << "," << "Z(mm)" << "," << "Sensitivity(%)" << endl;
0077
0078     for (int i_plane = 0; i_plane < numberOfSlices; i_plane++){
0079         Counts_per_plane[i_plane] = 0;
0080     }
0081 #ifdef ASCII
0082     cout<<"\nASCII coincidence list-mode data is being analised..."<<endl;
0083     while (InFile >> eventID0 >> blockID0 >> crystalID_axial0 >> crystalID_tangential0 >> DOI_ID0 >> timeStamp0 >> totalEdep0 >>
0084                      eventID1 >> blockID1 >> crystalID_axial1 >> crystalID_tangential1 >> DOI_ID1 >> timeStamp1 >> totalEdep1 >>
0085                      spositionX >> spositionY >> spositionZ){
0086
0087         planeNumber = int(spositionZ/planeWidth + numberOfSlices/2 - 0.5 + 0.5);
0088         Counts_per_plane[planeNumber]++;
0089         total_sensitivity++;
0090     }
0091 #endif
0092
0093
0094 #ifdef UseROOT
0095     cout<<"\nROOT coincidence list-mode data is being analised..."<<endl;
0096     TFile *f = new TFile("resultCoincidence.root","READ");
0097     TTree *Singles = (TTree*)gDirectory->Get("Coincidence");
0098
0099     Singles->SetBranchAddress("eventID0",&eventID0);
0100     Singles->SetBranchAddress("blockID0",&blockID0);
0101     Singles->SetBranchAddress("crystalID_axial0",&crystalID_axial0);
0102     Singles->SetBranchAddress("crystalID_tangential0",&crystalID_tangential0);
0103     Singles->SetBranchAddress("DOI_ID0",&DOI_ID0);
0104     Singles->SetBranchAddress("timeStamp0",&timeStamp0);
0105     Singles->SetBranchAddress("totalEdep0",&totalEdep0);
0106
0107     Singles->SetBranchAddress("eventID1",&eventID1);
0108     Singles->SetBranchAddress("blockID1",&blockID1);
0109     Singles->SetBranchAddress("crystalID_axial1",&crystalID_axial1);
0110     Singles->SetBranchAddress("crystalID_tangential1",&crystalID_tangential1);
0111     Singles->SetBranchAddress("DOI_ID1",&DOI_ID1);
0112     Singles->SetBranchAddress("timeStamp1",&timeStamp1);
0113     Singles->SetBranchAddress("totalEdep1",&totalEdep1);
0114
0115     Singles->SetBranchAddress("spositionX",&spositionX);
0116     Singles->SetBranchAddress("spositionY",&spositionY);
0117     Singles->SetBranchAddress("spositionZ",&spositionZ);
0118
0119     //

0120     int nentries = 0;
0121     nentries = Singles->GetEntries();
0122     for(int entry=0; entry<nentries; entry++){
0123         Singles->GetEntry(entry);
0124         planeNumber = int(spositionZ/planeWidth + numberOfSlices/2 - 0.5 + 0.5);
0125         Counts_per_plane[planeNumber]++;
0126         total_sensitivity++;
0127     }
0128 #endif
0129
0130
0131     //Save it into CSV file

0132     for (int i_plane = 0; i_plane < numberOfSlices; i_plane++){
0133         //Axial mid-position of the plane or slice

0134         z_offset = (i_plane - numberOfSlices/2 + 0.5)*planeWidth;
0135
0136         OutFile << i_plane << "," << z_offset << "," << (Counts_per_plane[i_plane]/NumberOfPrimaries)*100 << endl;
0137     }
0138     cout<<"\nSensitivity evaluation has completed."<<endl;
0139     cout<<"\nTotal Sensitivity (Number of coinsidence per total number of pair of photons): "<<(total_sensitivity/NumberOfPrimaries)*100 << "%"<<endl;
0140     InFile.close();
0141     OutFile.close();
0142
0143     return 0;
0144 }