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