From e8657815f4596a3349308006bf8981f9d004ca6b Mon Sep 17 00:00:00 2001 From: Stepan Date: Fri, 3 Jul 2020 22:54:31 +0200 Subject: [PATCH] apply suggestions --- .../plugins/HiFJRhoFlowModulationProducer.cc | 161 +++++----- RecoHI/HiJetAlgos/plugins/HiPuRhoProducer.cc | 282 ++++++++---------- RecoHI/HiJetAlgos/python/HiRecoPFJets_cff.py | 67 ++++- .../Configuration/python/RecoPFJets_cff.py | 7 +- .../JetProducers/plugins/CSJetProducer.cc | 32 +- RecoJets/JetProducers/plugins/CSJetProducer.h | 4 +- 6 files changed, 276 insertions(+), 277 deletions(-) diff --git a/RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer.cc b/RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer.cc index 53948fb73129e..3485a4c3808ff 100644 --- a/RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer.cc +++ b/RecoHI/HiJetAlgos/plugins/HiFJRhoFlowModulationProducer.cc @@ -26,6 +26,7 @@ #include "TF1.h" #include "TH1.h" #include "TMath.h" +#include "TMinuitMinimizer.h" #include #include @@ -34,52 +35,59 @@ #include #include +namespace { + double lineFunction(double* x, double* par) { return par[0]; } + double flowFunction(double* x, double* par) { + return par[0] * (1. + 2. * (par[1] * TMath::Cos(2. * (x[0] - par[2])) + par[3] * TMath::Cos(3. * (x[0] - par[4])))); + } +}; // namespace + class HiFJRhoFlowModulationProducer : public edm::stream::EDProducer<> { public: explicit HiFJRhoFlowModulationProducer(const edm::ParameterSet&); - static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); private: - void beginStream(edm::StreamID) override; void produce(edm::Event&, const edm::EventSetup&) override; - void endStream() override; - bool doEvtPlane_; - bool doFreePlaneFit_; - bool doJettyExclusion_; - int evtPlaneLevel_; + const bool doEvtPlane_; + const bool doFreePlaneFit_; + const bool doJettyExclusion_; + const int evtPlaneLevel_; edm::EDGetTokenT jetTag_; - edm::EDGetTokenT pfCandsToken_; - edm::EDGetTokenT evtPlaneToken_; - - float* hiEvtPlane; + const edm::EDGetTokenT pfCandsToken_; + const edm::EDGetTokenT evtPlaneToken_; + static constexpr int kMaxEvtPlane = 1000; + std::array hiEvtPlane_; + std::unique_ptr lineFit_p; + std::unique_ptr flowFit_p; + int nPhiBins; }; - HiFJRhoFlowModulationProducer::HiFJRhoFlowModulationProducer(const edm::ParameterSet& iConfig) : doEvtPlane_(iConfig.getParameter("doEvtPlane")), doFreePlaneFit_(iConfig.getParameter("doFreePlaneFit")), doJettyExclusion_(iConfig.getParameter("doJettyExclusion")), evtPlaneLevel_(iConfig.getParameter("evtPlaneLevel")), - jetTag_(consumes(iConfig.getParameter("jetTag"))), pfCandsToken_(consumes(iConfig.getParameter("pfCandSource"))), evtPlaneToken_(consumes(iConfig.getParameter("EvtPlane"))) { + if (doJettyExclusion_) + jetTag_ = consumes(iConfig.getParameter("jetTag")); produces>("rhoFlowFitParams"); + TMinuitMinimizer::UseStaticMinuit(false); + lineFit_p = std::unique_ptr(new TF1("lineFit", lineFunction, -TMath::Pi(), TMath::Pi())); + flowFit_p = std::unique_ptr(new TF1("flowFit", flowFunction, -TMath::Pi(), TMath::Pi())); } // ------------ method called to produce the data ------------ void HiFJRhoFlowModulationProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { - edm::Handle pfCands; - iEvent.getByToken(pfCandsToken_, pfCands); + auto const& pfCands = iEvent.get(pfCandsToken_); - edm::Handle evtPlanes; if (doEvtPlane_) { - iEvent.getByToken(evtPlaneToken_, evtPlanes); - if (evtPlanes.isValid()) { - for (unsigned int i = 0; i < evtPlanes->size(); ++i) { - hiEvtPlane[i] = (*evtPlanes)[i].angle(evtPlaneLevel_); - } - } + auto const& evtPlanes = iEvent.get(evtPlaneToken_); + assert(evtPlanes.size() < kMaxEvtPlane); + std::transform(evtPlanes.begin(), evtPlanes.end(), hiEvtPlane_.begin(), [this](auto const& ePlane) -> float { + return ePlane.angle(evtPlaneLevel_); + }); } edm::Handle jets; @@ -106,17 +114,28 @@ void HiFJRhoFlowModulationProducer::produce(edm::Event& iEvent, const edm::Event double eventPlane3Cos = 0; double eventPlane3Sin = 0; - for (auto const& pfCandidate : *pfCands) { - if (pfCandidate.eta() < -1.0) + std::vector pfcuts; + for (auto const& pfCandidate : pfCands) { + if (pfCandidate.particleId() != 1) { + pfcuts.push_back(false); continue; - if (pfCandidate.eta() > 1.0) + } + if (pfCandidate.eta() < -1.0) { + pfcuts.push_back(false); continue; - if (pfCandidate.pt() < .3) + } + if (pfCandidate.eta() > 1.0) { + pfcuts.push_back(false); continue; - if (pfCandidate.pt() > 3.) + } + if (pfCandidate.pt() < .3) { + pfcuts.push_back(false); continue; - if (pfCandidate.particleId() != 1) + } + if (pfCandidate.pt() > 3.) { + pfcuts.push_back(false); continue; + } if (doJettyExclusion_) { bool isGood = true; @@ -126,11 +145,14 @@ void HiFJRhoFlowModulationProducer::produce(edm::Event& iEvent, const edm::Event break; } } - if (!isGood) + if (!isGood) { + pfcuts.push_back(false); continue; + } } nFill++; + pfcuts.push_back(true); if (!doEvtPlane_) { eventPlane2Cos += std::cos(2 * pfCandidate.phi()); @@ -145,65 +167,35 @@ void HiFJRhoFlowModulationProducer::produce(edm::Event& iEvent, const edm::Event eventPlane2 = std::atan2(eventPlane2Sin, eventPlane2Cos) / 2.; eventPlane3 = std::atan2(eventPlane3Sin, eventPlane3Cos) / 3.; } else { - eventPlane2 = hiEvtPlane[8]; - eventPlane3 = hiEvtPlane[15]; + eventPlane2 = hiEvtPlane_[8]; + eventPlane3 = hiEvtPlane_[15]; } - + int pfcuts_count = 0; if (nFill >= 100 && eventPlane2 > -99) { - const int nPhiBins = std::max(10, nFill / 30); + nPhiBins = std::max(10, nFill / 30); std::string name = "phiTestIEta4_" + std::to_string(iEvent.id().event()) + "_h"; std::string nameFlat = "phiTestIEta4_Flat_" + std::to_string(iEvent.id().event()) + "_h"; - TH1F* phi_h = new TH1F(name.data(), "", nPhiBins, -TMath::Pi(), TMath::Pi()); - - for (auto const& pfCandidate : *pfCands) { - if (pfCandidate.eta() < -1.0) - continue; - if (pfCandidate.eta() > 1.0) - continue; - if (pfCandidate.pt() < .3) - continue; - if (pfCandidate.pt() > 3.) - continue; - if (pfCandidate.particleId() != 1) - continue; - - if (doJettyExclusion_) { - bool isGood = true; - for (auto const& jet : *jets) { - if (deltaR2(jet, pfCandidate) < .16) { - isGood = false; - break; - } - } - if (!isGood) - continue; - } - - phi_h->Fill(pfCandidate.phi()); + for (auto const& pfCandidate : pfCands) { + if (pfcuts.at(pfcuts_count)) + phi_h->Fill(pfCandidate.phi()); + pfcuts_count++; } - - std::string flowFitForm = "[0]*(1.+2.*([1]*TMath::Cos(2.*(x-[2]))+[3]*TMath::Cos(3.*(x-[4]))))"; - - TF1* flowFit_p = new TF1("flowFit", flowFitForm.c_str(), -TMath::Pi(), TMath::Pi()); flowFit_p->SetParameter(0, 10); flowFit_p->SetParameter(1, 0); flowFit_p->SetParameter(2, eventPlane2); flowFit_p->SetParameter(3, 0); flowFit_p->SetParameter(4, eventPlane3); - if (!doFreePlaneFit_) { flowFit_p->FixParameter(2, eventPlane2); flowFit_p->FixParameter(4, eventPlane3); } - TF1* lineFit_p = new TF1("lineFit", "[0]", -TMath::Pi(), TMath::Pi()); lineFit_p->SetParameter(0, 10); - phi_h->Fit(flowFit_p, "Q", "", -TMath::Pi(), TMath::Pi()); - phi_h->Fit(lineFit_p, "Q", "", -TMath::Pi(), TMath::Pi()); - + phi_h->Fit(flowFit_p.get(), "Q SERIAL", "", -TMath::Pi(), TMath::Pi()); + phi_h->Fit(lineFit_p.get(), "Q SERIAL", "", -TMath::Pi(), TMath::Pi()); rhoFlowFitParamsOut->at(0) = flowFit_p->GetParameter(0); rhoFlowFitParamsOut->at(1) = flowFit_p->GetParameter(1); rhoFlowFitParamsOut->at(2) = flowFit_p->GetParameter(2); @@ -217,35 +209,24 @@ void HiFJRhoFlowModulationProducer::produce(edm::Event& iEvent, const edm::Event rhoFlowFitParamsOut->at(8) = lineFit_p->GetNDF(); delete phi_h; - - delete flowFit_p; - delete lineFit_p; + pfcuts.clear(); } iEvent.put(std::move(rhoFlowFitParamsOut), "rhoFlowFitParams"); } -// ------------ method called once each stream just before starting event loop ------------ -void HiFJRhoFlowModulationProducer::beginStream(edm::StreamID) { - if (doEvtPlane_) { - constexpr int kMaxEvtPlanes = 1000; - hiEvtPlane = new float[kMaxEvtPlanes]; - } -} - -// ------------ method called once each stream just after ending the event loop ------------ -void HiFJRhoFlowModulationProducer::endStream() { - if (doEvtPlane_) - delete hiEvtPlane; -} - // ------------ method fills 'descriptions' with the allowed parameters for the module ------------ void HiFJRhoFlowModulationProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - // The following says we do not know what parameters are allowed so do no validation - // Please change this to state exactly what you do use, even if it is no parameters edm::ParameterSetDescription desc; - desc.setUnknown(); - descriptions.addDefault(desc); + desc.add("doEvtPlane", false); + desc.add("EvtPlane", edm::InputTag("hiEvtPlane")); + desc.add("doJettyExclusion", false); + desc.add("doFreePlaneFit", false); + desc.add("doFlatTest", false); + desc.add("jetTag", edm::InputTag("ak4PFJets")); + desc.add("pfCandSource", edm::InputTag("particleFlow")); + desc.add("evtPlaneLevel", 0); + descriptions.add("hiFJRhoFlowModulationProducer", desc); } DEFINE_FWK_MODULE(HiFJRhoFlowModulationProducer); diff --git a/RecoHI/HiJetAlgos/plugins/HiPuRhoProducer.cc b/RecoHI/HiJetAlgos/plugins/HiPuRhoProducer.cc index 9a95a02f3b11f..9f3c5ea400f2c 100644 --- a/RecoHI/HiJetAlgos/plugins/HiPuRhoProducer.cc +++ b/RecoHI/HiJetAlgos/plugins/HiPuRhoProducer.cc @@ -47,8 +47,6 @@ #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" -#include "TMath.h" - #include #include #include @@ -75,6 +73,26 @@ class HiPuRhoProducer : public edm::stream::EDProducer<> { virtual void calculateOrphanInput(std::vector& orphanInput); virtual void putRho(edm::Event& iEvent, const edm::EventSetup& iSetup); +private: + void produce(edm::Event&, const edm::EventSetup&) override; + + // This checks if the tower is anomalous (if a calo tower). + virtual void inputTowers(); + + bool postOrphan_; + bool setInitialValue_; + + const bool dropZeroTowers_; + const int medianWindowWidth_; + const double minimumTowersFraction_; + const double nSigmaPU_; // number of sigma for pileup + const double puPtMin_; + const double radiusPU_; // pileup radius + const double rParam_; // the R parameter to use + const double towSigmaCut_; + const edm::EDGetTokenT input_candidateview_token_; + edm::ESGetToken caloToken_; + constexpr static int nMaxJets_ = 200; float jteta[nMaxJets_]; @@ -89,8 +107,6 @@ class HiPuRhoProducer : public edm::stream::EDProducer<> { int vngeom[nEtaTow_]; int vntow[nEtaTow_]; - int vieta[nEtaTow_]; - float veta[nEtaTow_]; float vmean0[nEtaTow_]; float vrms0[nEtaTow_]; float vrho0[nEtaTow_]; @@ -102,6 +118,7 @@ class HiPuRhoProducer : public edm::stream::EDProducer<> { 0.957, 1.044, 1.131, 1.218, 1.305, 1.392, 1.479, 1.566, 1.653, 1.740, 1.830, 1.930, 2.043, 2.172, 2.322, 2.500, 2.650, 2.853, 3.000, 3.139, 3.314, 3.489, 3.664, 3.839, 4.013, 4.191, 4.363, 4.538, 4.716, 4.889, 5.191}; + const int initialValue = -99; std::vector etaEdgeLow_; std::vector etaEdgeHi_; @@ -116,30 +133,9 @@ class HiPuRhoProducer : public edm::stream::EDProducer<> { std::vector towExcludePhi_; std::vector towExcludeEta_; - int ieta(const reco::CandidatePtr& in) const; - int iphi(const reco::CandidatePtr& in) const; - -private: - void produce(edm::Event&, const edm::EventSetup&) override; - - // This checks if the tower is anomalous (if a calo tower). - virtual void inputTowers(); - - bool postOrphan_; - - bool dropZeroTowers_; - int medianWindowWidth_; - double minimumTowersFraction_; - double nSigmaPU_; // number of sigma for pileup - double puPtMin_; - double radiusPU_; // pileup radius - double rParam_; // the R parameter to use - edm::InputTag src_; // input constituent source - double towSigmaCut_; - - std::vector> inputs_; // input candidates - ClusterSequencePtr fjClusterSeq_; // fastjet cluster sequence - JetDefPtr fjJetDefinition_; // fastjet jet definition + std::vector inputs_; // input candidates + ClusterSequencePtr fjClusterSeq_; // fastjet cluster sequence + JetDefPtr fjJetDefinition_; // fastjet jet definition std::vector fjInputs_; // fastjet inputs std::vector fjJets_; // fastjet jets @@ -156,7 +152,8 @@ class HiPuRhoProducer : public edm::stream::EDProducer<> { std::map emean_; // energy mean std::map> eTop4_; // energy mean - edm::EDGetTokenT input_candidateview_token_; + typedef std::pair EtaPhi; + std::map towermap; }; HiPuRhoProducer::HiPuRhoProducer(const edm::ParameterSet& iConfig) @@ -167,10 +164,9 @@ HiPuRhoProducer::HiPuRhoProducer(const edm::ParameterSet& iConfig) puPtMin_(iConfig.getParameter("puPtMin")), radiusPU_(iConfig.getParameter("radiusPU")), rParam_(iConfig.getParameter("rParam")), - src_(iConfig.getParameter("src")), - towSigmaCut_(iConfig.getParameter("towSigmaCut")) { - input_candidateview_token_ = consumes(src_); - + towSigmaCut_(iConfig.getParameter("towSigmaCut")), + input_candidateview_token_(consumes(iConfig.getParameter("src"))), + caloToken_(esConsumes(edm::ESInputTag{})) { //register your products produces>("mapEtaEdges"); produces>("mapToRho"); @@ -187,21 +183,19 @@ HiPuRhoProducer::HiPuRhoProducer(const edm::ParameterSet& iConfig) void HiPuRhoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { setupGeometryMap(iEvent, iSetup); - inputs_.clear(); - fjInputs_.clear(); - fjJets_.clear(); - - // Get the particletowers - edm::Handle inputsHandle; - iEvent.getByToken(input_candidateview_token_, inputsHandle); + for (int i = ietamin_; i < ietamax_ + 1; i++) { + ntowersWithJets_[i] = 0; + } - for (size_t i = 0; i < inputsHandle->size(); ++i) - inputs_.push_back(inputsHandle->ptrAt(i)); + auto const& inputView = iEvent.get(input_candidateview_token_); + inputs_.reserve(inputView.size()); + for (auto const& input : inputView) + inputs_.push_back(&input); fjInputs_.reserve(inputs_.size()); inputTowers(); fjOriginalInputs_ = fjInputs_; - + setInitialValue_ = true; calculatePedestal(fjInputs_); subtractPedestal(fjInputs_); @@ -222,17 +216,21 @@ void HiPuRhoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) towExcludePhi_.clear(); towExcludeEta_.clear(); + setInitialValue_ = false; std::vector orphanInput; calculateOrphanInput(orphanInput); calculatePedestal(orphanInput); putRho(iEvent, iSetup); + + inputs_.clear(); + fjInputs_.clear(); + fjJets_.clear(); } void HiPuRhoProducer::inputTowers() { - auto inBegin = inputs_.begin(); - auto inEnd = inputs_.end(); - for (auto i = inBegin; i != inEnd; ++i) { - reco::CandidatePtr input = *i; + int index = -1; + for (auto const& input : inputs_) { + index++; if (edm::isNotFinite(input->pt())) continue; @@ -240,44 +238,38 @@ void HiPuRhoProducer::inputTowers() { continue; fjInputs_.push_back(fastjet::PseudoJet(input->px(), input->py(), input->pz(), input->energy())); - fjInputs_.back().set_user_index(i - inBegin); + fjInputs_.back().set_user_index(index); } } void HiPuRhoProducer::setupGeometryMap(edm::Event& iEvent, const edm::EventSetup& iSetup) { LogDebug("PileUpSubtractor") << "The subtractor setting up geometry...\n"; - - if (geo_ == nullptr) { - edm::ESHandle pG; - iSetup.get().get(pG); - geo_ = pG.product(); - std::vector alldid = geo_->getValidDetIds(); - - int ietaold = -10000; - ietamax_ = -10000; - ietamin_ = 10000; - for (auto const& did : alldid) { - if (did.det() == DetId::Hcal) { - HcalDetId hid = HcalDetId(did); - allgeomid_.push_back(did); - - if (hid.ieta() != ietaold) { - ietaold = hid.ieta(); - geomtowers_[hid.ieta()] = 1; - if (hid.ieta() > ietamax_) - ietamax_ = hid.ieta(); - if (hid.ieta() < ietamin_) - ietamin_ = hid.ieta(); - } else { - geomtowers_[hid.ieta()]++; - } + const auto& pG = iSetup.getData(caloToken_); + geo_ = &pG; + std::vector alldid = geo_->getValidDetIds(); + int ietaold = -10000; + ietamax_ = -10000; + ietamin_ = 10000; + towermap.clear(); + + for (auto const& did : alldid) { + if (did.det() == DetId::Hcal) { + HcalDetId hid = HcalDetId(did); + allgeomid_.push_back(did); + EtaPhi ep(geo_->getPosition((DetId)did).eta(), geo_->getPosition((DetId)did).phi()); + towermap[did] = ep; + if (hid.ieta() != ietaold) { + ietaold = hid.ieta(); + geomtowers_[hid.ieta()] = 1; + if (hid.ieta() > ietamax_) + ietamax_ = hid.ieta(); + if (hid.ieta() < ietamin_) + ietamin_ = hid.ieta(); + } else { + geomtowers_[hid.ieta()]++; } } } - - for (int i = ietamin_; i < ietamax_ + 1; i++) { - ntowersWithJets_[i] = 0; - } } void HiPuRhoProducer::calculatePedestal(std::vector const& coll) { @@ -294,21 +286,17 @@ void HiPuRhoProducer::calculatePedestal(std::vector const& c if (it > nEtaTow_ / 2) it = vi - nEtaTow_; - int sign = (it < 0) ? -1 : 1; - - vieta[vi] = it; - veta[vi] = sign * (etaedge[abs(it)] + etaedge[abs(it) - 1]) / 2.; - vngeom[vi] = -99; - vntow[vi] = -99; + vngeom[vi] = initialValue; + vntow[vi] = initialValue; - vmean1[vi] = -99; - vrms1[vi] = -99; - vrho1[vi] = -99; + vmean1[vi] = initialValue; + vrms1[vi] = initialValue; + vrho1[vi] = initialValue; - if ((*(ntowersWithJets_.find(it))).second == 0) { - vmean0[vi] = -99; - vrms0[vi] = -99; - vrho0[vi] = -99; + if (setInitialValue_) { + vmean0[vi] = initialValue; + vrms0[vi] = initialValue; + vrho0[vi] = initialValue; } } @@ -322,34 +310,33 @@ void HiPuRhoProducer::calculatePedestal(std::vector const& c } for (auto const& input_object : coll) { - const reco::CandidatePtr& originalTower = inputs_[input_object.user_index()]; - double Original_Et = originalTower->et(); - int ieta0 = ieta(originalTower); + const reco::Candidate* originalTower = inputs_[input_object.user_index()]; + double original_Et = originalTower->et(); + const CaloTower* ctc = dynamic_cast(originalTower); + int ieta0 = ctc->id().ieta(); - if (Original_Et > eTop4_[ieta0][0]) { + if (original_Et > eTop4_[ieta0][0]) { eTop4_[ieta0][3] = eTop4_[ieta0][2]; eTop4_[ieta0][2] = eTop4_[ieta0][1]; eTop4_[ieta0][1] = eTop4_[ieta0][0]; - eTop4_[ieta0][0] = Original_Et; - } else if (Original_Et > eTop4_[ieta0][1]) { + eTop4_[ieta0][0] = original_Et; + } else if (original_Et > eTop4_[ieta0][1]) { eTop4_[ieta0][3] = eTop4_[ieta0][2]; eTop4_[ieta0][2] = eTop4_[ieta0][1]; - eTop4_[ieta0][1] = Original_Et; - } else if (Original_Et > eTop4_[ieta0][2]) { + eTop4_[ieta0][1] = original_Et; + } else if (original_Et > eTop4_[ieta0][2]) { eTop4_[ieta0][3] = eTop4_[ieta0][2]; - eTop4_[ieta0][2] = Original_Et; - } else if (Original_Et > eTop4_[ieta0][3]) { - eTop4_[ieta0][3] = Original_Et; + eTop4_[ieta0][2] = original_Et; + } else if (original_Et > eTop4_[ieta0][3]) { + eTop4_[ieta0][3] = original_Et; } + emean_[ieta0] = emean_[ieta0] + original_Et; + emean2[ieta0] = emean2[ieta0] + original_Et * original_Et; if (ieta0 - ietaold != 0) { - emean_[ieta0] = emean_[ieta0] + Original_Et; - emean2[ieta0] = emean2[ieta0] + Original_Et * Original_Et; ntowers[ieta0] = 1; ietaold = ieta0; } else { - emean_[ieta0] = emean_[ieta0] + Original_Et; - emean2[ieta0] = emean2[ieta0] + Original_Et * Original_Et; ntowers[ieta0]++; } } @@ -362,9 +349,9 @@ void HiPuRhoProducer::calculatePedestal(std::vector const& c if (it < 0) vi = nEtaTow_ + it; - double e1 = (emean_.find(it))->second; - double e2 = (emean2.find(it))->second; - int nt = gt.second - (ntowersWithJets_.find(it))->second; + double e1 = emean_[it]; + double e2 = emean2[it]; + int nt = gt.second - ntowersWithJets_[it]; if (vi < nEtaTow_) { vngeom[vi] = gt.second; @@ -420,13 +407,13 @@ void HiPuRhoProducer::calculatePedestal(std::vector const& c if (vi < nEtaTow_) { vmean1[vi] = emean_[it]; - vrho1[vi] = emean_[it] / (etaWidth * (2. * TMath::Pi() / (double)vngeom[vi])); + vrho1[vi] = emean_[it] / (etaWidth * (2. * M_PI / (double)vngeom[vi])); rho_.push_back(vrho1[vi]); rhoM_.push_back(0); vrms1[vi] = esigma_[it]; if (vngeom[vi] == vntow[vi]) { vmean0[vi] = emean_[it]; - vrho0[vi] = emean_[it] / (etaWidth * (2. * TMath::Pi() / (double)vngeom[vi])); + vrho0[vi] = emean_[it] / (etaWidth * (2. * M_PI / (double)vngeom[vi])); vrms0[vi] = esigma_[it]; } rhoExtra_.push_back(vrho0[vi]); @@ -448,13 +435,13 @@ void HiPuRhoProducer::subtractPedestal(std::vector& coll) { std::vector newcoll; for (auto& input_object : coll) { int index = input_object.user_index(); - reco::CandidatePtr const& itow = inputs_[index]; + reco::Candidate const* itow = inputs_[index]; - int it = ieta(itow); - iphi(itow); + const CaloTower* ctc = dynamic_cast(itow); + int it = ctc->id().ieta(); - double Original_Et = itow->et(); - double etnew = Original_Et - (emean_.find(it))->second - (esigma_.find(it))->second; + double original_Et = itow->et(); + double etnew = original_Et - (emean_.find(it))->second - (esigma_.find(it))->second; float mScale = etnew / input_object.Et(); if (etnew < 0.) mScale = 0.; @@ -482,8 +469,8 @@ void HiPuRhoProducer::calculateOrphanInput(std::vector& orph std::vector> excludedTowers; // vector of excluded ieta, iphi values int32_t nref = 0; - for (auto const& pseudojetTMP : fjJets_) { + EtaPhi jet_etaphi(pseudojetTMP.eta(), pseudojetTMP.phi()); if (nref < nMaxJets_) { jtexngeom[nref] = 0; jtexntow[nref] = 0; @@ -497,13 +484,13 @@ void HiPuRhoProducer::calculateOrphanInput(std::vector& orph continue; for (auto const& im : allgeomid_) { - double dr = reco::deltaR(geo_->getPosition((DetId)im), pseudojetTMP); - std::vector>::const_iterator exclude = - std::find(excludedTowers.begin(), excludedTowers.end(), std::pair(im.ieta(), im.iphi())); - if (dr < radiusPU_ && exclude == excludedTowers.end() && + double dr2 = + reco::deltaR2(towermap[(DetId)im].first, towermap[(DetId)im].second, jet_etaphi.first, jet_etaphi.second); + auto exclude = std::find(excludedTowers.begin(), excludedTowers.end(), std::pair(im.ieta(), im.iphi())); + if (dr2 < radiusPU_ * radiusPU_ && exclude == excludedTowers.end() && (geomtowers_[im.ieta()] - ntowersWithJets_[im.ieta()]) > minimumTowersFraction_ * (geomtowers_[im.ieta()])) { ntowersWithJets_[im.ieta()]++; - excludedTowers.emplace_back(std::pair(im.ieta(), im.iphi())); + excludedTowers.emplace_back(im.ieta(), im.iphi()); if (nref < nMaxJets_) jtexngeom[nref]++; @@ -512,17 +499,17 @@ void HiPuRhoProducer::calculateOrphanInput(std::vector& orph for (auto const& input : fjInputs_) { int index = input.user_index(); - const reco::CandidatePtr& originalTower = inputs_[index]; - - int ie = ieta(originalTower); - int ip = iphi(originalTower); + const reco::Candidate* originalTower = inputs_[index]; + const CaloTower* ctc = dynamic_cast(originalTower); + int ie = ctc->id().ieta(); + int ip = ctc->id().iphi(); auto exclude = std::find(excludedTowers.begin(), excludedTowers.end(), std::pair(ie, ip)); if (exclude != excludedTowers.end()) { jettowers.push_back(index); } - double dr = reco::deltaR(input, pseudojetTMP); - if (dr < radiusPU_ && nref < nMaxJets_) { + double dr2 = reco::deltaR2(input.eta(), input.phi(), jet_etaphi.first, jet_etaphi.second); + if (dr2 < radiusPU_ * radiusPU_ && nref < nMaxJets_) { jtexntow[nref]++; jtexpt[nref] += originalTower->pt(); } @@ -531,11 +518,10 @@ void HiPuRhoProducer::calculateOrphanInput(std::vector& orph if (nref < nMaxJets_) nref++; } - // Create a new collections from the towers not included in jets for (auto const& input : fjInputs_) { int index = input.user_index(); - const reco::CandidatePtr& originalTower = inputs_[index]; + const reco::Candidate* originalTower = inputs_[index]; auto itjet = std::find(jettowers.begin(), jettowers.end(), index); if (itjet == jettowers.end()) { orphanInput.emplace_back(originalTower->px(), originalTower->py(), originalTower->pz(), originalTower->energy()); @@ -555,7 +541,7 @@ void HiPuRhoProducer::putRho(edm::Event& iEvent, const edm::EventSetup& iSetup) std::vector> order; for (std::size_t i = 0; i < size; ++i) { - order.emplace_back(std::make_pair(i, etaEdgeLow_[i])); + order.emplace_back(i, etaEdgeLow_[i]); } std::sort( @@ -612,33 +598,17 @@ void HiPuRhoProducer::putRho(edm::Event& iEvent, const edm::EventSetup& iSetup) // ------------ method fills 'descriptions' with the allowed parameters for the module ------------ void HiPuRhoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - //The following says we do not know what parameters are allowed so do no validation - // Please change this to state exactly what you do use, even if it is no parameters edm::ParameterSetDescription desc; - desc.setUnknown(); - descriptions.addDefault(desc); -} - -int HiPuRhoProducer::ieta(const reco::CandidatePtr& in) const { - const CaloTower* ctc = dynamic_cast(in.get()); - if (ctc) { - return ctc->id().ieta(); - } - - throw cms::Exception("Invalid Constituent") << "CaloJet constituent is not of CaloTower type"; - - return 0; -} - -int HiPuRhoProducer::iphi(const reco::CandidatePtr& in) const { - const CaloTower* ctc = dynamic_cast(in.get()); - if (ctc) { - return ctc->id().iphi(); - } - - throw cms::Exception("Invalid Constituent") << "CaloJet constituent is not of CaloTower type"; - - return 0; + desc.add("src", edm::InputTag("PFTowers")); + desc.add("medianWindowWidth", 2); + desc.add("towSigmaCut", 5.0); + desc.add("puPtMin", 15.0); + desc.add("rParam", 0.3); + desc.add("nSigmaPU", 1.0); + desc.add("radiusPU", 0.5); + desc.add("minimumTowersFraction", 0.7); + desc.add("dropZeroTowers", true); + descriptions.add("hiPuRhoProducer", desc); } //define this as a plug-in diff --git a/RecoHI/HiJetAlgos/python/HiRecoPFJets_cff.py b/RecoHI/HiJetAlgos/python/HiRecoPFJets_cff.py index 1deefa36f0de9..4ed7fc4151a6b 100644 --- a/RecoHI/HiJetAlgos/python/HiRecoPFJets_cff.py +++ b/RecoHI/HiJetAlgos/python/HiRecoPFJets_cff.py @@ -45,17 +45,33 @@ akPu6PFJets = akPu5PFJets.clone(rParam = cms.double(0.6), puPtMin = 30) akPu7PFJets = akPu5PFJets.clone(rParam = cms.double(0.7), puPtMin = 35) +hiPFCandCleanerforJets = cms.EDFilter('GenericPFCandidateSelector', + src = cms.InputTag('particleFlow'), + cut = cms.string("pt>5 && abs(eta)< 2") + ) + ak4PFJetsForFlow = akPu5PFJets.clone( - Ghost_EtaMax = cms.double(5.0), - Rho_EtaMax = cms.double(4.4), - doRhoFastjet = cms.bool(False), - jetPtMin = cms.double(15.0), - nSigmaPU = cms.double(1.0), - rParam = cms.double(0.4), - radiusPU = cms.double(0.5), - src = cms.InputTag("pfcandCleaner", "particleFlowCleaned"), + Ghost_EtaMax = 5.0, + Rho_EtaMax = 4.4, + doRhoFastjet = False, + jetPtMin = 15.0, + nSigmaPU = cms.double(1.0), + rParam = 0.4, + radiusPU = cms.double(0.5), + src = "hiPFCandCleanerforJets", ) +# ak4PFJetsForFlow = akPu5PFJets.clone( +# Ghost_EtaMax = cms.double(5.0), +# Rho_EtaMax = cms.double(4.4), +# doRhoFastjet = cms.bool(False), +# jetPtMin = cms.double(15.0), +# nSigmaPU = cms.double(1.0), +# rParam = cms.double(0.4), +# radiusPU = cms.double(0.5), +# src = cms.InputTag("hiPFCandCleanerforJets"), +# ) + kt4PFJetsForRho = cms.EDProducer( "FastjetJetProducer", HiPFJetParameters, @@ -70,9 +86,36 @@ kt4PFJetsForRho.GhostArea = cms.double(0.005) from RecoHI.HiJetAlgos.hiFJRhoProducer import hiFJRhoProducer -from RecoHI.HiJetAlgos.hiPuRhoProducer_cfi import hiPuRhoProducer -from RecoHI.HiJetAlgos.hiFJRhoFlowModulationProducer_cfi import hiFJRhoFlowModulationProducer -from RecoHI.HiJetAlgos.hiPFCandCleaner_cfi import hiPFCandCleaner +#from RecoHI.HiJetAlgos.hiPuRhoProducer_cfi import hiPuRhoProducer +#from RecoHI.HiJetAlgos.hiFJRhoFlowModulationProducer_cfi import hiFJRhoFlowModulationProducer +#from RecoHI.HiJetAlgos.hiPuRhoProducer import hiPuRhoProducer +#from RecoHI.HiJetAlgos.hiFJRhoFlowModulationProducer import hiFJRhoFlowModulationProducer +#from RecoHI.HiJetAlgos.hiPFCandCleaner_cfi import hiPFCandCleaner + +hiFJRhoFlowModulationProducer = cms.EDProducer( + 'HiFJRhoFlowModulationProducer', + pfCandSource = cms.InputTag('particleFlow'), + doJettyExclusion = cms.bool(False), + doFreePlaneFit = cms.bool(False), + doEvtPlane = cms.bool(False), + doFlatTest = cms.bool(False), + jetTag = cms.InputTag("ak4PFJets"), + EvtPlane = cms.InputTag("hiEvtPlane"), + evtPlaneLevel = cms.int32(0) + ) + +hiPuRhoProducer = cms.EDProducer( + 'HiPuRhoProducer', + dropZeroTowers = cms.bool(True), + medianWindowWidth = cms.int32(2), + minimumTowersFraction = cms.double(0.7), + nSigmaPU = cms.double(1.0), + puPtMin = cms.double(15.0), + rParam = cms.double(.3), + radiusPU = cms.double(.5), + src = cms.InputTag('PFTowers'), + towSigmaCut = cms.double(5.), + ) akCs4PFJets = cms.EDProducer( "CSJetProducer", @@ -103,7 +146,7 @@ akPu3PFJets, akPu4PFJets, akPu5PFJets, - hiPFCandCleaner, + hiPFCandCleanerforJets, kt4PFJetsForRho, ak4PFJetsForFlow, hiFJRhoProducer, diff --git a/RecoJets/Configuration/python/RecoPFJets_cff.py b/RecoJets/Configuration/python/RecoPFJets_cff.py index ba299241eb415..138f0471afb29 100644 --- a/RecoJets/Configuration/python/RecoPFJets_cff.py +++ b/RecoJets/Configuration/python/RecoPFJets_cff.py @@ -84,10 +84,10 @@ ) recoPFJetsWithSubstructure=cms.Sequence(recoPFJetsWithSubstructureTask) -from RecoHI.HiJetAlgos.HiRecoPFJets_cff import PFTowers, akPu3PFJets, akPu4PFJets, kt4PFJetsForRho, ak4PFJetsForFlow, akCs4PFJets, pfNoPileUpJMEHI +from RecoHI.HiJetAlgos.HiRecoPFJets_cff import PFTowers, akPu3PFJets, akPu4PFJets, kt4PFJetsForRho, ak4PFJetsForFlow, akCs4PFJets, pfNoPileUpJMEHI, hiFJRhoFlowModulationProducer, hiPuRhoProducer, hiPFCandCleanerforJets from RecoHI.HiJetAlgos.hiFJRhoProducer import hiFJRhoProducer -from RecoHI.HiJetAlgos.hiPuRhoProducer_cfi import hiPuRhoProducer -from RecoHI.HiJetAlgos.hiFJRhoFlowModulationProducer_cfi import hiFJRhoFlowModulationProducer + + recoPFJetsHITask =cms.Task(fixedGridRhoAll, fixedGridRhoFastjetAll, fixedGridRhoFastjetCentral, @@ -100,6 +100,7 @@ akPu3PFJets, akPu4PFJets, kt4PFJetsForRho, + hiPFCandCleanerforJets, ak4PFJetsForFlow, hiFJRhoProducer, hiPuRhoProducer, diff --git a/RecoJets/JetProducers/plugins/CSJetProducer.cc b/RecoJets/JetProducers/plugins/CSJetProducer.cc index 272cb63830745..5db800d412300 100644 --- a/RecoJets/JetProducers/plugins/CSJetProducer.cc +++ b/RecoJets/JetProducers/plugins/CSJetProducer.cc @@ -16,16 +16,18 @@ using namespace edm; using namespace cms; CSJetProducer::CSJetProducer(edm::ParameterSet const& conf) - : VirtualJetProducer(conf), csRParam_(-1.0), csAlpha_(0.), useModulatedRho_(false) { + : VirtualJetProducer(conf), + csRParam_(-1.0), + csAlpha_(0.), + useModulatedRho_(conf.getParameter("useModulatedRho")), + minFlowChi2Prob_(conf.getParameter("minFlowChi2Prob")), + maxFlowChi2Prob_(conf.getParameter("maxFlowChi2Prob")) { //get eta range, rho and rhom map etaToken_ = consumes>(conf.getParameter("etaMap")); rhoToken_ = consumes>(conf.getParameter("rho")); rhomToken_ = consumes>(conf.getParameter("rhom")); csRParam_ = conf.getParameter("csRParam"); csAlpha_ = conf.getParameter("csAlpha"); - useModulatedRho_ = conf.getParameter("useModulatedRho"); - minFlowChi2Prob_ = conf.getParameter("minFlowChi2Prob"); - maxFlowChi2Prob_ = conf.getParameter("maxFlowChi2Prob"); if (useModulatedRho_) rhoFlowFitParamsToken_ = consumes>(conf.getParameter("rhoFlowFitParams")); } @@ -60,7 +62,6 @@ void CSJetProducer::runAlgorithm(edm::Event& iEvent, edm::EventSetup const& iSet edm::Handle> rhoRanges; edm::Handle> rhomRanges; edm::Handle> rhoFlowFitParams; - iEvent.getByToken(etaToken_, etaRanges); iEvent.getByToken(rhoToken_, rhoRanges); iEvent.getByToken(rhomToken_, rhomRanges); @@ -77,6 +78,15 @@ void CSJetProducer::runAlgorithm(edm::Event& iEvent, edm::EventSetup const& iSet << "Size of etaRanges (" << bkgVecSize << ") and rhoRanges (" << rhoRanges->size() << ") and/or rhomRanges (" << rhomRanges->size() << ") vectors inconsistent\n"; } + bool minProb = false; + bool maxProb = false; + if (useModulatedRho_) { + if (!rhoFlowFitParams->empty()) { + double val = ROOT::Math::chisquared_cdf_c(rhoFlowFitParams->at(5), rhoFlowFitParams->at(6)); + minProb = val > minFlowChi2Prob_; + maxProb = val < maxFlowChi2Prob_; + } + } //Allow the background densities to change within the jet @@ -96,10 +106,6 @@ void CSJetProducer::runAlgorithm(edm::Event& iEvent, edm::EventSetup const& iSet if (useModulatedRho_) { if (!rhoFlowFitParams->empty()) { - double val = ROOT::Math::chisquared_cdf_c(rhoFlowFitParams->at(5), rhoFlowFitParams->at(6)); - bool minProb = val > minFlowChi2Prob_; - bool maxProb = val < maxFlowChi2Prob_; - if (minProb && maxProb) { rhoModulationFactor = getModulatedRhoFactor(ghostPhi, rhoFlowFitParams->at(2), @@ -171,10 +177,10 @@ void CSJetProducer::fillDescriptionsFromCSJetProducer(edm::ParameterSetDescripti desc.add("useModulatedRho", false); desc.add("minFlowChi2Prob", false); desc.add("maxFlowChi2Prob", false); - desc.add("etaMap", edm::InputTag("hiFJRhoProducer", "mapEtaEdges")); - desc.add("rho", edm::InputTag("hiFJRhoProducer", "mapToRho")); - desc.add("rhom", edm::InputTag("hiFJRhoProducer", "mapToRhoM")); - desc.add("rhoFlowFitParams", edm::InputTag("hiFJRhoFlowModulationProducer", "rhoFlowFitParams")); + desc.add("etaMap", {"hiFJRhoProducer", "mapEtaEdges"}); + desc.add("rho", {"hiFJRhoProducer", "mapToRho"}); + desc.add("rhom", {"hiFJRhoProducer", "mapToRhoM"}); + desc.add("rhoFlowFitParams", {"hiFJRhoFlowModulationProducer", "rhoFlowFitParams"}); } double CSJetProducer::getModulatedRhoFactor( diff --git a/RecoJets/JetProducers/plugins/CSJetProducer.h b/RecoJets/JetProducers/plugins/CSJetProducer.h index def133514d22c..308d6cb135c77 100644 --- a/RecoJets/JetProducers/plugins/CSJetProducer.h +++ b/RecoJets/JetProducers/plugins/CSJetProducer.h @@ -34,8 +34,7 @@ namespace cms { protected: void runAlgorithm(edm::Event& iEvent, const edm::EventSetup& iSetup) override; - double getModulatedRhoFactor( - const double phi, const double eventPlane2, const double eventPlane3, const double par1, const double par2); + double getModulatedRhoFactor(double phi, double eventPlane2, double eventPlane3, double par1, double par2); double csRParam_; /// for constituent subtraction : R parameter double csAlpha_; /// for HI constituent subtraction : alpha (power of pt in metric) @@ -43,7 +42,6 @@ namespace cms { bool useModulatedRho_; /// flag to turn on/off flow-modulated rho and rhom double minFlowChi2Prob_; /// flowFit chi2/ndof minimum compatability requirement double maxFlowChi2Prob_; /// flowFit chi2/ndof minimum compatability requirement - //input rho and rho_m + eta map edm::EDGetTokenT> etaToken_; edm::EDGetTokenT> rhoToken_;