doiPET/include/doiPETAnalysis.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025
0026 //GEANT4 - Depth-of-Interaction enabled Positron emission tomography (PET) advanced example
0027
0028 //Authors and contributors
0029
0030 // Author list to be updated, with names of co-authors and contributors from National Institute of Radiological Sciences (NIRS)
0031
0032 // Abdella M. Ahmed (1, 2), Andrew Chacon (1, 2), Harley Rutherford (1, 2),
0033 // Hideaki Tashima (3), Go Akamatsu (3), Akram Mohammadi (3), Eiji Yoshida (3), Taiga Yamaya (3)
0034 // Susanna Guatelli (2), and Mitra Safavi-Naeini (1, 2)
0035
0036 // (1) Australian Nuclear Science and Technology Organisation, Australia
0037 // (2) University of Wollongong, Australia
0038 // (3) National Institute of Radiological Sciences, Japan
0039
0040
0041 #ifndef doiPETAnalysis_h
0042 #define doiPETAnalysis_h  1
0043
0044 #include "doiPETGlobalParameters.hh"
0045 #include "globals.hh"
0046 #include <vector>
0047 #include <time.h>
0048 #include <map>
0049 #include <set>
0050 #include "G4ThreeVector.hh"
0051 #include <iostream>
0052 #include <fstream>
0053 #include <sstream>
0054 #include <iterator>
0055 #include <vector>
0056 #include <algorithm>
0057 #include "G4AnalysisManager.hh"
0058
0059 // Define the total number of columns in the ntuple
0060 const G4int MaxNtCol = 17;
0061
0062 class doiPETAnalysisMessenger;
0063
0064 //class InteractionInformation;
0065
0066 class doiPETAnalysis
0067 {
0068 private:
0069     doiPETAnalysis();
0070
0071 public:
0072     ~doiPETAnalysis();
0073     static doiPETAnalysis* GetInstance();
0074     void FindInteractingCrystal();
0075     void Open(G4String);
0076     void Close();
0077     void Delete();
0078     void ResetNumberOfHits();
0079     void Write(/*G4int, G4int, G4int, G4double*/);
0080     void WriteOutput();
0081
0082     //void GetIntractionInfomation(InteractionInformation*);
0083
0084     void GetParentParticleName(G4String);
0085     void GetSizeOfDetector (G4double, G4double, G4double);
0086     void SetScatterIndexInPhantom(G4int);
0087
0088     void SetSourcePosition(G4ThreeVector);//
0089     void SetEventID(G4int);
0090
0091     void BlurringParameters();
0092     void GetTimeOfAnnihilation(G4double);
0093
0094
0095     void PMTPosition();
0096     void AngerLogic(G4double, G4double, G4double, G4double, G4double, G4bool);
0097     void ReadReflectorPattern();
0098     void PrepareDOILookUpTable(G4String);
0099
0100     void SetActivity(G4double);
0101     void SetIsotopeHalfLife(G4double);
0102     void CrystalIDAfterAngerLogic(G4int, G4int, G4int);
0103     void TypeOfOutput(G4String);//Single or coincidence list-mode data
0104     void CalulateAcquisitionTime();
0105     //G4double QuantumEffifciency(G4double);
0106     G4double QuantumEffifciency(G4double, G4int, G4int);
0107     void ReadOut(G4int, G4int, G4double, G4double, G4ThreeVector, G4double);
0108
0109     //G4ROOT
0110     void book(); // booking the ROOT file
0111
0112
0113     void FillListModeEvent(); //Single or Coinsidence
0114     void finish();
0115     // Close the ROOT file with all the results stored in nutples
0116
0117
0118 private:
0119     static doiPETAnalysis* instance;
0120     doiPETAnalysisMessenger* fAnalysisMessenger;
0121     //std::multimap< G4int, InteractionInformation* > mapBlockInteraction;
0122     std::set<G4int> setBlockInteraction;
0123
0124     G4double upperThreshold, lowerThreshold;
0125     G4double triggerEnergy;
0126
0127     //G4ROOT
0128     G4bool factoryOn;
0129     G4int  fNtColId[MaxNtCol];
0130   //
0131
0132
0133     //G4ThreeVector sourcePosition;
0134
0135
0136
0137     //
0138     G4int scatterIndex;
0139
0140     G4String parentParticleName;//
0141
0142     //
0143     G4int numberofInteractions;
0144     G4int countCoincidence;
0145
0146     G4int numberOfBlocks_total;
0147
0148     G4double sizeOfDetector_DOI,sizeOfDetector_axial,sizeOfDetector_tangential;
0149
0150     //Virtual position of the PMT
0151     G4double signalPMT1, signalPMT2, signalPMT3, signalPMT4;
0152
0153     G4double posPMT1x, posPMT2x, posPMT3x, posPMT4x;
0154     G4double posPMT1y, posPMT2y, posPMT3y, posPMT4y;
0155     G4double posPMT1z, posPMT2z, posPMT3z, posPMT4z;
0156
0157     //
0158     G4double signalPMT1z, signalPMT2z, signalPMT3z, signalPMT4z;
0159     G4double signalPMT1y, signalPMT2y, signalPMT3y, signalPMT4y;
0160
0161     //
0162     G4double signalZplus, signalZminus;
0163     G4double signalYplus, signalYminus;
0164     //
0165
0166     G4double dist1z, dist2z, dist3z, dist4z, distz;
0167     G4double dist1y, dist2y, dist3y, dist4y, disty;
0168
0169     G4double shiftCoeff;
0170
0171     G4double PositionAngerZ, PositionAngerY;
0172
0173     //reflector pattern
0174     std::vector<G4int> ireflectorLayer1_Tangential;
0175     std::vector<G4int> ireflectorLayer1_Axial;
0176     std::vector<G4int> ireflectorLayer2_Tangential;
0177     std::vector<G4int> ireflectorLayer2_Axial;
0178     std::vector<G4int> ireflectorLayer3_Tangential;
0179     std::vector<G4int> ireflectorLayer3_Axial;
0180     std::vector<G4int> ireflectorLayer4_Tangential;
0181     std::vector<G4int> ireflectorLayer4_Axial;
0182     std::vector<G4int> doi_table;
0183     //
0184
0185     //The number of pixes for the 2D position histogram after Anger Logic calculation
0186     G4int numberOfPixel_axial;
0187     G4int numberOfPixel_tan;
0188
0189     //source position
0190     G4double spositionX;
0191     G4double spositionY;
0192     G4double spositionZ;
0193
0194     //interaction position with respect to the crystal axis
0195     G4ThreeVector interactionPos;
0196
0197     //interaction position
0198     G4double intPosX;
0199     G4double intPosY;
0200     G4double intPosZ;
0201
0202
0203     G4double interactionTime;
0204
0205     G4int crystalID;//contineous crystal ID in 3D
0206     G4int crystalID_2D;
0207
0208     G4int prev_eventID;
0209
0210     //Single output
0211     G4int eventID;
0212     G4int blockID;
0213     G4int crystalID_axial;
0214     G4int crystalID_tangential;
0215     G4int DOI_ID;
0216     G4double timeStamp;
0217     G4double totalEdep;
0218
0219     //coincidence output
0220     G4int eventID0,                 eventID1;
0221     G4int blockID0,                 blockID1;
0222     G4int crystalID_axial0,         crystalID_axial1;
0223     G4int crystalID_tangential0,    crystalID_tangential1;
0224     G4int DOI_ID0,                  DOI_ID1;
0225     G4double timeStamp0,            timeStamp1;
0226     G4double totalEdep0,            totalEdep1;
0227     G4double sposX, sposY, sposZ;
0228     //choice for the user
0229     G4bool getSinglesData;
0230     G4bool getCoincidenceData;
0231
0232     //
0233     G4bool ApplyAngerLogic;
0234
0235     G4double PMTblurring_tan;
0236     G4double PMTblurring_axial;
0237
0238     G4String outputData;
0239     G4int numberOfHit;
0240     std::vector<G4int> eventID_coin;
0241     std::vector<G4double> edep_coin;
0242     std::vector<G4int>blockID_coin;
0243     std::vector<G4int> cryID_axial_coin;
0244     std::vector<G4int> cryID_tan_coin;
0245     std::vector<G4int> cryDOI_coin;
0246     std::vector<G4double> time_coin;
0247
0248     //Crystal IDs after Anger Logic calculation
0249     G4int crystalIDNew_DOI, crystalIDNew_tan, crystalIDNew_axial;
0250
0251     //Crystal ID in the 2D position histogram along the axial and tangetial direction
0252     G4int crystalID_in2D_posHist_axial, crystalID_in2D_posHist_tan;
0253
0254     //continous crystal ID after after Anger Logic.
0255     G4int crystalID_in2D_posHist;
0256
0257
0258     //Crystal blurring
0259     G4double crystalResolution;
0260     G4double crystalResolutionMin;//
0261     G4double crystalResolutionMax;//
0262
0263     //G4bool variableResolution;
0264     G4bool fixedResolution;
0265     G4bool isDOIlookUpTablePrepared;
0266
0267     G4double energyResolution_fixed;
0268     std::vector<std::vector<G4double>> energyResolution_cryDependent;
0269
0270     G4double crystalEnergyRef;//This 511 keV
0271     G4double crystalQuantumEfficiency;//
0272     G4double edep_AfterCrystalBlurring;
0273     G4double crystalCoeff;
0274     G4double sigma_energyResolution;
0275
0276     G4double totalTime;
0277     G4double prev_totalTime;
0278     G4double timeInterval;
0279     G4double time_annihil;
0280     G4double time_tof;
0281     G4double block_DeadTime;
0282     G4double module_DeadTime;
0283
0284     G4double *blockTime;
0285     G4double *moduleTime;
0286
0287     //
0288     G4double activityNow;
0289     G4double InitialActivity;
0290     G4double halfLife;
0291
0292     G4String simulationType;
0293
0294     //
0295     //for output file to write results
0296     std::ofstream ofs;
0297     G4String asciiFileName;
0298     G4String rootFileName;
0299
0300 // #ifdef USEROOT
0301 //     TFile* file;
0302 //  TTree* tSingles;
0303 //  TTree* tCoincidence;
0304 //  //TH1F*hb;
0305 // #endif
0306
0307     //input file to read reflector pattern
0308     std::ifstream ifs;
0309 };
0310
0311 #endif
0312