Merge pull request #7577 from kreczko/pythia8-for-ttgenevent

Adding Pythia8 capability to TopDecaySubset

TopQuarkAnalysis/Configuration/test/TestForTTGenEvents.py

@@ -0,0 +1,88 @@
+import FWCore.ParameterSet.Config as cms
+import FWCore.ParameterSet.VarParsing as VarParsing
+import sys
+import os
+
+process = cms.Process("Analyzer")
+## enabling unscheduled mode for modules
+process.options = cms.untracked.PSet(
+    wantSummary = cms.untracked.bool(True),
+    allowUnscheduled = cms.untracked.bool(True),
+)
+options = VarParsing.VarParsing ('standard')
+options.register('runOnPythia8', True, 
+                 VarParsing.VarParsing.multiplicity.singleton, 
+                 VarParsing.VarParsing.varType.bool, 
+                 "decide to run on Pythia8 or Pythia6")
+options.register('useRun2Code', True, 
+                 VarParsing.VarParsing.multiplicity.singleton, 
+                 VarParsing.VarParsing.varType.bool, 
+                 "decide to use Run1 or Run2 code")
+
+# Get and parse the command line arguments
+if( hasattr(sys, "argv") ):
+    for args in sys.argv :
+        arg = args.split(',')
+        for val in arg:
+            val = val.split('=')
+            if(len(val)==2):
+                setattr(options,val[0], val[1])
+
+if options.runOnPythia8:
+    print 'Running on Pythia8'
+else:
+    print 'Running on Pythia6'
+
+## configure message logger
+process.load("FWCore.MessageLogger.MessageLogger_cfi")
+process.MessageLogger.cerr.threshold = 'INFO'
+process.MessageLogger.cerr.FwkReport.reportEvery = 100
+
+## define input
+if options.runOnPythia8:
+    process.source = cms.Source("PoolSource",
+       fileNames = cms.untracked.vstring(    
+         '/store/mc/Spring14dr/TTbarH_M-125_13TeV_amcatnlo-pythia8-tauola/AODSIM/PU20bx25_POSTLS170_V5-v1/00000/0E3D08A9-C610-E411-A862-0025B3E0657E.root',
+       ),
+       skipEvents = cms.untracked.uint32(0)
+    )
+else:
+    process.source = cms.Source("PoolSource",
+       fileNames = cms.untracked.vstring(
+        '/store/mc/Spring14dr/TTbarH_HToBB_M-125_13TeV_pythia6/AODSIM/PU20bx25_POSTLS170_V5-v1/00000/1CAB7E58-0BD0-E311-B688-00266CFFBC3C.root',
+       ),
+       skipEvents = cms.untracked.uint32(0)
+    )
+
+## define maximal number of events to loop over
+process.maxEvents = cms.untracked.PSet(
+    input = cms.untracked.int32(1000)
+)
+
+####################################################################
+genParticleCollection = 'genParticles'
+
+process.load("TopQuarkAnalysis.TopEventProducers.sequences.ttGenEvent_cff")
+process.initSubset.src = genParticleCollection
+process.decaySubset.src = genParticleCollection
+process.decaySubset.fillMode = "kME" # Status3, use kStable for Status2
+runMode = "Run2"
+if options.useRun2Code:
+    runMode = "Run2"
+else:
+    runMode = "Run1"
+process.decaySubset.runMode = runMode
+process.ttGenEventSequence = cms.Sequence(process.makeGenEvt)
+
+## module to store raw output from the processed modules into the ROOT file
+output_file = 'output.root'
+if options.runOnPythia8:
+    output_file = 'Pythia8_' + runMode + '.root'
+else:
+    output_file = 'Pythia6_' + runMode + '.root'
+
+process.out = cms.OutputModule("PoolOutputModule",
+    fileName = cms.untracked.string(output_file),
+    outputCommands = cms.untracked.vstring('drop *', 'keep *_*_*_Analyzer')
+)
+process.outpath = cms.EndPath(process.out)

TopQuarkAnalysis/TopEventProducers/interface/TopDecaySubset.h

@@ -11,6 +11,8 @@
 #include "FWCore/ParameterSet/interface/ParameterSet.h"
 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
 
