Par04/src/Par04EventAction.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 Par04EventAction.cc
0027 /// \brief Implementation of the Par04EventAction class
0028
0029 #include "Par04EventAction.hh"
0030
0031 #include "Par04DetectorConstruction.hh"  // for Par04DetectorConstruction
0032 #include "Par04Hit.hh"  // for Par04Hit, Par04HitsCollection
0033 #include "Par04ParallelFullWorld.hh"
0034
0035 #include "G4AnalysisManager.hh"  // for G4AnalysisManager
0036 #include "G4Event.hh"  // for G4Event
0037 #include "G4EventManager.hh"  // for G4EventManager
0038 #include "G4Exception.hh"  // for G4Exception, G4ExceptionDesc...
0039 #include "G4ExceptionSeverity.hh"  // for FatalException
0040 #include "G4GenericAnalysisManager.hh"  // for G4GenericAnalysisManager
0041 #include "G4HCofThisEvent.hh"  // for G4HCofThisEvent
0042 #include "G4PrimaryParticle.hh"  // for G4PrimaryParticle
0043 #include "G4PrimaryVertex.hh"  // for G4PrimaryVertex
0044 #include "G4SDManager.hh"  // for G4SDManager
0045 #include "G4SystemOfUnits.hh"  // for GeV
0046 #include "G4THitsCollection.hh"  // for G4THitsCollection
0047 #include "G4ThreeVector.hh"  // for G4ThreeVector
0048 #include "G4Timer.hh"  // for G4Timer
0049 #include "G4UserEventAction.hh"  // for G4UserEventAction
0050
0051 #include <CLHEP/Units/SystemOfUnits.h>  // for GeV
0052 #include <CLHEP/Vector/ThreeVector.h>  // for Hep3Vector
0053 #include <algorithm>  // for max
0054 #include <cmath>  // for log10
0055 #include <cstddef>  // for size_t
0056 #include <ostream>  // for basic_ostream::operator<<
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059
0060 Par04EventAction::Par04EventAction(Par04DetectorConstruction* aDetector,
0061                                    Par04ParallelFullWorld* aParallel)
0062   : G4UserEventAction(),
0063     fHitCollectionID(-1),
0064     fPhysicalFullHitCollectionID(-1),
0065     fPhysicalFastHitCollectionID(-1),
0066     fTimer(),
0067     fDetector(aDetector),
0068     fParallel(aParallel)
0069 {
0070   fCellNbRho = aDetector->GetMeshNbOfCells().x();
0071   fCellNbPhi = aDetector->GetMeshNbOfCells().y();
0072   fCellNbZ = aDetector->GetMeshNbOfCells().z();
0073   fCalEdep.reserve(fCellNbRho * fCellNbPhi * fCellNbZ);
0074   fCalRho.reserve(fCellNbRho * fCellNbPhi * fCellNbZ);
0075   fCalPhi.reserve(fCellNbRho * fCellNbPhi * fCellNbZ);
0076   fCalZ.reserve(fCellNbRho * fCellNbPhi * fCellNbZ);
0077   fCalPhysicalEdep.reserve(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices);
0078   fCalPhysicalLayer.reserve(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices);
0079   fCalPhysicalSlice.reserve(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices);
0080   fCalPhysicalRow.reserve(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices);
0081 }
0082
0083 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0084
0085 Par04EventAction::~Par04EventAction() = default;
0086
0087 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0088
0089 void Par04EventAction::BeginOfEventAction(const G4Event*)
0090 {
0091   StartTimer();
0092 }
0093
0094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0095
0096 void Par04EventAction::StartTimer()
0097 {
0098   fTimer.Start();
0099 }
0100
0101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0102
0103 void Par04EventAction::StopTimer()
0104 {
0105   fTimer.Stop();
0106 }
0107
0108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0109
0110 void Par04EventAction::EndOfEventAction(const G4Event* aEvent)
0111 {
0112   G4SDManager::GetSDMpointer()->GetHCtable();
0113   StopTimer();
0114
0115   // Get hits collection ID (only once)
0116   if (fHitCollectionID == -1) {
0117     fHitCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("hits");
0118   }
0119   if (fPhysicalFullHitCollectionID == -1) {
0120     fPhysicalFullHitCollectionID =
0121       G4SDManager::GetSDMpointer()->GetCollectionID("physicalCellsFullSim");
0122   }
0123   if (fPhysicalFastHitCollectionID == -1) {
0124     fPhysicalFastHitCollectionID =
0125       G4SDManager::GetSDMpointer()->GetCollectionID("physicalCellsFastSim");
0126   }
0127   // Get hits collection
0128   auto hitsCollection =
0129     static_cast<Par04HitsCollection*>(aEvent->GetHCofThisEvent()->GetHC(fHitCollectionID));
0130   auto physicalFullHitsCollection = static_cast<Par04HitsCollection*>(
0131     aEvent->GetHCofThisEvent()->GetHC(fPhysicalFullHitCollectionID));
0132   auto physicalFastHitsCollection = static_cast<Par04HitsCollection*>(
0133     aEvent->GetHCofThisEvent()->GetHC(fPhysicalFastHitCollectionID));
0134
0135   if (hitsCollection == nullptr) {
0136     G4ExceptionDescription msg;
0137     msg << "Cannot access hitsCollection ID " << fHitCollectionID;
0138     G4Exception("Par04EventAction::GetHitsCollection()", "MyCode0001", FatalException, msg);
0139   }
0140   if (physicalFullHitsCollection == nullptr) {
0141     G4ExceptionDescription msg;
0142     msg << "Cannot access physical full sim hitsCollection ID " << fPhysicalFullHitCollectionID;
0143     G4Exception("Par04EventAction::GetHitsCollection()", "MyCode0001", FatalException, msg);
0144   }
0145   if (physicalFastHitsCollection == nullptr) {
0146     G4ExceptionDescription msg;
0147     msg << "Cannot access physical fast sim hitsCollection ID " << fPhysicalFastHitCollectionID;
0148     G4Exception("Par04EventAction::GetHitsCollection()", "MyCode0001", FatalException, msg);
0149   }
0150   // Get analysis manager
0151   auto analysisManager = G4AnalysisManager::Instance();
0152   // Retrieve only once detector dimensions
0153   if (fCellSizeZ == 0) {
0154     fCellSizeZ = fDetector->GetMeshSizeOfCells().z();
0155     fCellSizePhi = fDetector->GetMeshSizeOfCells().y();
0156     fCellSizeRho = fDetector->GetMeshSizeOfCells().x();
0157     fCellNbRho = fDetector->GetMeshNbOfCells().x();
0158     fCellNbPhi = fDetector->GetMeshNbOfCells().y();
0159     fCellNbZ = fDetector->GetMeshNbOfCells().z();
0160   }
0161   if (fPhysicalNbLayers == 0) {
0162     fPhysicalNbLayers = fParallel->GetNbOfLayers();
0163     fPhysicalNbSlices = fParallel->GetNbOfSlices();
0164     fPhysicalNbRows = fParallel->GetNbOfRows();
0165   }
0166
0167   // Retrieve information from primary vertex and primary particle
0168   // To calculate shower axis and entry point to the detector
0169   auto primaryVertex =
0170     G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetPrimaryVertex();
0171   auto primaryParticle = primaryVertex->GetPrimary(0);
0172   G4double primaryEnergy = primaryParticle->GetTotalEnergy();
0173   // Estimate from vertex and particle direction the entry point to the detector
0174   // Calculate entrance point to the detector located at z = 0
0175   auto primaryDirection = primaryParticle->GetMomentumDirection();
0176   auto primaryEntrance =
0177     primaryVertex->GetPosition() - primaryVertex->GetPosition().z() * primaryDirection;
0178
0179   // Resize back to initial mesh size
0180   fCalEdep.resize(fCellNbRho * fCellNbPhi * fCellNbZ, 0);
0181   fCalRho.resize(fCellNbRho * fCellNbPhi * fCellNbZ, 0);
0182   fCalPhi.resize(fCellNbRho * fCellNbPhi * fCellNbZ, 0);
0183   fCalZ.resize(fCellNbRho * fCellNbPhi * fCellNbZ, 0);
0184   fCalPhysicalEdep.resize(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices, 0);
0185   fCalPhysicalLayer.resize(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices, 0);
0186   fCalPhysicalSlice.resize(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices, 0);
0187   fCalPhysicalRow.resize(fPhysicalNbLayers * fPhysicalNbRows * fPhysicalNbSlices, 0);
0188
0189   // Fill histograms
0190   Par04Hit* hit = nullptr;
0191   G4double hitEn = 0;
0192   G4double totalEnergy = 0;
0193   G4int hitNum = 0;
0194   G4int totalNum = 0;
0195   G4int hitZ = -1;
0196   G4int hitRho = -1;
0197   G4int hitPhi = -1;
0198   G4int hitType = -1;
0199   G4int numNonZeroThresholdCells = 0;
0200   G4double tDistance = 0., rDistance = 0., phiDistance = 0.;
0201   G4double tFirstMoment = 0., tSecondMoment = 0.;
0202   G4double rFirstMoment = 0., rSecondMoment = 0.;
0203   G4double phiMean = 0.;
0204   for (size_t iHit = 0; iHit < hitsCollection->entries(); iHit++) {
0205     hit = static_cast<Par04Hit*>(hitsCollection->GetHit(iHit));
0206     hitZ = hit->GetZid();
0207     hitRho = hit->GetRhoId();
0208     hitPhi = hit->GetPhiId();
0209     hitEn = hit->GetEdep();
0210     hitNum = hit->GetNdep();
0211     hitType = hit->GetType();
0212     if (hitEn > 0) {
0213       totalEnergy += hitEn;
0214       totalNum += hitNum;
0215       tDistance = hitZ * fCellSizeZ;
0216       rDistance = hitRho * fCellSizeRho;
0217       phiDistance = hitPhi * fCellSizePhi;
0218       tFirstMoment += hitEn * tDistance;
0219       rFirstMoment += hitEn * rDistance;
0220       phiMean += hitEn * phiDistance;
0221       analysisManager->FillH1(4, tDistance, hitEn);
0222       analysisManager->FillH1(5, rDistance, hitEn);
0223       analysisManager->FillH1(10, hitType);
0224       if (hitEn > 0.0005) {  // e > 0.5 keV
0225         fCalEdep[numNonZeroThresholdCells] = hitEn;
0226         fCalRho[numNonZeroThresholdCells] = hitRho;
0227         fCalPhi[numNonZeroThresholdCells] = hitPhi;
0228         fCalZ[numNonZeroThresholdCells] = hitZ;
0229         numNonZeroThresholdCells++;
0230         analysisManager->FillH1(13, std::log10(hitEn));
0231         analysisManager->FillH1(15, hitNum);
0232       }
0233     }
0234   }
0235   tFirstMoment /= totalEnergy;
0236   rFirstMoment /= totalEnergy;
0237   phiMean /= totalEnergy;
0238   analysisManager->FillH1(0, primaryEnergy / GeV);
0239   analysisManager->FillH1(1, totalEnergy / GeV);
0240   analysisManager->FillH1(2, totalEnergy / primaryEnergy);
0241   analysisManager->FillH1(3, fTimer.GetRealElapsed());
0242   analysisManager->FillH1(6, tFirstMoment);
0243   analysisManager->FillH1(7, rFirstMoment);
0244   analysisManager->FillH1(12, numNonZeroThresholdCells);
0245   analysisManager->FillH1(14, totalNum);
0246   // Resize to store only energy hits above threshold
0247   fCalEdep.resize(numNonZeroThresholdCells);
0248   fCalRho.resize(numNonZeroThresholdCells);
0249   fCalPhi.resize(numNonZeroThresholdCells);
0250   fCalZ.resize(numNonZeroThresholdCells);
0251   analysisManager->FillNtupleDColumn(0, 0, primaryEnergy);
0252   analysisManager->FillNtupleDColumn(0, 1, fTimer.GetRealElapsed());
0253   // Second loop over hits to calculate second moments
0254   for (size_t iHit = 0; iHit < hitsCollection->entries(); iHit++) {
0255     hit = static_cast<Par04Hit*>(hitsCollection->GetHit(iHit));
0256     hitEn = hit->GetEdep();
0257     hitZ = hit->GetZid();
0258     hitRho = hit->GetRhoId();
0259     hitPhi = hit->GetPhiId();
0260     if (hitEn > 0) {
0261       tDistance = hitZ * fCellSizeZ;
0262       rDistance = hitRho * fCellSizeRho;
0263       phiDistance = hitPhi * fCellSizePhi;
0264       tSecondMoment += hitEn * std::pow(tDistance - tFirstMoment, 2);
0265       rSecondMoment += hitEn * std::pow(rDistance - rFirstMoment, 2);
0266       analysisManager->FillH1(11, phiDistance - phiMean, hitEn);
0267     }
0268   }
0269   tSecondMoment /= totalEnergy;
0270   rSecondMoment /= totalEnergy;
0271   analysisManager->FillH1(8, tSecondMoment);
0272   analysisManager->FillH1(9, rSecondMoment);
0273
0274   // Fill ntuple with physical readout data
0275   G4double totalPhysicalEnergy = 0;
0276   totalNum = 0;
0277   hitEn = 0;
0278   hitNum = 0;
0279   G4int hitLayer = -1;
0280   G4int hitRow = -1;
0281   G4int hitSlice = -1;
0282   numNonZeroThresholdCells = 0;
0283   for (size_t iHit = 0; iHit < physicalFullHitsCollection->entries(); iHit++) {
0284     hit = static_cast<Par04Hit*>(physicalFullHitsCollection->GetHit(iHit));
0285     hitLayer = hit->GetRhoId();
0286     hitRow = hit->GetZid();
0287     hitSlice = hit->GetPhiId();
0288     hitEn = hit->GetEdep();
0289     hitNum = hit->GetNdep();
0290     if (hitEn > 0) {
0291       totalPhysicalEnergy += hitEn;
0292       totalNum += hitNum;
0293       if (hitEn > 0.0005) {  // e > 0.5 keV
0294         fCalPhysicalEdep[numNonZeroThresholdCells] = hitEn;
0295         fCalPhysicalLayer[numNonZeroThresholdCells] = hitLayer;
0296         fCalPhysicalRow[numNonZeroThresholdCells] = hitRow;
0297         fCalPhysicalSlice[numNonZeroThresholdCells] = hitSlice;
0298         numNonZeroThresholdCells++;
0299         analysisManager->FillH1(19, std::log10(hitEn));
0300         analysisManager->FillH1(21, hitNum);
0301       }
0302     }
0303   }
0304   for (size_t iHit = 0; iHit < physicalFastHitsCollection->entries(); iHit++) {
0305     hit = static_cast<Par04Hit*>(physicalFastHitsCollection->GetHit(iHit));
0306     hitLayer = hit->GetRhoId();
0307     hitRow = hit->GetZid();
0308     hitSlice = hit->GetPhiId();
0309     hitEn = hit->GetEdep();
0310     hitNum = hit->GetNdep();
0311     if (hitEn > 0) {
0312       totalPhysicalEnergy += hitEn;
0313       totalNum += hitNum;
0314       if (hitEn > 0.0005) {  // e > 0.5 keV
0315         fCalPhysicalEdep[numNonZeroThresholdCells] = hitEn;
0316         fCalPhysicalLayer[numNonZeroThresholdCells] = hitLayer;
0317         fCalPhysicalRow[numNonZeroThresholdCells] = hitRow;
0318         fCalPhysicalSlice[numNonZeroThresholdCells] = hitSlice;
0319         numNonZeroThresholdCells++;
0320         analysisManager->FillH1(19, std::log10(hitEn));
0321         analysisManager->FillH1(21, hitNum);
0322       }
0323     }
0324   }
0325   analysisManager->FillH1(16, totalPhysicalEnergy / GeV);
0326   analysisManager->FillH1(17, totalPhysicalEnergy / primaryEnergy);
0327   analysisManager->FillH1(18, numNonZeroThresholdCells);
0328   analysisManager->FillH1(20, totalNum);
0329   fCalPhysicalEdep.resize(numNonZeroThresholdCells);
0330   fCalPhysicalLayer.resize(numNonZeroThresholdCells);
0331   fCalPhysicalSlice.resize(numNonZeroThresholdCells);
0332   fCalPhysicalRow.resize(numNonZeroThresholdCells);
0333   analysisManager->AddNtupleRow(0);
0334   analysisManager->AddNtupleRow(1);
0335   analysisManager->AddNtupleRow(2);
0336 }