LXe/src/LXeEventAction.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 LXeEventAction.cc
0027 /// \brief Implementation of the LXeEventAction class
0028
0029 #include "LXeEventAction.hh"
0030
0031 #include "LXeDetectorConstruction.hh"
0032 #include "LXeHistoManager.hh"
0033 #include "LXePMTHit.hh"
0034 #include "LXeRun.hh"
0035 #include "LXeScintHit.hh"
0036 #include "LXeTrajectory.hh"
0037
0038 #include "G4Event.hh"
0039 #include "G4EventManager.hh"
0040 #include "G4RunManager.hh"
0041 #include "G4SDManager.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4Trajectory.hh"
0044 #include "G4TrajectoryContainer.hh"
0045 #include "G4UImanager.hh"
0046 #include "G4VVisManager.hh"
0047 #include "G4ios.hh"
0048
0049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0050
0051 LXeEventAction::LXeEventAction(const LXeDetectorConstruction* det) : fDetector(det)
0052 {
0053   fEventMessenger = new LXeEventMessenger(this);
0054 }
0055
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057
0058 LXeEventAction::~LXeEventAction()
0059 {
0060   delete fEventMessenger;
0061 }
0062
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064
0065 void LXeEventAction::BeginOfEventAction(const G4Event*)
0066 {
0067   fHitCount = 0;
0068   fPhotonCount_Scint = 0;
0069   fPhotonCount_Ceren = 0;
0070   fAbsorptionCount = 0;
0071   fBoundaryAbsorptionCount = 0;
0072   fTotE = 0.0;
0073
0074   fConvPosSet = false;
0075   fEdepMax = 0.0;
0076
0077   fPMTsAboveThreshold = 0;
0078
0079   G4SDManager* SDman = G4SDManager::GetSDMpointer();
0080   if (fScintCollID < 0) fScintCollID = SDman->GetCollectionID("scintCollection");
0081   if (fPMTCollID < 0) fPMTCollID = SDman->GetCollectionID("pmtHitCollection");
0082 }
0083
0084 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0085
0086 void LXeEventAction::EndOfEventAction(const G4Event* anEvent)
0087 {
0088   G4TrajectoryContainer* trajectoryContainer = anEvent->GetTrajectoryContainer();
0089
0090   G4int n_trajectories = 0;
0091   if (trajectoryContainer) n_trajectories = trajectoryContainer->entries();
0092
0093   // extract the trajectories and draw them
0094   if (G4VVisManager::GetConcreteInstance()) {
0095     for (G4int i = 0; i < n_trajectories; ++i) {
0096       auto trj = (LXeTrajectory*)((*(anEvent->GetTrajectoryContainer()))[i]);
0097       if (trj->GetParticleName() == "opticalphoton") {
0098         trj->SetForceDrawTrajectory(fForcedrawphotons);
0099         trj->SetForceNoDrawTrajectory(fForcenophotons);
0100       }
0101       trj->DrawTrajectory();
0102     }
0103   }
0104
0105   LXeScintHitsCollection* scintHC = nullptr;
0106   LXePMTHitsCollection* pmtHC = nullptr;
0107   G4HCofThisEvent* hitsCE = anEvent->GetHCofThisEvent();
0108
0109   // Get the hit collections
0110   if (hitsCE) {
0111     if (fScintCollID >= 0) {
0112       scintHC = (LXeScintHitsCollection*)(hitsCE->GetHC(fScintCollID));
0113     }
0114     if (fPMTCollID >= 0) {
0115       pmtHC = (LXePMTHitsCollection*)(hitsCE->GetHC(fPMTCollID));
0116     }
0117   }
0118
0119   // Hits in scintillator
0120   if (scintHC) {
0121     size_t n_hit = scintHC->entries();
0122     G4ThreeVector eWeightPos(0.);
0123     G4double edep;
0124     G4double edepMax = 0;
0125
0126     for (size_t i = 0; i < n_hit; ++i) {  // gather info on hits in scintillator
0127       edep = (*scintHC)[i]->GetEdep();
0128       fTotE += edep;
0129       eWeightPos += (*scintHC)[i]->GetPos() * edep;  // calculate energy weighted pos
0130       if (edep > edepMax) {
0131         edepMax = edep;  // store max energy deposit
0132         G4ThreeVector posMax = (*scintHC)[i]->GetPos();
0133         fPosMax = posMax;
0134         fEdepMax = edep;
0135       }
0136     }
0137
0138     G4AnalysisManager::Instance()->FillH1(7, fTotE);
0139
0140     if (fTotE == 0.) {
0141       if (fVerbose > 0) G4cout << "No hits in the scintillator this event." << G4endl;
0142     }
0143     else {
0144       // Finish calculation of energy weighted position
0145       eWeightPos /= fTotE;
0146       fEWeightPos = eWeightPos;
0147       if (fVerbose > 0) {
0148         G4cout << "\tEnergy weighted position of hits in LXe : " << eWeightPos / mm << G4endl;
0149       }
0150     }
0151     if (fVerbose > 0) {
0152       G4cout << "\tTotal energy deposition in scintillator : " << fTotE / keV << " (keV)" << G4endl;
0153     }
0154   }
0155
0156   if (pmtHC) {
0157     G4ThreeVector reconPos(0., 0., 0.);
0158     size_t pmts = pmtHC->entries();
0159     // Gather info from all PMTs
0160     for (size_t i = 0; i < pmts; ++i) {
0161       fHitCount += (*pmtHC)[i]->GetPhotonCount();
0162       reconPos += (*pmtHC)[i]->GetPMTPos() * (*pmtHC)[i]->GetPhotonCount();
0163       if ((*pmtHC)[i]->GetPhotonCount() >= fPMTThreshold) {
0164         ++fPMTsAboveThreshold;
0165       }
0166       else {  // wasn't above the threshold, turn it back off
0167         (*pmtHC)[i]->SetDrawit(false);
0168       }
0169     }
0170
0171     G4AnalysisManager::Instance()->FillH1(1, fHitCount);
0172     G4AnalysisManager::Instance()->FillH1(2, fPMTsAboveThreshold);
0173
0174     if (fHitCount > 0) {  // don't bother unless there were hits
0175       reconPos /= fHitCount;
0176       if (fVerbose > 0) {
0177         G4cout << "\tReconstructed position of hits in LXe : " << reconPos / mm << G4endl;
0178       }
0179       fReconPos = reconPos;
0180     }
0181     pmtHC->DrawAllHits();
0182   }
0183
0184   G4AnalysisManager::Instance()->FillH1(3, fPhotonCount_Scint);
0185   G4AnalysisManager::Instance()->FillH1(4, fPhotonCount_Ceren);
0186   G4AnalysisManager::Instance()->FillH1(5, fAbsorptionCount);
0187   G4AnalysisManager::Instance()->FillH1(6, fBoundaryAbsorptionCount);
0188
0189   if (fVerbose > 0) {
0190     // End of event output. later to be controlled by a verbose level
0191     G4cout << "\tNumber of photons that hit PMTs in this event : " << fHitCount << G4endl;
0192     G4cout << "\tNumber of PMTs above threshold(" << fPMTThreshold << ") : " << fPMTsAboveThreshold
0193            << G4endl;
0194     G4cout << "\tNumber of photons produced by scintillation in this event : " << fPhotonCount_Scint
0195            << G4endl;
0196     G4cout << "\tNumber of photons produced by cerenkov in this event : " << fPhotonCount_Ceren
0197            << G4endl;
0198     G4cout << "\tNumber of photons absorbed (OpAbsorption) in this event : " << fAbsorptionCount
0199            << G4endl;
0200     G4cout << "\tNumber of photons absorbed at boundaries (OpBoundary) in "
0201            << "this event : " << fBoundaryAbsorptionCount << G4endl;
0202     G4cout << "Unaccounted for photons in this event : "
0203            << (fPhotonCount_Scint + fPhotonCount_Ceren - fAbsorptionCount - fHitCount
0204                - fBoundaryAbsorptionCount)
0205            << G4endl;
0206   }
0207
0208   // update the run statistics
0209   auto run = static_cast<LXeRun*>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
0210
0211   run->IncHitCount(fHitCount);
0212   run->IncPhotonCount_Scint(fPhotonCount_Scint);
0213   run->IncPhotonCount_Ceren(fPhotonCount_Ceren);
0214   run->IncEDep(fTotE);
0215   run->IncAbsorption(fAbsorptionCount);
0216   run->IncBoundaryAbsorption(fBoundaryAbsorptionCount);
0217   run->IncHitsAboveThreshold(fPMTsAboveThreshold);
0218
0219   // If we have set the flag to save 'special' events, save here
0220   if (fPhotonCount_Scint + fPhotonCount_Ceren < fDetector->GetSaveThreshold()) {
0221     G4RunManager::GetRunManager()->rndmSaveThisEvent();
0222   }
0223 }
0224
0225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......