PATMHTProducer.cc

Go to the documentation of this file.

//
//

#include "PhysicsTools/PatAlgos/plugins/PATMHTProducer.h"

pat::PATMHTProducer::PATMHTProducer(const edm::ParameterSet& iConfig) {
  // Initialize the configurables
  verbose_ = iConfig.getParameter<double>("verbose");

  jetToken_ = consumes<edm::View<pat::Jet> >(iConfig.getUntrackedParameter<edm::InputTag>("jetTag"));
  eleToken_ = consumes<edm::View<pat::Electron> >(iConfig.getUntrackedParameter<edm::InputTag>("electronTag"));
  muoToken_ = consumes<edm::View<pat::Muon> >(iConfig.getUntrackedParameter<edm::InputTag>("muonTag"));
  tauToken_ = consumes<edm::View<pat::Tau> >(iConfig.getUntrackedParameter<edm::InputTag>("tauTag"));
  phoToken_ = consumes<edm::View<pat::Photon> >(iConfig.getUntrackedParameter<edm::InputTag>("photonTag"));

  uncertaintyScaleFactor_ = iConfig.getParameter<double>("uncertaintyScaleFactor");
  controlledUncertainty_ = iConfig.getParameter<bool>("controlledUncertainty");

  jetPtMin_ = iConfig.getParameter<double>("jetPtMin");
  jetEtaMax_ = iConfig.getParameter<double>("jetEtaMax");
  jetEMfracMax_ = iConfig.getParameter<double>("jetEMfracMax");
  elePtMin_ = iConfig.getParameter<double>("elePtMin");
  eleEtaMax_ = iConfig.getParameter<double>("eleEtaMax");
  muonPtMin_ = iConfig.getParameter<double>("muonPtMin");
  muonEtaMax_ = iConfig.getParameter<double>("muonEtaMax");

  jetEtUncertaintyParameter0_ = iConfig.getParameter<double>("jetEtUncertaintyParameter0");
  jetEtUncertaintyParameter1_ = iConfig.getParameter<double>("jetEtUncertaintyParameter1");
  jetEtUncertaintyParameter2_ = iConfig.getParameter<double>("jetEtUncertaintyParameter2");
  jetPhiUncertaintyParameter0_ = iConfig.getParameter<double>("jetPhiUncertaintyParameter0");
  jetPhiUncertaintyParameter1_ = iConfig.getParameter<double>("jetPhiUncertaintyParameter1");
  jetPhiUncertaintyParameter2_ = iConfig.getParameter<double>("jetPhiUncertaintyParameter2");

  eleEtUncertaintyParameter0_ = iConfig.getParameter<double>("eleEtUncertaintyParameter0");
  elePhiUncertaintyParameter0_ = iConfig.getParameter<double>("elePhiUncertaintyParameter0");

  muonEtUncertaintyParameter0_ = iConfig.getParameter<double>("muonEtUncertaintyParameter0");
  muonPhiUncertaintyParameter0_ = iConfig.getParameter<double>("muonPhiUncertaintyParameter0");

  CaloTowerTag_ = iConfig.getParameter<edm::InputTag>("CaloTowerTag");
  noHF_ = iConfig.getParameter<bool>("noHF");

  //  muonCalo_ = iConfig.getParameter<bool>("muonCalo");
  towerEtThreshold_ = iConfig.getParameter<double>("towerEtThreshold");
  useHO_ = iConfig.getParameter<bool>("useHO");

  setUncertaintyParameters();

  produces<pat::MHTCollection>();
}

pat::PATMHTProducer::~PATMHTProducer() {}

