Sub-event aware genJet and parton matching for HI miniAOD

Configuration/ProcessModifiers/python/genJetSubEvent_cff.py

@@ -0,0 +1,4 @@
+import FWCore.ParameterSet.Config as cms
+
+genJetSubEvent =  cms.Modifier()
+

Configuration/PyReleaseValidation/python/relval_standard.py

@@ -700,14 +700,14 @@
 workflows[148] = ['',['HydjetQ_MinBias_XeXe_5442GeV_2017','DIGIHI2017','RECOHI2017','HARVESTHI2017']]
 ### Run II cond., 2018
 workflows[150] = ['',['HydjetQ_B12_5020GeV_2018','DIGIHI2018','RECOHI2018','HARVESTHI2018']]
-workflows[158] = ['',['HydjetQ_B12_5020GeV_2018_ppReco','DIGIHI2018PPRECO','RECOHI2018PPRECO','ALCARECOHI2018PPRECO','HARVESTHI2018PPRECO']]
-workflows[158.01] = ['',['HydjetQ_reminiaodPbPb2018_INPUT','REMINIAODHI2018PPRECO','HARVESTHI2018PPRECOMINIAOD']]
+workflows[158] = ['',['HydjetQ_B12_5020GeV_2018_ppReco','DIGIHI2018PPRECO','RECOHI2018PPRECOMB','ALCARECOHI2018PPRECO','HARVESTHI2018PPRECO']]
+workflows[158.01] = ['',['HydjetQ_reminiaodPbPb2018_INPUT','REMINIAODHI2018PPRECOMB','HARVESTHI2018PPRECOMINIAOD']]
 workflows[158.1] = ['',['QCD_Pt_80_120_13_HI','DIGIHI2018PPRECO','RECOHI2018PPRECO','HARVESTHI2018PPRECO']]
 workflows[158.2] = ['',['PhotonJets_Pt_10_13_HI','DIGIHI2018PPRECO','RECOHI2018PPRECO','HARVESTHI2018PPRECO']]
 workflows[158.3] = ['',['ZEEMM_13_HI','DIGIHI2018PPRECO','RECOHI2018PPRECO','HARVESTHI2018PPRECO']]
 # 158.0 with multiple concurrent lumi sections
 workflows[158.181] = ['',['HydjetQ_B12_5020GeV_2018_ppReco_ml','DIGIHI2018PPRECOml','RECOHI2018PPRECOml','ALCARECOHI2018PPRECOml','HARVESTHI2018PPRECO']]
-workflows[159] = ['',['HydjetQ_B12_5020GeV_2021_ppReco','DIGIHI2021PPRECO','RECOHI2021PPRECO','ALCARECOHI2021PPRECO','HARVESTHI2021PPRECO']]
+workflows[159] = ['',['HydjetQ_B12_5020GeV_2021_ppReco','DIGIHI2021PPRECO','RECOHI2021PPRECOMB','ALCARECOHI2021PPRECO','HARVESTHI2021PPRECO']]
 workflows[159.1] = ['',['QCD_Pt_80_120_14_HI_2021','DIGIHI2021PPRECO','RECOHI2021PPRECO','HARVESTHI2021PPRECO']]
 workflows[159.2] = ['',['PhotonJets_Pt_10_14_HI_2021','DIGIHI2021PPRECO','RECOHI2021PPRECO','HARVESTHI2021PPRECO']]
 workflows[159.3] = ['',['ZMM_14_HI_2021','DIGIHI2021PPRECO','RECOHI2021PPRECO','HARVESTHI2021PPRECO']]

Configuration/PyReleaseValidation/python/relval_steps.py

@@ -2474,6 +2474,12 @@ def gen2021HiMix(fragment,howMuch):
                                                         '--datatier':'GEN-SIM-RECO,MINIAODSIM,DQMIO',
                                                         '--eventcontent':'RECOSIM,MINIAODSIM,DQM',
                                                         },step3Up2015Defaults])
