Par03/src/Par03EventAction.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 #include "Par03EventAction.hh"
0027
0028 #include "Par03DetectorConstruction.hh"
0029 #include "Par03Hit.hh"
0030
0031 #include "G4AnalysisManager.hh"
0032 #include "G4Event.hh"
0033 #include "G4EventManager.hh"
0034 #include "G4HCofThisEvent.hh"
0035 #include "G4SDManager.hh"
0036
0037 Par03EventAction::Par03EventAction(Par03DetectorConstruction* aDetector)
0038   : G4UserEventAction(), fHitCollectionID(-1), fTimer(), fDetector(aDetector)
0039 {}
0040
0041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0042
0043 Par03EventAction::~Par03EventAction() = default;
0044
0045 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0046
0047 void Par03EventAction::BeginOfEventAction(const G4Event*)
0048 {
0049   fTimer.Start();
0050 }
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053
0054 void Par03EventAction::EndOfEventAction(const G4Event* aEvent)
0055 {
0056   fTimer.Stop();
0057   // Get hits collection ID (only once)
0058   if (fHitCollectionID == -1) {
0059     fHitCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("hits");
0060   }
0061   // Get hits collection
0062   auto hitsCollection =
0063     static_cast<Par03HitsCollection*>(aEvent->GetHCofThisEvent()->GetHC(fHitCollectionID));
0064
0065   if (hitsCollection == nullptr) {
0066     G4ExceptionDescription msg;
0067     msg << "Cannot access hitsCollection ID " << fHitCollectionID;
0068     G4Exception("Par03EventAction::GetHitsCollection()", "MyCode0001", FatalException, msg);
0069   }
0070   // Get analysis manager
0071   auto analysisManager = G4AnalysisManager::Instance();
0072   // Retrieve only once detector dimensions
0073   if (fCellSizeZ == 0) {
0074     fCellSizeZ = fDetector->GetLength() / fDetector->GetNbOfLayers();
0075     fCellSizeRho = fDetector->GetRadius() / fDetector->GetNbOfRhoCells();
0076   }
0077
0078   // Retrieve information from primary vertex and primary particle
0079   // To calculate shower axis and entry point to the detector
0080   auto primaryVertex =
0081     G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetPrimaryVertex();
0082   auto primaryParticle = primaryVertex->GetPrimary(0);
0083   G4double primaryEnergy = primaryParticle->GetTotalEnergy();
0084   // Estimate from vertex and particle direction the entry point to the detector
0085   // Calculate entrance point to the detector located at z = 0
0086   auto primaryDirection = primaryParticle->GetMomentumDirection();
0087   auto primaryEntrance =
0088     primaryVertex->GetPosition() - primaryVertex->GetPosition().z() * primaryDirection;
0089   G4double cosDirection = std::cos(primaryDirection.theta());
0090   G4double sinDirection = std::sin(primaryDirection.theta());
0091
0092   // Fill histograms
0093   Par03Hit* hit = nullptr;
0094   G4double hitEn = 0;
0095   G4double totalEnergy = 0;
0096   G4int hitZ = -1;
0097   G4int hitRho = -1;
0098   G4int hitType = -1;
0099   G4double tDistance = 0., rDistance = 0.;
0100   G4double tFirstMoment = 0., tSecondMoment = 0.;
0101   G4double rFirstMoment = 0., rSecondMoment = 0.;
0102   for (size_t iHit = 0; iHit < hitsCollection->entries(); iHit++) {
0103     hit = static_cast<Par03Hit*>(hitsCollection->GetHit(iHit));
0104     hitZ = hit->GetZid();
0105     hitRho = hit->GetRhoId();
0106     hitEn = hit->GetEdep();
0107     hitType = hit->GetType();
0108     if (hitEn > 0) {
0109       totalEnergy += hitEn;
0110       tDistance = hitZ * fCellSizeZ * cosDirection
0111                   + (hitRho * fCellSizeRho - primaryEntrance.perp()) * sinDirection;
0112       rDistance = hitZ * fCellSizeZ * (-sinDirection)
0113                   + (hitRho * fCellSizeRho - primaryEntrance.perp()) * cosDirection;
0114       tFirstMoment += hitEn * tDistance;
0115       rFirstMoment += hitEn * rDistance;
0116       analysisManager->FillH1(4, tDistance, hitEn);
0117       analysisManager->FillH1(5, rDistance, hitEn);
0118       analysisManager->FillH1(10, hitType);
0119     }
0120   }
0121   tFirstMoment /= totalEnergy;
0122   rFirstMoment /= totalEnergy;
0123   analysisManager->FillH1(0, primaryEnergy / GeV);
0124   analysisManager->FillH1(1, totalEnergy / GeV);
0125   analysisManager->FillH1(2, totalEnergy / primaryEnergy);
0126   analysisManager->FillH1(3, fTimer.GetRealElapsed());
0127   analysisManager->FillH1(6, tFirstMoment);
0128   analysisManager->FillH1(7, rFirstMoment);
0129
0130   // Second loop over hits to calculate second moments
0131   for (size_t iHit = 0; iHit < hitsCollection->entries(); iHit++) {
0132     hit = static_cast<Par03Hit*>(hitsCollection->GetHit(iHit));
0133     hitEn = hit->GetEdep();
0134     hitZ = hit->GetZid();
0135     hitRho = hit->GetRhoId();
0136     if (hitEn > 0) {
0137       tDistance = hitZ * fCellSizeZ * cosDirection
0138                   + (hitRho * fCellSizeRho - primaryEntrance.r()) * sinDirection;
0139       rDistance = hitZ * fCellSizeZ * (-sinDirection)
0140                   + (hitRho * fCellSizeRho - primaryEntrance.r()) * cosDirection;
0141       tSecondMoment += hitEn * std::pow(tDistance - tFirstMoment, 2);
0142       rSecondMoment += hitEn * std::pow(rDistance - rFirstMoment, 2);
0143     }
0144   }
0145   tSecondMoment /= totalEnergy;
0146   rSecondMoment /= totalEnergy;
0147   analysisManager->FillH1(8, tSecondMoment);
0148   analysisManager->FillH1(9, rSecondMoment);
0149 }