Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:21:13

0001 // -*- C++ -*-
0002 //
0003 // Package:    L1Trigger/L1TNtuples
0004 // Class:      L1JetRecoTreeProducer
0005 //
0006 /**\class L1JetRecoTreeProducer L1JetRecoTreeProducer.cc L1Trigger/L1TNtuples/src/L1JetRecoTreeProducer.cc
0007 
0008  Description: Produces tree containing reco quantities
0009 
0010 
0011 */
0012 
0013 // system include files
0014 #include <memory>
0015 
0016 // framework
0017 #include "FWCore/Framework/interface/Frameworkfwd.h"
0018 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/Framework/interface/MakerMacros.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/Framework/interface/ESHandle.h"
0025 
0026 // cond formats
0027 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0028 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0029 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0030 
0031 // data formats
0032 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0033 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0034 #include "DataFormats/JetReco/interface/JetID.h"
0035 #include "DataFormats/METReco/interface/PFMETCollection.h"
0036 #include "DataFormats/MuonReco/interface/Muon.h"
0037 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0038 #include "DataFormats/METReco/interface/PFMET.h"
0039 #include "DataFormats/MuonReco/interface/Muon.h"
0040 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0041 #include "DataFormats/METReco/interface/CaloMET.h"
0042 #include "DataFormats/PatCandidates/interface/Jet.h"
0043 
0044 // ROOT output stuff
0045 #include "FWCore/ServiceRegistry/interface/Service.h"
0046 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0047 #include "TH1.h"
0048 #include "TTree.h"
0049 #include "TF1.h"
0050 #include <TVector2.h>
0051 
0052 //local  data formats
0053 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoJetDataFormat.h"
0054 #include "L1Trigger/L1TNtuples/interface/L1AnalysisRecoMetDataFormat.h"
0055 
0056 //
0057 // class declaration
0058 //
0059 
0060 class L1JetRecoTreeProducer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0061 public:
0062   explicit L1JetRecoTreeProducer(const edm::ParameterSet&);
0063   ~L1JetRecoTreeProducer() override;
0064 
0065 private:
0066   void beginJob(void) override;
0067   void analyze(const edm::Event&, const edm::EventSetup&) override;
0068   void endJob() override;
0069 
0070   void doPFJets(edm::Handle<reco::PFJetCollection> pfJets);
0071   void doPFJetCorr(edm::Handle<reco::PFJetCollection> pfJets, edm::Handle<reco::JetCorrector> pfJetCorr);
0072   void doPUPPIJets(edm::Handle<reco::PFJetCollection> puppiJets);
0073   void doCorrPUPPIJets(edm::Handle<std::vector<pat::Jet> > corrPuppiJets);
0074   void doCaloJets(edm::Handle<reco::CaloJetCollection> caloJets);
0075   void doCaloJetCorr(edm::Handle<reco::CaloJetCollection> caloJets, edm::Handle<reco::JetCorrector> caloJetCorr);
0076   void doCaloMet(edm::Handle<reco::CaloMETCollection> caloMet);
0077   void doCaloMetBE(edm::Handle<reco::CaloMETCollection> caloMetBE);
0078 
0079   void doPFMet(edm::Handle<reco::PFMETCollection> pfMet);
0080   void doPFMetNoMu(edm::Handle<reco::PFMETCollection> pfMet, edm::Handle<reco::MuonCollection>);
0081   void doPUPPIMetNoMu(edm::Handle<reco::PFMETCollection> puppiMet, edm::Handle<reco::MuonCollection>);
0082 
0083   void doZPt(edm::Handle<reco::MuonCollection> muons, edm::Handle<std::vector<pat::Jet> > corrPuppiJets);
0084 
0085   bool pfJetID(const reco::PFJet& jet);
0086   bool puppiJetID(const pat::Jet& jet);
0087   bool caloJetID(const reco::CaloJet& jet);
0088 
0089 public:
0090   L1Analysis::L1AnalysisRecoJetDataFormat* jet_data;
0091   L1Analysis::L1AnalysisRecoMetDataFormat* met_data;
0092 
0093 private:
0094   // tree
0095   TTree* tree_;
0096 
0097   // EDM input tags
0098   edm::EDGetTokenT<reco::PFJetCollection> pfJetToken_;
0099   edm::EDGetTokenT<reco::PFJetCollection> puppiJetToken_;
0100   edm::EDGetTokenT<std::vector<pat::Jet> > corrPuppiJetToken_;
0101   edm::EDGetTokenT<reco::CaloJetCollection> caloJetToken_;
0102   edm::EDGetTokenT<edm::ValueMap<reco::JetID> > caloJetIDToken_;
0103   edm::EDGetTokenT<reco::JetCorrector> pfJECToken_;
0104   edm::EDGetTokenT<reco::JetCorrector> caloJECToken_;
0105 
0106   edm::EDGetTokenT<reco::PFMETCollection> pfMetToken_;
0107   edm::EDGetTokenT<reco::PFMETCollection> puppiMetToken_;
0108   edm::EDGetTokenT<reco::CaloMETCollection> caloMetToken_;
0109   edm::EDGetTokenT<reco::CaloMETCollection> caloMetBEToken_;
0110 
0111   edm::EDGetTokenT<reco::MuonCollection> muonToken_;
0112 
0113   // debug stuff
0114   bool pfJetsMissing_;
0115   bool puppiJetsMissing_;
0116   bool corrPuppiJetsMissing_;
0117   double jetptThreshold_;
0118   double jetetaMax_;
0119   unsigned int maxCl_;
0120   unsigned int maxJet_;
0121   unsigned int maxVtx_;
0122   unsigned int maxTrk_;
0123 
0124   bool pfJetCorrMissing_;
0125   bool caloJetCorrMissing_;
0126   bool caloJetsMissing_;
0127   bool caloJetIDMissing_;
0128   bool pfMetMissing_;
0129   bool puppiMetMissing_;
0130   bool caloMetMissing_;
0131   bool caloMetBEMissing_;
0132   bool muonsMissing_;
0133 };
0134 
0135 L1JetRecoTreeProducer::L1JetRecoTreeProducer(const edm::ParameterSet& iConfig)
0136     : pfJetsMissing_(false),
0137       puppiJetsMissing_(false),
0138       corrPuppiJetsMissing_(false),
0139       pfJetCorrMissing_(false),
0140       caloJetCorrMissing_(false),
0141       caloJetsMissing_(false),
0142       caloJetIDMissing_(false),
0143       pfMetMissing_(false),
0144       puppiMetMissing_(false),
0145       caloMetMissing_(false),
0146       caloMetBEMissing_(false),
0147       muonsMissing_(false) {
0148   caloJetToken_ =
0149       consumes<reco::CaloJetCollection>(iConfig.getUntrackedParameter("caloJetToken", edm::InputTag("ak4CaloJets")));
0150   pfJetToken_ =
0151       consumes<reco::PFJetCollection>(iConfig.getUntrackedParameter("pfJetToken", edm::InputTag("ak4PFJetsCHS")));
0152   puppiJetToken_ = consumes<reco::PFJetCollection>(iConfig.getParameter<edm::InputTag>("puppiJetToken"));
0153   corrPuppiJetToken_ = consumes<std::vector<pat::Jet> >(
0154       iConfig.getUntrackedParameter("corrPuppiJetToken", edm::InputTag("patJetsCorrectedPuppiJets")));
0155   caloJetIDToken_ =
0156       consumes<edm::ValueMap<reco::JetID> >(iConfig.getUntrackedParameter("caloJetIDToken", edm::InputTag("ak4JetID")));
0157   pfJECToken_ = consumes<reco::JetCorrector>(
0158       iConfig.getUntrackedParameter<edm::InputTag>("pfJECToken", edm::InputTag("ak4PFCHSL1FastL2L3ResidualCorrector")));
0159   caloJECToken_ = consumes<reco::JetCorrector>(iConfig.getUntrackedParameter<edm::InputTag>(
0160       "caloJECToken", edm::InputTag("ak4CaloL1FastL2L3ResidualCorrector")));
0161 
0162   pfMetToken_ = consumes<reco::PFMETCollection>(iConfig.getUntrackedParameter("pfMetToken", edm::InputTag("pfMetT1")));
0163   puppiMetToken_ =
0164       consumes<reco::PFMETCollection>(iConfig.getUntrackedParameter("puppiMetToken", edm::InputTag("pfMetPuppi")));
0165   caloMetToken_ =
0166       consumes<reco::CaloMETCollection>(iConfig.getUntrackedParameter("caloMetToken", edm::InputTag("caloMet")));
0167   caloMetBEToken_ =
0168       consumes<reco::CaloMETCollection>(iConfig.getUntrackedParameter("caloMetBEToken", edm::InputTag("caloMetBE")));
0169 
0170   muonToken_ = consumes<reco::MuonCollection>(iConfig.getUntrackedParameter("muonToken", edm::InputTag("muons")));
0171 
0172   usesResource(TFileService::kSharedResource);
0173 
0174   jetptThreshold_ = iConfig.getParameter<double>("jetptThreshold");
0175   jetetaMax_ = iConfig.getParameter<double>("jetetaMax");
0176   maxJet_ = iConfig.getParameter<unsigned int>("maxJet");
0177 
0178   jet_data = new L1Analysis::L1AnalysisRecoJetDataFormat();
0179   met_data = new L1Analysis::L1AnalysisRecoMetDataFormat();
0180 
0181   // set up output
0182   edm::Service<TFileService> fs_;
0183   tree_ = fs_->make<TTree>("JetRecoTree", "JetRecoTree");
0184   tree_->Branch("Jet", "L1Analysis::L1AnalysisRecoJetDataFormat", &jet_data, 32000, 3);
0185   tree_->Branch("Sums", "L1Analysis::L1AnalysisRecoMetDataFormat", &met_data, 32000, 3);
0186 }
0187 
0188 L1JetRecoTreeProducer::~L1JetRecoTreeProducer() {
0189   // do anything here that needs to be done at desctruction time
0190   // (e.g. close files, deallocate resources etc.)
0191 }
0192 
0193 //
0194 // member functions
0195 //
0196 
0197 // ------------ method called to for each event  ------------
0198 void L1JetRecoTreeProducer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0199   jet_data->Reset();
0200   met_data->Reset();
0201 
0202   // get jets
0203   edm::Handle<reco::PFJetCollection> pfJets;
0204   iEvent.getByToken(pfJetToken_, pfJets);
0205 
0206   // get puppi jets
0207   edm::Handle<reco::PFJetCollection> puppiJets;
0208   iEvent.getByToken(puppiJetToken_, puppiJets);
0209 
0210   // get corrected puppi jets
0211   edm::Handle<std::vector<pat::Jet> > corrPuppiJets;
0212   iEvent.getByToken(corrPuppiJetToken_, corrPuppiJets);
0213 
0214   // get calo jets
0215   edm::Handle<reco::CaloJetCollection> caloJets;
0216   iEvent.getByToken(caloJetToken_, caloJets);
0217 
0218   //get sums
0219   edm::Handle<reco::PFMETCollection> pfMet;
0220   iEvent.getByToken(pfMetToken_, pfMet);
0221 
0222   //get sums
0223   edm::Handle<reco::PFMETCollection> puppiMet;
0224   iEvent.getByToken(puppiMetToken_, puppiMet);
0225 
0226   // get jet ID
0227   edm::Handle<edm::ValueMap<reco::JetID> > jetsID;
0228   iEvent.getByToken(caloJetIDToken_, jetsID);
0229 
0230   edm::Handle<reco::JetCorrector> pfJetCorr;
0231   iEvent.getByToken(pfJECToken_, pfJetCorr);
0232 
0233   edm::Handle<reco::JetCorrector> caloJetCorr;
0234   iEvent.getByToken(caloJECToken_, caloJetCorr);
0235 
0236   edm::Handle<reco::CaloMETCollection> caloMet;
0237   iEvent.getByToken(caloMetToken_, caloMet);
0238 
0239   edm::Handle<reco::CaloMETCollection> caloMetBE;
0240   iEvent.getByToken(caloMetBEToken_, caloMetBE);
0241 
0242   // get muons
0243   edm::Handle<reco::MuonCollection> muons;
0244   iEvent.getByToken(muonToken_, muons);
0245 
0246   if (pfJets.isValid()) {
0247     jet_data->nJets = 0;
0248 
0249     doPFJets(pfJets);
0250 
0251   } else {
0252     if (!pfJetsMissing_) {
0253       edm::LogWarning("MissingProduct") << "PFJets not found.  Branch will not be filled" << std::endl;
0254     }
0255     pfJetsMissing_ = true;
0256   }
0257 
0258   if (pfJetCorr.isValid()) {
0259     doPFJetCorr(pfJets, pfJetCorr);
0260 
0261   } else {
0262     if (!pfJetCorrMissing_) {
0263       edm::LogWarning("MissingProduct") << "PF Jet Corrector not found.  Branch will not be filled" << std::endl;
0264     }
0265     pfJetCorrMissing_ = true;
0266   }
0267 
0268   if (puppiJets.isValid()) {
0269     jet_data->puppi_nUncorrJets = 0;
0270 
0271     doPUPPIJets(puppiJets);
0272 
0273   } else {
0274     if (!puppiJetsMissing_) {
0275       edm::LogWarning("MissingProduct") << "PUPPIJets not found.  Branch will not be filled" << std::endl;
0276     }
0277     puppiJetsMissing_ = true;
0278   }
0279 
0280   if (corrPuppiJets.isValid()) {
0281     jet_data->puppi_nJets = 0;
0282 
0283     doCorrPUPPIJets(corrPuppiJets);
0284 
0285   } else {
0286     if (!corrPuppiJetsMissing_) {
0287       edm::LogWarning("MissingProduct") << "Corrected PUPPIJets not found.  Branch will not be filled" << std::endl;
0288     }
0289     corrPuppiJetsMissing_ = true;
0290   }
0291 
0292   if (caloJets.isValid()) {
0293     jet_data->nCaloJets = 0;
0294 
0295     doCaloJets(caloJets);
0296 
0297   } else {
0298     if (!caloJetsMissing_) {
0299       edm::LogWarning("MissingProduct") << "Calo Jets not found.  Branch will not be filled" << std::endl;
0300     }
0301     caloJetsMissing_ = true;
0302   }
0303 
0304   if (caloJetCorr.isValid()) {
0305     doCaloJetCorr(caloJets, caloJetCorr);
0306 
0307   } else {
0308     if (!caloJetCorrMissing_) {
0309       edm::LogWarning("MissingProduct") << "Calo Jet Corrector not found.  Branch will not be filled" << std::endl;
0310     }
0311     caloJetCorrMissing_ = true;
0312   }
0313 
0314   if (!jetsID.isValid()) {
0315     if (!caloJetIDMissing_) {
0316       edm::LogWarning("MissingProduct") << "Calo Jet ID not found.  Branch will not be filled" << std::endl;
0317     }
0318     caloJetIDMissing_ = true;
0319   }
0320 
0321   if (pfMet.isValid()) {
0322     doPFMet(pfMet);
0323 
0324     if (muons.isValid()) {
0325       doPFMetNoMu(pfMet, muons);
0326 
0327     } else {
0328       if (!muonsMissing_) {
0329         edm::LogWarning("MissingProduct") << "Muons not found.  PFMetNoMu branch will not be filled" << std::endl;
0330       }
0331       muonsMissing_ = true;
0332     }
0333   } else {
0334     if (!pfMetMissing_) {
0335       edm::LogWarning("MissingProduct") << "PFMet not found.  Branch will not be filled" << std::endl;
0336     }
0337     pfMetMissing_ = true;
0338   }
0339 
0340   if (puppiMet.isValid()) {
0341     if (muons.isValid()) {
0342       doPUPPIMetNoMu(puppiMet, muons);
0343 
0344     } else {
0345       if (!muonsMissing_) {
0346         edm::LogWarning("MissingProduct") << "Muons not found.  PUPPIMetNoMu branch will not be filled" << std::endl;
0347       }
0348       muonsMissing_ = true;
0349     }
0350   } else {
0351     if (!puppiMetMissing_) {
0352       edm::LogWarning("MissingProduct") << "PUPPIMet not found.  Branch will not be filled" << std::endl;
0353     }
0354     puppiMetMissing_ = true;
0355   }
0356 
0357   if (caloMet.isValid()) {
0358     doCaloMet(caloMet);
0359 
0360   } else {
0361     if (!caloMetMissing_) {
0362       edm::LogWarning("MissingProduct") << "CaloMet not found. Branch will not be filled" << std::endl;
0363     }
0364     caloMetMissing_ = true;
0365   }
0366 
0367   if (caloMetBE.isValid()) {
0368     doCaloMetBE(caloMetBE);
0369 
0370   } else {
0371     if (!caloMetBEMissing_) {
0372       edm::LogWarning("MissingProduct") << "CaloMetBE not found. Branch will not be filled" << std::endl;
0373     }
0374     caloMetBEMissing_ = true;
0375   }
0376 
0377   if (muons.isValid()) {
0378     if (puppiJets.isValid()) {
0379       doZPt(muons, corrPuppiJets);
0380 
0381     } else {
0382       if (!puppiJetsMissing_) {
0383         edm::LogWarning("MissingProduct") << "PUPPIJets not found.  Branch will not be filled" << std::endl;
0384       }
0385       puppiJetsMissing_ = true;
0386     }
0387   } else {
0388     if (!muonsMissing_) {
0389       edm::LogWarning("MissingProduct") << "Muons not found.  ZPt branch will not be filled" << std::endl;
0390     }
0391     muonsMissing_ = true;
0392   }
0393 
0394   tree_->Fill();
0395 }
0396 
0397 void L1JetRecoTreeProducer::doCaloJets(edm::Handle<reco::CaloJetCollection> caloJets) {
0398   for (auto it = caloJets->begin(); it != caloJets->end() && jet_data->nCaloJets < maxJet_; ++it) {
0399     if (!caloJetIDMissing_)
0400       if (!caloJetID(*it))
0401         continue;
0402 
0403     jet_data->caloEt.push_back(it->et());
0404     jet_data->caloEta.push_back(it->eta());
0405     jet_data->caloPhi.push_back(it->phi());
0406     jet_data->caloE.push_back(it->energy());
0407 
0408     jet_data->eEMF.push_back(it->emEnergyFraction());
0409     jet_data->eEmEB.push_back(it->emEnergyInEB());
0410     jet_data->eEmEE.push_back(it->emEnergyInEE());
0411     jet_data->eEmHF.push_back(it->emEnergyInHF());
0412     jet_data->eHadHB.push_back(it->hadEnergyInHB());
0413     jet_data->eHadHE.push_back(it->hadEnergyInHE());
0414     jet_data->eHadHO.push_back(it->hadEnergyInHO());
0415     jet_data->eHadHF.push_back(it->hadEnergyInHF());
0416     jet_data->eMaxEcalTow.push_back(it->maxEInEmTowers());
0417     jet_data->eMaxHcalTow.push_back(it->maxEInHadTowers());
0418     jet_data->towerArea.push_back(it->towersArea());
0419     jet_data->n60.push_back(it->n60());
0420 
0421     jet_data->nCaloJets++;
0422   }
0423 }
0424 
0425 void L1JetRecoTreeProducer::doPFJets(edm::Handle<reco::PFJetCollection> pfJets) {
0426   for (auto it = pfJets->begin(); it != pfJets->end() && jet_data->nJets < maxJet_; ++it) {
0427     if (!pfJetID(*it))
0428       continue;
0429 
0430     jet_data->et.push_back(it->et());
0431     jet_data->eta.push_back(it->eta());
0432     jet_data->phi.push_back(it->phi());
0433     jet_data->e.push_back(it->energy());
0434 
0435     jet_data->chef.push_back(it->chargedHadronEnergyFraction());
0436     jet_data->nhef.push_back(it->neutralHadronEnergyFraction());
0437     jet_data->pef.push_back(it->photonEnergyFraction());
0438     jet_data->eef.push_back(it->electronEnergyFraction());
0439     jet_data->mef.push_back(it->muonEnergyFraction());
0440     jet_data->hfhef.push_back(it->HFHadronEnergyFraction());
0441     jet_data->hfemef.push_back(it->HFEMEnergyFraction());
0442     jet_data->chMult.push_back(it->chargedHadronMultiplicity());
0443     jet_data->nhMult.push_back(it->neutralHadronMultiplicity());
0444     jet_data->phMult.push_back(it->photonMultiplicity());
0445     jet_data->elMult.push_back(it->electronMultiplicity());
0446     jet_data->muMult.push_back(it->muonMultiplicity());
0447     jet_data->hfhMult.push_back(it->HFHadronMultiplicity());
0448     jet_data->hfemMult.push_back(it->HFEMMultiplicity());
0449 
0450     jet_data->cemef.push_back(it->chargedEmEnergyFraction());
0451     jet_data->cmef.push_back(it->chargedMuEnergyFraction());
0452     jet_data->nemef.push_back(it->neutralEmEnergyFraction());
0453     jet_data->cMult.push_back(it->chargedMultiplicity());
0454     jet_data->nMult.push_back(it->neutralMultiplicity());
0455 
0456     jet_data->nJets++;
0457   }
0458 }
0459 
0460 void L1JetRecoTreeProducer::doPFJetCorr(edm::Handle<reco::PFJetCollection> pfJets,
0461                                         edm::Handle<reco::JetCorrector> pfJetCorr) {
0462   float corrFactor = 1.;
0463   unsigned int nJets = 0;
0464 
0465   float mHx = 0;
0466   float mHy = 0;
0467 
0468   met_data->Ht = 0;
0469   met_data->mHt = -999.;
0470   met_data->mHtPhi = -999.;
0471 
0472   for (auto it = pfJets->begin(); it != pfJets->end() && nJets < maxJet_; ++it) {
0473     if (!pfJetID(*it))
0474       continue;
0475 
0476     corrFactor = pfJetCorr.product()->correction(*it);
0477 
0478     jet_data->etCorr.push_back(it->et() * corrFactor);
0479     jet_data->corrFactor.push_back(corrFactor);
0480 
0481     nJets++;
0482 
0483     if (it->pt() * corrFactor > jetptThreshold_ && std::abs(it->eta()) < jetetaMax_) {
0484       mHx += -1. * it->px() * corrFactor;
0485       mHy += -1. * it->py() * corrFactor;
0486       met_data->Ht += it->pt() * corrFactor;
0487     }
0488   }
0489 
0490   TVector2 tv2 = TVector2(mHx, mHy);
0491   met_data->mHt = tv2.Mod();
0492   met_data->mHtPhi = tv2.Phi();
0493 }
0494 
0495 void L1JetRecoTreeProducer::doPUPPIJets(edm::Handle<reco::PFJetCollection> puppiJets) {
0496   for (auto it = puppiJets->begin(); it != puppiJets->end() && jet_data->puppi_nUncorrJets < maxJet_; ++it) {
0497     if (!puppiJetID(*it))
0498       continue;
0499     jet_data->puppi_et.push_back(it->et());
0500     jet_data->puppi_nUncorrJets++;
0501   }
0502 }
0503 
0504 void L1JetRecoTreeProducer::doCorrPUPPIJets(edm::Handle<std::vector<pat::Jet> > corrPuppiJets) {
0505   float mHx = 0;
0506   float mHy = 0;
0507 
0508   met_data->puppi_Ht = 0;
0509   met_data->puppi_mHt = -999.;
0510   met_data->puppi_mHtPhi = -999.;
0511 
0512   for (auto it = corrPuppiJets->begin(); it != corrPuppiJets->end() && jet_data->puppi_nJets < maxJet_; ++it) {
0513     if (!puppiJetID(*it))
0514       continue;
0515 
0516     jet_data->puppi_etCorr.push_back(it->et());
0517     jet_data->puppi_eta.push_back(it->eta());
0518     jet_data->puppi_phi.push_back(it->phi());
0519     jet_data->puppi_e.push_back(it->energy());
0520 
0521     jet_data->puppi_chef.push_back(it->chargedHadronEnergyFraction());
0522     jet_data->puppi_nhef.push_back(it->neutralHadronEnergyFraction());
0523     jet_data->puppi_pef.push_back(it->photonEnergyFraction());
0524     jet_data->puppi_eef.push_back(it->electronEnergyFraction());
0525     jet_data->puppi_mef.push_back(it->muonEnergyFraction());
0526     jet_data->puppi_hfhef.push_back(it->HFHadronEnergyFraction());
0527     jet_data->puppi_hfemef.push_back(it->HFEMEnergyFraction());
0528     jet_data->puppi_chMult.push_back(it->chargedHadronMultiplicity());
0529     jet_data->puppi_nhMult.push_back(it->neutralHadronMultiplicity());
0530     jet_data->puppi_phMult.push_back(it->photonMultiplicity());
0531     jet_data->puppi_elMult.push_back(it->electronMultiplicity());
0532     jet_data->puppi_muMult.push_back(it->muonMultiplicity());
0533     jet_data->puppi_hfhMult.push_back(it->HFHadronMultiplicity());
0534     jet_data->puppi_hfemMult.push_back(it->HFEMMultiplicity());
0535 
0536     jet_data->puppi_cemef.push_back(it->chargedEmEnergyFraction());
0537     jet_data->puppi_cmef.push_back(it->chargedMuEnergyFraction());
0538     jet_data->puppi_nemef.push_back(it->neutralEmEnergyFraction());
0539     jet_data->puppi_cMult.push_back(it->chargedMultiplicity());
0540     jet_data->puppi_nMult.push_back(it->neutralMultiplicity());
0541 
0542     jet_data->puppi_nJets++;
0543 
0544     if (it->pt() > jetptThreshold_ && std::abs(it->eta()) < jetetaMax_) {
0545       mHx += -1. * it->px();
0546       mHy += -1. * it->py();
0547       met_data->puppi_Ht += it->pt();
0548     }
0549   }
0550 
0551   TVector2 tv2 = TVector2(mHx, mHy);
0552   met_data->puppi_mHt = tv2.Mod();
0553   met_data->puppi_mHtPhi = tv2.Phi();
0554 }
0555 
0556 void L1JetRecoTreeProducer::doCaloJetCorr(edm::Handle<reco::CaloJetCollection> caloJets,
0557                                           edm::Handle<reco::JetCorrector> caloJetCorr) {
0558   float caloCorrFactor = 1.;
0559   unsigned int nCaloJets = 0;
0560 
0561   met_data->caloHt = 0;
0562 
0563   for (auto it = caloJets->begin(); it != caloJets->end() && nCaloJets < maxJet_; ++it) {
0564     if (!caloJetIDMissing_)
0565       if (!caloJetID(*it))
0566         continue;
0567 
0568     caloCorrFactor = caloJetCorr.product()->correction(*it);
0569 
0570     jet_data->caloEtCorr.push_back(it->et() * caloCorrFactor);
0571     jet_data->caloCorrFactor.push_back(caloCorrFactor);
0572 
0573     nCaloJets++;
0574 
0575     if (it->pt() * caloCorrFactor > jetptThreshold_ && std::abs(it->eta()) < jetetaMax_) {
0576       met_data->caloHt += it->pt() * caloCorrFactor;
0577     }
0578   }
0579 }
0580 
0581 void L1JetRecoTreeProducer::doPFMet(edm::Handle<reco::PFMETCollection> pfMet) {
0582   const reco::PFMETCollection* metCol = pfMet.product();
0583   const reco::PFMET theMet = metCol->front();
0584 
0585   met_data->met = theMet.et();
0586   met_data->metPhi = theMet.phi();
0587   met_data->sumEt = theMet.sumEt();
0588   met_data->metPx = theMet.px();
0589   met_data->metPy = theMet.py();
0590 }
0591 
0592 void L1JetRecoTreeProducer::doPFMetNoMu(edm::Handle<reco::PFMETCollection> pfMet,
0593                                         edm::Handle<reco::MuonCollection> muons) {
0594   const reco::PFMETCollection* metCol = pfMet.product();
0595   const reco::PFMET theMet = metCol->front();
0596   reco::PFMET thePFMetNoMu = metCol->front();
0597 
0598   double pfMetNoMuPx = theMet.px();
0599   double pfMetNoMuPy = theMet.py();
0600 
0601   double muPx(0.), muPy(0.);
0602 
0603   for (auto it = muons->begin(); it != muons->end(); ++it) {
0604     if (it->isPFMuon()) {
0605       muPx += it->px();
0606       muPy += it->py();
0607     }
0608   }
0609 
0610   pfMetNoMuPx += muPx;
0611   pfMetNoMuPy += muPy;
0612 
0613   math::XYZTLorentzVector pfMetNoMuP4(pfMetNoMuPx, pfMetNoMuPy, 0, hypot(pfMetNoMuPx, pfMetNoMuPy));
0614 
0615   thePFMetNoMu.setP4(pfMetNoMuP4);
0616 
0617   met_data->pfMetNoMu = thePFMetNoMu.et();
0618   met_data->pfMetNoMuPhi = thePFMetNoMu.phi();
0619   met_data->pfMetNoMuPx = thePFMetNoMu.px();
0620   met_data->pfMetNoMuPy = thePFMetNoMu.py();
0621 }
0622 
0623 void L1JetRecoTreeProducer::doPUPPIMetNoMu(edm::Handle<reco::PFMETCollection> puppiMet,
0624                                            edm::Handle<reco::MuonCollection> muons) {
0625   const reco::PFMETCollection* metCol = puppiMet.product();
0626   const reco::PFMET theMet = metCol->front();
0627   reco::PFMET thePUPPIMetNoMu = metCol->front();
0628 
0629   double puppiMetNoMuPx = theMet.px();
0630   double puppiMetNoMuPy = theMet.py();
0631 
0632   double muPx(0.), muPy(0.);
0633 
0634   for (auto it = muons->begin(); it != muons->end(); ++it) {
0635     if (it->isPFMuon()) {
0636       muPx += it->px();
0637       muPy += it->py();
0638     }
0639   }
0640 
0641   puppiMetNoMuPx += muPx;
0642   puppiMetNoMuPy += muPy;
0643 
0644   math::XYZTLorentzVector puppiMetNoMuP4(puppiMetNoMuPx, puppiMetNoMuPy, 0, hypot(puppiMetNoMuPx, puppiMetNoMuPy));
0645 
0646   thePUPPIMetNoMu.setP4(puppiMetNoMuP4);
0647 
0648   met_data->puppi_metNoMu = thePUPPIMetNoMu.et();
0649   met_data->puppi_metNoMuPhi = thePUPPIMetNoMu.phi();
0650   met_data->puppi_metNoMuPx = thePUPPIMetNoMu.px();
0651   met_data->puppi_metNoMuPy = thePUPPIMetNoMu.py();
0652 }
0653 
0654 void L1JetRecoTreeProducer::doCaloMet(edm::Handle<reco::CaloMETCollection> caloMet) {
0655   const reco::CaloMETCollection* metCol = caloMet.product();
0656   const reco::CaloMET theMet = metCol->front();
0657 
0658   met_data->caloMet = theMet.et();
0659   met_data->caloMetPhi = theMet.phi();
0660   met_data->caloSumEt = theMet.sumEt();
0661 }
0662 
0663 void L1JetRecoTreeProducer::doCaloMetBE(edm::Handle<reco::CaloMETCollection> caloMetBE) {
0664   const reco::CaloMETCollection* metCol = caloMetBE.product();
0665   const reco::CaloMET theMet = metCol->front();
0666 
0667   met_data->caloMetBE = theMet.et();
0668   met_data->caloMetPhiBE = theMet.phi();
0669   met_data->caloSumEtBE = theMet.sumEt();
0670 }
0671 
0672 void L1JetRecoTreeProducer::doZPt(edm::Handle<reco::MuonCollection> muons,
0673                                   edm::Handle<std::vector<pat::Jet> > corrPuppiJets) {
0674   bool passPuppiJetPtCut = false;
0675 
0676   for (auto it = corrPuppiJets->begin(); it != corrPuppiJets->end(); ++it) {
0677     if (!puppiJetID(*it))
0678       continue;
0679     if (it->muonEnergyFraction() > 0.5 || it->chargedEmEnergyFraction() > 0.5)
0680       continue;
0681     if (it->pt() > 30)
0682       passPuppiJetPtCut = true;
0683   }
0684 
0685   if (!passPuppiJetPtCut || muons->size() < 2) {
0686     met_data->zPt = -999;
0687     return;
0688   }
0689 
0690   reco::Muon muon1;
0691   reco::Muon muon2;
0692 
0693   float zMass = 91.2;
0694   float diMuMass = 0;
0695   float closestDiff = 999.;
0696   bool found2PFMuons = false;
0697 
0698   for (auto it1 = muons->begin(); it1 != muons->end(); ++it1) {
0699     if (!it1->isPFMuon())
0700       continue;
0701     for (auto it2 = std::next(it1); it2 != muons->end(); ++it2) {
0702       if (!it2->isPFMuon())
0703         continue;
0704       if (it1->charge() != (-1 * it2->charge()))
0705         continue;
0706 
0707       found2PFMuons = true;
0708       diMuMass = (it1->p4() + it2->p4()).M();
0709       float diff = abs(diMuMass - zMass);
0710       if (diff < closestDiff) {
0711         closestDiff = diff;
0712         muon1 = *it1;
0713         muon2 = *it2;
0714       }
0715     }
0716   }
0717 
0718   diMuMass = (muon1.p4() + muon2.p4()).M();
0719   if (abs(diMuMass - zMass) > 30 || !found2PFMuons) {
0720     met_data->zPt = -999;
0721     return;
0722   }
0723 
0724   float zPt = (muon1.p4() + muon2.p4()).pt();
0725   met_data->zPt = zPt;
0726 }
0727 
0728 bool L1JetRecoTreeProducer::pfJetID(const reco::PFJet& jet) {
0729   bool tmp = true;
0730   if (std::abs(jet.eta()) <= 2.6) {
0731     tmp &= jet.neutralHadronEnergyFraction() < 0.9;
0732     tmp &= jet.neutralEmEnergyFraction() < 0.9;
0733     tmp &= (jet.chargedMultiplicity() + jet.neutralMultiplicity()) > 1;
0734     tmp &= jet.muonEnergyFraction() < 0.8;
0735     tmp &= jet.chargedHadronEnergyFraction() > 0.01;
0736     tmp &= jet.chargedMultiplicity() > 0;
0737     tmp &= jet.chargedEmEnergyFraction() < 0.8;
0738   }
0739   if (std::abs(jet.eta()) > 2.6 && std::abs(jet.eta()) <= 2.7) {
0740     tmp &= jet.neutralHadronEnergyFraction() < 0.9;
0741     tmp &= jet.neutralEmEnergyFraction() < 0.99;
0742     tmp &= jet.muonEnergyFraction() < 0.8;
0743     tmp &= jet.chargedMultiplicity() > 0;
0744     tmp &= jet.chargedEmEnergyFraction() < 0.8;
0745   }
0746   if (std::abs(jet.eta()) > 2.7 && std::abs(jet.eta()) < 3.0) {
0747     tmp &= jet.neutralEmEnergyFraction() < 0.99;
0748     tmp &= jet.neutralMultiplicity() > 1;
0749   }
0750   if (std::abs(jet.eta()) > 3.0) {
0751     tmp &= jet.neutralHadronEnergyFraction() > 0.2;
0752     tmp &= jet.neutralEmEnergyFraction() < 0.9;
0753     tmp &= jet.neutralMultiplicity() > 10;
0754   }
0755 
0756   // our custom selection
0757   //tmp &= jet.muonMultiplicity() == 0;
0758   //tmp &= jet.electronMultiplicity() == 0;
0759 
0760   return tmp;
0761 }
0762 
0763 //https://twiki.cern.ch/twiki/bin/view/CMS/JetID13p6TeV
0764 bool L1JetRecoTreeProducer::puppiJetID(const pat::Jet& jet) {
0765   bool tmp = true;
0766   if (std::abs(jet.eta()) <= 2.6) {
0767     tmp &= jet.neutralHadronEnergyFraction() < 0.9;
0768     tmp &= jet.neutralEmEnergyFraction() < 0.9;
0769     tmp &= (jet.chargedMultiplicity() + jet.neutralMultiplicity()) > 1;
0770     tmp &= jet.muonEnergyFraction() < 0.8;
0771     tmp &= jet.chargedHadronEnergyFraction() > 0.01;
0772     tmp &= jet.chargedMultiplicity() > 0;
0773     tmp &= jet.chargedEmEnergyFraction() < 0.8;
0774   }
0775   if (std::abs(jet.eta()) > 2.6 && std::abs(jet.eta()) <= 2.7) {
0776     tmp &= jet.neutralHadronEnergyFraction() < 0.9;
0777     tmp &= jet.neutralEmEnergyFraction() < 0.99;
0778     tmp &= jet.muonEnergyFraction() < 0.8;
0779     tmp &= jet.chargedEmEnergyFraction() < 0.8;
0780   }
0781   if (std::abs(jet.eta()) > 2.7 && std::abs(jet.eta()) <= 3.0) {
0782     tmp &= jet.neutralHadronEnergyFraction() < 0.9999;
0783   }
0784   if (std::abs(jet.eta()) > 3.0) {
0785     tmp &= jet.neutralEmEnergyFraction() < 0.9;
0786     tmp &= jet.neutralMultiplicity() > 2;
0787   }
0788 
0789   return tmp;
0790 }
0791 
0792 bool L1JetRecoTreeProducer::caloJetID(const reco::CaloJet& jet) {
0793   bool tmp = true;
0794 
0795   return tmp;
0796 }
0797 
0798 // ------------ method called once each job just before starting event loop  ------------
0799 void L1JetRecoTreeProducer::beginJob(void) {}
0800 
0801 // ------------ method called once each job just after ending the event loop  ------------
0802 void L1JetRecoTreeProducer::endJob() {}
0803 
0804 //define this as a plug-in
0805 DEFINE_FWK_MODULE(L1JetRecoTreeProducer);