void pat::PATMHTProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  // make sure the SigInputObj container is empty
  while (!physobjvector_.empty()) {
    physobjvector_.erase(physobjvector_.begin(), physobjvector_.end());
  }

  // Clean the clustered towers
  s_clusteredTowers.clear();

  double number_of_jets = getJets(iEvent, iSetup);

  double number_of_electrons = getElectrons(iEvent, iSetup);

  double number_of_muons = getMuons(iEvent, iSetup);

  if (verbose_ == 1.) {
    std::cout << ">>>---> Number of jets: " << number_of_jets << std::endl;
    std::cout << ">>>---> Number of electrons: " << number_of_jets << std::endl;
    std::cout << ">>>---> Number of muons: " << number_of_muons << std::endl;
  }

  double met_x = 0;
  double met_y = 0;
  double met_et = 0;
  double met_phi = 0;
  double met_set = 0;

  auto themetsigcoll = std::make_unique<pat::MHTCollection>();

  if (!physobjvector_.empty()) {  // Only when the vector is not empty

    // calculate the MHT significance

    metsig::significanceAlgo signifAlgo;
    signifAlgo.addObjects(physobjvector_);
    double significance = signifAlgo.significance(met_et, met_phi, met_set);

    met_x = met_et * cos(met_phi);
    met_y = met_et * sin(met_phi);

    if (verbose_ == 1.) {
      std::cout << ">>>----> MHT Sgificance = " << significance << std::endl;
    }

    pat::MHT themetsigobj(reco::Particle::LorentzVector(met_x, met_y, 0, met_et), met_set, significance);

    // Store the number of jets, electrons, muons
    themetsigobj.setNumberOfJets(number_of_jets);
    themetsigobj.setNumberOfElectrons(number_of_electrons);
    themetsigobj.setNumberOfMuons(number_of_muons);

    themetsigcoll->push_back(themetsigobj);

  }  // If the vector is empty, just put empty product.

  iEvent.put(std::move(themetsigcoll));
}

// --------------------------------------------------
//  Fill Input Vector with Jets
// --------------------------------------------------
double pat::PATMHTProducer::getJets(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  std::string objectname = "jet";

  double number_of_jets_ = 0.0;

  edm::Handle<edm::View<pat::Jet> > jetHandle;
  iEvent.getByToken(jetToken_, jetHandle);
  edm::View<pat::Jet> jets = *jetHandle;

  // Fill Input Vector with Jets
  for (edm::View<pat::Jet>::const_iterator jet_iter = jets.begin(); jet_iter != jets.end(); ++jet_iter) {
    if ((jet_iter->pt() < jetPtMin_) || (TMath::Abs(jet_iter->eta()) > jetEtaMax_) ||
        (jet_iter->emEnergyFraction() > jetEMfracMax_))
      continue;

    double jet_et = jet_iter->et();
    double jet_phi = jet_iter->phi();

    if (verbose_ == 3.) {
      std::cout << "jet pt : " << jet_iter->pt() << " eta : " << jet_iter->eta()
                << " EMF: " << jet_iter->emEnergyFraction() << std::endl;
    }

    double sigma_et, sigma_phi;

    if (controlledUncertainty_) {
      sigma_et = jetUncertainty.etUncertainty->Eval(jet_et);
      sigma_phi = jetUncertainty.phiUncertainty->Eval(jet_et);
    } else {
      sigma_et = 0.0;   // jet_iter->resolutionEt();
      sigma_phi = 0.0;  //jet_iter->resolutionPhi();
    }

    if (verbose_ == 3.) {
      std::cout << "jet sigma_et : " << sigma_et << ", jet sigma_phi : " << sigma_phi << std::endl;
    }

    if (sigma_et <= 0 || sigma_phi <= 0)
      edm::LogWarning("PATMHTProducer") << " uncertainties for " << objectname << " are (et, phi): " << sigma_et << ","
                                        << sigma_phi << " (et,phi): " << jet_et << "," << jet_phi;
    // try to read out the jet resolution from the root file at PatUtils
    //-- Store jet for Significance Calculation --//

    if (uncertaintyScaleFactor_ != 1.0) {
      sigma_et = sigma_et * uncertaintyScaleFactor_;
      sigma_phi = sigma_phi * uncertaintyScaleFactor_;
      // edm::LogWarning("PATMHTProducer") << " using uncertainty scale factor: " << uncertaintyScaleFactor_ <<
      //" , uncertainties for " << objectname <<" changed to (et, phi): " << sigma_et << "," << sigma_phi;
    }

    if (verbose_ == 101.) {  // Study the Jets behavior

      std::cout << "v101> " << number_of_jets_ << "  " << jet_et << "  " << sigma_et << "  " << jet_phi << "  "
                << sigma_phi << std::endl;
    }

    metsig::SigInputObj tmp_jet(objectname, jet_et, jet_phi, sigma_et, sigma_phi);
    physobjvector_.push_back(tmp_jet);
    number_of_jets_++;

    //-- Store tower DetId's to be removed from Calo Tower sum later --//
    std::vector<CaloTowerPtr> v_towers = jet_iter->getCaloConstituents();
    //std::cout << "tower size = " << v_towers.size() << std::endl;

    for (unsigned int ii = 0; ii < v_towers.size(); ii++) {
      s_clusteredTowers.insert((*v_towers.at(ii)).id());
      //std::cout << "tower id = " << (*v_towers.at(ii)).id() << std::endl;
    }
  }

  if (verbose_ == 101.) {  // Study the Jets behavior - seperate events
    std::cout << "v101> --------------------------------------------" << std::endl;
  }

  return number_of_jets_;
}

