Process modifier for 2018 PbPb miniAOD reprocessing

CommonTools/ParticleFlow/python/pfEGammaToCandidateRemapper_cfi.py

@@ -7,3 +7,5 @@
     photon2pf = cms.InputTag("particleBasedIsolation","gedPhotons"), 
     pf2pf = cms.InputTag("FILLME")
 )
+
+pfEGammaToCandidateRemapperCleaned = pfEGammaToCandidateRemapper.clone(pf2pf = "cleanedParticleFlow")

Configuration/ProcessModifiers/python/run2_miniAOD_pp_on_AA_103X_Bugfix_cff.py

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

Configuration/PyReleaseValidation/python/relval_steps.py

@@ -1958,7 +1958,7 @@ def lhegensim2018ml(fragment,howMuch):
                           '--datatier':'MINIAOD,DQMIO',
                           '--eventcontent':'MINIAOD,DQM',
                           '--era':'Run2_2018_pp_on_AA',
-                          '--procModifiers':'run2_miniAOD_pp_on_AA_103X',
+                          '--procModifiers':'run2_miniAOD_pp_on_AA_103X,run2_miniAOD_pp_on_AA_103X_Bugfix',
                           '--data':'',
                           '--processName':'PAT',
                           '-n':'100'
@@ -2507,14 +2507,14 @@ def gen2021HiMix(fragment,howMuch):
                                        '--eventcontent':'MINIAODSIM,DQM',
                                        '-n':100,
                                        '--era':'Run2_2018_pp_on_AA',
-                                       '--procModifiers':'run2_miniAOD_pp_on_AA_103X',
+                                       '--procModifiers':'run2_miniAOD_pp_on_AA_103X,run2_miniAOD_pp_on_AA_103X_Bugfix',
                                    },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,run2_miniAOD_pp_on_AA_103X',
