added options needed for neutron background study, tuned exotica phys…

…ics list

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)

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)

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),

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<double>("RusRoElectronEnergyLimit")*MeV;
+  double energyLim = 
+    theParSet.getParameter<double>("RusRoElectronEnergyLimit")*MeV;
   if(energyLim > 0.0) {
     rrfact[0] = theParSet.getParameter<double>("RusRoEcalElectron");
     rrfact[1] = theParSet.getParameter<double>("RusRoHcalElectron");
@@ -166,4 +171,11 @@ void ParametrisedEMPhysics::ConstructProcess() {
       ph->RegisterProcess(theElectronLimiter, G4Positron::Positron());
     }
   }
+  // enable fluorescence
+  bool fluo = theParSet.getParameter<bool>("FlagFluo");
+  if(fluo) {
+    G4VAtomDeexcitation* de = new G4UAtomicDeexcitation();
+    G4LossTableManager::Instance()->SetAtomDeexcitation(de);
+    de->SetFluo(true);
+  }
 }

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)

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

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'),
         )

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));    

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<edm::FileInPath>("particlesDef");
   particleDefFilePath = fp.fullPath();
   edm::LogInfo("CustomPhysics")<<"Path for custom particle definition file: "
 			       <<particleDefFilePath;
-  myHelper = 0;
-
- }
+  myHelper = 0;  
+}
 
 CustomPhysicsList::~CustomPhysicsList() {
   delete myHelper;
@@ -47,11 +47,10 @@ void CustomPhysicsList::ConstructProcess() {
 }
 
 void CustomPhysicsList::addCustomPhysics(){
-  LogDebug("CustomPhysics") << " CustomPhysics: adding CustomPhysics processes";
+  LogDebug("CustomPhysics") << " CustomPhysicsList: adding CustomPhysics processes";
   aParticleIterator->reset();
 
-  while((*aParticleIterator)())    {
-    int i = 0;
+  while((*aParticleIterator)()) {
     G4ParticleDefinition* particle = aParticleIterator->value();
     CustomParticle* cp = dynamic_cast<CustomParticle*>(particle);
     if(CustomParticleFactory::isCustomParticle(particle)) {
@@ -60,14 +59,19 @@ void CustomPhysicsList::addCustomPhysics(){
 				<< " is Custom. Mass is "
 				<<particle->GetPDGMass()/GeV  <<" GeV.";
       if(cp->GetCloud()!=0) {
-	LogDebug("CustomPhysics")<<"Cloud mass is "
-				 <<cp->GetCloud()->GetPDGMass()/GeV
-				 <<" GeV. Spectator mass is "
-				 <<static_cast<CustomParticle*>(particle)->GetSpectator()->GetPDGMass()/GeV
-				 <<" GeV.";
+	LogDebug("CustomPhysics")
+	  <<"Cloud mass is "
+	  <<cp->GetCloud()->GetPDGMass()/GeV
+	  <<" GeV. Spectator mass is "
+	  <<static_cast<CustomParticle*>(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,53 +80,48 @@ 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";
     }
   }
 }
 
 
-void CustomPhysicsList::setupRHadronPhycis(G4ParticleDefinition* particle){
-
+void CustomPhysicsList::setupRHadronPhycis(G4ParticleDefinition* particle)
+{
   //    LogDebug("CustomPhysics")<<"Configuring rHadron: "
   //	<<cp->
 
   CustomParticle* cp = dynamic_cast<CustomParticle*>(particle);
-  if(cp->GetCloud()!=0) 
-    LogDebug("CustomPhysics")<<"Cloud mass is "
-			     <<cp->GetCloud()->GetPDGMass()/GeV
-			     <<" GeV. Spectator mass is "
-			     <<static_cast<CustomParticle*>(particle)->GetSpectator()->GetPDGMass()/GeV
-			     <<" GeV.";
-
+  if(cp->GetCloud()!=0) {
+    LogDebug("CustomPhysics")
+      <<"Cloud mass is "
+      <<cp->GetCloud()->GetPDGMass()/GeV
+      <<" GeV. Spectator mass is "
+      <<static_cast<CustomParticle*>(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<CustomParticle*>(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";
 }

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/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;
 };