From 9062c31dbaee8594cfb9e2a16d0f1af12e8d0590 Mon Sep 17 00:00:00 2001 From: Vladimir Date: Fri, 30 May 2014 11:03:41 +0200 Subject: [PATCH] added options needed for neutron background study, tuned exotica physics list --- .../python/NeutronBGforMuonsHP_cff.py | 3 + .../python/NeutronBGforMuonsXS_cff.py | 3 + SimG4Core/Application/python/g4SimHits_cfi.py | 5 +- .../Application/src/ParametrisedEMPhysics.cc | 18 ++++- .../test/python/runCuModel_cfg.py | 6 +- .../interface/CustomPhysicsList.h | 24 +++---- .../python/Exotica_HSCP_SIM_cfi.py | 8 +-- SimG4Core/CustomPhysics/src/CustomPhysics.cc | 6 +- .../CustomPhysics/src/CustomPhysicsList.cc | 71 +++++++++---------- .../Physics/interface/DDG4ProductionCuts.h | 5 +- SimG4Core/Physics/interface/PhysicsList.h | 2 +- SimG4Core/Physics/src/DDG4ProductionCuts.cc | 14 ++-- SimG4Core/Physics/src/PhysicsList.cc | 2 +- 13 files changed, 96 insertions(+), 71 deletions(-) diff --git a/SimG4Core/Application/python/NeutronBGforMuonsHP_cff.py b/SimG4Core/Application/python/NeutronBGforMuonsHP_cff.py index c6f9c6eb617fb..8d315fe34a5f1 100644 --- a/SimG4Core/Application/python/NeutronBGforMuonsHP_cff.py +++ b/SimG4Core/Application/python/NeutronBGforMuonsHP_cff.py @@ -10,6 +10,9 @@ def customise(process): process.common_maximum_time.DeadRegions = cms.vstring() # Physics List XS process.g4SimHits.Physics.type = cms.string('SimG4Core/Physics/FTFP_BERT_HP_EML') + process.g4SimHits.Physics.CutsOnProton = cms.bool(False) + process.g4SimHits.Physics.FlagMuNucl = cms.bool(True) + process.g4SimHits.Physics.FlagFluo = cms.bool(True) # Eta cut process.g4SimHits.Generator.MinEtaCut = cms.double(-7.0) process.g4SimHits.Generator.MaxEtaCut = cms.double(7.0) diff --git a/SimG4Core/Application/python/NeutronBGforMuonsXS_cff.py b/SimG4Core/Application/python/NeutronBGforMuonsXS_cff.py index c52b3e4bd07e0..11312d057a40b 100644 --- a/SimG4Core/Application/python/NeutronBGforMuonsXS_cff.py +++ b/SimG4Core/Application/python/NeutronBGforMuonsXS_cff.py @@ -10,6 +10,9 @@ def customise(process): process.common_maximum_time.DeadRegions = cms.vstring() # Physics List XS process.g4SimHits.Physics.type = cms.string('SimG4Core/Physics/FTFP_BERT_XS_EML') + process.g4SimHits.Physics.CutsOnProton = cms.bool(False) + process.g4SimHits.Physics.FlagMuNucl = cms.bool(True) + process.g4SimHits.Physics.FlagFluo = cms.bool(True) # Eta cut process.g4SimHits.Generator.MinEtaCut = cms.double(-7.0) process.g4SimHits.Generator.MaxEtaCut = cms.double(7.0) diff --git a/SimG4Core/Application/python/g4SimHits_cfi.py b/SimG4Core/Application/python/g4SimHits_cfi.py index ab425d58e968c..536efec5287b2 100644 --- a/SimG4Core/Application/python/g4SimHits_cfi.py +++ b/SimG4Core/Application/python/g4SimHits_cfi.py @@ -83,6 +83,7 @@ type = cms.string('SimG4Core/Physics/QGSP_FTFP_BERT_EML'), DummyEMPhysics = cms.bool(False), CutsPerRegion = cms.bool(True), + CutsOnProton = cms.bool(True), DefaultCutValue = cms.double(1.0), ## cuts in cm G4BremsstrahlungThreshold = cms.double(0.5), ## cut in GeV Verbosity = cms.untracked.int32(0), @@ -97,12 +98,10 @@ TrackingCut = cms.bool(False), SRType = cms.bool(True), FlagMuNucl = cms.bool(False), + FlagFluo = cms.bool(False), EMPhysics = cms.untracked.bool(True), HadPhysics = cms.untracked.bool(True), FlagBERT = cms.untracked.bool(False), - FlagFTF = cms.untracked.bool(False), - FlagGlauber = cms.untracked.bool(False), - FlagHP = cms.untracked.bool(False), GflashEcal = cms.bool(False), bField = cms.double(3.8), energyScaleEB = cms.double(1.032), diff --git a/SimG4Core/Application/src/ParametrisedEMPhysics.cc b/SimG4Core/Application/src/ParametrisedEMPhysics.cc index 6f24a68db0a5d..f037fde6535d7 100644 --- a/SimG4Core/Application/src/ParametrisedEMPhysics.cc +++ b/SimG4Core/Application/src/ParametrisedEMPhysics.cc @@ -25,6 +25,8 @@ #include "G4EmProcessOptions.hh" #include "G4PhysicsListHelper.hh" #include "G4SystemOfUnits.hh" +#include "G4UAtomicDeexcitation.hh" +#include "G4LossTableManager.hh" ParametrisedEMPhysics::ParametrisedEMPhysics(std::string name, const edm::ParameterSet & p) : G4VPhysicsConstructor(name), theParSet(p) @@ -79,7 +81,8 @@ void ParametrisedEMPhysics::ConstructProcess() { } if(gem) { - G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("EcalRegion"); + G4Region* aRegion = + G4RegionStore::GetInstance()->GetRegion("EcalRegion"); if(!aRegion){ edm::LogInfo("SimG4CoreApplication") @@ -95,7 +98,8 @@ void ParametrisedEMPhysics::ConstructProcess() { } } if(ghad) { - G4Region* aRegion = G4RegionStore::GetInstance()->GetRegion("HcalRegion"); + G4Region* aRegion = + G4RegionStore::GetInstance()->GetRegion("HcalRegion"); if(!aRegion) { edm::LogInfo("SimG4CoreApplication") << "ParametrisedEMPhysics::ConstructProcess: " @@ -119,7 +123,8 @@ void ParametrisedEMPhysics::ConstructProcess() { G4double rrfact[NREG] = { 1.0 }; // Russian roulette for e- - double energyLim = theParSet.getParameter("RusRoElectronEnergyLimit")*MeV; + double energyLim = + theParSet.getParameter("RusRoElectronEnergyLimit")*MeV; if(energyLim > 0.0) { rrfact[0] = theParSet.getParameter("RusRoEcalElectron"); rrfact[1] = theParSet.getParameter("RusRoHcalElectron"); @@ -166,4 +171,11 @@ void ParametrisedEMPhysics::ConstructProcess() { ph->RegisterProcess(theElectronLimiter, G4Positron::Positron()); } } + // enable fluorescence + bool fluo = theParSet.getParameter("FlagFluo"); + if(fluo) { + G4VAtomDeexcitation* de = new G4UAtomicDeexcitation(); + G4LossTableManager::Instance()->SetAtomDeexcitation(de); + de->SetFluo(true); + } } diff --git a/SimG4Core/CheckSecondary/test/python/runCuModel_cfg.py b/SimG4Core/CheckSecondary/test/python/runCuModel_cfg.py index 5487fa4c9a1ed..a056b85a4a330 100644 --- a/SimG4Core/CheckSecondary/test/python/runCuModel_cfg.py +++ b/SimG4Core/CheckSecondary/test/python/runCuModel_cfg.py @@ -102,14 +102,14 @@ G4BremsstrahlungThreshold = cms.double(0.5), DefaultCutValue = cms.double(1.0), CutsPerRegion = cms.bool(True), + CutsOnProton = cms.bool(True), Verbosity = cms.untracked.int32(0), EMPhysics = cms.untracked.bool(False), HadPhysics = cms.untracked.bool(True), QuasiElastic = cms.untracked.bool(True), FlagBERT = cms.untracked.bool(False), - FlagCHIPS = cms.untracked.bool(False), - FlagFTF = cms.untracked.bool(False), - FlagGlauber = cms.untracked.bool(False), + FlagMuNucl = cms.bool(False), + FlagFluo = cms.bool(False), Model = cms.untracked.string('LEP'), type = cms.string('SimG4Core/Physics/CMSModel'), DummyEMPhysics = cms.bool(False) diff --git a/SimG4Core/CustomPhysics/interface/CustomPhysicsList.h b/SimG4Core/CustomPhysics/interface/CustomPhysicsList.h index 695bd0057dd8b..11208326511b1 100644 --- a/SimG4Core/CustomPhysics/interface/CustomPhysicsList.h +++ b/SimG4Core/CustomPhysics/interface/CustomPhysicsList.h @@ -15,25 +15,25 @@ class CustomPhysicsList : public G4VPhysicsConstructor public: CustomPhysicsList(std::string name, const edm::ParameterSet & p); virtual ~CustomPhysicsList(); + + virtual void ConstructParticle(); + virtual void ConstructProcess(); + protected: - virtual void ConstructParticle(); - virtual void ConstructProcess(); - void addCustomPhysics(); -// void SetCuts(); -private: + void addCustomPhysics(); - void setupRHadronPhycis(G4ParticleDefinition* particle); - void setupSUSYPhycis(G4ParticleDefinition* particle); +private: - //HadronicProcessHelper *myHelper; - G4ProcessHelper *myHelper; + void setupRHadronPhycis(G4ParticleDefinition* particle); + void setupSUSYPhycis(G4ParticleDefinition* particle); - edm::ParameterSet myConfig; + G4ProcessHelper *myHelper; - std::string particleDefFilePath; - std::string processDefFilePath; + edm::ParameterSet myConfig; + std::string particleDefFilePath; + std::string processDefFilePath; }; #endif diff --git a/SimG4Core/CustomPhysics/python/Exotica_HSCP_SIM_cfi.py b/SimG4Core/CustomPhysics/python/Exotica_HSCP_SIM_cfi.py index a70d7a83ea8a3..86b7afae0f958 100644 --- a/SimG4Core/CustomPhysics/python/Exotica_HSCP_SIM_cfi.py +++ b/SimG4Core/CustomPhysics/python/Exotica_HSCP_SIM_cfi.py @@ -32,15 +32,14 @@ def customise(process): G4BremsstrahlungThreshold = cms.double(0.5), ## cut in GeV DefaultCutValue = cms.double(1.), ## cuts in cm,default 1cm CutsPerRegion = cms.bool(True), + CutsOnProton = cms.bool(True), Verbosity = cms.untracked.int32(0), type = cms.string('SimG4Core/Physics/CustomPhysics'), EMPhysics = cms.untracked.bool(True), ##G4 default true HadPhysics = cms.untracked.bool(True), ##G4 default true FlagBERT = cms.untracked.bool(False), - FlagCHIPS = cms.untracked.bool(False), - FlagFTF = cms.untracked.bool(False), - FlagGlauber = cms.untracked.bool(False), - FlagHP = cms.untracked.bool(False), + FlagMuNucl = cms.bool(False), + FlagFluo = cms.bool(False), GFlash = cms.PSet( GflashHistogram = cms.bool(False), GflashEMShowerModel = cms.bool(False), @@ -57,6 +56,7 @@ def customise(process): G4BremsstrahlungThreshold = cms.double(0.5), ## cut in GeV DefaultCutValue = cms.double(1.), ## cuts in cm,default 1cm CutsPerRegion = cms.bool(True), + CutsOnProton = cms.bool(True), Verbosity = cms.untracked.int32(0), type = cms.string('SimG4Core/Physics/CustomPhysics'), ) diff --git a/SimG4Core/CustomPhysics/src/CustomPhysics.cc b/SimG4Core/CustomPhysics/src/CustomPhysics.cc index 146f29b44911a..4ae9242a4d5ea 100644 --- a/SimG4Core/CustomPhysics/src/CustomPhysics.cc +++ b/SimG4Core/CustomPhysics/src/CustomPhysics.cc @@ -1,6 +1,7 @@ #include "SimG4Core/CustomPhysics/interface/CustomPhysics.h" #include "SimG4Core/CustomPhysics/interface/CustomPhysicsList.h" #include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysics95msc93.h" +#include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysics.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "G4DecayPhysics.hh" @@ -29,7 +30,8 @@ CustomPhysics::CustomPhysics(G4LogicalVolumeToDDLogicalPartMap& map, << hadPhys << "\n"; // EM Physics - RegisterPhysics(new CMSEmStandardPhysics95msc93("EM standard msc93",ver,"")); + RegisterPhysics(new CMSEmStandardPhysics(ver)); + //RegisterPhysics(new CMSEmStandardPhysics95msc93("EM standard msc93",ver,"")); // Synchroton Radiation & GN Physics RegisterPhysics(new G4EmExtraPhysics(ver)); @@ -50,7 +52,7 @@ CustomPhysics::CustomPhysics(G4LogicalVolumeToDDLogicalPartMap& map, RegisterPhysics(new G4IonPhysics(ver)); // Neutron tracking cut - RegisterPhysics( new G4NeutronTrackingCut(ver)); + RegisterPhysics(new G4NeutronTrackingCut(ver)); // Custom Physics RegisterPhysics(new CustomPhysicsList("custom",p)); diff --git a/SimG4Core/CustomPhysics/src/CustomPhysicsList.cc b/SimG4Core/CustomPhysics/src/CustomPhysicsList.cc index 164c9ad1b87ac..bcabea9818be2 100644 --- a/SimG4Core/CustomPhysics/src/CustomPhysicsList.cc +++ b/SimG4Core/CustomPhysics/src/CustomPhysicsList.cc @@ -23,16 +23,16 @@ using namespace CLHEP; -CustomPhysicsList::CustomPhysicsList(std::string name, const edm::ParameterSet & p) : G4VPhysicsConstructor(name) { - +CustomPhysicsList::CustomPhysicsList(std::string name, const edm::ParameterSet & p) + : G4VPhysicsConstructor(name) +{ myConfig = p; edm::FileInPath fp = p.getParameter("particlesDef"); particleDefFilePath = fp.fullPath(); edm::LogInfo("CustomPhysics")<<"Path for custom particle definition file: " <reset(); - while((*aParticleIterator)()) { - int i = 0; + while((*aParticleIterator)()) { G4ParticleDefinition* particle = aParticleIterator->value(); CustomParticle* cp = dynamic_cast(particle); if(CustomParticleFactory::isCustomParticle(particle)) { @@ -60,14 +59,19 @@ void CustomPhysicsList::addCustomPhysics(){ << " is Custom. Mass is " <GetPDGMass()/GeV <<" GeV."; if(cp->GetCloud()!=0) { - LogDebug("CustomPhysics")<<"Cloud mass is " - <GetCloud()->GetPDGMass()/GeV - <<" GeV. Spectator mass is " - <(particle)->GetSpectator()->GetPDGMass()/GeV - <<" GeV."; + LogDebug("CustomPhysics") + <<"Cloud mass is " + <GetCloud()->GetPDGMass()/GeV + <<" GeV. Spectator mass is " + <(particle)->GetSpectator()->GetPDGMass()/GeV + <<" GeV."; } G4ProcessManager* pmanager = particle->GetProcessManager(); if(pmanager) { + if(particle->GetPDGCharge()/eplus != 0) { + pmanager->AddProcess(new G4hMultipleScattering,-1, 1, 1); + pmanager->AddProcess(new G4hIonisation, -1, 2, 2); + } if(cp!=0) { if(particle->GetParticleType()=="rhadron" || particle->GetParticleType()=="mesonino" || @@ -76,10 +80,6 @@ void CustomPhysicsList::addCustomPhysics(){ pmanager->AddDiscreteProcess(new FullModelHadronicProcess(myHelper)); } } - if(particle->GetPDGCharge()/eplus != 0) { - pmanager->AddProcess(new G4hMultipleScattering,-1, 1,i+1); - pmanager->AddProcess(new G4hIonisation, -1, 2,i+2); - } } else LogDebug("CustomPhysics") << " No pmanager"; } @@ -87,42 +87,41 @@ void CustomPhysicsList::addCustomPhysics(){ } -void CustomPhysicsList::setupRHadronPhycis(G4ParticleDefinition* particle){ - +void CustomPhysicsList::setupRHadronPhycis(G4ParticleDefinition* particle) +{ // LogDebug("CustomPhysics")<<"Configuring rHadron: " // < CustomParticle* cp = dynamic_cast(particle); - if(cp->GetCloud()!=0) - LogDebug("CustomPhysics")<<"Cloud mass is " - <GetCloud()->GetPDGMass()/GeV - <<" GeV. Spectator mass is " - <(particle)->GetSpectator()->GetPDGMass()/GeV - <<" GeV."; - + if(cp->GetCloud()!=0) { + LogDebug("CustomPhysics") + <<"Cloud mass is " + <GetCloud()->GetPDGMass()/GeV + <<" GeV. Spectator mass is " + <(particle)->GetSpectator()->GetPDGMass()/GeV + <<" GeV."; + } G4ProcessManager* pmanager = particle->GetProcessManager(); if(pmanager){ if(!myHelper) myHelper = new G4ProcessHelper(myConfig); - pmanager->AddDiscreteProcess(new FullModelHadronicProcess(myHelper)); //GHEISHA if(particle->GetPDGCharge()/eplus != 0){ - pmanager->AddProcess(new G4hMultipleScattering,-1, 1,1); - pmanager->AddProcess(new G4hIonisation, -1, 2,2); + pmanager->AddProcess(new G4hMultipleScattering,-1, 1, 1); + pmanager->AddProcess(new G4hIonisation, -1, 2, 2); } + pmanager->AddDiscreteProcess(new FullModelHadronicProcess(myHelper)); //GHEISHA } else LogDebug("CustomPhysics") << " No pmanager"; } - -void CustomPhysicsList::setupSUSYPhycis(G4ParticleDefinition* particle){ - -// CustomParticle* cp = dynamic_cast(particle); +void CustomPhysicsList::setupSUSYPhycis(G4ParticleDefinition* particle) +{ G4ProcessManager* pmanager = particle->GetProcessManager(); if(pmanager){ - pmanager->AddProcess(new G4Decay,1, 1,1); if(particle->GetPDGCharge()/eplus != 0){ - pmanager->AddProcess(new G4hMultipleScattering,-1, 2,2); - pmanager->AddProcess(new G4hIonisation, -1, 3,3); + pmanager->AddProcess(new G4hMultipleScattering,-1, 1, 1); + pmanager->AddProcess(new G4hIonisation, -1, 2, 2); } + pmanager->AddProcess(new G4Decay, 1, -1, 3); } else LogDebug("CustomPhysics") << " No pmanager"; } diff --git a/SimG4Core/Physics/interface/DDG4ProductionCuts.h b/SimG4Core/Physics/interface/DDG4ProductionCuts.h index f366cd8de336c..5fc13cd2c4407 100644 --- a/SimG4Core/Physics/interface/DDG4ProductionCuts.h +++ b/SimG4Core/Physics/interface/DDG4ProductionCuts.h @@ -1,6 +1,7 @@ #ifndef SimG4Core_DDG4ProductionCuts_H #define SimG4Core_DDG4ProductionCuts_H +#include "FWCore/ParameterSet/interface/ParameterSet.h" #include "SimG4Core/Geometry/interface/G4LogicalVolumeToDDLogicalPartMap.h" #include @@ -14,7 +15,8 @@ class G4ProductionCuts; class DDG4ProductionCuts { public: - DDG4ProductionCuts(const G4LogicalVolumeToDDLogicalPartMap&, int); + DDG4ProductionCuts(const G4LogicalVolumeToDDLogicalPartMap&, int, + const edm::ParameterSet & p); ~DDG4ProductionCuts(); void update(); void SetVerbosity( int verb ) { m_Verbosity = verb; return ; } @@ -29,6 +31,7 @@ class DDG4ProductionCuts { G4LogicalVolumeToDDLogicalPartMap map_; std::string m_KeywordRegion; int m_Verbosity; + bool m_protonCut; G4LogicalVolumeToDDLogicalPartMap::Vector vec_; }; diff --git a/SimG4Core/Physics/interface/PhysicsList.h b/SimG4Core/Physics/interface/PhysicsList.h index a7d267c567078..7155f831def09 100644 --- a/SimG4Core/Physics/interface/PhysicsList.h +++ b/SimG4Core/Physics/interface/PhysicsList.h @@ -21,7 +21,7 @@ class PhysicsList : public G4VModularPhysicsList { virtual void SetCuts(); private: - edm::ParameterSet m_pPhysics; + const edm::ParameterSet m_pPhysics; DDG4ProductionCuts * prodCuts; int m_Verbosity; }; diff --git a/SimG4Core/Physics/src/DDG4ProductionCuts.cc b/SimG4Core/Physics/src/DDG4ProductionCuts.cc index 68ce87adbd865..c46ac7c65929e 100644 --- a/SimG4Core/Physics/src/DDG4ProductionCuts.cc +++ b/SimG4Core/Physics/src/DDG4ProductionCuts.cc @@ -7,9 +7,11 @@ #include -DDG4ProductionCuts::DDG4ProductionCuts(const G4LogicalVolumeToDDLogicalPartMap& map, int verb) : map_(map), m_Verbosity(verb) { - m_KeywordRegion = "CMSCutsRegion"; - initialize(); +DDG4ProductionCuts::DDG4ProductionCuts(const G4LogicalVolumeToDDLogicalPartMap& map, int verb, const edm::ParameterSet & p) + : map_(map), m_Verbosity(verb) { + m_KeywordRegion = "CMSCutsRegion"; + m_protonCut = p.getUntrackedParameter("CutsOnProton",true); + initialize(); } DDG4ProductionCuts::~DDG4ProductionCuts(){ @@ -108,11 +110,12 @@ void DDG4ProductionCuts::setProdCuts(const DDLogicalPart lpart, // // search for production cuts - // you must have three of them: e+ e- gamma + // you must have four of them: e+ e- gamma proton // double gammacut; double electroncut; double positroncut; + double protoncut = 0.0; int temp = map_.toDouble("ProdCutsForGamma",lpart,gammacut); if (temp != 1){ throw SimG4Exception("DDG4ProductionCuts: Problem with Region tags: no/more than one ProdCutsForGamma."); @@ -133,7 +136,8 @@ void DDG4ProductionCuts::setProdCuts(const DDLogicalPart lpart, prodCuts->SetProductionCut( electroncut, idxG4ElectronCut ); prodCuts->SetProductionCut( positroncut, idxG4PositronCut ); // For recoil use the same cut as for e- - prodCuts->SetProductionCut( electroncut, idxG4ProtonCut ); + if(m_protonCut) { protoncut = electroncut; } + prodCuts->SetProductionCut( protoncut, idxG4ProtonCut ); if ( m_Verbosity > 0 ) { LogDebug("Physics") << "DDG4ProductionCuts : Setting cuts for " << regionName << "\n Electrons: " << electroncut diff --git a/SimG4Core/Physics/src/PhysicsList.cc b/SimG4Core/Physics/src/PhysicsList.cc index 705cddb96c715..6bb935bd5682c 100644 --- a/SimG4Core/Physics/src/PhysicsList.cc +++ b/SimG4Core/Physics/src/PhysicsList.cc @@ -11,7 +11,7 @@ PhysicsList::PhysicsList(G4LogicalVolumeToDDLogicalPartMap & map, const edm::ParameterSet & p) : G4VModularPhysicsList(), m_pPhysics(p), prodCuts(0) { m_Verbosity = m_pPhysics.getUntrackedParameter("Verbosity",0); - prodCuts = new DDG4ProductionCuts(map, m_Verbosity); + prodCuts = new DDG4ProductionCuts(map, m_Verbosity, m_pPhysics); } PhysicsList::~PhysicsList() {