+                                         '--procModifiers':'genJetSubEvent,run2_miniAOD_pp_on_AA_103X,run2_miniAOD_pp_on_AA_103X_Bugfix',
                                      },hiDefaults2018_ppReco,step3Up2015Defaults])
 steps['ALCARECOHI2018PPRECO']=merge([hiDefaults2018_ppReco,{'-s':'ALCA:TkAlMinBias+SiStripCalMinBias',
                                                             '--datatier':'ALCARECO',

PhysicsTools/PatAlgos/python/slimming/slimming_cff.py

@@ -30,6 +30,7 @@
 from RecoLuminosity.LumiProducer.bunchSpacingProducer_cfi import bunchSpacingProducer
 from HeavyFlavorAnalysis.Onia2MuMu.OniaPhotonConversionProducer_cfi import PhotonCandidates as oniaPhotonCandidates
 from RecoLocalCalo.HcalRecProducers.HcalHitSelection_cfi import *
+from CommonTools.ParticleFlow.pfEGammaToCandidateRemapper_cfi import *
 
 slimmingTask = cms.Task(
     packedPFCandidatesTask,
@@ -82,6 +83,8 @@
     cms.Task(slimmingTask.copy(), packedCandidateMuonID, packedPFCandidateTrackChi2, lostTrackChi2, centralityBin, hiHFfilters))
 from Configuration.ProcessModifiers.run2_miniAOD_pp_on_AA_103X_cff import run2_miniAOD_pp_on_AA_103X
 run2_miniAOD_pp_on_AA_103X.toReplaceWith(slimmingTask,cms.Task(primaryVertexAssociationCleaned,slimmingTask.copy()))
+from Configuration.ProcessModifiers.run2_miniAOD_pp_on_AA_103X_Bugfix_cff import run2_miniAOD_pp_on_AA_103X_Bugfix
+run2_miniAOD_pp_on_AA_103X_Bugfix.toReplaceWith(slimmingTask,cms.Task(pfEGammaToCandidateRemapperCleaned,slimmingTask.copy()))
 
 from RecoHI.HiTracking.miniAODVertexRecovery_cff import offlinePrimaryVerticesRecovery, offlineSlimmedPrimaryVerticesRecovery
 pp_on_AA.toReplaceWith(

RecoEgamma/EgammaPhotonProducers/python/reducedEgamma_cfi.py

@@ -128,3 +128,10 @@
     hiPhotonIsolationMapInput = "photonIsolationHIProducerppGED",
     hiPhotonIsolationMapOutput = "photonIsolationHIProducerppGED"
 )
+
+from Configuration.ProcessModifiers.run2_miniAOD_pp_on_AA_103X_Bugfix_cff import run2_miniAOD_pp_on_AA_103X_Bugfix
+run2_miniAOD_pp_on_AA_103X_Bugfix.toModify(
+    reducedEgamma,
+    photonsPFValMap = "pfEGammaToCandidateRemapperCleaned:photons",
+    gsfElectronsPFValMap = "pfEGammaToCandidateRemapperCleaned:electrons"
+)

RecoHI/HiJetAlgos/plugins/HiBadParticleCleaner.cc

@@ -32,7 +32,7 @@ class HiBadParticleCleaner : public edm::global::EDProducer<> {
   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
   // ----------member data ---------------------------
 
-  edm::EDGetTokenT<edm::View<reco::PFCandidate> > tokenPFCandidates_;
+  edm::EDGetTokenT<edm::View<reco::PFCandidate>> tokenPFCandidates_;
   edm::EDGetTokenT<reco::VertexCollection> tokenPV_;
 
   const double minMuonPt_;
@@ -45,13 +45,14 @@ class HiBadParticleCleaner : public edm::global::EDProducer<> {
   const unsigned minPixelNHits_;
   const int minTrackerLayersForMuonLoose_;
   const int minTrackerLayersForMuonTight_;
+  const bool reMiniAODBugFix_;
 };
 
 //
 // constructors and destructor
 //
 HiBadParticleCleaner::HiBadParticleCleaner(const edm::ParameterSet& iConfig)
-    : tokenPFCandidates_(consumes<edm::View<reco::PFCandidate> >(iConfig.getParameter<edm::InputTag>("PFCandidates"))),
+    : tokenPFCandidates_(consumes<edm::View<reco::PFCandidate>>(iConfig.getParameter<edm::InputTag>("PFCandidates"))),
       tokenPV_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("offlinePV"))),
       minMuonPt_(iConfig.getParameter<double>("minMuonPt")),
       minChargedHadronPt_(iConfig.getParameter<double>("minChargedHadronPt")),
@@ -62,10 +63,12 @@ HiBadParticleCleaner::HiBadParticleCleaner(const edm::ParameterSet& iConfig)
       minTrackNHits_(iConfig.getParameter<uint>("minTrackNHits")),
       minPixelNHits_(iConfig.getParameter<uint>("minPixelNHits")),
       minTrackerLayersForMuonLoose_(iConfig.getParameter<int>("minTrackerLayersForMuonLoose")),
-      minTrackerLayersForMuonTight_(iConfig.getParameter<int>("minTrackerLayersForMuonTight")) {
+      minTrackerLayersForMuonTight_(iConfig.getParameter<int>("minTrackerLayersForMuonTight")),
+      reMiniAODBugFix_(iConfig.getParameter<bool>("reMiniAODBugFix")) {
   produces<bool>();
   produces<reco::PFCandidateCollection>();
   produces<reco::PFCandidateCollection>("removed");
+  produces<edm::ValueMap<reco::PFCandidateRef>>();
 }
 
 //
@@ -90,19 +93,27 @@ void HiBadParticleCleaner::produce(edm::StreamID, edm::Event& iEvent, const edm:
   auto pBadCandidateCollection = std::make_unique<reco::PFCandidateCollection>();
 
   bool foundBadCandidate = false;
+  size_t n = pfCandidates->size();
+  std::vector<int> candidateIndexMapper(n, 0);  // mapping between the original PF and post-cleaning PF
+  size_t iPF;
+  for (iPF = 0; iPF < n; iPF++) {
+    const reco::PFCandidate& pfCandidate = pfCandidates->at(iPF);
 
-  for (const reco::PFCandidate& pfCandidate : *pfCandidates) {
     if (pfCandidate.particleId() == reco::PFCandidate::ParticleType::mu)  // muon cleaning
     {
       if (pfCandidate.pt() > minMuonPt_) {
         if (!pfCandidate.muonRef()->isGlobalMuon() || !pfCandidate.muonRef()->isTrackerMuon() ||
             !pfCandidate.trackRef().isNonnull()) {
+          pBadCandidateCollection->push_back(pfCandidate);
+          candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
           foundBadCandidate = true;
           continue;
         }
         reco::TrackRef track = pfCandidate.trackRef();
 
         if (track->ptError() / track->pt() > minMuonTrackRelPtErr_ || track->pt() < pfCandidate.pt() / 2.) {
+          pBadCandidateCollection->push_back(pfCandidate);
+          candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
           foundBadCandidate = true;
           continue;
         }
@@ -128,19 +139,22 @@ void HiBadParticleCleaner::produce(edm::StreamID, edm::Event& iEvent, const edm:
 
           if (sig3d > maxSigLoose_) {
             pBadCandidateCollection->push_back(pfCandidate);
+            candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
             foundBadCandidate = true;
             continue;
           }
 
           if (track->pt() < pfCandidate.pt() / 1.5 || track->pt() > pfCandidate.pt() * 1.5) {
             foundBadCandidate = true;
             pBadCandidateCollection->push_back(pfCandidate);
+            candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
             continue;
           }
           if (track->originalAlgo() == reco::TrackBase::muonSeededStepOutIn &&
               track->hitPattern().trackerLayersWithMeasurement() < minTrackerLayersForMuonTight_) {
             foundBadCandidate = true;
             pBadCandidateCollection->push_back(pfCandidate);
+            candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
             continue;
           }
         }
@@ -156,6 +170,7 @@ void HiBadParticleCleaner::produce(edm::StreamID, edm::Event& iEvent, const edm:
         if ((nHits < minTrackNHits_ && nPixelHits < minPixelNHits_) || nHits == 3) {
           foundBadCandidate = true;
           pBadCandidateCollection->push_back(pfCandidate);
+          candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
           continue;
         }
 
@@ -176,12 +191,14 @@ void HiBadParticleCleaner::produce(edm::StreamID, edm::Event& iEvent, const edm:
         if (sig3d > maxSigLoose_) {
           foundBadCandidate = true;
           pBadCandidateCollection->push_back(pfCandidate);
+          candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
           continue;
         }
 
         if (sig3d > maxSigTight_ && nHits < minTrackNHits_) {
           foundBadCandidate = true;
           pBadCandidateCollection->push_back(pfCandidate);
+          candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
           continue;
         }
 
@@ -192,12 +209,14 @@ void HiBadParticleCleaner::produce(edm::StreamID, edm::Event& iEvent, const edm:
           if (sig3d > maxSigLoose_) {
             foundBadCandidate = true;
             pBadCandidateCollection->push_back(pfCandidate);
+            candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
             continue;
           }
 
           if (nHits < minTrackNHits_) {
             foundBadCandidate = true;
             pBadCandidateCollection->push_back(pfCandidate);
+            candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
             continue;
           }
         }
@@ -208,34 +227,56 @@ void HiBadParticleCleaner::produce(edm::StreamID, edm::Event& iEvent, const edm:
           if (sig3d > maxSigTight_) {
             foundBadCandidate = true;
             pBadCandidateCollection->push_back(pfCandidate);
+            candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
             continue;
           }
 
           if (nHits < minTrackNHits_) {
             foundBadCandidate = true;
             pBadCandidateCollection->push_back(pfCandidate);
+            candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
             continue;
           }
 
           if (nPixelHits < minPixelNHits_) {
             foundBadCandidate = true;
             pBadCandidateCollection->push_back(pfCandidate);
+            candidateIndexMapper[iPF] = -1 * (pBadCandidateCollection->size());
             continue;
           }
         }
       }
     }
 
     pOutputCandidateCollection->push_back(pfCandidate);
-
+    candidateIndexMapper[iPF] = (pOutputCandidateCollection->size());
   }  // end loop over pf candidates
 
   bool pass = !foundBadCandidate;
-
-  iEvent.put(std::move(pOutputCandidateCollection));
-  iEvent.put(std::move(pBadCandidateCollection), "removed");
+  edm::OrphanHandle<std::vector<reco::PFCandidate>> newpf = iEvent.put(std::move(pOutputCandidateCollection));
+  edm::OrphanHandle<std::vector<reco::PFCandidate>> badpf = iEvent.put(std::move(pBadCandidateCollection), "removed");
 
   iEvent.put(std::make_unique<bool>(pass));
+
+  if (reMiniAODBugFix_) {
+    std::unique_ptr<edm::ValueMap<reco::PFCandidateRef>> pf2pf(new edm::ValueMap<reco::PFCandidateRef>());
+    edm::ValueMap<reco::PFCandidateRef>::Filler filler(*pf2pf);
+
+    std::vector<reco::PFCandidateRef> refs;
+    refs.reserve(n);
+
+    for (iPF = 0; iPF < n; ++iPF) {
+      if (candidateIndexMapper[iPF] > 0) {
+        refs.push_back(reco::PFCandidateRef(newpf, candidateIndexMapper[iPF] - 1));
+      } else if (candidateIndexMapper[iPF] < 0) {
+        refs.push_back(reco::PFCandidateRef(badpf, -candidateIndexMapper[iPF] - 1));
+      }
+    }
+    filler.insert(pfCandidates, refs.begin(), refs.end());
+
+    filler.fill();
+    iEvent.put(std::move(pf2pf));
+  }
 }
 
 //define this as a plug-in

RecoHI/HiJetAlgos/plugins/HiPuRhoProducer.cc

@@ -90,6 +90,7 @@ class HiPuRhoProducer : public edm::stream::EDProducer<> {
   bool setInitialValue_;
 
   const bool dropZeroTowers_;
+  const bool reMiniAODBugFix_;
   const int medianWindowWidth_;
   const double minimumTowersFraction_;
   const double nSigmaPU_;  // number of sigma for pileup
@@ -149,6 +150,7 @@ class HiPuRhoProducer : public edm::stream::EDProducer<> {
 
 HiPuRhoProducer::HiPuRhoProducer(const edm::ParameterSet& iConfig)
     : dropZeroTowers_(iConfig.getParameter<bool>("dropZeroTowers")),
+      reMiniAODBugFix_(iConfig.getParameter<bool>("reMiniAODBugFix")),
       medianWindowWidth_(iConfig.getParameter<int>("medianWindowWidth")),
       minimumTowersFraction_(iConfig.getParameter<double>("minimumTowersFraction")),
       nSigmaPU_(iConfig.getParameter<double>("nSigmaPU")),
@@ -464,8 +466,11 @@ void HiPuRhoProducer::calculateOrphanInput(std::vector<fastjet::PseudoJet>& orph
 
     for (auto const& im : towermap_) {
       double dr2 = reco::deltaR2(im.eta, im.phi, jet_etaphi.first, jet_etaphi.second);
+      double nTowers = im.ieta;
+      if (reMiniAODBugFix_)
+        nTowers = geomtowers_[im.ieta];
       if (dr2 < radiusPU_ * radiusPU_ && !excludedTowers[std::pair(im.ieta, im.iphi)] &&
-          (im.ieta - ntowersWithJets_[im.ieta]) > minimumTowersFraction_ * im.ieta) {
+          (nTowers - ntowersWithJets_[im.ieta]) > minimumTowersFraction_ * nTowers) {
         ntowersWithJets_[im.ieta]++;
         excludedTowers[std::pair(im.ieta, im.iphi)] = 1;
       }
@@ -572,6 +577,7 @@ void HiPuRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descripti
   desc.add<double>("radiusPU", 0.5);
   desc.add<double>("minimumTowersFraction", 0.7);
   desc.add<bool>("dropZeroTowers", true);
+  desc.add<bool>("reMiniAODBugFix", false);
   descriptions.add("hiPuRhoProducer", desc);
 }
 

RecoHI/HiJetAlgos/python/HiBadParticleCleaner_cfi.py

@@ -14,5 +14,9 @@
     minTrackNHits = cms.uint32(10),
     minPixelNHits = cms.uint32(3),
     minTrackerLayersForMuonLoose  = cms.int32(7),  
-    minTrackerLayersForMuonTight  = cms.int32(10)  
+    minTrackerLayersForMuonTight  = cms.int32(10),
+    reMiniAODBugFix = cms.bool(False)
 )
+
+from Configuration.ProcessModifiers.run2_miniAOD_pp_on_AA_103X_Bugfix_cff import run2_miniAOD_pp_on_AA_103X_Bugfix
+run2_miniAOD_pp_on_AA_103X_Bugfix.toModify(cleanedParticleFlow,reMiniAODBugFix = True)

RecoHI/HiJetAlgos/python/HiRecoPFJets_cff.py

@@ -120,3 +120,5 @@
 from Configuration.ProcessModifiers.run2_miniAOD_pp_on_AA_103X_cff import run2_miniAOD_pp_on_AA_103X
 run2_miniAOD_pp_on_AA_103X.toModify(akCs4PFJets,src = 'cleanedParticleFlow')
 run2_miniAOD_pp_on_AA_103X.toModify(PFTowers,src = 'cleanedParticleFlow')
+from Configuration.ProcessModifiers.run2_miniAOD_pp_on_AA_103X_Bugfix_cff import run2_miniAOD_pp_on_AA_103X_Bugfix
+run2_miniAOD_pp_on_AA_103X.toModify(hiPuRho,reMiniAODBugFix = True)