+steps['RECOHI2021PPRECOMB']=merge([hiDefaults2021_ppReco,{'-s':'RAW2DIGI,L1Reco,RECO,EI,PAT,VALIDATION:@standardValidationNoHLT+@miniAODValidation,DQM:@standardDQMFakeHLT+@miniAODDQM',
+                                                          '--datatier':'GEN-SIM-RECO,MINIAODSIM,DQMIO',
+                                                          '--eventcontent':'RECOSIM,MINIAODSIM,DQM',
+                                                          '--era':'Run3_pp_on_PbPb',
+                                                          '--procModifiers':'genJetSubEvent',
+                                                        },step3Up2015Defaults])
 
 steps['ALCARECOHI2021PPRECO']=merge([hiDefaults2021_ppReco,{'-s':'ALCA:TkAlMinBias+SiStripCalMinBias',
                                                             '--datatier':'ALCARECO',
@@ -2484,13 +2490,26 @@ def gen2021HiMix(fragment,howMuch):
                                                         '--datatier':'GEN-SIM-RECO,MINIAODSIM,DQMIO,ALCARECO',
                                                         '--eventcontent':'RECOSIM,MINIAODSIM,DQM,ALCARECO',
                                                         },step3Up2015Defaults])
-steps['RECOHI2018PPRECOml']=merge([concurrentLumis,steps['RECOHI2018PPRECO']])
+steps['RECOHI2018PPRECOMB']=merge([hiDefaults2018_ppReco,{'-s':'RAW2DIGI,L1Reco,RECO,ALCA:SiStripCalZeroBias+SiPixelCalZeroBias,EI,PAT,VALIDATION:@standardValidationNoHLT+@miniAODValidation,DQM:@standardDQMFakeHLT+@miniAODDQM',
+                                                        '--datatier':'GEN-SIM-RECO,MINIAODSIM,DQMIO,ALCARECO',
+                                                        '--eventcontent':'RECOSIM,MINIAODSIM,DQM,ALCARECO',
+                                                          '--era':'Run2_2018_pp_on_AA',
+                                                          '--procModifiers':'genJetSubEvent',
+                               },step3Up2015Defaults])
+steps['RECOHI2018PPRECOml']=merge([concurrentLumis,steps['RECOHI2018PPRECOMB']])
 
 steps['REMINIAODHI2018PPRECO']=merge([{'-s':'PAT,VALIDATION:@miniAODValidation,DQM:@miniAODDQM',
                                                              '--datatier':'MINIAODSIM,DQMIO',
                                                              '--eventcontent':'MINIAODSIM,DQM',
                                                              '-n':100,
                                                          },hiDefaults2018_ppReco,step3Up2015Defaults])
