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,62 @@
+#include <memory>
+#include <vector>
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/global/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/HepMCCandidate/interface/GenParticle.h"
+
+class HiSignalGenJetProducer : public edm::global::EDProducer<> {
+public:
+  explicit HiSignalGenJetProducer(const edm::ParameterSet&);
+  ~HiSignalGenJetProducer() override = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions&);
+
+private:
+  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
+
+  edm::EDGetTokenT<edm::View<reco::GenJet> > jetSrc_;
+};
+
+HiSignalGenJetProducer::HiSignalGenJetProducer(const edm::ParameterSet& iConfig)
+    : jetSrc_(consumes<edm::View<reco::GenJet> >(iConfig.getParameter<edm::InputTag>("src"))) {
+  std::string alias = (iConfig.getParameter<edm::InputTag>("src")).label();
+  produces<reco::GenJetCollection>().setBranchAlias(alias);
+}
+
+void HiSignalGenJetProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
+  auto jets = std::make_unique<reco::GenJetCollection>();
+
+  edm::Handle<edm::View<reco::GenJet> > genjets;
+  iEvent.getByToken(jetSrc_, genjets);
+
+  for (const reco::GenJet& jet1 : *genjets) {
+    const reco::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);
+    }
+  }
+
+  iEvent.put(std::move(jets));
+}
+
+void HiSignalGenJetProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.setComment("Selects genJets from collision id = 0");
+  desc.add<edm::InputTag>("src", edm::InputTag("akHiGenJets"));
+  descriptions.add("HiSignalGenJetProducer", desc);
+}
+
+DEFINE_FWK_MODULE(HiSignalGenJetProducer);

RecoHI/HiJetAlgos/plugins/HiSignalParticleProducer.cc

@@ -0,0 +1,59 @@
+#include <memory>
+#include <vector>
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/global/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/HepMCCandidate/interface/GenParticle.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
+
+class HiSignalParticleProducer : public edm::global::EDProducer<> {
+public:
+  explicit HiSignalParticleProducer(const edm::ParameterSet&);
+  ~HiSignalParticleProducer() override = default;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions&);
+
+private:
+  void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
+
+  edm::EDGetTokenT<edm::View<reco::GenParticle> > genParticleSrc_;
+};
+
+HiSignalParticleProducer::HiSignalParticleProducer(const edm::ParameterSet& iConfig)
+    : genParticleSrc_(consumes<edm::View<reco::GenParticle> >(iConfig.getParameter<edm::InputTag>("src"))) {
+  std::string alias = (iConfig.getParameter<edm::InputTag>("src")).label();
+  produces<reco::GenParticleCollection>().setBranchAlias(alias);
+}
+
+void HiSignalParticleProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup&) const {
+  auto signalGenParticles = std::make_unique<reco::GenParticleCollection>();
+
+  edm::Handle<edm::View<reco::GenParticle> > genParticles;
+  iEvent.getByToken(genParticleSrc_, genParticles);
+
+  for (const reco::GenParticle& genParticle : *genParticles) {
+    if (genParticle.collisionId() == 0) {
+      signalGenParticles->push_back(genParticle);
+    }
+  }
+
+  iEvent.put(std::move(signalGenParticles));
+}
+
+void HiSignalParticleProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.setComment("Selects genParticles from collision id = 0");
+  desc.add<edm::InputTag>("src", edm::InputTag("genParticles"));
+  descriptions.add("HiSignalParticleProducer", desc);
+}
+
+DEFINE_FWK_MODULE(HiSignalParticleProducer);

RecoHI/HiJetAlgos/python/HiSignalGenJetProducer_cfi.py

@@ -0,0 +1,7 @@
+import FWCore.ParameterSet.Config as cms
+
+
+hiSignalGenJets = cms.EDProducer('HiSignalGenJetProducer',
+                                    src    = cms.InputTag('ak4HiGenJets')
+                                    )
+

RecoHI/HiJetAlgos/python/HiSignalParticleProducer_cfi.py

@@ -0,0 +1,7 @@
+import FWCore.ParameterSet.Config as cms
+
+
+hiSignalGenParticles = cms.EDProducer('HiSignalParticleProducer',
+                                    src    = cms.InputTag('genParticles')
+                                    )
+

RecoJets/JetProducers/plugins/SubEventGenJetProducer.cc

@@ -2,8 +2,6 @@
 #include "FWCore/Framework/interface/MakerMacros.h"
 #include "RecoJets/JetProducers/plugins/SubEventGenJetProducer.h"
 
-#include <memory>
-
 #include "FWCore/Utilities/interface/Exception.h"
 #include "FWCore/Utilities/interface/isFinite.h"
 #include "RecoJets/JetProducers/interface/JetSpecific.h"
@@ -36,9 +34,8 @@ namespace {
 }  // namespace
 
 SubEventGenJetProducer::SubEventGenJetProducer(edm::ParameterSet const& conf) : VirtualJetProducer(conf) {
-  //   mapSrc_ = conf.getParameter<edm::InputTag>( "srcMap");
   ignoreHydro_ = conf.getUntrackedParameter<bool>("ignoreHydro", true);
-  produces<reco::BasicJetCollection>();
+
   // the subjet collections are set through the config file in the "jetCollInstanceName" field.
 
   input_cand_token_ = consumes<reco::CandidateView>(src_);
@@ -104,8 +101,7 @@ void SubEventGenJetProducer::produce(edm::Event& iEvent, const edm::EventSetup&
 
   ////////////////
 
-  auto jets = std::make_unique<std::vector<GenJet>>();
-  subJets_ = jets.get();
+  jets_ = std::make_unique<std::vector<GenJet>>();
 
   LogDebug("VirtualJetProducer") << "Inputted towers\n";
 
@@ -123,7 +119,7 @@ void SubEventGenJetProducer::produce(edm::Event& iEvent, const edm::EventSetup&
   //Finalize
   LogDebug("SubEventJetProducer") << "Wrote jets\n";
 
-  iEvent.put(std::move(jets));
+  iEvent.put(std::move(jets_));
   return;
 }
 
@@ -156,7 +152,7 @@ void SubEventGenJetProducer::runAlgorithm(edm::Event& iEvent, edm::EventSetup co
     jet.setJetArea(jetArea);
     jet.setPileup(pu);
 
-    subJets_->push_back(jet);
+    jets_->push_back(jet);
   }
 }
 

RecoJets/JetProducers/plugins/SubEventGenJetProducer.h

@@ -10,6 +10,7 @@
 
  ************************************************************/
 
+#include <memory>
 #include <vector>
 #include "RecoJets/JetProducers/plugins/VirtualJetProducer.h"
 #include "DataFormats/JetReco/interface/GenJetCollection.h"
@@ -23,8 +24,8 @@ namespace cms {
     void runAlgorithm(edm::Event&, const edm::EventSetup&) override;
 
   protected:
-    std::vector<std::vector<fastjet::PseudoJet> > subInputs_;
-    std::vector<reco::GenJet>* subJets_;
+    std::vector<std::vector<fastjet::PseudoJet>> subInputs_;
+    std::unique_ptr<std::vector<reco::GenJet>> jets_;
     std::vector<int> hydroTag_;
     std::vector<int> nSubParticles_;
     bool ignoreHydro_;