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..c64983c9ce861 100644
--- a/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py
+++ b/PhysicsTools/PatAlgos/python/slimming/slimming_cff.py
@@ -67,7 +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
-(pp_on_AA_2018 | pp_on_PbPb_run3).toReplaceWith(slimmingTask, cms.Task(slimmingTask.copy(), hiPixelTracks))
+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,
+hiEvtPlane, hiEvtPlaneFlat))
from Configuration.Eras.Modifier_pp_on_PbPb_run3_cff import pp_on_PbPb_run3
from PhysicsTools.PatAlgos.packedCandidateMuonID_cfi import packedCandidateMuonID
@@ -76,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 0900da4ec6f0e..133149fb69b15 100644
--- a/RecoHI/HiEvtPlaneAlgos/BuildFile.xml
+++ b/RecoHI/HiEvtPlaneAlgos/BuildFile.xml
@@ -2,6 +2,7 @@
+
diff --git a/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h b/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h
new file mode 100644
index 0000000000000..b0e91e867634d
--- /dev/null
+++ b/RecoHI/HiEvtPlaneAlgos/interface/EPCuts.h
@@ -0,0 +1,168 @@
+#ifndef RecoHI_HiEvtPlaneAlgos_EPCuts_h
+#define RecoHI_HiEvtPlaneAlgos_EPCuts_h
+
+namespace hi {
+
+ enum class EP_ERA { ppReco, HIReco, Pixel, GenMC };
+
+ struct TrackStructure {
+ int centbin;
+ float eta;
+ float phi;
+ float et;
+ float pt;
+ int charge;
+ int pdgid;
+ int hits;
+ int algos;
+ int collection;
+ float dz;
+ float dxy;
+ float dzError;
+ float dxyError;
+ float ptError;
+ bool highPurity;
+ float dzSig;
+ float dxySig;
+ float normalizedChi2;
+ float dzError_Pix;
+ float chi2layer;
+ int numberOfValidHits;
+ int pixel;
+ };
+
+ class EPCuts {
+ 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;
+ }
+
+ bool isGoodHF(const TrackStructure& track) const {
+ if (track.pdgid != 1 && track.pdgid != 2)
+ return false;
+ if (std::abs(track.eta) < 3 || std::abs(track.eta) > 5)
+ return false;
+ return true;
+ }
+
+ bool isGoodCastor(const TrackStructure& track) const { return true; }
+
+ bool isGoodTrack(const TrackStructure& track) const {
+ 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(const TrackStructure& track) const {
+ if (track.charge == 0)
+ return false;
+ if (!track.highPurity)
+ return false;
+ if (track.ptError > pterror_ * track.pt)
+ return false;
+ if (track.numberOfValidHits < numberOfValidHits_)
+ return false;
+ if (track.chi2layer > chi2perlayer_)
+ return false;
+ if (std::abs(track.dxy) > dxyerror_ * track.dxyError)
+ return false;
+ if (std::abs(track.dz) > dzerror_ * track.dzError)
+ return false;
+ return true;
+ }
+
+ 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 > pterror_ * track.pt)
+ return false;
+ if (std::abs(track.dxy) > dxyerror_ * track.dxyError)
+ return false;
+ if (std::abs(track.dz) > dzerror_ * track.dzError)
+ return false;
+ if (track.chi2layer > chi2perlayer_)
+ return false;
+ //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(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 > pterror_ * track.pt)
+ return false;
+ if (track.pt < 2.4 and (nHits <= 6))
+ bPix = true;
+ if (not bPix) {
+ if (nHits < numberOfValidHits_)
+ return false;
+ if (track.chi2layer > chi2perlayer_)
+ return false;
+ if (track.ptError > pterror_ * track.pt)
+ return false;
+ int algo = track.algos;
+ if (track.pt > 2.4 && algo != reco::TrackBase::initialStep && algo != reco::TrackBase::lowPtTripletStep &&
+ algo != reco::TrackBase::pixelPairStep && algo != reco::TrackBase::detachedTripletStep)
+ return false;
+ if (std::abs(track.dxy) > dxyerror_ * track.dxyError)
+ return false;
+ if (std::abs(track.dz) > dzerror_ * track.dzError)
+ return false;
+ } else {
+ if (track.chi2layer > chi2Pix_)
+ return false;
+ if (std::abs(track.dz) > dzerror_Pix_ * track.dzError)
+ return false;
+ }
+ return true;
+ }
+
+ bool trackQuality_GenMC(const TrackStructure& track) const {
+ if (track.charge == 0)
+ return false;
+ if (std::abs(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 0e27e8a0b8342..01ca2ca215f24 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_;
@@ -65,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;
@@ -79,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;
@@ -99,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;
@@ -108,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;
@@ -117,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;
@@ -130,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);
@@ -140,20 +151,21 @@ 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;
if (ptval > 0)
- return pt * scale - pt2val / ptval;
+ return rv;
}
return 0.;
}
@@ -161,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));
@@ -172,23 +184,24 @@ class HiEvtPlaneFlatten {
return psi;
}
- double getSoffset(double s, double vtx, int centbin) const {
- int indx = getOffsetIndx(centbin, vtx);
- if (indx >= 0)
+ double soffset(double s, double vtx, int centbin) const {
+ int indx = offsetIndx(centbin, vtx);
+ if (indx >= 0) {
return s - yoffDB_[indx];
- else
+ } else {
return s;
+ }
}
- 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.;
@@ -197,6 +210,32 @@ class HiEvtPlaneFlatten {
return psi;
}
+ float minCent(int indx) const {
+ int ibin = indx / nvtxbins_;
+ return ibin * 100. / nbins_;
+ }
+
+ float maxCent(int indx) const {
+ int ibin = indx / nvtxbins_;
+ return (ibin + 1) * 100. / nbins_;
+ }
+
+ double minVtx(int indx) const {
+ int ivtx = indx - nvtxbins_ * (int)(indx / nvtxbins_);
+ return minvtx_ + ivtx * delvtx_;
+ }
+
+ double maxVtx(int indx) const {
+ int ivtx = indx - nvtxbins_ * (int)(indx / nvtxbins_);
+ return minvtx_ + (ivtx + 1) * delvtx_;
+ }
+
+ 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);
+ }
+
~HiEvtPlaneFlatten() {}
int getHBins() const { return hbins_; }
int getOBins() const { return obins_; }
@@ -421,9 +460,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 +506,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..ff634363248e9 100644
--- a/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h
+++ b/RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h
@@ -2,127 +2,111 @@
#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"};
-
- enum Detectors { Tracker, HF, Castor };
-
- 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 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 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 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 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 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 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 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 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 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 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 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};
-
- static const int NumEPNames = 29;
+ 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 std::array EPDet = {{
+ HF, HF, HF, Tracker, Tracker,
+ Tracker, HF, HF, HF, Tracker,
+ Tracker, Tracker
+ }};
+
+ const std::array EPOrder = {{
+ 2, 2, 2, 2, 2,
+ 2, 3, 3, 3, 3,
+ 3, 3
+ }};
+
+ 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 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 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 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 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 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::array ResCalcType = {{
+ "3sub", "3sub", "3sub", "3sub", "3sub",
+ "3sub", "3sub", "3sub", "3sub", "3sub",
+ "3sub", "3sub"
+ }};
+
+ const std::array MomConsWeight = {{
+ "no", "no", "no", "no", "no",
+ "no", "no", "no", "no", "no",
+ "no", "no"
+ }};
+
+ const std::array RCMate1 = {{
+ HFp2, HFm2, trackm2, HFm2, HFp2,
+ HFp2, HFp3, HFm3, trackm3, HFm3,
+ HFp3, HFp3
+ }};
+
+ const std::array RCMate2 = {{
+ trackmid2, trackmid2, trackp2, HFp2, HFm2,
+ HFm2, trackmid3, trackmid3, trackp3, HFp3,
+ HFm3, HFm3
+ }};
+
+ // 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..f4cc2524bd6f5 100644
--- a/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py
+++ b/RecoHI/HiEvtPlaneAlgos/python/HiEvtPlane_cfi.py
@@ -1,25 +1,50 @@
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 = "offlinePrimaryVertices",
+ trackTag = "packedPFCandidates",
+ caloTag = "particleFlow",
+ minet = 0.01,
+ minpt = 0.5,
+ dzdzerror_pix = 40.,
+ caloCentRef = -1,
+ caloCentRefWidth = -1,
+ cutEra = 0
+)
diff --git a/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py b/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py
index 33ac26edf1bdb..8d19848a88dea 100644
--- a/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py
+++ b/RecoHI/HiEvtPlaneAlgos/python/hiEvtPlaneFlat_cfi.py
@@ -1,19 +1,20 @@
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),
+ 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..2a7416577a59b 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;
@@ -116,8 +119,6 @@ namespace hi {
double &PtOrEt2,
uint &epmult) {
ang = -10;
- sv = 0;
- cv = 0;
sv = sumsin;
cv = sumcos;
svNoWgt = sumsinNoWgt;
@@ -175,49 +176,169 @@ class EvtPlaneProducer : public edm::stream::EDProducer<> {
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::InputTag castorTag_;
- edm::EDGetTokenT> castorToken;
+ edm::EDGetTokenT> castorToken_;
edm::Handle> castorCollection_;
edm::InputTag trackTag_;
- edm::EDGetTokenT trackToken;
+ edm::EDGetTokenT trackToken_;
+ edm::InputTag losttrackTag_;
edm::Handle trackCollection_;
+ bool bStrack_packedPFCandidates_;
+ bool bScalo_particleFlow_;
+ edm::EDGetTokenT packedToken_;
+ edm::EDGetTokenT lostToken_;
- edm::ESWatcher hiWatcher;
- edm::ESWatcher hirpWatcher;
+ edm::InputTag chi2MapTag_;
+ edm::EDGetTokenT> chi2MapToken_;
+ edm::InputTag chi2MapLostTag_;
+ edm::EDGetTokenT> chi2MapLostToken_;
bool loadDB_;
double minet_;
double maxet_;
double minpt_;
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];
+ TrackStructure track_;
+
+ edm::ESWatcher hiWatcher_;
+ edm::ESWatcher hirpWatcher_;
+
+ 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 (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;
+ if (loadDB_)
+ 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);
+ }
+ }
+ rp[i]->addParticle(w, track.et, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta);
+ }
+ };
+
+ 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 (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;
+ if (EPOrder[i] == 1) {
+ if (MomConsWeight[i][0] == 'y' && loadDB_) {
+ w = flat[i]->getW(track.et, vz, bin);
+ }
+ }
+ rp[i]->addParticle(w, track.et, sin(EPOrder[i] * track.phi), cos(EPOrder[i] * track.phi), track.eta);
+ }
+ }
+ }
+
+ 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 (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;
+ 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);
+ }
+ }
+ 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,20 +348,36 @@ 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")),
minpt_(iConfig.getParameter("minpt")),
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"))
+
+{
+ 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")) {
@@ -248,15 +385,27 @@ EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet &iConfig)
}
centralityLabel_ = centralityVariable_ + centralityMC_;
- centralityBinToken = consumes(centralityBinTag_);
+ centralityBinToken_ = consumes(centralityBinTag_);
- vertexToken = consumes>(vertexTag_);
+ vertexToken_ = consumes>(vertexTag_);
- caloToken = consumes(caloTag_);
+ 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_);
- castorToken = consumes>(castorTag_);
-
- trackToken = consumes(trackTag_);
+ } else {
+ if (bScalo_particleFlow_) {
+ caloTokenPF_ = consumes(caloTag_);
+ } else {
+ caloToken_ = consumes(caloTag_);
+ }
+ castorToken_ = consumes>(castorTag_);
+ trackToken_ = consumes(trackTag_);
+ }
produces();
for (int i = 0; i < NumEPNames; i++) {
@@ -264,7 +413,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 +434,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 +453,157 @@ 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 = iEvent.get(centralityBinToken_);
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;
- //
+ //best vertex
+ const reco::Vertex &vtx = iEvent.get(vertexToken_)[0];
+ double bestvz = vtx.z();
+ double bestvx = vtx.x();
+ double bestvy = vtx.y();
+ double bestvzError = vtx.zError();
+ math::XYZPoint bestvtx(bestvx, bestvy, bestvz);
+ math::Error<3>::type vtx_cov = vtx.covariance();
+
for (int i = 0; i < NumEPNames; i++)
rp[i]->reset();
- if (vzr_sell < minvtx_ or vzr_sell > maxvtx_)
- return;
-
- //calorimetry part
-
- double tower_eta, tower_phi;
- double tower_energyet, tower_energyet_e, tower_energyet_h;
-
- iEvent.getByToken(caloToken, caloCollection_);
-
- 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)
+ edm::Handle> chi2Map;
+ edm::Handle cands;
+ edm::Handle calocands;
+ if (bStrack_packedPFCandidates_) {
+ for (int idx = 1; idx < 3; idx++) {
+ if (idx == 1) {
+ iEvent.getByToken(packedToken_, cands);
+ iEvent.getByToken(chi2MapToken_, chi2Map);
+ }
+ if (idx == 2) {
+ iEvent.getByToken(lostToken_, cands);
+ iEvent.getByToken(chi2MapLostToken_, 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 ((idx == 1) and cuts_.isGoodHF(track_)) {
+ fillHF(track_, bestvz, bin);
+ }
+ if (!pf.hasTrackDetails())
continue;
- if (tower_energyet > maxet)
+ const reco::Track &trk = pf.pseudoTrack();
+ track_.highPurity = pf.trackHighPurity();
+ track_.charge = trk.charge();
+ if (!track_.highPurity || track_.charge == 0)
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);
+ 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)[pat::PackedCandidateRef(cands, i)];
+ track_.chi2layer = (*chi2Map)[pat::PackedCandidateRef(cands, 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;
- }
- rp[i]->addParticle(w, tower_energyet, sin(EPOrder[i] * tower_phi), cos(EPOrder[i] * tower_phi), tower_eta);
+ } else {
+ //calorimetry part
+ if (bScalo_particleFlow_) {
+ iEvent.getByToken(caloTokenPF_, calocands);
+ 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_);
+ 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);
+ }
}
- }
-
- //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;
- }
- 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);
- }
- }
- }
- } //end for
+ //Castor part
+ iEvent.getByToken(castorToken_, castorCollection_);
+ 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_);
+ 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);
+ }
}
auto evtplaneOutput = std::make_unique();
diff --git a/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc b/RecoHI/HiEvtPlaneAlgos/src/HiEvtPlaneFlatProducer.cc
index b2a106e0a9ed0..756fe2526a8b8 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"
@@ -29,17 +30,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
@@ -60,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;
@@ -80,15 +81,19 @@ class HiEvtPlaneFlatProducer : public edm::stream::EDProducer<> {
const int FlatOrder_;
int NumFlatBins_;
+ 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_;
};
-
+//
+// constructors and destructor
+//
HiEvtPlaneFlatProducer::HiEvtPlaneFlatProducer(const edm::ParameterSet& iConfig)
: centralityVariable_(iConfig.getParameter("centralityVariable")),
centralityBinTag_(iConfig.getParameter("centralityBinTag")),
@@ -98,66 +103,68 @@ HiEvtPlaneFlatProducer::HiEvtPlaneFlatProducer(const edm::ParameterSet& iConfig)
trackTag_(iConfig.getParameter("trackTag")),
FlatOrder_(iConfig.getParameter("FlatOrder")),
NumFlatBins_(iConfig.getParameter("NumFlatBins")),
+ 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")) {
+ nCentBins_ = 200.;
+
if (iConfig.exists("nonDefaultGlauberModel")) {
centralityMC_ = iConfig.getParameter("nonDefaultGlauberModel");
}
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();
//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,62 +174,59 @@ 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;
+ cbin = iEvent.get(centralityBinToken_);
+ bin = cbin / CentBinCompression_;
//
//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;
+
+ //best vertex
+ double bestvz = -999.9;
+ const reco::Vertex& vtx = iEvent.get(vertexToken_)[0];
+ bestvz = vtx.z();
//
//Get Event Planes
//
+ auto const& evtPlanes = iEvent.get(inputPlanesToken_);
+
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);
+ for (auto&& rp : (evtPlanes)) {
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]->soffset(s, bestvz, bin);
+ coff = flat[indx]->coffset(c, bestvz, bin);
+ psiOffset = flat[indx]->offsetPsi(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));
@@ -231,7 +235,7 @@ void HiEvtPlaneFlatProducer::produce(edm::Event& iEvent, const edm::EventSetup&
++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]);
}