// --------------------------------------------------
//  Fill Input Vector with Electrons
// --------------------------------------------------
double pat::PATMHTProducer::getElectrons(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  std::string objectname = "electron";

  double number_of_electrons_ = 0.0;

  // edm::ESHandle<CaloTowerConstituentsMap> cttopo;
  // iSetup.get<HcalRecNumberingRecord>().get(cttopo);
  // const CaloTowerConstituentsMap* caloTowerMap = cttopo.product();

  edm::Handle<edm::View<pat::Electron> > electronHandle;
  iEvent.getByToken(eleToken_, electronHandle);
  edm::View<pat::Electron> electrons = *electronHandle;

  // Fill Input Vector with Electrons
  for (edm::View<pat::Electron>::const_iterator electron_iter = electrons.begin(); electron_iter != electrons.end();
       ++electron_iter) {
    // Select electrons
    if (electron_iter->et() < elePtMin_ || TMath::Abs(electron_iter->eta()) > eleEtaMax_)
      continue;

    if (verbose_ == 3.) {
      std::cout << "electron pt = " << electron_iter->pt() << " eta : " << electron_iter->eta() << std::endl;
    }

    double electron_et = electron_iter->et();
    double electron_phi = electron_iter->phi();

    double sigma_et, sigma_phi;

    if (controlledUncertainty_) {
      sigma_et = eleUncertainty.etUncertainty->Eval(electron_et);
      sigma_phi = eleUncertainty.phiUncertainty->Eval(electron_et);
    } else {
      sigma_et = 0.0;   //electron_iter->resolutionEt();
      sigma_phi = 0.0;  // electron_iter->resolutionPhi();
    }

    if (verbose_ == 3.) {
      std::cout << "electron sigma_et : " << sigma_et << ", electron sigma_phi : " << sigma_phi << std::endl;
    }

    if (sigma_et < 0 || sigma_phi < 0)
      edm::LogWarning("PATMHTProducer") << " uncertainties for " << objectname << " are (et, phi): " << sigma_et << ","
                                        << sigma_phi << " (et,phi): " << electron_et << "," << electron_phi;

    if (uncertaintyScaleFactor_ != 1.0) {
      sigma_et = sigma_et * uncertaintyScaleFactor_;
      sigma_phi = sigma_phi * uncertaintyScaleFactor_;
    }

    metsig::SigInputObj tmp_electron(objectname, electron_et, electron_phi, sigma_et, sigma_phi);
    physobjvector_.push_back(tmp_electron);
    number_of_electrons_++;

    //-- Store tower DetId's to be removed from Calo Tower sum later --//
    /*
    const reco::SuperCluster& eleSC = *( electron_iter->superCluster() );

    std::vector<DetId> v_eleDetIds = eleSC.getHitsByDetId();

    //-- Convert cells to calo towers and add to set --//
    for( std::vector<DetId>::iterator cellId = v_eleDetIds.begin();
         cellId != v_eleDetIds.end();
         cellId++) {

      CaloTowerDetId towerId = caloTowerMap->towerOf(*cellId);
      if (towerId != nullDetId) {
    //std::cout << ">>> electron towerId: " << towerId << std::endl;
    std::pair<std::_Rb_tree_const_iterator<CaloTowerDetId>,bool> p1 = s_clusteredTowers.insert(towerId);
      }
      else
    std::cerr<<"No matching tower found for electron cell!\n";
    }

    */
  }

  return number_of_electrons_;
}

