Application/src/GFlashHadronShowerModel.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002
0003 #include "SimG4Core/Application/interface/GFlashHadronShowerModel.h"
0004 #include "SimG4Core/Application/interface/SteppingAction.h"
0005
0006 #include "SimGeneral/GFlash/interface/GflashHadronShowerProfile.h"
0007 #include "SimGeneral/GFlash/interface/GflashPiKShowerProfile.h"
0008 #include "SimGeneral/GFlash/interface/GflashKaonPlusShowerProfile.h"
0009 #include "SimGeneral/GFlash/interface/GflashKaonMinusShowerProfile.h"
0010 #include "SimGeneral/GFlash/interface/GflashProtonShowerProfile.h"
0011 #include "SimGeneral/GFlash/interface/GflashAntiProtonShowerProfile.h"
0012 #include "SimGeneral/GFlash/interface/GflashNameSpace.h"
0013 #include "SimGeneral/GFlash/interface/GflashHit.h"
0014
0015 #include "G4FastSimulationManager.hh"
0016 #include "G4TransportationManager.hh"
0017 #include "G4TouchableHandle.hh"
0018 #include "G4VSensitiveDetector.hh"
0019 #include "G4VPhysicalVolume.hh"
0020 #include "G4EventManager.hh"
0021
0022 #include "G4PionMinus.hh"
0023 #include "G4PionPlus.hh"
0024 #include "G4KaonMinus.hh"
0025 #include "G4KaonPlus.hh"
0026 #include "G4Proton.hh"
0027 #include "G4AntiProton.hh"
0028 #include "G4VProcess.hh"
0029
0030 #include <DD4hep/Filter.h>
0031 #include <vector>
0032
0033 using namespace CLHEP;
0034
0035 GFlashHadronShowerModel::GFlashHadronShowerModel(G4String modelName,
0036                                                  G4Region* envelope,
0037                                                  const edm::ParameterSet& parSet)
0038     : G4VFastSimulationModel(modelName, envelope), theParSet(parSet) {
0039   theWatcherOn = parSet.getParameter<bool>("watcherOn");
0040
0041   theProfile = nullptr;
0042   thePiKProfile = new GflashPiKShowerProfile(parSet);
0043   theKaonPlusProfile = new GflashKaonPlusShowerProfile(parSet);
0044   theKaonMinusProfile = new GflashKaonMinusShowerProfile(parSet);
0045   theProtonProfile = new GflashProtonShowerProfile(parSet);
0046   theAntiProtonProfile = new GflashAntiProtonShowerProfile(parSet);
0047
0048   theRegion = reinterpret_cast<G4Region*>(envelope);
0049
0050   theGflashStep = new G4Step();
0051   theGflashTouchableHandle = new G4TouchableHistory();
0052   theGflashNavigator = new G4Navigator();
0053 }
0054
0055 GFlashHadronShowerModel::~GFlashHadronShowerModel() {
0056   delete thePiKProfile;
0057   delete theKaonPlusProfile;
0058   delete theKaonMinusProfile;
0059   delete theProtonProfile;
0060   delete theAntiProtonProfile;
0061   delete theGflashStep;
0062 }
0063
0064 G4bool GFlashHadronShowerModel::IsApplicable(const G4ParticleDefinition& particleType) {
0065   return &particleType == G4PionMinus::PionMinusDefinition() || &particleType == G4PionPlus::PionPlusDefinition() ||
0066          &particleType == G4KaonMinus::KaonMinusDefinition() || &particleType == G4KaonPlus::KaonPlusDefinition() ||
0067          &particleType == G4AntiProton::AntiProtonDefinition() || &particleType == G4Proton::ProtonDefinition();
0068 }
0069
0070 G4bool GFlashHadronShowerModel::ModelTrigger(const G4FastTrack& fastTrack) {
0071   // ModelTrigger returns false for Gflash Hadronic Shower Model if it is not
0072   // tested from the corresponding wrapper process, GflashHadronWrapperProcess.
0073   // Temporarily the track status is set to fPostponeToNextEvent at the wrapper
0074   // process before ModelTrigger is really tested for the second time through
0075   // PostStepGPIL of enviaged parameterization processes.  The better
0076   // implmentation may be using via G4VUserTrackInformation of each track, which
0077   // requires to modify a geant source code of stepping (G4SteppingManager2)
0078
0079   // mininum energy cutoff to parameterize
0080   if (fastTrack.GetPrimaryTrack()->GetKineticEnergy() < CLHEP::GeV * Gflash::energyCutOff) {
0081     return false;
0082   }
0083
0084   // This will be changed accordingly when the way
0085   // dealing with CaloRegion changes later.
0086   const G4VPhysicalVolume* pCurrentVolume = fastTrack.GetPrimaryTrack()->GetTouchable()->GetVolume();
0087   if (pCurrentVolume == nullptr) {
0088     return false;
0089   }
0090
0091   const G4LogicalVolume* lv = pCurrentVolume->GetLogicalVolume();
0092   if (lv->GetRegion() != theRegion) {
0093     return false;
0094   }
0095   return !(excludeDetectorRegion(fastTrack));
0096 }
0097
0098 void GFlashHadronShowerModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) {
0099   // kill the particle
0100   fastStep.KillPrimaryTrack();
0101   fastStep.ProposePrimaryTrackPathLength(0.0);
0102
0103   // parameterize energy depostion by the particle type
0104   G4ParticleDefinition* particleType = fastTrack.GetPrimaryTrack()->GetDefinition();
0105
0106   theProfile = thePiKProfile;
0107   if (particleType == G4KaonMinus::KaonMinusDefinition())
0108     theProfile = theKaonMinusProfile;
0109   else if (particleType == G4KaonPlus::KaonPlusDefinition())
0110     theProfile = theKaonPlusProfile;
0111   else if (particleType == G4AntiProton::AntiProtonDefinition())
0112     theProfile = theAntiProtonProfile;
0113   else if (particleType == G4Proton::ProtonDefinition())
0114     theProfile = theProtonProfile;
0115
0116   //input variables for GflashHadronShowerProfile
0117   G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy() / GeV;
0118   G4double globalTime = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetGlobalTime();
0119   G4double charge = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint()->GetCharge();
0120   G4ThreeVector position = fastTrack.GetPrimaryTrack()->GetPosition() / cm;
0121   G4ThreeVector momentum = fastTrack.GetPrimaryTrack()->GetMomentum() / GeV;
0122   G4int showerType = Gflash::findShowerType(position);
0123
0124   theProfile->initialize(showerType, energy, globalTime, charge, position, momentum);
0125   theProfile->loadParameters();
0126   theProfile->hadronicParameterization();
0127
0128   // make hits
0129   makeHits(fastTrack);
0130 }
0131
0132 void GFlashHadronShowerModel::makeHits(const G4FastTrack& fastTrack) {
0133   std::vector<GflashHit>& gflashHitList = theProfile->getGflashHitList();
0134
0135   theGflashStep->SetTrack(const_cast<G4Track*>(fastTrack.GetPrimaryTrack()));
0136   theGflashStep->GetPostStepPoint()->SetProcessDefinedStep(
0137       const_cast<G4VProcess*>(fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()));
0138   theGflashNavigator->SetWorldVolume(
0139       G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
0140
0141   std::vector<GflashHit>::const_iterator spotIter = gflashHitList.begin();
0142   std::vector<GflashHit>::const_iterator spotIterEnd = gflashHitList.end();
0143
0144   for (; spotIter != spotIterEnd; spotIter++) {
0145     theGflashNavigator->LocateGlobalPointAndUpdateTouchableHandle(
0146         spotIter->getPosition(), G4ThreeVector(0, 0, 0), theGflashTouchableHandle, false);
0147     updateGflashStep(spotIter->getPosition(), spotIter->getTime());
0148
0149     // if there is a watcher defined in a job and the flag is turned on
0150     if (theWatcherOn) {
0151       theGflashStep->SetTotalEnergyDeposit(spotIter->getEnergy());
0152       SteppingAction* userSteppingAction = (SteppingAction*)G4EventManager::GetEventManager()->GetUserSteppingAction();
0153       userSteppingAction->m_g4StepSignal(theGflashStep);
0154     }
0155
0156     G4VPhysicalVolume* aCurrentVolume = theGflashStep->GetPreStepPoint()->GetPhysicalVolume();
0157     if (aCurrentVolume == nullptr) {
0158       continue;
0159     }
0160
0161     G4LogicalVolume* lv = aCurrentVolume->GetLogicalVolume();
0162     if (lv->GetRegion() != theRegion) {
0163       continue;
0164     }
0165
0166     theGflashStep->GetPreStepPoint()->SetSensitiveDetector(aCurrentVolume->GetLogicalVolume()->GetSensitiveDetector());
0167     G4VSensitiveDetector* aSensitive = theGflashStep->GetPreStepPoint()->GetSensitiveDetector();
0168
0169     if (aSensitive == nullptr)
0170       continue;
0171
0172     G4String nameCalor = (G4String)(dd4hep::dd::noNamespace(lv->GetName()));
0173     nameCalor.assign(nameCalor, 0, 2);
0174     G4double samplingWeight = 1.0;
0175     if (nameCalor == "HB") {
0176       samplingWeight = Gflash::scaleSensitiveHB;
0177     } else if (nameCalor == "HE" || nameCalor == "HT") {
0178       samplingWeight = Gflash::scaleSensitiveHE;
0179     }
0180     theGflashStep->SetTotalEnergyDeposit(spotIter->getEnergy() * samplingWeight);
0181
0182     aSensitive->Hit(theGflashStep);
0183   }
0184 }
0185
0186 void GFlashHadronShowerModel::updateGflashStep(const G4ThreeVector& spotPosition, G4double timeGlobal) {
0187   theGflashStep->GetPostStepPoint()->SetGlobalTime(timeGlobal);
0188   theGflashStep->GetPreStepPoint()->SetPosition(spotPosition);
0189   theGflashStep->GetPostStepPoint()->SetPosition(spotPosition);
0190   theGflashStep->GetPreStepPoint()->SetTouchableHandle(theGflashTouchableHandle);
0191 }
0192
0193 G4bool GFlashHadronShowerModel::isFirstInelasticInteraction(const G4FastTrack& fastTrack) {
0194   G4bool isFirst = false;
0195
0196   G4StepPoint* preStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint();
0197   G4StepPoint* postStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint();
0198
0199   G4String procName = postStep->GetProcessDefinedStep()->GetProcessName();
0200   G4ParticleDefinition* particleType = fastTrack.GetPrimaryTrack()->GetDefinition();
0201
0202   //@@@ this part is still temporary and the cut for the variable ratio should be optimized later
0203
0204   if ((particleType == G4PionPlus::PionPlusDefinition() && procName == "WrappedPionPlusInelastic") ||
0205       (particleType == G4PionMinus::PionMinusDefinition() && procName == "WrappedPionMinusInelastic") ||
0206       (particleType == G4KaonPlus::KaonPlusDefinition() && procName == "WrappedKaonPlusInelastic") ||
0207       (particleType == G4KaonMinus::KaonMinusDefinition() && procName == "WrappedKaonMinusInelastic") ||
0208       (particleType == G4AntiProton::AntiProtonDefinition() && procName == "WrappedAntiProtonInelastic") ||
0209       (particleType == G4Proton::ProtonDefinition() && procName == "WrappedProtonInelastic")) {
0210     //skip to the second interaction if the first inelastic is a quasi-elastic like interaction
0211     //@@@ the cut may be optimized later
0212
0213     const G4TrackVector* fSecondaryVector = fastTrack.GetPrimaryTrack()->GetStep()->GetSecondary();
0214     G4double leadingEnergy = 0.0;
0215
0216     //loop over 'all' secondaries including those produced by continuous processes.
0217     //@@@may require an additional condition only for hadron interaction with the process name,
0218     //but it will not change the result anyway
0219
0220     for (unsigned int isec = 0; isec < fSecondaryVector->size(); isec++) {
0221       G4Track* fSecondaryTrack = (*fSecondaryVector)[isec];
0222       G4double secondaryEnergy = fSecondaryTrack->GetKineticEnergy();
0223
0224       if (secondaryEnergy > leadingEnergy) {
0225         leadingEnergy = secondaryEnergy;
0226       }
0227     }
0228
0229     if ((preStep->GetTotalEnergy() != 0) && (leadingEnergy / preStep->GetTotalEnergy() < Gflash::QuasiElasticLike))
0230       isFirst = true;
0231   }
0232   return isFirst;
0233 }
0234
0235 G4bool GFlashHadronShowerModel::excludeDetectorRegion(const G4FastTrack& fastTrack) {
0236   const double invcm = 1.0 / CLHEP::cm;
0237   G4bool isExcluded = false;
0238   int verbosity = theParSet.getUntrackedParameter<int>("Verbosity");
0239
0240   //exclude regions where geometry are complicated
0241   //+- one supermodule around the EB/EE boundary: 1.479 +- 0.0174*5
0242   G4double eta = fastTrack.GetPrimaryTrack()->GetPosition().pseudoRapidity();
0243   if (std::fabs(eta) > 1.392 && std::fabs(eta) < 1.566) {
0244     if (verbosity > 0) {
0245       edm::LogVerbatim("SimG4CoreApplication") << "GFlashHadronShowerModel: excluding region of eta = " << eta;
0246     }
0247     return true;
0248   } else {
0249     G4StepPoint* postStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint();
0250
0251     Gflash::CalorimeterNumber kCalor = Gflash::getCalorimeterNumber(postStep->GetPosition() * invcm);
0252     G4double distOut = 9999.0;
0253
0254     //exclude the region where the shower starting point is inside the preshower
0255     if (std::fabs(eta) > Gflash::EtaMin[Gflash::kENCA] &&
0256         std::fabs((postStep->GetPosition()).getZ() * invcm) < Gflash::Zmin[Gflash::kENCA]) {
0257       return true;
0258     }
0259
0260     //<---the shower starting point is always inside envelopes
0261     //@@@exclude the region where the shower starting point is too close to the end of
0262     //the hadronic envelopes (may need to be optimized further!)
0263     //@@@if we extend parameterization including Magnet/HO, we need to change this strategy
0264     if (kCalor == Gflash::kHB) {
0265       distOut = Gflash::Rmax[Gflash::kHB] - postStep->GetPosition().getRho() * invcm;
0266       if (distOut < Gflash::MinDistanceToOut)
0267         isExcluded = true;
0268     } else if (kCalor == Gflash::kHE) {
0269       distOut = Gflash::Zmax[Gflash::kHE] - std::fabs(postStep->GetPosition().getZ() * invcm);
0270       if (distOut < Gflash::MinDistanceToOut)
0271         isExcluded = true;
0272     }
0273
0274     //@@@remove this print statement later
0275     if (isExcluded && verbosity > 0) {
0276       edm::LogVerbatim("SimG4CoreApplication")
0277           << "GFlashHadronShowerModel: skipping kCalor = " << kCalor << " DistanceToOut " << distOut << " from ("
0278           << (postStep->GetPosition()).getRho() * invcm << ":" << (postStep->GetPosition()).getZ() * invcm
0279           << ") of KE = " << fastTrack.GetPrimaryTrack()->GetKineticEnergy() / CLHEP::GeV;
0280     }
0281   }
0282
0283   return isExcluded;
0284 }