Application/src/GFlashEMShowerModel.cc

0001 //
0002 // initial setup : E.Barberio & Joanna Weng
0003 // big changes : Soon Jun & Dongwook Jang
0004 // V.Ivanchenko rename the class, cleanup, and move
0005 //              to SimG4Core/Application - 2012/08/14
0006 //
0007 #include "SimG4Core/Application/interface/SteppingAction.h"
0008 #include "SimG4Core/Application/interface/GFlashEMShowerModel.h"
0009
0010 #include "SimGeneral/GFlash/interface/GflashEMShowerProfile.h"
0011 #include "SimGeneral/GFlash/interface/GflashHit.h"
0012
0013 #include "G4Electron.hh"
0014 #include "G4Positron.hh"
0015 #include "G4VProcess.hh"
0016 #include "G4VPhysicalVolume.hh"
0017 #include "G4LogicalVolume.hh"
0018 #include "G4TransportationManager.hh"
0019 #include "G4EventManager.hh"
0020 #include "G4FastSimulationManager.hh"
0021 #include "G4TouchableHandle.hh"
0022 #include "G4VSensitiveDetector.hh"
0023 #include <CLHEP/Units/SystemOfUnits.h>
0024 using CLHEP::cm;
0025 using CLHEP::GeV;
0026
0027 GFlashEMShowerModel::GFlashEMShowerModel(const G4String& modelName,
0028                                          G4Envelope* envelope,
0029                                          const edm::ParameterSet& parSet)
0030     : G4VFastSimulationModel(modelName, envelope), theParSet(parSet) {
0031   theWatcherOn = parSet.getParameter<bool>("watcherOn");
0032
0033   theProfile = new GflashEMShowerProfile(parSet);
0034   theRegion = reinterpret_cast<G4Region*>(envelope);
0035
0036   theGflashStep = new G4Step();
0037   theGflashTouchableHandle = new G4TouchableHistory();
0038   theGflashNavigator = new G4Navigator();
0039 }
0040
0041 // --------------------------------------------------------------------------
0042
0043 GFlashEMShowerModel::~GFlashEMShowerModel() {
0044   delete theProfile;
0045   delete theGflashStep;
0046 }
0047
0048 G4bool GFlashEMShowerModel::IsApplicable(const G4ParticleDefinition& particleType) {
0049   return (&particleType == G4Electron::Electron() || &particleType == G4Positron::Positron());
0050 }
0051
0052 // --------------------------------------------------------------------------
0053 G4bool GFlashEMShowerModel::ModelTrigger(const G4FastTrack& fastTrack) {
0054   // Mininum energy cutoff to parameterize
0055   if (fastTrack.GetPrimaryTrack()->GetKineticEnergy() < Gflash::energyCutOff) {
0056     return false;
0057   }
0058   if (excludeDetectorRegion(fastTrack)) {
0059     return false;
0060   }
0061
0062   // This will be changed accordingly when the way
0063   // dealing with CaloRegion changes later.
0064   const G4TouchableHistory* touch = static_cast<const G4TouchableHistory*>(fastTrack.GetPrimaryTrack()->GetTouchable());
0065   const G4VPhysicalVolume* pCurrentVolume = touch->GetVolume();
0066   if (pCurrentVolume == nullptr) {
0067     return false;
0068   }
0069
0070   const G4LogicalVolume* lv = pCurrentVolume->GetLogicalVolume();
0071   if (lv->GetRegion() != theRegion) {
0072     return false;
0073   }
0074   return true;
0075 }
0076
0077 // ---------------------------------------------------------------------------
0078 void GFlashEMShowerModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) {
0079   // Kill the parameterised particle:
0080   fastStep.KillPrimaryTrack();
0081   fastStep.ProposePrimaryTrackPathLength(0.0);
0082
0083   // Input variables for GFlashEMShowerProfile with showerType = 1,5
0084   // Shower starts inside crystals
0085   G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy() / GeV;
0086   G4double globalTime = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetGlobalTime();
0087   G4double charge = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint()->GetCharge();
0088   G4ThreeVector position = fastTrack.GetPrimaryTrack()->GetPosition() / cm;
0089   G4ThreeVector momentum = fastTrack.GetPrimaryTrack()->GetMomentum() / GeV;
0090   G4int showerType = Gflash::findShowerType(position);
0091
0092   // Do actual parameterization
0093   // The result of parameterization is gflashHitList
0094   theProfile->initialize(showerType, energy, globalTime, charge, position, momentum);
0095   theProfile->parameterization();
0096
0097   // Make hits
0098   makeHits(fastTrack);
0099 }
0100
0101 // ---------------------------------------------------------------------------
0102 void GFlashEMShowerModel::makeHits(const G4FastTrack& fastTrack) {
0103   std::vector<GflashHit>& gflashHitList = theProfile->getGflashHitList();
0104
0105   theGflashStep->SetTrack(const_cast<G4Track*>(fastTrack.GetPrimaryTrack()));
0106
0107   theGflashStep->GetPostStepPoint()->SetProcessDefinedStep(
0108       const_cast<G4VProcess*>(fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()));
0109   theGflashNavigator->SetWorldVolume(
0110       G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
0111
0112   std::vector<GflashHit>::const_iterator spotIter = gflashHitList.begin();
0113   std::vector<GflashHit>::const_iterator spotIterEnd = gflashHitList.end();
0114
0115   for (; spotIter != spotIterEnd; ++spotIter) {
0116     // Put touchable for each hit so that touchable history
0117     //     keeps track of each step.
0118     theGflashNavigator->LocateGlobalPointAndUpdateTouchableHandle(
0119         spotIter->getPosition(), G4ThreeVector(0, 0, 0), theGflashTouchableHandle, false);
0120     updateGflashStep(spotIter->getPosition(), spotIter->getTime());
0121
0122     // If there is a watcher defined in a job and the flag is turned on
0123     if (theWatcherOn) {
0124       SteppingAction* userSteppingAction = (SteppingAction*)G4EventManager::GetEventManager()->GetUserSteppingAction();
0125       userSteppingAction->m_g4StepSignal(theGflashStep);
0126     }
0127
0128     // Send G4Step information to Hit/Digi if the volume is sensitive
0129     // Copied from G4SteppingManager.cc
0130
0131     G4VPhysicalVolume* aCurrentVolume = theGflashStep->GetPreStepPoint()->GetPhysicalVolume();
0132     if (aCurrentVolume == nullptr) {
0133       continue;
0134     }
0135
0136     G4LogicalVolume* lv = aCurrentVolume->GetLogicalVolume();
0137     if (lv->GetRegion() != theRegion) {
0138       continue;
0139     }
0140
0141     theGflashStep->GetPreStepPoint()->SetSensitiveDetector(aCurrentVolume->GetLogicalVolume()->GetSensitiveDetector());
0142     G4VSensitiveDetector* aSensitive = theGflashStep->GetPreStepPoint()->GetSensitiveDetector();
0143
0144     if (aSensitive == nullptr) {
0145       continue;
0146     }
0147
0148     theGflashStep->SetTotalEnergyDeposit(spotIter->getEnergy());
0149     aSensitive->Hit(theGflashStep);
0150   }
0151 }
0152
0153 // ---------------------------------------------------------------------------
0154 void GFlashEMShowerModel::updateGflashStep(const G4ThreeVector& spotPosition, G4double timeGlobal) {
0155   theGflashStep->GetPostStepPoint()->SetGlobalTime(timeGlobal);
0156   theGflashStep->GetPreStepPoint()->SetPosition(spotPosition);
0157   theGflashStep->GetPostStepPoint()->SetPosition(spotPosition);
0158   theGflashStep->GetPreStepPoint()->SetTouchableHandle(theGflashTouchableHandle);
0159 }
0160
0161 // ---------------------------------------------------------------------------
0162 G4bool GFlashEMShowerModel::excludeDetectorRegion(const G4FastTrack& fastTrack) {
0163   //exclude regions where geometry are complicated
0164   //+- one supermodule around the EB/EE boundary: 1.479 +- 0.0174*5
0165   G4double eta = fastTrack.GetPrimaryTrack()->GetPosition().pseudoRapidity();
0166   return (std::fabs(eta) > 1.392 && std::fabs(eta) < 1.566);
0167 }
0168
0169 // ---------------------------------------------------------------------------