HiFJGridEmptyAreaCalculator.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    HiJetBackground/HiFJGridEmptyAreaCalculator
// Class:      HiFJGridEmptyAreaCalculator
// Based on:   fastjet/tools/GridMedianBackgroundEstimator 
//
// Original Author:  Doga Gulhan
//         Created:  Wed Mar 16 14:00:04 CET 2016
//
//

#include "RecoHI/HiJetAlgos/plugins/HiFJGridEmptyAreaCalculator.h"
using namespace std;
using namespace edm;


HiFJGridEmptyAreaCalculator::HiFJGridEmptyAreaCalculator(const edm::ParameterSet& iConfig):
  gridWidth_(iConfig.getParameter<double>("gridWidth")),
  band_(iConfig.getParameter<double>("bandWidth")),
  hiBinCut_(iConfig.getParameter<int>("hiBinCut")),
  doCentrality_(iConfig.getParameter<bool>("doCentrality")),
  keepGridInfo_(iConfig.getParameter<bool>("keepGridInfo"))
{
  ymin_ = -99;
  ymax_ = -99;
  dy_ = -99;
  dphi_ = -99;
  tileArea_ = - 99;
  
  dyJet_ = 99;
  yminJet_ = - 99;
  ymaxJet_ = -99;
  totalInboundArea_ = -99;
  etaminJet_= -99;
  etamaxJet_ = -99;
  
  ny_ = 0;
  nphi_ = 0;
  ntotal_ = 0;
 
  ntotalJet_ = 0;
  nyJet_ = 0;
    
  pfCandsToken_ = consumes<reco::PFCandidateCollection>(iConfig.getParameter<edm::InputTag>( "pfCandSource" ));
  mapEtaToken_ = consumes<std::vector<double> >(iConfig.getParameter<edm::InputTag>( "mapEtaEdges" ));
  mapRhoToken_ = consumes<std::vector<double> >(iConfig.getParameter<edm::InputTag>( "mapToRho" ));
  mapRhoMToken_ = consumes<std::vector<double> >(iConfig.getParameter<edm::InputTag>( "mapToRhoM" ));
  jetsToken_ = consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>( "jetSource" ));
  centralityBinToken_ = consumes<int>(iConfig.getParameter<edm::InputTag>("CentralityBinSrc"));

  //register your products
  produces<std::vector<double > >("mapEmptyCorrFac"); 
  produces<std::vector<double > >("mapToRhoCorr");
  produces<std::vector<double > >("mapToRhoMCorr");
  produces<std::vector<double > >("mapToRhoCorr1Bin");
  produces<std::vector<double > >("mapToRhoMCorr1Bin");
  //rho calculation on a grid using median
  if(keepGridInfo_){
    produces<std::vector<double > >("mapRhoVsEtaGrid");
    produces<std::vector<double > >("mapMeanRhoVsEtaGrid");
    produces<std::vector<double > >("mapEtaMaxGrid");
    produces<std::vector<double > >("mapEtaMinGrid");
  }
}


HiFJGridEmptyAreaCalculator::~HiFJGridEmptyAreaCalculator()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}