// --------------------------------------------------
//  Fill Input Vector with Muons
// --------------------------------------------------

double pat::PATMHTProducer::getMuons(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  std::string objectname = "muon";
  edm::Handle<edm::View<pat::Muon> > muonHandle;
  iEvent.getByToken(muoToken_, muonHandle);
  edm::View<pat::Muon> muons = *muonHandle;

  if (!muonHandle.isValid()) {
    std::cout << ">>> PATMHTSelector not valid muon Handle!" << std::endl;
    return 0.0;
  }

  double number_of_muons_ = 0.0;

  for (edm::View<pat::Muon>::const_iterator muon_iter = muons.begin(); muon_iter != muons.end(); ++muon_iter) {
    if (muon_iter->pt() < muonPtMin_ || TMath::Abs(muon_iter->eta()) > muonEtaMax_)
      continue;

    if (verbose_ == 3.) {
      std::cout << "muon pt = " << muon_iter->pt() << " eta : " << muon_iter->eta() << std::endl;
    }

    double muon_pt = muon_iter->pt();
    double muon_phi = muon_iter->phi();

    double sigma_et, sigma_phi;

    if (controlledUncertainty_) {
      sigma_et = muonUncertainty.etUncertainty->Eval(muon_pt);
      sigma_phi = muonUncertainty.phiUncertainty->Eval(muon_pt);
    } else {
      sigma_et = 0.0;   //muon_iter->resolutionEt();
      sigma_phi = 0.0;  // muon_iter->resolutionPhi();
    }

    if (verbose_ == 3.) {
      std::cout << "muon sigma_et : " << sigma_et << ", muon sigma_phi : " << sigma_phi << std::endl;
    }

    if (sigma_et < 0 || sigma_phi < 0)
      edm::LogWarning("PATMHTProducer") << " uncertainties for " << objectname << " are (et, phi): " << sigma_et << ","
                                        << sigma_phi << " (pt,phi): " << muon_pt << "," << muon_phi;

    if (uncertaintyScaleFactor_ != 1.0) {
      sigma_et = sigma_et * uncertaintyScaleFactor_;
      sigma_phi = sigma_phi * uncertaintyScaleFactor_;
    }

    metsig::SigInputObj tmp_muon(objectname, muon_pt, muon_phi, sigma_et, sigma_phi);
    physobjvector_.push_back(tmp_muon);
    number_of_muons_++;

  }  // end Muon loop

  return number_of_muons_;
}