+steps['REMINIAODHI2018PPRECOMB']=merge([{'-s':'PAT,VALIDATION:@miniAODValidation,DQM:@miniAODDQM',
+                                                             '--datatier':'MINIAODSIM,DQMIO',
+                                                             '--eventcontent':'MINIAODSIM,DQM',
+                                                             '-n':100,
+                                         '--era':'Run2_2018_pp_on_AA',
+                                         '--procModifiers':'genJetSubEvent',
+                                                         },hiDefaults2018_ppReco,step3Up2015Defaults])
 
 steps['ALCARECOHI2018PPRECO']=merge([hiDefaults2018_ppReco,{'-s':'ALCA:TkAlMinBias+SiStripCalMinBias',
                                                             '--datatier':'ALCARECO',

PhysicsTools/PatAlgos/plugins/PATGenJetSlimmer.cc

@@ -76,6 +76,7 @@ void pat::PATGenJetSlimmer::produce(edm::Event& iEvent, const edm::EventSetup& i
   auto mapping = std::make_unique<std::vector<int> >();
   mapping->reserve(src->size());
 
+  unsigned nm = 0;  // number of mapped jets
   unsigned nl = 0;  // number of loose jets
   for (View<reco::GenJet>::const_iterator it = src->begin(), ed = src->end(); it != ed; ++it) {
     bool selectedLoose = false;
@@ -93,7 +94,7 @@ void pat::PATGenJetSlimmer::produce(edm::Event& iEvent, const edm::EventSetup& i
     out->push_back(*it);
     reco::GenJet& jet = out->back();
 
-    mapping->push_back(it - src->begin());
+    mapping->push_back(nm++);
 
     if (clearDaughters_) {
       jet.clearDaughters();

PhysicsTools/PatAlgos/python/mcMatchLayer0/jetMatch_cfi.py

@@ -30,13 +30,9 @@
 from Configuration.Eras.Modifier_pp_on_AA_2018_cff import pp_on_AA_2018
 from Configuration.Eras.Modifier_pp_on_PbPb_run3_cff import pp_on_PbPb_run3
 (pp_on_AA_2018 | pp_on_PbPb_run3).toModify(patJetGenJetMatch,
-                                           matched = 'ak4HiCleanedGenJets',
                                            maxDeltaR = 0.4,
                                            resolveByMatchQuality = True,
                                            src = "akCs4PFJets",
                                        )
 
-(pp_on_AA_2018 | pp_on_PbPb_run3).toModify(patJetPartonMatch,
-                                           matched = "cleanedPartons",
-                                           src = "akCs4PFJets",
-                                       )
+(pp_on_AA_2018 | pp_on_PbPb_run3).toModify(patJetPartonMatch, src = "akCs4PFJets")

PhysicsTools/PatAlgos/python/producersHeavyIons/heavyIonJets_cff.py

@@ -1,29 +1,40 @@
 import FWCore.ParameterSet.Config as cms
 
+from RecoHI.HiJetAlgos.HiSignalParticleProducer_cfi import hiSignalGenParticles
 from RecoJets.Configuration.GenJetParticles_cff import genParticlesForJets
 from RecoHI.HiJetAlgos.HiGenCleaner_cff import hiPartons
 from RecoHI.HiJetAlgos.HiGenJets_cff import ak4HiGenJets
 from RecoHI.HiJetAlgos.HiGenCleaner_cff import heavyIonCleanedGenJets
+from RecoHI.HiJetAlgos.HiSignalGenJetProducer_cfi import hiSignalGenJets
 
 allPartons = cms.EDProducer(
     "PartonSelector",
-    src = cms.InputTag('genParticles'),
+    src = cms.InputTag('hiSignalGenParticles'),
     withLeptons = cms.bool(False),
     )
 
+from Configuration.ProcessModifiers.genJetSubEvent_cff import genJetSubEvent
+genJetSubEvent.toModify(allPartons,src = "genParticles")
+
 cleanedPartons = hiPartons.clone(
     src = 'allPartons',
     )
 
-ak4HiGenJetsCleaned = heavyIonCleanedGenJets.clone(src = "ak4HiGenJets")
+ak4HiSignalGenJets = hiSignalGenJets.clone(src = "ak4HiGenJets")
 
-cleanedGenJetsTask = cms.Task(
+hiGenJetsTask = cms.Task(
+    hiSignalGenParticles,
     genParticlesForJets,
-    cleanedPartons,
+    allPartons,
     ak4HiGenJets,
-    ak4HiGenJetsCleaned
+    ak4HiSignalGenJets
 )
 
+ak4HiGenJetsCleaned = heavyIonCleanedGenJets.clone(src = "ak4HiGenJets")
+hiCleanedGenJetsTask_ = hiGenJetsTask.copyAndExclude([hiSignalGenParticles,ak4HiSignalGenJets])
+hiCleanedGenJetsTask_.add(cleanedPartons,ak4HiGenJetsCleaned)
+genJetSubEvent.toReplaceWith(hiGenJetsTask,hiCleanedGenJetsTask_)
+
 from RecoHI.HiJetAlgos.HiRecoPFJets_cff import PFTowers, pfEmptyCollection, ak4PFJetsForFlow, hiPuRho, hiFJRhoFlowModulation, akCs4PFJets
 from RecoHI.HiTracking.highPurityGeneralTracks_cfi import highPurityGeneralTracks
 from RecoJets.JetAssociationProducers.ak5JTA_cff import *
@@ -62,5 +73,5 @@
 recoJetsHIpostAODTask = cms.Task(
     recoPFJetsHIpostAODTask,
     allPartons,
-    cleanedGenJetsTask,
+    hiGenJetsTask,
     )

PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py

@@ -588,6 +588,9 @@ def miniAOD_customizeMC(process):
     process.tauGenJetsBoosted.GenParticles = "prunedGenParticles"
     process.patJetPartons.particles = "genParticles"
     process.patJetPartonMatch.matched = "prunedGenParticles"
+    _hiGeneral.toModify(process.patJetPartonMatch, matched =  "hiSignalGenParticles")
+    from Configuration.ProcessModifiers.genJetSubEvent_cff import genJetSubEvent
+    genJetSubEvent.toModify(process.patJetPartonMatch, matched =  "cleanedPartons")
     process.patJetPartonMatch.mcStatus = [ 3, 23 ]
     process.patJetGenJetMatch.matched = "slimmedGenJets"
     (~_hiGeneral).toModify(process, patJetGenJetMatchAK8Puppi = dict(matched =  "slimmedGenJetsAK8"))
@@ -603,6 +606,7 @@ def miniAOD_customizeMC(process):
 
     from PhysicsTools.PatAlgos.producersHeavyIons.heavyIonJetSetup import removeJECsForMC
     _hiGeneral.toModify(process, removeJECsForMC)
+    _hiGeneral.toReplaceWith(task,task.copyAndExclude([process.slimmedGenJetsFlavourInfos]))
 
 
 def miniAOD_customizeOutput(out):

PhysicsTools/PatAlgos/python/slimming/slimmedGenJets_cfi.py

@@ -22,4 +22,7 @@
 
 from Configuration.Eras.Modifier_pp_on_AA_2018_cff import pp_on_AA_2018
 from Configuration.Eras.Modifier_pp_on_PbPb_run3_cff import pp_on_PbPb_run3
+(pp_on_AA_2018 | pp_on_PbPb_run3).toModify(slimmedGenJets, src = "ak4HiSignalGenJets")
 (pp_on_AA_2018 | pp_on_PbPb_run3).toModify(slimmedGenJetsAK8, cut = 'pt>9999', nLoose = 0)
+from Configuration.ProcessModifiers.genJetSubEvent_cff import genJetSubEvent
+genJetSubEvent.toModify(slimmedGenJets, src = "ak4HiGenJetsCleaned")

RecoHI/HiJetAlgos/plugins/HiSignalGenJetProducer.cc

@@ -0,0 +1,128 @@
+// -*- C++ -*-
+//
+// Package:    HiSignalGenJetProducer
+// Class:      HiSignalGenJetProducer
+//
+/**\class HiSignalGenJetProducer HiSignalGenJetProducer.cc yetkin/HiSignalGenJetProducer/src/HiSignalGenJetProducer.cc
+
+ Description: <one line class summary>
+
+ Implementation:
+     <Notes on implementation>
+*/
+//
+// Original Author:  Yetkin Yilmaz
+//         Created:  Tue Jul 21 04:26:01 EDT 2009
+//
+//
+
+// system include files
+#include <memory>
+#include <vector>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/EDProducer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+
+#include "DataFormats/Common/interface/View.h"
+#include "DataFormats/JetReco/interface/GenJetCollection.h"
+#include "DataFormats/GeometryVector/interface/VectorUtil.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
+
+#include "DataFormats/Math/interface/Point3D.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
+
+using namespace std;
+using namespace edm;
+
+//
+// class decleration
+//
+
+class HiSignalGenJetProducer : public edm::EDProducer {
+public:
+  explicit HiSignalGenJetProducer(const edm::ParameterSet&);
+  ~HiSignalGenJetProducer() override;
+
+private:
+  void produce(edm::Event&, const edm::EventSetup&) override;
+  // ----------member data ---------------------------
+
+  edm::EDGetTokenT<edm::View<reco::GenJet> > jetSrc_;
+};
+
+//
+// constants, enums and typedefs
+//
+
+//
+// static data member definitions
+//
+
+//
+// constructors and destructor
+//
+
+HiSignalGenJetProducer::HiSignalGenJetProducer(const edm::ParameterSet& iConfig)
+    : jetSrc_(consumes<edm::View<reco::GenJet> >(iConfig.getParameter<edm::InputTag>("src"))) {
+  std::string alias = (iConfig.getParameter<InputTag>("src")).label();
+  produces<reco::GenJetCollection>().setBranchAlias(alias);
+}
+
+HiSignalGenJetProducer::~HiSignalGenJetProducer() {
+  // do anything here that needs to be done at desctruction time
+  // (e.g. close files, deallocate resources etc.)
+}
+
+//
+// member functions
+//
+
+// ------------ method called to produce the data  ------------
+
+void HiSignalGenJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  using namespace edm;
+  using namespace reco;
+
+  auto jets = std::make_unique<GenJetCollection>();
+
+  edm::Handle<edm::View<GenJet> > genjets;
+  iEvent.getByToken(jetSrc_, genjets);
+
+  int jetsize = genjets->size();
+
+  vector<int> selection;
+  selection.reserve(jetsize);
+  for (int ijet = 0; ijet < jetsize; ++ijet) {
+    selection.push_back(-1);
+  }
+
+  vector<int> selectedIndices;
+  vector<int> removedIndices;
+
+  for (int ijet = 0; ijet < jetsize; ++ijet) {
+    const GenJet* jet1 = &((*genjets)[ijet]);
+
+    const GenParticle* gencon = jet1->getGenConstituent(0);
+
+    if (gencon == nullptr)
+      throw cms::Exception("GenConstituent", "GenJet is missing its constituents");
+    else if (gencon->collisionId() == 0) {
+      jets->push_back(*jet1);
+      selection[ijet] = 1;
+    } else {
+      selection[ijet] = 0;
+      removedIndices.push_back(ijet);
+    }
+  }
+
+  iEvent.put(std::move(jets));
+}
+
+DEFINE_FWK_MODULE(HiSignalGenJetProducer);