void
HiFJGridEmptyAreaCalculator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  using namespace edm;
  //Values from rho calculator  
  edm::Handle<std::vector<double>> mapEtaRanges;
  iEvent.getByToken(mapEtaToken_, mapEtaRanges);
  
  edm::Handle<std::vector<double>> mapRho;
  iEvent.getByToken(mapRhoToken_, mapRho);

  edm::Handle<std::vector<double>> mapRhoM;  
  iEvent.getByToken(mapRhoMToken_, mapRhoM);

  // if doCentrality is set to true calculate empty area correction
  // for only events with hiBin > hiBinCut  
  int hiBin = -1;
  bool doEmptyArea = true; 
  if(doCentrality_){
    edm::Handle<int> cbin;
    iEvent.getByToken(centralityBinToken_,cbin);
    hiBin = *cbin;
   
    if(hiBin < hiBinCut_) doEmptyArea = false;
  }
  
  //Define output vectors
  int neta = (int)mapEtaRanges->size();
   
  auto mapToRhoCorrOut = std::make_unique<std::vector<double>>(neta-1,1e-6);
  auto mapToRhoMCorrOut = std::make_unique<std::vector<double>>(neta-1,1e-6);
  auto mapToRhoCorr1BinOut = std::make_unique<std::vector<double>>(neta-1,1e-6);
  auto mapToRhoMCorr1BinOut = std::make_unique<std::vector<double>>(neta-1,1e-6);

  setupGrid(mapEtaRanges->at(0), mapEtaRanges->at(neta-1));
  
  //calculate empty area correction over full acceptance leaving eta bands on the sides
  double allAcceptanceCorr = 1;
  if(doEmptyArea){
    etaminJet_ = mapEtaRanges->at(0) - band_;
    etamaxJet_ = mapEtaRanges->at(neta-1) + band_;
  
    calculateAreaFractionOfJets(iEvent, iSetup);
  
    allAcceptanceCorr =  totalInboundArea_;
  }
  
  //calculate empty area correction in each eta range
  for(int ieta = 0; ieta<(neta-1); ieta++) {
   
    double correctionKt = 1;   
    double rho = mapRho->at(ieta);
    double rhoM = mapRhoM->at(ieta);
    
    if(doEmptyArea){
      double etamin = mapEtaRanges->at(ieta);
      double etamax = mapEtaRanges->at(ieta+1);
      
      etaminJet_ = etamin + band_;  
      etamaxJet_ = etamax - band_;  
   
      calculateAreaFractionOfJets(iEvent, iSetup);
      correctionKt = totalInboundArea_;
    }
   
    mapToRhoCorrOut->at(ieta) = correctionKt*rho;
    mapToRhoMCorrOut->at(ieta) = correctionKt*rhoM;
   
    mapToRhoCorr1BinOut->at(ieta) = allAcceptanceCorr*rho;
    mapToRhoMCorr1BinOut->at(ieta) = allAcceptanceCorr*rhoM;
  }

  iEvent.put(std::move(mapToRhoCorrOut),"mapToRhoCorr");
  iEvent.put(std::move(mapToRhoMCorrOut),"mapToRhoMCorr");
  iEvent.put(std::move(mapToRhoCorr1BinOut),"mapToRhoCorr1Bin");
  iEvent.put(std::move(mapToRhoMCorr1BinOut),"mapToRhoMCorr1Bin");
  
  //calculate rho from grid as a function of eta over full range using PF candidates
  
  auto mapRhoVsEtaGridOut = std::make_unique<std::vector<double>>(ny_,0.);
  auto mapMeanRhoVsEtaGridOut = std::make_unique<std::vector<double>>(ny_,0.);
  auto mapEtaMaxGridOut = std::make_unique<std::vector<double>>(ny_,0.);
  auto mapEtaMinGridOut = std::make_unique<std::vector<double>>(ny_,0.);
  calculateGridRho(iEvent, iSetup);
  if(keepGridInfo_){
    for(int ieta = 0; ieta < ny_; ieta++) {
      mapRhoVsEtaGridOut->at(ieta) = rhoVsEta_[ieta];
      mapMeanRhoVsEtaGridOut->at(ieta) = meanRhoVsEta_[ieta];
      mapEtaMaxGridOut->at(ieta) = etaMaxGrid_[ieta];
      mapEtaMinGridOut->at(ieta) = etaMinGrid_[ieta];
    }

    iEvent.put(std::move(mapRhoVsEtaGridOut),"mapRhoVsEtaGrid");
    iEvent.put(std::move(mapMeanRhoVsEtaGridOut),"mapMeanRhoVsEtaGrid");
    iEvent.put(std::move(mapEtaMaxGridOut),"mapEtaMaxGrid");
    iEvent.put(std::move(mapEtaMinGridOut),"mapEtaMinGrid");
  }
}

