diff --git a/SimG4Core/Application/python/NeutronBGforMuonsHP_cff.py b/SimG4Core/Application/python/NeutronBGforMuonsHP_cff.py new file mode 100644 index 0000000000000..898447ea1844d --- /dev/null +++ b/SimG4Core/Application/python/NeutronBGforMuonsHP_cff.py @@ -0,0 +1,34 @@ +import FWCore.ParameterSet.Config as cms + +def customise(process): + + # fragment allowing to simulate neutron background in muon system + + if hasattr(process,'g4SimHits'): + # time window 10 millisecond + process.common_maximum_time.MaxTrackTime = cms.double(100000000.0) + process.common_maximum_time.MaxTimeNames = cms.vstring('ZDCRegion') + process.common_maximum_time.MaxTrackTimes = cms.vdouble(2000) + # Physics List XS + process.g4SimHits.Physics.type = cms.string('SimG4Core/Physics/FTFP_BERT_HP_EML') + process.g4SimHits.Physics.CutsOnProton = cms.untracked.bool(False) + process.g4SimHits.Physics.FlagFluo = cms.bool(True) + process.g4SimHits.Physics.RusRoGammaEnergyLimit = cms.double(0.0) + # Eta cut + process.g4SimHits.Generator.MinEtaCut = cms.double(-7.0) + process.g4SimHits.Generator.MaxEtaCut = cms.double(7.0) + # stacking action + process.g4SimHits.StackingAction.MaxTrackTime = cms.double(100000000.0) + process.g4SimHits.StackingAction.MaxTimeNames = cms.vstring('ZDCRegion') + process.g4SimHits.StackingAction.MaxTrackTimes = cms.vdouble(2000) + # stepping action + process.g4SimHits.SteppingAction.MaxTrackTime = cms.double(100000000.0) + process.g4SimHits.SteppingAction.MaxTimeNames = cms.vstring('ZDCRegion') + process.g4SimHits.SteppingAction.MaxTrackTimes = cms.vdouble(2000) + # Russian roulette disabled + process.g4SimHits.StackingAction.RusRoGammaEnergyLimit = cms.double(0.0) + process.g4SimHits.StackingAction.RusRoNeutronEnergyLimit = cms.double(0.0) + # full simulation of HF + process.g4SimHits.HFShower.UseHFGflash = cms.bool(False) + + return(process) diff --git a/SimG4Core/Application/python/NeutronBGforMuonsXS_cff.py b/SimG4Core/Application/python/NeutronBGforMuonsXS_cff.py new file mode 100644 index 0000000000000..9b31cb15b33a1 --- /dev/null +++ b/SimG4Core/Application/python/NeutronBGforMuonsXS_cff.py @@ -0,0 +1,34 @@ +import FWCore.ParameterSet.Config as cms + +def customise(process): + + # fragment allowing to simulate neutron background in muon system + + if hasattr(process,'g4SimHits'): + # time window 10 millisecond + process.common_maximum_time.MaxTrackTime = cms.double(100000000.0) + process.common_maximum_time.MaxTimeNames = cms.vstring('ZDCRegion') + process.common_maximum_time.MaxTrackTimes = cms.vdouble(2000) + # Physics List XS + process.g4SimHits.Physics.type = cms.string('SimG4Core/Physics/FTFP_BERT_XS_EML') + process.g4SimHits.Physics.CutsOnProton = cms.untracked.bool(False) + process.g4SimHits.Physics.FlagFluo = cms.bool(True) + process.g4SimHits.Physics.RusRoGammaEnergyLimit = cms.double(0.0) + # Eta cut + process.g4SimHits.Generator.MinEtaCut = cms.double(-7.0) + process.g4SimHits.Generator.MaxEtaCut = cms.double(7.0) + # stacking action + process.g4SimHits.StackingAction.MaxTrackTime = cms.double(100000000.0) + process.g4SimHits.StackingAction.MaxTimeNames = cms.vstring('ZDCRegion') + process.g4SimHits.StackingAction.MaxTrackTimes = cms.vdouble(2000) + # stepping action + process.g4SimHits.SteppingAction.MaxTrackTime = cms.double(100000000.0) + process.g4SimHits.SteppingAction.MaxTimeNames = cms.vstring('ZDCRegion') + process.g4SimHits.SteppingAction.MaxTrackTimes = cms.vdouble(2000) + # Russian roulette disabled + process.g4SimHits.StackingAction.RusRoGammaEnergyLimit = cms.double(0.0) + process.g4SimHits.StackingAction.RusRoNeutronEnergyLimit = cms.double(0.0) + # full simulation of HF + process.g4SimHits.HFShower.UseHFGflash = cms.bool(False) + + return(process) diff --git a/SimG4Core/Application/python/g4SimHits_cfi.py b/SimG4Core/Application/python/g4SimHits_cfi.py index f70ab05d5742d..6df2644eb7a06 100644 --- a/SimG4Core/Application/python/g4SimHits_cfi.py +++ b/SimG4Core/Application/python/g4SimHits_cfi.py @@ -77,6 +77,7 @@ delta = cms.double(1.0) ), Physics = cms.PSet( + common_maximum_time, # NOTE : if you want EM Physics only, # please select "SimG4Core/Physics/DummyPhysics" for type # and turn ON DummyEMPhysics @@ -84,6 +85,7 @@ type = cms.string('SimG4Core/Physics/QGSP_FTFP_BERT_EML'), DummyEMPhysics = cms.bool(False), CutsPerRegion = cms.bool(True), + CutsOnProton = cms.untracked.bool(True), DefaultCutValue = cms.double(1.0), ## cuts in cm G4BremsstrahlungThreshold = cms.double(0.5), ## cut in GeV Verbosity = cms.untracked.int32(0), @@ -97,6 +99,7 @@ Region = cms.string(' '), TrackingCut = cms.bool(True), SRType = cms.bool(True), + FlagFluo = cms.bool(False), EMPhysics = cms.untracked.bool(True), HadPhysics = cms.untracked.bool(True), FlagBERT = cms.untracked.bool(False), diff --git a/SimG4Core/Application/src/ParametrisedEMPhysics.cc b/SimG4Core/Application/src/ParametrisedEMPhysics.cc index 3047902b0d604..38d140f792d68 100644 --- a/SimG4Core/Application/src/ParametrisedEMPhysics.cc +++ b/SimG4Core/Application/src/ParametrisedEMPhysics.cc @@ -20,6 +20,8 @@ #include "G4RegionStore.hh" #include "G4EmProcessOptions.hh" +#include "G4UAtomicDeexcitation.hh" +#include "G4LossTableManager.hh" ParametrisedEMPhysics::ParametrisedEMPhysics(std::string name, const edm::ParameterSet & p) : G4VPhysicsConstructor(name), theParSet(p) @@ -162,4 +164,11 @@ void ParametrisedEMPhysics::ConstructProcess() { } } } + // enable fluorescence + bool fluo = theParSet.getParameter("FlagFluo"); + if(fluo) { + G4VAtomDeexcitation* de = new G4UAtomicDeexcitation(); + G4LossTableManager::Instance()->SetAtomDeexcitation(de); + de->SetFluo(true); + } } 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/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 16157f966e106..6bb935bd5682c 100644 --- a/SimG4Core/Physics/src/PhysicsList.cc +++ b/SimG4Core/Physics/src/PhysicsList.cc @@ -3,6 +3,7 @@ #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "G4LossTableManager.hh" +#include "G4SystemOfUnits.hh" PhysicsList::PhysicsList(G4LogicalVolumeToDDLogicalPartMap & map, const HepPDT::ParticleDataTable * table_, @@ -10,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() { diff --git a/SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsXS.h b/SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsXS.h new file mode 100644 index 0000000000000..b4298834a1937 --- /dev/null +++ b/SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsXS.h @@ -0,0 +1,26 @@ +#ifndef SimG4Core_PhysicsLists_CMSEmStandardPhysicsXS_h +#define SimG4Core_PhysicsLists_CMSEmStandardPhysicsXS_h + +#include "G4VPhysicsConstructor.hh" +#include "globals.hh" + +class CMSEmStandardPhysicsXS : public G4VPhysicsConstructor { + +public: + CMSEmStandardPhysicsXS(G4int ver); + virtual ~CMSEmStandardPhysicsXS(); + + virtual void ConstructParticle(); + virtual void ConstructProcess(); + +private: + G4int verbose; +}; + +#endif + + + + + + diff --git a/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_HP_EML.cc b/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_HP_EML.cc new file mode 100644 index 0000000000000..af8047a510235 --- /dev/null +++ b/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_HP_EML.cc @@ -0,0 +1,71 @@ +#include "FTFPCMS_BERT_HP_EML.hh" +#include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsXS.h" +#include "SimG4Core/PhysicsLists/interface/CMSMonopolePhysics.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "G4DecayPhysics.hh" +#include "G4EmExtraPhysics.hh" +#include "G4IonPhysics.hh" +#include "G4StoppingPhysics.hh" +#include "G4HadronElasticPhysics.hh" +#include "G4NeutronTrackingCut.hh" +#include "G4HadronicProcessStore.hh" + +#include "G4DataQuestionaire.hh" +#include "HadronPhysicsFTFP_BERT_HP.hh" + +FTFPCMS_BERT_HP_EML::FTFPCMS_BERT_HP_EML(G4LogicalVolumeToDDLogicalPartMap& map, + const HepPDT::ParticleDataTable * table_, + sim::FieldBuilder *fieldBuilder_, + const edm::ParameterSet & p) : PhysicsList(map, table_, fieldBuilder_, p) { + + G4DataQuestionaire it(photon); + + int ver = p.getUntrackedParameter("Verbosity",0); + bool emPhys = p.getUntrackedParameter("EMPhysics",true); + bool hadPhys = p.getUntrackedParameter("HadPhysics",true); + bool tracking= p.getParameter("TrackingCut"); + //bool munucl = p.getParameter("FlagMuNucl"); + edm::LogInfo("PhysicsList") << "You are using the simulation engine: " + << "FTFP_BERT_EML with Flags for EM Physics " + << emPhys << ", for Hadronic Physics " + << hadPhys << " and tracking cut " << tracking; + + if (emPhys) { + // EM Physics + RegisterPhysics( new CMSEmStandardPhysicsXS(ver)); + + // Synchroton Radiation & GN Physics + G4EmExtraPhysics* gn = new G4EmExtraPhysics(ver); + //if(munucl) { G4String yes = "on"; gn->MuonNuclear(yes); } + RegisterPhysics(gn); + } + + // Decays + this->RegisterPhysics( new G4DecayPhysics(ver) ); + + if (hadPhys) { + G4HadronicProcessStore::Instance()->SetVerbose(ver); + + // Hadron Elastic scattering + RegisterPhysics( new G4HadronElasticPhysics(ver)); + + // Hadron Physics + RegisterPhysics( new HadronPhysicsFTFP_BERT_HP(ver)); + + // Stopping Physics + RegisterPhysics( new G4StoppingPhysics(ver)); + + // Ion Physics + RegisterPhysics( new G4IonPhysics(ver)); + + // Neutron tracking cut + if (tracking) { + RegisterPhysics( new G4NeutronTrackingCut(ver)); + } + } + + // Monopoles + RegisterPhysics( new CMSMonopolePhysics(table_,fieldBuilder_,p)); +} + diff --git a/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_HP_EML.hh b/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_HP_EML.hh new file mode 100644 index 0000000000000..04a0b7c7d71ec --- /dev/null +++ b/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_HP_EML.hh @@ -0,0 +1,16 @@ +#ifndef SimG4Core_PhysicsLists_FTFPCMS_BERT_HP_EML_H +#define SimG4Core_PhysicsLists_FTFPCMS_BERT_HP_EML_H + +#include "SimG4Core/Physics/interface/PhysicsList.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +class FTFPCMS_BERT_HP_EML: public PhysicsList { + +public: + FTFPCMS_BERT_HP_EML(G4LogicalVolumeToDDLogicalPartMap& map, const HepPDT::ParticleDataTable * table_, sim::FieldBuilder *fieldBuilder_, const edm::ParameterSet & p); +}; + +#endif + + + diff --git a/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_XS_EML.cc b/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_XS_EML.cc new file mode 100644 index 0000000000000..693bffdf0427d --- /dev/null +++ b/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_XS_EML.cc @@ -0,0 +1,74 @@ +#include "FTFPCMS_BERT_XS_EML.hh" +#include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsXS.h" +#include "SimG4Core/PhysicsLists/interface/CMSMonopolePhysics.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "G4DecayPhysics.hh" +#include "G4EmExtraPhysics.hh" +#include "G4IonPhysics.hh" +#include "G4StoppingPhysics.hh" +#include "G4HadronElasticPhysicsXS.hh" +#include "G4NeutronCrossSectionXS.hh" +#include "G4NeutronTrackingCut.hh" +#include "G4HadronicProcessStore.hh" + +#include "G4DataQuestionaire.hh" +#include "HadronPhysicsFTFP_BERT.hh" + +FTFPCMS_BERT_XS_EML::FTFPCMS_BERT_XS_EML(G4LogicalVolumeToDDLogicalPartMap& map, + const HepPDT::ParticleDataTable * table_, + sim::FieldBuilder *fieldBuilder_, + const edm::ParameterSet & p) : PhysicsList(map, table_, fieldBuilder_, p) { + + G4DataQuestionaire it(photon); + + int ver = p.getUntrackedParameter("Verbosity",0); + bool emPhys = p.getUntrackedParameter("EMPhysics",true); + bool hadPhys = p.getUntrackedParameter("HadPhysics",true); + bool tracking= p.getParameter("TrackingCut"); + //bool munucl = p.getParameter("FlagMuNucl"); + edm::LogInfo("PhysicsList") << "You are using the simulation engine: " + << "FTFP_BERT_EML with Flags for EM Physics " + << emPhys << ", for Hadronic Physics " + << hadPhys << " and tracking cut " << tracking; + + if (emPhys) { + // EM Physics + RegisterPhysics( new CMSEmStandardPhysicsXS(ver)); + + // Synchroton Radiation & GN Physics + G4EmExtraPhysics* gn = new G4EmExtraPhysics(ver); + //if(munucl) { G4String yes = "on"; gn->MuonNuclear(yes); } + RegisterPhysics(gn); + } + + // Decays + this->RegisterPhysics( new G4DecayPhysics(ver) ); + + if (hadPhys) { + G4HadronicProcessStore::Instance()->SetVerbose(ver); + + // Hadron Elastic scattering + RegisterPhysics( new G4HadronElasticPhysicsXS(ver)); + + // Hadron Physics + RegisterPhysics( new HadronPhysicsFTFP_BERT(ver)); + + // Stopping Physics + RegisterPhysics( new G4StoppingPhysics(ver)); + + // Ion Physics + RegisterPhysics( new G4IonPhysics(ver)); + + RegisterPhysics( new G4NeutronCrossSectionXS(ver)); + + // Neutron tracking cut + if (tracking) { + RegisterPhysics( new G4NeutronTrackingCut(ver)); + } + } + + // Monopoles + RegisterPhysics( new CMSMonopolePhysics(table_,fieldBuilder_,p)); +} + diff --git a/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_XS_EML.hh b/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_XS_EML.hh new file mode 100644 index 0000000000000..0c3a5e06ecc07 --- /dev/null +++ b/SimG4Core/PhysicsLists/plugins/FTFPCMS_BERT_XS_EML.hh @@ -0,0 +1,16 @@ +#ifndef SimG4Core_PhysicsLists_FTFPCMS_BERT_XS_EML_H +#define SimG4Core_PhysicsLists_FTFPCMS_BERT_XS_EML_H + +#include "SimG4Core/Physics/interface/PhysicsList.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +class FTFPCMS_BERT_XS_EML: public PhysicsList { + +public: + FTFPCMS_BERT_XS_EML(G4LogicalVolumeToDDLogicalPartMap& map, const HepPDT::ParticleDataTable * table_, sim::FieldBuilder *fieldBuilder_, const edm::ParameterSet & p); +}; + +#endif + + + diff --git a/SimG4Core/PhysicsLists/plugins/module.cc b/SimG4Core/PhysicsLists/plugins/module.cc index aa9fe1d63d807..50f9346ada431 100644 --- a/SimG4Core/PhysicsLists/plugins/module.cc +++ b/SimG4Core/PhysicsLists/plugins/module.cc @@ -24,6 +24,9 @@ #include "QGSPCMS_FTFP_BERT_EML95.hh" #include "QGSPCMS_FTFP_BERT_EML_New.hh" #include "QGSPCMS_FTFP_BERT_EML95msc93.hh" +#include "FTFPCMS_BERT_HP_EML.hh" +#include "FTFPCMS_BERT_XS_EML.hh" + DEFINE_PHYSICSLIST(CMSModel); DEFINE_PHYSICSLIST(DummyPhysics); @@ -33,6 +36,10 @@ typedef FTFPCMS_BERT FTFP_BERT; DEFINE_PHYSICSLIST(FTFP_BERT); typedef FTFPCMS_BERT_EML FTFP_BERT_EML; DEFINE_PHYSICSLIST(FTFP_BERT_EML); +typedef FTFPCMS_BERT_HP_EML FTFP_BERT_HP_EML; +DEFINE_PHYSICSLIST(FTFP_BERT_HP_EML); +typedef FTFPCMS_BERT_XS_EML FTFP_BERT_XS_EML; +DEFINE_PHYSICSLIST(FTFP_BERT_XS_EML); typedef FTFPCMS_BERT_EML95 FTFP_BERT_EML95; DEFINE_PHYSICSLIST(FTFP_BERT_EML95); typedef FTFPCMS_BERT_EMV FTFP_BERT_EMV; diff --git a/SimG4Core/PhysicsLists/src/CMSEmStandardPhysicsXS.cc b/SimG4Core/PhysicsLists/src/CMSEmStandardPhysicsXS.cc new file mode 100644 index 0000000000000..c4a3288bb5acd --- /dev/null +++ b/SimG4Core/PhysicsLists/src/CMSEmStandardPhysicsXS.cc @@ -0,0 +1,323 @@ +#include "SimG4Core/PhysicsLists/interface/CMSEmStandardPhysicsXS.h" + +#include "G4ParticleDefinition.hh" +#include "G4LossTableManager.hh" +#include "G4EmProcessOptions.hh" + +#include "G4ComptonScattering.hh" +#include "G4GammaConversion.hh" +#include "G4PhotoElectricEffect.hh" +#include "G4LivermorePhotoElectricModel.hh" +#include "G4KleinNishinaModel.hh" + +#include "G4hMultipleScattering.hh" +#include "G4eMultipleScattering.hh" +#include "G4MuMultipleScattering.hh" +#include "G4CoulombScattering.hh" +#include "G4eCoulombScatteringModel.hh" +#include "G4WentzelVIModel.hh" +#include "G4UrbanMscModel93.hh" + +#include "G4eIonisation.hh" +#include "G4eBremsstrahlung.hh" +#include "G4eplusAnnihilation.hh" + +#include "G4MuIonisation.hh" +#include "G4MuBremsstrahlung.hh" +#include "G4MuPairProduction.hh" + +#include "G4hIonisation.hh" +#include "G4ionIonisation.hh" +#include "G4hBremsstrahlung.hh" +#include "G4hPairProduction.hh" + +#include "G4Gamma.hh" +#include "G4Electron.hh" +#include "G4Positron.hh" +#include "G4MuonPlus.hh" +#include "G4MuonMinus.hh" +#include "G4TauMinus.hh" +#include "G4TauPlus.hh" +#include "G4PionPlus.hh" +#include "G4PionMinus.hh" +#include "G4KaonPlus.hh" +#include "G4KaonMinus.hh" +#include "G4BMesonMinus.hh" +#include "G4BMesonPlus.hh" +#include "G4DMesonMinus.hh" +#include "G4DMesonPlus.hh" +#include "G4Proton.hh" +#include "G4AntiProton.hh" +#include "G4SigmaMinus.hh" +#include "G4AntiSigmaMinus.hh" +#include "G4SigmaPlus.hh" +#include "G4AntiSigmaPlus.hh" +#include "G4XiMinus.hh" +#include "G4AntiXiMinus.hh" +#include "G4OmegaMinus.hh" +#include "G4AntiOmegaMinus.hh" +#include "G4LambdacPlus.hh" +#include "G4AntiLambdacPlus.hh" +#include "G4XicPlus.hh" +#include "G4AntiXicPlus.hh" +#include "G4Deuteron.hh" +#include "G4Triton.hh" +#include "G4He3.hh" +#include "G4Alpha.hh" +#include "G4GenericIon.hh" + +#include "G4PhysicsListHelper.hh" +#include "G4BuilderType.hh" + +#include "G4SystemOfUnits.hh" + +CMSEmStandardPhysicsXS::CMSEmStandardPhysicsXS(G4int ver) : + G4VPhysicsConstructor("CMSEmStandardXS_opt1"), verbose(ver) { + G4LossTableManager::Instance(); + SetPhysicsType(bElectromagnetic); +} + +CMSEmStandardPhysicsXS::~CMSEmStandardPhysicsXS() {} + +void CMSEmStandardPhysicsXS::ConstructParticle() { + // gamma + G4Gamma::Gamma(); + + // leptons + G4Electron::Electron(); + G4Positron::Positron(); + G4MuonPlus::MuonPlus(); + G4MuonMinus::MuonMinus(); + G4TauMinus::TauMinusDefinition(); + G4TauPlus::TauPlusDefinition(); + + // mesons + G4PionPlus::PionPlusDefinition(); + G4PionMinus::PionMinusDefinition(); + G4KaonPlus::KaonPlusDefinition(); + G4KaonMinus::KaonMinusDefinition(); + G4DMesonMinus::DMesonMinusDefinition(); + G4DMesonPlus::DMesonPlusDefinition(); + G4BMesonMinus::BMesonMinusDefinition(); + G4BMesonPlus::BMesonPlusDefinition(); + + // barions + G4Proton::Proton(); + G4AntiProton::AntiProton(); + G4SigmaMinus::SigmaMinusDefinition(); + G4AntiSigmaMinus::AntiSigmaMinusDefinition(); + G4SigmaPlus::SigmaPlusDefinition(); + G4AntiSigmaPlus::AntiSigmaPlusDefinition(); + G4XiMinus::XiMinusDefinition(); + G4AntiXiMinus::AntiXiMinusDefinition(); + G4OmegaMinus::OmegaMinusDefinition(); + G4AntiOmegaMinus::AntiOmegaMinusDefinition(); + G4LambdacPlus::LambdacPlusDefinition(); + G4AntiLambdacPlus::AntiLambdacPlusDefinition(); + G4XicPlus::XicPlusDefinition(); + G4AntiXicPlus::AntiXicPlusDefinition(); + + // ions + G4Deuteron::Deuteron(); + G4Triton::Triton(); + G4He3::He3(); + G4Alpha::Alpha(); + G4GenericIon::GenericIonDefinition(); +} + +void CMSEmStandardPhysicsXS::ConstructProcess() { + + if(verbose > 0) { + G4cout << "### " << GetPhysicsName() << " Construct Processes " << G4endl; + } + + // This EM builder takes default models of Geant4 10 EMV. + // Multiple scattering by Urban for all particles + // except e+e- below 100 MeV for which the Urban93 model is used + + G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper(); + + // muon & hadron bremsstrahlung and pair production + G4MuBremsstrahlung* mub = new G4MuBremsstrahlung(); + G4MuPairProduction* mup = new G4MuPairProduction(); + G4hBremsstrahlung* pib = new G4hBremsstrahlung(); + G4hPairProduction* pip = new G4hPairProduction(); + G4hBremsstrahlung* kb = new G4hBremsstrahlung(); + G4hPairProduction* kp = new G4hPairProduction(); + G4hBremsstrahlung* pb = new G4hBremsstrahlung(); + G4hPairProduction* pp = new G4hPairProduction(); + + // muon & hadron multiple scattering + G4MuMultipleScattering* mumsc = new G4MuMultipleScattering(); + mumsc->AddEmModel(0, new G4WentzelVIModel()); + G4MuMultipleScattering* pimsc = new G4MuMultipleScattering(); + pimsc->AddEmModel(0, new G4WentzelVIModel()); + G4MuMultipleScattering* kmsc = new G4MuMultipleScattering(); + kmsc->AddEmModel(0, new G4WentzelVIModel()); + G4MuMultipleScattering* pmsc = new G4MuMultipleScattering(); + pmsc->AddEmModel(0, new G4WentzelVIModel()); + G4hMultipleScattering* hmsc = new G4hMultipleScattering("ionmsc"); + + // high energy limit for e+- scattering models and bremsstrahlung + G4double highEnergyLimit = 100*MeV; + + theParticleIterator->reset(); + while( (*theParticleIterator)() ){ + G4ParticleDefinition* particle = theParticleIterator->value(); + G4String particleName = particle->GetParticleName(); + + if (particleName == "gamma") { + + G4PhotoElectricEffect* photo = new G4PhotoElectricEffect(); + photo->SetEmModel(new G4LivermorePhotoElectricModel(),1); + ph->RegisterProcess(photo, particle); + G4ComptonScattering* compt = new G4ComptonScattering(); + compt->SetEmModel(new G4KleinNishinaModel(), 1); + ph->RegisterProcess(compt, particle); + ph->RegisterProcess(new G4GammaConversion(), particle); + + } else if (particleName == "e-") { + + G4eIonisation* eioni = new G4eIonisation(); + eioni->SetStepFunction(0.8, 1.0*mm); + + G4eMultipleScattering* msc = new G4eMultipleScattering; + msc->SetStepLimitType(fMinimal); + G4UrbanMscModel93* msc1 = new G4UrbanMscModel93(); + G4WentzelVIModel* msc2 = new G4WentzelVIModel(); + msc1->SetHighEnergyLimit(highEnergyLimit); + msc2->SetLowEnergyLimit(highEnergyLimit); + msc->AddEmModel(0, msc1); + msc->AddEmModel(0, msc2); + + G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel(); + G4CoulombScattering* ss = new G4CoulombScattering(); + ss->SetEmModel(ssm, 1); + ss->SetMinKinEnergy(highEnergyLimit); + ssm->SetLowEnergyLimit(highEnergyLimit); + ssm->SetActivationLowEnergyLimit(highEnergyLimit); + + ph->RegisterProcess(msc, particle); + ph->RegisterProcess(eioni, particle); + ph->RegisterProcess(new G4eBremsstrahlung(), particle); + ph->RegisterProcess(ss, particle); + + } else if (particleName == "e+") { + + G4eIonisation* eioni = new G4eIonisation(); + eioni->SetStepFunction(0.8, 1.0*mm); + + G4eMultipleScattering* msc = new G4eMultipleScattering; + msc->SetStepLimitType(fMinimal); + G4UrbanMscModel93* msc1 = new G4UrbanMscModel93(); + G4WentzelVIModel* msc2 = new G4WentzelVIModel(); + msc1->SetHighEnergyLimit(highEnergyLimit); + msc2->SetLowEnergyLimit(highEnergyLimit); + msc->AddEmModel(0, msc1); + msc->AddEmModel(0, msc2); + + G4eCoulombScatteringModel* ssm = new G4eCoulombScatteringModel(); + G4CoulombScattering* ss = new G4CoulombScattering(); + ss->SetEmModel(ssm, 1); + ss->SetMinKinEnergy(highEnergyLimit); + ssm->SetLowEnergyLimit(highEnergyLimit); + ssm->SetActivationLowEnergyLimit(highEnergyLimit); + + ph->RegisterProcess(msc, particle); + ph->RegisterProcess(eioni, particle); + ph->RegisterProcess(new G4eBremsstrahlung(), particle); + ph->RegisterProcess(new G4eplusAnnihilation(), particle); + ph->RegisterProcess(ss, particle); + + } else if (particleName == "mu+" || + particleName == "mu-" ) { + + ph->RegisterProcess(mumsc, particle); + ph->RegisterProcess(new G4MuIonisation(), particle); + ph->RegisterProcess(mub, particle); + ph->RegisterProcess(mup, particle); + ph->RegisterProcess(new G4CoulombScattering(), particle); + + } else if (particleName == "alpha" || + particleName == "He3" ) { + + //ph->RegisterProcess(hmsc, particle); + ph->RegisterProcess(new G4hMultipleScattering(), particle); + ph->RegisterProcess(new G4ionIonisation(), particle); + + } else if (particleName == "GenericIon") { + + ph->RegisterProcess(hmsc, particle); + ph->RegisterProcess(new G4ionIonisation(), particle); + + } else if (particleName == "pi+" || + particleName == "pi-" ) { + + //G4hMultipleScattering* pimsc = new G4hMultipleScattering(); + ph->RegisterProcess(pimsc, particle); + ph->RegisterProcess(new G4hIonisation(), particle); + ph->RegisterProcess(pib, particle); + ph->RegisterProcess(pip, particle); + ph->RegisterProcess(new G4CoulombScattering(), particle); + + } else if (particleName == "kaon+" || + particleName == "kaon-" ) { + + //G4hMultipleScattering* kmsc = new G4hMultipleScattering(); + ph->RegisterProcess(kmsc, particle); + ph->RegisterProcess(new G4hIonisation(), particle); + ph->RegisterProcess(kb, particle); + ph->RegisterProcess(kp, particle); + ph->RegisterProcess(new G4CoulombScattering(), particle); + + // } else if (particleName == "proton" ) { + } else if (particleName == "proton" || + particleName == "anti_proton") { + + //G4hMultipleScattering* pmsc = new G4hMultipleScattering(); + ph->RegisterProcess(pmsc, particle); + ph->RegisterProcess(new G4hIonisation(), particle); + ph->RegisterProcess(pb, particle); + ph->RegisterProcess(pp, particle); + ph->RegisterProcess(new G4CoulombScattering(), particle); + + } else if (particleName == "B+" || + particleName == "B-" || + particleName == "D+" || + particleName == "D-" || + particleName == "Ds+" || + particleName == "Ds-" || + particleName == "anti_He3" || + particleName == "anti_alpha" || + particleName == "anti_deuteron" || + particleName == "anti_lambda_c+" || + particleName == "anti_omega-" || + particleName == "anti_sigma_c+" || + particleName == "anti_sigma_c++" || + particleName == "anti_sigma+" || + particleName == "anti_sigma-" || + particleName == "anti_triton" || + particleName == "anti_xi_c+" || + particleName == "anti_xi-" || + particleName == "deuteron" || + particleName == "lambda_c+" || + particleName == "omega-" || + particleName == "sigma_c+" || + particleName == "sigma_c++" || + particleName == "sigma+" || + particleName == "sigma-" || + particleName == "tau+" || + particleName == "tau-" || + particleName == "triton" || + particleName == "xi_c+" || + particleName == "xi-" ) { + + ph->RegisterProcess(hmsc, particle); + ph->RegisterProcess(new G4hIonisation(), particle); + } + } + G4EmProcessOptions opt; + opt.SetVerbose(verbose); + opt.SetPolarAngleLimit(CLHEP::pi); + opt.SetApplyCuts(true); +}