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