From 9dcbc9251ca81c0f06e8ea0e60458fd13676cacf Mon Sep 17 00:00:00 2001 From: Vladimir Date: Mon, 15 Dec 2014 18:28:17 +0100 Subject: [PATCH] add optional hadronic GFlash to any Physics List --- .../interface/GFlashHadronShowerModel.h | 65 ++++ .../interface/ParametrisedEMPhysics.h | 7 +- .../Application/src/GFlashEMShowerModel.cc | 10 +- .../src/GFlashHadronShowerModel.cc | 315 ++++++++++++++++++ .../Application/src/ParametrisedEMPhysics.cc | 131 +++++--- SimG4Core/Physics/src/PhysicsList.cc | 13 +- 6 files changed, 472 insertions(+), 69 deletions(-) create mode 100644 SimG4Core/Application/interface/GFlashHadronShowerModel.h create mode 100644 SimG4Core/Application/src/GFlashHadronShowerModel.cc diff --git a/SimG4Core/Application/interface/GFlashHadronShowerModel.h b/SimG4Core/Application/interface/GFlashHadronShowerModel.h new file mode 100644 index 0000000000000..e659b7617f16a --- /dev/null +++ b/SimG4Core/Application/interface/GFlashHadronShowerModel.h @@ -0,0 +1,65 @@ +#ifndef GflashHadronShowerModel_H +#define GflashHadronShowerModel_H + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "G4VFastSimulationModel.hh" + +#include "G4TouchableHandle.hh" +#include "G4Navigator.hh" +#include "G4Step.hh" + +class GflashHadronShowerProfile; +class GflashPiKShowerProfile; +class GflashKaonPlusShowerProfile; +class GflashKaonMinusShowerProfile; +class GflashProtonShowerProfile; +class GflashAntiProtonShowerProfile; +class G4Region; +//class GflashHistogram; + +class GFlashHadronShowerModel : public G4VFastSimulationModel +{ +public: + //------------------------- + // Constructor, destructor + //------------------------- + GFlashHadronShowerModel (G4String modelName, G4Region* envelope, + const edm::ParameterSet& parSet); + ~GFlashHadronShowerModel (); + + //------------------------------------------------------------------------ + // Virtual methods that should be implemented for this hadron shower model + //------------------------------------------------------------------------ + + G4bool IsApplicable(const G4ParticleDefinition&); + G4bool ModelTrigger(const G4FastTrack &); + void DoIt(const G4FastTrack&, G4FastStep&); + +private: + G4bool isFirstInelasticInteraction(const G4FastTrack& fastTrack); + G4bool excludeDetectorRegion(const G4FastTrack& fastTrack); + void makeHits(const G4FastTrack& fastTrack); + void updateGflashStep(const G4ThreeVector& position, G4double time); + +private: + + G4bool theWatcherOn; + edm::ParameterSet theParSet; + GflashHadronShowerProfile *theProfile; + GflashPiKShowerProfile *thePiKProfile; + GflashKaonPlusShowerProfile *theKaonPlusProfile; + GflashKaonMinusShowerProfile *theKaonMinusProfile; + GflashProtonShowerProfile *theProtonProfile; + GflashAntiProtonShowerProfile *theAntiProtonProfile; + + const G4Region* theRegion; + + G4Step *theGflashStep; + G4Navigator *theGflashNavigator; + G4TouchableHandle theGflashTouchableHandle; + + //debugging histograms + //GflashHistogram* theHisto; +}; + +#endif diff --git a/SimG4Core/Application/interface/ParametrisedEMPhysics.h b/SimG4Core/Application/interface/ParametrisedEMPhysics.h index db7f47e12c86d..8ef93e19929cf 100644 --- a/SimG4Core/Application/interface/ParametrisedEMPhysics.h +++ b/SimG4Core/Application/interface/ParametrisedEMPhysics.h @@ -12,6 +12,7 @@ #include "G4VPhysicsConstructor.hh" class GFlashEMShowerModel; +class GFlashHadronShowerModel; class ElectronLimiter; class ParametrisedEMPhysics : public G4VPhysicsConstructor @@ -28,8 +29,10 @@ class ParametrisedEMPhysics : public G4VPhysicsConstructor edm::ParameterSet theParSet; - GFlashEMShowerModel *theEMShowerModel; - GFlashEMShowerModel *theHadShowerModel; + GFlashEMShowerModel *theEcalEMShowerModel; + GFlashEMShowerModel *theHcalEMShowerModel; + GFlashHadronShowerModel *theEcalHadShowerModel; + GFlashHadronShowerModel *theHcalHadShowerModel; ElectronLimiter *theElectronLimiter; ElectronLimiter *thePositronLimiter; diff --git a/SimG4Core/Application/src/GFlashEMShowerModel.cc b/SimG4Core/Application/src/GFlashEMShowerModel.cc index 9b690e07a4723..5a45dfcc4d55e 100644 --- a/SimG4Core/Application/src/GFlashEMShowerModel.cc +++ b/SimG4Core/Application/src/GFlashEMShowerModel.cc @@ -57,7 +57,8 @@ GFlashEMShowerModel::IsApplicable(const G4ParticleDefinition& particleType) G4bool GFlashEMShowerModel::ModelTrigger(const G4FastTrack & fastTrack ) { // Mininum energy cutoff to parameterize - if(fastTrack.GetPrimaryTrack()->GetKineticEnergy() < GeV) { return false; } + if(fastTrack.GetPrimaryTrack()->GetKineticEnergy() < Gflash::energyCutOff) + { return false; } if(excludeDetectorRegion(fastTrack)) { return false; } // This will be changed accordingly when the way @@ -133,12 +134,15 @@ void GFlashEMShowerModel::makeHits(const G4FastTrack& fastTrack) // Put touchable for each hit so that touchable history // keeps track of each step. - theGflashNavigator->LocateGlobalPointAndUpdateTouchableHandle(spotIter->getPosition(),G4ThreeVector(0,0,0),theGflashTouchableHandle, false); + theGflashNavigator->LocateGlobalPointAndUpdateTouchableHandle(spotIter->getPosition(), + G4ThreeVector(0,0,0), + theGflashTouchableHandle, false); updateGflashStep(spotIter->getPosition(),spotIter->getTime()); // If there is a watcher defined in a job and the flag is turned on if(theWatcherOn) { - SteppingAction* userSteppingAction = (SteppingAction*) G4EventManager::GetEventManager()->GetUserSteppingAction(); + SteppingAction* userSteppingAction = + (SteppingAction*) G4EventManager::GetEventManager()->GetUserSteppingAction(); userSteppingAction->m_g4StepSignal(theGflashStep); } diff --git a/SimG4Core/Application/src/GFlashHadronShowerModel.cc b/SimG4Core/Application/src/GFlashHadronShowerModel.cc new file mode 100644 index 0000000000000..0958ff824e198 --- /dev/null +++ b/SimG4Core/Application/src/GFlashHadronShowerModel.cc @@ -0,0 +1,315 @@ +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "SimG4Core/Application/interface/GFlashHadronShowerModel.h" +#include "SimG4Core/Application/interface/SteppingAction.h" + +#include "SimGeneral/GFlash/interface/GflashHadronShowerProfile.h" +#include "SimGeneral/GFlash/interface/GflashPiKShowerProfile.h" +#include "SimGeneral/GFlash/interface/GflashKaonPlusShowerProfile.h" +#include "SimGeneral/GFlash/interface/GflashKaonMinusShowerProfile.h" +#include "SimGeneral/GFlash/interface/GflashProtonShowerProfile.h" +#include "SimGeneral/GFlash/interface/GflashAntiProtonShowerProfile.h" +#include "SimGeneral/GFlash/interface/GflashNameSpace.h" +//#include "SimGeneral/GFlash/interface/GflashHistogram.h" +#include "SimGeneral/GFlash/interface/GflashHit.h" + +#include "G4FastSimulationManager.hh" +#include "G4TransportationManager.hh" +#include "G4TouchableHandle.hh" +#include "G4VSensitiveDetector.hh" +#include "G4VPhysicalVolume.hh" +#include "G4EventManager.hh" + +#include "G4PionMinus.hh" +#include "G4PionPlus.hh" +#include "G4KaonMinus.hh" +#include "G4KaonPlus.hh" +#include "G4Proton.hh" +#include "G4AntiProton.hh" +#include "G4VProcess.hh" + +#include + +using namespace CLHEP; + +GFlashHadronShowerModel::GFlashHadronShowerModel(G4String modelName, G4Region* envelope, + const edm::ParameterSet& parSet) + : G4VFastSimulationModel(modelName, envelope), theParSet(parSet) +{ + theWatcherOn = parSet.getParameter("watcherOn"); + + theProfile = 0; + thePiKProfile = new GflashPiKShowerProfile(parSet); + theKaonPlusProfile = new GflashKaonPlusShowerProfile(parSet); + theKaonMinusProfile = new GflashKaonMinusShowerProfile(parSet); + theProtonProfile = new GflashProtonShowerProfile(parSet); + theAntiProtonProfile = new GflashAntiProtonShowerProfile(parSet); + //theHisto = GflashHistogram::instance(); + + theRegion = const_cast(envelope); + + theGflashStep = new G4Step(); + theGflashTouchableHandle = new G4TouchableHistory(); + theGflashNavigator = new G4Navigator(); + +} + +GFlashHadronShowerModel::~GFlashHadronShowerModel() +{ + delete thePiKProfile; + delete theKaonPlusProfile; + delete theKaonMinusProfile; + delete theProtonProfile; + delete theAntiProtonProfile; + delete theGflashStep; +} + +G4bool GFlashHadronShowerModel::IsApplicable(const G4ParticleDefinition& particleType) +{ + return + &particleType == G4PionMinus::PionMinusDefinition() || + &particleType == G4PionPlus::PionPlusDefinition() || + &particleType == G4KaonMinus::KaonMinusDefinition() || + &particleType == G4KaonPlus::KaonPlusDefinition() || + &particleType == G4AntiProton::AntiProtonDefinition() || + &particleType == G4Proton::ProtonDefinition() ; +} + +G4bool GFlashHadronShowerModel::ModelTrigger(const G4FastTrack& fastTrack) +{ + // ModelTrigger returns false for Gflash Hadronic Shower Model if it is not + // tested from the corresponding wrapper process, GflashHadronWrapperProcess. + // Temporarily the track status is set to fPostponeToNextEvent at the wrapper + // process before ModelTrigger is really tested for the second time through + // PostStepGPIL of enviaged parameterization processes. The better + // implmentation may be using via G4VUserTrackInformation of each track, which + // requires to modify a geant source code of stepping (G4SteppingManager2) + + // mininum energy cutoff to parameterize + if (fastTrack.GetPrimaryTrack()->GetKineticEnergy() < GeV*Gflash::energyCutOff ) + { return false; } + + // This will be changed accordingly when the way + // dealing with CaloRegion changes later. + G4TouchableHistory* touch = + (G4TouchableHistory*)(fastTrack.GetPrimaryTrack()->GetTouchable()); + G4VPhysicalVolume* pCurrentVolume = touch->GetVolume(); + if( pCurrentVolume == 0) { return false; } + + G4LogicalVolume* lv = pCurrentVolume->GetLogicalVolume(); + if(lv->GetRegion() != theRegion) { return false; } + + // check whether this is called from the normal GPIL or the wrapper process GPIL + //if(fastTrack.GetPrimaryTrack()->GetTrackStatus() == fPostponeToNextEvent ) { + + // Shower pameterization start at the first inelastic interaction point + //if(isFirstInelasticInteraction(fastTrack) && excludeDetectorRegion(fastTrack)) + if(excludeDetectorRegion(fastTrack)) + { return false; } + + return true; +} + +void GFlashHadronShowerModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) +{ + // kill the particle + fastStep.KillPrimaryTrack(); + fastStep.ProposePrimaryTrackPathLength(0.0); + + // parameterize energy depostion by the particle type + G4ParticleDefinition* particleType = fastTrack.GetPrimaryTrack()->GetDefinition(); + + theProfile = thePiKProfile; + if(particleType == G4KaonMinus::KaonMinusDefinition()) theProfile = theKaonMinusProfile; + else if(particleType == G4KaonPlus::KaonPlusDefinition()) theProfile = theKaonPlusProfile; + else if(particleType == G4AntiProton::AntiProtonDefinition()) theProfile = theAntiProtonProfile; + else if(particleType == G4Proton::ProtonDefinition()) theProfile = theProtonProfile; + + //input variables for GflashHadronShowerProfile + G4double energy = fastTrack.GetPrimaryTrack()->GetKineticEnergy()/GeV; + G4double globalTime = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetGlobalTime(); + G4double charge = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint()->GetCharge(); + G4ThreeVector position = fastTrack.GetPrimaryTrack()->GetPosition()/cm; + G4ThreeVector momentum = fastTrack.GetPrimaryTrack()->GetMomentum()/GeV; + G4int showerType = Gflash::findShowerType(position); + + theProfile->initialize(showerType,energy,globalTime,charge,position,momentum); + theProfile->loadParameters(); + theProfile->hadronicParameterization(); + + // make hits + makeHits(fastTrack); + +} + +void GFlashHadronShowerModel::makeHits(const G4FastTrack& fastTrack) { + + std::vector& gflashHitList = theProfile->getGflashHitList(); + + theGflashStep->SetTrack(const_cast(fastTrack.GetPrimaryTrack())); + theGflashStep->GetPostStepPoint()->SetProcessDefinedStep(const_cast + (fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint()->GetProcessDefinedStep())); + theGflashNavigator->SetWorldVolume(G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume()); + + std::vector::const_iterator spotIter = gflashHitList.begin(); + std::vector::const_iterator spotIterEnd = gflashHitList.end(); + + for( ; spotIter != spotIterEnd; spotIter++){ + + theGflashNavigator->LocateGlobalPointAndUpdateTouchableHandle(spotIter->getPosition(), + G4ThreeVector(0,0,0), + theGflashTouchableHandle, false); + updateGflashStep(spotIter->getPosition(),spotIter->getTime()); + + // if there is a watcher defined in a job and the flag is turned on + if(theWatcherOn) { + theGflashStep->SetTotalEnergyDeposit(spotIter->getEnergy()); + SteppingAction* userSteppingAction = + (SteppingAction*) G4EventManager::GetEventManager()->GetUserSteppingAction(); + userSteppingAction->m_g4StepSignal(theGflashStep); + } + + G4VPhysicalVolume* aCurrentVolume = theGflashStep->GetPreStepPoint()->GetPhysicalVolume(); + if( aCurrentVolume == 0 ) { continue; } + + G4LogicalVolume* lv = aCurrentVolume->GetLogicalVolume(); + if(lv->GetRegion() != theRegion) { continue; } + + theGflashStep->GetPreStepPoint()->SetSensitiveDetector(aCurrentVolume->GetLogicalVolume()->GetSensitiveDetector()); + G4VSensitiveDetector* aSensitive = theGflashStep->GetPreStepPoint()->GetSensitiveDetector(); + + if( aSensitive == 0 ) continue; + + G4String nameCalor = aCurrentVolume->GetName(); + nameCalor.assign(nameCalor,0,2); + G4double samplingWeight = 1.0; + if(nameCalor == "HB" ) { + samplingWeight = Gflash::scaleSensitiveHB; + } + else if(nameCalor=="HE" || nameCalor == "HT") { + samplingWeight = Gflash::scaleSensitiveHE; + } + theGflashStep->SetTotalEnergyDeposit(spotIter->getEnergy()*samplingWeight); + + aSensitive->Hit(theGflashStep); + + } +} + +void GFlashHadronShowerModel::updateGflashStep(const G4ThreeVector& spotPosition, G4double timeGlobal) +{ + theGflashStep->GetPostStepPoint()->SetGlobalTime(timeGlobal); + theGflashStep->GetPreStepPoint()->SetPosition(spotPosition); + theGflashStep->GetPostStepPoint()->SetPosition(spotPosition); + theGflashStep->GetPreStepPoint()->SetTouchableHandle(theGflashTouchableHandle); +} + +G4bool GFlashHadronShowerModel::isFirstInelasticInteraction(const G4FastTrack& fastTrack) +{ + G4bool isFirst=false; + + G4StepPoint* preStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPreStepPoint(); + G4StepPoint* postStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint(); + + G4String procName = postStep->GetProcessDefinedStep()->GetProcessName(); + G4ParticleDefinition* particleType = fastTrack.GetPrimaryTrack()->GetDefinition(); + + //@@@ this part is still temporary and the cut for the variable ratio should be optimized later + + if((particleType == G4PionPlus::PionPlusDefinition() && procName == "WrappedPionPlusInelastic") || + (particleType == G4PionMinus::PionMinusDefinition() && procName == "WrappedPionMinusInelastic") || + (particleType == G4KaonPlus::KaonPlusDefinition() && procName == "WrappedKaonPlusInelastic") || + (particleType == G4KaonMinus::KaonMinusDefinition() && procName == "WrappedKaonMinusInelastic") || + (particleType == G4AntiProton::AntiProtonDefinition() && procName == "WrappedAntiProtonInelastic") || + (particleType == G4Proton::ProtonDefinition() && procName == "WrappedProtonInelastic")) { + + //skip to the second interaction if the first inelastic is a quasi-elastic like interaction + //@@@ the cut may be optimized later + + const G4TrackVector* fSecondaryVector = fastTrack.GetPrimaryTrack()->GetStep()->GetSecondary(); + G4double leadingEnergy = 0.0; + + //loop over 'all' secondaries including those produced by continuous processes. + //@@@may require an additional condition only for hadron interaction with the process name, + //but it will not change the result anyway + + for (unsigned int isec = 0 ; isec < fSecondaryVector->size() ; isec++) { + G4Track* fSecondaryTrack = (*fSecondaryVector)[isec]; + G4double secondaryEnergy = fSecondaryTrack->GetKineticEnergy(); + + if(secondaryEnergy > leadingEnergy ) { + leadingEnergy = secondaryEnergy; + } + } + + if((preStep->GetTotalEnergy()!=0) && + (leadingEnergy/preStep->GetTotalEnergy() < Gflash::QuasiElasticLike)) isFirst = true; + + //Fill debugging histograms and check information on secondaries - + //remove after final implimentation + /* + if(theHisto->getStoreFlag()) { + theHisto->preStepPosition->Fill(preStep->GetPosition().getRho()/cm); + theHisto->postStepPosition->Fill(postStep->GetPosition().getRho()/cm); + theHisto->deltaStep->Fill((postStep->GetPosition() - preStep->GetPosition()).getRho()/cm); + theHisto->kineticEnergy->Fill(fastTrack.GetPrimaryTrack()->GetKineticEnergy()/GeV); + theHisto->energyLoss->Fill(fabs(fastTrack.GetPrimaryTrack()->GetStep()->GetDeltaEnergy()/GeV)); + theHisto->energyRatio->Fill(leadingEnergy/preStep->GetTotalEnergy()); + } + */ + } + return isFirst; +} + +G4bool GFlashHadronShowerModel::excludeDetectorRegion(const G4FastTrack& fastTrack) +{ + G4bool isExcluded=false; + int verbosity = theParSet.getUntrackedParameter("Verbosity"); + + //exclude regions where geometry are complicated + //+- one supermodule around the EB/EE boundary: 1.479 +- 0.0174*5 + G4double eta = fastTrack.GetPrimaryTrack()->GetPosition().pseudoRapidity() ; + if(std::fabs(eta) > 1.392 && std::fabs(eta) < 1.566) { + if(verbosity>0) { + edm::LogInfo("SimG4CoreApplication") + << "GFlashHadronShowerModel: excluding region of eta = " << eta; + } + return true; + } + else { + G4StepPoint* postStep = fastTrack.GetPrimaryTrack()->GetStep()->GetPostStepPoint(); + + Gflash::CalorimeterNumber kCalor = Gflash::getCalorimeterNumber(postStep->GetPosition()/cm); + G4double distOut = 9999.0; + + //exclude the region where the shower starting point is inside the preshower + if( std::fabs(eta) > Gflash::EtaMin[Gflash::kENCA] && + std::fabs((postStep->GetPosition()).getZ()/cm) < Gflash::Zmin[Gflash::kENCA]) { + return true; + } + + //<---the shower starting point is always inside envelopes + //@@@exclude the region where the shower starting point is too close to the end of + //the hadronic envelopes (may need to be optimized further!) + //@@@if we extend parameterization including Magnet/HO, we need to change this strategy + if(kCalor == Gflash::kHB) { + distOut = Gflash::Rmax[Gflash::kHB] - postStep->GetPosition().getRho()/cm; + if (distOut < Gflash::MinDistanceToOut ) isExcluded = true; + } + else if(kCalor == Gflash::kHE) { + distOut = Gflash::Zmax[Gflash::kHE] - std::fabs(postStep->GetPosition().getZ()/cm); + if (distOut < Gflash::MinDistanceToOut ) isExcluded = true; + } + + //@@@remove this print statement later + if(isExcluded && verbosity > 0) { + edm::LogInfo("SimG4CoreApplication") + << "GFlashHadronShowerModel: skipping kCalor = " << kCalor << + " DistanceToOut " << distOut << " from (" << (postStep->GetPosition()).getRho()/cm << + ":" << (postStep->GetPosition()).getZ()/cm << ") of KE = " + << fastTrack.GetPrimaryTrack()->GetKineticEnergy()/GeV; + } + } + + return isExcluded; +} diff --git a/SimG4Core/Application/src/ParametrisedEMPhysics.cc b/SimG4Core/Application/src/ParametrisedEMPhysics.cc index 3811be089616c..3b576c7a59d54 100644 --- a/SimG4Core/Application/src/ParametrisedEMPhysics.cc +++ b/SimG4Core/Application/src/ParametrisedEMPhysics.cc @@ -7,6 +7,7 @@ #include "SimG4Core/Application/interface/ParametrisedEMPhysics.h" #include "SimG4Core/Application/interface/GFlashEMShowerModel.h" +#include "SimG4Core/Application/interface/GFlashHadronShowerModel.h" #include "SimG4Core/Application/interface/ElectronLimiter.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" @@ -21,6 +22,12 @@ #include "G4RegionStore.hh" #include "G4Electron.hh" #include "G4Positron.hh" +#include "G4PionMinus.hh" +#include "G4PionPlus.hh" +#include "G4KaonMinus.hh" +#include "G4KaonPlus.hh" +#include "G4Proton.hh" +#include "G4AntiProton.hh" #include "G4EmProcessOptions.hh" #include "G4PhysicsListHelper.hh" @@ -31,13 +38,17 @@ ParametrisedEMPhysics::ParametrisedEMPhysics(std::string name, const edm::ParameterSet & p) : G4VPhysicsConstructor(name), theParSet(p) { - theEMShowerModel = 0; - theHadShowerModel = 0; + theEcalEMShowerModel = 0; + theEcalHadShowerModel = 0; + theHcalEMShowerModel = 0; + theHcalHadShowerModel = 0; } ParametrisedEMPhysics::~ParametrisedEMPhysics() { - delete theEMShowerModel; - delete theHadShowerModel; + delete theEcalEMShowerModel; + delete theEcalHadShowerModel; + delete theHcalEMShowerModel; + delete theHcalHadShowerModel; } void ParametrisedEMPhysics::ConstructParticle() @@ -63,41 +74,58 @@ void ParametrisedEMPhysics::ConstructProcess() { // GFlash part bool gem = theParSet.getParameter("GflashEcal"); bool ghad = theParSet.getParameter("GflashHcal"); + bool gemHad = theParSet.getParameter("GflashEcalHad"); + bool ghadHad = theParSet.getParameter("GflashHcalHad"); - if(gem || ghad) { + G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper(); + if(gem || ghad || gemHad || ghadHad) { edm::LogInfo("SimG4CoreApplication") - << "ParametrisedEMPhysics: GFlash Construct: " - << gem << " " << ghad; + << "ParametrisedEMPhysics: GFlash Construct for e+-: " + << gem << " " << ghad << " for hadrons: " << gemHad << " " << ghadHad; + G4FastSimulationManagerProcess * theFastSimulationManagerProcess = new G4FastSimulationManagerProcess(); - aParticleIterator->reset(); - while ((*aParticleIterator)()) { - G4ParticleDefinition * particle = aParticleIterator->value(); - G4ProcessManager * pmanager = particle->GetProcessManager(); - G4String pname = particle->GetParticleName(); - if(pname == "e-" || pname == "e+") { - pmanager->AddProcess(theFastSimulationManagerProcess, -1, -1, 1); - } + + if(gem || ghad) { + ph->RegisterProcess(theFastSimulationManagerProcess, G4Electron::Electron()); + ph->RegisterProcess(theFastSimulationManagerProcess, G4Positron::Positron()); + } + if(gemHad || ghadHad) { + ph->RegisterProcess(theFastSimulationManagerProcess, G4Proton::Proton()); + ph->RegisterProcess(theFastSimulationManagerProcess, G4AntiProton::AntiProton()); + ph->RegisterProcess(theFastSimulationManagerProcess, G4PionPlus::PionPlus()); + ph->RegisterProcess(theFastSimulationManagerProcess, G4PionMinus::PionMinus()); + ph->RegisterProcess(theFastSimulationManagerProcess, G4KaonPlus::KaonPlus()); + ph->RegisterProcess(theFastSimulationManagerProcess, G4KaonMinus::KaonMinus()); } - if(gem) { + if(gem || gemHad) { G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("EcalRegion"); - + if(!aRegion){ edm::LogInfo("SimG4CoreApplication") << "ParametrisedEMPhysics::ConstructProcess: " << "EcalRegion is not defined, GFlash will not be enabled for ECAL!"; } else { - - //Electromagnetic Shower Model for ECAL - theEMShowerModel = - new GFlashEMShowerModel("GflashEMShowerModel",aRegion,theParSet); - //std::cout << "GFlash is defined for EcalRegion" << std::endl; - } + if(gem) { + + //Electromagnetic Shower Model for ECAL + theEcalEMShowerModel = + new GFlashEMShowerModel("GflashEcalEMShowerModel",aRegion,theParSet); + //std::cout << "GFlash is defined for EcalRegion" << std::endl; + } + if(gemHad) { + + //Electromagnetic Shower Model for ECAL + theEcalHadShowerModel = + new GFlashHadronShowerModel("GflashEcalHadShowerModel",aRegion,theParSet); + //std::cout << "GFlash is defined for EcalRegion" << std::endl; + } + } } - if(ghad) { + if(ghad || ghadHad) { G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion"); if(!aRegion) { @@ -106,15 +134,24 @@ void ParametrisedEMPhysics::ConstructProcess() { << "HcalRegion is not defined, GFlash will not be enabled for HCAL!"; } else { - - //Electromagnetic Shower Model for HCAL - theHadShowerModel = - new GFlashEMShowerModel("GflashHadShowerModel",aRegion,theParSet); - //std::cout << "GFlash is defined for HcalRegion" << std::endl; + if(ghad) { + + //Electromagnetic Shower Model for HCAL + theHcalEMShowerModel = + new GFlashEMShowerModel("GflashHcalEMShowerModel",aRegion,theParSet); + //std::cout << "GFlash is defined for HcalRegion" << std::endl; + } + if(ghadHad) { + + //Electromagnetic Shower Model for ECAL + theHcalHadShowerModel = + new GFlashHadronShowerModel("GflashHcalHadShowerModel",aRegion,theParSet); + //std::cout << "GFlash is defined for EcalRegion" << std::endl; + } } } } - // bremsstrahlung threshold + // bremsstrahlung threshold and EM verbosity G4EmProcessOptions opt; G4int verb = theParSet.getUntrackedParameter("Verbosity",0); opt.SetVerbose(verb - 1); @@ -124,15 +161,14 @@ void ParametrisedEMPhysics::ConstructProcess() { << "ParametrisedEMPhysics::ConstructProcess: bremsstrahlung threshold Eth= " << bremth/GeV << " GeV"; opt.SetBremsstrahlungTh(bremth); - // (m_pPhysics.getParameter("G4BremsstrahlungThreshold")*GeV); - // Russian Roulette part + + // Russian roulette and tracking cut for e+- const G4int NREG = 6; const G4String rname[NREG] = {"EcalRegion", "HcalRegion", "MuonIron", "PreshowerRegion","CastorRegion", "DefaultRegionForTheWorld"}; G4double rrfact[NREG] = { 1.0 }; - // Russian roulette for e- double energyLim = theParSet.getParameter("RusRoElectronEnergyLimit")*MeV; if(energyLim > 0.0) { @@ -146,11 +182,6 @@ void ParametrisedEMPhysics::ConstructProcess() { if(rrfact[i] < 1.0) { opt.ActivateSecondaryBiasing("eIoni",rname[i],rrfact[i],energyLim); opt.ActivateSecondaryBiasing("hIoni",rname[i],rrfact[i],energyLim); - //opt.ActivateSecondaryBiasing("muIoni",rname[i],rrfact[i],energyLim); - //opt.ActivateSecondaryBiasing("ionIoni",rname[i],rrfact[i],energyLim); - //opt.ActivateSecondaryBiasingForGamma("phot",rname[i],rrfact[i],energyLim); - //opt.ActivateSecondaryBiasingForGamma("compt",rname[i],rrfact[i],energyLim); - //opt.ActivateSecondaryBiasingForGamma("conv",rname[i],rrfact[i],energyLim); edm::LogInfo("SimG4CoreApplication") << "ParametrisedEMPhysics: Russian Roulette" << " for e- Prob= " << rrfact[i] @@ -165,21 +196,17 @@ void ParametrisedEMPhysics::ConstructProcess() { bool rLimiter = theParSet.getParameter("ElectronRangeTest"); bool pLimiter = theParSet.getParameter("PositronStepLimit"); - if(eLimiter || rLimiter || pLimiter) { - G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper(); - - if(eLimiter || rLimiter) { - theElectronLimiter = new ElectronLimiter(theParSet); - theElectronLimiter->SetRangeCheckFlag(rLimiter); - theElectronLimiter->SetFieldCheckFlag(eLimiter); - ph->RegisterProcess(theElectronLimiter, G4Electron::Electron()); - } + if(eLimiter || rLimiter) { + theElectronLimiter = new ElectronLimiter(theParSet); + theElectronLimiter->SetRangeCheckFlag(rLimiter); + theElectronLimiter->SetFieldCheckFlag(eLimiter); + ph->RegisterProcess(theElectronLimiter, G4Electron::Electron()); + } - if(pLimiter){ - thePositronLimiter = new ElectronLimiter(theParSet); - thePositronLimiter->SetFieldCheckFlag(pLimiter); - ph->RegisterProcess(theElectronLimiter, G4Positron::Positron()); - } + if(pLimiter){ + thePositronLimiter = new ElectronLimiter(theParSet); + thePositronLimiter->SetFieldCheckFlag(pLimiter); + ph->RegisterProcess(theElectronLimiter, G4Positron::Positron()); } // enable fluorescence bool fluo = theParSet.getParameter("FlagFluo"); diff --git a/SimG4Core/Physics/src/PhysicsList.cc b/SimG4Core/Physics/src/PhysicsList.cc index a7207351ab9e7..15729ddf833f6 100644 --- a/SimG4Core/Physics/src/PhysicsList.cc +++ b/SimG4Core/Physics/src/PhysicsList.cc @@ -15,13 +15,7 @@ PhysicsList::PhysicsList(G4LogicalVolumeToDDLogicalPartMap & map, } PhysicsList::~PhysicsList() { - /* - if (m_Verbosity > 1) - LogDebug("Physics") << " G4BremsstrahlungThreshold was " - << G4LossTableManager::Instance()->BremsstrahlungTh()/GeV - << " GeV "; - */ - if (prodCuts!=0) delete prodCuts; + delete prodCuts; } void PhysicsList::SetCuts() { @@ -29,16 +23,11 @@ void PhysicsList::SetCuts() { SetDefaultCutValue(m_pPhysics.getParameter("DefaultCutValue")*cm); SetCutsWithDefault(); - // VI potentially unsafe to set bremsstrahlung threashold here - //G4LossTableManager::Instance()->SetBremsstrahlungTh - // (m_pPhysics.getParameter("G4BremsstrahlungThreshold")*GeV); - if ( m_pPhysics.getParameter("CutsPerRegion") ) { prodCuts->update(); } if ( m_Verbosity > 1) { - //G4LossTableManager::Instance()->SetVerbose(m_Verbosity-1); G4VUserPhysicsList::DumpCutValuesTable(); }