File indexing completed on 2024-04-06 12:30:24
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 if (theProfile)
0055 delete theProfile;
0056 if (theGflashStep)
0057 delete theGflashStep;
0058 }
0059
0060 G4bool GflashHadronShowerModel::IsApplicable(const G4ParticleDefinition &particleType) {
0061 return &particleType == G4PionMinus::PionMinusDefinition() || &particleType == G4PionPlus::PionPlusDefinition() ||
0062 &particleType == G4KaonMinus::KaonMinusDefinition() || &particleType == G4KaonPlus::KaonPlusDefinition() ||
0063 &particleType == G4AntiProton::AntiProtonDefinition() || &particleType == G4Proton::ProtonDefinition();
0064 }
0065
0066 G4bool GflashHadronShowerModel::ModelTrigger(const G4FastTrack &fastTrack) {
0067
0068
0069
0070
0071
0072
0073
0074
0075 G4bool trigger = false;
0076
0077
0078 if (fastTrack.GetPrimaryTrack()->GetKineticEnergy() < Gflash::energyCutOff * CLHEP::GeV)
0079 return trigger;
0080
0081
0082
0083 if (fastTrack.GetPrimaryTrack()->GetTrackStatus() == fPostponeToNextEvent) {
0084
0085 G4bool isInelastic = isFirstInelasticInteraction(fastTrack);
0086
0087
0088 if (isInelastic) {
0089 trigger = (!excludeDetectorRegion(fastTrack));
0090 }
0091 }
0092
0093 return trigger;
0094 }
0095
0096 void GflashHadronShowerModel::DoIt(const G4FastTrack &fastTrack, G4FastStep &fastStep) {
0097
0098 fastStep.KillPrimaryTrack();
0099 fastStep.ProposePrimaryTrackPathLength(0.0);
0100
0101
0102 G4ParticleDefinition *particleType = fastTrack.GetPrimaryTrack()->GetDefinition();
0103
0104 theProfile = thePiKProfile;
0105 if (particleType == G4KaonMinus::KaonMinusDefinition())
0106 theProfile = theKaonMinusProfile;
0107 else if (particleType == G4KaonPlus::KaonPlusDefinition())
0108 theProfile = theKaonPlusProfile;
0109 else if (particleType == G4AntiProton::AntiProtonDefinition())
0110 theProfile = theAntiProtonProfile;
0111 else if (particleType == G4Proton::ProtonDefinition())
0112 theProfile = theProtonProfile;
0113
0114
0115 G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy() / GeV;
0116 G4double globalTime = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetGlobalTime();
0117 G4double charge = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint()->GetCharge();
0118 G4ThreeVector position = fastTrack.GetPrimaryTrack()->GetPosition() / cm;
0119 G4ThreeVector momentum = fastTrack.GetPrimaryTrack()->GetMomentum() / GeV;
0120 G4int showerType = Gflash::findShowerType(position);
0121
0122 theProfile->initialize(showerType, energy, globalTime, charge, position, momentum);
0123 theProfile->loadParameters();
0124 theProfile->hadronicParameterization();
0125
0126
0127 makeHits(fastTrack);
0128 }
0129
0130 void GflashHadronShowerModel::makeHits(const G4FastTrack &fastTrack) {
0131 std::vector<GflashHit> &gflashHitList = theProfile->getGflashHitList();
0132
0133 theGflashStep->SetTrack(const_cast<G4Track *>(fastTrack.GetPrimaryTrack()));
0134 theGflashStep->GetPostStepPoint()->SetProcessDefinedStep(
0135 const_cast<G4VProcess *>(fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()));
0136 theGflashNavigator->SetWorldVolume(
0137 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume());
0138
0139 std::vector<GflashHit>::const_iterator spotIter = gflashHitList.begin();
0140 std::vector<GflashHit>::const_iterator spotIterEnd = gflashHitList.end();
0141
0142 for (; spotIter != spotIterEnd; spotIter++) {
0143 theGflashNavigator->LocateGlobalPointAndUpdateTouchableHandle(
0144 spotIter->getPosition(), G4ThreeVector(0, 0, 0), theGflashTouchableHandle, false);
0145 updateGflashStep(spotIter->getPosition(), spotIter->getTime());
0146
0147 const G4VPhysicalVolume *aCurrentVolume = theGflashStep->GetPreStepPoint()->GetPhysicalVolume();
0148 if (aCurrentVolume == nullptr)
0149 continue;
0150
0151 const G4LogicalVolume *lv = aCurrentVolume->GetLogicalVolume();
0152 if (lv->GetRegion()->GetName() != "CaloRegion")
0153 continue;
0154
0155 theGflashStep->GetPreStepPoint()->SetSensitiveDetector(aCurrentVolume->GetLogicalVolume()->GetSensitiveDetector());
0156 G4VSensitiveDetector *aSensitive = theGflashStep->GetPreStepPoint()->GetSensitiveDetector();
0157
0158 if (aSensitive == nullptr)
0159 continue;
0160
0161 G4String nameCalor = aCurrentVolume->GetName();
0162 nameCalor.assign(nameCalor, 0, 2);
0163 G4double samplingWeight = 1.0;
0164 if (nameCalor == "HB") {
0165 samplingWeight = Gflash::scaleSensitiveHB;
0166 } else if (nameCalor == "HE" || nameCalor == "HT") {
0167 samplingWeight = Gflash::scaleSensitiveHE;
0168 }
0169 theGflashStep->SetTotalEnergyDeposit(spotIter->getEnergy() * samplingWeight);
0170
0171 aSensitive->Hit(theGflashStep);
0172 }
0173 }
0174
0175 void GflashHadronShowerModel::updateGflashStep(const G4ThreeVector &spotPosition, G4double timeGlobal) {
0176 theGflashStep->GetPostStepPoint()->SetGlobalTime(timeGlobal);
0177 theGflashStep->GetPreStepPoint()->SetPosition(spotPosition);
0178 theGflashStep->GetPostStepPoint()->SetPosition(spotPosition);
0179 theGflashStep->GetPreStepPoint()->SetTouchableHandle(theGflashTouchableHandle);
0180 }
0181
0182 G4bool GflashHadronShowerModel::isFirstInelasticInteraction(const G4FastTrack &fastTrack) {
0183 G4bool isFirst = false;
0184
0185 G4StepPoint *preStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint();
0186 G4StepPoint *postStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint();
0187
0188 G4String procName = postStep->GetProcessDefinedStep()->GetProcessName();
0189 G4ParticleDefinition *particleType = fastTrack.GetPrimaryTrack()->GetDefinition();
0190
0191
0192
0193
0194 if ((particleType == G4PionPlus::PionPlusDefinition() && procName == "WrappedPionPlusInelastic") ||
0195 (particleType == G4PionMinus::PionMinusDefinition() && procName == "WrappedPionMinusInelastic") ||
0196 (particleType == G4KaonPlus::KaonPlusDefinition() && procName == "WrappedKaonPlusInelastic") ||
0197 (particleType == G4KaonMinus::KaonMinusDefinition() && procName == "WrappedKaonMinusInelastic") ||
0198 (particleType == G4AntiProton::AntiProtonDefinition() && procName == "WrappedAntiProtonInelastic") ||
0199 (particleType == G4Proton::ProtonDefinition() && procName == "WrappedProtonInelastic")) {
0200
0201
0202
0203
0204 const G4TrackVector *fSecondaryVector = fastTrack.GetPrimaryTrack()->GetStep()->GetSecondary();
0205 G4double leadingEnergy = 0.0;
0206
0207
0208
0209
0210
0211
0212 for (unsigned int isec = 0; isec < fSecondaryVector->size(); isec++) {
0213 G4Track *fSecondaryTrack = (*fSecondaryVector)[isec];
0214 G4double secondaryEnergy = fSecondaryTrack->GetKineticEnergy();
0215
0216 if (secondaryEnergy > leadingEnergy) {
0217 leadingEnergy = secondaryEnergy;
0218 }
0219 }
0220
0221 if ((preStep->GetTotalEnergy() != 0) && (leadingEnergy / preStep->GetTotalEnergy() < Gflash::QuasiElasticLike))
0222 isFirst = true;
0223
0224
0225
0226
0227 if (theHisto->getStoreFlag()) {
0228 theHisto->preStepPosition->Fill(preStep->GetPosition().getRho() / cm);
0229 theHisto->postStepPosition->Fill(postStep->GetPosition().getRho() / cm);
0230 theHisto->deltaStep->Fill((postStep->GetPosition() - preStep->GetPosition()).getRho() / cm);
0231 theHisto->kineticEnergy->Fill(fastTrack.GetPrimaryTrack()->GetKineticEnergy() / GeV);
0232 theHisto->energyLoss->Fill(fabs(fastTrack.GetPrimaryTrack()->GetStep()->GetDeltaEnergy() / GeV));
0233 theHisto->energyRatio->Fill(leadingEnergy / preStep->GetTotalEnergy());
0234 }
0235 }
0236 return isFirst;
0237 }
0238
0239 G4bool GflashHadronShowerModel::excludeDetectorRegion(const G4FastTrack &fastTrack) {
0240 const double invcm = 1.0 / CLHEP::cm;
0241 G4bool isExcluded = false;
0242 int verbosity = theParSet.getUntrackedParameter<int>("Verbosity");
0243
0244
0245
0246 G4double eta = fastTrack.GetPrimaryTrack()->GetPosition().pseudoRapidity();
0247 if (std::fabs(eta) > 1.392 && std::fabs(eta) < 1.566) {
0248 if (verbosity > 0) {
0249 edm::LogVerbatim("SimGeneralGFlash") << "GflashHadronShowerModel: excluding region of eta = " << eta;
0250 }
0251 return true;
0252 } else {
0253 const G4StepPoint *postStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint();
0254
0255 Gflash::CalorimeterNumber kCalor = Gflash::getCalorimeterNumber(postStep->GetPosition() * invcm);
0256 G4double distOut = 9999.0;
0257
0258
0259
0260 if (std::fabs(eta) > Gflash::EtaMin[Gflash::kENCA] &&
0261 std::fabs((postStep->GetPosition()).getZ() / CLHEP::cm) < Gflash::Zmin[Gflash::kENCA]) {
0262 return true;
0263 }
0264
0265
0266
0267
0268
0269
0270 if (kCalor == Gflash::kHB) {
0271 distOut = Gflash::Rmax[Gflash::kHB] - postStep->GetPosition().getRho() * invcm;
0272 if (distOut < Gflash::MinDistanceToOut)
0273 isExcluded = true;
0274 } else if (kCalor == Gflash::kHE) {
0275 distOut = Gflash::Zmax[Gflash::kHE] - std::fabs(postStep->GetPosition().getZ() * invcm);
0276 if (distOut < Gflash::MinDistanceToOut)
0277 isExcluded = true;
0278 }
0279
0280
0281 if (isExcluded && verbosity > 0) {
0282 G4cout << "GflashHadronShowerModel: skipping kCalor = " << kCalor << " DistanceToOut " << distOut << " from ("
0283 << (postStep->GetPosition()).getRho() * invcm << ":" << (postStep->GetPosition()).getZ() * invcm
0284 << ") of KE = " << fastTrack.GetPrimaryTrack()->GetKineticEnergy() / CLHEP::GeV << G4endl;
0285 }
0286 }
0287
0288 return isExcluded;
0289 }