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