Merge pull request #5146 from civanch/tau-decay-bug-exotica

Backport of tau-decay bug fix and exotica tuning

SimG4Core/Application/python/g4SimHits_cfi.py

@@ -136,6 +136,7 @@
         ApplyPhiCuts = cms.bool(False),
         MinPhiCut = cms.double(-3.14159265359), ## (radians)
         MaxPhiCut = cms.double(3.14159265359), ## according to CMS conventions
+        ApplyLumiMonitorCuts = cms.bool(False), ## primary for lumi monitors
         Verbosity = cms.untracked.int32(0)
     ),
     RunAction = cms.PSet(
@@ -296,6 +297,7 @@
         ParametrizeLast   = cms.untracked.bool(False)
     ),
     HFShowerLibrary = cms.PSet(
+      #  FileName        = cms.FileInPath('SimG4CMS/Calo/data/HFShowerLibrary_oldpmt_noatt_eta4_16en.root'),
         FileName        = cms.FileInPath('SimG4CMS/Calo/data/hfshowerlibrary_lhep_140_edm.root'),
         BackProbability = cms.double(0.2),
         TreeEMID        = cms.string('emParticles'),

SimG4Core/CustomPhysics/python/Exotica_HSCP_SIM_cfi.py

@@ -2,68 +2,47 @@
 
 def customise(process):
 
-	FLAVOR = process.generator.hscpFlavor.value()
-	MASS_POINT = process.generator.massPoint.value()
-	SLHA_FILE = process.generator.slhaFile.value()
-	PROCESS_FILE = process.generator.processFile.value()
-	PARTICLE_FILE = process.generator.particleFile.value()
-	USE_REGGE = process.generator.useregge.value()
+    FLAVOR = process.generator.hscpFlavor.value()
+    MASS_POINT = process.generator.massPoint.value()
+    SLHA_FILE = process.generator.slhaFile.value()
+    PROCESS_FILE = process.generator.processFile.value()
+    PARTICLE_FILE = process.generator.particleFile.value()
+    USE_REGGE = process.generator.useregge.value()
 
-	process.load("SimG4Core.CustomPhysics.CustomPhysics_cfi")
-	process.customPhysicsSetup.particlesDef = PARTICLE_FILE
-	process.customPhysicsSetup.reggeModel = USE_REGGE
-	process.g4SimHits.Watchers = cms.VPSet (
-		cms.PSet(
-		type = cms.string('RHStopTracer'),
+    process.load("SimG4Core.CustomPhysics.CustomPhysics_cfi")
+    process.customPhysicsSetup.particlesDef = PARTICLE_FILE
+    process.customPhysicsSetup.reggeModel = USE_REGGE
+
+    if hasattr(process,'g4SimHits'):
+        # defined watches
+        process.g4SimHits.Watchers = cms.VPSet (
+	    cms.PSet(
+	        type = cms.string('RHStopTracer'),
 		RHStopTracer = cms.PSet(
 		verbose = cms.untracked.bool (False),
 		traceParticle = cms.string ("(~|tau1).*"),
 		stopRegularParticles = cms.untracked.bool (False)
 		)        
-		)
-		)
-
-
+	    )
+	)
+	# defined custom Physics List
+	process.g4SimHits.Physics.type = cms.string('SimG4Core/Physics/CustomPhysics')
+	# add verbosity
+	process.g4SimHits.Physics.Verbosity = cms.untracked.int32(0)
+	process.g4SimHits.G4Commands = cms.vstring('')
+	# check flavor of exotics and choose exotica Physics List
 	if FLAVOR=="gluino" or FLAVOR=="stop":
-		process.customPhysicsSetup.processesDef = PROCESS_FILE
-		process.g4SimHits.Physics = cms.PSet(
-        process.customPhysicsSetup,
-        DummyEMPhysics = cms.bool(True),
-        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),
-        FlagMuNucl  = cms.bool(False),
-        FlagFluo    = cms.bool(False),
-        GFlash = cms.PSet(
-        GflashHistogram = cms.bool(False),
-        GflashEMShowerModel = cms.bool(False),
-        GflashHadronPhysics = cms.string('QGSP_BERT_EMV'),
-        GflashHadronShowerModel = cms.bool(False)
-        )
-        )
-		process.g4SimHits.G4Commands = cms.vstring('/tracking/verbose 1')
-
+          process.customPhysicsSetup.processesDef = PROCESS_FILE
+          process.g4SimHits.Physics.ExoticaPhysicsSS = cms.untracked.bool(False)
 	elif FLAVOR =="stau":