//=== Uncertainty Functions ===============================================
void pat::PATMHTProducer::setUncertaintyParameters() {
  // set the various functions here:

  //-- For Et functions, [0]= par_n, [1]=par_s, [2]= par_c ---//
  //-- Ecal Uncertainty Functions ------------------------------------//
  //-- From: FastSimulation/Calorimetry/data/HcalResponse.cfi --//
  //-- Ecal Barrel --//
  ecalEBUncertainty.etUncertainty.reset(new TF1("ecalEBEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3));
  ecalEBUncertainty.etUncertainty->SetParameter(0, 0.2);
  ecalEBUncertainty.etUncertainty->SetParameter(1, 0.03);
  ecalEBUncertainty.etUncertainty->SetParameter(2, 0.005);

  ecalEBUncertainty.phiUncertainty.reset(new TF1("ecalEBphiFunc", "[0]*x", 1));
  ecalEBUncertainty.phiUncertainty->SetParameter(0, 0.0174);

  //-- Ecal Endcap --//
  ecalEEUncertainty.etUncertainty.reset(new TF1("ecalEEEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3));
  ecalEEUncertainty.etUncertainty->SetParameter(0, 0.2);
  ecalEEUncertainty.etUncertainty->SetParameter(1, 0.03);
  ecalEEUncertainty.etUncertainty->SetParameter(2, 0.005);

  ecalEEUncertainty.phiUncertainty.reset(new TF1("ecalEEphiFunc", "[0]*x", 1));
  ecalEEUncertainty.phiUncertainty->SetParameter(0, 0.087);

  //-- Hcal Uncertainty Functions --------------------------------------//
  //-- From: FastSimulation/Calorimetry/data/HcalResponse.cfi --//
  //-- Hcal Barrel --//
  hcalHBUncertainty.etUncertainty.reset(new TF1("hcalHBEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3));
  hcalHBUncertainty.etUncertainty->SetParameter(0, 0.);
  hcalHBUncertainty.etUncertainty->SetParameter(1, 1.22);
  hcalHBUncertainty.etUncertainty->SetParameter(2, 0.05);

  hcalHBUncertainty.phiUncertainty.reset(new TF1("ecalHBphiFunc", "[0]*x", 1));
  hcalHBUncertainty.phiUncertainty->SetParameter(0, 0.087);

  //-- Hcal Endcap --//
  hcalHEUncertainty.etUncertainty.reset(new TF1("hcalHEEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3));
  hcalHEUncertainty.etUncertainty->SetParameter(0, 0.);
  hcalHEUncertainty.etUncertainty->SetParameter(1, 1.3);
  hcalHEUncertainty.etUncertainty->SetParameter(2, 0.05);

  hcalHEUncertainty.phiUncertainty.reset(new TF1("ecalHEphiFunc", "[0]*x", 1));
  hcalHEUncertainty.phiUncertainty->SetParameter(0, 0.087);

  //-- Hcal Outer --//
  hcalHOUncertainty.etUncertainty.reset(new TF1("hcalHOEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3));
  hcalHOUncertainty.etUncertainty->SetParameter(0, 0.);
  hcalHOUncertainty.etUncertainty->SetParameter(1, 1.82);
  hcalHOUncertainty.etUncertainty->SetParameter(2, 0.09);

  hcalHOUncertainty.phiUncertainty.reset(new TF1("ecalHOphiFunc", "[0]*x", 1));
  hcalHOUncertainty.phiUncertainty->SetParameter(0, 0.087);

  //-- Hcal Forward --//
  hcalHFUncertainty.etUncertainty.reset(new TF1("hcalHFEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3));
  hcalHFUncertainty.etUncertainty->SetParameter(0, 0.);
  hcalHFUncertainty.etUncertainty->SetParameter(1, 1.82);
  hcalHFUncertainty.etUncertainty->SetParameter(2, 0.09);

  hcalHFUncertainty.phiUncertainty.reset(new TF1("ecalHFphiFunc", "[0]*x", 1));
  hcalHFUncertainty.phiUncertainty->SetParameter(0, 0.174);

  //--- Jet Uncertainty Functions --------------------------------------//
  jetUncertainty.etUncertainty.reset(new TF1("jetEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3));
  //-- values from PTDR 1, ch 11.4 --//
  jetUncertainty.etUncertainty->SetParameter(0, jetEtUncertaintyParameter0_);
  jetUncertainty.etUncertainty->SetParameter(1, jetEtUncertaintyParameter1_);
  jetUncertainty.etUncertainty->SetParameter(2, jetEtUncertaintyParameter2_);

  //-- phi value from our own fits --//
  //jetUncertainty.phiUncertainty.reset( new TF1("jetPhiFunc","[0]*x",1) );
  //jetUncertainty.phiUncertainty->SetParameter(0, jetPhiUncertaintyParameter0_);

  //-- phi Functions and values from
  // http://indico.cern.ch/getFile.py/access?contribId=9&sessionId=0&resId=0&materialId=slides&confId=46394
  jetUncertainty.phiUncertainty.reset(new TF1("jetPhiFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3));
  jetUncertainty.phiUncertainty->SetParameter(0, jetPhiUncertaintyParameter0_);
  jetUncertainty.phiUncertainty->SetParameter(1, jetPhiUncertaintyParameter1_);
  jetUncertainty.phiUncertainty->SetParameter(2, jetPhiUncertaintyParameter2_);

  //-- Jet corrections are assumed not to have an error --//
  /*jetCorrUncertainty.etUncertainty.reset( new TF1("jetCorrEtFunc","[0]*x",1) );
  jetCorrUncertainty.etUncertainty->SetParameter(0,0.0);
  jetCorrUncertainty.phiUncertainty.reset( new TF1("jetCorrPhiFunc","[0]*x",1) );
  jetCorrUncertainty.phiUncertainty->SetParameter(0,0.0*(3.14159/180.));*/

  //--- Electron Uncertainty Functions ---------------------------------//
  // completely ambiguious values for electron-like jets...
  // the egamma group keeps track of these here:
  // https://twiki.cern.ch/twiki/bin/view/CMS/EgammaCMSSWVal
  // electron resolution in energy is around 3.4%, measured for 10 < pT < 50 at realistic events with pile-up.

  eleUncertainty.etUncertainty.reset(new TF1("eleEtFunc", "[0] * x", 1));
  //  eleUncertainty.etUncertainty->SetParameter(0,0.034);
  eleUncertainty.etUncertainty->SetParameter(0, eleEtUncertaintyParameter0_);

  eleUncertainty.phiUncertainty.reset(new TF1("elePhiFunc", "[0] * x", 1));
  //  eleUncertainty.phiUncertainty->SetParameter(0,1*(3.14159/180.));
  eleUncertainty.phiUncertainty->SetParameter(0, elePhiUncertaintyParameter0_);

  //--- Muon Uncertainty Functions ------------------------------------//
  // and ambiguious values for the muons...

  muonUncertainty.etUncertainty.reset(new TF1("muonEtFunc", "[0] * x", 1));
  //  muonUncertainty.etUncertainty->SetParameter(0,0.01);
  muonUncertainty.etUncertainty->SetParameter(0, muonEtUncertaintyParameter0_);
  muonUncertainty.phiUncertainty.reset(new TF1("muonPhiFunc", "[0] * x", 1));
  //  muonUncertainty.phiUncertainty->SetParameter(0,1*(3.14159/180.));
  muonUncertainty.phiUncertainty->SetParameter(0, muonPhiUncertaintyParameter0_);

  //-- Muon calo deposites are assumed not to have an error --//
  /*muonCorrUncertainty.etUncertainty.reset( new TF1("muonCorrEtFunc","[0] * x",1) );
  muonCorrUncertainty.etUncertainty->SetParameter(0,0.0);
  muonCorrUncertainty.phiUncertainty.reset( new TF1("muonCorrPhiFunc","[0] * x",1) );
  muonCorrUncertainty.phiUncertainty->SetParameter(0,0.0*(3.14159/180.)); */
}

using namespace pat;
DEFINE_FWK_MODULE(PATMHTProducer);