+#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
+
 /**
    \class   TopDecaySubset TopDecaySubset.h "TopQuarkAnalysis/TopEventProducers/interface/TopDecaySubset.h"
 
@@ -31,7 +33,9 @@ class TopDecaySubset : public edm::EDProducer {
   /// of gen particles
   enum  FillMode {kStable, kME};
   /// classification of potential shower types
-  enum ShowerModel{kStart=-1, kNone, kPythia, kHerwig};
+  enum ShowerModel{kStart=-1, kNone, kPythia, kHerwig, kPythia8, kSherpa};
+  /// supported modes to run the code
+  enum RunMode {kRun1, kRun2};
 
   /// default constructor
   explicit TopDecaySubset(const edm::ParameterSet& cfg);
@@ -43,12 +47,34 @@ class TopDecaySubset : public edm::EDProducer {
  private:
   /// find top quarks in list of input particles
   std::vector<const reco::GenParticle*> findTops(const reco::GenParticleCollection& parts);
+  /// find primal top quarks (top quarks from the hard interaction)
+  /// for Pythia6 this is identical to findDecayingTops
+  std::vector<const reco::GenParticle*> findPrimalTops(const reco::GenParticleCollection& parts);
+  /// find decaying top quarks (quarks that decay to qW)
+  /// for Pythia6 this is identical to findPrimalTops
+  std::vector<const reco::GenParticle*> findDecayingTops(const reco::GenParticleCollection& parts);
+  /// find W bosons that come from top quark decays
+  /// for Pythia6 this is identical to findDecayingW
+  const reco::GenParticle* findPrimalW(const reco::GenParticle* top);
+  /// find W bosons that come from top quark decays and decay themselves (end of the MC chain)
+  /// for Pythia6 this is identical to findPrimalW
+//  const reco::GenParticle* findDecayingW(const reco::GenParticle* top);
+  /// find the last particle in a (potentially) long chain of state transitions
+  /// e.g. top[status==22]-> top[status==44 -> top[status==44] ->
+  /// top[status==44] -> top[status==62]
+  /// this function would pick the top with status 62
+  const reco::GenParticle* findLastParticleInChain(const reco::GenParticle* p);
   /// check the decay chain for the used shower model
   ShowerModel checkShowerModel(const std::vector<const reco::GenParticle*>& tops) const;
+  ShowerModel checkShowerModel(edm::Event& event);
   /// check whether W bosons are contained in the original gen particle listing
   void checkWBosons(std::vector<const reco::GenParticle*>& tops) const;
   /// fill output vector for full decay chain
   void fillListing(const std::vector<const reco::GenParticle*>& tops, reco::GenParticleCollection& target);
+  /// fill output vector for full decay chain
+  void fillListing(const std::vector<const reco::GenParticle*>& primalTops,
+				   const std::vector<const reco::GenParticle*>& decayingTops,
+				   reco::GenParticleCollection& target);
 
   /// clear references
   void clearReferences();
@@ -62,17 +88,22 @@ class TopDecaySubset : public edm::EDProducer {
   void addDaughters(int& idx, const reco::GenParticle::const_iterator part, reco::GenParticleCollection& target, bool recursive=true);
   /// fill vector including all radiations from quarks originating from W/top
   void addRadiation(int& idx, const reco::GenParticle::const_iterator part, reco::GenParticleCollection& target);
+  void addRadiation(int& idx, const reco::GenParticle* part, reco::GenParticleCollection& target);
 
  private:
   /// input tag for the genParticle source
   edm::EDGetTokenT<reco::GenParticleCollection> srcToken_;
+  /// input tag for the genEventInfo source
+  edm::EDGetTokenT<GenEventInfoProduct> genEventInfo_srcToken_;
   /// add radiation or not?
   bool addRadiation_;
   /// print the whole list of input particles or not?
   /// mode of decaySubset creation
   FillMode fillMode_;
   /// parton shower mode (filled in checkShowerModel)
   ShowerModel showerModel_;
+  /// run mode (Run1 || Run2)
+  RunMode runMode_;
 
   /// index in new evt listing of parts with daughters;
   /// has to be set to -1 in produce to deliver consistent

TopQuarkAnalysis/TopEventProducers/python/producers/TopDecaySubset_cfi.py

@@ -13,6 +13,10 @@
   ## 'kME'     : status 3 particles (from matrix element, before
   ##             parton shower) are saved (as status 3 particles)
   fillMode = cms.string("kStable"),
+  ## define run mode. The following run modes are available:
+  ## 'Run1' : will run the same code as during Run1
+  ## 'Run2' : will run new code that can deal with Pythia8
+  runMode = cms.string("Run1"),
   ## choose whether to save additionally radiated gluons in the
   ## decay chain or not
   addRadiation = cms.bool(True)