//----------------------------------------------------------------------
// setting a new event
//----------------------------------------------------------------------
// tell the background estimator that it has a new event, composed
// of the specified particles.
void
HiFJGridEmptyAreaCalculator::calculateGridRho(const edm::Event& iEvent, const edm::EventSetup& iSetup) {

  vector<vector<double>> scalarPt(ny_, vector<double>(nphi_, 0.0));
  
  edm::Handle<reco::PFCandidateCollection> pfCands;
  iEvent.getByToken(pfCandsToken_, pfCands);
  const reco::PFCandidateCollection *pfCandidateColl = pfCands.product();
  for(unsigned icand = 0; icand < pfCandidateColl->size(); icand++) {
   const reco::PFCandidate pfCandidate = pfCandidateColl->at(icand);
   //use ony the particles within the eta range
   if (pfCandidate.eta() < ymin_ || pfCandidate.eta() > ymax_ ) continue;
   int jeta = tileIndexEta(&pfCandidate);
   int jphi = tileIndexPhi(&pfCandidate);
   scalarPt[jeta][jphi] += pfCandidate.pt();
  }
  
  rhoVsEta_.resize(ny_);
  meanRhoVsEta_.resize(ny_);
  for(int jeta = 0; jeta < ny_; jeta++){
  
    rhoVsEta_[jeta] = 0;
    meanRhoVsEta_[jeta] = 0;
    vector<double> rhoVsPhi;
    int nEmpty = 0;
    
    for(int jphi = 0; jphi < nphi_; jphi++){
      double binpt = scalarPt[jeta][jphi];
      meanRhoVsEta_[jeta] += binpt;
      //fill in the vector for median calculation
      if(binpt > 0) rhoVsPhi.push_back(binpt);
      else nEmpty++;
    }
    meanRhoVsEta_[jeta] /= ((double)nphi_);
    meanRhoVsEta_[jeta] /= tileArea_;

    //median calculation
    sort(rhoVsPhi.begin(), rhoVsPhi.end());
    //use only the nonzero grid cells for median calculation;
    int nFull = nphi_ - nEmpty;
    if(nFull == 0){
      rhoVsEta_[jeta] = 0;
      continue;
    }
    if (nFull  % 2 == 0)
      {
        rhoVsEta_[jeta] = (rhoVsPhi[(int)(nFull / 2 - 1)] + rhoVsPhi[(int)(nFull / 2)]) / 2;
      }
    else 
      {
        rhoVsEta_[jeta] = rhoVsPhi[(int)(nFull / 2)];
      }
    //correct for empty cells 
    rhoVsEta_[jeta] *= (((double) nFull)/((double) nphi_));
    //normalize to area
    rhoVsEta_[jeta] /= tileArea_;
  }
}

void
HiFJGridEmptyAreaCalculator::calculateAreaFractionOfJets(const edm::Event& iEvent, const edm::EventSetup& iSetup){
  edm::Handle<edm::View<reco::Jet> > jets;
  iEvent.getByToken(jetsToken_, jets);
  
  //calculate jet kt area fraction inside boundary by grid
  totalInboundArea_ = 0;
  
  for(auto jet = jets->begin(); jet != jets->end(); ++jet) {
    if (jet->eta() < etaminJet_ || jet->eta() > etamaxJet_) continue;
   
    double areaKt = jet->jetArea(); 
    setupGridJet(&*jet);
    std::vector<std::pair<int, int> > pfIndicesJet;
    std::vector<std::pair<int, int> > pfIndicesJetInbound;
    int nConstitJet = 0;
    int nConstitJetInbound = 0;
    for(auto daughter : jet->getJetConstituentsQuick()){
      auto pfCandidate = static_cast<const reco::PFCandidate*>(daughter);

      int jeta = tileIndexEtaJet(&*pfCandidate);
      int jphi = tileIndexPhi(&*pfCandidate);
      pfIndicesJet.push_back(std::make_pair(jphi, jeta));
      nConstitJet++;
      if (pfCandidate->eta() < etaminJet_ && pfCandidate->eta() > etamaxJet_) continue;
      pfIndicesJetInbound.push_back(std::make_pair(jphi, jeta));
      nConstitJetInbound++;
   }   
   
   //if the jet is well within the eta range just add the area
   if(nConstitJet == nConstitJetInbound){
    totalInboundArea_ += areaKt;
    continue;
   }
   
   //for jets that fall outside of eta range calculate fraction of area
   //inside the range with a grid
   int nthis = 0;
   if (nConstitJet > 0) nthis = numJetGridCells(pfIndicesJet);

   int nthisInbound = 0;
   if (nConstitJetInbound > 0) nthisInbound = numJetGridCells(pfIndicesJetInbound);
   
  
   double fractionArea = ((double)nthisInbound)/((double)nthis);
   totalInboundArea_ += areaKt*fractionArea;    
  } 
  
  //divide by the total area in that range
  totalInboundArea_ /= ((etamaxJet_ - etaminJet_)*twopi_);
  
  //the fraction can still be greater than 1 because kt area fraction inside 
  //the range can differ from what we calculated with the grid
  if (totalInboundArea_ > 1) totalInboundArea_ = 1;
}


