PATMHTProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    PATMHTProducer
// Class:      PATMHTProducer
//
//
// Original Author:  Xin Shi & Freya Blekman, Cornell University
//         Created:  Fri Sep 12 17:58:29 CEST 2008
//
//

#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/Common/interface/View.h"
#include "DataFormats/METReco/interface/SigInputObj.h"
#include "DataFormats/Math/interface/LorentzVector.h"
#include "DataFormats/PatCandidates/interface/Electron.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/MHT.h"
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/Photon.h"
#include "DataFormats/PatCandidates/interface/Tau.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "RecoMET/METAlgorithms/interface/SignAlgoResolutions.h"
#include "RecoMET/METAlgorithms/interface/significanceAlgo.h"

namespace pat {
  class PATMHTProducer : public edm::stream::EDProducer<> {
  public:
    explicit PATMHTProducer(const edm::ParameterSet&);
    ~PATMHTProducer() override;

  private:
    void produce(edm::Event&, const edm::EventSetup&) override;

    double getJets(edm::Event&, const edm::EventSetup&);
    double getElectrons(edm::Event&, const edm::EventSetup&);
    double getMuons(edm::Event&, const edm::EventSetup&);
    void getTowers(edm::Event&, const edm::EventSetup&);

    // ----------member data ---------------------------

    double verbose_;

    // input tags.
    edm::InputTag mhtLabel_;
    edm::EDGetTokenT<edm::View<pat::Jet> > jetToken_;
    edm::EDGetTokenT<edm::View<pat::Electron> > eleToken_;
    edm::EDGetTokenT<edm::View<pat::Muon> > muoToken_;
    edm::EDGetTokenT<edm::View<pat::Tau> > tauToken_;
    edm::EDGetTokenT<edm::View<pat::Photon> > phoToken_;

    std::vector<metsig::SigInputObj> physobjvector_;

    double uncertaintyScaleFactor_;  // scale factor for the uncertainty parameters.
    bool controlledUncertainty_;     // use controlled uncertainty parameters.

    //--- test the uncertainty parameters ---//

    class uncertaintyFunctions {
    public:
      std::unique_ptr<TF1> etUncertainty;
      std::unique_ptr<TF1> phiUncertainty;
    };

    void setUncertaintyParameters();  // fills the following uncertaintyFunctions objects:
    uncertaintyFunctions ecalEBUncertainty;
    uncertaintyFunctions ecalEEUncertainty;
    uncertaintyFunctions hcalHBUncertainty;
    uncertaintyFunctions hcalHEUncertainty;
    uncertaintyFunctions hcalHOUncertainty;
    uncertaintyFunctions hcalHFUncertainty;

    uncertaintyFunctions jetUncertainty;
    uncertaintyFunctions jetCorrUncertainty;
    uncertaintyFunctions eleUncertainty;
    uncertaintyFunctions muonUncertainty;
    uncertaintyFunctions muonCorrUncertainty;

    //--- tags and parameters ---//

    bool useCaloTowers_;
    bool useJets_;
    bool useElectrons_;
    bool useMuons_;
    std::set<CaloTowerDetId> s_clusteredTowers;

    bool noHF_;

    double jetPtMin_;
    double jetEtaMax_;
    double jetEMfracMax_;

    double elePtMin_;
    double eleEtaMax_;

    double muonPtMin_;
    double muonEtaMax_;
    double muonTrackD0Max_;
    double muonTrackDzMax_;
    int muonNHitsMin_;
    double muonDPtMax_;
    double muonChiSqMax_;

    //  double uncertaintyScaleFactor_; // scale factor for the uncertainty parameters.

    double jetEtUncertaintyParameter0_;
    double jetEtUncertaintyParameter1_;
    double jetEtUncertaintyParameter2_;

    double jetPhiUncertaintyParameter0_;
    double jetPhiUncertaintyParameter1_;
    double jetPhiUncertaintyParameter2_;

    double eleEtUncertaintyParameter0_;
    double elePhiUncertaintyParameter0_;

    double muonEtUncertaintyParameter0_;
    double muonPhiUncertaintyParameter0_;

    edm::InputTag CaloJetAlgorithmTag_;
    edm::InputTag CorJetAlgorithmTag_;
    std::string JetCorrectionService_;
    edm::InputTag MuonTag_;
    edm::InputTag ElectronTag_;
    edm::InputTag CaloTowerTag_;
    std::string metCollectionLabel_;
    std::string significanceLabel_;

    //--- For Muon Calo Deposits ---//
    //TrackDetectorAssociator   trackAssociator_;
    //TrackAssociatorParameters trackAssociatorParameters_;

    double towerEtThreshold_;
    bool useHO_;
  };
  //define this as a plug-in

}  // namespace pat

#include <memory>

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 =
      std::make_unique<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 = std::make_unique<TF1>("ecalEBphiFunc", "[0]*x", 1);
  ecalEBUncertainty.phiUncertainty->SetParameter(0, 0.0174);

  //-- Ecal Endcap --//
  ecalEEUncertainty.etUncertainty =
      std::make_unique<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 = std::make_unique<TF1>("ecalEEphiFunc", "[0]*x", 1);
  ecalEEUncertainty.phiUncertainty->SetParameter(0, 0.087);

  //-- Hcal Uncertainty Functions --------------------------------------//
  //-- From: FastSimulation/Calorimetry/data/HcalResponse.cfi --//
  //-- Hcal Barrel --//
  hcalHBUncertainty.etUncertainty =
      std::make_unique<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 = std::make_unique<TF1>("ecalHBphiFunc", "[0]*x", 1);
  hcalHBUncertainty.phiUncertainty->SetParameter(0, 0.087);

  //-- Hcal Endcap --//
  hcalHEUncertainty.etUncertainty =
      std::make_unique<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 = std::make_unique<TF1>("ecalHEphiFunc", "[0]*x", 1);
  hcalHEUncertainty.phiUncertainty->SetParameter(0, 0.087);

  //-- Hcal Outer --//
  hcalHOUncertainty.etUncertainty =
      std::make_unique<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 = std::make_unique<TF1>("ecalHOphiFunc", "[0]*x", 1);
  hcalHOUncertainty.phiUncertainty->SetParameter(0, 0.087);

  //-- Hcal Forward --//
  hcalHFUncertainty.etUncertainty =
      std::make_unique<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 = std::make_unique<TF1>("ecalHFphiFunc", "[0]*x", 1);
  hcalHFUncertainty.phiUncertainty->SetParameter(0, 0.174);

  //--- Jet Uncertainty Functions --------------------------------------//
  jetUncertainty.etUncertainty = std::make_unique<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 =
      std::make_unique<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 = std::make_unique<TF1>("eleEtFunc", "[0] * x", 1);
  //  eleUncertainty.etUncertainty->SetParameter(0,0.034);
  eleUncertainty.etUncertainty->SetParameter(0, eleEtUncertaintyParameter0_);

  eleUncertainty.phiUncertainty = std::make_unique<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 = std::make_unique<TF1>("muonEtFunc", "[0] * x", 1);
  //  muonUncertainty.etUncertainty->SetParameter(0,0.01);
  muonUncertainty.etUncertainty->SetParameter(0, muonEtUncertaintyParameter0_);
  muonUncertainty.phiUncertainty = std::make_unique<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.)); */
}

#include "FWCore/Framework/interface/MakerMacros.h"
using namespace pat;
DEFINE_FWK_MODULE(PATMHTProducer);