Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-02-16 06:12:09 

0001 //
0002 // DiJetAnalyzer.cc
0003 //
0004 // description: Create an ntuple necessary for dijet balance calibration for the HCAL
0005 //
0006 
0007 // system include files
0008 #include <memory>
0009 #include <string>
0010 #include <vector>
0011 #include <set>
0012 #include <map>
0013 
0014 #include "TTree.h"
0015 #include "TFile.h"
0016 #include "TH1D.h"
0017 #include "TH2D.h"
0018 #include "TClonesArray.h"
0019 
0020 // user include files
0021 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0022 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0023 
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/EventSetup.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/ServiceRegistry/interface/Service.h"
0031 #include "FWCore/Utilities/interface/EDMException.h"
0032 
0033 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0034 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0035 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0036 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
0037 #include "DataFormats/HcalRecHit/interface/HFRecHit.h"
0038 #include "DataFormats/HcalRecHit/interface/HORecHit.h"
0039 #include "DataFormats/ParticleFlowReco/interface/PFBlockFwd.h"
0040 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0041 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0042 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0043 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
0044 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0045 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0046 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0047 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0048 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0049 
0050 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0051 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0052 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0053 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0054 
0055 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
0056 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0057 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0058 
0059 //
0060 // class declarations
0061 //
0062 
0063 class JetCorretPair : protected std::pair<const reco::PFJet*, double> {
0064 public:
0065   JetCorretPair() {
0066     first = 0;
0067     second = 1.0;
0068   }
0069   JetCorretPair(const reco::PFJet* j, double s) {
0070     first = j;
0071     second = s;
0072   }
0073   ~JetCorretPair() {}
0074 
0075   inline const reco::PFJet* jet(void) const { return first; }
0076   inline void jet(const reco::PFJet* j) {
0077     first = j;
0078     return;
0079   }
0080   inline double scale(void) const { return second; }
0081   inline void scale(double d) {
0082     second = d;
0083     return;
0084   }
0085 
0086 private:
0087 };
0088 
0089 class DiJetAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0090 public:
0091   explicit DiJetAnalyzer(const edm::ParameterSet&);
0092   ~DiJetAnalyzer() override = default;
0093 
0094 private:
0095   void beginJob() override;  //(const edm::EventSetup&);
0096   void analyze(const edm::Event&, const edm::EventSetup&) override;
0097   void endJob() override;
0098 
0099   // parameters
0100   const std::string pfJetCollName_;     // label for the PF jet collection
0101   const std::string pfJetCorrName_;     // label for the PF jet correction service
0102   const std::string hbheRecHitName_;    // label for HBHERecHits collection
0103   const std::string hfRecHitName_;      // label for HFRecHit collection
0104   const std::string hoRecHitName_;      // label for HORecHit collection
0105   const std::string pvCollName_;        // label for primary vertex collection
0106   const std::string rootHistFilename_;  // name of the histogram file
0107   const double maxDeltaEta_;            // maximum delta-|Eta| between Jets
0108   const double minTagJetEta_;           // minimum |eta| of the tag jet
0109   const double maxTagJetEta_;           // maximum |eta| of the tag jet
0110   const double minSumJetEt_;            // minimum Sum of the tag and probe jet Et
0111   const double minJetEt_;               // minimum Jet Et
0112   const double maxThirdJetEt_;          // maximum 3rd jet Et
0113   const bool debug_;                    // print debug statements
0114 
0115   const edm::EDGetTokenT<reco::PFJetCollection> tok_PFJet_;
0116   const edm::EDGetTokenT<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>> tok_HBHE_;
0117   const edm::EDGetTokenT<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>> tok_HF_;
0118   const edm::EDGetTokenT<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>> tok_HO_;
0119   const edm::EDGetTokenT<reco::VertexCollection> tok_Vertex_;
0120 
0121   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0122 
0123   // root file/histograms
0124   TH1D* h_PassSelPF_;
0125   TTree* tree_;
0126 
0127   float tpfjet_pt_, tpfjet_p_, tpfjet_E_, tpfjet_eta_, tpfjet_phi_, tpfjet_EMfrac_, tpfjet_hadEcalEfrac_, tpfjet_scale_,
0128       tpfjet_area_;
0129   int tpfjet_jetID_;
0130   float tpfjet_EBE_, tpfjet_EEE_, tpfjet_HBE_, tpfjet_HEE_, tpfjet_HFE_;
0131   float tpfjet_unkown_E_, tpfjet_unkown_px_, tpfjet_unkown_py_, tpfjet_unkown_pz_, tpfjet_unkown_EcalE_;
0132   float tpfjet_electron_E_, tpfjet_electron_px_, tpfjet_electron_py_, tpfjet_electron_pz_, tpfjet_electron_EcalE_;
0133   float tpfjet_muon_E_, tpfjet_muon_px_, tpfjet_muon_py_, tpfjet_muon_pz_, tpfjet_muon_EcalE_;
0134   float tpfjet_photon_E_, tpfjet_photon_px_, tpfjet_photon_py_, tpfjet_photon_pz_, tpfjet_photon_EcalE_;
0135   int tpfjet_unkown_n_, tpfjet_electron_n_, tpfjet_muon_n_, tpfjet_photon_n_;
0136   int tpfjet_had_n_;
0137   std::vector<float> tpfjet_had_E_, tpfjet_had_px_, tpfjet_had_py_, tpfjet_had_pz_, tpfjet_had_EcalE_,
0138       tpfjet_had_rawHcalE_, tpfjet_had_emf_;
0139   std::vector<int> tpfjet_had_id_, tpfjet_had_candtrackind_, tpfjet_had_ntwrs_;
0140   int tpfjet_ntwrs_;
0141   std::vector<int> tpfjet_twr_ieta_, tpfjet_twr_iphi_, tpfjet_twr_depth_, tpfjet_twr_subdet_, tpfjet_twr_candtrackind_,
0142       tpfjet_twr_hadind_, tpfjet_twr_elmttype_, tpfjet_twr_clusterind_;
0143   std::vector<float> tpfjet_twr_hade_, tpfjet_twr_frac_, tpfjet_twr_dR_;
0144   int tpfjet_cluster_n_;
0145   std::vector<float> tpfjet_cluster_eta_, tpfjet_cluster_phi_, tpfjet_cluster_dR_;
0146   int tpfjet_ncandtracks_;
0147   std::vector<float> tpfjet_candtrack_px_, tpfjet_candtrack_py_, tpfjet_candtrack_pz_, tpfjet_candtrack_EcalE_;
0148   float ppfjet_pt_, ppfjet_p_, ppfjet_E_, ppfjet_eta_, ppfjet_phi_, ppfjet_EMfrac_, ppfjet_hadEcalEfrac_, ppfjet_scale_,
0149       ppfjet_area_;
0150   int ppfjet_jetID_;
0151   float ppfjet_EBE_, ppfjet_EEE_, ppfjet_HBE_, ppfjet_HEE_, ppfjet_HFE_;
0152   float ppfjet_unkown_E_, ppfjet_unkown_px_, ppfjet_unkown_py_, ppfjet_unkown_pz_, ppfjet_unkown_EcalE_;
0153   float ppfjet_electron_E_, ppfjet_electron_px_, ppfjet_electron_py_, ppfjet_electron_pz_, ppfjet_electron_EcalE_;
0154   float ppfjet_muon_E_, ppfjet_muon_px_, ppfjet_muon_py_, ppfjet_muon_pz_, ppfjet_muon_EcalE_;
0155   float ppfjet_photon_E_, ppfjet_photon_px_, ppfjet_photon_py_, ppfjet_photon_pz_, ppfjet_photon_EcalE_;
0156   int ppfjet_unkown_n_, ppfjet_electron_n_, ppfjet_muon_n_, ppfjet_photon_n_;
0157   int ppfjet_had_n_;
0158   std::vector<float> ppfjet_had_E_, ppfjet_had_px_, ppfjet_had_py_, ppfjet_had_pz_, ppfjet_had_EcalE_,
0159       ppfjet_had_rawHcalE_, ppfjet_had_emf_;
0160   std::vector<int> ppfjet_had_id_, ppfjet_had_candtrackind_, ppfjet_had_ntwrs_;
0161   int ppfjet_ntwrs_;
0162   std::vector<int> ppfjet_twr_ieta_, ppfjet_twr_iphi_, ppfjet_twr_depth_, ppfjet_twr_subdet_, ppfjet_twr_candtrackind_,
0163       ppfjet_twr_hadind_, ppfjet_twr_elmttype_, ppfjet_twr_clusterind_;
0164   std::vector<float> ppfjet_twr_hade_, ppfjet_twr_frac_, ppfjet_twr_dR_;
0165   int ppfjet_cluster_n_;
0166   std::vector<float> ppfjet_cluster_eta_, ppfjet_cluster_phi_, ppfjet_cluster_dR_;
0167   int ppfjet_ncandtracks_;
0168   std::vector<float> ppfjet_candtrack_px_, ppfjet_candtrack_py_, ppfjet_candtrack_pz_, ppfjet_candtrack_EcalE_;
0169   float pf_dijet_deta_, pf_dijet_dphi_, pf_dijet_balance_;
0170   float pf_thirdjet_px_, pf_thirdjet_py_, pf_realthirdjet_px_, pf_realthirdjet_py_, pf_realthirdjet_scale_;
0171   int pf_Run_, pf_Lumi_, pf_Event_;
0172   int pf_NPV_;
0173 
0174   // helper functions
0175   double deltaR(const reco::Jet* j1, const reco::Jet* j2);
0176   double deltaR(const double eta1, const double phi1, const double eta2, const double phi2);
0177   int getEtaPhi(const DetId id);
0178   int getEtaPhi(const HcalDetId id);
0179 
0180   struct JetCorretPairComp {
0181     inline bool operator()(const JetCorretPair& a, const JetCorretPair& b) const {
0182       return (a.jet()->pt() * a.scale()) > (b.jet()->pt() * b.scale());
0183     }
0184   };
0185 };
0186 
0187 DiJetAnalyzer::DiJetAnalyzer(const edm::ParameterSet& iConfig)
0188     : pfJetCollName_(iConfig.getParameter<std::string>("pfJetCollName")),
0189       pfJetCorrName_(iConfig.getParameter<std::string>("pfJetCorrName")),
0190       hbheRecHitName_(iConfig.getParameter<std::string>("hbheRecHitName")),
0191       hfRecHitName_(iConfig.getParameter<std::string>("hfRecHitName")),
0192       hoRecHitName_(iConfig.getParameter<std::string>("hoRecHitName")),
0193       pvCollName_(iConfig.getParameter<std::string>("pvCollName")),
0194       rootHistFilename_(iConfig.getParameter<std::string>("rootHistFilename")),
0195       maxDeltaEta_(iConfig.getParameter<double>("maxDeltaEta")),
0196       minTagJetEta_(iConfig.getParameter<double>("minTagJetEta")),
0197       maxTagJetEta_(iConfig.getParameter<double>("maxTagJetEta")),
0198       minSumJetEt_(iConfig.getParameter<double>("minSumJetEt")),
0199       minJetEt_(iConfig.getParameter<double>("minJetEt")),
0200       maxThirdJetEt_(iConfig.getParameter<double>("maxThirdJetEt")),
0201       debug_(iConfig.getUntrackedParameter<bool>("debug", false)),
0202       tok_PFJet_(consumes<reco::PFJetCollection>(pfJetCollName_)),
0203       tok_HBHE_(consumes<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>>(hbheRecHitName_)),
0204       tok_HF_(consumes<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>>(hfRecHitName_)),
0205       tok_HO_(consumes<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>>(hoRecHitName_)),
0206       tok_Vertex_(consumes<reco::VertexCollection>(pvCollName_)),
0207       tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()) {
0208   usesResource(TFileService::kSharedResource);
0209 }
0210 
0211 //
0212 // member functions
0213 //
0214 
0215 // ------------ method called to for each event  ------------
0216 void DiJetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& evSetup) {
0217   pf_Run_ = iEvent.id().run();
0218   pf_Lumi_ = iEvent.id().luminosityBlock();
0219   pf_Event_ = iEvent.id().event();
0220 
0221   // Get PFJets
0222   const edm::Handle<reco::PFJetCollection> pfjets = iEvent.getHandle(tok_PFJet_);
0223   if (!pfjets.isValid()) {
0224     throw edm::Exception(edm::errors::ProductNotFound)
0225         << " could not find PFJetCollection named " << pfJetCollName_ << ".\n";
0226     return;
0227   }
0228 
0229   // Get RecHits in HB and HE
0230   const edm::Handle<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>> hbhereco =
0231       iEvent.getHandle(tok_HBHE_);
0232   if (!hbhereco.isValid()) {
0233     throw edm::Exception(edm::errors::ProductNotFound)
0234         << " could not find HBHERecHit named " << hbheRecHitName_ << ".\n";
0235     return;
0236   }
0237 
0238   // Get RecHits in HF
0239   const edm::Handle<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>> hfreco =
0240       iEvent.getHandle(tok_HF_);
0241   if (!hfreco.isValid()) {
0242     throw edm::Exception(edm::errors::ProductNotFound) << " could not find HFRecHit named " << hfRecHitName_ << ".\n";
0243     return;
0244   }
0245 
0246   // Get RecHits in HO
0247   const edm::Handle<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>> horeco =
0248       iEvent.getHandle(tok_HO_);
0249   if (!horeco.isValid()) {
0250     throw edm::Exception(edm::errors::ProductNotFound) << " could not find HORecHit named " << hoRecHitName_ << ".\n";
0251     return;
0252   }
0253 
0254   // Get geometry
0255   const CaloGeometry* geo = &evSetup.getData(tok_geom_);
0256   const HcalGeometry* HBGeom = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, 1));
0257   const HcalGeometry* HEGeom = static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, 2));
0258   const CaloSubdetectorGeometry* HOGeom = geo->getSubdetectorGeometry(DetId::Hcal, 3);
0259   const CaloSubdetectorGeometry* HFGeom = geo->getSubdetectorGeometry(DetId::Hcal, 4);
0260 
0261   int HBHE_n = 0;
0262   for (edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>::const_iterator ith = hbhereco->begin();
0263        ith != hbhereco->end();
0264        ++ith) {
0265     HBHE_n++;
0266   }
0267   int HF_n = 0;
0268   for (edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith = hfreco->begin();
0269        ith != hfreco->end();
0270        ++ith) {
0271     HF_n++;
0272   }
0273   int HO_n = 0;
0274   for (edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>::const_iterator ith = horeco->begin();
0275        ith != horeco->end();
0276        ++ith) {
0277     HO_n++;
0278   }
0279 
0280   // Get primary vertices
0281   const edm::Handle<std::vector<reco::Vertex>> pv = iEvent.getHandle(tok_Vertex_);
0282   if (!pv.isValid()) {
0283     throw edm::Exception(edm::errors::ProductNotFound) << " could not find Vertex named " << pvCollName_ << ".\n";
0284     return;
0285   }
0286   pf_NPV_ = 0;
0287   for (std::vector<reco::Vertex>::const_iterator it = pv->begin(); it != pv->end(); ++it) {
0288     if (!it->isFake() && it->ndof() > 4)
0289       ++pf_NPV_;
0290   }
0291 
0292   // Get jet corrections
0293   const JetCorrector* correctorPF = JetCorrector::getJetCorrector(pfJetCorrName_, evSetup);
0294 
0295   //////////////////////////////
0296   // Event Selection
0297   //////////////////////////////
0298 
0299   // determine which cut results in failure
0300   int passSelPF = 0;
0301 
0302   // sort jets by corrected et
0303   std::set<JetCorretPair, JetCorretPairComp> pfjetcorretpairset;
0304   for (reco::PFJetCollection::const_iterator it = pfjets->begin(); it != pfjets->end(); ++it) {
0305     const reco::PFJet* jet = &(*it);
0306     double jec = correctorPF->correction(*it, iEvent, evSetup);
0307     pfjetcorretpairset.insert(JetCorretPair(jet, jec));
0308   }
0309 
0310   JetCorretPair pf_tag, pf_probe;
0311   pf_thirdjet_px_ = pf_thirdjet_py_ = 0.0;
0312   pf_realthirdjet_px_ = pf_realthirdjet_py_ = 0.0;
0313   pf_realthirdjet_px_ = 1;
0314   int cntr = 0;
0315   for (std::set<JetCorretPair, JetCorretPairComp>::const_iterator it = pfjetcorretpairset.begin();
0316        it != pfjetcorretpairset.end();
0317        ++it) {
0318     JetCorretPair jet = (*it);
0319     ++cntr;
0320     if (cntr == 1)
0321       pf_tag = jet;
0322     else if (cntr == 2)
0323       pf_probe = jet;
0324     else {
0325       pf_thirdjet_px_ += jet.scale() * jet.jet()->px();
0326       pf_thirdjet_py_ += jet.scale() * jet.jet()->py();
0327       if (cntr == 3) {
0328         pf_realthirdjet_px_ = jet.jet()->px();
0329         pf_realthirdjet_py_ = jet.jet()->py();
0330         pf_realthirdjet_scale_ = jet.scale();
0331       }
0332     }
0333   }
0334 
0335   if (pf_tag.jet() && pf_probe.jet()) {
0336     // require that the first two jets are above some minimum,
0337     // and the rest are below some maximum
0338     if ((pf_tag.jet()->et() + pf_probe.jet()->et()) < minSumJetEt_)
0339       passSelPF |= 0x1;
0340     if (pf_tag.jet()->et() < minJetEt_ || pf_probe.jet()->et() < minJetEt_)
0341       passSelPF |= 0x2;
0342     if (sqrt(pf_thirdjet_px_ * pf_thirdjet_px_ + pf_thirdjet_py_ * pf_thirdjet_py_) > maxThirdJetEt_)
0343       passSelPF |= 0x4;
0344 
0345     // force the tag jet to have the smaller |eta|
0346     if (std::fabs(pf_tag.jet()->eta()) > std::fabs(pf_probe.jet()->eta())) {
0347       JetCorretPair temp = pf_tag;
0348       pf_tag = pf_probe;
0349       pf_probe = temp;
0350     }
0351 
0352     // eta cuts
0353     double dAbsEta = std::fabs(std::fabs(pf_tag.jet()->eta()) - std::fabs(pf_probe.jet()->eta()));
0354     if (dAbsEta > maxDeltaEta_)
0355       passSelPF |= 0x8;
0356     if (fabs(pf_tag.jet()->eta()) < minTagJetEta_)
0357       passSelPF |= 0x10;
0358     if (fabs(pf_tag.jet()->eta()) > maxTagJetEta_)
0359       passSelPF |= 0x10;
0360   } else {
0361     passSelPF = 0x40;
0362   }
0363 
0364   h_PassSelPF_->Fill(passSelPF);
0365   if (passSelPF)
0366     return;
0367   // dump
0368   if (debug_) {
0369     edm::LogVerbatim("HcalCalib") << "Run: " << iEvent.id().run() << "; Event: " << iEvent.id().event();
0370     for (reco::PFJetCollection::const_iterator it = pfjets->begin(); it != pfjets->end(); ++it) {
0371       const reco::PFJet* jet = &(*it);
0372       edm::LogVerbatim("HcalCalib") << "istag=" << (jet == pf_tag.jet()) << "; isprobe=" << (jet == pf_probe.jet())
0373                                     << "; et=" << jet->et() << "; eta=" << jet->eta();
0374     }
0375   }
0376 
0377   // Reset particle variables
0378   tpfjet_unkown_E_ = tpfjet_unkown_px_ = tpfjet_unkown_py_ = tpfjet_unkown_pz_ = tpfjet_unkown_EcalE_ = 0.0;
0379   tpfjet_electron_E_ = tpfjet_electron_px_ = tpfjet_electron_py_ = tpfjet_electron_pz_ = tpfjet_electron_EcalE_ = 0.0;
0380   tpfjet_muon_E_ = tpfjet_muon_px_ = tpfjet_muon_py_ = tpfjet_muon_pz_ = tpfjet_muon_EcalE_ = 0.0;
0381   tpfjet_photon_E_ = tpfjet_photon_px_ = tpfjet_photon_py_ = tpfjet_photon_pz_ = tpfjet_photon_EcalE_ = 0.0;
0382   tpfjet_unkown_n_ = tpfjet_electron_n_ = tpfjet_muon_n_ = tpfjet_photon_n_ = 0;
0383   tpfjet_had_n_ = 0;
0384   tpfjet_cluster_n_ = 0;
0385   ppfjet_unkown_E_ = ppfjet_unkown_px_ = ppfjet_unkown_py_ = ppfjet_unkown_pz_ = ppfjet_unkown_EcalE_ = 0.0;
0386   ppfjet_electron_E_ = ppfjet_electron_px_ = ppfjet_electron_py_ = ppfjet_electron_pz_ = ppfjet_electron_EcalE_ = 0.0;
0387   ppfjet_muon_E_ = ppfjet_muon_px_ = ppfjet_muon_py_ = ppfjet_muon_pz_ = ppfjet_muon_EcalE_ = 0.0;
0388   ppfjet_photon_E_ = ppfjet_photon_px_ = ppfjet_photon_py_ = ppfjet_photon_pz_ = ppfjet_photon_EcalE_ = 0.0;
0389   ppfjet_unkown_n_ = ppfjet_electron_n_ = ppfjet_muon_n_ = ppfjet_photon_n_ = 0;
0390   ppfjet_had_n_ = 0;
0391   ppfjet_cluster_n_ = 0;
0392 
0393   tpfjet_had_E_.clear();
0394   tpfjet_had_px_.clear();
0395   tpfjet_had_py_.clear();
0396   tpfjet_had_pz_.clear();
0397   tpfjet_had_EcalE_.clear();
0398   tpfjet_had_rawHcalE_.clear();
0399   tpfjet_had_emf_.clear();
0400   tpfjet_had_id_.clear();
0401   tpfjet_had_candtrackind_.clear();
0402   tpfjet_had_ntwrs_.clear();
0403   tpfjet_twr_ieta_.clear();
0404   tpfjet_twr_iphi_.clear();
0405   tpfjet_twr_depth_.clear();
0406   tpfjet_twr_subdet_.clear();
0407   tpfjet_twr_candtrackind_.clear();
0408   tpfjet_twr_hadind_.clear();
0409   tpfjet_twr_elmttype_.clear();
0410   tpfjet_twr_hade_.clear();
0411   tpfjet_twr_frac_.clear();
0412   tpfjet_twr_dR_.clear();
0413   tpfjet_twr_clusterind_.clear();
0414   tpfjet_cluster_eta_.clear();
0415   tpfjet_cluster_phi_.clear();
0416   tpfjet_cluster_dR_.clear();
0417   tpfjet_candtrack_px_.clear();
0418   tpfjet_candtrack_py_.clear();
0419   tpfjet_candtrack_pz_.clear();
0420   tpfjet_candtrack_EcalE_.clear();
0421   ppfjet_had_E_.clear();
0422   ppfjet_had_px_.clear();
0423   ppfjet_had_py_.clear();
0424   ppfjet_had_pz_.clear();
0425   ppfjet_had_EcalE_.clear();
0426   ppfjet_had_rawHcalE_.clear();
0427   ppfjet_had_emf_.clear();
0428   ppfjet_had_id_.clear();
0429   ppfjet_had_candtrackind_.clear();
0430   ppfjet_had_ntwrs_.clear();
0431   ppfjet_twr_ieta_.clear();
0432   ppfjet_twr_iphi_.clear();
0433   ppfjet_twr_depth_.clear();
0434   ppfjet_twr_subdet_.clear();
0435   ppfjet_twr_candtrackind_.clear();
0436   ppfjet_twr_hadind_.clear();
0437   ppfjet_twr_elmttype_.clear();
0438   ppfjet_twr_hade_.clear();
0439   ppfjet_twr_frac_.clear();
0440   ppfjet_twr_dR_.clear();
0441   ppfjet_twr_clusterind_.clear();
0442   ppfjet_cluster_eta_.clear();
0443   ppfjet_cluster_phi_.clear();
0444   ppfjet_cluster_dR_.clear();
0445   ppfjet_candtrack_px_.clear();
0446   ppfjet_candtrack_py_.clear();
0447   ppfjet_candtrack_pz_.clear();
0448   ppfjet_candtrack_EcalE_.clear();
0449 
0450   std::map<int, std::pair<int, std::set<float>>> tpfjet_rechits;
0451   std::map<int, std::pair<int, std::set<float>>> ppfjet_rechits;
0452   std::map<float, int> tpfjet_clusters;
0453   std::map<float, int> ppfjet_clusters;
0454 
0455   // fill tag jet variables
0456   tpfjet_pt_ = pf_tag.jet()->pt();
0457   tpfjet_p_ = pf_tag.jet()->p();
0458   tpfjet_E_ = pf_tag.jet()->energy();
0459   tpfjet_eta_ = pf_tag.jet()->eta();
0460   tpfjet_phi_ = pf_tag.jet()->phi();
0461   tpfjet_scale_ = pf_tag.scale();
0462   tpfjet_area_ = pf_tag.jet()->jetArea();
0463   tpfjet_ntwrs_ = 0;
0464   tpfjet_ncandtracks_ = 0;
0465 
0466   tpfjet_jetID_ = 0;  // Not a loose, medium, or tight jet
0467   if (fabs(pf_tag.jet()->eta()) < 2.4) {
0468     if (pf_tag.jet()->chargedHadronEnergyFraction() > 0 && pf_tag.jet()->chargedMultiplicity() > 0 &&
0469         pf_tag.jet()->chargedEmEnergyFraction() < 0.99 &&
0470         (pf_tag.jet()->chargedMultiplicity() + pf_tag.jet()->neutralMultiplicity()) > 1) {
0471       if (pf_tag.jet()->neutralHadronEnergyFraction() < 0.9 && pf_tag.jet()->neutralEmEnergyFraction() < 0.9) {
0472         tpfjet_jetID_ = 3;  // Tight jet
0473       } else if (pf_tag.jet()->neutralHadronEnergyFraction() < 0.95 && pf_tag.jet()->neutralEmEnergyFraction() < 0.95) {
0474         tpfjet_jetID_ = 2;  // Medium jet
0475       } else if (pf_tag.jet()->neutralHadronEnergyFraction() < 0.99 && pf_tag.jet()->neutralEmEnergyFraction() < 0.99) {
0476         tpfjet_jetID_ = 1;  // Loose jet
0477       }
0478     }
0479   } else if ((pf_tag.jet()->chargedMultiplicity() + pf_tag.jet()->neutralMultiplicity()) > 1) {
0480     if (pf_tag.jet()->neutralHadronEnergyFraction() < 0.9 && pf_tag.jet()->neutralEmEnergyFraction() < 0.9) {
0481       tpfjet_jetID_ = 3;  // Tight jet
0482     } else if (pf_tag.jet()->neutralHadronEnergyFraction() < 0.95 && pf_tag.jet()->neutralEmEnergyFraction() < 0.95) {
0483       tpfjet_jetID_ = 2;  // Medium jet
0484     } else if (pf_tag.jet()->neutralHadronEnergyFraction() < 0.99 && pf_tag.jet()->neutralEmEnergyFraction() < 0.99) {
0485       tpfjet_jetID_ = 1;  // Loose jet
0486     }
0487   }
0488 
0489   /////////////////////////////////////////////
0490   // Get PF constituents and fill HCAL towers
0491   /////////////////////////////////////////////
0492 
0493   // Get tag PFCandidates
0494   std::vector<reco::PFCandidatePtr> tagconst = pf_tag.jet()->getPFConstituents();
0495   for (std::vector<reco::PFCandidatePtr>::const_iterator it = tagconst.begin(); it != tagconst.end(); ++it) {
0496     bool hasTrack = false;
0497     // Test PFCandidate type
0498     reco::PFCandidate::ParticleType candidateType = (*it)->particleId();
0499     switch (candidateType) {
0500       case reco::PFCandidate::X:
0501         tpfjet_unkown_E_ += (*it)->energy();
0502         tpfjet_unkown_px_ += (*it)->px();
0503         tpfjet_unkown_py_ += (*it)->py();
0504         tpfjet_unkown_pz_ += (*it)->pz();
0505         tpfjet_unkown_EcalE_ += (*it)->ecalEnergy();
0506         tpfjet_unkown_n_++;
0507         continue;
0508       case reco::PFCandidate::h: {
0509         tpfjet_had_E_.push_back((*it)->energy());
0510         tpfjet_had_px_.push_back((*it)->px());
0511         tpfjet_had_py_.push_back((*it)->py());
0512         tpfjet_had_pz_.push_back((*it)->pz());
0513         tpfjet_had_EcalE_.push_back((*it)->ecalEnergy());
0514         tpfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
0515         tpfjet_had_id_.push_back(0);
0516         tpfjet_had_ntwrs_.push_back(0);
0517         tpfjet_had_n_++;
0518 
0519         reco::TrackRef trackRef = (*it)->trackRef();
0520         if (trackRef.isNonnull()) {
0521           reco::Track track = *trackRef;
0522           tpfjet_candtrack_px_.push_back(track.px());
0523           tpfjet_candtrack_py_.push_back(track.py());
0524           tpfjet_candtrack_pz_.push_back(track.pz());
0525           tpfjet_candtrack_EcalE_.push_back((*it)->ecalEnergy());
0526           tpfjet_had_candtrackind_.push_back(tpfjet_ncandtracks_);
0527           hasTrack = true;
0528           tpfjet_ncandtracks_++;
0529         } else {
0530           tpfjet_had_candtrackind_.push_back(-2);
0531         }
0532       } break;
0533       case reco::PFCandidate::e:
0534         tpfjet_electron_E_ += (*it)->energy();
0535         tpfjet_electron_px_ += (*it)->px();
0536         tpfjet_electron_py_ += (*it)->py();
0537         tpfjet_electron_pz_ += (*it)->pz();
0538         tpfjet_electron_EcalE_ += (*it)->ecalEnergy();
0539         tpfjet_electron_n_++;
0540         continue;
0541       case reco::PFCandidate::mu:
0542         tpfjet_muon_E_ += (*it)->energy();
0543         tpfjet_muon_px_ += (*it)->px();
0544         tpfjet_muon_py_ += (*it)->py();
0545         tpfjet_muon_pz_ += (*it)->pz();
0546         tpfjet_muon_EcalE_ += (*it)->ecalEnergy();
0547         tpfjet_muon_n_++;
0548         continue;
0549       case reco::PFCandidate::gamma:
0550         tpfjet_photon_E_ += (*it)->energy();
0551         tpfjet_photon_px_ += (*it)->px();
0552         tpfjet_photon_py_ += (*it)->py();
0553         tpfjet_photon_pz_ += (*it)->pz();
0554         tpfjet_photon_EcalE_ += (*it)->ecalEnergy();
0555         tpfjet_photon_n_++;
0556         continue;
0557       case reco::PFCandidate::h0: {
0558         tpfjet_had_E_.push_back((*it)->energy());
0559         tpfjet_had_px_.push_back((*it)->px());
0560         tpfjet_had_py_.push_back((*it)->py());
0561         tpfjet_had_pz_.push_back((*it)->pz());
0562         tpfjet_had_EcalE_.push_back((*it)->ecalEnergy());
0563         tpfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
0564         tpfjet_had_id_.push_back(1);
0565         tpfjet_had_candtrackind_.push_back(-1);
0566         tpfjet_had_ntwrs_.push_back(0);
0567         tpfjet_had_n_++;
0568         break;
0569       }
0570       case reco::PFCandidate::h_HF: {
0571         tpfjet_had_E_.push_back((*it)->energy());
0572         tpfjet_had_px_.push_back((*it)->px());
0573         tpfjet_had_py_.push_back((*it)->py());
0574         tpfjet_had_pz_.push_back((*it)->pz());
0575         tpfjet_had_EcalE_.push_back((*it)->ecalEnergy());
0576         tpfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
0577         tpfjet_had_id_.push_back(2);
0578         tpfjet_had_candtrackind_.push_back(-1);
0579         tpfjet_had_ntwrs_.push_back(0);
0580         tpfjet_had_n_++;
0581         break;
0582       }
0583       case reco::PFCandidate::egamma_HF: {
0584         tpfjet_had_E_.push_back((*it)->energy());
0585         tpfjet_had_px_.push_back((*it)->px());
0586         tpfjet_had_py_.push_back((*it)->py());
0587         tpfjet_had_pz_.push_back((*it)->pz());
0588         tpfjet_had_EcalE_.push_back((*it)->ecalEnergy());
0589         tpfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
0590         tpfjet_had_id_.push_back(3);
0591         tpfjet_had_candtrackind_.push_back(-1);
0592         tpfjet_had_ntwrs_.push_back(0);
0593         tpfjet_had_n_++;
0594         break;
0595       }
0596     }
0597 
0598     std::map<int, int> twrietas;
0599     float HFHAD_E = 0;
0600     float HFEM_E = 0;
0601     int maxElement = (*it)->elementsInBlocks().size();
0602     for (int e = 0; e < maxElement; ++e) {
0603       // Get elements from block
0604       reco::PFBlockRef blockRef = (*it)->elementsInBlocks()[e].first;
0605       const edm::OwnVector<reco::PFBlockElement>& elements = blockRef->elements();
0606       for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
0607         if (elements[iEle].index() == (*it)->elementsInBlocks()[e].second) {
0608           if (elements[iEle].type() == reco::PFBlockElement::HCAL) {  // Element is HB or HE
0609             // Get cluster and hits
0610             reco::PFClusterRef clusterref = elements[iEle].clusterRef();
0611             reco::PFCluster cluster = *clusterref;
0612             double cluster_dR = deltaR(tpfjet_eta_, tpfjet_phi_, cluster.eta(), cluster.phi());
0613             if (tpfjet_clusters.count(cluster_dR) == 0) {
0614               tpfjet_clusters[cluster_dR] = tpfjet_cluster_n_;
0615               tpfjet_cluster_eta_.push_back(cluster.eta());
0616               tpfjet_cluster_phi_.push_back(cluster.phi());
0617               tpfjet_cluster_dR_.push_back(cluster_dR);
0618               tpfjet_cluster_n_++;
0619             }
0620             int cluster_ind = tpfjet_clusters[cluster_dR];
0621 
0622             std::vector<std::pair<DetId, float>> hitsAndFracs = cluster.hitsAndFractions();
0623 
0624             // Run over hits and match
0625             int nHits = hitsAndFracs.size();
0626             for (int iHit = 0; iHit < nHits; iHit++) {
0627               int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
0628 
0629               for (edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>::const_iterator ith =
0630                        hbhereco->begin();
0631                    ith != hbhereco->end();
0632                    ++ith) {
0633                 int etaPhiRecHit = getEtaPhi((*ith).id());
0634                 if (etaPhiPF == etaPhiRecHit) {
0635                   tpfjet_had_ntwrs_.at(tpfjet_had_n_ - 1)++;
0636                   if (tpfjet_rechits.count((*ith).id()) == 0) {
0637                     tpfjet_twr_ieta_.push_back((*ith).id().ieta());
0638                     tpfjet_twr_iphi_.push_back((*ith).id().iphi());
0639                     tpfjet_twr_depth_.push_back((*ith).id().depth());
0640                     tpfjet_twr_subdet_.push_back((*ith).id().subdet());
0641                     if (hitsAndFracs[iHit].second > 0.05 && (*ith).energy() > 0.0)
0642                       twrietas[(*ith).id().ieta()]++;
0643                     tpfjet_twr_hade_.push_back((*ith).energy());
0644                     tpfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
0645                     tpfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
0646                     tpfjet_twr_hadind_.push_back(tpfjet_had_n_ - 1);
0647                     tpfjet_twr_elmttype_.push_back(0);
0648                     tpfjet_twr_clusterind_.push_back(cluster_ind);
0649                     if (hasTrack) {
0650                       tpfjet_twr_candtrackind_.push_back(tpfjet_ncandtracks_ - 1);
0651                     } else {
0652                       tpfjet_twr_candtrackind_.push_back(-1);
0653                     }
0654                     switch ((*ith).id().subdet()) {
0655                       case HcalSubdetector::HcalBarrel: {
0656                         CaloCellGeometry::CornersVec cv = HBGeom->getCorners((*ith).id());
0657                         float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
0658                         float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
0659                         if (cv[0].phi() < cv[2].phi())
0660                           avgphi = (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) +
0661                                     static_cast<double>(cv[2].phi())) /
0662                                    2.0;
0663                         tpfjet_twr_dR_.push_back(deltaR(tpfjet_eta_, tpfjet_phi_, avgeta, avgphi));
0664                         break;
0665                       }
0666                       case HcalSubdetector::HcalEndcap: {
0667                         CaloCellGeometry::CornersVec cv = HEGeom->getCorners((*ith).id());
0668                         float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
0669                         float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
0670                         if (cv[0].phi() < cv[2].phi())
0671                           avgphi = (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) +
0672                                     static_cast<double>(cv[2].phi())) /
0673                                    2.0;
0674                         tpfjet_twr_dR_.push_back(deltaR(tpfjet_eta_, tpfjet_phi_, avgeta, avgphi));
0675                         break;
0676                       }
0677                       default:
0678                         tpfjet_twr_dR_.push_back(-1);
0679                         break;
0680                     }
0681                     tpfjet_rechits[(*ith).id()].first = tpfjet_ntwrs_;
0682                     ++tpfjet_ntwrs_;
0683                   } else if (tpfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0) {
0684                     tpfjet_twr_frac_.at(tpfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
0685                     if (cluster_dR <
0686                         tpfjet_cluster_dR_.at(tpfjet_twr_clusterind_.at(tpfjet_rechits[(*ith).id()].first))) {
0687                       tpfjet_twr_clusterind_.at(tpfjet_rechits[(*ith).id()].first) = cluster_ind;
0688                     }
0689                     tpfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
0690                   }
0691                 }                                                           // Test if ieta,iphi matches
0692               }                                                             // Loop over rechits
0693             }                                                               // Loop over hits
0694           }                                                                 // Test if element is from HCAL
0695           else if (elements[iEle].type() == reco::PFBlockElement::HFHAD) {  // Element is HF
0696             for (edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith =
0697                      hfreco->begin();
0698                  ith != hfreco->end();
0699                  ++ith) {
0700               if ((*ith).id().depth() == 1)
0701                 continue;  // Remove long fibers
0702               auto thisCell = HFGeom->getGeometry((*ith).id().rawId());
0703               const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
0704 
0705               bool passMatch = false;
0706               if ((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()) {
0707                 if ((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi())
0708                   passMatch = true;
0709                 else if (cv[0].phi() < cv[2].phi()) {
0710                   if ((*it)->phi() < cv[0].phi())
0711                     passMatch = true;
0712                   else if ((*it)->phi() > cv[2].phi())
0713                     passMatch = true;
0714                 }
0715               }
0716 
0717               if (passMatch) {
0718                 tpfjet_had_ntwrs_.at(tpfjet_had_n_ - 1)++;
0719                 tpfjet_twr_ieta_.push_back((*ith).id().ieta());
0720                 tpfjet_twr_iphi_.push_back((*ith).id().iphi());
0721                 tpfjet_twr_depth_.push_back((*ith).id().depth());
0722                 tpfjet_twr_subdet_.push_back((*ith).id().subdet());
0723                 tpfjet_twr_hade_.push_back((*ith).energy());
0724                 tpfjet_twr_frac_.push_back(1.0);
0725                 tpfjet_twr_hadind_.push_back(tpfjet_had_n_ - 1);
0726                 tpfjet_twr_elmttype_.push_back(1);
0727                 tpfjet_twr_clusterind_.push_back(-1);
0728                 tpfjet_twr_candtrackind_.push_back(-1);
0729                 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
0730                 float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
0731                 if (cv[0].phi() < cv[2].phi())
0732                   avgphi =
0733                       (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
0734                 tpfjet_twr_dR_.push_back(deltaR(tpfjet_eta_, tpfjet_phi_, avgeta, avgphi));
0735                 ++tpfjet_ntwrs_;
0736                 HFHAD_E += (*ith).energy();
0737               }
0738             }
0739           } else if (elements[iEle].type() == reco::PFBlockElement::HFEM) {  // Element is HF
0740             for (edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith =
0741                      hfreco->begin();
0742                  ith != hfreco->end();
0743                  ++ith) {
0744               if ((*ith).id().depth() == 2)
0745                 continue;  // Remove short fibers
0746               auto thisCell = HFGeom->getGeometry((*ith).id().rawId());
0747               const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
0748 
0749               bool passMatch = false;
0750               if ((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()) {
0751                 if ((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi())
0752                   passMatch = true;
0753                 else if (cv[0].phi() < cv[2].phi()) {
0754                   if ((*it)->phi() < cv[0].phi())
0755                     passMatch = true;
0756                   else if ((*it)->phi() > cv[2].phi())
0757                     passMatch = true;
0758                 }
0759               }
0760 
0761               if (passMatch) {
0762                 tpfjet_had_ntwrs_.at(tpfjet_had_n_ - 1)++;
0763                 tpfjet_twr_ieta_.push_back((*ith).id().ieta());
0764                 tpfjet_twr_iphi_.push_back((*ith).id().iphi());
0765                 tpfjet_twr_depth_.push_back((*ith).id().depth());
0766                 tpfjet_twr_subdet_.push_back((*ith).id().subdet());
0767                 tpfjet_twr_hade_.push_back((*ith).energy());
0768                 tpfjet_twr_frac_.push_back(1.0);
0769                 tpfjet_twr_hadind_.push_back(tpfjet_had_n_ - 1);
0770                 tpfjet_twr_elmttype_.push_back(2);
0771                 tpfjet_twr_clusterind_.push_back(-1);
0772                 tpfjet_twr_candtrackind_.push_back(-1);
0773                 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
0774                 float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
0775                 if (cv[0].phi() < cv[2].phi())
0776                   avgphi =
0777                       (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
0778                 tpfjet_twr_dR_.push_back(deltaR(tpfjet_eta_, tpfjet_phi_, avgeta, avgphi));
0779                 ++tpfjet_ntwrs_;
0780                 HFEM_E += (*ith).energy();
0781               }
0782             }
0783           } else if (elements[iEle].type() == reco::PFBlockElement::HO) {  // Element is HO
0784             reco::PFClusterRef clusterref = elements[iEle].clusterRef();
0785             reco::PFCluster cluster = *clusterref;
0786             double cluster_dR = deltaR(tpfjet_eta_, tpfjet_phi_, cluster.eta(), cluster.phi());
0787             if (tpfjet_clusters.count(cluster_dR) == 0) {
0788               tpfjet_clusters[cluster_dR] = tpfjet_cluster_n_;
0789               tpfjet_cluster_eta_.push_back(cluster.eta());
0790               tpfjet_cluster_phi_.push_back(cluster.phi());
0791               tpfjet_cluster_dR_.push_back(cluster_dR);
0792               tpfjet_cluster_n_++;
0793             }
0794             int cluster_ind = tpfjet_clusters[cluster_dR];
0795 
0796             std::vector<std::pair<DetId, float>> hitsAndFracs = cluster.hitsAndFractions();
0797             int nHits = hitsAndFracs.size();
0798             for (int iHit = 0; iHit < nHits; iHit++) {
0799               int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
0800 
0801               for (edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>::const_iterator ith =
0802                        horeco->begin();
0803                    ith != horeco->end();
0804                    ++ith) {
0805                 int etaPhiRecHit = getEtaPhi((*ith).id());
0806                 if (etaPhiPF == etaPhiRecHit) {
0807                   tpfjet_had_ntwrs_.at(tpfjet_had_n_ - 1)++;
0808                   if (tpfjet_rechits.count((*ith).id()) == 0) {
0809                     tpfjet_twr_ieta_.push_back((*ith).id().ieta());
0810                     tpfjet_twr_iphi_.push_back((*ith).id().iphi());
0811                     tpfjet_twr_depth_.push_back((*ith).id().depth());
0812                     tpfjet_twr_subdet_.push_back((*ith).id().subdet());
0813                     if (hitsAndFracs[iHit].second > 0.05 && (*ith).energy() > 0.0)
0814                       twrietas[(*ith).id().ieta()]++;
0815                     tpfjet_twr_hade_.push_back((*ith).energy());
0816                     tpfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
0817                     tpfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
0818                     tpfjet_twr_hadind_.push_back(tpfjet_had_n_ - 1);
0819                     tpfjet_twr_elmttype_.push_back(3);
0820                     tpfjet_twr_clusterind_.push_back(cluster_ind);
0821                     if (hasTrack) {
0822                       tpfjet_twr_candtrackind_.push_back(tpfjet_ncandtracks_ - 1);
0823                     } else {
0824                       tpfjet_twr_candtrackind_.push_back(-1);
0825                     }
0826                     auto thisCell = HOGeom->getGeometry((*ith).id().rawId());
0827                     const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
0828                     float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
0829                     float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
0830                     if (cv[0].phi() < cv[2].phi())
0831                       avgphi =
0832                           (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) /
0833                           2.0;
0834                     tpfjet_twr_dR_.push_back(deltaR(tpfjet_eta_, tpfjet_phi_, avgeta, avgphi));
0835                     tpfjet_rechits[(*ith).id()].first = tpfjet_ntwrs_;
0836                     ++tpfjet_ntwrs_;
0837                   } else if (tpfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0) {
0838                     tpfjet_twr_frac_.at(tpfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
0839                     if (cluster_dR <
0840                         tpfjet_cluster_dR_.at(tpfjet_twr_clusterind_.at(tpfjet_rechits[(*ith).id()].first))) {
0841                       tpfjet_twr_clusterind_.at(tpfjet_rechits[(*ith).id()].first) = cluster_ind;
0842                     }
0843                     tpfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
0844                   }
0845                 }  // Test if ieta,iphi match
0846               }    // Loop over rechits
0847             }      // Loop over hits
0848           }        // Test if element is HO
0849         }          // Test for right element index
0850       }            // Loop over elements
0851     }              // Loop over elements in blocks
0852 
0853     switch (candidateType) {
0854       case reco::PFCandidate::h_HF:
0855         tpfjet_had_emf_.push_back(HFEM_E / (HFEM_E + HFHAD_E));
0856         break;
0857       case reco::PFCandidate::egamma_HF:
0858         tpfjet_had_emf_.push_back(-1);
0859         break;
0860       default:
0861         tpfjet_had_emf_.push_back(-1);
0862         break;
0863     }
0864   }  // Loop over PF constitutents
0865 
0866   int tag_had_EcalE = 0;
0867   int tag_had_rawHcalE = 0;
0868   for (int i = 0; i < tpfjet_had_n_; i++) {
0869     tag_had_EcalE += tpfjet_had_EcalE_[i];
0870     tag_had_rawHcalE += tpfjet_had_rawHcalE_[i];
0871   }
0872   tpfjet_EMfrac_ = 1.0 - tag_had_rawHcalE / (tag_had_rawHcalE + tag_had_EcalE + tpfjet_unkown_E_ + tpfjet_electron_E_ +
0873                                              tpfjet_muon_E_ + tpfjet_photon_E_);
0874   tpfjet_hadEcalEfrac_ = tag_had_EcalE / (tag_had_rawHcalE + tag_had_EcalE + tpfjet_unkown_E_ + tpfjet_electron_E_ +
0875                                           tpfjet_muon_E_ + tpfjet_photon_E_);
0876 
0877   // fill probe jet variables
0878   ppfjet_pt_ = pf_probe.jet()->pt();
0879   ppfjet_p_ = pf_probe.jet()->p();
0880   ppfjet_E_ = pf_probe.jet()->energy();
0881   ppfjet_eta_ = pf_probe.jet()->eta();
0882   ppfjet_phi_ = pf_probe.jet()->phi();
0883   ppfjet_scale_ = pf_probe.scale();
0884   ppfjet_area_ = pf_probe.jet()->jetArea();
0885   ppfjet_ntwrs_ = 0;
0886   ppfjet_ncandtracks_ = 0;
0887 
0888   ppfjet_jetID_ = 0;  // Not a loose, medium, or tight jet
0889   if (fabs(pf_probe.jet()->eta()) < 2.4) {
0890     if (pf_probe.jet()->chargedHadronEnergyFraction() > 0 && pf_probe.jet()->chargedMultiplicity() > 0 &&
0891         pf_probe.jet()->chargedEmEnergyFraction() < 0.99 &&
0892         (pf_probe.jet()->chargedMultiplicity() + pf_probe.jet()->neutralMultiplicity()) > 1) {
0893       if (pf_probe.jet()->neutralHadronEnergyFraction() < 0.9 && pf_probe.jet()->neutralEmEnergyFraction() < 0.9) {
0894         ppfjet_jetID_ = 3;  // Tight jet
0895       } else if (pf_probe.jet()->neutralHadronEnergyFraction() < 0.95 &&
0896                  pf_probe.jet()->neutralEmEnergyFraction() < 0.95) {
0897         ppfjet_jetID_ = 2;  // Medium jet
0898       } else if (pf_probe.jet()->neutralHadronEnergyFraction() < 0.99 &&
0899                  pf_probe.jet()->neutralEmEnergyFraction() < 0.99) {
0900         ppfjet_jetID_ = 1;  // Loose jet
0901       }
0902     }
0903   } else if ((pf_probe.jet()->chargedMultiplicity() + pf_probe.jet()->neutralMultiplicity()) > 1) {
0904     if (pf_probe.jet()->neutralHadronEnergyFraction() < 0.9 && pf_probe.jet()->neutralEmEnergyFraction() < 0.9) {
0905       ppfjet_jetID_ = 3;  // Tight jet
0906     } else if (pf_probe.jet()->neutralHadronEnergyFraction() < 0.95 &&
0907                pf_probe.jet()->neutralEmEnergyFraction() < 0.95) {
0908       ppfjet_jetID_ = 2;  // Medium jet
0909     } else if (pf_probe.jet()->neutralHadronEnergyFraction() < 0.99 &&
0910                pf_probe.jet()->neutralEmEnergyFraction() < 0.99) {
0911       ppfjet_jetID_ = 1;  // Loose jet
0912     }
0913   }
0914 
0915   // Get PF constituents and fill HCAL towers
0916   std::vector<reco::PFCandidatePtr> probeconst = pf_probe.jet()->getPFConstituents();
0917   for (std::vector<reco::PFCandidatePtr>::const_iterator it = probeconst.begin(); it != probeconst.end(); ++it) {
0918     bool hasTrack = false;
0919     reco::PFCandidate::ParticleType candidateType = (*it)->particleId();
0920     switch (candidateType) {
0921       case reco::PFCandidate::X:
0922         ppfjet_unkown_E_ += (*it)->energy();
0923         ppfjet_unkown_px_ += (*it)->px();
0924         ppfjet_unkown_py_ += (*it)->py();
0925         ppfjet_unkown_pz_ += (*it)->pz();
0926         ppfjet_unkown_EcalE_ += (*it)->ecalEnergy();
0927         ppfjet_unkown_n_++;
0928         continue;
0929       case reco::PFCandidate::h: {
0930         ppfjet_had_E_.push_back((*it)->energy());
0931         ppfjet_had_px_.push_back((*it)->px());
0932         ppfjet_had_py_.push_back((*it)->py());
0933         ppfjet_had_pz_.push_back((*it)->pz());
0934         ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
0935         ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
0936         ppfjet_had_id_.push_back(0);
0937         ppfjet_had_ntwrs_.push_back(0);
0938         ppfjet_had_n_++;
0939 
0940         reco::TrackRef trackRef = (*it)->trackRef();
0941         if (trackRef.isNonnull()) {
0942           reco::Track track = *trackRef;
0943           ppfjet_candtrack_px_.push_back(track.px());
0944           ppfjet_candtrack_py_.push_back(track.py());
0945           ppfjet_candtrack_pz_.push_back(track.pz());
0946           ppfjet_candtrack_EcalE_.push_back((*it)->ecalEnergy());
0947           ppfjet_had_candtrackind_.push_back(ppfjet_ncandtracks_);
0948           hasTrack = true;
0949           ppfjet_ncandtracks_++;
0950         } else {
0951           ppfjet_had_candtrackind_.push_back(-2);
0952         }
0953       } break;
0954       case reco::PFCandidate::e:
0955         ppfjet_electron_E_ += (*it)->energy();
0956         ppfjet_electron_px_ += (*it)->px();
0957         ppfjet_electron_py_ += (*it)->py();
0958         ppfjet_electron_pz_ += (*it)->pz();
0959         ppfjet_electron_EcalE_ += (*it)->ecalEnergy();
0960         ppfjet_electron_n_++;
0961         continue;
0962       case reco::PFCandidate::mu:
0963         ppfjet_muon_E_ += (*it)->energy();
0964         ppfjet_muon_px_ += (*it)->px();
0965         ppfjet_muon_py_ += (*it)->py();
0966         ppfjet_muon_pz_ += (*it)->pz();
0967         ppfjet_muon_EcalE_ += (*it)->ecalEnergy();
0968         ppfjet_muon_n_++;
0969         continue;
0970       case reco::PFCandidate::gamma:
0971         ppfjet_photon_E_ += (*it)->energy();
0972         ppfjet_photon_px_ += (*it)->px();
0973         ppfjet_photon_py_ += (*it)->py();
0974         ppfjet_photon_pz_ += (*it)->pz();
0975         ppfjet_photon_EcalE_ += (*it)->ecalEnergy();
0976         ppfjet_photon_n_++;
0977         continue;
0978       case reco::PFCandidate::h0: {
0979         ppfjet_had_E_.push_back((*it)->energy());
0980         ppfjet_had_px_.push_back((*it)->px());
0981         ppfjet_had_py_.push_back((*it)->py());
0982         ppfjet_had_pz_.push_back((*it)->pz());
0983         ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
0984         ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
0985         ppfjet_had_id_.push_back(1);
0986         ppfjet_had_candtrackind_.push_back(-1);
0987         ppfjet_had_ntwrs_.push_back(0);
0988         ppfjet_had_n_++;
0989         break;
0990       }
0991       case reco::PFCandidate::h_HF: {
0992         ppfjet_had_E_.push_back((*it)->energy());
0993         ppfjet_had_px_.push_back((*it)->px());
0994         ppfjet_had_py_.push_back((*it)->py());
0995         ppfjet_had_pz_.push_back((*it)->pz());
0996         ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
0997         ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
0998         ppfjet_had_id_.push_back(2);
0999         ppfjet_had_candtrackind_.push_back(-1);
1000         ppfjet_had_ntwrs_.push_back(0);
1001         ppfjet_had_n_++;
1002         break;
1003       }
1004       case reco::PFCandidate::egamma_HF: {
1005         ppfjet_had_E_.push_back((*it)->energy());
1006         ppfjet_had_px_.push_back((*it)->px());
1007         ppfjet_had_py_.push_back((*it)->py());
1008         ppfjet_had_pz_.push_back((*it)->pz());
1009         ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
1010         ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
1011         ppfjet_had_id_.push_back(3);
1012         ppfjet_had_candtrackind_.push_back(-1);
1013         ppfjet_had_ntwrs_.push_back(0);
1014         ppfjet_had_n_++;
1015         break;
1016       }
1017     }
1018 
1019     float HFHAD_E = 0;
1020     float HFEM_E = 0;
1021     int maxElement = (*it)->elementsInBlocks().size();
1022     for (int e = 0; e < maxElement; ++e) {
1023       // Get elements from block
1024       reco::PFBlockRef blockRef = (*it)->elementsInBlocks()[e].first;
1025       const edm::OwnVector<reco::PFBlockElement>& elements = blockRef->elements();
1026       for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
1027         if (elements[iEle].index() == (*it)->elementsInBlocks()[e].second) {
1028           if (elements[iEle].type() == reco::PFBlockElement::HCAL) {  // Element is HB or HE
1029             // Get cluster and hits
1030             reco::PFClusterRef clusterref = elements[iEle].clusterRef();
1031             reco::PFCluster cluster = *clusterref;
1032             double cluster_dR = deltaR(ppfjet_eta_, ppfjet_phi_, cluster.eta(), cluster.phi());
1033             if (ppfjet_clusters.count(cluster_dR) == 0) {
1034               ppfjet_clusters[cluster_dR] = ppfjet_cluster_n_;
1035               ppfjet_cluster_eta_.push_back(cluster.eta());
1036               ppfjet_cluster_phi_.push_back(cluster.phi());
1037               ppfjet_cluster_dR_.push_back(cluster_dR);
1038               ppfjet_cluster_n_++;
1039             }
1040             int cluster_ind = ppfjet_clusters[cluster_dR];
1041             std::vector<std::pair<DetId, float>> hitsAndFracs = cluster.hitsAndFractions();
1042 
1043             // Run over hits and match
1044             int nHits = hitsAndFracs.size();
1045             for (int iHit = 0; iHit < nHits; iHit++) {
1046               int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
1047 
1048               for (edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>::const_iterator ith =
1049                        hbhereco->begin();
1050                    ith != hbhereco->end();
1051                    ++ith) {
1052                 int etaPhiRecHit = getEtaPhi((*ith).id());
1053                 if (etaPhiPF == etaPhiRecHit) {
1054                   ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1055                   if (ppfjet_rechits.count((*ith).id()) == 0) {
1056                     ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1057                     ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1058                     ppfjet_twr_depth_.push_back((*ith).id().depth());
1059                     ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1060                     ppfjet_twr_hade_.push_back((*ith).energy());
1061                     ppfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
1062                     ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1063                     ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1064                     ppfjet_twr_elmttype_.push_back(0);
1065                     ppfjet_twr_clusterind_.push_back(cluster_ind);
1066                     if (hasTrack) {
1067                       ppfjet_twr_candtrackind_.push_back(ppfjet_ncandtracks_ - 1);
1068                     } else {
1069                       ppfjet_twr_candtrackind_.push_back(-1);
1070                     }
1071                     switch ((*ith).id().subdet()) {
1072                       case HcalSubdetector::HcalBarrel: {
1073                         CaloCellGeometry::CornersVec cv = HBGeom->getCorners((*ith).id());
1074                         float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1075                         float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1076                         if (cv[0].phi() < cv[2].phi())
1077                           avgphi = (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) +
1078                                     static_cast<double>(cv[2].phi())) /
1079                                    2.0;
1080                         ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1081                         break;
1082                       }
1083                       case HcalSubdetector::HcalEndcap: {
1084                         CaloCellGeometry::CornersVec cv = HEGeom->getCorners((*ith).id());
1085                         float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1086                         float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1087                         if (cv[0].phi() < cv[2].phi())
1088                           avgphi = (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) +
1089                                     static_cast<double>(cv[2].phi())) /
1090                                    2.0;
1091                         ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1092                         break;
1093                       }
1094                       default:
1095                         ppfjet_twr_dR_.push_back(-1);
1096                         break;
1097                     }
1098                     ppfjet_rechits[(*ith).id()].first = ppfjet_ntwrs_;
1099                     ++ppfjet_ntwrs_;
1100                   } else if (ppfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0) {
1101                     ppfjet_twr_frac_.at(ppfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
1102                     if (cluster_dR <
1103                         ppfjet_cluster_dR_.at(ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first))) {
1104                       ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first) = cluster_ind;
1105                     }
1106                     ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1107                   }
1108                 }                                                           // Test if ieta,iphi matches
1109               }                                                             // Loop over rechits
1110             }                                                               // Loop over hits
1111           }                                                                 // Test if element is from HCAL
1112           else if (elements[iEle].type() == reco::PFBlockElement::HFHAD) {  // Element is HF
1113             for (edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith =
1114                      hfreco->begin();
1115                  ith != hfreco->end();
1116                  ++ith) {
1117               if ((*ith).id().depth() == 1)
1118                 continue;  // Remove long fibers
1119               auto thisCell = HFGeom->getGeometry((*ith).id().rawId());
1120               const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
1121 
1122               bool passMatch = false;
1123               if ((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()) {
1124                 if ((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi())
1125                   passMatch = true;
1126                 else if (cv[0].phi() < cv[2].phi()) {
1127                   if ((*it)->phi() < cv[0].phi())
1128                     passMatch = true;
1129                   else if ((*it)->phi() > cv[2].phi())
1130                     passMatch = true;
1131                 }
1132               }
1133 
1134               if (passMatch) {
1135                 ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1136                 ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1137                 ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1138                 ppfjet_twr_depth_.push_back((*ith).id().depth());
1139                 ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1140                 ppfjet_twr_hade_.push_back((*ith).energy());
1141                 ppfjet_twr_frac_.push_back(1.0);
1142                 ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1143                 ppfjet_twr_elmttype_.push_back(1);
1144                 ppfjet_twr_clusterind_.push_back(-1);
1145                 ppfjet_twr_candtrackind_.push_back(-1);
1146                 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1147                 float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1148                 if (cv[0].phi() < cv[2].phi())
1149                   avgphi =
1150                       (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1151                 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1152                 ++ppfjet_ntwrs_;
1153                 HFHAD_E += (*ith).energy();
1154               }
1155             }
1156           } else if (elements[iEle].type() == reco::PFBlockElement::HFEM) {  // Element is HF
1157             for (edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith =
1158                      hfreco->begin();
1159                  ith != hfreco->end();
1160                  ++ith) {
1161               if ((*ith).id().depth() == 2)
1162                 continue;  // Remove short fibers
1163               auto thisCell = HFGeom->getGeometry((*ith).id().rawId());
1164               const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
1165 
1166               bool passMatch = false;
1167               if ((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()) {
1168                 if ((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi())
1169                   passMatch = true;
1170                 else if (cv[0].phi() < cv[2].phi()) {
1171                   if ((*it)->phi() < cv[0].phi())
1172                     passMatch = true;
1173                   else if ((*it)->phi() > cv[2].phi())
1174                     passMatch = true;
1175                 }
1176               }
1177 
1178               if (passMatch) {
1179                 ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1180                 ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1181                 ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1182                 ppfjet_twr_depth_.push_back((*ith).id().depth());
1183                 ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1184                 ppfjet_twr_hade_.push_back((*ith).energy());
1185                 ppfjet_twr_frac_.push_back(1.0);
1186                 ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1187                 ppfjet_twr_elmttype_.push_back(2);
1188                 ppfjet_twr_clusterind_.push_back(-1);
1189                 ppfjet_twr_candtrackind_.push_back(-1);
1190                 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1191                 float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1192                 if (cv[0].phi() < cv[2].phi())
1193                   avgphi =
1194                       (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1195                 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1196                 ++ppfjet_ntwrs_;
1197                 HFEM_E += (*ith).energy();
1198               }
1199             }
1200           } else if (elements[iEle].type() == reco::PFBlockElement::HO) {  // Element is HO
1201             reco::PFClusterRef clusterref = elements[iEle].clusterRef();
1202             reco::PFCluster cluster = *clusterref;
1203             double cluster_dR = deltaR(ppfjet_eta_, ppfjet_phi_, cluster.eta(), cluster.phi());
1204             if (ppfjet_clusters.count(cluster_dR) == 0) {
1205               ppfjet_clusters[cluster_dR] = ppfjet_cluster_n_;
1206               ppfjet_cluster_eta_.push_back(cluster.eta());
1207               ppfjet_cluster_phi_.push_back(cluster.phi());
1208               ppfjet_cluster_dR_.push_back(cluster_dR);
1209               ppfjet_cluster_n_++;
1210             }
1211             int cluster_ind = ppfjet_clusters[cluster_dR];
1212 
1213             std::vector<std::pair<DetId, float>> hitsAndFracs = cluster.hitsAndFractions();
1214             int nHits = hitsAndFracs.size();
1215             for (int iHit = 0; iHit < nHits; iHit++) {
1216               int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
1217 
1218               for (edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>::const_iterator ith =
1219                        horeco->begin();
1220                    ith != horeco->end();
1221                    ++ith) {
1222                 int etaPhiRecHit = getEtaPhi((*ith).id());
1223                 if (etaPhiPF == etaPhiRecHit) {
1224                   ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1225                   if (ppfjet_rechits.count((*ith).id()) == 0) {
1226                     ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1227                     ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1228                     ppfjet_twr_depth_.push_back((*ith).id().depth());
1229                     ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1230                     ppfjet_twr_hade_.push_back((*ith).energy());
1231                     ppfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
1232                     ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1233                     ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1234                     ppfjet_twr_elmttype_.push_back(3);
1235                     ppfjet_twr_clusterind_.push_back(cluster_ind);
1236                     if (hasTrack) {
1237                       ppfjet_twr_candtrackind_.push_back(ppfjet_ncandtracks_ - 1);
1238                     } else {
1239                       ppfjet_twr_candtrackind_.push_back(-1);
1240                     }
1241                     auto thisCell = HOGeom->getGeometry((*ith).id().rawId());
1242                     const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
1243                     float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1244                     float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1245                     if (cv[0].phi() < cv[2].phi())
1246                       avgphi =
1247                           (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) /
1248                           2.0;
1249                     ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1250                     ppfjet_rechits[(*ith).id()].first = ppfjet_ntwrs_;
1251                     ++ppfjet_ntwrs_;
1252                   } else if (ppfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0) {
1253                     ppfjet_twr_frac_.at(ppfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
1254                     if (cluster_dR <
1255                         ppfjet_cluster_dR_.at(ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first))) {
1256                       ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first) = cluster_ind;
1257                     }
1258                     ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1259                   }
1260                 }  // Test if ieta,iphi match
1261               }    // Loop over rechits
1262             }      // Loop over hits
1263           }        // Test if element is from HO
1264         }          // Test for right element index
1265       }            // Loop over elements
1266     }              // Loop over elements in blocks
1267     switch (candidateType) {
1268       case reco::PFCandidate::h_HF:
1269         ppfjet_had_emf_.push_back(HFEM_E / (HFEM_E + HFHAD_E));
1270         break;
1271       case reco::PFCandidate::egamma_HF:
1272         ppfjet_had_emf_.push_back(-1);
1273         break;
1274       default:
1275         ppfjet_had_emf_.push_back(-1);
1276         break;
1277     }
1278   }  // Loop over PF constitutents
1279 
1280   int probe_had_EcalE = 0;
1281   int probe_had_rawHcalE = 0;
1282   for (int i = 0; i < ppfjet_had_n_; i++) {
1283     probe_had_EcalE += ppfjet_had_EcalE_[i];
1284     probe_had_rawHcalE += ppfjet_had_rawHcalE_[i];
1285   }
1286   ppfjet_EMfrac_ = 1.0 - probe_had_rawHcalE / (probe_had_rawHcalE + probe_had_EcalE + ppfjet_unkown_E_ +
1287                                                ppfjet_electron_E_ + ppfjet_muon_E_ + ppfjet_photon_E_);
1288   ppfjet_hadEcalEfrac_ = probe_had_EcalE / (probe_had_rawHcalE + probe_had_EcalE + ppfjet_unkown_E_ +
1289                                             ppfjet_electron_E_ + ppfjet_muon_E_ + ppfjet_photon_E_);
1290 
1291   // fill dijet variables
1292   pf_dijet_deta_ = std::fabs(std::fabs(pf_tag.jet()->eta()) - std::fabs(pf_probe.jet()->eta()));
1293   pf_dijet_dphi_ = pf_tag.jet()->phi() - pf_probe.jet()->phi();
1294   if (pf_dijet_dphi_ > 3.1415)
1295     pf_dijet_dphi_ = 6.2832 - pf_dijet_dphi_;
1296   pf_dijet_balance_ = (tpfjet_pt_ - ppfjet_pt_) / (tpfjet_pt_ + ppfjet_pt_);
1297 
1298   tree_->Fill();
1299 
1300   return;
1301 }
1302 
1303 // ------------ method called once each job just before starting event loop  ------------
1304 void DiJetAnalyzer::beginJob() {
1305   // book histograms
1306 
1307   edm::Service<TFileService> fs;
1308   h_PassSelPF_ = fs->make<TH1D>("h_PassSelectionPF", "Selection Pass Failures PFJets", 200, -0.5, 199.5);
1309 
1310   tree_ = fs->make<TTree>("dijettree", "tree for dijet balancing");
1311 
1312   tree_->Branch("tpfjet_pt", &tpfjet_pt_, "tpfjet_pt/F");
1313   tree_->Branch("tpfjet_p", &tpfjet_p_, "tpfjet_p/F");
1314   tree_->Branch("tpfjet_E", &tpfjet_E_, "tpfjet_E/F");
1315   tree_->Branch("tpfjet_eta", &tpfjet_eta_, "tpfjet_eta/F");
1316   tree_->Branch("tpfjet_phi", &tpfjet_phi_, "tpfjet_phi/F");
1317   tree_->Branch("tpfjet_EMfrac", &tpfjet_EMfrac_, "tpfjet_EMfrac/F");
1318   tree_->Branch("tpfjet_hadEcalEfrac", &tpfjet_hadEcalEfrac_, "tpfjet_hadEcalEfrac/F");
1319   tree_->Branch("tpfjet_scale", &tpfjet_scale_, "tpfjet_scale/F");
1320   tree_->Branch("tpfjet_area", &tpfjet_area_, "tpfjet_area/F");
1321   tree_->Branch("tpfjet_jetID", &tpfjet_jetID_, "tpfjet_jetID/I");
1322   tree_->Branch("tpfjet_unkown_E", &tpfjet_unkown_E_, "tpfjet_unkown_E/F");
1323   tree_->Branch("tpfjet_electron_E", &tpfjet_electron_E_, "tpfjet_electron_E/F");
1324   tree_->Branch("tpfjet_muon_E", &tpfjet_muon_E_, "tpfjet_muon_E/F");
1325   tree_->Branch("tpfjet_photon_E", &tpfjet_photon_E_, "tpfjet_photon_E/F");
1326   tree_->Branch("tpfjet_unkown_px", &tpfjet_unkown_px_, "tpfjet_unkown_px/F");
1327   tree_->Branch("tpfjet_electron_px", &tpfjet_electron_px_, "tpfjet_electron_px/F");
1328   tree_->Branch("tpfjet_muon_px", &tpfjet_muon_px_, "tpfjet_muon_px/F");
1329   tree_->Branch("tpfjet_photon_px", &tpfjet_photon_px_, "tpfjet_photon_px/F");
1330   tree_->Branch("tpfjet_unkown_py", &tpfjet_unkown_py_, "tpfjet_unkown_py/F");
1331   tree_->Branch("tpfjet_electron_py", &tpfjet_electron_py_, "tpfjet_electron_py/F");
1332   tree_->Branch("tpfjet_muon_py", &tpfjet_muon_py_, "tpfjet_muon_py/F");
1333   tree_->Branch("tpfjet_photon_py", &tpfjet_photon_py_, "tpfjet_photon_py/F");
1334   tree_->Branch("tpfjet_unkown_pz", &tpfjet_unkown_pz_, "tpfjet_unkown_pz/F");
1335   tree_->Branch("tpfjet_electron_pz", &tpfjet_electron_pz_, "tpfjet_electron_pz/F");
1336   tree_->Branch("tpfjet_muon_pz", &tpfjet_muon_pz_, "tpfjet_muon_pz/F");
1337   tree_->Branch("tpfjet_photon_pz", &tpfjet_photon_pz_, "tpfjet_photon_pz/F");
1338   tree_->Branch("tpfjet_unkown_EcalE", &tpfjet_unkown_EcalE_, "tpfjet_unkown_EcalE/F");
1339   tree_->Branch("tpfjet_electron_EcalE", &tpfjet_electron_EcalE_, "tpfjet_electron_EcalE/F");
1340   tree_->Branch("tpfjet_muon_EcalE", &tpfjet_muon_EcalE_, "tpfjet_muon_EcalE/F");
1341   tree_->Branch("tpfjet_photon_EcalE", &tpfjet_photon_EcalE_, "tpfjet_photon_EcalE/F");
1342   tree_->Branch("tpfjet_unkown_n", &tpfjet_unkown_n_, "tpfjet_unkown_n/I");
1343   tree_->Branch("tpfjet_electron_n", &tpfjet_electron_n_, "tpfjet_electron_n/I");
1344   tree_->Branch("tpfjet_muon_n", &tpfjet_muon_n_, "tpfjet_muon_n/I");
1345   tree_->Branch("tpfjet_photon_n", &tpfjet_photon_n_, "tpfjet_photon_n/I");
1346   tree_->Branch("tpfjet_had_n", &tpfjet_had_n_, "tpfjet_had_n/I");
1347   tree_->Branch("tpfjet_had_E", &tpfjet_had_E_);
1348   tree_->Branch("tpfjet_had_px", &tpfjet_had_px_);
1349   tree_->Branch("tpfjet_had_py", &tpfjet_had_py_);
1350   tree_->Branch("tpfjet_had_pz", &tpfjet_had_pz_);
1351   tree_->Branch("tpfjet_had_EcalE", &tpfjet_had_EcalE_);
1352   tree_->Branch("tpfjet_had_rawHcalE", &tpfjet_had_rawHcalE_);
1353   tree_->Branch("tpfjet_had_emf", &tpfjet_had_emf_);
1354   tree_->Branch("tpfjet_had_id", &tpfjet_had_id_);
1355   tree_->Branch("tpfjet_had_candtrackind", &tpfjet_had_candtrackind_);
1356   tree_->Branch("tpfjet_had_ntwrs", &tpfjet_had_ntwrs_);
1357   tree_->Branch("tpfjet_ntwrs", &tpfjet_ntwrs_, "tpfjet_ntwrs/I");
1358   tree_->Branch("tpfjet_twr_ieta", &tpfjet_twr_ieta_);
1359   tree_->Branch("tpfjet_twr_iphi", &tpfjet_twr_iphi_);
1360   tree_->Branch("tpfjet_twr_depth", &tpfjet_twr_depth_);
1361   tree_->Branch("tpfjet_twr_subdet", &tpfjet_twr_subdet_);
1362   tree_->Branch("tpfjet_twr_hade", &tpfjet_twr_hade_);
1363   tree_->Branch("tpfjet_twr_frac", &tpfjet_twr_frac_);
1364   tree_->Branch("tpfjet_twr_candtrackind", &tpfjet_twr_candtrackind_);
1365   tree_->Branch("tpfjet_twr_hadind", &tpfjet_twr_hadind_);
1366   tree_->Branch("tpfjet_twr_elmttype", &tpfjet_twr_elmttype_);
1367   tree_->Branch("tpfjet_twr_dR", &tpfjet_twr_dR_);
1368   tree_->Branch("tpfjet_twr_clusterind", &tpfjet_twr_clusterind_);
1369   tree_->Branch("tpfjet_cluster_n", &tpfjet_cluster_n_, "tpfjet_cluster_n/I");
1370   tree_->Branch("tpfjet_cluster_eta", &tpfjet_cluster_eta_);
1371   tree_->Branch("tpfjet_cluster_phi", &tpfjet_cluster_phi_);
1372   tree_->Branch("tpfjet_cluster_dR", &tpfjet_cluster_dR_);
1373   tree_->Branch("tpfjet_ncandtracks", &tpfjet_ncandtracks_, "tpfjet_ncandtracks/I");
1374   tree_->Branch("tpfjet_candtrack_px", &tpfjet_candtrack_px_);
1375   tree_->Branch("tpfjet_candtrack_py", &tpfjet_candtrack_py_);
1376   tree_->Branch("tpfjet_candtrack_pz", &tpfjet_candtrack_pz_);
1377   tree_->Branch("tpfjet_candtrack_EcalE", &tpfjet_candtrack_EcalE_);
1378   tree_->Branch("ppfjet_pt", &ppfjet_pt_, "ppfjet_pt/F");
1379   tree_->Branch("ppfjet_p", &ppfjet_p_, "ppfjet_p/F");
1380   tree_->Branch("ppfjet_E", &ppfjet_E_, "ppfjet_E/F");
1381   tree_->Branch("ppfjet_eta", &ppfjet_eta_, "ppfjet_eta/F");
1382   tree_->Branch("ppfjet_phi", &ppfjet_phi_, "ppfjet_phi/F");
1383   tree_->Branch("ppfjet_EMfrac", &ppfjet_EMfrac_, "ppfjet_EMfrac/F");
1384   tree_->Branch("ppfjet_hadEcalEfrac", &ppfjet_hadEcalEfrac_, "ppfjet_hadEcalEfrac/F");
1385   tree_->Branch("ppfjet_scale", &ppfjet_scale_, "ppfjet_scale/F");
1386   tree_->Branch("ppfjet_area", &ppfjet_area_, "ppfjet_area/F");
1387   tree_->Branch("ppfjet_jetID", &ppfjet_jetID_, "ppfjet_jetID/I");
1388   tree_->Branch("ppfjet_unkown_E", &ppfjet_unkown_E_, "ppfjet_unkown_E/F");
1389   tree_->Branch("ppfjet_electron_E", &ppfjet_electron_E_, "ppfjet_electron_E/F");
1390   tree_->Branch("ppfjet_muon_E", &ppfjet_muon_E_, "ppfjet_muon_E/F");
1391   tree_->Branch("ppfjet_photon_E", &ppfjet_photon_E_, "ppfjet_photon_E/F");
1392   tree_->Branch("ppfjet_unkown_px", &ppfjet_unkown_px_, "ppfjet_unkown_px/F");
1393   tree_->Branch("ppfjet_electron_px", &ppfjet_electron_px_, "ppfjet_electron_px/F");
1394   tree_->Branch("ppfjet_muon_px", &ppfjet_muon_px_, "ppfjet_muon_px/F");
1395   tree_->Branch("ppfjet_photon_px", &ppfjet_photon_px_, "ppfjet_photon_px/F");
1396   tree_->Branch("ppfjet_unkown_py", &ppfjet_unkown_py_, "ppfjet_unkown_py/F");
1397   tree_->Branch("ppfjet_electron_py", &ppfjet_electron_py_, "ppfjet_electron_py/F");
1398   tree_->Branch("ppfjet_muon_py", &ppfjet_muon_py_, "ppfjet_muon_py/F");
1399   tree_->Branch("ppfjet_photon_py", &ppfjet_photon_py_, "ppfjet_photon_py/F");
1400   tree_->Branch("ppfjet_unkown_pz", &ppfjet_unkown_pz_, "ppfjet_unkown_pz/F");
1401   tree_->Branch("ppfjet_electron_pz", &ppfjet_electron_pz_, "ppfjet_electron_pz/F");
1402   tree_->Branch("ppfjet_muon_pz", &ppfjet_muon_pz_, "ppfjet_muon_pz/F");
1403   tree_->Branch("ppfjet_photon_pz", &ppfjet_photon_pz_, "ppfjet_photon_pz/F");
1404   tree_->Branch("ppfjet_unkown_EcalE", &ppfjet_unkown_EcalE_, "ppfjet_unkown_EcalE/F");
1405   tree_->Branch("ppfjet_electron_EcalE", &ppfjet_electron_EcalE_, "ppfjet_electron_EcalE/F");
1406   tree_->Branch("ppfjet_muon_EcalE", &ppfjet_muon_EcalE_, "ppfjet_muon_EcalE/F");
1407   tree_->Branch("ppfjet_photon_EcalE", &ppfjet_photon_EcalE_, "ppfjet_photon_EcalE/F");
1408   tree_->Branch("ppfjet_unkown_n", &ppfjet_unkown_n_, "ppfjet_unkown_n/I");
1409   tree_->Branch("ppfjet_electron_n", &ppfjet_electron_n_, "ppfjet_electron_n/I");
1410   tree_->Branch("ppfjet_muon_n", &ppfjet_muon_n_, "ppfjet_muon_n/I");
1411   tree_->Branch("ppfjet_photon_n", &ppfjet_photon_n_, "ppfjet_photon_n/I");
1412   tree_->Branch("ppfjet_had_n", &ppfjet_had_n_, "ppfjet_had_n/I");
1413   tree_->Branch("ppfjet_had_E", &ppfjet_had_E_);
1414   tree_->Branch("ppfjet_had_px", &ppfjet_had_px_);
1415   tree_->Branch("ppfjet_had_py", &ppfjet_had_py_);
1416   tree_->Branch("ppfjet_had_pz", &ppfjet_had_pz_);
1417   tree_->Branch("ppfjet_had_EcalE", &ppfjet_had_EcalE_);
1418   tree_->Branch("ppfjet_had_rawHcalE", &ppfjet_had_rawHcalE_);
1419   tree_->Branch("ppfjet_had_emf", &ppfjet_had_emf_);
1420   tree_->Branch("ppfjet_had_id", &ppfjet_had_id_);
1421   tree_->Branch("ppfjet_had_candtrackind", &ppfjet_had_candtrackind_);
1422   tree_->Branch("ppfjet_had_ntwrs", &ppfjet_had_ntwrs_);
1423   tree_->Branch("ppfjet_ntwrs", &ppfjet_ntwrs_, "ppfjet_ntwrs/I");
1424   tree_->Branch("ppfjet_twr_ieta", &ppfjet_twr_ieta_);
1425   tree_->Branch("ppfjet_twr_iphi", &ppfjet_twr_iphi_);
1426   tree_->Branch("ppfjet_twr_depth", &ppfjet_twr_depth_);
1427   tree_->Branch("ppfjet_twr_subdet", &ppfjet_twr_subdet_);
1428   tree_->Branch("ppfjet_twr_hade", &ppfjet_twr_hade_);
1429   tree_->Branch("ppfjet_twr_frac", &ppfjet_twr_frac_);
1430   tree_->Branch("ppfjet_twr_candtrackind", &ppfjet_twr_candtrackind_);
1431   tree_->Branch("ppfjet_twr_hadind", &ppfjet_twr_hadind_);
1432   tree_->Branch("ppfjet_twr_elmttype", &ppfjet_twr_elmttype_);
1433   tree_->Branch("ppfjet_twr_dR", &ppfjet_twr_dR_);
1434   tree_->Branch("ppfjet_twr_clusterind", &ppfjet_twr_clusterind_);
1435   tree_->Branch("ppfjet_cluster_n", &ppfjet_cluster_n_, "ppfjet_cluster_n/I");
1436   tree_->Branch("ppfjet_cluster_eta", &ppfjet_cluster_eta_);
1437   tree_->Branch("ppfjet_cluster_phi", &ppfjet_cluster_phi_);
1438   tree_->Branch("ppfjet_cluster_dR", &ppfjet_cluster_dR_);
1439   tree_->Branch("ppfjet_ncandtracks", &ppfjet_ncandtracks_, "ppfjet_ncandtracks/I");
1440   tree_->Branch("ppfjet_candtrack_px", &ppfjet_candtrack_px_);
1441   tree_->Branch("ppfjet_candtrack_py", &ppfjet_candtrack_py_);
1442   tree_->Branch("ppfjet_candtrack_pz", &ppfjet_candtrack_pz_);
1443   tree_->Branch("ppfjet_candtrack_EcalE", &ppfjet_candtrack_EcalE_);
1444   tree_->Branch("pf_dijet_deta", &pf_dijet_deta_, "pf_dijet_deta/F");
1445   tree_->Branch("pf_dijet_dphi", &pf_dijet_dphi_, "pf_dijet_dphi/F");
1446   tree_->Branch("pf_dijet_balance", &pf_dijet_balance_, "pf_dijet_balance/F");
1447   tree_->Branch("pf_thirdjet_px", &pf_thirdjet_px_, "pf_thirdjet_px/F");
1448   tree_->Branch("pf_thirdjet_py", &pf_thirdjet_py_, "pf_thirdjet_py/F");
1449   tree_->Branch("pf_realthirdjet_px", &pf_realthirdjet_px_, "pf_realthirdjet_px/F");
1450   tree_->Branch("pf_realthirdjet_py", &pf_realthirdjet_py_, "pf_realthirdjet_py/F");
1451   tree_->Branch("pf_realthirdjet_scale", &pf_realthirdjet_scale_, "pf_realthirdjet_scale/F");
1452   tree_->Branch("pf_Run", &pf_Run_, "pf_Run/I");
1453   tree_->Branch("pf_Lumi", &pf_Lumi_, "pf_Lumi/I");
1454   tree_->Branch("pf_Event", &pf_Event_, "pf_Event/I");
1455   tree_->Branch("pf_NPV", &pf_NPV_, "pf_NPV/I");
1456 
1457   return;
1458 }
1459 
1460 // ------------ method called once each job just after ending the event loop  ------------
1461 void DiJetAnalyzer::endJob() {}
1462 
1463 // helper function
1464 
1465 double DiJetAnalyzer::deltaR(const reco::Jet* j1, const reco::Jet* j2) {
1466   double deta = j1->eta() - j2->eta();
1467   double dphi = std::fabs(j1->phi() - j2->phi());
1468   if (dphi > 3.1415927)
1469     dphi = 2 * 3.1415927 - dphi;
1470   return std::sqrt(deta * deta + dphi * dphi);
1471 }
1472 
1473 double DiJetAnalyzer::deltaR(const double eta1, const double phi1, const double eta2, const double phi2) {
1474   double deta = eta1 - eta2;
1475   double dphi = std::fabs(phi1 - phi2);
1476   if (dphi > 3.1415927)
1477     dphi = 2 * 3.1415927 - dphi;
1478   return std::sqrt(deta * deta + dphi * dphi);
1479 }
1480 
1481 int DiJetAnalyzer::getEtaPhi(const DetId id) {
1482   return id.rawId() & 0x3FFF;  // Get 14 least-significant digits
1483 }
1484 
1485 int DiJetAnalyzer::getEtaPhi(const HcalDetId id) {
1486   return id.rawId() & 0x3FFF;  // Get 14 least-significant digits
1487 }
1488 
1489 //define this as a plug-in
1490 DEFINE_FWK_MODULE(DiJetAnalyzer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!