File indexing completed on 2024-04-06 12:30:14
0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002
0003 #include "SimG4Core/Application/interface/GFlashHadronShowerModel.h"
0004 #include "SimG4Core/Application/interface/SteppingAction.h"
0005 #include "SimG4Core/Geometry/interface/DD4hep2DDDName.h"
0006
0007 #include "SimGeneral/GFlash/interface/GflashHadronShowerProfile.h"
0008 #include "SimGeneral/GFlash/interface/GflashPiKShowerProfile.h"
0009 #include "SimGeneral/GFlash/interface/GflashKaonPlusShowerProfile.h"
0010 #include "SimGeneral/GFlash/interface/GflashKaonMinusShowerProfile.h"
0011 #include "SimGeneral/GFlash/interface/GflashProtonShowerProfile.h"
0012 #include "SimGeneral/GFlash/interface/GflashAntiProtonShowerProfile.h"
0013 #include "SimGeneral/GFlash/interface/GflashNameSpace.h"
0014 #include "SimGeneral/GFlash/interface/GflashHit.h"
0015
0016 #include "G4FastSimulationManager.hh"
0017 #include "G4TransportationManager.hh"
0018 #include "G4TouchableHandle.hh"
0019 #include "G4VSensitiveDetector.hh"
0020 #include "G4VPhysicalVolume.hh"
0021 #include "G4EventManager.hh"
0022
0023 #include "G4PionMinus.hh"
0024 #include "G4PionPlus.hh"
0025 #include "G4KaonMinus.hh"
0026 #include "G4KaonPlus.hh"
0027 #include "G4Proton.hh"
0028 #include "G4AntiProton.hh"
0029 #include "G4VProcess.hh"
0030
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
0072
0073
0074
0075
0076
0077
0078
0079
0080 if (fastTrack.GetPrimaryTrack()->GetKineticEnergy() < CLHEP::GeV * Gflash::energyCutOff) {
0081 return false;
0082 }
0083
0084
0085
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
0100 fastStep.KillPrimaryTrack();
0101 fastStep.ProposePrimaryTrackPathLength(0.0);
0102
0103
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
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
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
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)(DD4hep2DDDName::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
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
0211
0212
0213 const G4TrackVector* fSecondaryVector = fastTrack.GetPrimaryTrack()->GetStep()->GetSecondary();
0214 G4double leadingEnergy = 0.0;
0215
0216
0217
0218
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
0241
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
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
0261
0262
0263
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
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 }