// #ifndef FASTJET_GMBGE_USEFJGRID
//----------------------------------------------------------------------
// protected material
//----------------------------------------------------------------------
// configure the grid
void
HiFJGridEmptyAreaCalculator::setupGrid(double etamin, double etamax) {

  // since we've exchanged the arguments of the grid constructor,
  // there's a danger of calls with exchanged ymax,spacing arguments -- 
  // the following check should catch most such situations.
  ymin_ = etamin;
  ymax_ = etamax;
  
  assert(ymax_ - ymin_ >= gridWidth_);

  // this grid-definition code is becoming repetitive -- it should
  // probably be moved somewhere central...
  double nyDouble = (ymax_ - ymin_) / gridWidth_;
  ny_ = int(nyDouble+0.5);
  dy_ = (ymax_ - ymin_) / ny_;
  
  nphi_ = int (twopi_ / gridWidth_ + 0.5);
  dphi_ = twopi_ / nphi_;

  // some sanity checking (could throw a fastjet::Error)
  assert(ny_ >= 1 && nphi_ >= 1);

  ntotal_ = nphi_ * ny_;
  //_scalar_pt.resize(_ntotal);
  tileArea_ = dy_ * dphi_;
  
  
  etaMaxGrid_.resize(ny_);
  etaMinGrid_.resize(ny_);
  for(int jeta = 0; jeta < ny_; jeta++){
   etaMinGrid_[jeta] = etamin + dy_*((double)jeta);
   etaMaxGrid_[jeta] = etamin + dy_*((double)jeta + 1.);
  }
}

//----------------------------------------------------------------------
// retrieve the grid tile index for a given PseudoJet
int
HiFJGridEmptyAreaCalculator::tileIndexPhi(const reco::PFCandidate *pfCand)  {
  // directly taking int does not work for values between -1 and 0
  // so use floor instead
  // double iy_double = (p.rap() - _ymin) / _dy;
  // if (iy_double < 0.0) return -1;
  // int iy = int(iy_double);
  // if (iy >= _ny) return -1;

  // writing it as below gives a huge speed gain (factor two!). Even
  // though answers are identical and the routine here is not the
  // speed-critical step. It's not at all clear why.
 
  int iphi = int( (pfCand->phi() + (twopi_/2.))/dphi_ );
  assert(iphi >= 0 && iphi <= nphi_);
  if (iphi == nphi_) iphi = 0; // just in case of rounding errors

  return iphi;
}

//----------------------------------------------------------------------
// retrieve the grid tile index for a given PseudoJet
int
HiFJGridEmptyAreaCalculator::tileIndexEta(const reco::PFCandidate *pfCand)  {
  // directly taking int does not work for values between -1 and 0
  // so use floor instead
  // double iy_double = (p.rap() - _ymin) / _dy;
  // if (iy_double < 0.0) return -1;
  // int iy = int(iy_double);
  // if (iy >= _ny) return -1;

  // writing it as below gives a huge speed gain (factor two!). Even
  // though answers are identical and the routine here is not the
  // speed-critical step. It's not at all clear why.
  int iy = int(floor( (pfCand->eta() - ymin_) / dy_ ));
  if (iy < 0 || iy >= ny_) return -1;
  
  assert (iy < ny_ && iy >= 0);

  return iy;
}

