CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
HiFJGridEmptyAreaCalculator.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: HiJetBackground/HiFJGridEmptyAreaCalculator
4 // Class: HiFJGridEmptyAreaCalculator
5 // Based on: fastjet/tools/GridMedianBackgroundEstimator
6 //
7 // Original Author: Doga Gulhan
8 // Created: Wed Mar 16 14:00:04 CET 2016
9 //
10 //
11 
13 using namespace std;
14 using namespace edm;
15 
17  : gridWidth_(iConfig.getParameter<double>("gridWidth")),
18  band_(iConfig.getParameter<double>("bandWidth")),
19  hiBinCut_(iConfig.getParameter<int>("hiBinCut")),
20  doCentrality_(iConfig.getParameter<bool>("doCentrality")),
21  keepGridInfo_(iConfig.getParameter<bool>("keepGridInfo")) {
22  ymin_ = -99;
23  ymax_ = -99;
24  dy_ = -99;
25  dphi_ = -99;
26  tileArea_ = -99;
27 
28  dyJet_ = 99;
29  yminJet_ = -99;
30  ymaxJet_ = -99;
31  totalInboundArea_ = -99;
32  etaminJet_ = -99;
33  etamaxJet_ = -99;
34 
35  ny_ = 0;
36  nphi_ = 0;
37  ntotal_ = 0;
38 
39  ntotalJet_ = 0;
40  nyJet_ = 0;
41 
42  pfCandsToken_ = consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfCandSource"));
43  mapEtaToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapEtaEdges"));
44  mapRhoToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapToRho"));
45  mapRhoMToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapToRhoM"));
46  jetsToken_ = consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jetSource"));
47  centralityBinToken_ = consumes<int>(iConfig.getParameter<edm::InputTag>("CentralityBinSrc"));
48 
49  //register your products
50  produces<std::vector<double>>("mapEmptyCorrFac");
51  produces<std::vector<double>>("mapToRhoCorr");
52  produces<std::vector<double>>("mapToRhoMCorr");
53  produces<std::vector<double>>("mapToRhoCorr1Bin");
54  produces<std::vector<double>>("mapToRhoMCorr1Bin");
55  //rho calculation on a grid using median
56  if (keepGridInfo_) {
57  produces<std::vector<double>>("mapRhoVsEtaGrid");
58  produces<std::vector<double>>("mapMeanRhoVsEtaGrid");
59  produces<std::vector<double>>("mapEtaMaxGrid");
60  produces<std::vector<double>>("mapEtaMinGrid");
61  }
62 }
63 
65  // do anything here that needs to be done at desctruction time
66  // (e.g. close files, deallocate resources etc.)
67 }
68 
70  using namespace edm;
71  //Values from rho calculator
72  edm::Handle<std::vector<double>> mapEtaRanges;
73  iEvent.getByToken(mapEtaToken_, mapEtaRanges);
74 
76  iEvent.getByToken(mapRhoToken_, mapRho);
77 
79  iEvent.getByToken(mapRhoMToken_, mapRhoM);
80 
81  // if doCentrality is set to true calculate empty area correction
82  // for only events with hiBin > hiBinCut
83  int hiBin = -1;
84  bool doEmptyArea = true;
85  if (doCentrality_) {
86  edm::Handle<int> cbin;
87  iEvent.getByToken(centralityBinToken_, cbin);
88  hiBin = *cbin;
89 
90  if (hiBin < hiBinCut_)
91  doEmptyArea = false;
92  }
93 
94  //Define output vectors
95  int neta = (int)mapEtaRanges->size();
96 
97  auto mapToRhoCorrOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
98  auto mapToRhoMCorrOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
99  auto mapToRhoCorr1BinOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
100  auto mapToRhoMCorr1BinOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
101 
102  setupGrid(mapEtaRanges->at(0), mapEtaRanges->at(neta - 1));
103 
104  //calculate empty area correction over full acceptance leaving eta bands on the sides
105  double allAcceptanceCorr = 1;
106  if (doEmptyArea) {
107  etaminJet_ = mapEtaRanges->at(0) - band_;
108  etamaxJet_ = mapEtaRanges->at(neta - 1) + band_;
109 
111 
112  allAcceptanceCorr = totalInboundArea_;
113  }
114 
115  //calculate empty area correction in each eta range
116  for (int ieta = 0; ieta < (neta - 1); ieta++) {
117  double correctionKt = 1;
118  double rho = mapRho->at(ieta);
119  double rhoM = mapRhoM->at(ieta);
120 
121  if (doEmptyArea) {
122  double etamin = mapEtaRanges->at(ieta);
123  double etamax = mapEtaRanges->at(ieta + 1);
124 
125  etaminJet_ = etamin + band_;
126  etamaxJet_ = etamax - band_;
127 
129  correctionKt = totalInboundArea_;
130  }
131 
132  mapToRhoCorrOut->at(ieta) = correctionKt * rho;
133  mapToRhoMCorrOut->at(ieta) = correctionKt * rhoM;
134 
135  mapToRhoCorr1BinOut->at(ieta) = allAcceptanceCorr * rho;
136  mapToRhoMCorr1BinOut->at(ieta) = allAcceptanceCorr * rhoM;
137  }
138 
139  iEvent.put(std::move(mapToRhoCorrOut), "mapToRhoCorr");
140  iEvent.put(std::move(mapToRhoMCorrOut), "mapToRhoMCorr");
141  iEvent.put(std::move(mapToRhoCorr1BinOut), "mapToRhoCorr1Bin");
142  iEvent.put(std::move(mapToRhoMCorr1BinOut), "mapToRhoMCorr1Bin");
143 
144  //calculate rho from grid as a function of eta over full range using PF candidates
145 
146  auto mapRhoVsEtaGridOut = std::make_unique<std::vector<double>>(ny_, 0.);
147  auto mapMeanRhoVsEtaGridOut = std::make_unique<std::vector<double>>(ny_, 0.);
148  auto mapEtaMaxGridOut = std::make_unique<std::vector<double>>(ny_, 0.);
149  auto mapEtaMinGridOut = std::make_unique<std::vector<double>>(ny_, 0.);
150  calculateGridRho(iEvent, iSetup);
151  if (keepGridInfo_) {
152  for (int ieta = 0; ieta < ny_; ieta++) {
153  mapRhoVsEtaGridOut->at(ieta) = rhoVsEta_[ieta];
154  mapMeanRhoVsEtaGridOut->at(ieta) = meanRhoVsEta_[ieta];
155  mapEtaMaxGridOut->at(ieta) = etaMaxGrid_[ieta];
156  mapEtaMinGridOut->at(ieta) = etaMinGrid_[ieta];
157  }
158 
159  iEvent.put(std::move(mapRhoVsEtaGridOut), "mapRhoVsEtaGrid");
160  iEvent.put(std::move(mapMeanRhoVsEtaGridOut), "mapMeanRhoVsEtaGrid");
161  iEvent.put(std::move(mapEtaMaxGridOut), "mapEtaMaxGrid");
162  iEvent.put(std::move(mapEtaMinGridOut), "mapEtaMinGrid");
163  }
164 }
165 
166 //----------------------------------------------------------------------
167 // setting a new event
168 //----------------------------------------------------------------------
169 // tell the background estimator that it has a new event, composed
170 // of the specified particles.
172  vector<vector<double>> scalarPt(ny_, vector<double>(nphi_, 0.0));
173 
175  iEvent.getByToken(pfCandsToken_, pfCands);
176  const reco::PFCandidateCollection* pfCandidateColl = pfCands.product();
177  for (unsigned icand = 0; icand < pfCandidateColl->size(); icand++) {
178  const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
179  //use ony the particles within the eta range
180  if (pfCandidate.eta() < ymin_ || pfCandidate.eta() > ymax_)
181  continue;
182  int jeta = tileIndexEta(&pfCandidate);
183  int jphi = tileIndexPhi(&pfCandidate);
184  scalarPt[jeta][jphi] += pfCandidate.pt();
185  }
186 
187  rhoVsEta_.resize(ny_);
188  meanRhoVsEta_.resize(ny_);
189  for (int jeta = 0; jeta < ny_; jeta++) {
190  rhoVsEta_[jeta] = 0;
191  meanRhoVsEta_[jeta] = 0;
192  vector<double> rhoVsPhi;
193  int nEmpty = 0;
194 
195  for (int jphi = 0; jphi < nphi_; jphi++) {
196  double binpt = scalarPt[jeta][jphi];
197  meanRhoVsEta_[jeta] += binpt;
198  //fill in the vector for median calculation
199  if (binpt > 0)
200  rhoVsPhi.push_back(binpt);
201  else
202  nEmpty++;
203  }
204  meanRhoVsEta_[jeta] /= ((double)nphi_);
206 
207  //median calculation
208  sort(rhoVsPhi.begin(), rhoVsPhi.end());
209  //use only the nonzero grid cells for median calculation;
210  int nFull = nphi_ - nEmpty;
211  if (nFull == 0) {
212  rhoVsEta_[jeta] = 0;
213  continue;
214  }
215  if (nFull % 2 == 0) {
216  rhoVsEta_[jeta] = (rhoVsPhi[(int)(nFull / 2 - 1)] + rhoVsPhi[(int)(nFull / 2)]) / 2;
217  } else {
218  rhoVsEta_[jeta] = rhoVsPhi[(int)(nFull / 2)];
219  }
220  //correct for empty cells
221  rhoVsEta_[jeta] *= (((double)nFull) / ((double)nphi_));
222  //normalize to area
224  }
225 }
226 
229  iEvent.getByToken(jetsToken_, jets);
230 
231  //calculate jet kt area fraction inside boundary by grid
232  totalInboundArea_ = 0;
233 
234  for (auto jet = jets->begin(); jet != jets->end(); ++jet) {
235  if (jet->eta() < etaminJet_ || jet->eta() > etamaxJet_)
236  continue;
237 
238  double areaKt = jet->jetArea();
239  setupGridJet(&*jet);
240  std::vector<std::pair<int, int>> pfIndicesJet;
241  std::vector<std::pair<int, int>> pfIndicesJetInbound;
242  int nConstitJet = 0;
243  int nConstitJetInbound = 0;
244  for (auto daughter : jet->getJetConstituentsQuick()) {
245  auto pfCandidate = static_cast<const reco::PFCandidate*>(daughter);
246 
247  int jeta = tileIndexEtaJet(&*pfCandidate);
248  int jphi = tileIndexPhi(&*pfCandidate);
249  pfIndicesJet.push_back(std::make_pair(jphi, jeta));
250  nConstitJet++;
251  if (pfCandidate->eta() < etaminJet_ && pfCandidate->eta() > etamaxJet_)
252  continue;
253  pfIndicesJetInbound.push_back(std::make_pair(jphi, jeta));
254  nConstitJetInbound++;
255  }
256 
257  //if the jet is well within the eta range just add the area
258  if (nConstitJet == nConstitJetInbound) {
259  totalInboundArea_ += areaKt;
260  continue;
261  }
262 
263  //for jets that fall outside of eta range calculate fraction of area
264  //inside the range with a grid
265  int nthis = 0;
266  if (nConstitJet > 0)
267  nthis = numJetGridCells(pfIndicesJet);
268 
269  int nthisInbound = 0;
270  if (nConstitJetInbound > 0)
271  nthisInbound = numJetGridCells(pfIndicesJetInbound);
272 
273  double fractionArea = ((double)nthisInbound) / ((double)nthis);
274  totalInboundArea_ += areaKt * fractionArea;
275  }
276 
277  //divide by the total area in that range
279 
280  //the fraction can still be greater than 1 because kt area fraction inside
281  //the range can differ from what we calculated with the grid
282  if (totalInboundArea_ > 1)
283  totalInboundArea_ = 1;
284 }
285 
286 // #ifndef FASTJET_GMBGE_USEFJGRID
287 //----------------------------------------------------------------------
288 // protected material
289 //----------------------------------------------------------------------
290 // configure the grid
292  // since we've exchanged the arguments of the grid constructor,
293  // there's a danger of calls with exchanged ymax,spacing arguments --
294  // the following check should catch most such situations.
295  ymin_ = etamin;
296  ymax_ = etamax;
297 
299 
300  // this grid-definition code is becoming repetitive -- it should
301  // probably be moved somewhere central...
302  double nyDouble = (ymax_ - ymin_) / gridWidth_;
303  ny_ = int(nyDouble + 0.5);
304  dy_ = (ymax_ - ymin_) / ny_;
305 
306  nphi_ = int(twopi_ / gridWidth_ + 0.5);
307  dphi_ = twopi_ / nphi_;
308 
309  // some sanity checking (could throw a fastjet::Error)
310  assert(ny_ >= 1 && nphi_ >= 1);
311 
312  ntotal_ = nphi_ * ny_;
313  //_scalar_pt.resize(_ntotal);
314  tileArea_ = dy_ * dphi_;
315 
316  etaMaxGrid_.resize(ny_);
317  etaMinGrid_.resize(ny_);
318  for (int jeta = 0; jeta < ny_; jeta++) {
319  etaMinGrid_[jeta] = etamin + dy_ * ((double)jeta);
320  etaMaxGrid_[jeta] = etamin + dy_ * ((double)jeta + 1.);
321  }
322 }
323 
324 //----------------------------------------------------------------------
325 // retrieve the grid tile index for a given PseudoJet
327  // directly taking int does not work for values between -1 and 0
328  // so use floor instead
329  // double iy_double = (p.rap() - _ymin) / _dy;
330  // if (iy_double < 0.0) return -1;
331  // int iy = int(iy_double);
332  // if (iy >= _ny) return -1;
333 
334  // writing it as below gives a huge speed gain (factor two!). Even
335  // though answers are identical and the routine here is not the
336  // speed-critical step. It's not at all clear why.
337 
338  int iphi = int((pfCand->phi() + (twopi_ / 2.)) / dphi_);
339  assert(iphi >= 0 && iphi <= nphi_);
340  if (iphi == nphi_)
341  iphi = 0; // just in case of rounding errors
342 
343  return iphi;
344 }
345 
346 //----------------------------------------------------------------------
347 // retrieve the grid tile index for a given PseudoJet
349  // directly taking int does not work for values between -1 and 0
350  // so use floor instead
351  // double iy_double = (p.rap() - _ymin) / _dy;
352  // if (iy_double < 0.0) return -1;
353  // int iy = int(iy_double);
354  // if (iy >= _ny) return -1;
355 
356  // writing it as below gives a huge speed gain (factor two!). Even
357  // though answers are identical and the routine here is not the
358  // speed-critical step. It's not at all clear why.
359  int iy = int(floor((pfCand->eta() - ymin_) / dy_));
360  if (iy < 0 || iy >= ny_)
361  return -1;
362 
363  assert(iy < ny_ && iy >= 0);
364 
365  return iy;
366 }
367 
368 // #ifndef FASTJET_GMBGE_USEFJGRID
369 //----------------------------------------------------------------------
370 // protected material
371 //----------------------------------------------------------------------
372 // configure the grid
374  // since we've exchanged the arguments of the grid constructor,
375  // there's a danger of calls with exchanged ymax,spacing arguments --
376  // the following check should catch most such situations.
377  yminJet_ = jet->eta() - 0.6;
378  ymaxJet_ = jet->eta() + 0.6;
379 
381 
382  // this grid-definition code is becoming repetitive -- it should
383  // probably be moved somewhere central...
384  double nyDouble = (ymaxJet_ - yminJet_) / gridWidth_;
385  nyJet_ = int(nyDouble + 0.5);
386  dyJet_ = (ymaxJet_ - yminJet_) / nyJet_;
387 
388  assert(nyJet_ >= 1);
389 
390  ntotalJet_ = nphi_ * nyJet_;
391  //_scalar_pt.resize(_ntotal);
392 }
393 
394 //----------------------------------------------------------------------
395 // retrieve the grid tile index for a given PseudoJet
397  // directly taking int does not work for values between -1 and 0
398  // so use floor instead
399  // double iy_double = (p.rap() - _ymin) / _dy;
400  // if (iy_double < 0.0) return -1;
401  // int iy = int(iy_double);
402  // if (iy >= _ny) return -1;
403 
404  // writing it as below gives a huge speed gain (factor two!). Even
405  // though answers are identical and the routine here is not the
406  // speed-critical step. It's not at all clear why.
407  int iyjet = int(floor((pfCand->eta() - yminJet_) / dy_));
408  if (iyjet < 0 || iyjet >= nyJet_)
409  return -1;
410 
411  assert(iyjet < nyJet_ && iyjet >= 0);
412 
413  return iyjet;
414 }
415 
417  int ngrid = 0;
418  //sort phi eta grid indices in phi
419  std::sort(indices.begin(), indices.end());
420  int lowestJPhi = indices[0].first;
421  int lowestJEta = indices[0].second;
422  int highestJEta = lowestJEta;
423 
424  //for each fixed phi value calculate the number of grids in eta
425  for (unsigned int iconst = 1; iconst < indices.size(); iconst++) {
426  int jphi = indices[iconst].first;
427  int jeta = indices[iconst].second;
428  if (jphi == lowestJPhi) {
429  if (jeta < lowestJEta)
430  lowestJEta = jeta;
431  if (jeta > highestJEta)
432  highestJEta = jeta;
433  } else {
434  lowestJPhi = jphi;
435  ngrid += highestJEta - lowestJEta + 1;
436  lowestJEta = jeta;
437  highestJEta = jeta;
438  }
439  }
440  ngrid += highestJEta - lowestJEta + 1;
441  return ngrid;
442 }
443 
445 
447 
450  desc.add<edm::InputTag>("jetSource", edm::InputTag("kt4PFJets"));
451  desc.add<edm::InputTag>("CentralityBinSrc", edm::InputTag("centralityBin"));
452  desc.add<edm::InputTag>("mapEtaEdges", edm::InputTag("mapEtaEdges"));
453  desc.add<edm::InputTag>("mapToRho", edm::InputTag("mapToRho"));
454  desc.add<edm::InputTag>("mapToRhoM", edm::InputTag("mapToRhoM"));
455  desc.add<edm::InputTag>("pfCandSource", edm::InputTag("particleFlow"));
456  desc.add<double>("gridWidth", 0.05);
457  desc.add<double>("bandWidth", 0.2);
458  desc.add<bool>("doCentrality", true);
459  desc.add<int>("hiBinCut", 100);
460  desc.add<bool>("keepGridInfo", false);
461  descriptions.add("hiFJGridEmptyAreaCalculator", desc);
462 }
463 
edm::EDGetTokenT< std::vector< double > > mapRhoToken_
edm::EDGetTokenT< int > centralityBinToken_
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
int numJetGridCells(std::vector< std::pair< int, int >> &indices)
number of grid cells that overlap with jet constituents filling in the in between area ...
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
double ymin_
internal parameters for grid
double pt() const final
transverse momentum
edm::EDGetTokenT< edm::View< reco::Jet > > jetsToken_
input tokens
int tileIndexEta(const reco::PFCandidate *pfCand)
edm::EDGetTokenT< reco::PFCandidateCollection > pfCandsToken_
Base class for all types of Jets.
Definition: Jet.h:20
T const * product() const
Definition: Handle.h:70
void setupGrid(double eta_min, double eta_max)
configure the grid
edm::EDGetTokenT< std::vector< double > > mapRhoMToken_
int tileIndexEtaJet(const reco::PFCandidate *pfCand)
assert(be >=bs)
void setupGridJet(const reco::Jet *jet)
edm::EDGetTokenT< std::vector< double > > mapEtaToken_
int tileIndexPhi(const reco::PFCandidate *pfCand)
int iEvent
Definition: GenABIO.cc:224
const int neta
void calculateAreaFractionOfJets(const edm::Event &iEvent, const edm::EventSetup &iSetup)
const double twopi_
information about the grid
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
void beginStream(edm::StreamID) override
std::vector< reco::PFCandidate > PFCandidateCollection
collection of PFCandidates
void calculateGridRho(const edm::Event &iEvent, const edm::EventSetup &iSetup)
HiFJGridEmptyAreaCalculator(const edm::ParameterSet &)
void produce(edm::Event &, const edm::EventSetup &) override
void add(std::string const &label, ParameterSetDescription const &psetDescription)
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:41
HLT enums.
double phi() const final
momentum azimuthal angle
def move(src, dest)
Definition: eostools.py:511
double eta() const final
momentum pseudorapidity