Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:18 

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::CandidateView>(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::CandidateView> pfCands;
0175   iEvent.getByToken(pfCandsToken_, pfCands);
0176   const auto& pfCandidateColl = pfCands.product();
0177   for (const auto& pfCandidate : *pfCandidateColl) {
0178     //use ony the particles within the eta range
0179     if (pfCandidate.eta() < ymin_ || pfCandidate.eta() > ymax_)
0180       continue;
0181     int jeta = tileIndexEta(&pfCandidate);
0182     int jphi = tileIndexPhi(&pfCandidate);
0183     scalarPt[jeta][jphi] += pfCandidate.pt();
0184   }
0185 
0186   rhoVsEta_.resize(ny_);
0187   meanRhoVsEta_.resize(ny_);
0188   for (int jeta = 0; jeta < ny_; jeta++) {
0189     rhoVsEta_[jeta] = 0;
0190     meanRhoVsEta_[jeta] = 0;
0191     vector<double> rhoVsPhi;
0192     int nEmpty = 0;
0193 
0194     for (int jphi = 0; jphi < nphi_; jphi++) {
0195       double binpt = scalarPt[jeta][jphi];
0196       meanRhoVsEta_[jeta] += binpt;
0197       //fill in the vector for median calculation
0198       if (binpt > 0)
0199         rhoVsPhi.push_back(binpt);
0200       else
0201         nEmpty++;
0202     }
0203     meanRhoVsEta_[jeta] /= ((double)nphi_);
0204     meanRhoVsEta_[jeta] /= tileArea_;
0205 
0206     //median calculation
0207     sort(rhoVsPhi.begin(), rhoVsPhi.end());
0208     //use only the nonzero grid cells for median calculation;
0209     int nFull = nphi_ - nEmpty;
0210     if (nFull == 0) {
0211       rhoVsEta_[jeta] = 0;
0212       continue;
0213     }
0214     if (nFull % 2 == 0) {
0215       rhoVsEta_[jeta] = (rhoVsPhi[(int)(nFull / 2 - 1)] + rhoVsPhi[(int)(nFull / 2)]) / 2;
0216     } else {
0217       rhoVsEta_[jeta] = rhoVsPhi[(int)(nFull / 2)];
0218     }
0219     //correct for empty cells
0220     rhoVsEta_[jeta] *= (((double)nFull) / ((double)nphi_));
0221     //normalize to area
0222     rhoVsEta_[jeta] /= tileArea_;
0223   }
0224 }
0225 
0226 void HiFJGridEmptyAreaCalculator::calculateAreaFractionOfJets(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0227   edm::Handle<edm::View<reco::Jet>> jets;
0228   iEvent.getByToken(jetsToken_, jets);
0229 
0230   //calculate jet kt area fraction inside boundary by grid
0231   totalInboundArea_ = 0;
0232 
0233   for (const auto& jet : *jets) {
0234     if (jet.eta() < etaminJet_ || jet.eta() > etamaxJet_)
0235       continue;
0236 
0237     double areaKt = jet.jetArea();
0238     setupGridJet(&jet);
0239     std::vector<std::pair<int, int>> pfIndicesJet;
0240     std::vector<std::pair<int, int>> pfIndicesJetInbound;
0241     int nConstitJet = 0;
0242     int nConstitJetInbound = 0;
0243     for (const auto& daughter : jet.getJetConstituentsQuick()) {
0244       auto pfCandidate = static_cast<const reco::PFCandidate*>(daughter);
0245 
0246       int jeta = tileIndexEtaJet(&*pfCandidate);
0247       int jphi = tileIndexPhi(&*pfCandidate);
0248       pfIndicesJet.push_back(std::make_pair(jphi, jeta));
0249       nConstitJet++;
0250       if (pfCandidate->eta() < etaminJet_ && pfCandidate->eta() > etamaxJet_)
0251         continue;
0252       pfIndicesJetInbound.push_back(std::make_pair(jphi, jeta));
0253       nConstitJetInbound++;
0254     }
0255 
0256     //if the jet is well within the eta range just add the area
0257     if (nConstitJet == nConstitJetInbound) {
0258       totalInboundArea_ += areaKt;
0259       continue;
0260     }
0261 
0262     //for jets that fall outside of eta range calculate fraction of area
0263     //inside the range with a grid
0264     int nthis = 0;
0265     if (nConstitJet > 0)
0266       nthis = numJetGridCells(pfIndicesJet);
0267 
0268     int nthisInbound = 0;
0269     if (nConstitJetInbound > 0)
0270       nthisInbound = numJetGridCells(pfIndicesJetInbound);
0271 
0272     double fractionArea = ((double)nthisInbound) / ((double)nthis);
0273     totalInboundArea_ += areaKt * fractionArea;
0274   }
0275 
0276   //divide by the total area in that range
0277   totalInboundArea_ /= ((etamaxJet_ - etaminJet_) * twopi_);
0278 
0279   //the fraction can still be greater than 1 because kt area fraction inside
0280   //the range can differ from what we calculated with the grid
0281   if (totalInboundArea_ > 1)
0282     totalInboundArea_ = 1;
0283 }
0284 
0285 // #ifndef FASTJET_GMBGE_USEFJGRID
0286 //----------------------------------------------------------------------
0287 // protected material
0288 //----------------------------------------------------------------------
0289 // configure the grid
0290 void HiFJGridEmptyAreaCalculator::setupGrid(double etamin, double etamax) {
0291   // since we've exchanged the arguments of the grid constructor,
0292   // there's a danger of calls with exchanged ymax,spacing arguments --
0293   // the following check should catch most such situations.
0294   ymin_ = etamin;
0295   ymax_ = etamax;
0296 
0297   assert(ymax_ - ymin_ >= gridWidth_);
0298 
0299   // this grid-definition code is becoming repetitive -- it should
0300   // probably be moved somewhere central...
0301   double nyDouble = (ymax_ - ymin_) / gridWidth_;
0302   ny_ = int(nyDouble + 0.5);
0303   dy_ = (ymax_ - ymin_) / ny_;
0304 
0305   nphi_ = int(twopi_ / gridWidth_ + 0.5);
0306   dphi_ = twopi_ / nphi_;
0307 
0308   // some sanity checking (could throw a fastjet::Error)
0309   assert(ny_ >= 1 && nphi_ >= 1);
0310 
0311   ntotal_ = nphi_ * ny_;
0312   //_scalar_pt.resize(_ntotal);
0313   tileArea_ = dy_ * dphi_;
0314 
0315   etaMaxGrid_.resize(ny_);
0316   etaMinGrid_.resize(ny_);
0317   for (int jeta = 0; jeta < ny_; jeta++) {
0318     etaMinGrid_[jeta] = etamin + dy_ * ((double)jeta);
0319     etaMaxGrid_[jeta] = etamin + dy_ * ((double)jeta + 1.);
0320   }
0321 }
0322 
0323 //----------------------------------------------------------------------
0324 // retrieve the grid tile index for a given PseudoJet
0325 int HiFJGridEmptyAreaCalculator::tileIndexPhi(const reco::Candidate* pfCand) {
0326   // directly taking int does not work for values between -1 and 0
0327   // so use floor instead
0328   // double iy_double = (p.rap() - _ymin) / _dy;
0329   // if (iy_double < 0.0) return -1;
0330   // int iy = int(iy_double);
0331   // if (iy >= _ny) return -1;
0332 
0333   // writing it as below gives a huge speed gain (factor two!). Even
0334   // though answers are identical and the routine here is not the
0335   // speed-critical step. It's not at all clear why.
0336 
0337   int iphi = int((pfCand->phi() + (twopi_ / 2.)) / dphi_);
0338   assert(iphi >= 0 && iphi <= nphi_);
0339   if (iphi == nphi_)
0340     iphi = 0;  // just in case of rounding errors
0341 
0342   return iphi;
0343 }
0344 
0345 //----------------------------------------------------------------------
0346 // retrieve the grid tile index for a given PseudoJet
0347 int HiFJGridEmptyAreaCalculator::tileIndexEta(const reco::Candidate* pfCand) {
0348   // directly taking int does not work for values between -1 and 0
0349   // so use floor instead
0350   // double iy_double = (p.rap() - _ymin) / _dy;
0351   // if (iy_double < 0.0) return -1;
0352   // int iy = int(iy_double);
0353   // if (iy >= _ny) return -1;
0354 
0355   // writing it as below gives a huge speed gain (factor two!). Even
0356   // though answers are identical and the routine here is not the
0357   // speed-critical step. It's not at all clear why.
0358   int iy = int(floor((pfCand->eta() - ymin_) / dy_));
0359   if (iy < 0 || iy >= ny_)
0360     return -1;
0361 
0362   assert(iy < ny_ && iy >= 0);
0363 
0364   return iy;
0365 }
0366 
0367 // #ifndef FASTJET_GMBGE_USEFJGRID
0368 //----------------------------------------------------------------------
0369 // protected material
0370 //----------------------------------------------------------------------
0371 // configure the grid
0372 void HiFJGridEmptyAreaCalculator::setupGridJet(const reco::Jet* jet) {
0373   // since we've exchanged the arguments of the grid constructor,
0374   // there's a danger of calls with exchanged ymax,spacing arguments --
0375   // the following check should catch most such situations.
0376   yminJet_ = jet->eta() - 0.6;
0377   ymaxJet_ = jet->eta() + 0.6;
0378 
0379   assert(ymaxJet_ - yminJet_ >= gridWidth_);
0380 
0381   // this grid-definition code is becoming repetitive -- it should
0382   // probably be moved somewhere central...
0383   double nyDouble = (ymaxJet_ - yminJet_) / gridWidth_;
0384   nyJet_ = int(nyDouble + 0.5);
0385   dyJet_ = (ymaxJet_ - yminJet_) / nyJet_;
0386 
0387   assert(nyJet_ >= 1);
0388 
0389   ntotalJet_ = nphi_ * nyJet_;
0390   //_scalar_pt.resize(_ntotal);
0391 }
0392 
0393 //----------------------------------------------------------------------
0394 // retrieve the grid tile index for a given PseudoJet
0395 int HiFJGridEmptyAreaCalculator::tileIndexEtaJet(const reco::Candidate* pfCand) {
0396   // directly taking int does not work for values between -1 and 0
0397   // so use floor instead
0398   // double iy_double = (p.rap() - _ymin) / _dy;
0399   // if (iy_double < 0.0) return -1;
0400   // int iy = int(iy_double);
0401   // if (iy >= _ny) return -1;
0402 
0403   // writing it as below gives a huge speed gain (factor two!). Even
0404   // though answers are identical and the routine here is not the
0405   // speed-critical step. It's not at all clear why.
0406   int iyjet = int(floor((pfCand->eta() - yminJet_) / dy_));
0407   if (iyjet < 0 || iyjet >= nyJet_)
0408     return -1;
0409 
0410   assert(iyjet < nyJet_ && iyjet >= 0);
0411 
0412   return iyjet;
0413 }
0414 
0415 int HiFJGridEmptyAreaCalculator::numJetGridCells(std::vector<std::pair<int, int>>& indices) {
0416   int ngrid = 0;
0417   //sort phi eta grid indices in phi
0418   std::sort(indices.begin(), indices.end());
0419   int lowestJPhi = indices[0].first;
0420   int lowestJEta = indices[0].second;
0421   int highestJEta = lowestJEta;
0422 
0423   //for each fixed phi value calculate the number of grids in eta
0424   for (const auto& index : indices) {
0425     int jphi = index.first;
0426     int jeta = index.second;
0427     if (jphi == lowestJPhi) {
0428       if (jeta < lowestJEta)
0429         lowestJEta = jeta;
0430       if (jeta > highestJEta)
0431         highestJEta = jeta;
0432     } else {
0433       lowestJPhi = jphi;
0434       ngrid += highestJEta - lowestJEta + 1;
0435       lowestJEta = jeta;
0436       highestJEta = jeta;
0437     }
0438   }
0439   ngrid += highestJEta - lowestJEta + 1;
0440   return ngrid;
0441 }
0442 
0443 void HiFJGridEmptyAreaCalculator::beginStream(edm::StreamID) {}
0444 
0445 void HiFJGridEmptyAreaCalculator::endStream() {}
0446 
0447 void HiFJGridEmptyAreaCalculator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0448   edm::ParameterSetDescription desc;
0449   desc.add<edm::InputTag>("jetSource", edm::InputTag("kt4PFJets"));
0450   desc.add<edm::InputTag>("CentralityBinSrc", edm::InputTag("centralityBin"));
0451   desc.add<edm::InputTag>("mapEtaEdges", edm::InputTag("mapEtaEdges"));
0452   desc.add<edm::InputTag>("mapToRho", edm::InputTag("mapToRho"));
0453   desc.add<edm::InputTag>("mapToRhoM", edm::InputTag("mapToRhoM"));
0454   desc.add<edm::InputTag>("pfCandSource", edm::InputTag("particleFlow"));
0455   desc.add<double>("gridWidth", 0.05);
0456   desc.add<double>("bandWidth", 0.2);
0457   desc.add<bool>("doCentrality", true);
0458   desc.add<int>("hiBinCut", 100);
0459   desc.add<bool>("keepGridInfo", false);
0460   descriptions.add("hiFJGridEmptyAreaCalculator", desc);
0461 }
0462 
0463 DEFINE_FWK_MODULE(HiFJGridEmptyAreaCalculator);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!