// #ifndef FASTJET_GMBGE_USEFJGRID
//----------------------------------------------------------------------
// protected material
//----------------------------------------------------------------------
// configure the grid
void
HiFJGridEmptyAreaCalculator::setupGridJet(const reco::Jet *jet) {

  // since we've exchanged the arguments of the grid constructor,
  // there's a danger of calls with exchanged ymax,spacing arguments -- 
  // the following check should catch most such situations.
  yminJet_ = jet->eta()-0.6;
  ymaxJet_ = jet->eta()+0.6;
  
  assert(ymaxJet_ - yminJet_ >= gridWidth_);

  // this grid-definition code is becoming repetitive -- it should
  // probably be moved somewhere central...
  double nyDouble = (ymaxJet_ - yminJet_) / gridWidth_;
  nyJet_ = int(nyDouble+0.5);
  dyJet_ = (ymaxJet_ - yminJet_) / nyJet_;
  
  assert(nyJet_ >= 1);

  ntotalJet_ = nphi_ * nyJet_;
  //_scalar_pt.resize(_ntotal);
}


//----------------------------------------------------------------------
// retrieve the grid tile index for a given PseudoJet
int
HiFJGridEmptyAreaCalculator::tileIndexEtaJet(const reco::PFCandidate *pfCand) {
  // directly taking int does not work for values between -1 and 0
  // so use floor instead
  // double iy_double = (p.rap() - _ymin) / _dy;
  // if (iy_double < 0.0) return -1;
  // int iy = int(iy_double);
  // if (iy >= _ny) return -1;

  // writing it as below gives a huge speed gain (factor two!). Even
  // though answers are identical and the routine here is not the
  // speed-critical step. It's not at all clear why.
  int iyjet = int(floor( (pfCand->eta() - yminJet_) / dy_ ));
  if (iyjet < 0 || iyjet >= nyJet_) return -1;
  
  assert (iyjet < nyJet_ && iyjet >= 0);

  return iyjet;
}

int 
HiFJGridEmptyAreaCalculator::numJetGridCells( std::vector<std::pair<int, int> >& indices )
{
  int ngrid = 0;
  //sort phi eta grid indices in phi
  std::sort(indices.begin(),indices.end());
  int lowestJPhi = indices[0].first;
  int lowestJEta = indices[0].second;
  int highestJEta = lowestJEta;
    
  //for each fixed phi value calculate the number of grids in eta
  for(unsigned int iconst = 1; iconst < indices.size(); iconst++){
     int jphi = indices[iconst].first;
     int jeta = indices[iconst].second;
     if (jphi == lowestJPhi){
      if (jeta < lowestJEta) lowestJEta = jeta;
      if (jeta > highestJEta) highestJEta = jeta;
     }else{
      lowestJPhi = jphi;
      ngrid += highestJEta - lowestJEta + 1;
      lowestJEta = jeta;
      highestJEta = jeta;
    }
  }
  ngrid += highestJEta - lowestJEta + 1;
  return ngrid;
}

void
HiFJGridEmptyAreaCalculator::beginStream(edm::StreamID)
{
}

void
HiFJGridEmptyAreaCalculator::endStream() {
}

void HiFJGridEmptyAreaCalculator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<edm::InputTag>("jetSource",edm::InputTag("kt4PFJets"));
  desc.add<edm::InputTag>("CentralityBinSrc",edm::InputTag("centralityBin"));
  desc.add<edm::InputTag>("mapEtaEdges",edm::InputTag("mapEtaEdges"));
  desc.add<edm::InputTag>("mapToRho",edm::InputTag("mapToRho"));
  desc.add<edm::InputTag>("mapToRhoM",edm::InputTag("mapToRhoM"));
  desc.add<edm::InputTag>("pfCandSource",edm::InputTag("particleFlow"));
  desc.add<double>("gridWidth",0.05);
  desc.add<double>("bandWidth",0.2);
  desc.add<bool>("doCentrality", true);
  desc.add<int>("hiBinCut",100);
  desc.add<bool>("keepGridInfo",false);
  descriptions.add("hiFJGridEmptyAreaCalculator",desc);
}

DEFINE_FWK_MODULE(HiFJGridEmptyAreaCalculator);