From 2c00969b61552be12ec47d62458f40eeeaa70505 Mon Sep 17 00:00:00 2001 From: Quan Wang Date: Fri, 18 Sep 2020 21:16:06 +0200 Subject: [PATCH 1/5] [+] update RP for mini-AOD --- DQM/Physics/src/CentralityDQM.cc | 13 - DQM/Physics/src/CentralityDQM.h | 7 - DQM/Physics/src/CentralitypADQM.h | 17 - .../python/slimming/MicroEventContent_cff.py | 8 +- .../PatAlgos/python/slimming/slimming_cff.py | 6 +- RecoHI/HiEvtPlaneAlgos/BuildFile.xml | 6 + RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h | 166 +++++ .../interface/HiEvtPlaneFlatten.h | 59 +- .../interface/HiEvtPlaneList.h | 185 +++--- RecoHI/HiEvtPlaneAlgos/interface/LoadEPDB.h | 6 +- .../HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py | 48 +- .../python/RecoHiEvtPlane_EventContent_cff.py | 29 +- .../python/hiEvtPlaneFlat_cfi.py | 21 +- .../HiEvtPlaneAlgos/src/EvtPlaneProducer.cc | 579 ++++++++++++------ .../src/HiEvtPlaneFlatProducer.cc | 167 +++-- 15 files changed, 872 insertions(+), 445 deletions(-) create mode 100644 RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h diff --git a/DQM/Physics/src/CentralityDQM.cc b/DQM/Physics/src/CentralityDQM.cc index e8c2ad6889a22..f49b9d8df38d0 100644 --- a/DQM/Physics/src/CentralityDQM.cc +++ b/DQM/Physics/src/CentralityDQM.cc @@ -94,17 +94,11 @@ void CentralityDQM::bookHistograms(DQMStore::IBooker& bei, edm::Run const&, edm: Double_t psirange = 4; bei.setCurrentFolder("Physics/Centrality/EventPlane/"); - h_ep_HFm1 = bei.book1D("h_ep_HFm1", "h_ep_HFm1", 800, -psirange, psirange); - h_ep_HFp1 = bei.book1D("h_ep_HFp1", "h_ep_HFp1", 800, -psirange, psirange); - h_ep_trackm1 = bei.book1D("h_ep_trackm1", "h_ep_trackm1", 800, -psirange, psirange); - h_ep_trackp1 = bei.book1D("h_ep_trackp1", "h_ep_trackp1", 800, -psirange, psirange); - h_ep_castor1 = bei.book1D("h_ep_castor1", "h_ep_castor1", 800, -psirange, psirange); h_ep_HFm2 = bei.book1D("h_ep_HFm2", "h_ep_HFm2", 800, -psirange, psirange); h_ep_HFp2 = bei.book1D("h_ep_HFp2", "h_ep_HFp2", 800, -psirange, psirange); h_ep_trackmid2 = bei.book1D("h_ep_trackmid2", "h_ep_trackmid2", 800, -psirange, psirange); h_ep_trackm2 = bei.book1D("h_ep_trackm2", "h_ep_trackm2", 800, -psirange, psirange); h_ep_trackp2 = bei.book1D("h_ep_trackp2", "h_ep_trackp2", 800, -psirange, psirange); - h_ep_castor2 = bei.book1D("h_ep_castor2", "h_ep_castor2", 800, -psirange, psirange); h_ep_HFm3 = bei.book1D("h_ep_HFm3", "h_ep_HFm3", 800, -psirange, psirange); h_ep_HFp3 = bei.book1D("h_ep_HFp3", "h_ep_HFp3", 800, -psirange, psirange); h_ep_trackmid3 = bei.book1D("h_ep_trackmid3", "h_ep_trackmid3", 800, -psirange, psirange); @@ -172,18 +166,11 @@ void CentralityDQM::analyze(const edm::Event& iEvent, const edm::EventSetup& iSe if (ep.isValid()) { EvtPlaneCollection::const_iterator rp = ep->begin(); - h_ep_HFm1->Fill(rp[HFm1].angle(0)); - h_ep_HFp1->Fill(rp[HFp1].angle(0)); - h_ep_trackm1->Fill(rp[trackm1].angle(0)); - h_ep_trackp1->Fill(rp[trackp1].angle(0)); - h_ep_castor1->Fill(rp[Castor1].angle(0)); - h_ep_HFm2->Fill(rp[HFm2].angle(0)); h_ep_HFp2->Fill(rp[HFp2].angle(0)); h_ep_trackmid2->Fill(rp[trackmid2].angle(0)); h_ep_trackm2->Fill(rp[trackm2].angle(0)); h_ep_trackp2->Fill(rp[trackp2].angle(0)); - h_ep_castor2->Fill(rp[Castor2].angle(0)); h_ep_HFm3->Fill(rp[HFm3].angle(0)); h_ep_HFp3->Fill(rp[HFp3].angle(0)); diff --git a/DQM/Physics/src/CentralityDQM.h b/DQM/Physics/src/CentralityDQM.h index ba0a730ef56ed..02a4c90a6dd03 100644 --- a/DQM/Physics/src/CentralityDQM.h +++ b/DQM/Physics/src/CentralityDQM.h @@ -83,18 +83,11 @@ class CentralityDQM : public DQMEDAnalyzer { MonitorElement* h_cent_bin; - MonitorElement* h_ep_HFm1; - MonitorElement* h_ep_HFp1; - MonitorElement* h_ep_trackm1; - MonitorElement* h_ep_trackp1; - MonitorElement* h_ep_castor1; - MonitorElement* h_ep_HFm2; MonitorElement* h_ep_HFp2; MonitorElement* h_ep_trackmid2; MonitorElement* h_ep_trackm2; MonitorElement* h_ep_trackp2; - MonitorElement* h_ep_castor2; MonitorElement* h_ep_HFm3; MonitorElement* h_ep_HFp3; diff --git a/DQM/Physics/src/CentralitypADQM.h b/DQM/Physics/src/CentralitypADQM.h index 979d3e7b7c505..0dfd1fd7aba26 100644 --- a/DQM/Physics/src/CentralitypADQM.h +++ b/DQM/Physics/src/CentralitypADQM.h @@ -79,23 +79,6 @@ class CentralitypADQM : public DQMEDAnalyzer { MonitorElement* h_vertex_z; MonitorElement* h_cent_bin; - - MonitorElement* h_ep_HFm1; - MonitorElement* h_ep_HFp1; - MonitorElement* h_ep_trackm1; - MonitorElement* h_ep_trackp1; - MonitorElement* h_ep_castor1; - - MonitorElement* h_ep_HFm2; - MonitorElement* h_ep_HFp2; - MonitorElement* h_ep_trackmid2; - MonitorElement* h_ep_trackm2; - MonitorElement* h_ep_trackp2; - MonitorElement* h_ep_castor2; - - MonitorElement* h_ep_HFm3; - MonitorElement* h_ep_HFp3; - MonitorElement* h_ep_trackmid3; }; #endif diff --git a/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py b/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py index db6e283aa8a3f..a5e86d3f8cca5 100644 --- a/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py +++ b/PhysicsTools/PatAlgos/python/slimming/MicroEventContent_cff.py @@ -66,7 +66,7 @@ 'keep *_l1extraParticles_*_*', 'keep L1GlobalTriggerReadoutRecord_gtDigis_*_*', # stage 2 L1 trigger - 'keep *_gtStage2Digis__*', + 'keep *_gtStage2Digis__*', 'keep *_gmtStage2Digis_Muon_*', 'keep *_caloStage2Digis_Jet_*', 'keep *_caloStage2Digis_Tau_*', @@ -100,7 +100,7 @@ 'keep GenLumiInfoHeader_generator_*_*', 'keep GenLumiInfoProduct_*_*_*', 'keep GenEventInfoProduct_generator_*_*', - 'keep recoGenParticles_genPUProtons_*_*', + 'keep recoGenParticles_genPUProtons_*_*', 'keep *_slimmedGenJetsFlavourInfos_*_*', 'keep *_slimmedGenJets__*', 'keep *_slimmedGenJetsAK8__*', @@ -136,7 +136,7 @@ _run3_common_extraCommands = ["drop *_packedPFCandidates_hcalDepthEnergyFractions_*"] from Configuration.Eras.Modifier_run3_common_cff import run3_common run3_common.toModify(MicroEventContent, outputCommands = MicroEventContent.outputCommands + _run3_common_extraCommands) -# --- +# --- _pp_on_AA_extraCommands = [ 'keep patPackedCandidates_hiPixelTracks_*_*', @@ -147,6 +147,8 @@ 'keep recoCentrality_hiCentrality_*_*', 'keep int_centralityBin_*_*', 'keep recoHFFilterInfo_hiHFfilters_*_*', + 'keep *_hiEvtPlane_*_*', + 'keep *_hiEvtPlaneFlat_*_*' ] 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 diff --git a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py index ca2b48f515bf1..8c1e497f0bf92 100644 --- a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py +++ b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py @@ -67,7 +67,11 @@ pp_on_AA_2018.toReplaceWith(slimmingTask, slimmingTask.copyAndExclude([slimmedOOTPhotons])) from Configuration.Eras.Modifier_pp_on_PbPb_run3_cff import pp_on_PbPb_run3 from PhysicsTools.PatAlgos.slimming.hiPixelTracks_cfi import hiPixelTracks -(pp_on_AA_2018 | pp_on_PbPb_run3).toReplaceWith(slimmingTask, cms.Task(slimmingTask.copy(), hiPixelTracks)) +from RecoHI.HiCentralityAlgos.CentralityBin_cfi import centralityBin +from RecoHI.HiEvtPlaneAlgos.HiEvtPlane_cfi import hiEvtPlane +from RecoHI.HiEvtPlaneAlgos.hiEvtPlaneFlat_cfi import hiEvtPlaneFlat +(pp_on_AA_2018 | pp_on_PbPb_run3).toReplaceWith(slimmingTask, cms.Task(slimmingTask.copy(), hiPixelTracks, +centralityBin, hiEvtPlane, hiEvtPlaneFlat)) from Configuration.Eras.Modifier_pp_on_PbPb_run3_cff import pp_on_PbPb_run3 from PhysicsTools.PatAlgos.packedCandidateMuonID_cfi import packedCandidateMuonID diff --git a/RecoHI/HiEvtPlaneAlgos/BuildFile.xml b/RecoHI/HiEvtPlaneAlgos/BuildFile.xml index 0900da4ec6f0e..e64ad6baf50de 100644 --- a/RecoHI/HiEvtPlaneAlgos/BuildFile.xml +++ b/RecoHI/HiEvtPlaneAlgos/BuildFile.xml @@ -1,14 +1,20 @@ + + + + + + diff --git a/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h b/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h new file mode 100644 index 0000000000000..7445ef49f4483 --- /dev/null +++ b/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h @@ -0,0 +1,166 @@ +#ifndef __EPCuts__ +#define __EPCuts__ + +namespace hi { + + enum class EP_ERA { ppReco, HIReco, Pixel, GenMC }; + + struct TrackStructure { + int centbin; + double eta; + double phi; + double et; + double pt; + int charge; + int pdgid; + int hits; + int algos; + int collection; + double dz; + double dxy; + double dzError; + double dxyError; + double ptError; + bool highPurity; + double dzSig; + double dxySig; + double normalizedChi2; + double dzError_Pix; + double chi2layer; + int numberOfValidHits; + int pixel; + }; + + class EPCuts { + private: + EP_ERA cutera_; + double pterror_; + double dzerror_; + double dxyerror_; + double chi2perlayer_; + double dzerror_Pix_; + double chi2Pix_; + int numberOfValidHits_; + + public: + explicit EPCuts(EP_ERA cutEra = EP_ERA::ppReco, + double pterror = 0.1, + double dzerror = 3.0, + double dxyerror = 3.0, + double chi2perlayer = 0.18, + double dzError_Pix = 10.0, + double chi2Pix = 40., + int numberOfValidHits = 11) { + cutera_ = cutEra; + pterror_ = pterror; + dzerror_ = dzerror; + dxyerror_ = dxyerror; + chi2perlayer_ = chi2perlayer; + dzerror_Pix_ = dzError_Pix; + chi2Pix_ = chi2Pix; + numberOfValidHits_ = numberOfValidHits; + } + ~EPCuts() {} + + bool isGoodHF(TrackStructure track) { + if (track.pdgid != 1 && track.pdgid != 2) + return false; + if (fabs(track.eta) < 3 || fabs(track.eta) > 5) + return false; + return true; + } + + bool isGoodCastor(TrackStructure track) { return true; } + + bool isGoodTrack(TrackStructure track) { + if (cutera_ == EP_ERA::ppReco) + return TrackQuality_ppReco(track); + if (cutera_ == EP_ERA::HIReco) + return TrackQuality_HIReco(track); + if (cutera_ == EP_ERA::Pixel) + return TrackQuality_Pixel(track); + return false; + } + + bool TrackQuality_ppReco(TrackStructure track) { + if (track.charge == 0) + return false; + if (!track.highPurity) + return false; + if (track.ptError / track.pt > pterror_) + return false; + if (track.numberOfValidHits < numberOfValidHits_) + return false; + if (track.chi2layer > chi2perlayer_) + return false; + if (fabs(track.dxy / track.dxyError) > dxyerror_) + return false; + if (fabs(track.dz / track.dzError) > dzerror_) + return false; + return true; + } + + bool TrackQuality_HIReco(TrackStructure track) { + if (track.charge == 0) + return false; + if (!track.highPurity) + return false; + if (track.numberOfValidHits < numberOfValidHits_) + return false; + if (track.ptError / track.pt > pterror_) + return false; + if (fabs(track.dxy / track.dxyError) > dxyerror_) + return false; + if (fabs(track.dz / track.dzError) > dzerror_) + return false; + if (track.chi2layer > chi2perlayer_) + return false; + if (track.algos != 4 && track.algos != 5 && track.algos != 6 && track.algos != 7) + return false; + return true; + } + + bool TrackQuality_Pixel(TrackStructure track) { + if (track.charge == 0) + return false; + if (!track.highPurity) + return false; + bool bPix = false; + int nHits = track.numberOfValidHits; + if (track.ptError / track.pt > pterror_) + return false; + if (track.pt < 2.4 and (nHits == 3 or nHits == 4 or nHits == 5 or nHits == 6)) + bPix = true; + if (not bPix) { + if (nHits < numberOfValidHits_) + return false; + if (track.chi2layer > chi2perlayer_) + return false; + if (track.ptError / track.pt > pterror_) + return false; + int algo = track.algos; + if (track.pt > 2.4 && algo != 4 && algo != 5 && algo != 6 && algo != 7) + return false; + if (fabs(track.dxy / track.dxyError) > dxyerror_) + return false; + if (fabs(track.dz / track.dzError) > dzerror_) + return false; + } else { + if (track.chi2layer > chi2Pix_) + return false; + if (fabs(track.dz / track.dzError) > dzerror_Pix_) + return false; + } + return true; + } + + bool TrackQuality_GenMC(TrackStructure track) { + if (track.charge == 0) + return false; + if (fabs(track.eta) > 2.4) + return false; + return true; + } + }; +} // namespace hi +#endif diff --git a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h index 0e27e8a0b8342..4dd857e297c6f 100644 --- a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h +++ b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h @@ -32,9 +32,20 @@ class HiEvtPlaneFlatten { vorder_ = 2; //sets default order of event plane } - void init(int order, int nbins, std::string tag, int vord) { + void init(int order, + int nbins, + int nvtxbins = 10, + double minvtx = -25, + double delvtx = 5, + std::string tag = "", + int vord = 2) { hOrder_ = order; //order of flattening vorder_ = vord; //1(v1), 2(v2), 3(v3), 4(v4) + nvtxbins_ = nvtxbins; + nbins_ = nbins; + minvtx_ = minvtx; + delvtx_ = delvtx; + caloCentRefMinBin_ = -1; caloCentRefMaxBin_ = -1; hbins_ = nbins * nvtxbins_ * hOrder_; @@ -152,8 +163,9 @@ class HiEvtPlaneFlatten { double scale = getEtScale(vtx, centbin); double ptval = getPtDB(indx) * scale; double pt2val = getPt2DB(indx) * pow(scale, 2); + double rv = pt * scale - pt2val / ptval; if (ptval > 0) - return pt * scale - pt2val / ptval; + return rv; } return 0.; } @@ -174,10 +186,11 @@ class HiEvtPlaneFlatten { double getSoffset(double s, double vtx, int centbin) const { int indx = getOffsetIndx(centbin, vtx); - if (indx >= 0) + if (indx >= 0) { return s - yoffDB_[indx]; - else + } else { return s; + } } double getCoffset(double c, double vtx, int centbin) const { @@ -197,6 +210,37 @@ class HiEvtPlaneFlatten { return psi; } + float getMinCent(int indx) { + int ibin = (int)(indx / nvtxbins_); + return ibin * 100. / nbins_; + } + + float getMaxCent(int indx) { + int ibin = (int)(indx / nvtxbins_); + return (ibin + 1) * 100. / nbins_; + } + + double getMinVtx(int indx) { + int ivtx = indx - nvtxbins_ * (int)(indx / nvtxbins_); + return minvtx_ + ivtx * delvtx_; + } + + double getMaxVtx(int indx) { + int ivtx = indx - nvtxbins_ * (int)(indx / nvtxbins_); + return minvtx_ + (ivtx + 1) * delvtx_; + } + + std::string getRangeString(int indx) { + char buf[120]; + sprintf(buf, + "%5.1f < cent < %5.1f; %4.1f < vtx < %4.1f", + getMinCent(indx), + getMaxCent(indx), + getMinVtx(indx), + getMaxVtx(indx)); + return std::string(buf); + } + ~HiEvtPlaneFlatten() {} int getHBins() const { return hbins_; } int getOBins() const { return obins_; } @@ -421,9 +465,6 @@ class HiEvtPlaneFlatten { } private: - static constexpr int nvtxbins_ = 10; - static constexpr double minvtx_ = -25.; - static constexpr double delvtx_ = 5.; static const int MAXCUT = 10000; static const int MAXCUTOFF = 1000; @@ -470,6 +511,10 @@ class HiEvtPlaneFlatten { double centRes40_[2]; double centResErr40_[2]; + int nvtxbins_; + int nbins_; + double minvtx_; + double delvtx_; int hOrder_; //flattening order int hbins_; //number of bins needed for flattening int obins_; //number of (x,y) offset bins diff --git a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h index 784f708198d3f..d3d10a6788088 100644 --- a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h +++ b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h @@ -2,127 +2,110 @@ #define __HiEvtPlaneList__ /* Index Name Detector Order hmin1 hmax1 hmin2 hmax2 minpt maxpt nsub mcw rmate1 rmate2 - 0 HFm1 HF 1 -5.00 -3.00 0.00 0.00 0.01 30.00 3sub no HFp1 trackp1 - 1 HFp1 HF 1 3.00 5.00 0.00 0.00 0.01 30.00 3sub no HFm1 trackm1 - 2 HF1 HF 1 -5.00 -3.00 3.00 5.00 0.01 30.00 3sub no trackm1 trackp1 - 3 trackm1 Tracker 1 -2.00 -1.00 0.00 0.00 0.30 3.00 3sub no HFm1 HFp1 - 4 trackp1 Tracker 1 1.00 2.00 0.00 0.00 0.30 3.00 3sub no HFm1 HFp1 - 5 Castor1 Castor 1 -6.55 -5.10 0.00 0.00 0.01 50.00 3sub no HFp1 trackp1 - 6 HFm2 HF 2 -5.00 -3.00 0.00 0.00 0.01 30.00 3sub no HFp2 trackmid2 - 7 HFp2 HF 2 3.00 5.00 0.00 0.00 0.01 30.00 3sub no HFm2 trackmid2 - 8 HF2 HF 2 -5.00 -3.00 3.00 5.00 0.01 30.00 3sub no trackm2 trackp2 - 9 trackmid2 Tracker 2 -0.75 0.75 0.00 0.00 0.30 3.00 3sub no HFm2 HFp2 - 10 trackm2 Tracker 2 -2.00 -1.00 0.00 0.00 0.30 3.00 3sub no HFm2 HFp2 - 11 trackp2 Tracker 2 1.00 2.00 0.00 0.00 0.30 3.00 3sub no HFm2 HFp2 - 12 Castor2 Castor 2 -6.55 -5.10 0.00 0.00 0.01 50.00 3sub no trackmid2 HFp2 - 13 HFm3 HF 3 -5.00 -3.00 0.00 0.00 0.01 30.00 3sub no HFp3 trackmid3 - 14 HFp3 HF 3 3.00 5.00 0.00 0.00 0.01 30.00 3sub no HFm3 trackmid3 - 15 HF3 HF 3 -5.00 -3.00 3.00 5.00 0.01 30.00 3sub no trackm3 trackp3 - 16 trackmid3 Tracker 3 -0.75 0.75 0.00 0.00 0.30 3.00 3sub no HFm3 HFp3 - 17 trackm3 Tracker 3 -2.00 -1.00 0.00 0.00 0.30 3.00 3sub no HFm3 HFp3 - 18 trackp3 Tracker 3 1.00 2.00 0.00 0.00 0.30 3.00 3sub no HFm3 HFp3 - 19 HFm4 HF 4 -5.00 -3.00 0.00 0.00 0.01 30.00 3sub no HFp4 trackmid4 - 20 HFp4 HF 4 3.00 5.00 0.00 0.00 0.01 30.00 3sub no HFm4 trackmid4 - 21 HF4 HF 4 -5.00 -3.00 3.00 5.00 0.01 30.00 3sub no trackm4 trackp4 - 22 trackmid4 Tracker 4 -0.75 0.75 0.00 0.00 0.30 3.00 3sub no HFm4 HFp4 - 23 trackm4 Tracker 4 -2.00 -1.00 0.00 0.00 0.30 3.00 3sub no HFm4 HFp4 - 24 trackp4 Tracker 4 1.00 2.00 0.00 0.00 0.30 3.00 3sub no HFm4 HFp4 - 25 HFm1mc HF 1 -5.00 -3.00 0.00 0.00 0.01 30.00 3sub yes HFp1mc trackp1mc - 26 HFp1mc HF 1 3.00 5.00 0.00 0.00 0.01 30.00 3sub yes HFm1mc trackm1mc - 27 trackm1mc Tracker 1 -2.20 -1.40 0.00 0.00 0.30 3.00 3sub yes HFm1mc HFp1mc - 28 trackp1mc Tracker 1 1.40 2.20 0.00 0.00 0.30 3.00 3sub yes HFm1mc HFp1mc - 29 Castor1mc Castor 1 -6.55 -5.10 0.00 0.00 0.01 50.00 3sub yes HFp1mc trackp1mc + 0 HFm2 HF 2 -5.00 -3.00 0.00 0.00 0.01 30.00 3sub no HFp2 trackmid2 + 1 HFp2 HF 2 3.00 5.00 0.00 0.00 0.01 30.00 3sub no HFm2 trackmid2 + 2 HF2 HF 2 -5.00 -3.00 3.00 5.00 0.01 30.00 3sub no trackm2 trackp2 + 3 trackmid2 Tracker 2 -0.50 0.50 0.00 0.00 0.50 3.00 3sub no HFm2 HFp2 + 4 trackm2 Tracker 2 -2.00 -1.00 0.00 0.00 0.50 3.00 3sub no HFp2 HFm2 + 5 trackp2 Tracker 2 1.00 2.00 0.00 0.00 0.50 3.00 3sub no HFp2 HFm2 + 6 HFm3 HF 3 -5.00 -3.00 0.00 0.00 0.01 30.00 3sub no HFp3 trackmid3 + 7 HFp3 HF 3 3.00 5.00 0.00 0.00 0.01 30.00 3sub no HFm3 trackmid3 + 8 HF3 HF 3 -5.00 -3.00 3.00 5.00 0.01 30.00 3sub no trackm3 trackp3 + 9 trackmid3 Tracker 3 -0.50 0.50 0.00 0.00 0.50 3.00 3sub no HFm3 HFp3 + 10 trackm3 Tracker 3 -2.00 -1.00 0.00 0.00 0.50 3.00 3sub no HFp3 HFm3 + 11 trackp3 Tracker 3 1.00 2.00 0.00 0.00 0.50 3.00 3sub no HFp3 HFm3 */ #include - namespace hi { + // clang-format off enum EPNamesInd { - HFm1, - HFp1, - HF1, - trackm1, - trackp1, - Castor1, - HFm2, - HFp2, - HF2, - trackmid2, - trackm2, - trackp2, - Castor2, - HFm3, - HFp3, - HF3, - trackmid3, - trackm3, - trackp3, - HFm4, - HFp4, - HF4, - trackmid4, - trackm4, - trackp4, - HFm1mc, - HFp1mc, - trackm1mc, - trackp1mc, - EPBLANK + HFm2, HFp2, HF2, trackmid2, trackm2, + trackp2, HFm3, HFp3, HF3, trackmid3, + trackm3, trackp3, EPBLANK }; - const std::string EPNames[] = {"HFm1", "HFp1", "HF1", "trackm1", "trackp1", "Castor1", - "HFm2", "HFp2", "HF2", "trackmid2", "trackm2", "trackp2", - "Castor2", "HFm3", "HFp3", "HF3", "trackmid3", "trackm3", - "trackp3", "HFm4", "HFp4", "HF4", "trackmid4", "trackm4", - "trackp4", "HFm1mc", "HFp1mc", "trackm1mc", "trackp1mc"}; + const std::string EPNames[] = { + "HFm2", "HFp2", "HF2", "trackmid2", "trackm2", + "trackp2", "HFm3", "HFp3", "HF3", "trackmid3", + "trackm3", "trackp3" + }; - enum Detectors { Tracker, HF, Castor }; + enum Detectors { Tracker, HF, Castor, RPD }; - const int EPDet[] = {HF, HF, HF, Tracker, Tracker, Castor, HF, HF, HF, Tracker, - Tracker, Tracker, Castor, HF, HF, HF, Tracker, Tracker, Tracker, HF, - HF, HF, Tracker, Tracker, Tracker, HF, HF, Tracker, Tracker}; + const int EPDet[] = { + HF, HF, HF, Tracker, Tracker, + Tracker, HF, HF, HF, Tracker, + Tracker, Tracker + }; - const int EPOrder[] = {1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 1, 1, 1, 1}; + const int EPOrder[] = { + 2, 2, 2, 2, 2, + 2, 3, 3, 3, 3, + 3, 3 + }; - const double EPEtaMin1[] = {-5.00, 3.00, -5.00, -2.00, 1.00, -6.55, -5.00, 3.00, -5.00, -0.75, - -2.00, 1.00, -6.55, -5.00, 3.00, -5.00, -0.75, -2.00, 1.00, -5.00, - 3.00, -5.00, -0.75, -2.00, 1.00, -5.00, 3.00, -2.20, 1.40}; + const double EPEtaMin1[] = { + -5.00, 3.00, -5.00, -0.50, -2.00, + 1.00, -5.00, 3.00, -5.00, -0.50, + -2.00, 1.00 + }; - const double EPEtaMax1[] = {-3.00, 5.00, -3.00, -1.00, 2.00, -5.10, -3.00, 5.00, -3.00, 0.75, - -1.00, 2.00, -5.10, -3.00, 5.00, -3.00, 0.75, -1.00, 2.00, -3.00, - 5.00, -3.00, 0.75, -1.00, 2.00, -3.00, 5.00, -1.40, 2.20}; + const double EPEtaMax1[] = { + -3.00, 5.00, -3.00, 0.50, -1.00, + 2.00, -3.00, 5.00, -3.00, 0.50, + -1.00, 2.00 + }; - const double EPEtaMin2[] = {0.00, 0.00, 3.00, 0.00, 0.00, 0.00, 0.00, 0.00, 3.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, - 3.00, 0.00, 0.00, 0.00, 0.00, 0.00, 3.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}; + const double EPEtaMin2[] = { + 0.00, 0.00, 3.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 3.00, 0.00, + 0.00, 0.00 + }; - const double EPEtaMax2[] = {0.00, 0.00, 5.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, - 5.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00}; + const double EPEtaMax2[] = { + 0.00, 0.00, 5.00, 0.00, 0.00, + 0.00, 0.00, 0.00, 5.00, 0.00, + 0.00, 0.00 + }; - const double minTransverse[] = {0.01, 0.01, 0.01, 0.30, 0.30, 0.01, 0.01, 0.01, 0.01, 0.30, - 0.30, 0.30, 0.01, 0.01, 0.01, 0.01, 0.30, 0.30, 0.30, 0.01, - 0.01, 0.01, 0.30, 0.30, 0.30, 0.01, 0.01, 0.30, 0.30}; + const double minTransverse[] = { + 0.01, 0.01, 0.01, 0.50, 0.50, + 0.50, 0.01, 0.01, 0.01, 0.50, + 0.50, 0.50 + }; - const double maxTransverse[] = {30.00, 30.00, 30.00, 3.00, 3.00, 50.00, 30.00, 30.00, 30.00, 3.00, - 3.00, 3.00, 50.00, 30.00, 30.00, 30.00, 3.00, 3.00, 3.00, 30.00, - 30.00, 30.00, 3.00, 3.00, 3.00, 30.00, 30.00, 3.00, 3.00}; + const double maxTransverse[] = { + 30.00, 30.00, 30.00, 3.00, 3.00, + 3.00, 30.00, 30.00, 30.00, 3.00, + 3.00, 3.00 + }; - const std::string ResCalcType[] = {"3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", - "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", - "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub"}; + const std::string ResCalcType[] = { + "3sub", "3sub", "3sub", "3sub", "3sub", + "3sub", "3sub", "3sub", "3sub", "3sub", + "3sub", "3sub" + }; - const std::string MomConsWeight[] = {"no", "no", "no", "no", "no", "no", "no", "no", "no", "no", - "no", "no", "no", "no", "no", "no", "no", "no", "no", "no", - "no", "no", "no", "no", "no", "yes", "yes", "yes", "yes"}; + const std::string MomConsWeight[] = { + "no", "no", "no", "no", "no", + "no", "no", "no", "no", "no", + "no", "no" + }; - const int RCMate1[] = {HFp1, HFm1, trackm1, HFm1, HFm1, HFp1, HFp2, HFm2, trackm2, HFm2, - HFm2, HFm2, trackmid2, HFp3, HFm3, trackm3, HFm3, HFm3, HFm3, HFp4, - HFm4, trackm4, HFm4, HFm4, HFm4, HFp1mc, HFm1mc, HFm1mc, HFm1mc}; + const int RCMate1[] = { + HFp2, HFm2, trackm2, HFm2, HFp2, + HFp2, HFp3, HFm3, trackm3, HFm3, + HFp3, HFp3 + }; - const int RCMate2[] = {trackp1, trackm1, trackp1, HFp1, HFp1, trackp1, trackmid2, trackmid2, - trackp2, HFp2, HFp2, HFp2, HFp2, trackmid3, trackmid3, trackp3, - HFp3, HFp3, HFp3, trackmid4, trackmid4, trackp4, HFp4, HFp4, - HFp4, trackp1mc, trackm1mc, HFp1mc, HFp1mc}; + const int RCMate2[] = { + trackmid2, trackmid2, trackp2, HFp2, HFm2, + HFm2, trackmid3, trackmid3, trackp3, HFp3, + HFm3, HFm3 + }; - static const int NumEPNames = 29; + static const int NumEPNames = 12; + // clang-format on } // namespace hi #endif diff --git a/RecoHI/HiEvtPlaneAlgos/interface/LoadEPDB.h b/RecoHI/HiEvtPlaneAlgos/interface/LoadEPDB.h index eae35f4434468..c32598bb230e8 100644 --- a/RecoHI/HiEvtPlaneAlgos/interface/LoadEPDB.h +++ b/RecoHI/HiEvtPlaneAlgos/interface/LoadEPDB.h @@ -25,8 +25,6 @@ #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h" #include -//using namespace hi; - class LoadEPDB { public: explicit LoadEPDB(const edm::ESHandle flatparmsDB_, HiEvtPlaneFlatten** flat) { @@ -55,7 +53,6 @@ class LoadEPDB { } else if (i >= Hbins && i < Hbins + Obins) { flat[indx]->setXoffDB(i - Hbins, thisBin->x[j]); flat[indx]->setYoffDB(i - Hbins, thisBin->y[j]); - } else if (i >= Hbins + Obins && i < Hbins + 2 * Obins) { flat[indx]->setPtDB(i - Hbins - Obins, thisBin->x[j]); flat[indx]->setPt2DB(i - Hbins - Obins, thisBin->y[j]); @@ -73,8 +70,9 @@ class LoadEPDB { for (int j = 0; j < ncentbins; j++) { const RPFlatParams::EP* thisBin = &(flatparmsDB_->m_table[Hbins + 2 * Obins + cbins + j + 1]); if (fabs(centbinning - 1.) < 0.01) { - for (int i = 0; i < hi::NumEPNames; i++) + for (int i = 0; i < hi::NumEPNames; i++) { flat[i]->setCentRes1(j, thisBin->x[i], thisBin->y[i]); + } } if (fabs(centbinning - 2.) < 0.01) { for (int i = 0; i < hi::NumEPNames; i++) diff --git a/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py b/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py index 5cbdd824ea845..13efd5a09736d 100644 --- a/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py +++ b/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py @@ -1,25 +1,53 @@ import FWCore.ParameterSet.Config as cms hiEvtPlane = cms.EDProducer("EvtPlaneProducer", + centralityVariable = cms.string("HFtowers"), + centralityBinTag = cms.InputTag("centralityBin","HFtowers"), vertexTag = cms.InputTag("hiSelectedVertex"), caloTag = cms.InputTag("towerMaker"), castorTag = cms.InputTag("CastorTowerReco"), trackTag = cms.InputTag("hiGeneralTracks"), - centralityBinTag = cms.InputTag("centralityBin","HFtowers"), - centralityVariable = cms.string("HFtowers"), + lostTag = cms.InputTag("lostTracks"), + chi2MapTag = cms.InputTag("packedPFCandidateTrackChi2"), + chi2MapLostTag = cms.InputTag("lostTrackChi2"), nonDefaultGlauberModel = cms.string(""), - FlatOrder = cms.int32(9), - NumFlatBins = cms.int32(40), - CentBinCompression = cms.int32(5), - caloCentRef = cms.double(80.), - caloCentRefWidth = cms.double(5.0), loadDB = cms.bool(False), minet = cms.double(-1.), - maxet = cms.double(-1.), + maxet = cms.double(-1), minpt = cms.double(0.3), maxpt = cms.double(3.0), minvtx = cms.double(-25.), maxvtx = cms.double(25.), - dzerr = cms.double(10.), - chi2 = cms.double(40.) + flatnvtxbins = cms.int32(10), + flatminvtx = cms.double(-15.0), + flatdelvtx = cms.double(3.0), + dzdzerror = cms.double(3.0), + d0d0error = cms.double(3.0), + pterror = cms.double(0.1), + chi2perlayer = cms.double(0.18), + dzdzerror_pix = cms.double(10.), + chi2 = cms.double(40.), + nhitsValid = cms.int32(11), + FlatOrder = cms.int32(9), + NumFlatBins = cms.int32(40), + caloCentRef = cms.double(80.), + caloCentRefWidth = cms.double(5.0), + CentBinCompression = cms.int32(5), + cutEra = cms.int32(2) # 0:ppReco, 1:HIReco, 2:Pixel, 3: GenMC ) + +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(hiEvtPlane, + vertexTag = cms.InputTag("offlinePrimaryVertices"), + trackTag = cms.InputTag("generalTracks"), + minet = cms.double(0.01), + minpt = cms.double(0.5), + minvtx = cms.double(-15.), + maxvtx = cms.double(15.), + dzdzerror_pix = cms.double(40.), + caloCentRef = cms.double(-1), + caloCentRefWidth = cms.double(-1), + cutEra = cms.int32(0) +) diff --git a/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py b/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py index 33e7b661b7c8f..c32f5b744ec8a 100644 --- a/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py +++ b/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py @@ -1,22 +1,17 @@ import FWCore.ParameterSet.Config as cms -# AOD content -RecoHiEvtPlaneAOD = cms.PSet( - outputCommands = cms.untracked.vstring( - 'keep recoEvtPlanes_hiEvtPlane_*_*', - 'keep ZDCRecHitsSorted_zdcreco_*_*', - 'keep ZDCDataFramesSorted_hcalDigis_*_*', - 'keep HFRecHitsSorted_hfreco_*_*') -) +RecoHiEvtPlaneFEVT = cms.PSet( + outputCommands = cms.untracked.vstring('keep recoEvtPlanes_hiEvtPlane_*_*') + ) -# RECO content RecoHiEvtPlaneRECO = cms.PSet( - outputCommands = cms.untracked.vstring() -) -RecoHiEvtPlaneRECO.outputCommands.extend(RecoHiEvtPlaneAOD.outputCommands) + outputCommands = cms.untracked.vstring('keep recoEvtPlanes_hiEvtPlane_*_*') + ) -# FEVT content -RecoHiEvtPlaneFEVT = cms.PSet( - outputCommands = cms.untracked.vstring() -) -RecoHiEvtPlaneFEVT.outputCommands.extend(RecoHiEvtPlaneRECO.outputCommands) +RecoHiEvtPlaneAOD = cms.PSet( + outputCommands = cms.untracked.vstring('keep recoEvtPlanes_hiEvtPlane_*_*', + 'keep ZDCRecHitsSorted_zdcreco_*_*', + 'keep ZDCDataFramesSorted_hcalDigis_*_*', + 'keep HFRecHitsSorted_hfreco_*_*' + ) + ) diff --git a/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py b/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py index 33ac26edf1bdb..29cc0f551572e 100644 --- a/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py +++ b/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py @@ -1,19 +1,22 @@ import FWCore.ParameterSet.Config as cms hiEvtPlaneFlat = cms.EDProducer('HiEvtPlaneFlatProducer', - vertexTag = cms.InputTag("hiSelectedVertex"), - centralityTag = cms.InputTag("hiCentrality"), - centralityBinTag = cms.InputTag("centralityBin","HFtowers"), centralityVariable = cms.string("HFtowers"), - nonDefaultGlauberModel = cms.string(""), + centralityBinTag = cms.InputTag("centralityBin","HFtowers"), + centralityTag = cms.InputTag("hiCentrality"), + vertexTag = cms.InputTag("offlinePrimaryVertices"), inputPlanesTag = cms.InputTag("hiEvtPlane"), - trackTag = cms.InputTag("hiGeneralTracks"), + nonDefaultGlauberModel = cms.string(""), + trackTag = cms.InputTag("generalTracks"), FlatOrder = cms.int32(9), NumFlatBins = cms.int32(40), - Noffmin = cms.int32 (-1), - Noffmax = cms.int32 (10000), + minvtx = cms.double(-15.), + maxvtx = cms.double(15.), + flatnvtxbins = cms.int32(10), + flatminvtx = cms.double(-15.0), + flatdelvtx = cms.double(3.0), + caloCentRef = cms.double(-1.), + caloCentRefWidth = cms.double(-1.), CentBinCompression = cms.int32(5), - caloCentRef = cms.double(80.), - caloCentRefWidth = cms.double(5.0), useOffsetPsi = cms.bool(True) ) diff --git a/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc index fe27144182afe..30834b29e53bd 100644 --- a/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc +++ b/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc @@ -4,7 +4,7 @@ // Class: EvtPlaneProducer // /**\class EvtPlaneProducer EvtPlaneProducer.cc RecoHI/EvtPlaneProducer/src/EvtPlaneProducer.cc - + Description: Implementation: @@ -21,6 +21,7 @@ Description: #include #include #include +#include // user include files #include "FWCore/Framework/interface/Frameworkfwd.h" @@ -50,7 +51,6 @@ Description: #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/VertexReco/interface/Vertex.h" #include "DataFormats/VertexReco/interface/VertexFwd.h" -#include #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h" #include "CondFormats/HIObjects/interface/RPFlatParams.h" #include "CondFormats/DataRecord/interface/HeavyIonRPRcd.h" @@ -62,8 +62,11 @@ Description: #include "DataFormats/Common/interface/Ref.h" #include "DataFormats/Common/interface/RefVector.h" +#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h" +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h" #include "RecoHI/HiEvtPlaneAlgos/interface/LoadEPDB.h" +#include "RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h" using namespace std; using namespace hi; @@ -171,6 +174,7 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { private: GenPlane *rp[NumEPNames]; + EPCuts *cuts; void produce(edm::Event &, const edm::EventSetup &) override; @@ -190,6 +194,7 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { edm::InputTag caloTag_; edm::EDGetTokenT caloToken; edm::Handle caloCollection_; + edm::EDGetTokenT caloTokenPF; edm::InputTag castorTag_; edm::EDGetTokenT> castorToken; @@ -197,10 +202,19 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { edm::InputTag trackTag_; edm::EDGetTokenT trackToken; + edm::InputTag losttrackTag_; + edm::EDGetTokenT losttrackToken; edm::Handle trackCollection_; + string strack; + string scalo; + edm::EDGetTokenT> packedToken; + edm::EDGetTokenT> lostToken; - edm::ESWatcher hiWatcher; - edm::ESWatcher hirpWatcher; + edm::InputTag chi2MapTag_; + edm::EDGetTokenT> chi2MapToken; + edm::InputTag chi2MapLostTag_; + edm::EDGetTokenT> chi2MapLostToken; + std::vector trkNormChi2; bool loadDB_; double minet_; @@ -209,15 +223,149 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { double maxpt_; double minvtx_; double maxvtx_; + int flatnvtxbins_; + double flatminvtx_; + double flatdelvtx_; + double dzdzerror_; + double d0d0error_; + double pterror_; + double chi2perlayer_; double dzerr_; + double dzdzerror_pix_; double chi2_; + int nhitsValid_; int FlatOrder_; int NumFlatBins_; double nCentBins_; double caloCentRef_; double caloCentRefWidth_; int CentBinCompression_; + int cutEra_; HiEvtPlaneFlatten *flat[NumEPNames]; + int evtCount; + TrackStructure track; + int pcnt = 0; + + edm::ESWatcher hiWatcher; + edm::ESWatcher hirpWatcher; + + void FillHF(TrackStructure track, double vz, int bin) { + double minet = minet_; + double maxet = maxet_; + for (int i = 0; i < NumEPNames; i++) { + if (EPDet[i] != HF) + continue; + if (minet_ < 0) + minet = minTransverse[i]; + if (maxet_ < 0) + maxet = maxTransverse[i]; + if (track.et < minet) + continue; + if (track.et > maxet) + continue; + if (EPEtaMin2[i] == EPEtaMax2[i]) { + if (track.eta < EPEtaMin1[i]) + continue; + if (track.eta > EPEtaMax1[i]) + continue; + } else { + if (track.eta < EPEtaMin1[i]) + continue; + if (track.eta > EPEtaMax2[i]) + continue; + if (track.eta > EPEtaMax1[i] && track.eta < EPEtaMin2[i]) + continue; + } + double w = track.et; + if (loadDB_) + w = track.et * flat[i]->getEtScale(vz, bin); + if (EPOrder[i] == 1) { + if (MomConsWeight[i][0] == 'y' && loadDB_) { + w = flat[i]->getW(track.et, vz, bin); + } + //if(track.eta<0 ) w=-w; + } + rp[i]->addParticle(w, track.et, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta); + } + }; + + void FillCastor(TrackStructure track, double vz, int bin) { + double minet = minet_; + double maxet = maxet_; + for (int i = 0; i < NumEPNames; i++) { + if (EPDet[i] == Castor) { + if (minet_ < 0) + minet = minTransverse[i]; + if (maxet_ < 0) + maxet = maxTransverse[i]; + if (track.et < minet) + continue; + if (track.et > maxet) + continue; + if (EPEtaMin2[i] == EPEtaMax2[i]) { + if (track.eta < EPEtaMin1[i]) + continue; + if (track.eta > EPEtaMax1[i]) + continue; + } else { + if (track.eta < EPEtaMin1[i]) + continue; + if (track.eta > EPEtaMax2[i]) + continue; + if (track.eta > EPEtaMax1[i] && track.eta < EPEtaMin2[i]) + continue; + } + double w = track.et; + if (EPOrder[i] == 1) { + if (MomConsWeight[i][0] == 'y' && loadDB_) { + w = flat[i]->getW(track.et, vz, bin); + } + //if(track.eta<0 ) w=-w; + } + rp[i]->addParticle(w, track.et, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta); + } + } + } + + void FillTracker(TrackStructure track, double vz, int bin) { + double minpt = minpt_; + double maxpt = maxpt_; + for (int i = 0; i < NumEPNames; i++) { + if (EPDet[i] == Tracker) { + if (minpt_ < 0) + minpt = minTransverse[i]; + if (maxpt_ < 0) + maxpt = maxTransverse[i]; + if (track.pt < minpt) + continue; + if (track.pt > maxpt) + continue; + if (EPEtaMin2[i] == EPEtaMax2[i]) { + if (track.eta < EPEtaMin1[i]) + continue; + if (track.eta > EPEtaMax1[i]) + continue; + } else { + if (track.eta < EPEtaMin1[i]) + continue; + if (track.eta > EPEtaMax2[i]) + continue; + if (track.eta > EPEtaMax1[i] && track.eta < EPEtaMin2[i]) + continue; + } + double w = track.pt; + if (w > 2.5) + w = 2.0; //v2 starts decreasing above ~2.5 GeV/c + if (EPOrder[i] == 1) { + if (MomConsWeight[i][0] == 'y' && loadDB_) { + w = flat[i]->getW(track.pt, vz, bin); + } + //if(track.eta<0) w=-w; + } + rp[i]->addParticle(w, track.pt, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta); + } + } + }; }; EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig) @@ -227,6 +375,9 @@ EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig) caloTag_(iConfig.getParameter("caloTag")), castorTag_(iConfig.getParameter("castorTag")), trackTag_(iConfig.getParameter("trackTag")), + losttrackTag_(iConfig.getParameter("lostTag")), + chi2MapTag_(iConfig.getParameter("chi2MapTag")), + chi2MapLostTag_(iConfig.getParameter("chi2MapLostTag")), loadDB_(iConfig.getParameter("loadDB")), minet_(iConfig.getParameter("minet")), maxet_(iConfig.getParameter("maxet")), @@ -234,13 +385,44 @@ EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig) maxpt_(iConfig.getParameter("maxpt")), minvtx_(iConfig.getParameter("minvtx")), maxvtx_(iConfig.getParameter("maxvtx")), - dzerr_(iConfig.getParameter("dzerr")), + flatnvtxbins_(iConfig.getParameter("flatnvtxbins")), + flatminvtx_(iConfig.getParameter("flatminvtx")), + flatdelvtx_(iConfig.getParameter("flatdelvtx")), + dzdzerror_(iConfig.getParameter("dzdzerror")), + d0d0error_(iConfig.getParameter("d0d0error")), + pterror_(iConfig.getParameter("pterror")), + chi2perlayer_(iConfig.getParameter("chi2perlayer")), + dzdzerror_pix_(iConfig.getParameter("dzdzerror_pix")), chi2_(iConfig.getParameter("chi2")), + nhitsValid_(iConfig.getParameter("nhitsValid")), FlatOrder_(iConfig.getParameter("FlatOrder")), NumFlatBins_(iConfig.getParameter("NumFlatBins")), caloCentRef_(iConfig.getParameter("caloCentRef")), caloCentRefWidth_(iConfig.getParameter("caloCentRefWidth")), - CentBinCompression_(iConfig.getParameter("CentBinCompression")) { + CentBinCompression_(iConfig.getParameter("CentBinCompression")), + cutEra_(iConfig.getParameter("cutEra")) + +{ + switch (cutEra_) { + case 0: + cuts = new EPCuts( + EP_ERA::ppReco, pterror_, dzdzerror_, d0d0error_, chi2perlayer_, dzdzerror_pix_, chi2_, nhitsValid_); + break; + case 1: + cuts = new EPCuts( + EP_ERA::HIReco, pterror_, dzdzerror_, d0d0error_, chi2perlayer_, dzdzerror_pix_, chi2_, nhitsValid_); + break; + case 2: + cuts = new EPCuts( + EP_ERA::Pixel, pterror_, dzdzerror_, d0d0error_, chi2perlayer_, dzdzerror_pix_, chi2_, nhitsValid_); + break; + case 3: + cuts = new EPCuts( + EP_ERA::GenMC, pterror_, dzdzerror_, d0d0error_, chi2perlayer_, dzdzerror_pix_, chi2_, nhitsValid_); + break; + default: + cuts = nullptr; + } nCentBins_ = 200.; if (iConfig.exists("nonDefaultGlauberModel")) { @@ -252,11 +434,23 @@ EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig) vertexToken = consumes>(vertexTag_); - caloToken = consumes(caloTag_); - - castorToken = consumes>(castorTag_); - - trackToken = consumes(trackTag_); + strack = trackTag_.label(); + scalo = caloTag_.label(); + if (strack.find("packedPFCandidates") != std::string::npos) { + packedToken = consumes>(trackTag_); + lostToken = consumes>(losttrackTag_); + chi2MapToken = consumes>(chi2MapTag_); + chi2MapLostToken = consumes>(chi2MapLostTag_); + + } else { + if (scalo.find("particleFlow") != std::string::npos) { + caloTokenPF = consumes(caloTag_); + } else { + caloToken = consumes(caloTag_); + } + castorToken = consumes>(castorTag_); + trackToken = consumes(trackTag_); + } produces(); for (int i = 0; i < NumEPNames; i++) { @@ -264,7 +458,7 @@ EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig) } for (int i = 0; i < NumEPNames; i++) { flat[i] = new HiEvtPlaneFlatten(); - flat[i]->init(FlatOrder_, NumFlatBins_, EPNames[i], EPOrder[i]); + flat[i]->init(FlatOrder_, NumFlatBins_, flatnvtxbins_, flatminvtx_, flatdelvtx_, EPNames[i], EPOrder[i]); } } @@ -285,8 +479,7 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup using namespace edm; using namespace std; using namespace reco; - - if (loadDB_ && (hiWatcher.check(iSetup) || hirpWatcher.check(iSetup))) { + if (hiWatcher.check(iSetup)) { // //Get Size of Centrality Table // @@ -305,214 +498,207 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup } } } - - // - //Get flattening parameter file. - // - if (loadDB_) { - edm::ESHandle flatparmsDB_; - iSetup.get().get(flatparmsDB_); - LoadEPDB db(flatparmsDB_, flat); - if (!db.IsSuccess()) { - loadDB_ = kFALSE; - } + } + // + //Get flattening parameter file. + // + if (loadDB_ && hirpWatcher.check(iSetup)) { + edm::ESHandle flatparmsDB_; + iSetup.get().get(flatparmsDB_); + LoadEPDB db(flatparmsDB_, flat); + if (!db.IsSuccess()) { + loadDB_ = kFALSE; } - - } //rp record change - + } // //Get Centrality // int bin = 0; + int cbin = 0; if (loadDB_) { edm::Handle cbin_; iEvent.getByToken(centralityBinToken, cbin_); - int cbin = *cbin_; + cbin = *cbin_; bin = cbin / CentBinCompression_; } // //Get Vertex // - int vs_sell = 0.; - float vzr_sell; - iEvent.getByToken(vertexToken, vertex_); - const reco::VertexCollection *vertices3 = nullptr; - if (vertex_.isValid()) { - vertices3 = vertex_.product(); - vs_sell = vertices3->size(); - } - if (vs_sell > 0) { - vzr_sell = vertices3->begin()->z(); - } else - vzr_sell = -999.9; - // - for (int i = 0; i < NumEPNames; i++) - rp[i]->reset(); - if (vzr_sell < minvtx_ or vzr_sell > maxvtx_) + edm::Handle vertices; + iEvent.getByToken(vertexToken, vertices); + + //best vertex + const reco::Vertex &vtx = (*vertices)[0]; + double bestvz = vtx.z(); + double bestvx = vtx.x(); + double bestvy = vtx.y(); + double bestvzError = vtx.zError(); + double bestvxError = vtx.xError(); + double bestvyError = vtx.yError(); + math::XYZPoint bestvtx(bestvx, bestvy, bestvz); + if (bestvz < minvtx_ || bestvz > maxvtx_) return; - //calorimetry part - - double tower_eta, tower_phi; - double tower_energyet, tower_energyet_e, tower_energyet_h; - - iEvent.getByToken(caloToken, caloCollection_); + for (int i = 0; i < NumEPNames; i++) + rp[i]->reset(); + edm::Handle> chi2Map; + edm::Handle> cands; + edm::Handle calocands; + if (strack.find("packedPFCandidates") != std::string::npos) { + iEvent.getByToken(packedToken, cands); + iEvent.getByToken(chi2MapToken, chi2Map); + for (unsigned int i = 0, n = cands->size(); i < n; ++i) { + track = {}; + track.centbin = cbin; + const pat::PackedCandidate &pf = (*cands)[i]; + track.et = pf.et(); + track.eta = pf.eta(); + track.phi = pf.phi(); + track.pdgid = pf.pdgId(); + if (cuts->isGoodHF(track)) { + FillHF(track, bestvz, bin); + } + if (!pf.hasTrackDetails()) + continue; + const reco::Track &trk = pf.pseudoTrack(); + track.highPurity = pf.trackHighPurity(); + track.charge = trk.charge(); + if (!track.highPurity || track.charge == 0) + continue; + track.collection = 1; + track.eta = trk.eta(); + track.phi = trk.phi(); + track.pt = trk.pt(); + track.ptError = trk.ptError(); + track.numberOfValidHits = trk.numberOfValidHits(); + track.algos = trk.algo(); + track.dz = trk.dz(bestvtx); + track.dxy = -1. * trk.dxy(bestvtx); + track.dzError = sqrt(pow(trk.dzError(), 2) + pow(bestvzError, 2)); + track.dxyError = sqrt(pow(trk.dxyError(), 2) + bestvxError * bestvyError); + track.dzSig = track.dz / track.dzError; + track.dxySig = track.dxy / track.dxyError; + const reco::HitPattern &hit_pattern = trk.hitPattern(); + track.normalizedChi2 = (*chi2Map)[cands->ptrAt(i)]; + track.chi2layer = (*chi2Map)[cands->ptrAt(i)] / hit_pattern.trackerLayersWithMeasurement(); + if (cuts->isGoodTrack(track)) { + FillTracker(track, bestvz, bin); + } + } - if (caloCollection_.isValid()) { - for (CaloTowerCollection::const_iterator j = caloCollection_->begin(); j != caloCollection_->end(); j++) { - tower_eta = j->eta(); - tower_phi = j->phi(); - tower_energyet_e = j->emEt(); - tower_energyet_h = j->hadEt(); - tower_energyet = tower_energyet_e + tower_energyet_h; - double minet = minet_; - double maxet = maxet_; - for (int i = 0; i < NumEPNames; i++) { - if (minet_ < 0) - minet = minTransverse[i]; - if (maxet_ < 0) - maxet = maxTransverse[i]; - if (tower_energyet < minet) - continue; - if (tower_energyet > maxet) - continue; - if (EPDet[i] == HF) { - double w = tower_energyet; - if (loadDB_) - w = tower_energyet * flat[i]->getEtScale(vzr_sell, bin); - if (EPOrder[i] == 1) { - if (MomConsWeight[i][0] == 'y' && loadDB_) { - w = flat[i]->getW(tower_energyet, vzr_sell, bin); - } - if (tower_eta < 0) - w = -w; - } - rp[i]->addParticle(w, tower_energyet, sin(EPOrder[i] * tower_phi), cos(EPOrder[i] * tower_phi), tower_eta); - } + iEvent.getByToken(lostToken, cands); + iEvent.getByToken(chi2MapLostToken, chi2Map); + for (unsigned int i = 0, n = cands->size(); i < n; ++i) { + track = {}; + track.centbin = cbin; + if (cuts->isGoodHF(track)) { + FillHF(track, bestvz, bin); + } + const pat::PackedCandidate &pf = (*cands)[i]; + if (!pf.hasTrackDetails()) + continue; + const reco::Track &trk = pf.pseudoTrack(); + track.highPurity = pf.trackHighPurity(); + track.charge = trk.charge(); + if (!track.highPurity || track.charge == 0) + continue; + track.collection = 2; + track.eta = pf.eta(); + track.phi = pf.phi(); + track.et = pf.et(); + track.pdgid = pf.pdgId(); + track.pt = trk.pt(); + track.ptError = trk.ptError(); + track.numberOfValidHits = trk.numberOfValidHits(); + track.algos = trk.algo(); + track.dz = trk.dz(bestvtx); + track.dxy = -1. * trk.dxy(bestvtx); + track.dzError = sqrt(pow(trk.dzError(), 2) + pow(bestvzError, 2)); + track.dxyError = sqrt(pow(trk.dxyError(), 2) + bestvxError * bestvyError); + track.dzSig = track.dz / track.dzError; + track.dxySig = track.dxy / track.dxyError; + const reco::HitPattern &hit_pattern = trk.hitPattern(); + track.normalizedChi2 = (*chi2Map)[cands->ptrAt(i)]; + track.chi2layer = (*chi2Map)[cands->ptrAt(i)] / hit_pattern.trackerLayersWithMeasurement(); + if (cuts->isGoodTrack(track)) { + FillTracker(track, bestvz, bin); } } - } - //Castor part - - iEvent.getByToken(castorToken, castorCollection_); - - if (castorCollection_.isValid()) { - for (std::vector::const_iterator j = castorCollection_->begin(); j != castorCollection_->end(); - j++) { - tower_eta = j->eta(); - tower_phi = j->phi(); - tower_energyet = j->et(); - double minet = minet_; - double maxet = maxet_; - for (int i = 0; i < NumEPNames; i++) { - if (EPDet[i] == Castor) { - if (minet_ < 0) - minet = minTransverse[i]; - if (maxet_ < 0) - maxet = maxTransverse[i]; - if (tower_energyet < minet) - continue; - if (tower_energyet > maxet) - continue; - double w = tower_energyet; - if (EPOrder[i] == 1) { - if (MomConsWeight[i][0] == 'y' && loadDB_) { - w = flat[i]->getW(tower_energyet, vzr_sell, bin); - } - if (tower_eta < 0) - w = -w; + } else { + //calorimetry part + if (scalo.find("particleFlow") != std::string::npos) { + iEvent.getByToken(caloTokenPF, calocands); + if (cands.isValid()) { + for (unsigned int i = 0, n = calocands->size(); i < n; ++i) { + track = {}; + track.centbin = cbin; + const reco::PFCandidate &pf = (*calocands)[i]; + track.et = pf.et(); + track.eta = pf.eta(); + track.phi = pf.phi(); + track.pdgid = pf.pdgId(); + if (cuts->isGoodHF(track)) { + FillHF(track, bestvz, bin); } - rp[i]->addParticle(w, tower_energyet, sin(EPOrder[i] * tower_phi), cos(EPOrder[i] * tower_phi), tower_eta); + } + } + } else { + iEvent.getByToken(caloToken, caloCollection_); + if (caloCollection_.isValid()) { + for (CaloTowerCollection::const_iterator j = caloCollection_->begin(); j != caloCollection_->end(); j++) { + track.eta = j->eta(); + track.phi = j->phi(); + track.et = j->emEt() + j->hadEt(); + track.pdgid = 1; + if (cuts->isGoodHF(track)) + FillHF(track, bestvz, bin); } } } - } - - //Tracking part - double track_eta; - double track_phi; - double track_pt; - - double vzErr2 = 0.0, vxyErr = 0.0; - math::XYZPoint vtxPoint(0.0, 0.0, 0.0); - if (vertex_.isValid() && !vertex_->empty()) { - vtxPoint = vertex_->begin()->position(); - vzErr2 = (vertex_->begin()->zError()) * (vertex_->begin()->zError()); - vxyErr = vertex_->begin()->xError() * vertex_->begin()->yError(); - } - - iEvent.getByToken(trackToken, trackCollection_); - if (trackCollection_.isValid()) { - for (reco::TrackCollection::const_iterator j = trackCollection_->begin(); j != trackCollection_->end(); j++) { - bool accepted = true; - bool isPixel = false; - // determine if the track is a pixel track - if (j->numberOfValidHits() < 7) - isPixel = true; - - // determine the vertex significance - double d0 = 0.0, dz = 0.0, d0sigma = 0.0, dzsigma = 0.0; - d0 = -1. * j->dxy(vtxPoint); - dz = j->dz(vtxPoint); - d0sigma = sqrt(j->d0Error() * j->d0Error() + vxyErr); - dzsigma = sqrt(j->dzError() * j->dzError() + vzErr2); - - // cuts for pixel tracks - if (isPixel) { - // dz significance cut - if (fabs(dz / dzsigma) > dzerr_) - accepted = false; - // chi2/ndof cut - if (j->normalizedChi2() > chi2_) - accepted = false; - } - // cuts for full tracks - if (!isPixel) { - // dz and d0 significance cuts - if (fabs(dz / dzsigma) > 3) - accepted = false; - if (fabs(d0 / d0sigma) > 3) - accepted = false; - // pt resolution cut - if (j->ptError() / j->pt() > 0.1) - accepted = false; - // number of valid hits cut - if (j->numberOfValidHits() < 12) - accepted = false; + //Castor part + iEvent.getByToken(castorToken, castorCollection_); + if (castorCollection_.isValid()) { + for (std::vector::const_iterator j = castorCollection_->begin(); j != castorCollection_->end(); + j++) { + track.eta = j->eta(); + track.phi = j->phi(); + track.et = j->et(); + track.pdgid = 1; + if (cuts->isGoodCastor(track)) + FillCastor(track, bestvz, bin); } - if (accepted) { - track_eta = j->eta(); - track_phi = j->phi(); - track_pt = j->pt(); - double minpt = minpt_; - double maxpt = maxpt_; - for (int i = 0; i < NumEPNames; i++) { - if (minpt_ < 0) - minpt = minTransverse[i]; - if (maxpt_ < 0) - maxpt = maxTransverse[i]; - if (track_pt < minpt) - continue; - if (track_pt > maxpt) - continue; - if (EPDet[i] == Tracker) { - double w = track_pt; - if (w > 2.5) - w = 2.0; //v2 starts decreasing above ~2.5 GeV/c - if (EPOrder[i] == 1) { - if (MomConsWeight[i][0] == 'y' && loadDB_) { - w = flat[i]->getW(track_pt, vzr_sell, bin); - } - if (track_eta < 0) - w = -w; - } - rp[i]->addParticle(w, track_pt, sin(EPOrder[i] * track_phi), cos(EPOrder[i] * track_phi), track_eta); - } - } + } + //Tracking part + iEvent.getByToken(trackToken, trackCollection_); + if (trackCollection_.isValid()) { + for (reco::TrackCollection::const_iterator j = trackCollection_->begin(); j != trackCollection_->end(); j++) { + track.highPurity = j->quality(reco::TrackBase::highPurity); + track.charge = j->charge(); + if (!track.highPurity || track.charge == 0) + continue; + track.centbin = cbin; + track.collection = 0; + track.eta = j->eta(); + track.phi = j->phi(); + track.pt = j->pt(); + track.ptError = j->ptError(); + track.numberOfValidHits = j->numberOfValidHits(); + track.algos = j->algo(); + track.dz = j->dz(bestvtx); + track.dxy = -1. * j->dxy(bestvtx); + track.dzError = sqrt(pow(j->dzError(), 2) + pow(bestvzError, 2)); + track.dxyError = sqrt(pow(j->d0Error(), 2) + bestvxError * bestvyError); + track.dzSig = track.dz / track.dzError; + track.dxySig = track.dxy / track.dxyError; + track.normalizedChi2 = j->normalizedChi2(); + track.chi2layer = track.normalizedChi2 / j->hitPattern().trackerLayersWithMeasurement(); + + if (cuts->isGoodTrack(track)) + FillTracker(track, bestvz, bin); } - } //end for + } } auto evtplaneOutput = std::make_unique(); @@ -536,6 +722,7 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup } iEvent.put(std::move(evtplaneOutput)); + ++pcnt; } //define this as a plug-in diff --git a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc index b2a106e0a9ed0..ab3dca1f4bdd9 100644 --- a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc +++ b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc @@ -1,3 +1,4 @@ +#include #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/stream/EDProducer.h" @@ -19,6 +20,17 @@ #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/VertexReco/interface/Vertex.h" #include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "SimDataFormats/Track/interface/SimTrack.h" +#include "SimDataFormats/Track/interface/SimTrackContainer.h" +#include "SimDataFormats/Track/interface/CoreSimTrack.h" +#include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h" +#include "SimDataFormats/Vertex/interface/SimVertex.h" +#include "SimDataFormats/Vertex/interface/SimVertexContainer.h" +#include "SimDataFormats/TrackingHit/interface/PSimHit.h" +#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" +#include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h" +#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h" +#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h" #include "FWCore/ServiceRegistry/interface/Service.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" @@ -29,17 +41,17 @@ #include "CondFormats/HIObjects/interface/RPFlatParams.h" #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h" +#include +#include #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h" #include "RecoHI/HiEvtPlaneAlgos/interface/LoadEPDB.h" -#include "TList.h" -#include "TString.h" +using namespace std; +using namespace hi; #include -#include -#include -#include +using std::vector; // // class declaration @@ -80,15 +92,32 @@ class HiEvtPlaneFlatProducer : public edm::stream::EDProducer<> { const int FlatOrder_; int NumFlatBins_; + double minvtx_; + double maxvtx_; + int flatnvtxbins_; + double flatminvtx_; + double flatdelvtx_; double caloCentRef_; double caloCentRefWidth_; int CentBinCompression_; - int Noffmin_; - int Noffmax_; - HiEvtPlaneFlatten* flat[hi::NumEPNames]; + HiEvtPlaneFlatten* flat[NumEPNames]; bool useOffsetPsi_; + double nCentBins_; }; +// +// constants, enums and typedefs +// + +typedef std::vector TrackingParticleCollection; +typedef TrackingParticleRefVector::iterator tp_iterator; +// +// static data member definitions +// + +// +// constructors and destructor +// HiEvtPlaneFlatProducer::HiEvtPlaneFlatProducer(const edm::ParameterSet& iConfig) : centralityVariable_(iConfig.getParameter("centralityVariable")), centralityBinTag_(iConfig.getParameter("centralityBinTag")), @@ -98,12 +127,18 @@ HiEvtPlaneFlatProducer::HiEvtPlaneFlatProducer(const edm::ParameterSet& iConfig) trackTag_(iConfig.getParameter("trackTag")), FlatOrder_(iConfig.getParameter("FlatOrder")), NumFlatBins_(iConfig.getParameter("NumFlatBins")), + minvtx_(iConfig.getParameter("minvtx")), + maxvtx_(iConfig.getParameter("maxvtx")), + flatnvtxbins_(iConfig.getParameter("flatnvtxbins")), + flatminvtx_(iConfig.getParameter("flatminvtx")), + flatdelvtx_(iConfig.getParameter("flatdelvtx")), caloCentRef_(iConfig.getParameter("caloCentRef")), caloCentRefWidth_(iConfig.getParameter("caloCentRefWidth")), CentBinCompression_(iConfig.getParameter("CentBinCompression")), - Noffmin_(iConfig.getParameter("Noffmin")), - Noffmax_(iConfig.getParameter("Noffmax")), useOffsetPsi_(iConfig.getParameter("useOffsetPsi")) { + // UseEtHF = kFALSE; + nCentBins_ = 200.; + if (iConfig.exists("nonDefaultGlauberModel")) { centralityMC_ = iConfig.getParameter("nonDefaultGlauberModel"); } @@ -122,42 +157,41 @@ HiEvtPlaneFlatProducer::HiEvtPlaneFlatProducer(const edm::ParameterSet& iConfig) //register your products produces(); //now do what ever other initialization is needed - for (int i = 0; i < hi::NumEPNames; i++) { + for (int i = 0; i < NumEPNames; i++) { flat[i] = new HiEvtPlaneFlatten(); - flat[i]->init(FlatOrder_, NumFlatBins_, hi::EPNames[i], hi::EPOrder[i]); + flat[i]->init(FlatOrder_, NumFlatBins_, flatnvtxbins_, flatminvtx_, flatdelvtx_, EPNames[i], EPOrder[i]); } } HiEvtPlaneFlatProducer::~HiEvtPlaneFlatProducer() { // do anything here that needs to be done at desctruction time // (e.g. close files, deallocate resources etc.) - for (int i = 0; i < hi::NumEPNames; i++) { + for (int i = 0; i < NumEPNames; i++) { delete flat[i]; } } +// +// member functions +// + // ------------ method called to produce the data ------------ void HiEvtPlaneFlatProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { using namespace edm; using namespace std; using namespace reco; - using namespace hi; - // - //Get Flattening Parameters - // if (hiWatcher.check(iSetup)) { // //Get Size of Centrality Table // edm::ESHandle centDB_; iSetup.get().get(centralityLabel_, centDB_); - int nCentBins = centDB_->m_table.size(); - for (int i = 0; i < hi::NumEPNames; i++) { - flat[i]->setCaloCentRefBins(-1, -1); + nCentBins_ = centDB_->m_table.size(); + for (int i = 0; i < NumEPNames; i++) { if (caloCentRef_ > 0) { - int minbin = (caloCentRef_ - caloCentRefWidth_ / 2.) * nCentBins / 100.; - int maxbin = (caloCentRef_ + caloCentRefWidth_ / 2.) * nCentBins / 100.; + int minbin = (caloCentRef_ - caloCentRefWidth_ / 2.) * nCentBins_ / 100.; + int maxbin = (caloCentRef_ + caloCentRefWidth_ / 2.) * nCentBins_ / 100.; minbin /= CentBinCompression_; maxbin /= CentBinCompression_; if (minbin > 0 && maxbin >= minbin) { @@ -167,71 +201,84 @@ void HiEvtPlaneFlatProducer::produce(edm::Event& iEvent, const edm::EventSetup& } } } - + // + //Get flattening parameter file. + // if (hirpWatcher.check(iSetup)) { edm::ESHandle flatparmsDB_; iSetup.get().get(flatparmsDB_); LoadEPDB db(flatparmsDB_, flat); - if (!db.IsSuccess()) - return; - } + } //rp record change + // //Get Centrality // - - int bin = iEvent.get(centralityBinToken) / CentBinCompression_; - - if (Noffmin_ >= 0) { - int nOff = iEvent.get(centralityToken).Ntracks(); - if ((nOff < Noffmin_) or (nOff >= Noffmax_)) { - return; - } - } + int bin = 0; + int cbin = 0; + edm::Handle cbin_; + iEvent.getByToken(centralityBinToken, cbin_); + cbin = *cbin_; + bin = cbin / CentBinCompression_; + //double centval = cbin*100./(double) nCentBins_; // //Get Vertex // - int vs_sell; // vertex collection size - float vzr_sell; - auto const& vertices3 = iEvent.get(vertexToken); - vs_sell = vertices3.size(); - if (vs_sell > 0) { - vzr_sell = vertices3.begin()->z(); - } else - vzr_sell = -999.9; + edm::Handle vertices; + iEvent.getByToken(vertexToken, vertices); + + //best vertex + double bestvz = -999.9, bestvx = -999.9, bestvy = -999.9; + const reco::Vertex& vtx = (*vertices)[0]; + bestvz = vtx.z(); + bestvx = vtx.x(); + bestvy = vtx.y(); + math::XYZPoint bestvtx(bestvx, bestvy, bestvz); + if (bestvz < minvtx_ || bestvz > maxvtx_) + return; // //Get Event Planes // + edm::Handle evtPlanes_; + iEvent.getByToken(inputPlanesToken, evtPlanes_); + + if (!evtPlanes_.isValid()) { + return; + } + auto evtplaneOutput = std::make_unique(); - EvtPlane* ep[hi::NumEPNames]; - for (int i = 0; i < hi::NumEPNames; i++) { + EvtPlane* ep[NumEPNames]; + for (int i = 0; i < NumEPNames; i++) { ep[i] = nullptr; } int indx = 0; - for (auto const& rp : iEvent.get(inputPlanesToken)) { - double psiOffset = rp.angle(0); - double s = rp.sumSin(0); - double c = rp.sumCos(0); - uint m = rp.mult(); - + for (EvtPlaneCollection::const_iterator rp = evtPlanes_->begin(); rp != evtPlanes_->end(); rp++) { + double s = rp->sumSin(0); + double c = rp->sumCos(0); + uint m = rp->mult(); double soff = s; double coff = c; - if (useOffsetPsi_) { - soff = flat[indx]->getSoffset(s, vzr_sell, bin); - coff = flat[indx]->getCoffset(c, vzr_sell, bin); - psiOffset = flat[indx]->getOffsetPsi(soff, coff); + double psiOffset = -10; + double psiFlat = -10; + if (rp->angle(0) > -5) { + if (useOffsetPsi_) { + soff = flat[indx]->getSoffset(s, bestvz, bin); + coff = flat[indx]->getCoffset(c, bestvz, bin); + psiOffset = flat[indx]->getOffsetPsi(soff, coff); + } + psiFlat = flat[indx]->getFlatPsi(psiOffset, bestvz, bin); } - double psiFlat = flat[indx]->getFlatPsi(psiOffset, vzr_sell, bin); - ep[indx] = new EvtPlane(indx, 2, psiFlat, soff, coff, rp.sumw(), rp.sumw2(), rp.sumPtOrEt(), rp.sumPtOrEt2(), m); - ep[indx]->addLevel(0, rp.angle(0), s, c); - ep[indx]->addLevel(3, 0., rp.sumSin(3), rp.sumCos(3)); + ep[indx] = + new EvtPlane(indx, 2, psiFlat, soff, coff, rp->sumw(), rp->sumw2(), rp->sumPtOrEt(), rp->sumPtOrEt2(), m); + ep[indx]->addLevel(0, rp->angle(0), s, c); + ep[indx]->addLevel(3, 0., rp->sumSin(3), rp->sumCos(3)); if (useOffsetPsi_) ep[indx]->addLevel(1, psiOffset, soff, coff); ++indx; } - for (int i = 0; i < hi::NumEPNames; i++) { + for (int i = 0; i < NumEPNames; i++) { if (ep[i] != nullptr) evtplaneOutput->push_back(*ep[i]); } From 5075c815467bb37913a6890f60f02e7196fb5116 Mon Sep 17 00:00:00 2001 From: Quan Wang Date: Fri, 2 Oct 2020 06:25:57 +0200 Subject: [PATCH 2/5] [+] update for code review --- .../PatAlgos/python/slimming/slimming_cff.py | 5 +- RecoHI/HiEvtPlaneAlgos/BuildFile.xml | 3 - RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h | 88 ++++++------- .../interface/HiEvtPlaneFlatten.h | 53 ++++---- .../interface/HiEvtPlaneList.h | 55 ++++---- .../HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py | 20 +-- .../HiEvtPlaneAlgos/src/EvtPlaneProducer.cc | 117 ++++++++---------- .../src/HiEvtPlaneFlatProducer.cc | 91 +++++--------- 8 files changed, 190 insertions(+), 242 deletions(-) diff --git a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py index 8c1e497f0bf92..c64983c9ce861 100644 --- a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py +++ b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py @@ -67,11 +67,10 @@ pp_on_AA_2018.toReplaceWith(slimmingTask, slimmingTask.copyAndExclude([slimmedOOTPhotons])) from Configuration.Eras.Modifier_pp_on_PbPb_run3_cff import pp_on_PbPb_run3 from PhysicsTools.PatAlgos.slimming.hiPixelTracks_cfi import hiPixelTracks -from RecoHI.HiCentralityAlgos.CentralityBin_cfi import centralityBin from RecoHI.HiEvtPlaneAlgos.HiEvtPlane_cfi import hiEvtPlane from RecoHI.HiEvtPlaneAlgos.hiEvtPlaneFlat_cfi import hiEvtPlaneFlat (pp_on_AA_2018 | pp_on_PbPb_run3).toReplaceWith(slimmingTask, cms.Task(slimmingTask.copy(), hiPixelTracks, -centralityBin, hiEvtPlane, hiEvtPlaneFlat)) +hiEvtPlane, hiEvtPlaneFlat)) from Configuration.Eras.Modifier_pp_on_PbPb_run3_cff import pp_on_PbPb_run3 from PhysicsTools.PatAlgos.packedCandidateMuonID_cfi import packedCandidateMuonID @@ -80,7 +79,7 @@ from RecoHI.HiCentralityAlgos.hiHFfilters_cfi import hiHFfilters lostTrackChi2 = packedPFCandidateTrackChi2.clone(candidates = "lostTracks", doLostTracks = True) (pp_on_AA_2018 | pp_on_PbPb_run3).toReplaceWith( - slimmingTask, + slimmingTask, cms.Task(slimmingTask.copy(), packedCandidateMuonID, packedPFCandidateTrackChi2, lostTrackChi2, centralityBin, hiHFfilters)) from Configuration.Eras.Modifier_phase2_timing_cff import phase2_timing diff --git a/RecoHI/HiEvtPlaneAlgos/BuildFile.xml b/RecoHI/HiEvtPlaneAlgos/BuildFile.xml index e64ad6baf50de..9b54f07c134c9 100644 --- a/RecoHI/HiEvtPlaneAlgos/BuildFile.xml +++ b/RecoHI/HiEvtPlaneAlgos/BuildFile.xml @@ -6,9 +6,6 @@ - - - diff --git a/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h b/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h index 7445ef49f4483..f739101891bb1 100644 --- a/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h +++ b/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h @@ -1,5 +1,5 @@ -#ifndef __EPCuts__ -#define __EPCuts__ +#ifndef RecoHI_HiEvtPlaneAlgos_EPCuts_h +#define RecoHI_HiEvtPlaneAlgos_EPCuts_h namespace hi { @@ -7,41 +7,31 @@ namespace hi { struct TrackStructure { int centbin; - double eta; - double phi; - double et; - double pt; + float eta; + float phi; + float et; + float pt; int charge; int pdgid; int hits; int algos; int collection; - double dz; - double dxy; - double dzError; - double dxyError; - double ptError; + float dz; + float dxy; + float dzError; + float dxyError; + float ptError; bool highPurity; - double dzSig; - double dxySig; - double normalizedChi2; - double dzError_Pix; - double chi2layer; + float dzSig; + float dxySig; + float normalizedChi2; + float dzError_Pix; + float chi2layer; int numberOfValidHits; int pixel; }; class EPCuts { - private: - EP_ERA cutera_; - double pterror_; - double dzerror_; - double dxyerror_; - double chi2perlayer_; - double dzerror_Pix_; - double chi2Pix_; - int numberOfValidHits_; - public: explicit EPCuts(EP_ERA cutEra = EP_ERA::ppReco, double pterror = 0.1, @@ -60,9 +50,8 @@ namespace hi { chi2Pix_ = chi2Pix; numberOfValidHits_ = numberOfValidHits; } - ~EPCuts() {} - bool isGoodHF(TrackStructure track) { + bool isGoodHF(const TrackStructure& track) const { if (track.pdgid != 1 && track.pdgid != 2) return false; if (fabs(track.eta) < 3 || fabs(track.eta) > 5) @@ -70,24 +59,24 @@ namespace hi { return true; } - bool isGoodCastor(TrackStructure track) { return true; } + bool isGoodCastor(const TrackStructure& track) const { return true; } - bool isGoodTrack(TrackStructure track) { + bool isGoodTrack(const TrackStructure& track) const { if (cutera_ == EP_ERA::ppReco) - return TrackQuality_ppReco(track); + return trackQuality_ppReco(track); if (cutera_ == EP_ERA::HIReco) - return TrackQuality_HIReco(track); + return trackQuality_HIReco(track); if (cutera_ == EP_ERA::Pixel) - return TrackQuality_Pixel(track); + return trackQuality_Pixel(track); return false; } - bool TrackQuality_ppReco(TrackStructure track) { + bool trackQuality_ppReco(const TrackStructure& track) const { if (track.charge == 0) return false; if (!track.highPurity) return false; - if (track.ptError / track.pt > pterror_) + if (track.ptError > pterror_ * track.pt) return false; if (track.numberOfValidHits < numberOfValidHits_) return false; @@ -100,14 +89,14 @@ namespace hi { return true; } - bool TrackQuality_HIReco(TrackStructure track) { + bool trackQuality_HIReco(const TrackStructure& track) const { if (track.charge == 0) return false; if (!track.highPurity) return false; if (track.numberOfValidHits < numberOfValidHits_) return false; - if (track.ptError / track.pt > pterror_) + if (track.ptError > pterror_ * track.pt) return false; if (fabs(track.dxy / track.dxyError) > dxyerror_) return false; @@ -115,19 +104,21 @@ namespace hi { return false; if (track.chi2layer > chi2perlayer_) return false; - if (track.algos != 4 && track.algos != 5 && track.algos != 6 && track.algos != 7) + //if (track.algos != 4 && track.algos != 5 && track.algos != 6 && track.algos != 7) + if (track.algos != reco::TrackBase::initialStep && track.algos != reco::TrackBase::lowPtTripletStep && + track.algos != reco::TrackBase::pixelPairStep && track.algos != reco::TrackBase::detachedTripletStep) return false; return true; } - bool TrackQuality_Pixel(TrackStructure track) { + bool trackQuality_Pixel(const TrackStructure& track) const { if (track.charge == 0) return false; if (!track.highPurity) return false; bool bPix = false; int nHits = track.numberOfValidHits; - if (track.ptError / track.pt > pterror_) + if (track.ptError > pterror_ * track.pt) return false; if (track.pt < 2.4 and (nHits == 3 or nHits == 4 or nHits == 5 or nHits == 6)) bPix = true; @@ -136,10 +127,11 @@ namespace hi { return false; if (track.chi2layer > chi2perlayer_) return false; - if (track.ptError / track.pt > pterror_) + if (track.ptError > pterror_ * track.pt) return false; int algo = track.algos; - if (track.pt > 2.4 && algo != 4 && algo != 5 && algo != 6 && algo != 7) + if (track.pt > 2.4 && algo != reco::TrackBase::initialStep && algo != reco::TrackBase::lowPtTripletStep && + algo != reco::TrackBase::pixelPairStep && algo != reco::TrackBase::detachedTripletStep) return false; if (fabs(track.dxy / track.dxyError) > dxyerror_) return false; @@ -154,13 +146,23 @@ namespace hi { return true; } - bool TrackQuality_GenMC(TrackStructure track) { + bool TrackQuality_GenMC(const TrackStructure& track) const { if (track.charge == 0) return false; if (fabs(track.eta) > 2.4) return false; return true; } + + private: + EP_ERA cutera_; + double pterror_; + double dzerror_; + double dxyerror_; + double chi2perlayer_; + double dzerror_Pix_; + double chi2Pix_; + int numberOfValidHits_; }; } // namespace hi #endif diff --git a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h index 4dd857e297c6f..ab7d96c46e7b5 100644 --- a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h +++ b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h @@ -76,7 +76,7 @@ class HiEvtPlaneFlatten { } } - int getCutIndx(int centbin, double vtx, int iord) const { + int cutIndx(int centbin, double vtx, int iord) const { int cut; if (centbin < 0) return -1; @@ -90,7 +90,7 @@ class HiEvtPlaneFlatten { return cut; } - int getOffsetIndx(int centbin, double vtx) const { + int offsetIndx(int centbin, double vtx) const { int cut; if (centbin < 0) return -1; @@ -110,7 +110,7 @@ class HiEvtPlaneFlatten { for (int k = 0; k < hOrder_; k++) { double fsin = sin(vorder_ * (k + 1) * psi); double fcos = cos(vorder_ * (k + 1) * psi); - int indx = getCutIndx(centbin, vtx, k); + int indx = cutIndx(centbin, vtx, k); if (indx >= 0) { flatX_[indx] += fcos; flatY_[indx] += fsin; @@ -119,7 +119,7 @@ class HiEvtPlaneFlatten { } } void fillOffset(double s, double c, uint m, double vtx, int centbin) { - int indx = getOffsetIndx(centbin, vtx); + int indx = offsetIndx(centbin, vtx); if (indx >= 0) { xoff_[indx] += c; yoff_[indx] += s; @@ -128,7 +128,7 @@ class HiEvtPlaneFlatten { } } void fillPt(double ptval, double vtx, int centbin) { - int indx = getOffsetIndx(centbin, vtx); + int indx = offsetIndx(centbin, vtx); if (indx >= 0) { pt_[indx] += ptval; pt2_[indx] += ptval * ptval; @@ -141,9 +141,9 @@ class HiEvtPlaneFlatten { caloCentRefMaxBin_ = caloCentRefMaxBin; } - double getEtScale(double vtx, int centbin) const { - int refmin = getOffsetIndx(caloCentRefMinBin_, vtx); - int refmax = getOffsetIndx(caloCentRefMaxBin_, vtx); + double etScale(double vtx, int centbin) const { + int refmin = offsetIndx(caloCentRefMinBin_, vtx); + int refmax = offsetIndx(caloCentRefMaxBin_, vtx); double caloCentRefVal_ = 0; for (int i = refmin; i <= refmax; i++) { caloCentRefVal_ += getPtDB(i); @@ -151,16 +151,16 @@ class HiEvtPlaneFlatten { caloCentRefVal_ /= refmax - refmin + 1.; if (caloCentRefMinBin_ < 0) return 1.; - int indx = getOffsetIndx(centbin, vtx); + int indx = offsetIndx(centbin, vtx); if (indx < 0 || caloCentRefVal_ == 0 || getPtDB(indx) == 0) return 1.; return caloCentRefVal_ / getPtDB(indx); } double getW(double pt, double vtx, int centbin) const { - int indx = getOffsetIndx(centbin, vtx); + int indx = offsetIndx(centbin, vtx); if (indx >= 0) { - double scale = getEtScale(vtx, centbin); + double scale = etScale(vtx, centbin); double ptval = getPtDB(indx) * scale; double pt2val = getPt2DB(indx) * pow(scale, 2); double rv = pt * scale - pt2val / ptval; @@ -173,7 +173,7 @@ class HiEvtPlaneFlatten { double getFlatPsi(double psi, double vtx, int centbin) const { double correction = 0; for (int k = 0; k < hOrder_; k++) { - int indx = getCutIndx(centbin, vtx, k); + int indx = cutIndx(centbin, vtx, k); if (indx >= 0) correction += (2. / (double)((k + 1) * vorder_)) * (flatXDB_[indx] * sin(vorder_ * (k + 1) * psi) - flatYDB_[indx] * cos(vorder_ * (k + 1) * psi)); @@ -184,8 +184,8 @@ class HiEvtPlaneFlatten { return psi; } - double getSoffset(double s, double vtx, int centbin) const { - int indx = getOffsetIndx(centbin, vtx); + double soffset(double s, double vtx, int centbin) const { + int indx = offsetIndx(centbin, vtx); if (indx >= 0) { return s - yoffDB_[indx]; } else { @@ -193,15 +193,15 @@ class HiEvtPlaneFlatten { } } - double getCoffset(double c, double vtx, int centbin) const { - int indx = getOffsetIndx(centbin, vtx); + double coffset(double c, double vtx, int centbin) const { + int indx = offsetIndx(centbin, vtx); if (indx >= 0) return c - xoffDB_[indx]; else return c; } - double getOffsetPsi(double s, double c) const { + double offsetPsi(double s, double c) const { double psi = atan2(s, c) / vorder_; if ((fabs(s) < 1e-4) && (fabs(c) < 1e-4)) psi = 0.; @@ -210,34 +210,29 @@ class HiEvtPlaneFlatten { return psi; } - float getMinCent(int indx) { - int ibin = (int)(indx / nvtxbins_); + float minCent(int indx) const { + int ibin = indx / nvtxbins_; return ibin * 100. / nbins_; } - float getMaxCent(int indx) { - int ibin = (int)(indx / nvtxbins_); + float maxCent(int indx) const { + int ibin = indx / nvtxbins_; return (ibin + 1) * 100. / nbins_; } - double getMinVtx(int indx) { + double minVtx(int indx) const { int ivtx = indx - nvtxbins_ * (int)(indx / nvtxbins_); return minvtx_ + ivtx * delvtx_; } - double getMaxVtx(int indx) { + double maxVtx(int indx) const { int ivtx = indx - nvtxbins_ * (int)(indx / nvtxbins_); return minvtx_ + (ivtx + 1) * delvtx_; } std::string getRangeString(int indx) { char buf[120]; - sprintf(buf, - "%5.1f < cent < %5.1f; %4.1f < vtx < %4.1f", - getMinCent(indx), - getMaxCent(indx), - getMinVtx(indx), - getMaxVtx(indx)); + sprintf(buf, "%5.1f < cent < %5.1f; %4.1f < vtx < %4.1f", minCent(indx), maxCent(indx), minVtx(indx), maxVtx(indx)); return std::string(buf); } diff --git a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h index d3d10a6788088..ff634363248e9 100644 --- a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h +++ b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h @@ -25,87 +25,88 @@ namespace hi { trackm3, trackp3, EPBLANK }; - const std::string EPNames[] = { + static const int NumEPNames = 12; + + const std::array EPNames = {{ "HFm2", "HFp2", "HF2", "trackmid2", "trackm2", "trackp2", "HFm3", "HFp3", "HF3", "trackmid3", "trackm3", "trackp3" - }; + }}; enum Detectors { Tracker, HF, Castor, RPD }; - const int EPDet[] = { + const std::array EPDet = {{ HF, HF, HF, Tracker, Tracker, Tracker, HF, HF, HF, Tracker, Tracker, Tracker - }; + }}; - const int EPOrder[] = { + const std::array EPOrder = {{ 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3 - }; + }}; - const double EPEtaMin1[] = { + const std::array EPEtaMin1 = {{ -5.00, 3.00, -5.00, -0.50, -2.00, 1.00, -5.00, 3.00, -5.00, -0.50, -2.00, 1.00 - }; + }}; - const double EPEtaMax1[] = { + const std::array EPEtaMax1 = {{ -3.00, 5.00, -3.00, 0.50, -1.00, 2.00, -3.00, 5.00, -3.00, 0.50, -1.00, 2.00 - }; + }}; - const double EPEtaMin2[] = { + const std::array EPEtaMin2 = {{ 0.00, 0.00, 3.00, 0.00, 0.00, 0.00, 0.00, 0.00, 3.00, 0.00, 0.00, 0.00 - }; + }}; - const double EPEtaMax2[] = { + const std::array EPEtaMax2 = {{ 0.00, 0.00, 5.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00, 0.00, 0.00 - }; + }}; - const double minTransverse[] = { + const std::array minTransverse = {{ 0.01, 0.01, 0.01, 0.50, 0.50, 0.50, 0.01, 0.01, 0.01, 0.50, 0.50, 0.50 - }; + }}; - const double maxTransverse[] = { + const std::array maxTransverse = {{ 30.00, 30.00, 30.00, 3.00, 3.00, 3.00, 30.00, 30.00, 30.00, 3.00, 3.00, 3.00 - }; + }}; - const std::string ResCalcType[] = { + const std::array ResCalcType = {{ "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub", "3sub" - }; + }}; - const std::string MomConsWeight[] = { + const std::array MomConsWeight = {{ "no", "no", "no", "no", "no", "no", "no", "no", "no", "no", "no", "no" - }; + }}; - const int RCMate1[] = { + const std::array RCMate1 = {{ HFp2, HFm2, trackm2, HFm2, HFp2, HFp2, HFp3, HFm3, trackm3, HFm3, HFp3, HFp3 - }; + }}; - const int RCMate2[] = { + const std::array RCMate2 = {{ trackmid2, trackmid2, trackp2, HFp2, HFm2, HFm2, trackmid3, trackmid3, trackp3, HFp3, HFm3, HFm3 - }; + }}; - static const int NumEPNames = 12; // clang-format on } // namespace hi #endif diff --git a/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py b/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py index 13efd5a09736d..5f77961194ebf 100644 --- a/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py +++ b/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py @@ -40,14 +40,14 @@ from Configuration.Eras.Modifier_pp_on_PbPb_run3_cff import pp_on_PbPb_run3 (pp_on_AA_2018 | pp_on_PbPb_run3).toModify(hiEvtPlane, - vertexTag = cms.InputTag("offlinePrimaryVertices"), - trackTag = cms.InputTag("generalTracks"), - minet = cms.double(0.01), - minpt = cms.double(0.5), - minvtx = cms.double(-15.), - maxvtx = cms.double(15.), - dzdzerror_pix = cms.double(40.), - caloCentRef = cms.double(-1), - caloCentRefWidth = cms.double(-1), - cutEra = cms.int32(0) + vertexTag = "offlinePrimaryVertices", + trackTag = "generalTracks", + minet = 0.01, + minpt = 0.5, + minvtx = -15., + maxvtx = 15., + dzdzerror_pix = 40., + caloCentRef = -1, + caloCentRefWidth = -1, + cutEra = 0 ) diff --git a/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc index 30834b29e53bd..c1cc7aea99334 100644 --- a/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc +++ b/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc @@ -174,47 +174,45 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { private: GenPlane *rp[NumEPNames]; - EPCuts *cuts; void produce(edm::Event &, const edm::EventSetup &) override; // ----------member data --------------------------- + EPCuts *cuts; std::string centralityVariable_; std::string centralityLabel_; std::string centralityMC_; edm::InputTag centralityBinTag_; - edm::EDGetTokenT centralityBinToken; + edm::EDGetTokenT centralityBinToken_; edm::InputTag vertexTag_; - edm::EDGetTokenT> vertexToken; + edm::EDGetTokenT> vertexToken_; edm::Handle> vertex_; edm::InputTag caloTag_; - edm::EDGetTokenT caloToken; + edm::EDGetTokenT caloToken_; edm::Handle caloCollection_; - edm::EDGetTokenT caloTokenPF; + edm::EDGetTokenT caloTokenPF_; edm::InputTag castorTag_; - edm::EDGetTokenT> castorToken; + edm::EDGetTokenT> castorToken_; edm::Handle> castorCollection_; edm::InputTag trackTag_; - edm::EDGetTokenT trackToken; + edm::EDGetTokenT trackToken_; edm::InputTag losttrackTag_; - edm::EDGetTokenT losttrackToken; edm::Handle trackCollection_; - string strack; - string scalo; - edm::EDGetTokenT> packedToken; - edm::EDGetTokenT> lostToken; + bool bStrack_packedPFCandidates_; + bool bScalo_particleFlow_; + edm::EDGetTokenT> packedToken_; + edm::EDGetTokenT> lostToken_; edm::InputTag chi2MapTag_; - edm::EDGetTokenT> chi2MapToken; + edm::EDGetTokenT> chi2MapToken_; edm::InputTag chi2MapLostTag_; - edm::EDGetTokenT> chi2MapLostToken; - std::vector trkNormChi2; + edm::EDGetTokenT> chi2MapLostToken_; bool loadDB_; double minet_; @@ -242,14 +240,12 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { int CentBinCompression_; int cutEra_; HiEvtPlaneFlatten *flat[NumEPNames]; - int evtCount; TrackStructure track; - int pcnt = 0; edm::ESWatcher hiWatcher; edm::ESWatcher hirpWatcher; - void FillHF(TrackStructure track, double vz, int bin) { + void fillHF(TrackStructure track, double vz, int bin) { double minet = minet_; double maxet = maxet_; for (int i = 0; i < NumEPNames; i++) { @@ -278,12 +274,11 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { } double w = track.et; if (loadDB_) - w = track.et * flat[i]->getEtScale(vz, bin); + w = track.et * flat[i]->etScale(vz, bin); if (EPOrder[i] == 1) { if (MomConsWeight[i][0] == 'y' && loadDB_) { w = flat[i]->getW(track.et, vz, bin); } - //if(track.eta<0 ) w=-w; } rp[i]->addParticle(w, track.et, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta); } @@ -320,7 +315,6 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { if (MomConsWeight[i][0] == 'y' && loadDB_) { w = flat[i]->getW(track.et, vz, bin); } - //if(track.eta<0 ) w=-w; } rp[i]->addParticle(w, track.et, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta); } @@ -360,7 +354,6 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { if (MomConsWeight[i][0] == 'y' && loadDB_) { w = flat[i]->getW(track.pt, vz, bin); } - //if(track.eta<0) w=-w; } rp[i]->addParticle(w, track.pt, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta); } @@ -430,26 +423,26 @@ EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig) } centralityLabel_ = centralityVariable_ + centralityMC_; - centralityBinToken = consumes(centralityBinTag_); + centralityBinToken_ = consumes(centralityBinTag_); - vertexToken = consumes>(vertexTag_); + vertexToken_ = consumes>(vertexTag_); - strack = trackTag_.label(); - scalo = caloTag_.label(); - if (strack.find("packedPFCandidates") != std::string::npos) { - packedToken = consumes>(trackTag_); - lostToken = consumes>(losttrackTag_); - chi2MapToken = consumes>(chi2MapTag_); - chi2MapLostToken = consumes>(chi2MapLostTag_); + bStrack_packedPFCandidates_ = (trackTag_.label().find("packedPFCandidates") != std::string::npos); + bScalo_particleFlow_ = (caloTag_.label().find("particleFlow") != std::string::npos); + if (bStrack_packedPFCandidates_) { + packedToken_ = consumes>(trackTag_); + lostToken_ = consumes>(losttrackTag_); + chi2MapToken_ = consumes>(chi2MapTag_); + chi2MapLostToken_ = consumes>(chi2MapLostTag_); } else { - if (scalo.find("particleFlow") != std::string::npos) { - caloTokenPF = consumes(caloTag_); + if (bScalo_particleFlow_) { + caloTokenPF_ = consumes(caloTag_); } else { - caloToken = consumes(caloTag_); + caloToken_ = consumes(caloTag_); } - castorToken = consumes>(castorTag_); - trackToken = consumes(trackTag_); + castorToken_ = consumes>(castorTag_); + trackToken_ = consumes(trackTag_); } produces(); @@ -516,19 +509,14 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup int bin = 0; int cbin = 0; if (loadDB_) { - edm::Handle cbin_; - iEvent.getByToken(centralityBinToken, cbin_); - cbin = *cbin_; + cbin = iEvent.get(centralityBinToken_); bin = cbin / CentBinCompression_; } // //Get Vertex // - edm::Handle vertices; - iEvent.getByToken(vertexToken, vertices); - //best vertex - const reco::Vertex &vtx = (*vertices)[0]; + const reco::Vertex &vtx = iEvent.get(vertexToken_)[0]; double bestvz = vtx.z(); double bestvx = vtx.x(); double bestvy = vtx.y(); @@ -536,17 +524,19 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup double bestvxError = vtx.xError(); double bestvyError = vtx.yError(); math::XYZPoint bestvtx(bestvx, bestvy, bestvz); - if (bestvz < minvtx_ || bestvz > maxvtx_) - return; + math::Error<3>::type vtx_cov = vtx.covariance(); + // Produce the EP regardless of vz position + // if (bestvz < minvtx_ || bestvz > maxvtx_) + // return; for (int i = 0; i < NumEPNames; i++) rp[i]->reset(); edm::Handle> chi2Map; edm::Handle> cands; edm::Handle calocands; - if (strack.find("packedPFCandidates") != std::string::npos) { - iEvent.getByToken(packedToken, cands); - iEvent.getByToken(chi2MapToken, chi2Map); + if (bStrack_packedPFCandidates_) { + iEvent.getByToken(packedToken_, cands); + iEvent.getByToken(chi2MapToken_, chi2Map); for (unsigned int i = 0, n = cands->size(); i < n; ++i) { track = {}; track.centbin = cbin; @@ -556,7 +546,7 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup track.phi = pf.phi(); track.pdgid = pf.pdgId(); if (cuts->isGoodHF(track)) { - FillHF(track, bestvz, bin); + fillHF(track, bestvz, bin); } if (!pf.hasTrackDetails()) continue; @@ -572,10 +562,10 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup track.ptError = trk.ptError(); track.numberOfValidHits = trk.numberOfValidHits(); track.algos = trk.algo(); - track.dz = trk.dz(bestvtx); - track.dxy = -1. * trk.dxy(bestvtx); - track.dzError = sqrt(pow(trk.dzError(), 2) + pow(bestvzError, 2)); - track.dxyError = sqrt(pow(trk.dxyError(), 2) + bestvxError * bestvyError); + track.dz = std::abs(trk.dz(bestvtx)); + track.dxy = std::abs(trk.dxy(bestvtx)); + track.dzError = std::hypot(trk.dzError(), bestvzError); + track.dxyError = trk.dxyError(bestvtx, vtx_cov); track.dzSig = track.dz / track.dzError; track.dxySig = track.dxy / track.dxyError; const reco::HitPattern &hit_pattern = trk.hitPattern(); @@ -586,13 +576,13 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup } } - iEvent.getByToken(lostToken, cands); - iEvent.getByToken(chi2MapLostToken, chi2Map); + iEvent.getByToken(lostToken_, cands); + iEvent.getByToken(chi2MapLostToken_, chi2Map); for (unsigned int i = 0, n = cands->size(); i < n; ++i) { track = {}; track.centbin = cbin; if (cuts->isGoodHF(track)) { - FillHF(track, bestvz, bin); + fillHF(track, bestvz, bin); } const pat::PackedCandidate &pf = (*cands)[i]; if (!pf.hasTrackDetails()) @@ -627,8 +617,8 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup } else { //calorimetry part - if (scalo.find("particleFlow") != std::string::npos) { - iEvent.getByToken(caloTokenPF, calocands); + if (bScalo_particleFlow_) { + iEvent.getByToken(caloTokenPF_, calocands); if (cands.isValid()) { for (unsigned int i = 0, n = calocands->size(); i < n; ++i) { track = {}; @@ -639,12 +629,12 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup track.phi = pf.phi(); track.pdgid = pf.pdgId(); if (cuts->isGoodHF(track)) { - FillHF(track, bestvz, bin); + fillHF(track, bestvz, bin); } } } } else { - iEvent.getByToken(caloToken, caloCollection_); + iEvent.getByToken(caloToken_, caloCollection_); if (caloCollection_.isValid()) { for (CaloTowerCollection::const_iterator j = caloCollection_->begin(); j != caloCollection_->end(); j++) { track.eta = j->eta(); @@ -652,13 +642,13 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup track.et = j->emEt() + j->hadEt(); track.pdgid = 1; if (cuts->isGoodHF(track)) - FillHF(track, bestvz, bin); + fillHF(track, bestvz, bin); } } } //Castor part - iEvent.getByToken(castorToken, castorCollection_); + iEvent.getByToken(castorToken_, castorCollection_); if (castorCollection_.isValid()) { for (std::vector::const_iterator j = castorCollection_->begin(); j != castorCollection_->end(); j++) { @@ -671,7 +661,7 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup } } //Tracking part - iEvent.getByToken(trackToken, trackCollection_); + iEvent.getByToken(trackToken_, trackCollection_); if (trackCollection_.isValid()) { for (reco::TrackCollection::const_iterator j = trackCollection_->begin(); j != trackCollection_->end(); j++) { track.highPurity = j->quality(reco::TrackBase::highPurity); @@ -722,7 +712,6 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup } iEvent.put(std::move(evtplaneOutput)); - ++pcnt; } //define this as a plug-in diff --git a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc index ab3dca1f4bdd9..711cca39a6107 100644 --- a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc +++ b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc @@ -20,17 +20,6 @@ #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/VertexReco/interface/Vertex.h" #include "DataFormats/VertexReco/interface/VertexFwd.h" -#include "SimDataFormats/Track/interface/SimTrack.h" -#include "SimDataFormats/Track/interface/SimTrackContainer.h" -#include "SimDataFormats/Track/interface/CoreSimTrack.h" -#include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h" -#include "SimDataFormats/Vertex/interface/SimVertex.h" -#include "SimDataFormats/Vertex/interface/SimVertexContainer.h" -#include "SimDataFormats/TrackingHit/interface/PSimHit.h" -#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h" -#include "SimDataFormats/TrackingHit/interface/UpdatablePSimHit.h" -#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h" -#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h" #include "FWCore/ServiceRegistry/interface/Service.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" @@ -72,19 +61,19 @@ class HiEvtPlaneFlatProducer : public edm::stream::EDProducer<> { std::string centralityMC_; edm::InputTag centralityBinTag_; - edm::EDGetTokenT centralityBinToken; + edm::EDGetTokenT centralityBinToken_; edm::InputTag centralityTag_; - edm::EDGetTokenT centralityToken; + edm::EDGetTokenT centralityToken_; edm::InputTag vertexTag_; - edm::EDGetTokenT> vertexToken; + edm::EDGetTokenT> vertexToken_; edm::InputTag inputPlanesTag_; - edm::EDGetTokenT inputPlanesToken; + edm::EDGetTokenT inputPlanesToken_; edm::InputTag trackTag_; - edm::EDGetTokenT trackToken; + edm::EDGetTokenT trackToken_; edm::Handle trackCollection_; edm::ESWatcher hiWatcher; @@ -104,17 +93,6 @@ class HiEvtPlaneFlatProducer : public edm::stream::EDProducer<> { bool useOffsetPsi_; double nCentBins_; }; -// -// constants, enums and typedefs -// - -typedef std::vector TrackingParticleCollection; -typedef TrackingParticleRefVector::iterator tp_iterator; - -// -// static data member definitions -// - // // constructors and destructor // @@ -136,7 +114,6 @@ HiEvtPlaneFlatProducer::HiEvtPlaneFlatProducer(const edm::ParameterSet& iConfig) caloCentRefWidth_(iConfig.getParameter("caloCentRefWidth")), CentBinCompression_(iConfig.getParameter("CentBinCompression")), useOffsetPsi_(iConfig.getParameter("useOffsetPsi")) { - // UseEtHF = kFALSE; nCentBins_ = 200.; if (iConfig.exists("nonDefaultGlauberModel")) { @@ -144,15 +121,15 @@ HiEvtPlaneFlatProducer::HiEvtPlaneFlatProducer(const edm::ParameterSet& iConfig) } centralityLabel_ = centralityVariable_ + centralityMC_; - centralityBinToken = consumes(centralityBinTag_); + centralityBinToken_ = consumes(centralityBinTag_); - centralityToken = consumes(centralityTag_); + centralityToken_ = consumes(centralityTag_); - vertexToken = consumes>(vertexTag_); + vertexToken_ = consumes>(vertexTag_); - trackToken = consumes(trackTag_); + trackToken_ = consumes(trackTag_); - inputPlanesToken = consumes(inputPlanesTag_); + inputPlanesToken_ = consumes(inputPlanesTag_); //register your products produces(); @@ -215,37 +192,26 @@ void HiEvtPlaneFlatProducer::produce(edm::Event& iEvent, const edm::EventSetup& // int bin = 0; int cbin = 0; - edm::Handle cbin_; - iEvent.getByToken(centralityBinToken, cbin_); - cbin = *cbin_; + cbin = iEvent.get(centralityBinToken_); bin = cbin / CentBinCompression_; - //double centval = cbin*100./(double) nCentBins_; // //Get Vertex // - edm::Handle vertices; - iEvent.getByToken(vertexToken, vertices); //best vertex - double bestvz = -999.9, bestvx = -999.9, bestvy = -999.9; - const reco::Vertex& vtx = (*vertices)[0]; + double bestvz = -999.9; + const reco::Vertex& vtx = iEvent.get(vertexToken_)[0]; bestvz = vtx.z(); - bestvx = vtx.x(); - bestvy = vtx.y(); - math::XYZPoint bestvtx(bestvx, bestvy, bestvz); - if (bestvz < minvtx_ || bestvz > maxvtx_) - return; + // Produce the EP regardless of vz position + // if (bestvz < minvtx_ || bestvz > maxvtx_) + // return; // //Get Event Planes // edm::Handle evtPlanes_; - iEvent.getByToken(inputPlanesToken, evtPlanes_); - - if (!evtPlanes_.isValid()) { - return; - } + iEvent.getByToken(inputPlanesToken_, evtPlanes_); auto evtplaneOutput = std::make_unique(); EvtPlane* ep[NumEPNames]; @@ -253,26 +219,25 @@ void HiEvtPlaneFlatProducer::produce(edm::Event& iEvent, const edm::EventSetup& ep[i] = nullptr; } int indx = 0; - for (EvtPlaneCollection::const_iterator rp = evtPlanes_->begin(); rp != evtPlanes_->end(); rp++) { - double s = rp->sumSin(0); - double c = rp->sumCos(0); - uint m = rp->mult(); + for (auto&& rp : (*evtPlanes_)) { + double s = rp.sumSin(0); + double c = rp.sumCos(0); + uint m = rp.mult(); double soff = s; double coff = c; double psiOffset = -10; double psiFlat = -10; - if (rp->angle(0) > -5) { + if (rp.angle(0) > -5) { if (useOffsetPsi_) { - soff = flat[indx]->getSoffset(s, bestvz, bin); - coff = flat[indx]->getCoffset(c, bestvz, bin); - psiOffset = flat[indx]->getOffsetPsi(soff, coff); + soff = flat[indx]->soffset(s, bestvz, bin); + coff = flat[indx]->coffset(c, bestvz, bin); + psiOffset = flat[indx]->offsetPsi(soff, coff); } psiFlat = flat[indx]->getFlatPsi(psiOffset, bestvz, bin); } - ep[indx] = - new EvtPlane(indx, 2, psiFlat, soff, coff, rp->sumw(), rp->sumw2(), rp->sumPtOrEt(), rp->sumPtOrEt2(), m); - ep[indx]->addLevel(0, rp->angle(0), s, c); - ep[indx]->addLevel(3, 0., rp->sumSin(3), rp->sumCos(3)); + ep[indx] = new EvtPlane(indx, 2, psiFlat, soff, coff, rp.sumw(), rp.sumw2(), rp.sumPtOrEt(), rp.sumPtOrEt2(), m); + ep[indx]->addLevel(0, rp.angle(0), s, c); + ep[indx]->addLevel(3, 0., rp.sumSin(3), rp.sumCos(3)); if (useOffsetPsi_) ep[indx]->addLevel(1, psiOffset, soff, coff); ++indx; From 85242d576df06bdcc53d81e45b37f2a34889bc7d Mon Sep 17 00:00:00 2001 From: Quan Wang Date: Sun, 18 Oct 2020 22:02:30 +0200 Subject: [PATCH 3/5] [+] update for code review --- RecoHI/HiEvtPlaneAlgos/BuildFile.xml | 2 - RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h | 20 +- .../interface/HiEvtPlaneFlatten.h | 2 +- .../python/RecoHiEvtPlane_EventContent_cff.py | 29 +- .../HiEvtPlaneAlgos/src/EvtPlaneProducer.cc | 356 +++++++----------- .../src/HiEvtPlaneFlatProducer.cc | 5 +- 6 files changed, 170 insertions(+), 244 deletions(-) diff --git a/RecoHI/HiEvtPlaneAlgos/BuildFile.xml b/RecoHI/HiEvtPlaneAlgos/BuildFile.xml index 9b54f07c134c9..847d75c7055d8 100644 --- a/RecoHI/HiEvtPlaneAlgos/BuildFile.xml +++ b/RecoHI/HiEvtPlaneAlgos/BuildFile.xml @@ -1,6 +1,5 @@ - @@ -11,7 +10,6 @@ - diff --git a/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h b/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h index f739101891bb1..65ba8b014a472 100644 --- a/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h +++ b/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h @@ -54,7 +54,7 @@ namespace hi { bool isGoodHF(const TrackStructure& track) const { if (track.pdgid != 1 && track.pdgid != 2) return false; - if (fabs(track.eta) < 3 || fabs(track.eta) > 5) + if (std::abs(track.eta) < 3 || std::abs(track.eta) > 5) return false; return true; } @@ -82,9 +82,9 @@ namespace hi { return false; if (track.chi2layer > chi2perlayer_) return false; - if (fabs(track.dxy / track.dxyError) > dxyerror_) + if (std::abs(track.dxy) > dxyerror_ * track.dxyError) return false; - if (fabs(track.dz / track.dzError) > dzerror_) + if (std::abs(track.dz) > dzerror_ * track.dzError) return false; return true; } @@ -98,9 +98,9 @@ namespace hi { return false; if (track.ptError > pterror_ * track.pt) return false; - if (fabs(track.dxy / track.dxyError) > dxyerror_) + if (std::abs(track.dxy) > dxyerror_ * track.dxyError) return false; - if (fabs(track.dz / track.dzError) > dzerror_) + if (std::abs(track.dz) > dzerror_ * track.dzError) return false; if (track.chi2layer > chi2perlayer_) return false; @@ -133,23 +133,23 @@ namespace hi { if (track.pt > 2.4 && algo != reco::TrackBase::initialStep && algo != reco::TrackBase::lowPtTripletStep && algo != reco::TrackBase::pixelPairStep && algo != reco::TrackBase::detachedTripletStep) return false; - if (fabs(track.dxy / track.dxyError) > dxyerror_) + if (std::abs(track.dxy) > dxyerror_ * track.dxyError) return false; - if (fabs(track.dz / track.dzError) > dzerror_) + if (std::abs(track.dz) > dzerror_ * track.dzError) return false; } else { if (track.chi2layer > chi2Pix_) return false; - if (fabs(track.dz / track.dzError) > dzerror_Pix_) + if (std::abs(track.dz) > dzerror_Pix_ * track.dzError) return false; } return true; } - bool TrackQuality_GenMC(const TrackStructure& track) const { + bool trackQuality_GenMC(const TrackStructure& track) const { if (track.charge == 0) return false; - if (fabs(track.eta) > 2.4) + if (std::abs(track.eta) > 2.4) return false; return true; } diff --git a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h index ab7d96c46e7b5..01ca2ca215f24 100644 --- a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h +++ b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneFlatten.h @@ -230,7 +230,7 @@ class HiEvtPlaneFlatten { return minvtx_ + (ivtx + 1) * delvtx_; } - std::string getRangeString(int indx) { + std::string rangeString(int indx) { char buf[120]; sprintf(buf, "%5.1f < cent < %5.1f; %4.1f < vtx < %4.1f", minCent(indx), maxCent(indx), minVtx(indx), maxVtx(indx)); return std::string(buf); diff --git a/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py b/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py index c32f5b744ec8a..33e7b661b7c8f 100644 --- a/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py +++ b/RecoHI/HiEvtPlaneAlgos/python/RecoHiEvtPlane_EventContent_cff.py @@ -1,17 +1,22 @@ import FWCore.ParameterSet.Config as cms -RecoHiEvtPlaneFEVT = cms.PSet( - outputCommands = cms.untracked.vstring('keep recoEvtPlanes_hiEvtPlane_*_*') - ) +# AOD content +RecoHiEvtPlaneAOD = cms.PSet( + outputCommands = cms.untracked.vstring( + 'keep recoEvtPlanes_hiEvtPlane_*_*', + 'keep ZDCRecHitsSorted_zdcreco_*_*', + 'keep ZDCDataFramesSorted_hcalDigis_*_*', + 'keep HFRecHitsSorted_hfreco_*_*') +) +# RECO content RecoHiEvtPlaneRECO = cms.PSet( - outputCommands = cms.untracked.vstring('keep recoEvtPlanes_hiEvtPlane_*_*') - ) + outputCommands = cms.untracked.vstring() +) +RecoHiEvtPlaneRECO.outputCommands.extend(RecoHiEvtPlaneAOD.outputCommands) -RecoHiEvtPlaneAOD = cms.PSet( - outputCommands = cms.untracked.vstring('keep recoEvtPlanes_hiEvtPlane_*_*', - 'keep ZDCRecHitsSorted_zdcreco_*_*', - 'keep ZDCDataFramesSorted_hcalDigis_*_*', - 'keep HFRecHitsSorted_hfreco_*_*' - ) - ) +# FEVT content +RecoHiEvtPlaneFEVT = cms.PSet( + outputCommands = cms.untracked.vstring() +) +RecoHiEvtPlaneFEVT.outputCommands.extend(RecoHiEvtPlaneRECO.outputCommands) diff --git a/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc index c1cc7aea99334..dfbfcb5997397 100644 --- a/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc +++ b/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc @@ -178,7 +178,7 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { void produce(edm::Event &, const edm::EventSetup &) override; // ----------member data --------------------------- - EPCuts *cuts; + EPCuts cuts_; std::string centralityVariable_; std::string centralityLabel_; @@ -206,8 +206,8 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { edm::Handle trackCollection_; bool bStrack_packedPFCandidates_; bool bScalo_particleFlow_; - edm::EDGetTokenT> packedToken_; - edm::EDGetTokenT> lostToken_; + edm::EDGetTokenT packedToken_; + edm::EDGetTokenT lostToken_; edm::InputTag chi2MapTag_; edm::EDGetTokenT> chi2MapToken_; @@ -240,38 +240,27 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { int CentBinCompression_; int cutEra_; HiEvtPlaneFlatten *flat[NumEPNames]; - TrackStructure track; + TrackStructure track_; - edm::ESWatcher hiWatcher; - edm::ESWatcher hirpWatcher; + edm::ESWatcher hiWatcher_; + edm::ESWatcher hirpWatcher_; - void fillHF(TrackStructure track, double vz, int bin) { + void fillHF(const TrackStructure &track, double vz, int bin) { double minet = minet_; double maxet = maxet_; for (int i = 0; i < NumEPNames; i++) { if (EPDet[i] != HF) continue; + if (track.et < minet) + continue; + if (track.et > maxet) + continue; if (minet_ < 0) minet = minTransverse[i]; if (maxet_ < 0) maxet = maxTransverse[i]; - if (track.et < minet) - continue; - if (track.et > maxet) + if (not passEta(track.eta, i)) continue; - if (EPEtaMin2[i] == EPEtaMax2[i]) { - if (track.eta < EPEtaMin1[i]) - continue; - if (track.eta > EPEtaMax1[i]) - continue; - } else { - if (track.eta < EPEtaMin1[i]) - continue; - if (track.eta > EPEtaMax2[i]) - continue; - if (track.eta > EPEtaMax1[i] && track.eta < EPEtaMin2[i]) - continue; - } double w = track.et; if (loadDB_) w = track.et * flat[i]->etScale(vz, bin); @@ -284,32 +273,21 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { } }; - void FillCastor(TrackStructure track, double vz, int bin) { + void fillCastor(const TrackStructure &track, double vz, int bin) { double minet = minet_; double maxet = maxet_; for (int i = 0; i < NumEPNames; i++) { if (EPDet[i] == Castor) { + if (track.et < minet) + continue; + if (track.et > maxet) + continue; if (minet_ < 0) minet = minTransverse[i]; if (maxet_ < 0) maxet = maxTransverse[i]; - if (track.et < minet) - continue; - if (track.et > maxet) + if (not passEta(track.eta, i)) continue; - if (EPEtaMin2[i] == EPEtaMax2[i]) { - if (track.eta < EPEtaMin1[i]) - continue; - if (track.eta > EPEtaMax1[i]) - continue; - } else { - if (track.eta < EPEtaMin1[i]) - continue; - if (track.eta > EPEtaMax2[i]) - continue; - if (track.eta > EPEtaMax1[i] && track.eta < EPEtaMin2[i]) - continue; - } double w = track.et; if (EPOrder[i] == 1) { if (MomConsWeight[i][0] == 'y' && loadDB_) { @@ -321,32 +299,38 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { } } - void FillTracker(TrackStructure track, double vz, int bin) { + bool passEta(float eta, int i) { + if (EPEtaMin2[i] == EPEtaMax2[i]) { + if (eta < EPEtaMin1[i]) + return false; + if (eta > EPEtaMax1[i]) + return false; + } else { + if (eta < EPEtaMin1[i]) + return false; + if (eta > EPEtaMax2[i]) + return false; + if (eta > EPEtaMax1[i] && eta < EPEtaMin2[i]) + return false; + } + return true; + } + + void fillTracker(const TrackStructure &track, double vz, int bin) { double minpt = minpt_; double maxpt = maxpt_; for (int i = 0; i < NumEPNames; i++) { if (EPDet[i] == Tracker) { + if (track.pt < minpt) + continue; + if (track.pt > maxpt) + continue; if (minpt_ < 0) minpt = minTransverse[i]; if (maxpt_ < 0) maxpt = maxTransverse[i]; - if (track.pt < minpt) + if (not passEta(track.eta, i)) continue; - if (track.pt > maxpt) - continue; - if (EPEtaMin2[i] == EPEtaMax2[i]) { - if (track.eta < EPEtaMin1[i]) - continue; - if (track.eta > EPEtaMax1[i]) - continue; - } else { - if (track.eta < EPEtaMin1[i]) - continue; - if (track.eta > EPEtaMax2[i]) - continue; - if (track.eta > EPEtaMax1[i] && track.eta < EPEtaMin2[i]) - continue; - } double w = track.pt; if (w > 2.5) w = 2.0; //v2 starts decreasing above ~2.5 GeV/c @@ -396,26 +380,10 @@ EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig) cutEra_(iConfig.getParameter("cutEra")) { - switch (cutEra_) { - case 0: - cuts = new EPCuts( - EP_ERA::ppReco, pterror_, dzdzerror_, d0d0error_, chi2perlayer_, dzdzerror_pix_, chi2_, nhitsValid_); - break; - case 1: - cuts = new EPCuts( - EP_ERA::HIReco, pterror_, dzdzerror_, d0d0error_, chi2perlayer_, dzdzerror_pix_, chi2_, nhitsValid_); - break; - case 2: - cuts = new EPCuts( - EP_ERA::Pixel, pterror_, dzdzerror_, d0d0error_, chi2perlayer_, dzdzerror_pix_, chi2_, nhitsValid_); - break; - case 3: - cuts = new EPCuts( - EP_ERA::GenMC, pterror_, dzdzerror_, d0d0error_, chi2perlayer_, dzdzerror_pix_, chi2_, nhitsValid_); - break; - default: - cuts = nullptr; - } + if (cutEra_ > 3) + throw edm::Exception(edm::errors::Configuration) << "wrong range in cutEra parameter"; + cuts_ = EPCuts( + static_cast(cutEra_), pterror_, dzdzerror_, d0d0error_, chi2perlayer_, dzdzerror_pix_, chi2_, nhitsValid_); nCentBins_ = 200.; if (iConfig.exists("nonDefaultGlauberModel")) { @@ -430,8 +398,8 @@ EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig) bStrack_packedPFCandidates_ = (trackTag_.label().find("packedPFCandidates") != std::string::npos); bScalo_particleFlow_ = (caloTag_.label().find("particleFlow") != std::string::npos); if (bStrack_packedPFCandidates_) { - packedToken_ = consumes>(trackTag_); - lostToken_ = consumes>(losttrackTag_); + packedToken_ = consumes(trackTag_); + lostToken_ = consumes(losttrackTag_); chi2MapToken_ = consumes>(chi2MapTag_); chi2MapLostToken_ = consumes>(chi2MapLostTag_); @@ -472,7 +440,7 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup using namespace edm; using namespace std; using namespace reco; - if (hiWatcher.check(iSetup)) { + if (hiWatcher_.check(iSetup)) { // //Get Size of Centrality Table // @@ -495,7 +463,7 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup // //Get flattening parameter file. // - if (loadDB_ && hirpWatcher.check(iSetup)) { + if (loadDB_ && hirpWatcher_.check(iSetup)) { edm::ESHandle flatparmsDB_; iSetup.get().get(flatparmsDB_); LoadEPDB db(flatparmsDB_, flat); @@ -521,8 +489,6 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup double bestvx = vtx.x(); double bestvy = vtx.y(); double bestvzError = vtx.zError(); - double bestvxError = vtx.xError(); - double bestvyError = vtx.yError(); math::XYZPoint bestvtx(bestvx, bestvy, bestvz); math::Error<3>::type vtx_cov = vtx.covariance(); // Produce the EP regardless of vz position @@ -535,159 +501,117 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup edm::Handle> cands; edm::Handle calocands; if (bStrack_packedPFCandidates_) { - iEvent.getByToken(packedToken_, cands); - iEvent.getByToken(chi2MapToken_, chi2Map); - for (unsigned int i = 0, n = cands->size(); i < n; ++i) { - track = {}; - track.centbin = cbin; - const pat::PackedCandidate &pf = (*cands)[i]; - track.et = pf.et(); - track.eta = pf.eta(); - track.phi = pf.phi(); - track.pdgid = pf.pdgId(); - if (cuts->isGoodHF(track)) { - fillHF(track, bestvz, bin); - } - if (!pf.hasTrackDetails()) - continue; - const reco::Track &trk = pf.pseudoTrack(); - track.highPurity = pf.trackHighPurity(); - track.charge = trk.charge(); - if (!track.highPurity || track.charge == 0) - continue; - track.collection = 1; - track.eta = trk.eta(); - track.phi = trk.phi(); - track.pt = trk.pt(); - track.ptError = trk.ptError(); - track.numberOfValidHits = trk.numberOfValidHits(); - track.algos = trk.algo(); - track.dz = std::abs(trk.dz(bestvtx)); - track.dxy = std::abs(trk.dxy(bestvtx)); - track.dzError = std::hypot(trk.dzError(), bestvzError); - track.dxyError = trk.dxyError(bestvtx, vtx_cov); - track.dzSig = track.dz / track.dzError; - track.dxySig = track.dxy / track.dxyError; - const reco::HitPattern &hit_pattern = trk.hitPattern(); - track.normalizedChi2 = (*chi2Map)[cands->ptrAt(i)]; - track.chi2layer = (*chi2Map)[cands->ptrAt(i)] / hit_pattern.trackerLayersWithMeasurement(); - if (cuts->isGoodTrack(track)) { - FillTracker(track, bestvz, bin); + for (int idx = 1; idx < 3; idx++) { + if (idx == 1) { + iEvent.getByToken(packedToken_, cands); + iEvent.getByToken(chi2MapToken_, chi2Map); } - } - - iEvent.getByToken(lostToken_, cands); - iEvent.getByToken(chi2MapLostToken_, chi2Map); - for (unsigned int i = 0, n = cands->size(); i < n; ++i) { - track = {}; - track.centbin = cbin; - if (cuts->isGoodHF(track)) { - fillHF(track, bestvz, bin); + if (idx == 2) { + iEvent.getByToken(lostToken_, cands); + iEvent.getByToken(chi2MapLostToken_, chi2Map); } - const pat::PackedCandidate &pf = (*cands)[i]; - if (!pf.hasTrackDetails()) - continue; - const reco::Track &trk = pf.pseudoTrack(); - track.highPurity = pf.trackHighPurity(); - track.charge = trk.charge(); - if (!track.highPurity || track.charge == 0) - continue; - track.collection = 2; - track.eta = pf.eta(); - track.phi = pf.phi(); - track.et = pf.et(); - track.pdgid = pf.pdgId(); - track.pt = trk.pt(); - track.ptError = trk.ptError(); - track.numberOfValidHits = trk.numberOfValidHits(); - track.algos = trk.algo(); - track.dz = trk.dz(bestvtx); - track.dxy = -1. * trk.dxy(bestvtx); - track.dzError = sqrt(pow(trk.dzError(), 2) + pow(bestvzError, 2)); - track.dxyError = sqrt(pow(trk.dxyError(), 2) + bestvxError * bestvyError); - track.dzSig = track.dz / track.dzError; - track.dxySig = track.dxy / track.dxyError; - const reco::HitPattern &hit_pattern = trk.hitPattern(); - track.normalizedChi2 = (*chi2Map)[cands->ptrAt(i)]; - track.chi2layer = (*chi2Map)[cands->ptrAt(i)] / hit_pattern.trackerLayersWithMeasurement(); - if (cuts->isGoodTrack(track)) { - FillTracker(track, bestvz, bin); + for (unsigned int i = 0, n = cands->size(); i < n; ++i) { + track_ = {}; + track_.centbin = cbin; + const pat::PackedCandidate &pf = (*cands)[i]; + track_.et = pf.et(); + track_.eta = pf.eta(); + track_.phi = pf.phi(); + track_.pdgid = pf.pdgId(); + if ((idx == 1) and cuts_.isGoodHF(track_)) { + fillHF(track_, bestvz, bin); + } + if (!pf.hasTrackDetails()) + continue; + const reco::Track &trk = pf.pseudoTrack(); + track_.highPurity = pf.trackHighPurity(); + track_.charge = trk.charge(); + if (!track_.highPurity || track_.charge == 0) + continue; + track_.collection = idx; + track_.eta = trk.eta(); + track_.phi = trk.phi(); + track_.pt = trk.pt(); + track_.ptError = trk.ptError(); + track_.numberOfValidHits = trk.numberOfValidHits(); + track_.algos = trk.algo(); + track_.dz = std::abs(trk.dz(bestvtx)); + track_.dxy = std::abs(trk.dxy(bestvtx)); + track_.dzError = std::hypot(trk.dzError(), bestvzError); + track_.dxyError = trk.dxyError(bestvtx, vtx_cov); + track_.dzSig = track_.dz / track_.dzError; + track_.dxySig = track_.dxy / track_.dxyError; + const reco::HitPattern &hit_pattern = trk.hitPattern(); + track_.normalizedChi2 = (*chi2Map)[cands->ptrAt(i)]; + track_.chi2layer = (*chi2Map)[cands->ptrAt(i)] / hit_pattern.trackerLayersWithMeasurement(); + if (cuts_.isGoodTrack(track_)) { + fillTracker(track_, bestvz, bin); + } } } - } else { //calorimetry part if (bScalo_particleFlow_) { iEvent.getByToken(caloTokenPF_, calocands); - if (cands.isValid()) { - for (unsigned int i = 0, n = calocands->size(); i < n; ++i) { - track = {}; - track.centbin = cbin; - const reco::PFCandidate &pf = (*calocands)[i]; - track.et = pf.et(); - track.eta = pf.eta(); - track.phi = pf.phi(); - track.pdgid = pf.pdgId(); - if (cuts->isGoodHF(track)) { - fillHF(track, bestvz, bin); - } + for (unsigned int i = 0, n = calocands->size(); i < n; ++i) { + track_ = {}; + track_.centbin = cbin; + const reco::PFCandidate &pf = (*calocands)[i]; + track_.et = pf.et(); + track_.eta = pf.eta(); + track_.phi = pf.phi(); + track_.pdgid = pf.pdgId(); + if (cuts_.isGoodHF(track_)) { + fillHF(track_, bestvz, bin); } } } else { iEvent.getByToken(caloToken_, caloCollection_); - if (caloCollection_.isValid()) { - for (CaloTowerCollection::const_iterator j = caloCollection_->begin(); j != caloCollection_->end(); j++) { - track.eta = j->eta(); - track.phi = j->phi(); - track.et = j->emEt() + j->hadEt(); - track.pdgid = 1; - if (cuts->isGoodHF(track)) - fillHF(track, bestvz, bin); - } + for (const auto &tower : *caloCollection_) { + track_.eta = tower.eta(); + track_.phi = tower.phi(); + track_.et = tower.emEt() + tower.hadEt(); + track_.pdgid = 1; + if (cuts_.isGoodHF(track_)) + fillHF(track_, bestvz, bin); } } //Castor part iEvent.getByToken(castorToken_, castorCollection_); - if (castorCollection_.isValid()) { - for (std::vector::const_iterator j = castorCollection_->begin(); j != castorCollection_->end(); - j++) { - track.eta = j->eta(); - track.phi = j->phi(); - track.et = j->et(); - track.pdgid = 1; - if (cuts->isGoodCastor(track)) - FillCastor(track, bestvz, bin); - } + for (const auto &tower : *castorCollection_) { + track_.eta = tower.eta(); + track_.phi = tower.phi(); + track_.et = tower.et(); + track_.pdgid = 1; + if (cuts_.isGoodCastor(track_)) + fillCastor(track_, bestvz, bin); } //Tracking part iEvent.getByToken(trackToken_, trackCollection_); - if (trackCollection_.isValid()) { - for (reco::TrackCollection::const_iterator j = trackCollection_->begin(); j != trackCollection_->end(); j++) { - track.highPurity = j->quality(reco::TrackBase::highPurity); - track.charge = j->charge(); - if (!track.highPurity || track.charge == 0) - continue; - track.centbin = cbin; - track.collection = 0; - track.eta = j->eta(); - track.phi = j->phi(); - track.pt = j->pt(); - track.ptError = j->ptError(); - track.numberOfValidHits = j->numberOfValidHits(); - track.algos = j->algo(); - track.dz = j->dz(bestvtx); - track.dxy = -1. * j->dxy(bestvtx); - track.dzError = sqrt(pow(j->dzError(), 2) + pow(bestvzError, 2)); - track.dxyError = sqrt(pow(j->d0Error(), 2) + bestvxError * bestvyError); - track.dzSig = track.dz / track.dzError; - track.dxySig = track.dxy / track.dxyError; - track.normalizedChi2 = j->normalizedChi2(); - track.chi2layer = track.normalizedChi2 / j->hitPattern().trackerLayersWithMeasurement(); - - if (cuts->isGoodTrack(track)) - FillTracker(track, bestvz, bin); - } + for (const auto &trk : *trackCollection_) { + track_.highPurity = trk.quality(reco::TrackBase::highPurity); + track_.charge = trk.charge(); + if (!track_.highPurity || track_.charge == 0) + continue; + track_.centbin = cbin; + track_.collection = 0; + track_.eta = trk.eta(); + track_.phi = trk.phi(); + track_.pt = trk.pt(); + track_.ptError = trk.ptError(); + track_.numberOfValidHits = trk.numberOfValidHits(); + track_.algos = trk.algo(); + track_.dz = std::abs(trk.dz(bestvtx)); + track_.dxy = std::abs(trk.dxy(bestvtx)); + track_.dzError = std::hypot(trk.dzError(), bestvzError); + track_.dxyError = trk.dxyError(bestvtx, vtx_cov); + track_.dzSig = track_.dz / track_.dzError; + track_.dxySig = track_.dxy / track_.dxyError; + track_.normalizedChi2 = trk.normalizedChi2(); + track_.chi2layer = track_.normalizedChi2 / trk.hitPattern().trackerLayersWithMeasurement(); + if (cuts_.isGoodTrack(track_)) + fillTracker(track_, bestvz, bin); } } diff --git a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc index 711cca39a6107..0353e14e28d97 100644 --- a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc +++ b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc @@ -210,8 +210,7 @@ void HiEvtPlaneFlatProducer::produce(edm::Event& iEvent, const edm::EventSetup& //Get Event Planes // - edm::Handle evtPlanes_; - iEvent.getByToken(inputPlanesToken_, evtPlanes_); + auto const& evtPlanes = iEvent.get(inputPlanesToken_); auto evtplaneOutput = std::make_unique(); EvtPlane* ep[NumEPNames]; @@ -219,7 +218,7 @@ void HiEvtPlaneFlatProducer::produce(edm::Event& iEvent, const edm::EventSetup& ep[i] = nullptr; } int indx = 0; - for (auto&& rp : (*evtPlanes_)) { + for (auto&& rp : (evtPlanes)) { double s = rp.sumSin(0); double c = rp.sumCos(0); uint m = rp.mult(); From 695650005225b0e3ac4fd2cbfceb4fa77d57f083 Mon Sep 17 00:00:00 2001 From: Quan Wang Date: Tue, 20 Oct 2020 05:56:54 +0200 Subject: [PATCH 4/5] [+] update --- RecoHI/HiEvtPlaneAlgos/BuildFile.xml | 2 +- RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h | 2 +- RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc | 6 +++--- RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc | 3 --- 4 files changed, 5 insertions(+), 8 deletions(-) diff --git a/RecoHI/HiEvtPlaneAlgos/BuildFile.xml b/RecoHI/HiEvtPlaneAlgos/BuildFile.xml index 847d75c7055d8..133149fb69b15 100644 --- a/RecoHI/HiEvtPlaneAlgos/BuildFile.xml +++ b/RecoHI/HiEvtPlaneAlgos/BuildFile.xml @@ -2,7 +2,7 @@ - + diff --git a/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h b/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h index 65ba8b014a472..b0e91e867634d 100644 --- a/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h +++ b/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h @@ -120,7 +120,7 @@ namespace hi { int nHits = track.numberOfValidHits; if (track.ptError > pterror_ * track.pt) return false; - if (track.pt < 2.4 and (nHits == 3 or nHits == 4 or nHits == 5 or nHits == 6)) + if (track.pt < 2.4 and (nHits <= 6)) bPix = true; if (not bPix) { if (nHits < numberOfValidHits_) diff --git a/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc index dfbfcb5997397..4e7e95f477406 100644 --- a/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc +++ b/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc @@ -498,7 +498,7 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup for (int i = 0; i < NumEPNames; i++) rp[i]->reset(); edm::Handle> chi2Map; - edm::Handle> cands; + edm::Handle cands; edm::Handle calocands; if (bStrack_packedPFCandidates_) { for (int idx = 1; idx < 3; idx++) { @@ -542,8 +542,8 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup track_.dzSig = track_.dz / track_.dzError; track_.dxySig = track_.dxy / track_.dxyError; const reco::HitPattern &hit_pattern = trk.hitPattern(); - track_.normalizedChi2 = (*chi2Map)[cands->ptrAt(i)]; - track_.chi2layer = (*chi2Map)[cands->ptrAt(i)] / hit_pattern.trackerLayersWithMeasurement(); + track_.normalizedChi2 = (*chi2Map)[pat::PackedCandidateRef(cands, i)]; + track_.chi2layer = (*chi2Map)[pat::PackedCandidateRef(cands, i)] / hit_pattern.trackerLayersWithMeasurement(); if (cuts_.isGoodTrack(track_)) { fillTracker(track_, bestvz, bin); } diff --git a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc index 0353e14e28d97..d6f642dd64038 100644 --- a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc +++ b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc @@ -202,9 +202,6 @@ void HiEvtPlaneFlatProducer::produce(edm::Event& iEvent, const edm::EventSetup& double bestvz = -999.9; const reco::Vertex& vtx = iEvent.get(vertexToken_)[0]; bestvz = vtx.z(); - // Produce the EP regardless of vz position - // if (bestvz < minvtx_ || bestvz > maxvtx_) - // return; // //Get Event Planes From cce57527aae44680bf3098db7b488f7164570015 Mon Sep 17 00:00:00 2001 From: Quan Wang Date: Wed, 21 Oct 2020 22:45:46 +0200 Subject: [PATCH 5/5] [!] bug fix and code clean --- .../HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py | 7 ++-- .../python/hiEvtPlaneFlat_cfi.py | 2 -- .../HiEvtPlaneAlgos/src/EvtPlaneProducer.cc | 33 +++++++------------ .../src/HiEvtPlaneFlatProducer.cc | 4 --- 4 files changed, 14 insertions(+), 32 deletions(-) diff --git a/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py b/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py index 5f77961194ebf..f4cc2524bd6f5 100644 --- a/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py +++ b/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py @@ -16,8 +16,6 @@ maxet = cms.double(-1), minpt = cms.double(0.3), maxpt = cms.double(3.0), - minvtx = cms.double(-25.), - maxvtx = cms.double(25.), flatnvtxbins = cms.int32(10), flatminvtx = cms.double(-15.0), flatdelvtx = cms.double(3.0), @@ -41,11 +39,10 @@ (pp_on_AA_2018 | pp_on_PbPb_run3).toModify(hiEvtPlane, vertexTag = "offlinePrimaryVertices", - trackTag = "generalTracks", + trackTag = "packedPFCandidates", + caloTag = "particleFlow", minet = 0.01, minpt = 0.5, - minvtx = -15., - maxvtx = 15., dzdzerror_pix = 40., caloCentRef = -1, caloCentRefWidth = -1, diff --git a/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py b/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py index 29cc0f551572e..8d19848a88dea 100644 --- a/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py +++ b/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py @@ -10,8 +10,6 @@ trackTag = cms.InputTag("generalTracks"), FlatOrder = cms.int32(9), NumFlatBins = cms.int32(40), - minvtx = cms.double(-15.), - maxvtx = cms.double(15.), flatnvtxbins = cms.int32(10), flatminvtx = cms.double(-15.0), flatdelvtx = cms.double(3.0), diff --git a/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc index 4e7e95f477406..2a7416577a59b 100644 --- a/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc +++ b/RecoHI/HiEvtPlaneAlgos/src/EvtPlaneProducer.cc @@ -119,8 +119,6 @@ namespace hi { double &PtOrEt2, uint &epmult) { ang = -10; - sv = 0; - cv = 0; sv = sumsin; cv = sumcos; svNoWgt = sumsinNoWgt; @@ -219,8 +217,6 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { double maxet_; double minpt_; double maxpt_; - double minvtx_; - double maxvtx_; int flatnvtxbins_; double flatminvtx_; double flatdelvtx_; @@ -251,14 +247,14 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { for (int i = 0; i < NumEPNames; i++) { if (EPDet[i] != HF) continue; - if (track.et < minet) - continue; - if (track.et > maxet) - continue; if (minet_ < 0) minet = minTransverse[i]; if (maxet_ < 0) maxet = maxTransverse[i]; + if (track.et < minet) + continue; + if (track.et > maxet) + continue; if (not passEta(track.eta, i)) continue; double w = track.et; @@ -278,14 +274,14 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { double maxet = maxet_; for (int i = 0; i < NumEPNames; i++) { if (EPDet[i] == Castor) { - if (track.et < minet) - continue; - if (track.et > maxet) - continue; if (minet_ < 0) minet = minTransverse[i]; if (maxet_ < 0) maxet = maxTransverse[i]; + if (track.et < minet) + continue; + if (track.et > maxet) + continue; if (not passEta(track.eta, i)) continue; double w = track.et; @@ -321,14 +317,14 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> { double maxpt = maxpt_; for (int i = 0; i < NumEPNames; i++) { if (EPDet[i] == Tracker) { - if (track.pt < minpt) - continue; - if (track.pt > maxpt) - continue; if (minpt_ < 0) minpt = minTransverse[i]; if (maxpt_ < 0) maxpt = maxTransverse[i]; + if (track.pt < minpt) + continue; + if (track.pt > maxpt) + continue; if (not passEta(track.eta, i)) continue; double w = track.pt; @@ -360,8 +356,6 @@ EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig) maxet_(iConfig.getParameter("maxet")), minpt_(iConfig.getParameter("minpt")), maxpt_(iConfig.getParameter("maxpt")), - minvtx_(iConfig.getParameter("minvtx")), - maxvtx_(iConfig.getParameter("maxvtx")), flatnvtxbins_(iConfig.getParameter("flatnvtxbins")), flatminvtx_(iConfig.getParameter("flatminvtx")), flatdelvtx_(iConfig.getParameter("flatdelvtx")), @@ -491,9 +485,6 @@ void EvtPlaneProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSetup double bestvzError = vtx.zError(); math::XYZPoint bestvtx(bestvx, bestvy, bestvz); math::Error<3>::type vtx_cov = vtx.covariance(); - // Produce the EP regardless of vz position - // if (bestvz < minvtx_ || bestvz > maxvtx_) - // return; for (int i = 0; i < NumEPNames; i++) rp[i]->reset(); diff --git a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc index d6f642dd64038..756fe2526a8b8 100644 --- a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc +++ b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc @@ -81,8 +81,6 @@ class HiEvtPlaneFlatProducer : public edm::stream::EDProducer<> { const int FlatOrder_; int NumFlatBins_; - double minvtx_; - double maxvtx_; int flatnvtxbins_; double flatminvtx_; double flatdelvtx_; @@ -105,8 +103,6 @@ HiEvtPlaneFlatProducer::HiEvtPlaneFlatProducer(const edm::ParameterSet& iConfig) trackTag_(iConfig.getParameter("trackTag")), FlatOrder_(iConfig.getParameter("FlatOrder")), NumFlatBins_(iConfig.getParameter("NumFlatBins")), - minvtx_(iConfig.getParameter("minvtx")), - maxvtx_(iConfig.getParameter("maxvtx")), flatnvtxbins_(iConfig.getParameter("flatnvtxbins")), flatminvtx_(iConfig.getParameter("flatminvtx")), flatdelvtx_(iConfig.getParameter("flatdelvtx")),