-		process.g4SimHits.Physics = cms.PSet(
-        process.customPhysicsSetup,
-        DummyEMPhysics = cms.bool(True),
-        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'),
-        )
-		process.g4SimHits.G4Commands = cms.vstring('/tracking/verbose 1')
-
-	else:
-		print "Wrong flavor %s. Only accepted are gluino, stau, stop." % FLAVOR
+          process.g4SimHits.Physics.ExoticaPhysicsSS = cms.untracked.bool(False)
+        else:
+          print "Wrong flavor %s. Only accepted are gluino, stau, stop." % FLAVOR
+        # add custom options
+        process.g4SimHits.Physics = cms.PSet(
+            process.g4SimHits.Physics, #keep all default value and add others
+            process.customPhysicsSetup,
+        )	
 
 	return process
 

SimG4Core/CustomPhysics/src/CustomPhysics.cc

@@ -23,13 +23,9 @@ CustomPhysics::CustomPhysics(G4LogicalVolumeToDDLogicalPartMap& map,
   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 ssPhys  = p.getUntrackedParameter<bool>("ExoticaPhysicsSS",false);
   edm::LogInfo("PhysicsList") << "You are using the simulation engine: "
-			      << "QQGSP_FTFP_BERT_EML with Flags for EM Physics "
-			      << emPhys << " and for Hadronic Physics "
-			      << hadPhys << "\n";
+			      << "QGSP_FTFP_BERT_EML for regular particles";
 
   // EM Physics
   RegisterPhysics(new CMSEmStandardPhysics(ver));

SimG4Core/CustomPhysics/src/CustomPhysicsList.cc

@@ -48,32 +48,34 @@ void CustomPhysicsList::ConstructProcess() {
 }
 
 void CustomPhysicsList::addCustomPhysics(){
-  LogDebug("CustomPhysics") << " CustomPhysicsList: adding CustomPhysics processes";
+
+  edm::LogInfo("CustomPhysics") << " CustomPhysicsList: adding CustomPhysics processes "
+				<< "for the list of particles: \n";
   aParticleIterator->reset();
 
   while((*aParticleIterator)()) {
     G4ParticleDefinition* particle = aParticleIterator->value();
-    CustomParticle* cp = dynamic_cast<CustomParticle*>(particle);
     if(CustomParticleFactory::isCustomParticle(particle)) {
-      LogDebug("CustomPhysics") << particle->GetParticleName()
-				<<", "<<particle->GetPDGEncoding()
-				<< " is Custom. Mass is "
-				<<particle->GetPDGMass()/GeV  <<" GeV.";
+      CustomParticle* cp = dynamic_cast<CustomParticle*>(particle);
+      edm::LogInfo("CustomPhysics") << particle->GetParticleName()
+				    <<"  PDGcode= "<<particle->GetPDGEncoding()
+				    << "  Mass= "
+				    <<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.";
+	edm::LogInfo("CustomPhysics") << particle->GetParticleName()
+				      <<" CloudMass= "
+				      <<cp->GetCloud()->GetPDGMass()/GeV
+				      <<" GeV; SpectatorMass= "
+				      << cp->GetSpectator()->GetPDGMass()/GeV
+				      <<" GeV.";
       }
       G4ProcessManager* pmanager = particle->GetProcessManager();
       if(pmanager) {
-	if(particle->GetPDGCharge()/eplus != 0) {
+	if(particle->GetPDGCharge() != 0.0) {
 	  pmanager->AddProcess(new G4hMultipleScattering,-1, 1, 1);
 	  pmanager->AddProcess(new G4hIonisation,        -1, 2, 2);
 	}
-	if(cp!=0) {
+	if(cp != 0) {
 	  if(particle->GetParticleType()=="rhadron" || 
 	     particle->GetParticleType()=="mesonino" || 
 	     particle->GetParticleType() == "sbaryon"){

SimG4Core/CustomPhysics/src/CustomPhysicsListSS.cc

@@ -48,24 +48,24 @@ void CustomPhysicsListSS::ConstructProcess() {
 }
 
 void CustomPhysicsListSS::addCustomPhysics(){
-  LogDebug("CustomPhysics") << " CustomPhysicsListSS: adding CustomPhysics processes";
+  edm::LogInfo("CustomPhysics") << " CustomPhysicsListSS: adding CustomPhysics processes";
   aParticleIterator->reset();
 
   while((*aParticleIterator)()) {
     G4ParticleDefinition* particle = aParticleIterator->value();
-    CustomParticle* cp = dynamic_cast<CustomParticle*>(particle);
     if(CustomParticleFactory::isCustomParticle(particle)) {
-      LogDebug("CustomPhysics") << particle->GetParticleName()
-				<<", "<<particle->GetPDGEncoding()
-				<< " is Custom. Mass is "
-				<<particle->GetPDGMass()/GeV  <<" GeV.";
+      CustomParticle* cp = dynamic_cast<CustomParticle*>(particle);
+      edm::LogInfo("CustomPhysics") << particle->GetParticleName()
+				    <<" PDGcode= "<<particle->GetPDGEncoding()
+				    << " Mass= "
+				    <<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.";
+	edm::LogInfo("CustomPhysics") <<  particle->GetParticleName()
+				      << " CloudMass= "
+				      <<cp->GetCloud()->GetPDGMass()/GeV
+				      <<" GeV; SpectatorMass= "
+				      <<cp->GetSpectator()->GetPDGMass()/GeV
+				      <<" GeV.";
       }
       G4ProcessManager* pmanager = particle->GetProcessManager();
       if(pmanager) {

SimG4Core/Generators/interface/Generator.h

@@ -14,6 +14,7 @@
 
 class G4Event;
 class G4PrimaryParticle;
+class LumiMonitorFilter;
 
 class Generator
 {
@@ -53,6 +54,7 @@ class Generator
   double theEtaCutForHector; 
   double theDecLenCut;
   int verbose;
+  LumiMonitorFilter* fLumiFilter; 
   HepMC::GenEvent*  evt_;
   math::XYZTLorentzVector* vtx_;
   double weight_;    

SimG4Core/Generators/interface/LumiMonitorFilter.h

@@ -0,0 +1,19 @@
+#ifndef SimG4Core_LumiMonitorFilter_H
+#define SimG4Core_LumiMonitorFilter_H
+
+#include "HepMC/GenParticle.h"
+
+class LumiMonitorFilter
+{
+public:
+  LumiMonitorFilter();
+  virtual ~LumiMonitorFilter();
+
+  void Describe() const; 
+  bool isGoodForLumiMonitor(const HepMC::GenParticle*) const;
+
+private:
+
+};
+
+#endif

SimG4Core/Generators/src/Generator.cc

@@ -1,5 +1,6 @@
 #include "SimG4Core/Generators/interface/Generator.h"
 #include "SimG4Core/Generators/interface/HepMCParticle.h"
+#include "SimG4Core/Generators/interface/LumiMonitorFilter.h"
 
 #include "SimG4Core/Notification/interface/SimG4Exception.h"
 
@@ -31,6 +32,7 @@ Generator::Generator(const ParameterSet & p) :
   theMaxPCut(p.getParameter<double>("MaxPCut")),   
   theEtaCutForHector(p.getParameter<double>("EtaCutForHector")),
   verbose(p.getUntrackedParameter<int>("Verbosity",0)),
+  fLumiFilter(0),
   evt_(0),
   vtx_(0),
   weight_(0),
@@ -40,6 +42,9 @@ Generator::Generator(const ParameterSet & p) :
   pdgFilterSel(false), 
   fPDGFilter(false) 
 {
+  bool lumi = p.getParameter<bool>("ApplyLumiMonitorCuts");
+  if(lumi) { fLumiFilter = new LumiMonitorFilter(); }
+
   double theRDecLenCut = p.getParameter<double>("RDecLenCut")*cm;
   theDecRCut2 = theRDecLenCut*theRDecLenCut;
 
@@ -60,8 +65,8 @@ Generator::Generator(const ParameterSet & p) :
 	  edm::LogWarning("SimG4CoreGenerator") 
 	    << " *** Selecting only PDG ID = " << pdgFilter[ii];
 	} else {
-	  edm::LogWarning("SimG4CoreGenerator") << " *** Filtering out PDG ID = " 
-						<< pdgFilter[ii];
+	  edm::LogWarning("SimG4CoreGenerator") 
+	    << " *** Filtering out PDG ID = " << pdgFilter[ii];
 	}
       }
     }
@@ -82,11 +87,15 @@ Generator::Generator(const ParameterSet & p) :
     << " cm;  Z_hector = " << Z_hector << " cm\n"
     << "ApplyPCuts: " << fPCuts << "  ApplyPtransCut: " << fPtransCut
     << "  ApplyEtaCuts: " << fEtaCuts 
-    << "  ApplyPhiCuts: " << fPhiCuts;
+    << "  ApplyPhiCuts: " << fPhiCuts
+    << "  ApplyLumiMonitorCuts: " << lumi;
+  if(lumi) { fLumiFilter->Describe(); }
 }
 
 Generator::~Generator() 
-{}
+{
+  delete fLumiFilter;
+}
 
 void Generator::HepMC2G4(const HepMC::GenEvent * evt, G4Event * g4evt)
 {
@@ -302,6 +311,8 @@ void Generator::HepMC2G4(const HepMC::GenEvent * evt, G4Event * g4evt)
 	      continue;
 	    }
 	  }
+	  if(fLumiFilter && !fLumiFilter->isGoodForLumiMonitor(*pitr)) 
+	    { continue; }
 	  toBeAdded = true;
 	  if ( verbose > 2 ) LogDebug("SimG4CoreGenerator") 
 	    << "GenParticle barcode = " << (*pitr)->barcode() 
@@ -436,8 +447,8 @@ void Generator::particleAssignDaughters( G4PrimaryParticle* g4p,
         double dd = std::sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
         particleAssignDaughters(g4daught,*vpdec,dd);
       }
-
-    (*vpdec)->set_status(1000+(*vpdec)->status()); 
+    // this line is the source of the "tau-bug" in 7_2_0_pre5
+    //(*vpdec)->set_status(1000+(*vpdec)->status()); 
     g4p->SetDaughter(g4daught);
 
     if ( verbose > 1 ) g4daught->Print();

SimG4Core/Generators/src/LumiMonitorFilter.cc

@@ -0,0 +1,24 @@
+#include "SimG4Core/Generators/interface/LumiMonitorFilter.h"
+
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+
+using namespace edm;
+//using std::cout;
+//using std::endl;
+
+LumiMonitorFilter::LumiMonitorFilter() 
+{} 
+
+LumiMonitorFilter::~LumiMonitorFilter() 
+{}
+
+void LumiMonitorFilter::Describe() const
+{
+  edm::LogInfo("LumiMonitorFilter") 
+    << " is active ";
+}
+
+bool LumiMonitorFilter::isGoodForLumiMonitor(const HepMC::GenParticle*) const
+{
+  return true;
+}

SimG4Core/Physics/src/G4ProcessTypeEnumerator.cc

@@ -1,7 +1,7 @@
 #include "SimG4Core/Physics/interface/G4ProcessTypeEnumerator.h"
 
-static const int nprocesses = 47;
-static const std::string g4processes[nprocesses] = { 
+static const int nprocesses = 49;
+static const std::string g4processes[nprocesses] = { "Primary", 
   "Transportation",  "CoupleTrans",  "CoulombScat",  "Ionisation",  "Brems",
   "PairProdCharged",  "Annih",  "AnnihToMuMu",  "AnnihToHad",  "NuclearStopp",
   "Msc",  "Rayleigh",  "PhotoElectric",  "Compton",  "Conv",
@@ -10,9 +10,10 @@ static const std::string g4processes[nprocesses] = {
   "DNAElastic",  "DNAExcit",  "DNAIonisation",  "DNAVibExcit",  "DNAAttachment",
   "DNAChargeDec",  "DNAChargeInc",  "HadElastic",  "HadInelastic",  "HadCapture",
   "HadFission",  "HadAtRest",  "HadCEX",  "Decay",  "DecayWSpin",
-  "DecayPiWSpin",  "DecayRadio",  "DecayUnKnown",  "DecayExt",  "StepLimiter",
+  "DecayPiWSpin",  "DecayRadio",  "DecayUnKnown",  "DecayExt", "GFlash", "StepLimiter",
   "UsrSpecCuts", "NeutronKiller"}; 
 static const int g4subtype[nprocesses] = { 
+  0,   // Primary generator
   91,  // Transportation
   92,  // CoupleTrans
   1,   // CoulombScat
@@ -57,6 +58,7 @@ static const int g4subtype[nprocesses] = {
   210, // DecayRadio
   211, // DecayUnKnown
   231, // DecayExt
+  301, // GFlash
   401, // StepLimiter
   402,
   403  // NeutronKiller