Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:43

0001 // -*- C++ -*-
0002 //
0003 // Package:    HiJetBackground/HiFJGridEmptyAreaCalculator
0004 // Class:      HiFJGridEmptyAreaCalculator
0005 // Based on:   fastjet/tools/GridMedianBackgroundEstimator
0006 //
0007 // Original Author:  Doga Gulhan
0008 //         Created:  Wed Mar 16 14:00:04 CET 2016
0009 //
0010 //
0011 
0012 #include "RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.h"
0013 using namespace std;
0014 using namespace edm;
0015 
0016 HiFJGridEmptyAreaCalculator::HiFJGridEmptyAreaCalculator(const edm::ParameterSet& iConfig)
0017     : gridWidth_(iConfig.getParameter<double>("gridWidth")),
0018       band_(iConfig.getParameter<double>("bandWidth")),
0019       hiBinCut_(iConfig.getParameter<int>("hiBinCut")),
0020       doCentrality_(iConfig.getParameter<bool>("doCentrality")),
0021       keepGridInfo_(iConfig.getParameter<bool>("keepGridInfo")) {
0022   ymin_ = -99;
0023   ymax_ = -99;
0024   dy_ = -99;
0025   dphi_ = -99;
0026   tileArea_ = -99;
0027 
0028   dyJet_ = 99;
0029   yminJet_ = -99;
0030   ymaxJet_ = -99;
0031   totalInboundArea_ = -99;
0032   etaminJet_ = -99;
0033   etamaxJet_ = -99;
0034 
0035   ny_ = 0;
0036   nphi_ = 0;
0037   ntotal_ = 0;
0038 
0039   ntotalJet_ = 0;
0040   nyJet_ = 0;
0041 
0042   pfCandsToken_ = consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>("pfCandSource"));
0043   mapEtaToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapEtaEdges"));
0044   mapRhoToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapToRho"));
0045   mapRhoMToken_ = consumes<std::vector<double>>(iConfig.getParameter<edm::InputTag>("mapToRhoM"));
0046   jetsToken_ = consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jetSource"));
0047   centralityBinToken_ = consumes<int>(iConfig.getParameter<edm::InputTag>("CentralityBinSrc"));
0048 
0049   //register your products
0050   produces<std::vector<double>>("mapEmptyCorrFac");
0051   produces<std::vector<double>>("mapToRhoCorr");
0052   produces<std::vector<double>>("mapToRhoMCorr");
0053   produces<std::vector<double>>("mapToRhoCorr1Bin");
0054   produces<std::vector<double>>("mapToRhoMCorr1Bin");
0055   //rho calculation on a grid using median
0056   if (keepGridInfo_) {
0057     produces<std::vector<double>>("mapRhoVsEtaGrid");
0058     produces<std::vector<double>>("mapMeanRhoVsEtaGrid");
0059     produces<std::vector<double>>("mapEtaMaxGrid");
0060     produces<std::vector<double>>("mapEtaMinGrid");
0061   }
0062 }
0063 
0064 HiFJGridEmptyAreaCalculator::~HiFJGridEmptyAreaCalculator() {
0065   // do anything here that needs to be done at desctruction time
0066   // (e.g. close files, deallocate resources etc.)
0067 }
0068 
0069 void HiFJGridEmptyAreaCalculator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0070   using namespace edm;
0071   //Values from rho calculator
0072   edm::Handle<std::vector<double>> mapEtaRanges;
0073   iEvent.getByToken(mapEtaToken_, mapEtaRanges);
0074 
0075   edm::Handle<std::vector<double>> mapRho;
0076   iEvent.getByToken(mapRhoToken_, mapRho);
0077 
0078   edm::Handle<std::vector<double>> mapRhoM;
0079   iEvent.getByToken(mapRhoMToken_, mapRhoM);
0080 
0081   // if doCentrality is set to true calculate empty area correction
0082   // for only events with hiBin > hiBinCut
0083   int hiBin = -1;
0084   bool doEmptyArea = true;
0085   if (doCentrality_) {
0086     edm::Handle<int> cbin;
0087     iEvent.getByToken(centralityBinToken_, cbin);
0088     hiBin = *cbin;
0089 
0090     if (hiBin < hiBinCut_)
0091       doEmptyArea = false;
0092   }
0093 
0094   //Define output vectors
0095   int neta = (int)mapEtaRanges->size();
0096 
0097   auto mapToRhoCorrOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
0098   auto mapToRhoMCorrOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
0099   auto mapToRhoCorr1BinOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
0100   auto mapToRhoMCorr1BinOut = std::make_unique<std::vector<double>>(neta - 1, 1e-6);
0101 
0102   setupGrid(mapEtaRanges->at(0), mapEtaRanges->at(neta - 1));
0103 
0104   //calculate empty area correction over full acceptance leaving eta bands on the sides
0105   double allAcceptanceCorr = 1;
0106   if (doEmptyArea) {
0107     etaminJet_ = mapEtaRanges->at(0) - band_;
0108     etamaxJet_ = mapEtaRanges->at(neta - 1) + band_;
0109 
0110     calculateAreaFractionOfJets(iEvent, iSetup);
0111 
0112     allAcceptanceCorr = totalInboundArea_;
0113   }
0114 
0115   //calculate empty area correction in each eta range
0116   for (int ieta = 0; ieta < (neta - 1); ieta++) {
0117     double correctionKt = 1;
0118     double rho = mapRho->at(ieta);
0119     double rhoM = mapRhoM->at(ieta);
0120 
0121     if (doEmptyArea) {
0122       double etamin = mapEtaRanges->at(ieta);
0123       double etamax = mapEtaRanges->at(ieta + 1);
0124 
0125       etaminJet_ = etamin + band_;
0126       etamaxJet_ = etamax - band_;
0127 
0128       calculateAreaFractionOfJets(iEvent, iSetup);
0129       correctionKt = totalInboundArea_;
0130     }
0131 
0132     mapToRhoCorrOut->at(ieta) = correctionKt * rho;
0133     mapToRhoMCorrOut->at(ieta) = correctionKt * rhoM;
0134 
0135     mapToRhoCorr1BinOut->at(ieta) = allAcceptanceCorr * rho;
0136     mapToRhoMCorr1BinOut->at(ieta) = allAcceptanceCorr * rhoM;
0137   }
0138 
0139   iEvent.put(std::move(mapToRhoCorrOut), "mapToRhoCorr");
0140   iEvent.put(std::move(mapToRhoMCorrOut), "mapToRhoMCorr");
0141   iEvent.put(std::move(mapToRhoCorr1BinOut), "mapToRhoCorr1Bin");
0142   iEvent.put(std::move(mapToRhoMCorr1BinOut), "mapToRhoMCorr1Bin");
0143 
0144   //calculate rho from grid as a function of eta over full range using PF candidates
0145 
0146   auto mapRhoVsEtaGridOut = std::make_unique<std::vector<double>>(ny_, 0.);
0147   auto mapMeanRhoVsEtaGridOut = std::make_unique<std::vector<double>>(ny_, 0.);
0148   auto mapEtaMaxGridOut = std::make_unique<std::vector<double>>(ny_, 0.);
0149   auto mapEtaMinGridOut = std::make_unique<std::vector<double>>(ny_, 0.);
0150   calculateGridRho(iEvent, iSetup);
0151   if (keepGridInfo_) {
0152     for (int ieta = 0; ieta < ny_; ieta++) {
0153       mapRhoVsEtaGridOut->at(ieta) = rhoVsEta_[ieta];
0154       mapMeanRhoVsEtaGridOut->at(ieta) = meanRhoVsEta_[ieta];
0155       mapEtaMaxGridOut->at(ieta) = etaMaxGrid_[ieta];
0156       mapEtaMinGridOut->at(ieta) = etaMinGrid_[ieta];
0157     }
0158 
0159     iEvent.put(std::move(mapRhoVsEtaGridOut), "mapRhoVsEtaGrid");
0160     iEvent.put(std::move(mapMeanRhoVsEtaGridOut), "mapMeanRhoVsEtaGrid");
0161     iEvent.put(std::move(mapEtaMaxGridOut), "mapEtaMaxGrid");
0162     iEvent.put(std::move(mapEtaMinGridOut), "mapEtaMinGrid");
0163   }
0164 }
0165 
0166 //----------------------------------------------------------------------
0167 // setting a new event
0168 //----------------------------------------------------------------------
0169 // tell the background estimator that it has a new event, composed
0170 // of the specified particles.
0171 void HiFJGridEmptyAreaCalculator::calculateGridRho(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0172   vector<vector<double>> scalarPt(ny_, vector<double>(nphi_, 0.0));
0173 
0174   edm::Handle<reco::PFCandidateCollection> pfCands;
0175   iEvent.getByToken(pfCandsToken_, pfCands);
0176   const reco::PFCandidateCollection* pfCandidateColl = pfCands.product();
0177   for (unsigned icand = 0; icand < pfCandidateColl->size(); icand++) {
0178     const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
0179     //use ony the particles within the eta range
0180     if (pfCandidate.eta() < ymin_ || pfCandidate.eta() > ymax_)
0181       continue;
0182     int jeta = tileIndexEta(&pfCandidate);
0183     int jphi = tileIndexPhi(&pfCandidate);
0184     scalarPt[jeta][jphi] += pfCandidate.pt();
0185   }
0186 
0187   rhoVsEta_.resize(ny_);
0188   meanRhoVsEta_.resize(ny_);
0189   for (int jeta = 0; jeta < ny_; jeta++) {
0190     rhoVsEta_[jeta] = 0;
0191     meanRhoVsEta_[jeta] = 0;
0192     vector<double> rhoVsPhi;
0193     int nEmpty = 0;
0194 
0195     for (int jphi = 0; jphi < nphi_; jphi++) {
0196       double binpt = scalarPt[jeta][jphi];
0197       meanRhoVsEta_[jeta] += binpt;
0198       //fill in the vector for median calculation
0199       if (binpt > 0)
0200         rhoVsPhi.push_back(binpt);
0201       else
0202         nEmpty++;
0203     }
0204     meanRhoVsEta_[jeta] /= ((double)nphi_);
0205     meanRhoVsEta_[jeta] /= tileArea_;
0206 
0207     //median calculation
0208     sort(rhoVsPhi.begin(), rhoVsPhi.end());
0209     //use only the nonzero grid cells for median calculation;
0210     int nFull = nphi_ - nEmpty;
0211     if (nFull == 0) {
0212       rhoVsEta_[jeta] = 0;
0213       continue;
0214     }
0215     if (nFull % 2 == 0) {
0216       rhoVsEta_[jeta] = (rhoVsPhi[(int)(nFull / 2 - 1)] + rhoVsPhi[(int)(nFull / 2)]) / 2;
0217     } else {
0218       rhoVsEta_[jeta] = rhoVsPhi[(int)(nFull / 2)];
0219     }
0220     //correct for empty cells
0221     rhoVsEta_[jeta] *= (((double)nFull) / ((double)nphi_));
0222     //normalize to area
0223     rhoVsEta_[jeta] /= tileArea_;
0224   }
0225 }
0226 
0227 void HiFJGridEmptyAreaCalculator::calculateAreaFractionOfJets(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0228   edm::Handle<edm::View<reco::Jet>> jets;
0229   iEvent.getByToken(jetsToken_, jets);
0230 
0231   //calculate jet kt area fraction inside boundary by grid
0232   totalInboundArea_ = 0;
0233 
0234   for (auto jet = jets->begin(); jet != jets->end(); ++jet) {
0235     if (jet->eta() < etaminJet_ || jet->eta() > etamaxJet_)
0236       continue;
0237 
0238     double areaKt = jet->jetArea();
0239     setupGridJet(&*jet);
0240     std::vector<std::pair<int, int>> pfIndicesJet;
0241     std::vector<std::pair<int, int>> pfIndicesJetInbound;
0242     int nConstitJet = 0;
0243     int nConstitJetInbound = 0;
0244     for (auto daughter : jet->getJetConstituentsQuick()) {
0245       auto pfCandidate = static_cast<const reco::PFCandidate*>(daughter);
0246 
0247       int jeta = tileIndexEtaJet(&*pfCandidate);
0248       int jphi = tileIndexPhi(&*pfCandidate);
0249       pfIndicesJet.push_back(std::make_pair(jphi, jeta));
0250       nConstitJet++;
0251       if (pfCandidate->eta() < etaminJet_ && pfCandidate->eta() > etamaxJet_)
0252         continue;
0253       pfIndicesJetInbound.push_back(std::make_pair(jphi, jeta));
0254       nConstitJetInbound++;
0255     }
0256 
0257     //if the jet is well within the eta range just add the area
0258     if (nConstitJet == nConstitJetInbound) {
0259       totalInboundArea_ += areaKt;
0260       continue;
0261     }
0262 
0263     //for jets that fall outside of eta range calculate fraction of area
0264     //inside the range with a grid
0265     int nthis = 0;
0266     if (nConstitJet > 0)
0267       nthis = numJetGridCells(pfIndicesJet);
0268 
0269     int nthisInbound = 0;
0270     if (nConstitJetInbound > 0)
0271       nthisInbound = numJetGridCells(pfIndicesJetInbound);
0272 
0273     double fractionArea = ((double)nthisInbound) / ((double)nthis);
0274     totalInboundArea_ += areaKt * fractionArea;
0275   }
0276 
0277   //divide by the total area in that range
0278   totalInboundArea_ /= ((etamaxJet_ - etaminJet_) * twopi_);
0279 
0280   //the fraction can still be greater than 1 because kt area fraction inside
0281   //the range can differ from what we calculated with the grid
0282   if (totalInboundArea_ > 1)
0283     totalInboundArea_ = 1;
0284 }
0285 
0286 // #ifndef FASTJET_GMBGE_USEFJGRID
0287 //----------------------------------------------------------------------
0288 // protected material
0289 //----------------------------------------------------------------------
0290 // configure the grid
0291 void HiFJGridEmptyAreaCalculator::setupGrid(double etamin, double etamax) {
0292   // since we've exchanged the arguments of the grid constructor,
0293   // there's a danger of calls with exchanged ymax,spacing arguments --
0294   // the following check should catch most such situations.
0295   ymin_ = etamin;
0296   ymax_ = etamax;
0297 
0298   assert(ymax_ - ymin_ >= gridWidth_);
0299 
0300   // this grid-definition code is becoming repetitive -- it should
0301   // probably be moved somewhere central...
0302   double nyDouble = (ymax_ - ymin_) / gridWidth_;
0303   ny_ = int(nyDouble + 0.5);
0304   dy_ = (ymax_ - ymin_) / ny_;
0305 
0306   nphi_ = int(twopi_ / gridWidth_ + 0.5);
0307   dphi_ = twopi_ / nphi_;
0308 
0309   // some sanity checking (could throw a fastjet::Error)
0310   assert(ny_ >= 1 && nphi_ >= 1);
0311 
0312   ntotal_ = nphi_ * ny_;
0313   //_scalar_pt.resize(_ntotal);
0314   tileArea_ = dy_ * dphi_;
0315 
0316   etaMaxGrid_.resize(ny_);
0317   etaMinGrid_.resize(ny_);
0318   for (int jeta = 0; jeta < ny_; jeta++) {
0319     etaMinGrid_[jeta] = etamin + dy_ * ((double)jeta);
0320     etaMaxGrid_[jeta] = etamin + dy_ * ((double)jeta + 1.);
0321   }
0322 }
0323 
0324 //----------------------------------------------------------------------
0325 // retrieve the grid tile index for a given PseudoJet
0326 int HiFJGridEmptyAreaCalculator::tileIndexPhi(const reco::PFCandidate* pfCand) {
0327   // directly taking int does not work for values between -1 and 0
0328   // so use floor instead
0329   // double iy_double = (p.rap() - _ymin) / _dy;
0330   // if (iy_double < 0.0) return -1;
0331   // int iy = int(iy_double);
0332   // if (iy >= _ny) return -1;
0333 
0334   // writing it as below gives a huge speed gain (factor two!). Even
0335   // though answers are identical and the routine here is not the
0336   // speed-critical step. It's not at all clear why.
0337 
0338   int iphi = int((pfCand->phi() + (twopi_ / 2.)) / dphi_);
0339   assert(iphi >= 0 && iphi <= nphi_);
0340   if (iphi == nphi_)
0341     iphi = 0;  // just in case of rounding errors
0342 
0343   return iphi;
0344 }
0345 
0346 //----------------------------------------------------------------------
0347 // retrieve the grid tile index for a given PseudoJet
0348 int HiFJGridEmptyAreaCalculator::tileIndexEta(const reco::PFCandidate* pfCand) {
0349   // directly taking int does not work for values between -1 and 0
0350   // so use floor instead
0351   // double iy_double = (p.rap() - _ymin) / _dy;
0352   // if (iy_double < 0.0) return -1;
0353   // int iy = int(iy_double);
0354   // if (iy >= _ny) return -1;
0355 
0356   // writing it as below gives a huge speed gain (factor two!). Even
0357   // though answers are identical and the routine here is not the
0358   // speed-critical step. It's not at all clear why.
0359   int iy = int(floor((pfCand->eta() - ymin_) / dy_));
0360   if (iy < 0 || iy >= ny_)
0361     return -1;
0362 
0363   assert(iy < ny_ && iy >= 0);
0364 
0365   return iy;
0366 }
0367 
0368 // #ifndef FASTJET_GMBGE_USEFJGRID
0369 //----------------------------------------------------------------------
0370 // protected material
0371 //----------------------------------------------------------------------
0372 // configure the grid
0373 void HiFJGridEmptyAreaCalculator::setupGridJet(const reco::Jet* jet) {
0374   // since we've exchanged the arguments of the grid constructor,
0375   // there's a danger of calls with exchanged ymax,spacing arguments --
0376   // the following check should catch most such situations.
0377   yminJet_ = jet->eta() - 0.6;
0378   ymaxJet_ = jet->eta() + 0.6;
0379 
0380   assert(ymaxJet_ - yminJet_ >= gridWidth_);
0381 
0382   // this grid-definition code is becoming repetitive -- it should
0383   // probably be moved somewhere central...
0384   double nyDouble = (ymaxJet_ - yminJet_) / gridWidth_;
0385   nyJet_ = int(nyDouble + 0.5);
0386   dyJet_ = (ymaxJet_ - yminJet_) / nyJet_;
0387 
0388   assert(nyJet_ >= 1);
0389 
0390   ntotalJet_ = nphi_ * nyJet_;
0391   //_scalar_pt.resize(_ntotal);
0392 }
0393 
0394 //----------------------------------------------------------------------
0395 // retrieve the grid tile index for a given PseudoJet
0396 int HiFJGridEmptyAreaCalculator::tileIndexEtaJet(const reco::PFCandidate* pfCand) {
0397   // directly taking int does not work for values between -1 and 0
0398   // so use floor instead
0399   // double iy_double = (p.rap() - _ymin) / _dy;
0400   // if (iy_double < 0.0) return -1;
0401   // int iy = int(iy_double);
0402   // if (iy >= _ny) return -1;
0403 
0404   // writing it as below gives a huge speed gain (factor two!). Even
0405   // though answers are identical and the routine here is not the
0406   // speed-critical step. It's not at all clear why.
0407   int iyjet = int(floor((pfCand->eta() - yminJet_) / dy_));
0408   if (iyjet < 0 || iyjet >= nyJet_)
0409     return -1;
0410 
0411   assert(iyjet < nyJet_ && iyjet >= 0);
0412 
0413   return iyjet;
0414 }
0415 
0416 int HiFJGridEmptyAreaCalculator::numJetGridCells(std::vector<std::pair<int, int>>& indices) {
0417   int ngrid = 0;
0418   //sort phi eta grid indices in phi
0419   std::sort(indices.begin(), indices.end());
0420   int lowestJPhi = indices[0].first;
0421   int lowestJEta = indices[0].second;
0422   int highestJEta = lowestJEta;
0423 
0424   //for each fixed phi value calculate the number of grids in eta
0425   for (unsigned int iconst = 1; iconst < indices.size(); iconst++) {
0426     int jphi = indices[iconst].first;
0427     int jeta = indices[iconst].second;
0428     if (jphi == lowestJPhi) {
0429       if (jeta < lowestJEta)
0430         lowestJEta = jeta;
0431       if (jeta > highestJEta)
0432         highestJEta = jeta;
0433     } else {
0434       lowestJPhi = jphi;
0435       ngrid += highestJEta - lowestJEta + 1;
0436       lowestJEta = jeta;
0437       highestJEta = jeta;
0438     }
0439   }
0440   ngrid += highestJEta - lowestJEta + 1;
0441   return ngrid;
0442 }
0443 
0444 void HiFJGridEmptyAreaCalculator::beginStream(edm::StreamID) {}
0445 
0446 void HiFJGridEmptyAreaCalculator::endStream() {}
0447 
0448 void HiFJGridEmptyAreaCalculator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0449   edm::ParameterSetDescription desc;
0450   desc.add<edm::InputTag>("jetSource", edm::InputTag("kt4PFJets"));
0451   desc.add<edm::InputTag>("CentralityBinSrc", edm::InputTag("centralityBin"));
0452   desc.add<edm::InputTag>("mapEtaEdges", edm::InputTag("mapEtaEdges"));
0453   desc.add<edm::InputTag>("mapToRho", edm::InputTag("mapToRho"));
0454   desc.add<edm::InputTag>("mapToRhoM", edm::InputTag("mapToRhoM"));
0455   desc.add<edm::InputTag>("pfCandSource", edm::InputTag("particleFlow"));
0456   desc.add<double>("gridWidth", 0.05);
0457   desc.add<double>("bandWidth", 0.2);
0458   desc.add<bool>("doCentrality", true);
0459   desc.add<int>("hiBinCut", 100);
0460   desc.add<bool>("keepGridInfo", false);
0461   descriptions.add("hiFJGridEmptyAreaCalculator", desc);
0462 }
0463 
0464 DEFINE_FWK_MODULE(HiFJGridEmptyAreaCalculator);