GFlash/src/GflashHadronShowerModel.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002
0003 #include "SimG4Core/GFlash/interface/GflashHadronShowerModel.h"
0004
0005 #include "SimGeneral/GFlash/interface/GflashAntiProtonShowerProfile.h"
0006 #include "SimGeneral/GFlash/interface/GflashHadronShowerProfile.h"
0007 #include "SimGeneral/GFlash/interface/GflashHistogram.h"
0008 #include "SimGeneral/GFlash/interface/GflashHit.h"
0009 #include "SimGeneral/GFlash/interface/GflashKaonMinusShowerProfile.h"
0010 #include "SimGeneral/GFlash/interface/GflashKaonPlusShowerProfile.h"
0011 #include "SimGeneral/GFlash/interface/GflashNameSpace.h"
0012 #include "SimGeneral/GFlash/interface/GflashPiKShowerProfile.h"
0013 #include "SimGeneral/GFlash/interface/GflashProtonShowerProfile.h"
0014
0015 #include "G4EventManager.hh"
0016 #include "G4FastSimulationManager.hh"
0017 #include "G4TouchableHandle.hh"
0018 #include "G4TransportationManager.hh"
0019 #include "G4VPhysicalVolume.hh"
0020 #include "G4VSensitiveDetector.hh"
0021
0022 #include "G4AntiProton.hh"
0023 #include "G4KaonMinus.hh"
0024 #include "G4KaonPlus.hh"
0025 #include "G4PionMinus.hh"
0026 #include "G4PionPlus.hh"
0027 #include "G4Proton.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 = new GflashHadronShowerProfile(parSet);
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   theHisto = GflashHistogram::instance();
0047
0048   theGflashStep = new G4Step();
0049   theGflashTouchableHandle = new G4TouchableHistory();
0050   theGflashNavigator = new G4Navigator();
0051 }
0052
0053 GflashHadronShowerModel::~GflashHadronShowerModel() {
0054   delete theProfile;
0055   delete theGflashStep;
0056 }
0057
0058 G4bool GflashHadronShowerModel::IsApplicable(const G4ParticleDefinition &particleType) {
0059   return &particleType == G4PionMinus::PionMinusDefinition() || &particleType == G4PionPlus::PionPlusDefinition() ||
0060          &particleType == G4KaonMinus::KaonMinusDefinition() || &particleType == G4KaonPlus::KaonPlusDefinition() ||
0061          &particleType == G4AntiProton::AntiProtonDefinition() || &particleType == G4Proton::ProtonDefinition();
0062 }
0063
0064 G4bool GflashHadronShowerModel::ModelTrigger(const G4FastTrack &fastTrack) {
0065   // ModelTrigger returns false for Gflash Hadronic Shower Model if it is not
0066   // tested from the corresponding wrapper process, GflashHadronWrapperProcess.
0067   // Temporarily the track status is set to fPostponeToNextEvent at the wrapper
0068   // process before ModelTrigger is really tested for the second time through
0069   // PostStepGPIL of enviaged parameterization processes.  The better
0070   // implmentation may be using via G4VUserTrackInformation of each track, which
0071   // requires to modify a geant source code of stepping (G4SteppingManager2)
0072
0073   G4bool trigger = false;
0074
0075   // mininum energy cutoff to parameterize
0076   if (fastTrack.GetPrimaryTrack()->GetKineticEnergy() < Gflash::energyCutOff * CLHEP::GeV)
0077     return trigger;
0078
0079   // check whether this is called from the normal GPIL or the wrapper process
0080   // GPIL
0081   if (fastTrack.GetPrimaryTrack()->GetTrackStatus() == fPostponeToNextEvent) {
0082     // Shower pameterization start at the first inelastic interaction point
0083     G4bool isInelastic = isFirstInelasticInteraction(fastTrack);
0084
0085     // Other conditions
0086     if (isInelastic) {
0087       trigger = (!excludeDetectorRegion(fastTrack));
0088     }
0089   }
0090
0091   return trigger;
0092 }
0093
0094 void GflashHadronShowerModel::DoIt(const G4FastTrack &fastTrack, G4FastStep &fastStep) {
0095   // kill the particle
0096   fastStep.KillPrimaryTrack();
0097   fastStep.ProposePrimaryTrackPathLength(0.0);
0098
0099   // parameterize energy depostion by the particle type
0100   G4ParticleDefinition *particleType = fastTrack.GetPrimaryTrack()->GetDefinition();
0101
0102   theProfile = thePiKProfile;
0103   if (particleType == G4KaonMinus::KaonMinusDefinition())
0104     theProfile = theKaonMinusProfile;
0105   else if (particleType == G4KaonPlus::KaonPlusDefinition())
0106     theProfile = theKaonPlusProfile;
0107   else if (particleType == G4AntiProton::AntiProtonDefinition())
0108     theProfile = theAntiProtonProfile;
0109   else if (particleType == G4Proton::ProtonDefinition())
0110     theProfile = theProtonProfile;
0111
0112   // input variables for GflashHadronShowerProfile
0113   G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy() / GeV;
0114   G4double globalTime = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetGlobalTime();
0115   G4double charge = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint()->GetCharge();
0116   G4ThreeVector position = fastTrack.GetPrimaryTrack()->GetPosition() / cm;
0117   G4ThreeVector momentum = fastTrack.GetPrimaryTrack()->GetMomentum() / GeV;
0118   G4int showerType = Gflash::findShowerType(position);
0119
0120   theProfile->initialize(showerType, energy, globalTime, charge, position, momentum);
0121   theProfile->loadParameters();
0122   theProfile->hadronicParameterization();
0123
0124   // make hits
0125   makeHits(fastTrack);
0126 }
0127
0128 void GflashHadronShowerModel::makeHits(const G4FastTrack &fastTrack) {
0129   std::vector<GflashHit> &gflashHitList = theProfile->getGflashHitList();
0130
0131   theGflashStep->SetTrack(const_cast<G4Track *>(fastTrack.GetPrimaryTrack()));
0132   theGflashStep->GetPostStepPoint()->SetProcessDefinedStep(
0133       const_cast<G4VProcess *>(fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()));
0134   theGflashNavigator->SetWorldVolume(
0135       G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
0136
0137   for (auto &hit : gflashHitList) {
0138     theGflashNavigator->LocateGlobalPointAndUpdateTouchableHandle(
0139         hit.getPosition(), G4ThreeVector(0, 0, 0), theGflashTouchableHandle, false);
0140     updateGflashStep(hit.getPosition(), hit.getTime());
0141
0142     const G4VPhysicalVolume *aCurrentVolume = theGflashStep->GetPreStepPoint()->GetPhysicalVolume();
0143     if (aCurrentVolume == nullptr)
0144       continue;
0145
0146     const G4LogicalVolume *lv = aCurrentVolume->GetLogicalVolume();
0147     if (lv->GetRegion()->GetName() != "CaloRegion")
0148       continue;
0149
0150     theGflashStep->GetPreStepPoint()->SetSensitiveDetector(aCurrentVolume->GetLogicalVolume()->GetSensitiveDetector());
0151     G4VSensitiveDetector *aSensitive = theGflashStep->GetPreStepPoint()->GetSensitiveDetector();
0152
0153     if (aSensitive == nullptr)
0154       continue;
0155
0156     G4String nameCalor = aCurrentVolume->GetName();
0157     nameCalor.assign(nameCalor, 0, 2);
0158     G4double samplingWeight = 1.0;
0159     if (nameCalor == "HB") {
0160       samplingWeight = Gflash::scaleSensitiveHB;
0161     } else if (nameCalor == "HE" || nameCalor == "HT") {
0162       samplingWeight = Gflash::scaleSensitiveHE;
0163     }
0164     theGflashStep->SetTotalEnergyDeposit(hit.getEnergy() * samplingWeight);
0165
0166     aSensitive->Hit(theGflashStep);
0167   }
0168 }
0169
0170 void GflashHadronShowerModel::updateGflashStep(const G4ThreeVector &spotPosition, G4double timeGlobal) {
0171   theGflashStep->GetPostStepPoint()->SetGlobalTime(timeGlobal);
0172   theGflashStep->GetPreStepPoint()->SetPosition(spotPosition);
0173   theGflashStep->GetPostStepPoint()->SetPosition(spotPosition);
0174   theGflashStep->GetPreStepPoint()->SetTouchableHandle(theGflashTouchableHandle);
0175 }
0176
0177 G4bool GflashHadronShowerModel::isFirstInelasticInteraction(const G4FastTrack &fastTrack) {
0178   G4bool isFirst = false;
0179
0180   G4StepPoint *preStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint();
0181   G4StepPoint *postStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint();
0182
0183   G4String procName = postStep->GetProcessDefinedStep()->GetProcessName();
0184   G4ParticleDefinition *particleType = fastTrack.GetPrimaryTrack()->GetDefinition();
0185
0186   //@@@ this part is still temporary and the cut for the variable ratio should
0187   // be optimized later
0188
0189   if ((particleType == G4PionPlus::PionPlusDefinition() && procName == "WrappedPionPlusInelastic") ||
0190       (particleType == G4PionMinus::PionMinusDefinition() && procName == "WrappedPionMinusInelastic") ||
0191       (particleType == G4KaonPlus::KaonPlusDefinition() && procName == "WrappedKaonPlusInelastic") ||
0192       (particleType == G4KaonMinus::KaonMinusDefinition() && procName == "WrappedKaonMinusInelastic") ||
0193       (particleType == G4AntiProton::AntiProtonDefinition() && procName == "WrappedAntiProtonInelastic") ||
0194       (particleType == G4Proton::ProtonDefinition() && procName == "WrappedProtonInelastic")) {
0195     // skip to the second interaction if the first inelastic is a quasi-elastic
0196     // like interaction
0197     //@@@ the cut may be optimized later
0198
0199     const G4TrackVector *fSecondaryVector = fastTrack.GetPrimaryTrack()->GetStep()->GetSecondary();
0200     G4double leadingEnergy = 0.0;
0201
0202     // loop over 'all' secondaries including those produced by continuous
0203     // processes.
0204     //@@@may require an additional condition only for hadron interaction with
0205     // the process name, but it will not change the result anyway
0206
0207     for (unsigned int isec = 0; isec < fSecondaryVector->size(); isec++) {
0208       G4Track *fSecondaryTrack = (*fSecondaryVector)[isec];
0209       G4double secondaryEnergy = fSecondaryTrack->GetKineticEnergy();
0210
0211       if (secondaryEnergy > leadingEnergy) {
0212         leadingEnergy = secondaryEnergy;
0213       }
0214     }
0215
0216     if ((preStep->GetTotalEnergy() != 0) && (leadingEnergy / preStep->GetTotalEnergy() < Gflash::QuasiElasticLike))
0217       isFirst = true;
0218
0219     // Fill debugging histograms and check information on secondaries -
0220     // remove after final implimentation
0221
0222     if (theHisto->getStoreFlag()) {
0223       theHisto->preStepPosition->Fill(preStep->GetPosition().getRho() / cm);
0224       theHisto->postStepPosition->Fill(postStep->GetPosition().getRho() / cm);
0225       theHisto->deltaStep->Fill((postStep->GetPosition() - preStep->GetPosition()).getRho() / cm);
0226       theHisto->kineticEnergy->Fill(fastTrack.GetPrimaryTrack()->GetKineticEnergy() / GeV);
0227       theHisto->energyLoss->Fill(fabs(fastTrack.GetPrimaryTrack()->GetStep()->GetDeltaEnergy() / GeV));
0228       theHisto->energyRatio->Fill(leadingEnergy / preStep->GetTotalEnergy());
0229     }
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("SimGeneralGFlash") << "GflashHadronShowerModel: excluding region of eta = " << eta;
0245     }
0246     return true;
0247   } else {
0248     const 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
0254     // preshower
0255     if (std::fabs(eta) > Gflash::EtaMin[Gflash::kENCA] &&
0256         std::fabs((postStep->GetPosition()).getZ() / CLHEP::cm) < 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
0262     // end of the hadronic envelopes (may need to be optimized further!)
0263     //@@@if we extend parameterization including Magnet/HO, we need to change
0264     // this strategy
0265     if (kCalor == Gflash::kHB) {
0266       distOut = Gflash::Rmax[Gflash::kHB] - postStep->GetPosition().getRho() * invcm;
0267       if (distOut < Gflash::MinDistanceToOut)
0268         isExcluded = true;
0269     } else if (kCalor == Gflash::kHE) {
0270       distOut = Gflash::Zmax[Gflash::kHE] - std::fabs(postStep->GetPosition().getZ() * invcm);
0271       if (distOut < Gflash::MinDistanceToOut)
0272         isExcluded = true;
0273     }
0274
0275     //@@@remove this print statement later
0276     if (isExcluded && verbosity > 0) {
0277       G4cout << "GflashHadronShowerModel: skipping kCalor = " << kCalor << " DistanceToOut " << distOut << " from ("
0278              << (postStep->GetPosition()).getRho() * invcm << ":" << (postStep->GetPosition()).getZ() * invcm
0279              << ") of KE = " << fastTrack.GetPrimaryTrack()->GetKineticEnergy() / CLHEP::GeV << G4endl;
0280     }
0281   }
0282
0283   return isExcluded;
0284 }