Merge pull request #4054 from civanch/neutron_backgound

Neutron backgound: backport from 7_1_X

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)

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)

SimG4Core/Application/python/g4SimHits_cfi.py

@@ -77,13 +77,15 @@
         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
         #
         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),

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<bool>("FlagFluo");
+  if(fluo) {
+    G4VAtomDeexcitation* de = new G4UAtomicDeexcitation();
+    G4LossTableManager::Instance()->SetAtomDeexcitation(de);
+    de->SetFluo(true);
+  }
 }

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 <string>
@@ -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_;    
 };
 

SimG4Core/Physics/src/DDG4ProductionCuts.cc

@@ -7,9 +7,11 @@
 
 #include <algorithm>
 
-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<bool>("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

SimG4Core/Physics/src/PhysicsList.cc

@@ -3,14 +3,15 @@
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 
 #include "G4LossTableManager.hh"
+#include "G4SystemOfUnits.hh"
 
 PhysicsList::PhysicsList(G4LogicalVolumeToDDLogicalPartMap & map,
 			 const HepPDT::ParticleDataTable * table_,
 			 sim::FieldBuilder *fieldBuilder_,
 			 const edm::ParameterSet & p) 
   : G4VModularPhysicsList(), m_pPhysics(p),  prodCuts(0) {
   m_Verbosity = m_pPhysics.getUntrackedParameter<int>("Verbosity",0);
-  prodCuts = new DDG4ProductionCuts(map, m_Verbosity);	
+  prodCuts = new DDG4ProductionCuts(map, m_Verbosity, m_pPhysics);	
 }
 
 PhysicsList::~PhysicsList() {

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
+
+
+
+
+
+

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<int>("Verbosity",0);
+  bool emPhys  = p.getUntrackedParameter<bool>("EMPhysics",true);
+  bool hadPhys = p.getUntrackedParameter<bool>("HadPhysics",true);
+  bool tracking= p.getParameter<bool>("TrackingCut");
+  //bool munucl  = p.getParameter<bool>("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));
+}
+

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
+
+
+

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<int>("Verbosity",0);
+  bool emPhys  = p.getUntrackedParameter<bool>("EMPhysics",true);
+  bool hadPhys = p.getUntrackedParameter<bool>("HadPhysics",true);
+  bool tracking= p.getParameter<bool>("TrackingCut");
+  //bool munucl  = p.getParameter<bool>("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));
+}
+