doiPET/src/doiPETRun.cc

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 /// \file electromagnetic/TestEm11/src/doiPETRun.cc
0027 /// \brief Implementation of the doiPETRun class
0028 //
0029 // $Id: doiPETRun.cc 71376 2013-06-14 07:44:50Z maire $
0030 //
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0033
0034 #include "doiPETRun.hh"
0035 #include "doiPETDetectorConstruction.hh"
0036 #include "doiPETPrimaryGeneratorAction.hh"
0037
0038 #include "G4UnitsTable.hh"
0039 #include "G4SystemOfUnits.hh"
0040 #include "G4Event.hh"
0041 #include "doiPETAnalysis.hh"
0042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0043
0044 doiPETRun::doiPETRun()
0045     : G4Run()
0046 {
0047     //
0048     totalTime = 0 * s;
0049     prev_totalTime = 0 * s;
0050     totalEdep = 0*keV;
0051 }
0052
0053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0054
0055 doiPETRun::~doiPETRun()
0056 {}
0057
0058 //
0059 void doiPETRun::GetIntractionInfomation(InteractionInformation*step)
0060 {
0061     //The photon can have many interactions in one block. This stores all the interaction (including multiple copies of the same block.)
0062     //Save the interacting block in the associative container including multiple copies. Since there will be multiple interactions in the crystals wthin a single block
0063     mapBlockInteraction.insert( std::multimap< G4int, InteractionInformation* >::value_type( step->GetBlockNo(),  step) );
0064
0065     //The following stores the block number uniquely. Save the interacting block in the associative container without saving multiple copies
0066     setBlockInteraction.insert( step->GetBlockNo());
0067 }
0068
0069 ///////// InteractingCrystal ///////////////////////////////////
0070
0071 void doiPETRun::FindInteractingCrystal()
0072 {
0073     std::multimap< G4int, InteractionInformation* >::iterator mitr;
0074     std::set<G4int>::iterator sitr;
0075     G4int cystalIDTemp;
0076
0077     for(sitr=setBlockInteraction.begin(); sitr!=setBlockInteraction.end(); sitr++){
0078         blockID = (*sitr);
0079         mitr = mapBlockInteraction.find(blockID);
0080         numberofInteractions = mapBlockInteraction.count(blockID); //number of interactions in which the interaction deposite energy in a block.
0081         totalEdep=edepMax=0.;
0082         for(G4int i=0; i<numberofInteractions; i++, mitr++){
0083             edep = (*mitr).second->GetEdep();
0084             cystalIDTemp = (*mitr).second->GetCrystalNo();
0085             edep_AfterCrystalBlurring = doiPETAnalysis::GetInstance()->QuantumEffifciency(edep, blockID, cystalIDTemp);
0086             totalEdep += edep_AfterCrystalBlurring;
0087
0088             crystalID_vec.push_back(cystalIDTemp);
0089             posInter_vec.push_back(((*mitr).second->GetInteractionPositionInCrystal()));
0090             edepInCry_vec.push_back(edep_AfterCrystalBlurring);
0091
0092             //This is  to identify the crystal with the highest energy deposition.
0093             if(edepMax<edep_AfterCrystalBlurring){
0094                 edepMax=edep_AfterCrystalBlurring;
0095                 crystalID = cystalIDTemp;
0096                 interactionTime = (*mitr).second->GetGlobalTime();
0097
0098                 //position of interaction based on highest edep. See below for energy weighted position of interaction
0099                 //interactionPos = (*mitr).second->GetInteractionPositionInCrystal();
0100             }
0101         }
0102         if(totalEdep==0)continue;
0103
0104         //The interaction position is based on energy weighted center of mass calculation
0105         interactionPos = CenterOfMassInteractionPos(crystalID_vec, edepInCry_vec, totalEdep, posInter_vec);
0106         doiPETAnalysis::GetInstance()->ReadOut(blockID,crystalID,interactionTime,timeOfAnnihil,interactionPos,totalEdep);
0107
0108         crystalID_vec.clear();
0109         posInter_vec.clear();
0110         edepInCry_vec.clear();
0111     }
0112
0113 }
0114
0115 //Apply center of mass calculation in determining the interaction position
0116 G4ThreeVector doiPETRun::CenterOfMassInteractionPos(const std::vector<G4int> &cryID, const std::vector<G4double> &energyDep, G4double edepTotal, const std::vector<G4ThreeVector> &pos){
0117     posInterInCrystal = G4ThreeVector();
0118     for(std::size_t i=0; i<cryID.size();i++){
0119         posInterInCrystal += pos[i] * energyDep[i];
0120     }
0121     return (posInterInCrystal)/edepTotal;
0122 }
0123
0124
0125 void doiPETRun::SetAnnihilationTime(G4double t){
0126     timeOfAnnihil = t;
0127 }
0128 //
0129 void doiPETRun::Clear(){
0130     //clear the content of the associative containers, the multimap and the set
0131     std::multimap< G4int, InteractionInformation* >::iterator mitr2;
0132     for( mitr2= mapBlockInteraction.begin(); mitr2!=mapBlockInteraction.end(); mitr2++){
0133         delete (*mitr2).second;
0134     }
0135
0136     mapBlockInteraction.clear();
0137     setBlockInteraction.clear();
0138     //
0139 }
0140 //
0141 void doiPETRun::GetEventIDRun(G4int evt){
0142     eventID = evt;
0143 }
0144 //
0145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0146
0147 void doiPETRun::Merge(const G4Run* run)
0148 {
0149     G4Run::Merge(run);
0150 }