PatAlgos/plugins/PATMHTProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PATMHTProducer
0004 // Class:      PATMHTProducer
0005 //
0006 /**\class PATMHTProducer
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Xin Shi & Freya Blekman, Cornell University
0015 //         Created:  Fri Sep 12 17:58:29 CEST 2008
0016 //
0017 //
0018
0019 #include "DataFormats/Candidate/interface/Candidate.h"
0020 #include "DataFormats/Common/interface/View.h"
0021 #include "DataFormats/METReco/interface/SigInputObj.h"
0022 #include "DataFormats/Math/interface/LorentzVector.h"
0023 #include "DataFormats/PatCandidates/interface/Electron.h"
0024 #include "DataFormats/PatCandidates/interface/Jet.h"
0025 #include "DataFormats/PatCandidates/interface/MHT.h"
0026 #include "DataFormats/PatCandidates/interface/Muon.h"
0027 #include "DataFormats/PatCandidates/interface/Photon.h"
0028 #include "DataFormats/PatCandidates/interface/Tau.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/stream/EDProducer.h"
0031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0032 #include "FWCore/ParameterSet/interface/FileInPath.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 #include "FWCore/Utilities/interface/InputTag.h"
0035 #include "RecoMET/METAlgorithms/interface/SignAlgoResolutions.h"
0036 #include "RecoMET/METAlgorithms/interface/significanceAlgo.h"
0037
0038 namespace pat {
0039   class PATMHTProducer : public edm::stream::EDProducer<> {
0040   public:
0041     explicit PATMHTProducer(const edm::ParameterSet&);
0042     ~PATMHTProducer() override;
0043
0044   private:
0045     void produce(edm::Event&, const edm::EventSetup&) override;
0046
0047     double getJets(edm::Event&, const edm::EventSetup&);
0048     double getElectrons(edm::Event&, const edm::EventSetup&);
0049     double getMuons(edm::Event&, const edm::EventSetup&);
0050     void getTowers(edm::Event&, const edm::EventSetup&);
0051
0052     // ----------member data ---------------------------
0053
0054     double verbose_;
0055
0056     // input tags.
0057     edm::InputTag mhtLabel_;
0058     edm::EDGetTokenT<edm::View<pat::Jet> > jetToken_;
0059     edm::EDGetTokenT<edm::View<pat::Electron> > eleToken_;
0060     edm::EDGetTokenT<edm::View<pat::Muon> > muoToken_;
0061     edm::EDGetTokenT<edm::View<pat::Tau> > tauToken_;
0062     edm::EDGetTokenT<edm::View<pat::Photon> > phoToken_;
0063
0064     std::vector<metsig::SigInputObj> physobjvector_;
0065
0066     double uncertaintyScaleFactor_;  // scale factor for the uncertainty parameters.
0067     bool controlledUncertainty_;     // use controlled uncertainty parameters.
0068
0069     //--- test the uncertainty parameters ---//
0070
0071     class uncertaintyFunctions {
0072     public:
0073       std::unique_ptr<TF1> etUncertainty;
0074       std::unique_ptr<TF1> phiUncertainty;
0075     };
0076
0077     void setUncertaintyParameters();  // fills the following uncertaintyFunctions objects:
0078     uncertaintyFunctions ecalEBUncertainty;
0079     uncertaintyFunctions ecalEEUncertainty;
0080     uncertaintyFunctions hcalHBUncertainty;
0081     uncertaintyFunctions hcalHEUncertainty;
0082     uncertaintyFunctions hcalHOUncertainty;
0083     uncertaintyFunctions hcalHFUncertainty;
0084
0085     uncertaintyFunctions jetUncertainty;
0086     uncertaintyFunctions jetCorrUncertainty;
0087     uncertaintyFunctions eleUncertainty;
0088     uncertaintyFunctions muonUncertainty;
0089     uncertaintyFunctions muonCorrUncertainty;
0090
0091     //--- tags and parameters ---//
0092
0093     bool useCaloTowers_;
0094     bool useJets_;
0095     bool useElectrons_;
0096     bool useMuons_;
0097     std::set<CaloTowerDetId> s_clusteredTowers;
0098
0099     bool noHF_;
0100
0101     double jetPtMin_;
0102     double jetEtaMax_;
0103     double jetEMfracMax_;
0104
0105     double elePtMin_;
0106     double eleEtaMax_;
0107
0108     double muonPtMin_;
0109     double muonEtaMax_;
0110     double muonTrackD0Max_;
0111     double muonTrackDzMax_;
0112     int muonNHitsMin_;
0113     double muonDPtMax_;
0114     double muonChiSqMax_;
0115
0116     //  double uncertaintyScaleFactor_; // scale factor for the uncertainty parameters.
0117
0118     double jetEtUncertaintyParameter0_;
0119     double jetEtUncertaintyParameter1_;
0120     double jetEtUncertaintyParameter2_;
0121
0122     double jetPhiUncertaintyParameter0_;
0123     double jetPhiUncertaintyParameter1_;
0124     double jetPhiUncertaintyParameter2_;
0125
0126     double eleEtUncertaintyParameter0_;
0127     double elePhiUncertaintyParameter0_;
0128
0129     double muonEtUncertaintyParameter0_;
0130     double muonPhiUncertaintyParameter0_;
0131
0132     edm::InputTag CaloJetAlgorithmTag_;
0133     edm::InputTag CorJetAlgorithmTag_;
0134     std::string JetCorrectionService_;
0135     edm::InputTag MuonTag_;
0136     edm::InputTag ElectronTag_;
0137     edm::InputTag CaloTowerTag_;
0138     std::string metCollectionLabel_;
0139     std::string significanceLabel_;
0140
0141     //--- For Muon Calo Deposits ---//
0142     //TrackDetectorAssociator   trackAssociator_;
0143     //TrackAssociatorParameters trackAssociatorParameters_;
0144
0145     double towerEtThreshold_;
0146     bool useHO_;
0147   };
0148   //define this as a plug-in
0149
0150 }  // namespace pat
0151
0152 #include <memory>
0153
0154 pat::PATMHTProducer::PATMHTProducer(const edm::ParameterSet& iConfig) {
0155   // Initialize the configurables
0156   verbose_ = iConfig.getParameter<double>("verbose");
0157
0158   jetToken_ = consumes<edm::View<pat::Jet> >(iConfig.getUntrackedParameter<edm::InputTag>("jetTag"));
0159   eleToken_ = consumes<edm::View<pat::Electron> >(iConfig.getUntrackedParameter<edm::InputTag>("electronTag"));
0160   muoToken_ = consumes<edm::View<pat::Muon> >(iConfig.getUntrackedParameter<edm::InputTag>("muonTag"));
0161   tauToken_ = consumes<edm::View<pat::Tau> >(iConfig.getUntrackedParameter<edm::InputTag>("tauTag"));
0162   phoToken_ = consumes<edm::View<pat::Photon> >(iConfig.getUntrackedParameter<edm::InputTag>("photonTag"));
0163
0164   uncertaintyScaleFactor_ = iConfig.getParameter<double>("uncertaintyScaleFactor");
0165   controlledUncertainty_ = iConfig.getParameter<bool>("controlledUncertainty");
0166
0167   jetPtMin_ = iConfig.getParameter<double>("jetPtMin");
0168   jetEtaMax_ = iConfig.getParameter<double>("jetEtaMax");
0169   jetEMfracMax_ = iConfig.getParameter<double>("jetEMfracMax");
0170   elePtMin_ = iConfig.getParameter<double>("elePtMin");
0171   eleEtaMax_ = iConfig.getParameter<double>("eleEtaMax");
0172   muonPtMin_ = iConfig.getParameter<double>("muonPtMin");
0173   muonEtaMax_ = iConfig.getParameter<double>("muonEtaMax");
0174
0175   jetEtUncertaintyParameter0_ = iConfig.getParameter<double>("jetEtUncertaintyParameter0");
0176   jetEtUncertaintyParameter1_ = iConfig.getParameter<double>("jetEtUncertaintyParameter1");
0177   jetEtUncertaintyParameter2_ = iConfig.getParameter<double>("jetEtUncertaintyParameter2");
0178   jetPhiUncertaintyParameter0_ = iConfig.getParameter<double>("jetPhiUncertaintyParameter0");
0179   jetPhiUncertaintyParameter1_ = iConfig.getParameter<double>("jetPhiUncertaintyParameter1");
0180   jetPhiUncertaintyParameter2_ = iConfig.getParameter<double>("jetPhiUncertaintyParameter2");
0181
0182   eleEtUncertaintyParameter0_ = iConfig.getParameter<double>("eleEtUncertaintyParameter0");
0183   elePhiUncertaintyParameter0_ = iConfig.getParameter<double>("elePhiUncertaintyParameter0");
0184
0185   muonEtUncertaintyParameter0_ = iConfig.getParameter<double>("muonEtUncertaintyParameter0");
0186   muonPhiUncertaintyParameter0_ = iConfig.getParameter<double>("muonPhiUncertaintyParameter0");
0187
0188   CaloTowerTag_ = iConfig.getParameter<edm::InputTag>("CaloTowerTag");
0189   noHF_ = iConfig.getParameter<bool>("noHF");
0190
0191   //  muonCalo_ = iConfig.getParameter<bool>("muonCalo");
0192   towerEtThreshold_ = iConfig.getParameter<double>("towerEtThreshold");
0193   useHO_ = iConfig.getParameter<bool>("useHO");
0194
0195   setUncertaintyParameters();
0196
0197   produces<pat::MHTCollection>();
0198 }
0199
0200 pat::PATMHTProducer::~PATMHTProducer() {}
0201
0202 void pat::PATMHTProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0203   // make sure the SigInputObj container is empty
0204   while (!physobjvector_.empty()) {
0205     physobjvector_.erase(physobjvector_.begin(), physobjvector_.end());
0206   }
0207
0208   // Clean the clustered towers
0209   s_clusteredTowers.clear();
0210
0211   double number_of_jets = getJets(iEvent, iSetup);
0212
0213   double number_of_electrons = getElectrons(iEvent, iSetup);
0214
0215   double number_of_muons = getMuons(iEvent, iSetup);
0216
0217   if (verbose_ == 1.) {
0218     std::cout << ">>>---> Number of jets: " << number_of_jets << std::endl;
0219     std::cout << ">>>---> Number of electrons: " << number_of_jets << std::endl;
0220     std::cout << ">>>---> Number of muons: " << number_of_muons << std::endl;
0221   }
0222
0223   double met_x = 0;
0224   double met_y = 0;
0225   double met_et = 0;
0226   double met_phi = 0;
0227   double met_set = 0;
0228
0229   auto themetsigcoll = std::make_unique<pat::MHTCollection>();
0230
0231   if (!physobjvector_.empty()) {  // Only when the vector is not empty
0232
0233     // calculate the MHT significance
0234
0235     metsig::significanceAlgo signifAlgo;
0236     signifAlgo.addObjects(physobjvector_);
0237     double significance = signifAlgo.significance(met_et, met_phi, met_set);
0238
0239     met_x = met_et * cos(met_phi);
0240     met_y = met_et * sin(met_phi);
0241
0242     if (verbose_ == 1.) {
0243       std::cout << ">>>----> MHT Sgificance = " << significance << std::endl;
0244     }
0245
0246     pat::MHT themetsigobj(reco::Particle::LorentzVector(met_x, met_y, 0, met_et), met_set, significance);
0247
0248     // Store the number of jets, electrons, muons
0249     themetsigobj.setNumberOfJets(number_of_jets);
0250     themetsigobj.setNumberOfElectrons(number_of_electrons);
0251     themetsigobj.setNumberOfMuons(number_of_muons);
0252
0253     themetsigcoll->push_back(themetsigobj);
0254
0255   }  // If the vector is empty, just put empty product.
0256
0257   iEvent.put(std::move(themetsigcoll));
0258 }
0259
0260 // --------------------------------------------------
0261 //  Fill Input Vector with Jets
0262 // --------------------------------------------------
0263 double pat::PATMHTProducer::getJets(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0264   std::string objectname = "jet";
0265
0266   double number_of_jets_ = 0.0;
0267
0268   edm::Handle<edm::View<pat::Jet> > jetHandle;
0269   iEvent.getByToken(jetToken_, jetHandle);
0270   edm::View<pat::Jet> jets = *jetHandle;
0271
0272   // Fill Input Vector with Jets
0273   for (edm::View<pat::Jet>::const_iterator jet_iter = jets.begin(); jet_iter != jets.end(); ++jet_iter) {
0274     if ((jet_iter->pt() < jetPtMin_) || (TMath::Abs(jet_iter->eta()) > jetEtaMax_) ||
0275         (jet_iter->emEnergyFraction() > jetEMfracMax_))
0276       continue;
0277
0278     double jet_et = jet_iter->et();
0279     double jet_phi = jet_iter->phi();
0280
0281     if (verbose_ == 3.) {
0282       std::cout << "jet pt : " << jet_iter->pt() << " eta : " << jet_iter->eta()
0283                 << " EMF: " << jet_iter->emEnergyFraction() << std::endl;
0284     }
0285
0286     double sigma_et, sigma_phi;
0287
0288     if (controlledUncertainty_) {
0289       sigma_et = jetUncertainty.etUncertainty->Eval(jet_et);
0290       sigma_phi = jetUncertainty.phiUncertainty->Eval(jet_et);
0291     } else {
0292       sigma_et = 0.0;   // jet_iter->resolutionEt();
0293       sigma_phi = 0.0;  //jet_iter->resolutionPhi();
0294     }
0295
0296     if (verbose_ == 3.) {
0297       std::cout << "jet sigma_et : " << sigma_et << ", jet sigma_phi : " << sigma_phi << std::endl;
0298     }
0299
0300     if (sigma_et <= 0 || sigma_phi <= 0)
0301       edm::LogWarning("PATMHTProducer") << " uncertainties for " << objectname << " are (et, phi): " << sigma_et << ","
0302                                         << sigma_phi << " (et,phi): " << jet_et << "," << jet_phi;
0303     // try to read out the jet resolution from the root file at PatUtils
0304     //-- Store jet for Significance Calculation --//
0305
0306     if (uncertaintyScaleFactor_ != 1.0) {
0307       sigma_et = sigma_et * uncertaintyScaleFactor_;
0308       sigma_phi = sigma_phi * uncertaintyScaleFactor_;
0309       // edm::LogWarning("PATMHTProducer") << " using uncertainty scale factor: " << uncertaintyScaleFactor_ <<
0310       //" , uncertainties for " << objectname <<" changed to (et, phi): " << sigma_et << "," << sigma_phi;
0311     }
0312
0313     if (verbose_ == 101.) {  // Study the Jets behavior
0314
0315       std::cout << "v101> " << number_of_jets_ << "  " << jet_et << "  " << sigma_et << "  " << jet_phi << "  "
0316                 << sigma_phi << std::endl;
0317     }
0318
0319     metsig::SigInputObj tmp_jet(objectname, jet_et, jet_phi, sigma_et, sigma_phi);
0320     physobjvector_.push_back(tmp_jet);
0321     number_of_jets_++;
0322
0323     //-- Store tower DetId's to be removed from Calo Tower sum later --//
0324     std::vector<CaloTowerPtr> v_towers = jet_iter->getCaloConstituents();
0325     //std::cout << "tower size = " << v_towers.size() << std::endl;
0326
0327     for (unsigned int ii = 0; ii < v_towers.size(); ii++) {
0328       s_clusteredTowers.insert((*v_towers.at(ii)).id());
0329       //std::cout << "tower id = " << (*v_towers.at(ii)).id() << std::endl;
0330     }
0331   }
0332
0333   if (verbose_ == 101.) {  // Study the Jets behavior - seperate events
0334     std::cout << "v101> --------------------------------------------" << std::endl;
0335   }
0336
0337   return number_of_jets_;
0338 }
0339
0340 // --------------------------------------------------
0341 //  Fill Input Vector with Electrons
0342 // --------------------------------------------------
0343 double pat::PATMHTProducer::getElectrons(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0344   std::string objectname = "electron";
0345
0346   double number_of_electrons_ = 0.0;
0347
0348   // edm::ESHandle<CaloTowerConstituentsMap> cttopo;
0349   // iSetup.get<HcalRecNumberingRecord>().get(cttopo);
0350   // const CaloTowerConstituentsMap* caloTowerMap = cttopo.product();
0351
0352   edm::Handle<edm::View<pat::Electron> > electronHandle;
0353   iEvent.getByToken(eleToken_, electronHandle);
0354   edm::View<pat::Electron> electrons = *electronHandle;
0355
0356   // Fill Input Vector with Electrons
0357   for (edm::View<pat::Electron>::const_iterator electron_iter = electrons.begin(); electron_iter != electrons.end();
0358        ++electron_iter) {
0359     // Select electrons
0360     if (electron_iter->et() < elePtMin_ || TMath::Abs(electron_iter->eta()) > eleEtaMax_)
0361       continue;
0362
0363     if (verbose_ == 3.) {
0364       std::cout << "electron pt = " << electron_iter->pt() << " eta : " << electron_iter->eta() << std::endl;
0365     }
0366
0367     double electron_et = electron_iter->et();
0368     double electron_phi = electron_iter->phi();
0369
0370     double sigma_et, sigma_phi;
0371
0372     if (controlledUncertainty_) {
0373       sigma_et = eleUncertainty.etUncertainty->Eval(electron_et);
0374       sigma_phi = eleUncertainty.phiUncertainty->Eval(electron_et);
0375     } else {
0376       sigma_et = 0.0;   //electron_iter->resolutionEt();
0377       sigma_phi = 0.0;  // electron_iter->resolutionPhi();
0378     }
0379
0380     if (verbose_ == 3.) {
0381       std::cout << "electron sigma_et : " << sigma_et << ", electron sigma_phi : " << sigma_phi << std::endl;
0382     }
0383
0384     if (sigma_et < 0 || sigma_phi < 0)
0385       edm::LogWarning("PATMHTProducer") << " uncertainties for " << objectname << " are (et, phi): " << sigma_et << ","
0386                                         << sigma_phi << " (et,phi): " << electron_et << "," << electron_phi;
0387
0388     if (uncertaintyScaleFactor_ != 1.0) {
0389       sigma_et = sigma_et * uncertaintyScaleFactor_;
0390       sigma_phi = sigma_phi * uncertaintyScaleFactor_;
0391     }
0392
0393     metsig::SigInputObj tmp_electron(objectname, electron_et, electron_phi, sigma_et, sigma_phi);
0394     physobjvector_.push_back(tmp_electron);
0395     number_of_electrons_++;
0396
0397     //-- Store tower DetId's to be removed from Calo Tower sum later --//
0398     /*
0399     const reco::SuperCluster& eleSC = *( electron_iter->superCluster() );
0400
0401     std::vector<DetId> v_eleDetIds = eleSC.getHitsByDetId();
0402
0403     //-- Convert cells to calo towers and add to set --//
0404     for( std::vector<DetId>::iterator cellId = v_eleDetIds.begin();
0405          cellId != v_eleDetIds.end();
0406          cellId++) {
0407
0408       CaloTowerDetId towerId = caloTowerMap->towerOf(*cellId);
0409       if (towerId != nullDetId) {
0410     //std::cout << ">>> electron towerId: " << towerId << std::endl;
0411     std::pair<std::_Rb_tree_const_iterator<CaloTowerDetId>,bool> p1 = s_clusteredTowers.insert(towerId);
0412       }
0413       else
0414     std::cerr<<"No matching tower found for electron cell!\n";
0415     }
0416
0417     */
0418   }
0419
0420   return number_of_electrons_;
0421 }
0422
0423 // --------------------------------------------------
0424 //  Fill Input Vector with Muons
0425 // --------------------------------------------------
0426
0427 double pat::PATMHTProducer::getMuons(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0428   std::string objectname = "muon";
0429   edm::Handle<edm::View<pat::Muon> > muonHandle;
0430   iEvent.getByToken(muoToken_, muonHandle);
0431   edm::View<pat::Muon> muons = *muonHandle;
0432
0433   if (!muonHandle.isValid()) {
0434     std::cout << ">>> PATMHTSelector not valid muon Handle!" << std::endl;
0435     return 0.0;
0436   }
0437
0438   double number_of_muons_ = 0.0;
0439
0440   for (edm::View<pat::Muon>::const_iterator muon_iter = muons.begin(); muon_iter != muons.end(); ++muon_iter) {
0441     if (muon_iter->pt() < muonPtMin_ || TMath::Abs(muon_iter->eta()) > muonEtaMax_)
0442       continue;
0443
0444     if (verbose_ == 3.) {
0445       std::cout << "muon pt = " << muon_iter->pt() << " eta : " << muon_iter->eta() << std::endl;
0446     }
0447
0448     double muon_pt = muon_iter->pt();
0449     double muon_phi = muon_iter->phi();
0450
0451     double sigma_et, sigma_phi;
0452
0453     if (controlledUncertainty_) {
0454       sigma_et = muonUncertainty.etUncertainty->Eval(muon_pt);
0455       sigma_phi = muonUncertainty.phiUncertainty->Eval(muon_pt);
0456     } else {
0457       sigma_et = 0.0;   //muon_iter->resolutionEt();
0458       sigma_phi = 0.0;  // muon_iter->resolutionPhi();
0459     }
0460
0461     if (verbose_ == 3.) {
0462       std::cout << "muon sigma_et : " << sigma_et << ", muon sigma_phi : " << sigma_phi << std::endl;
0463     }
0464
0465     if (sigma_et < 0 || sigma_phi < 0)
0466       edm::LogWarning("PATMHTProducer") << " uncertainties for " << objectname << " are (et, phi): " << sigma_et << ","
0467                                         << sigma_phi << " (pt,phi): " << muon_pt << "," << muon_phi;
0468
0469     if (uncertaintyScaleFactor_ != 1.0) {
0470       sigma_et = sigma_et * uncertaintyScaleFactor_;
0471       sigma_phi = sigma_phi * uncertaintyScaleFactor_;
0472     }
0473
0474     metsig::SigInputObj tmp_muon(objectname, muon_pt, muon_phi, sigma_et, sigma_phi);
0475     physobjvector_.push_back(tmp_muon);
0476     number_of_muons_++;
0477
0478   }  // end Muon loop
0479
0480   return number_of_muons_;
0481 }
0482
0483 //=== Uncertainty Functions ===============================================
0484 void pat::PATMHTProducer::setUncertaintyParameters() {
0485   // set the various functions here:
0486
0487   //-- For Et functions, [0]= par_n, [1]=par_s, [2]= par_c ---//
0488   //-- Ecal Uncertainty Functions ------------------------------------//
0489   //-- From: FastSimulation/Calorimetry/data/HcalResponse.cfi --//
0490   //-- Ecal Barrel --//
0491   ecalEBUncertainty.etUncertainty =
0492       std::make_unique<TF1>("ecalEBEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3);
0493   ecalEBUncertainty.etUncertainty->SetParameter(0, 0.2);
0494   ecalEBUncertainty.etUncertainty->SetParameter(1, 0.03);
0495   ecalEBUncertainty.etUncertainty->SetParameter(2, 0.005);
0496
0497   ecalEBUncertainty.phiUncertainty = std::make_unique<TF1>("ecalEBphiFunc", "[0]*x", 1);
0498   ecalEBUncertainty.phiUncertainty->SetParameter(0, 0.0174);
0499
0500   //-- Ecal Endcap --//
0501   ecalEEUncertainty.etUncertainty =
0502       std::make_unique<TF1>("ecalEEEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3);
0503   ecalEEUncertainty.etUncertainty->SetParameter(0, 0.2);
0504   ecalEEUncertainty.etUncertainty->SetParameter(1, 0.03);
0505   ecalEEUncertainty.etUncertainty->SetParameter(2, 0.005);
0506
0507   ecalEEUncertainty.phiUncertainty = std::make_unique<TF1>("ecalEEphiFunc", "[0]*x", 1);
0508   ecalEEUncertainty.phiUncertainty->SetParameter(0, 0.087);
0509
0510   //-- Hcal Uncertainty Functions --------------------------------------//
0511   //-- From: FastSimulation/Calorimetry/data/HcalResponse.cfi --//
0512   //-- Hcal Barrel --//
0513   hcalHBUncertainty.etUncertainty =
0514       std::make_unique<TF1>("hcalHBEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3);
0515   hcalHBUncertainty.etUncertainty->SetParameter(0, 0.);
0516   hcalHBUncertainty.etUncertainty->SetParameter(1, 1.22);
0517   hcalHBUncertainty.etUncertainty->SetParameter(2, 0.05);
0518
0519   hcalHBUncertainty.phiUncertainty = std::make_unique<TF1>("ecalHBphiFunc", "[0]*x", 1);
0520   hcalHBUncertainty.phiUncertainty->SetParameter(0, 0.087);
0521
0522   //-- Hcal Endcap --//
0523   hcalHEUncertainty.etUncertainty =
0524       std::make_unique<TF1>("hcalHEEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3);
0525   hcalHEUncertainty.etUncertainty->SetParameter(0, 0.);
0526   hcalHEUncertainty.etUncertainty->SetParameter(1, 1.3);
0527   hcalHEUncertainty.etUncertainty->SetParameter(2, 0.05);
0528
0529   hcalHEUncertainty.phiUncertainty = std::make_unique<TF1>("ecalHEphiFunc", "[0]*x", 1);
0530   hcalHEUncertainty.phiUncertainty->SetParameter(0, 0.087);
0531
0532   //-- Hcal Outer --//
0533   hcalHOUncertainty.etUncertainty =
0534       std::make_unique<TF1>("hcalHOEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3);
0535   hcalHOUncertainty.etUncertainty->SetParameter(0, 0.);
0536   hcalHOUncertainty.etUncertainty->SetParameter(1, 1.82);
0537   hcalHOUncertainty.etUncertainty->SetParameter(2, 0.09);
0538
0539   hcalHOUncertainty.phiUncertainty = std::make_unique<TF1>("ecalHOphiFunc", "[0]*x", 1);
0540   hcalHOUncertainty.phiUncertainty->SetParameter(0, 0.087);
0541
0542   //-- Hcal Forward --//
0543   hcalHFUncertainty.etUncertainty =
0544       std::make_unique<TF1>("hcalHFEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3);
0545   hcalHFUncertainty.etUncertainty->SetParameter(0, 0.);
0546   hcalHFUncertainty.etUncertainty->SetParameter(1, 1.82);
0547   hcalHFUncertainty.etUncertainty->SetParameter(2, 0.09);
0548
0549   hcalHFUncertainty.phiUncertainty = std::make_unique<TF1>("ecalHFphiFunc", "[0]*x", 1);
0550   hcalHFUncertainty.phiUncertainty->SetParameter(0, 0.174);
0551
0552   //--- Jet Uncertainty Functions --------------------------------------//
0553   jetUncertainty.etUncertainty = std::make_unique<TF1>("jetEtFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3);
0554   //-- values from PTDR 1, ch 11.4 --//
0555   jetUncertainty.etUncertainty->SetParameter(0, jetEtUncertaintyParameter0_);
0556   jetUncertainty.etUncertainty->SetParameter(1, jetEtUncertaintyParameter1_);
0557   jetUncertainty.etUncertainty->SetParameter(2, jetEtUncertaintyParameter2_);
0558
0559   //-- phi value from our own fits --//
0560   //jetUncertainty.phiUncertainty.reset( new TF1("jetPhiFunc","[0]*x",1) );
0561   //jetUncertainty.phiUncertainty->SetParameter(0, jetPhiUncertaintyParameter0_);
0562
0563   //-- phi Functions and values from
0564   // http://indico.cern.ch/getFile.py/access?contribId=9&sessionId=0&resId=0&materialId=slides&confId=46394
0565   jetUncertainty.phiUncertainty =
0566       std::make_unique<TF1>("jetPhiFunc", "x*sqrt(([0]*[0]/(x*x))+([1]*[1]/x)+([2]*[2]))", 3);
0567   jetUncertainty.phiUncertainty->SetParameter(0, jetPhiUncertaintyParameter0_);
0568   jetUncertainty.phiUncertainty->SetParameter(1, jetPhiUncertaintyParameter1_);
0569   jetUncertainty.phiUncertainty->SetParameter(2, jetPhiUncertaintyParameter2_);
0570
0571   //-- Jet corrections are assumed not to have an error --//
0572   /*jetCorrUncertainty.etUncertainty.reset( new TF1("jetCorrEtFunc","[0]*x",1) );
0573   jetCorrUncertainty.etUncertainty->SetParameter(0,0.0);
0574   jetCorrUncertainty.phiUncertainty.reset( new TF1("jetCorrPhiFunc","[0]*x",1) );
0575   jetCorrUncertainty.phiUncertainty->SetParameter(0,0.0*(3.14159/180.));*/
0576
0577   //--- Electron Uncertainty Functions ---------------------------------//
0578   // completely ambiguious values for electron-like jets...
0579   // the egamma group keeps track of these here:
0580   // https://twiki.cern.ch/twiki/bin/view/CMS/EgammaCMSSWVal
0581   // electron resolution in energy is around 3.4%, measured for 10 < pT < 50 at realistic events with pile-up.
0582
0583   eleUncertainty.etUncertainty = std::make_unique<TF1>("eleEtFunc", "[0] * x", 1);
0584   //  eleUncertainty.etUncertainty->SetParameter(0,0.034);
0585   eleUncertainty.etUncertainty->SetParameter(0, eleEtUncertaintyParameter0_);
0586
0587   eleUncertainty.phiUncertainty = std::make_unique<TF1>("elePhiFunc", "[0] * x", 1);
0588   //  eleUncertainty.phiUncertainty->SetParameter(0,1*(3.14159/180.));
0589   eleUncertainty.phiUncertainty->SetParameter(0, elePhiUncertaintyParameter0_);
0590
0591   //--- Muon Uncertainty Functions ------------------------------------//
0592   // and ambiguious values for the muons...
0593
0594   muonUncertainty.etUncertainty = std::make_unique<TF1>("muonEtFunc", "[0] * x", 1);
0595   //  muonUncertainty.etUncertainty->SetParameter(0,0.01);
0596   muonUncertainty.etUncertainty->SetParameter(0, muonEtUncertaintyParameter0_);
0597   muonUncertainty.phiUncertainty = std::make_unique<TF1>("muonPhiFunc", "[0] * x", 1);
0598   //  muonUncertainty.phiUncertainty->SetParameter(0,1*(3.14159/180.));
0599   muonUncertainty.phiUncertainty->SetParameter(0, muonPhiUncertaintyParameter0_);
0600
0601   //-- Muon calo deposites are assumed not to have an error --//
0602   /*muonCorrUncertainty.etUncertainty.reset( new TF1("muonCorrEtFunc","[0] * x",1) );
0603   muonCorrUncertainty.etUncertainty->SetParameter(0,0.0);
0604   muonCorrUncertainty.phiUncertainty.reset( new TF1("muonCorrPhiFunc","[0] * x",1) );
0605   muonCorrUncertainty.phiUncertainty->SetParameter(0,0.0*(3.14159/180.)); */
0606 }
0607
0608 #include "FWCore/Framework/interface/MakerMacros.h"
0609 using namespace pat;
0610 DEFINE_FWK_MODULE(PATMHTProducer);