File indexing completed on 2023-01-15 23:49:19
0001
0002 #include <iostream>
0003 #include <map>
0004 #include <memory>
0005 #include <regex>
0006 #include <set>
0007 #include <vector>
0008 #include <string>
0009
0010 #include "TTree.h"
0011 #include "TFile.h"
0012 #include "TClonesArray.h"
0013 #include "TObjString.h"
0014
0015
0016 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0017 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0018 #include "FWCore/Common/interface/TriggerNames.h"
0019 #include "FWCore/Framework/interface/Frameworkfwd.h"
0020 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/Framework/interface/EventSetup.h"
0023 #include "FWCore/Framework/interface/MakerMacros.h"
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "FWCore/ServiceRegistry/interface/Service.h"
0027 #include "FWCore/Utilities/interface/EDGetToken.h"
0028 #include "FWCore/Utilities/interface/EDMException.h"
0029 #include "DataFormats/Common/interface/TriggerResults.h"
0030 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0031 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0032 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0033 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0034 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0035 #include "DataFormats/JetReco/interface/PFJetCollection.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/METReco/interface/METCollection.h"
0040 #include "DataFormats/METReco/interface/PFMET.h"
0041 #include "DataFormats/METReco/interface/PFMETCollection.h"
0042 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0043 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0044 #include "DataFormats/ParticleFlowReco/interface/PFBlockFwd.h"
0045 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0046 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0047 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0048 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
0049 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0050 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0051 #include "DataFormats/TrackReco/interface/Track.h"
0052 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0053 #include "DataFormats/Common/interface/TriggerResults.h"
0054 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0055 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0056 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0057 #include "DataFormats/VertexReco/interface/Vertex.h"
0058 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0059
0060 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0061 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0062 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0063 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0064
0065 #include "HLTrigger/HLTcore/interface/HLTPrescaleProvider.h"
0066 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0067 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0068
0069 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0070
0071
0072
0073
0074
0075 class PhotonPair : protected std::pair<const reco::Photon*, double> {
0076 public:
0077 PhotonPair() {
0078 first = 0;
0079 second = 0.0;
0080 fIdx = -1;
0081 }
0082 PhotonPair(const reco::Photon* ph, double pt, int setIdx = -1) {
0083 first = ph;
0084 second = pt;
0085 fIdx = setIdx;
0086 }
0087 ~PhotonPair() = default;
0088
0089 inline const reco::Photon* photon(void) const { return first; }
0090 inline void photon(const reco::Photon* ph) {
0091 first = ph;
0092 return;
0093 }
0094 inline double pt(void) const { return second; }
0095 inline void pt(double d) {
0096 second = d;
0097 return;
0098 }
0099 void idx(int set_idx) { fIdx = set_idx; };
0100 int idx() const { return fIdx; }
0101 bool isValid() const { return (first != NULL) ? true : false; }
0102
0103 private:
0104 int fIdx;
0105 };
0106
0107 class PFJetCorretPair : protected std::pair<const reco::PFJet*, double> {
0108 public:
0109 PFJetCorretPair() {
0110 first = 0;
0111 second = 1.0;
0112 }
0113 PFJetCorretPair(const reco::PFJet* j, double s) {
0114 first = j;
0115 second = s;
0116 }
0117 ~PFJetCorretPair() = default;
0118
0119 inline const reco::PFJet* jet(void) const { return first; }
0120 inline void jet(const reco::PFJet* j) {
0121 first = j;
0122 return;
0123 }
0124 inline double scale(void) const { return second; }
0125 inline void scale(double d) {
0126 second = d;
0127 return;
0128 }
0129 double scaledEt() const { return first->et() * second; }
0130 bool isValid() const { return (first != NULL) ? true : false; }
0131
0132 private:
0133 };
0134
0135
0136
0137
0138
0139 class GammaJetAnalysis : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0140 public:
0141 explicit GammaJetAnalysis(const edm::ParameterSet&);
0142 ~GammaJetAnalysis() override = default;
0143
0144 float pfEcalIso(const reco::Photon* localPho1,
0145 edm::Handle<reco::PFCandidateCollection> pfHandle,
0146 float dRmax,
0147 float dRVetoBarrel,
0148 float dRVetoEndcap,
0149 float etaStripBarrel,
0150 float etaStripEndcap,
0151 float energyBarrel,
0152 float energyEndcap,
0153 reco::PFCandidate::ParticleType pfToUse);
0154
0155 float pfHcalIso(const reco::Photon* localPho,
0156 edm::Handle<reco::PFCandidateCollection> pfHandle,
0157 float dRmax,
0158 float dRveto,
0159 reco::PFCandidate::ParticleType pfToUse);
0160
0161 std::vector<float> pfTkIsoWithVertex(const reco::Photon* localPho1,
0162 edm::Handle<reco::PFCandidateCollection> pfHandle,
0163 edm::Handle<reco::VertexCollection> vtxHandle,
0164 float dRmax,
0165 float dRvetoBarrel,
0166 float dRvetoEndcap,
0167 float ptMin,
0168 float dzMax,
0169 float dxyMax,
0170 reco::PFCandidate::ParticleType pfToUse);
0171
0172 private:
0173 void beginJob() override;
0174 void analyze(const edm::Event&, const edm::EventSetup&) override;
0175 void endJob() override;
0176 void beginRun(const edm::Run&, const edm::EventSetup&) override;
0177 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0178
0179
0180 const int debug_;
0181 const unsigned int debugEvent;
0182 int debugHLTTrigNames;
0183
0184 const edm::InputTag rhoCollection_;
0185 const edm::InputTag pfType1METColl, pfMETColl;
0186
0187 const std::string photonCollName_;
0188 const std::string pfJetCollName_;
0189 const std::string genJetCollName_;
0190 const std::string genParticleCollName_;
0191 const std::string genEventInfoName_;
0192 const std::string hbheRecHitName_;
0193 const std::string hfRecHitName_;
0194 const std::string hoRecHitName_;
0195 const std::string rootHistFilename_;
0196 const std::string pvCollName_;
0197 const std::string prodProcess_;
0198
0199 const bool allowNoPhoton_;
0200 const double photonPtMin_;
0201 const double photonJetDPhiMin_;
0202 const double jetEtMin_;
0203 const double jet2EtMax_;
0204 const double jet3EtMax_;
0205 std::vector<std::string> photonTrigNamesV_;
0206 std::vector<std::string> jetTrigNamesV_;
0207 const bool writeTriggerPrescale_;
0208
0209 bool doPFJets_;
0210 bool doGenJets_;
0211 int workOnAOD_;
0212 bool ignoreHLT_;
0213
0214
0215 edm::EDGetTokenT<reco::PhotonCollection> tok_Photon_;
0216 edm::EDGetTokenT<reco::PFJetCollection> tok_PFJet_;
0217 edm::EDGetTokenT<std::vector<reco::GenJet>> tok_GenJet_;
0218 edm::EDGetTokenT<std::vector<reco::GenParticle>> tok_GenPart_;
0219 edm::EDGetTokenT<GenEventInfoProduct> tok_GenEvInfo_;
0220 edm::EDGetTokenT<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>> tok_HBHE_;
0221 edm::EDGetTokenT<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>> tok_HF_;
0222 edm::EDGetTokenT<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>> tok_HO_;
0223 edm::EDGetTokenT<edm::ValueMap<Bool_t>> tok_loosePhoton_;
0224 edm::EDGetTokenT<edm::ValueMap<Bool_t>> tok_tightPhoton_;
0225 edm::EDGetTokenT<std::vector<Bool_t>> tok_loosePhotonV_;
0226 edm::EDGetTokenT<std::vector<Bool_t>> tok_tightPhotonV_;
0227 edm::EDGetTokenT<reco::PFCandidateCollection> tok_PFCand_;
0228 edm::EDGetTokenT<reco::VertexCollection> tok_Vertex_;
0229 edm::EDGetTokenT<reco::GsfElectronCollection> tok_GsfElec_;
0230 edm::EDGetTokenT<double> tok_Rho_;
0231 edm::EDGetTokenT<reco::ConversionCollection> tok_Conv_;
0232 edm::EDGetTokenT<reco::BeamSpot> tok_BS_;
0233 edm::EDGetTokenT<std::vector<reco::Vertex>> tok_PV_;
0234 edm::EDGetTokenT<reco::PFMETCollection> tok_PFMET_;
0235 edm::EDGetTokenT<reco::PFMETCollection> tok_PFType1MET_;
0236 edm::EDGetTokenT<edm::TriggerResults> tok_TrigRes_;
0237 edm::EDGetTokenT<reco::JetCorrector> jetCorrectorToken_;
0238 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0239
0240
0241 TTree* misc_tree_;
0242 TTree* calo_tree_;
0243 TTree* pf_tree_;
0244
0245
0246 HLTPrescaleProvider hltPrescaleProvider_;
0247 std::vector<int> photonTrigFired_;
0248 std::vector<double> photonTrigPrescale_;
0249 std::vector<int> jetTrigFired_;
0250 std::vector<double> jetTrigPrescale_;
0251
0252
0253 int runNumber_, lumiBlock_, eventNumber_;
0254 float eventWeight_, eventPtHat_;
0255 int nPhotons_, nGenJets_;
0256 int nPFJets_;
0257 ULong64_t nProcessed_;
0258 int pf_NPV_;
0259
0260
0261 float met_value_, met_phi_, met_sumEt_;
0262 float metType1_value_, metType1_phi_, metType1_sumEt_;
0263
0264
0265 float rho2012_;
0266 float tagPho_pt_, pho_2nd_pt_, tagPho_energy_, tagPho_eta_, tagPho_phi_, tagPho_sieie_;
0267 float tagPho_HoE_, tagPho_r9_, tagPho_EcalIsoDR04_, tagPho_HcalIsoDR04_, tagPho_HcalIsoDR0412_,
0268 tagPho_TrkIsoHollowDR04_, tagPho_pfiso_myphoton03_;
0269 float tagPho_pfiso_myneutral03_;
0270 std::vector<std::vector<float>> tagPho_pfiso_mycharged03;
0271 int tagPho_pixelSeed_;
0272 int tagPho_ConvSafeEleVeto_;
0273 int tagPho_idTight_, tagPho_idLoose_;
0274 float tagPho_genPt_, tagPho_genEnergy_, tagPho_genEta_, tagPho_genPhi_;
0275 float tagPho_genDeltaR_;
0276
0277
0278
0279 float ppfjet_pt_, ppfjet_p_, ppfjet_E_, ppfjet_eta_, ppfjet_phi_, ppfjet_scale_;
0280 float ppfjet_area_, ppfjet_E_NPVcorr_;
0281 float ppfjet_NeutralHadronFrac_, ppfjet_NeutralEMFrac_;
0282 int ppfjet_nConstituents_;
0283 float ppfjet_ChargedHadronFrac_, ppfjet_ChargedMultiplicity_, ppfjet_ChargedEMFrac_;
0284 float ppfjet_gendr_, ppfjet_genpt_, ppfjet_genp_, ppfjet_genE_;
0285 float ppfjet_unkown_E_, ppfjet_unkown_px_, ppfjet_unkown_py_, ppfjet_unkown_pz_, ppfjet_unkown_EcalE_;
0286 float ppfjet_electron_E_, ppfjet_electron_px_, ppfjet_electron_py_, ppfjet_electron_pz_, ppfjet_electron_EcalE_;
0287 float ppfjet_muon_E_, ppfjet_muon_px_, ppfjet_muon_py_, ppfjet_muon_pz_, ppfjet_muon_EcalE_;
0288 float ppfjet_photon_E_, ppfjet_photon_px_, ppfjet_photon_py_, ppfjet_photon_pz_, ppfjet_photon_EcalE_;
0289 int ppfjet_unkown_n_, ppfjet_electron_n_, ppfjet_muon_n_, ppfjet_photon_n_;
0290 int ppfjet_had_n_;
0291 std::vector<float> ppfjet_had_E_, ppfjet_had_px_, ppfjet_had_py_, ppfjet_had_pz_, ppfjet_had_EcalE_,
0292 ppfjet_had_rawHcalE_, ppfjet_had_emf_, ppfjet_had_E_mctruth_;
0293 std::vector<int> ppfjet_had_id_, ppfjet_had_candtrackind_, ppfjet_had_mcpdgId_, ppfjet_had_ntwrs_;
0294 int ppfjet_ntwrs_;
0295 std::vector<int> ppfjet_twr_ieta_, ppfjet_twr_iphi_, ppfjet_twr_depth_, ppfjet_twr_subdet_, ppfjet_twr_candtrackind_,
0296 ppfjet_twr_hadind_, ppfjet_twr_elmttype_, ppfjet_twr_clusterind_;
0297 std::vector<float> ppfjet_twr_hade_, ppfjet_twr_frac_, ppfjet_twr_dR_;
0298 int ppfjet_cluster_n_;
0299 std::vector<float> ppfjet_cluster_eta_, ppfjet_cluster_phi_, ppfjet_cluster_dR_;
0300 int ppfjet_ncandtracks_;
0301 std::vector<float> ppfjet_candtrack_px_, ppfjet_candtrack_py_, ppfjet_candtrack_pz_, ppfjet_candtrack_EcalE_;
0302
0303
0304 float pfjet2_pt_, pfjet2_p_, pfjet2_E_, pfjet2_eta_, pfjet2_phi_, pfjet2_scale_;
0305 float pfjet2_area_, pfjet2_E_NPVcorr_;
0306 float pfjet2_NeutralHadronFrac_, pfjet2_NeutralEMFrac_;
0307 int pfjet2_nConstituents_;
0308 float pfjet2_ChargedHadronFrac_, pfjet2_ChargedMultiplicity_, pfjet2_ChargedEMFrac_;
0309 float pfjet2_gendr_, pfjet2_genpt_, pfjet2_genp_, pfjet2_genE_;
0310 float pfjet2_unkown_E_, pfjet2_unkown_px_, pfjet2_unkown_py_, pfjet2_unkown_pz_, pfjet2_unkown_EcalE_;
0311 float pfjet2_electron_E_, pfjet2_electron_px_, pfjet2_electron_py_, pfjet2_electron_pz_, pfjet2_electron_EcalE_;
0312 float pfjet2_muon_E_, pfjet2_muon_px_, pfjet2_muon_py_, pfjet2_muon_pz_, pfjet2_muon_EcalE_;
0313 float pfjet2_photon_E_, pfjet2_photon_px_, pfjet2_photon_py_, pfjet2_photon_pz_, pfjet2_photon_EcalE_;
0314 int pfjet2_unkown_n_, pfjet2_electron_n_, pfjet2_muon_n_, pfjet2_photon_n_;
0315 int pfjet2_had_n_;
0316 std::vector<float> pfjet2_had_E_, pfjet2_had_px_, pfjet2_had_py_, pfjet2_had_pz_, pfjet2_had_EcalE_,
0317 pfjet2_had_rawHcalE_, pfjet2_had_emf_, pfjet2_had_E_mctruth_;
0318 std::vector<int> pfjet2_had_id_, pfjet2_had_candtrackind_, pfjet2_had_mcpdgId_, pfjet2_had_ntwrs_;
0319 int pfjet2_ntwrs_;
0320 std::vector<int> pfjet2_twr_ieta_, pfjet2_twr_iphi_, pfjet2_twr_depth_, pfjet2_twr_subdet_, pfjet2_twr_candtrackind_,
0321 pfjet2_twr_hadind_, pfjet2_twr_elmttype_, pfjet2_twr_clusterind_;
0322 std::vector<float> pfjet2_twr_hade_, pfjet2_twr_frac_, pfjet2_twr_dR_;
0323 int pfjet2_cluster_n_;
0324 std::vector<float> pfjet2_cluster_eta_, pfjet2_cluster_phi_, pfjet2_cluster_dR_;
0325 int pfjet2_ncandtracks_;
0326 std::vector<float> pfjet2_candtrack_px_, pfjet2_candtrack_py_, pfjet2_candtrack_pz_, pfjet2_candtrack_EcalE_;
0327
0328 float pf_thirdjet_et_;
0329 float pf_thirdjet_pt_, pf_thirdjet_p_, pf_thirdjet_px_, pf_thirdjet_py_;
0330 float pf_thirdjet_E_, pf_thirdjet_eta_, pf_thirdjet_phi_, pf_thirdjet_scale_;
0331
0332
0333 template <class JetPair_type>
0334 float calc_dPhi(const PhotonPair& pho, const JetPair_type& jet) {
0335 if (!pho.isValid() || !jet.isValid())
0336 return 9999.;
0337 float phi1 = pho.photon()->phi();
0338 float phi2 = jet.jet()->phi();
0339 float dphi = fabs(phi1 - phi2);
0340 const float cPi = 4 * atan(1);
0341 while (dphi > cPi)
0342 dphi = fabs(2 * cPi - dphi);
0343 return dphi;
0344 }
0345
0346 double deltaR(const reco::Jet* j1, const reco::Jet* j2);
0347 double deltaR(const double eta1, const double phi1, const double eta2, const double phi2);
0348 int getEtaPhi(const DetId id);
0349 int getEtaPhi(const HcalDetId id);
0350
0351 void clear_leadingPfJetVars();
0352 void copy_leadingPfJetVars_to_pfJet2();
0353
0354 template <class Jet_type>
0355 double deltaR(const PhotonPair& photon, const Jet_type* jet) {
0356 if (!photon.isValid())
0357 return 9999.;
0358 return deltaR(photon.photon()->eta(), photon.photon()->phi(), jet->eta(), jet->phi());
0359 }
0360
0361 struct PFJetCorretPairComp {
0362 inline bool operator()(const PFJetCorretPair& a, const PFJetCorretPair& b) const {
0363 return (a.jet()->pt() * a.scale()) > (b.jet()->pt() * b.scale());
0364 }
0365 };
0366
0367 struct PhotonPairComp {
0368 inline bool operator()(const PhotonPair& a, const PhotonPair& b) const {
0369 return ((a.photon()->pt()) > (b.photon()->pt()));
0370 }
0371 };
0372 };
0373
0374 inline void HERE(const char* msg) {
0375 if (0 && msg)
0376 edm::LogWarning("GammaJetAnalysis") << msg;
0377 }
0378
0379 double getNeutralPVCorr(double eta, int intNPV, double area, bool isMC_) {
0380 double NPV = static_cast<double>(intNPV);
0381 double etaArray[101] = {-5, -4.9, -4.8, -4.7, -4.6, -4.5, -4.4, -4.3, -4.2, -4.1, -4, -3.9, -3.8, -3.7, -3.6,
0382 -3.5, -3.4, -3.3, -3.2, -3.1, -3, -2.9, -2.8, -2.7, -2.6, -2.5, -2.4, -2.3, -2.2, -2.1,
0383 -2, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1, -0.9, -0.8, -0.7, -0.6,
0384 -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
0385 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.1, 2.2, 2.3, 2.4,
0386 2.5, 2.6, 2.7, 2.8, 2.9, 3, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9,
0387 4, 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5};
0388 int ind = -1;
0389 for (int i = 0; i < 100; ++i) {
0390 if (eta > etaArray[i] && eta < etaArray[i + 1]) {
0391 ind = i;
0392 break;
0393 }
0394 }
0395 if (ind < 0)
0396 return 0;
0397 double pt_density;
0398 if (isMC_) {
0399 double p0[100] = {0.08187, 0.096718, 0.11565, 0.12115, 0.12511, 0.12554, 0.13858, 0.14282, 0.14302, 0.15054,
0400 0.14136, 0.14992, 0.13812, 0.13771, 0.13165, 0.12609, 0.12446, 0.11311, 0.13771, 0.16401,
0401 0.20454, 0.27899, 0.34242, 0.43096, 0.50742, 0.59683, 0.66877, 0.68664, 0.69541, 0.66873,
0402 0.64175, 0.61097, 0.58524, 0.5712, 0.55752, 0.54869, 0.31073, 0.22667, 0.55614, 0.55962,
0403 0.54348, 0.53206, 0.51594, 0.49928, 0.49139, 0.48766, 0.49157, 0.49587, 0.50109, 0.5058,
0404 0.51279, 0.51515, 0.51849, 0.52607, 0.52362, 0.52169, 0.53579, 0.54821, 0.56262, 0.58355,
0405 0.58809, 0.57525, 0.52539, 0.53505, 0.52307, 0.52616, 0.52678, 0.53536, 0.55141, 0.58107,
0406 0.60556, 0.62601, 0.60897, 0.59018, 0.49593, 0.40462, 0.32052, 0.24436, 0.18867, 0.12591,
0407 0.095421, 0.090578, 0.078767, 0.11797, 0.14057, 0.14614, 0.15232, 0.14742, 0.15647, 0.14947,
0408 0.15805, 0.14467, 0.14526, 0.14081, 0.1262, 0.12429, 0.11951, 0.11146, 0.095677, 0.083126};
0409 double p1[100] = {0.26831, 0.30901, 0.37017, 0.38747, 0.41547, 0.45237, 0.49963, 0.54074, 0.54949, 0.5937,
0410 0.56904, 0.60766, 0.58042, 0.59823, 0.58535, 0.54594, 0.58403, 0.601, 0.65401, 0.65049,
0411 0.65264, 0.6387, 0.60646, 0.59669, 0.55561, 0.5053, 0.42889, 0.37264, 0.36456, 0.36088,
0412 0.36728, 0.37439, 0.38779, 0.40133, 0.40989, 0.41722, 0.47539, 0.49848, 0.42642, 0.42431,
0413 0.42113, 0.41285, 0.41003, 0.41116, 0.41231, 0.41634, 0.41795, 0.41806, 0.41786, 0.41765,
0414 0.41779, 0.41961, 0.42144, 0.42192, 0.4209, 0.41885, 0.4163, 0.4153, 0.41864, 0.4257,
0415 0.43018, 0.43218, 0.43798, 0.42723, 0.42185, 0.41349, 0.40553, 0.39132, 0.3779, 0.37055,
0416 0.36522, 0.37057, 0.38058, 0.43259, 0.51052, 0.55918, 0.60178, 0.60995, 0.64087, 0.65554,
0417 0.65308, 0.65654, 0.60466, 0.58678, 0.54392, 0.58277, 0.59651, 0.57916, 0.60744, 0.56882,
0418 0.59323, 0.5499, 0.54003, 0.49938, 0.4511, 0.41499, 0.38676, 0.36955, 0.30803, 0.26659};
0419 double p2[100] = {
0420 0.00080918, 0.00083447, 0.0011378, 0.0011221, 0.0013613, 0.0016362, 0.0015854, 0.0019131,
0421 0.0017474, 0.0020078, 0.001856, 0.0020331, 0.0020823, 0.001898, 0.0020096, 0.0016464,
0422 0.0032413, 0.0045615, 0.0054495, 0.0057584, 0.0058982, 0.0058956, 0.0055109, 0.0051433,
0423 0.0042098, 0.0032096, 0.00044089, -0.0003884, -0.0007059, -0.00092769, -0.001116, -0.0010437,
0424 -0.00080318, -0.00044142, 6.7232e-05, 0.00055265, -0.0014486, -0.0020432, 0.0015121, 0.0016343,
0425 0.0015638, 0.0015707, 0.0014403, 0.0012886, 0.0011684, 0.00099089, 0.00091497, 0.00087915,
0426 0.00084703, 0.00084542, 0.00087419, 0.00088013, 0.00090493, 0.00095853, 0.0010389, 0.0011191,
0427 0.0012643, 0.0013833, 0.001474, 0.0015401, 0.0015582, 0.0014265, 0.00087453, 0.00086639,
0428 0.00042986, -5.0257e-06, -0.00053124, -0.00086417, -0.0011228, -0.0011749, -0.0010068, -0.00083012,
0429 -0.00062906, 0.00021515, 0.0028714, 0.0038835, 0.0047212, 0.0051427, 0.0055762, 0.0055872,
0430 0.0054989, 0.0053033, 0.0044519, 0.0032223, 0.0017641, 0.0021165, 0.0019909, 0.0021061,
0431 0.0020322, 0.0018357, 0.0019829, 0.001683, 0.0018553, 0.0015304, 0.0015822, 0.0013119,
0432 0.0010745, 0.0010808, 0.00080678, 0.00079756};
0433 pt_density = p0[ind] + p1[ind] * (NPV - 1) + p2[ind] * (NPV - 1) * (NPV - 1);
0434 } else {
0435 double p0[100] = {0.12523, 0.14896, 0.17696, 0.19376, 0.20038, 0.21353, 0.25069, 0.27089, 0.29124, 0.31947,
0436 0.31781, 0.35453, 0.35424, 0.38159, 0.39453, 0.4003, 0.34798, 0.26303, 0.24824, 0.22857,
0437 0.22609, 0.26793, 0.30096, 0.37637, 0.44461, 0.55692, 0.70328, 0.72458, 0.75065, 0.73569,
0438 0.72485, 0.69933, 0.69804, 0.70775, 0.70965, 0.71439, 0.72189, 0.73691, 0.74847, 0.74968,
0439 0.73467, 0.70115, 0.6732, 0.65971, 0.65724, 0.67751, 0.69569, 0.70905, 0.71815, 0.72119,
0440 0.72128, 0.71645, 0.70588, 0.68958, 0.66978, 0.65959, 0.66889, 0.68713, 0.71063, 0.74283,
0441 0.75153, 0.74733, 0.73335, 0.71346, 0.70168, 0.69445, 0.68841, 0.67761, 0.67654, 0.6957,
0442 0.70276, 0.71057, 0.68176, 0.64651, 0.49156, 0.38366, 0.31375, 0.24127, 0.21395, 0.17783,
0443 0.19026, 0.21486, 0.24689, 0.3434, 0.40184, 0.39876, 0.3873, 0.36462, 0.36337, 0.32777,
0444 0.328, 0.29868, 0.28087, 0.25713, 0.22466, 0.20784, 0.19798, 0.18054, 0.15022, 0.12811};
0445 double p1[100] = {0.26829, 0.30825, 0.37034, 0.38736, 0.41645, 0.45985, 0.51433, 0.56215, 0.5805, 0.63926,
0446 0.62007, 0.67895, 0.66015, 0.68817, 0.67975, 0.64161, 0.70887, 0.74454, 0.80197, 0.78873,
0447 0.77892, 0.74943, 0.70034, 0.6735, 0.60774, 0.53312, 0.42132, 0.36279, 0.3547, 0.35014,
0448 0.35655, 0.3646, 0.37809, 0.38922, 0.39599, 0.40116, 0.40468, 0.40645, 0.40569, 0.4036,
0449 0.39874, 0.39326, 0.39352, 0.39761, 0.40232, 0.40729, 0.41091, 0.41247, 0.413, 0.41283,
0450 0.41289, 0.4134, 0.41322, 0.41185, 0.40769, 0.40193, 0.39707, 0.39254, 0.39274, 0.3989,
0451 0.40474, 0.40758, 0.40788, 0.40667, 0.40433, 0.40013, 0.39371, 0.38154, 0.36723, 0.3583,
0452 0.35148, 0.35556, 0.36172, 0.41073, 0.50629, 0.57068, 0.62972, 0.65188, 0.69954, 0.72967,
0453 0.74333, 0.76148, 0.71418, 0.69062, 0.63065, 0.67117, 0.68278, 0.66028, 0.68147, 0.62494,
0454 0.64452, 0.58685, 0.57076, 0.52387, 0.47132, 0.42637, 0.39554, 0.37989, 0.31825, 0.27969};
0455 double p2[100] = {
0456 -0.0014595, -0.0014618, -0.0011988, -0.00095404, -5.3893e-05, 0.00018901, 0.00012553, 0.0004172,
0457 0.00020229, 0.00051942, 0.00052088, 0.00076727, 0.0010407, 0.0010184, 0.0013442, 0.0011271,
0458 0.0032841, 0.0045259, 0.0051803, 0.0054477, 0.0055691, 0.0056668, 0.0053084, 0.0050978,
0459 0.0042061, 0.003321, 0.00045155, 0.00021376, 0.0001178, -2.6836e-05, -0.00017689, -0.00014723,
0460 0.00016887, 0.00067322, 0.0012952, 0.0019229, 0.0024702, 0.0028854, 0.0031576, 0.003284,
0461 0.0032643, 0.0031061, 0.0028377, 0.0025386, 0.0022583, 0.0020448, 0.001888, 0.0017968,
0462 0.0017286, 0.0016989, 0.0017014, 0.0017302, 0.0017958, 0.0018891, 0.0020609, 0.0022876,
0463 0.0025391, 0.0028109, 0.0030294, 0.0031867, 0.0032068, 0.0030755, 0.0028181, 0.0023893,
0464 0.0018359, 0.0012192, 0.00061654, 0.00016088, -0.00015204, -0.00019503, -3.7236e-05, 0.00016663,
0465 0.00033833, 0.00082988, 0.0034005, 0.0042941, 0.0050884, 0.0052612, 0.0055901, 0.0054357,
0466 0.0052671, 0.0049174, 0.0042236, 0.0031138, 0.0011733, 0.0014057, 0.0010843, 0.0010992,
0467 0.0007966, 0.00052196, 0.00053029, 0.00021273, 0.00041664, 0.00010455, 0.00015173, -9.7827e-05,
0468 -0.0010859, -0.0013748, -0.0016641, -0.0016887};
0469 pt_density = p0[ind] + p1[ind] * (NPV - 1) + p2[ind] * (NPV - 1) * (NPV - 1);
0470 }
0471 double ECorr = pt_density * area * cosh(eta);
0472 return ECorr;
0473 }
0474
0475
0476
0477 inline unsigned int helper_findTrigger(const std::vector<std::string>& list, const std::string& name) {
0478 std::regex re(std::string("^(") + name + "|" + name + "_v\\d*)$");
0479 for (unsigned int i = 0, n = list.size(); i < n; ++i) {
0480 if (std::regex_match(list[i], re))
0481 return i;
0482 }
0483 return list.size();
0484 }
0485
0486
0487
0488 GammaJetAnalysis::GammaJetAnalysis(const edm::ParameterSet& iConfig)
0489 : debug_(iConfig.getUntrackedParameter<int>("debug", 0)),
0490 debugEvent(iConfig.getUntrackedParameter<unsigned int>("debugEvent", 0)),
0491 debugHLTTrigNames(iConfig.getUntrackedParameter<int>("debugHLTTrigNames", 1)),
0492 rhoCollection_(iConfig.getParameter<edm::InputTag>("rhoColl")),
0493 pfType1METColl(iConfig.getParameter<edm::InputTag>("PFMETTYPE1Coll")),
0494 pfMETColl(iConfig.getParameter<edm::InputTag>("PFMETColl")),
0495 photonCollName_(iConfig.getParameter<std::string>("photonCollName")),
0496 pfJetCollName_(iConfig.getParameter<std::string>("pfJetCollName")),
0497 genJetCollName_(iConfig.getParameter<std::string>("genJetCollName")),
0498 genParticleCollName_(iConfig.getParameter<std::string>("genParticleCollName")),
0499 genEventInfoName_(iConfig.getParameter<std::string>("genEventInfoName")),
0500 hbheRecHitName_(iConfig.getParameter<std::string>("hbheRecHitName")),
0501 hfRecHitName_(iConfig.getParameter<std::string>("hfRecHitName")),
0502 hoRecHitName_(iConfig.getParameter<std::string>("hoRecHitName")),
0503 rootHistFilename_(iConfig.getParameter<std::string>("rootHistFilename")),
0504 pvCollName_(iConfig.getParameter<std::string>("pvCollName")),
0505 prodProcess_((iConfig.exists("prodProcess")) ? iConfig.getUntrackedParameter<std::string>("prodProcess")
0506 : "MYGAMMA"),
0507 allowNoPhoton_(iConfig.getParameter<bool>("allowNoPhoton")),
0508 photonPtMin_(iConfig.getParameter<double>("photonPtMin")),
0509 photonJetDPhiMin_(iConfig.getParameter<double>("photonJetDPhiMin")),
0510 jetEtMin_(iConfig.getParameter<double>("jetEtMin")),
0511 jet2EtMax_(iConfig.getParameter<double>("jet2EtMax")),
0512 jet3EtMax_(iConfig.getParameter<double>("jet3EtMax")),
0513 photonTrigNamesV_(iConfig.getParameter<std::vector<std::string>>("photonTriggers")),
0514 jetTrigNamesV_(iConfig.getParameter<std::vector<std::string>>("jetTriggers")),
0515 writeTriggerPrescale_(iConfig.getParameter<bool>("writeTriggerPrescale")),
0516 doPFJets_(iConfig.getParameter<bool>("doPFJets")),
0517 doGenJets_(iConfig.getParameter<bool>("doGenJets")),
0518 workOnAOD_(iConfig.getParameter<int>("workOnAOD")),
0519 ignoreHLT_(iConfig.getUntrackedParameter<bool>("ignoreHLT", false)),
0520 jetCorrectorToken_(consumes<reco::JetCorrector>(iConfig.getParameter<edm::InputTag>("JetCorrections"))),
0521 tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0522 hltPrescaleProvider_(iConfig, consumesCollector(), *this) {
0523 usesResource(TFileService::kSharedResource);
0524
0525
0526 eventWeight_ = 1.0;
0527 eventPtHat_ = 0.;
0528 nProcessed_ = 0;
0529
0530
0531
0532 if (workOnAOD_ < 2) {
0533 tok_Photon_ = consumes<reco::PhotonCollection>(photonCollName_);
0534 tok_PFJet_ = consumes<reco::PFJetCollection>(pfJetCollName_);
0535 tok_GenJet_ = consumes<std::vector<reco::GenJet>>(genJetCollName_);
0536 tok_GenPart_ = consumes<std::vector<reco::GenParticle>>(genParticleCollName_);
0537 tok_GenEvInfo_ = consumes<GenEventInfoProduct>(genEventInfoName_);
0538 tok_HBHE_ = consumes<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>>(hbheRecHitName_);
0539 tok_HF_ = consumes<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>>(hfRecHitName_);
0540 tok_HO_ = consumes<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>>(hoRecHitName_);
0541 tok_loosePhoton_ = consumes<edm::ValueMap<Bool_t>>(edm::InputTag("PhotonIDProdGED", "PhotonCutBasedIDLoose"));
0542 tok_tightPhoton_ = consumes<edm::ValueMap<Bool_t>>(edm::InputTag("PhotonIDProdGED:PhotonCutBasedIDTight"));
0543 tok_PFCand_ = consumes<reco::PFCandidateCollection>(edm::InputTag("particleFlow"));
0544 tok_Vertex_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0545 tok_GsfElec_ = consumes<reco::GsfElectronCollection>(edm::InputTag("gsfElectrons"));
0546 tok_Rho_ = consumes<double>(rhoCollection_);
0547 tok_Conv_ = consumes<reco::ConversionCollection>(edm::InputTag("allConversions"));
0548 tok_BS_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
0549 tok_PV_ = consumes<std::vector<reco::Vertex>>(pvCollName_);
0550 tok_PFMET_ = consumes<reco::PFMETCollection>(pfMETColl);
0551 tok_PFType1MET_ = consumes<reco::PFMETCollection>(pfType1METColl);
0552 tok_TrigRes_ = consumes<edm::TriggerResults>(edm::InputTag("TriggerResults::HLT"));
0553
0554 } else {
0555
0556 const char* prod = "GammaJetProd";
0557 if (prodProcess_.size() == 0) {
0558 edm::LogError("GammaJetAnalysis") << "prodProcess needs to be defined";
0559 throw edm::Exception(edm::errors::ProductNotFound);
0560 }
0561 const char* an = prodProcess_.c_str();
0562 edm::LogWarning("GammaJetAnalysis") << "FAST FIX: changing " << photonCollName_ << " to"
0563 << edm::InputTag(prod, photonCollName_, an);
0564 tok_Photon_ = consumes<reco::PhotonCollection>(edm::InputTag(prod, photonCollName_, an));
0565 tok_PFJet_ = consumes<reco::PFJetCollection>(edm::InputTag(prod, pfJetCollName_, an));
0566 tok_GenJet_ = consumes<std::vector<reco::GenJet>>(edm::InputTag(prod, genJetCollName_, an));
0567 tok_GenPart_ = consumes<std::vector<reco::GenParticle>>(edm::InputTag(prod, genParticleCollName_, an));
0568 tok_GenEvInfo_ = consumes<GenEventInfoProduct>(edm::InputTag(prod, genEventInfoName_, an));
0569 tok_HBHE_ = consumes<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>>(
0570 edm::InputTag(prod, hbheRecHitName_, an));
0571 tok_HF_ = consumes<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>>(
0572 edm::InputTag(prod, hfRecHitName_, an));
0573 tok_HO_ = consumes<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>>(
0574 edm::InputTag(prod, hoRecHitName_, an));
0575
0576
0577 tok_loosePhotonV_ = consumes<std::vector<Bool_t>>(edm::InputTag(prod, "PhotonIDProdGED:PhotonCutBasedIDLoose", an));
0578 tok_tightPhotonV_ = consumes<std::vector<Bool_t>>(edm::InputTag(prod, "PhotonIDProdGED:PhotonCutBasedIDTight", an));
0579 tok_PFCand_ = consumes<reco::PFCandidateCollection>(edm::InputTag(prod, "particleFlow", an));
0580 tok_Vertex_ = consumes<reco::VertexCollection>(edm::InputTag(prod, "offlinePrimaryVertices", an));
0581 tok_GsfElec_ = consumes<reco::GsfElectronCollection>(edm::InputTag(prod, "gedGsfElectrons", an));
0582 tok_Rho_ = consumes<double>(edm::InputTag(prod, rhoCollection_.label(), an));
0583 tok_Conv_ = consumes<reco::ConversionCollection>(edm::InputTag(prod, "allConversions", an));
0584 tok_BS_ = consumes<reco::BeamSpot>(edm::InputTag(prod, "offlineBeamSpot", an));
0585 tok_PV_ = consumes<std::vector<reco::Vertex>>(edm::InputTag(prod, pvCollName_, an));
0586 tok_PFMET_ = consumes<reco::PFMETCollection>(edm::InputTag(prod, pfMETColl.label(), an));
0587 tok_PFType1MET_ = consumes<reco::PFMETCollection>(edm::InputTag(prod, pfType1METColl.label(), an));
0588 TString HLTlabel = "TriggerResults::HLT";
0589 if (prodProcess_.find("reRECO") != std::string::npos)
0590 HLTlabel.ReplaceAll("HLT", "reHLT");
0591 tok_TrigRes_ = consumes<edm::TriggerResults>(edm::InputTag(prod, HLTlabel.Data(), an));
0592 }
0593 }
0594
0595
0596
0597
0598
0599
0600 void GammaJetAnalysis::analyze(const edm::Event& iEvent, const edm::EventSetup& evSetup) {
0601 nProcessed_++;
0602
0603 edm::LogVerbatim("GammaJetAnalysis") << "nProcessed=" << nProcessed_ << "\n";
0604
0605
0606 const edm::Handle<reco::PhotonCollection> photons = iEvent.getHandle(tok_Photon_);
0607 if (!photons.isValid()) {
0608 edm::LogWarning("GammaJetAnalysis") << "Could not find PhotonCollection named " << photonCollName_;
0609 return;
0610 }
0611
0612 if ((photons->size() == 0) && !allowNoPhoton_) {
0613 if (debug_ > 0)
0614 edm::LogVerbatim("GammaJetAnalysis") << "No photons in the event";
0615 return;
0616 }
0617
0618 nPhotons_ = photons->size();
0619 edm::LogVerbatim("GammaJetAnalysis") << "nPhotons_=" << nPhotons_;
0620
0621
0622 edm::Handle<edm::ValueMap<Bool_t>> loosePhotonQual, tightPhotonQual;
0623 edm::Handle<std::vector<Bool_t>> loosePhotonQual_Vec, tightPhotonQual_Vec;
0624 if (workOnAOD_ < 2) {
0625 iEvent.getByToken(tok_loosePhoton_, loosePhotonQual);
0626 iEvent.getByToken(tok_tightPhoton_, tightPhotonQual);
0627 if (!loosePhotonQual.isValid() || !tightPhotonQual.isValid()) {
0628 edm::LogWarning("GammaJetAnalysis") << "Failed to get photon quality flags";
0629 return;
0630 }
0631 } else {
0632 iEvent.getByToken(tok_loosePhotonV_, loosePhotonQual_Vec);
0633 iEvent.getByToken(tok_tightPhotonV_, tightPhotonQual_Vec);
0634 if (!loosePhotonQual_Vec.isValid() || !tightPhotonQual_Vec.isValid()) {
0635 edm::LogWarning("GammaJetAnalysis") << "Failed to get photon quality flags (vec)";
0636 return;
0637 }
0638 }
0639
0640
0641
0642 std::set<PhotonPair, PhotonPairComp> photonpairset;
0643 int counter = 0;
0644 for (reco::PhotonCollection::const_iterator it = photons->begin(); it != photons->end(); ++it) {
0645
0646 const reco::Photon* photon = &(*it);
0647
0648 photonpairset.insert(PhotonPair(photon, photon->pt(), counter));
0649 counter++;
0650
0651 }
0652
0653 HERE(Form("photonpairset.size=%d", int(photonpairset.size())));
0654
0655 if ((photonpairset.size() == 0) && !allowNoPhoton_) {
0656 if (debug_ > 0)
0657 edm::LogVerbatim("GammaJetAnalysis") << "No good quality photons in the event";
0658 return;
0659 }
0660
0661
0662
0663
0664
0665
0666 PhotonPair photon_tag;
0667 PhotonPair photon_2nd;
0668 counter = 0;
0669 for (std::set<PhotonPair, PhotonPairComp>::const_iterator it = photonpairset.begin(); it != photonpairset.end();
0670 ++it) {
0671 PhotonPair photon = (*it);
0672 ++counter;
0673 if (counter == 1)
0674 photon_tag = photon;
0675 else if (counter == 2)
0676 photon_2nd = photon;
0677 else
0678 break;
0679 }
0680
0681 if (counter == 0) {
0682 edm::LogWarning("GammaJetAnalysis") << "Code bug";
0683 return;
0684 }
0685
0686 HERE(Form("counter=%d", counter));
0687
0688
0689 if (photon_tag.isValid() && (photon_tag.pt() < photonPtMin_)) {
0690 if (debug_ > 0)
0691 edm::LogVerbatim("GammaJetAnalysis") << "largest photonPt=" << photon_tag.pt();
0692 return;
0693 }
0694
0695 HERE("aa");
0696
0697
0698 edm::Handle<reco::PFJetCollection> pfjets;
0699 nPFJets_ = 0;
0700 nGenJets_ = 0;
0701
0702 unsigned int anyJetCount = 0;
0703
0704 if (doPFJets_) {
0705 iEvent.getByToken(tok_PFJet_, pfjets);
0706 if (!pfjets.isValid()) {
0707 edm::LogWarning("GammaJetAnalysis") << "Could not find PFJetCollection named " << pfJetCollName_;
0708 return;
0709 }
0710 anyJetCount += pfjets->size();
0711 nPFJets_ = pfjets->size();
0712 }
0713
0714 HERE(Form("anyJetCount=%d", anyJetCount));
0715
0716 if (anyJetCount == 0) {
0717 if (debug_ > 0)
0718 edm::LogVerbatim("GammaJetAnalysis") << "Event contains no jets";
0719 return;
0720 }
0721 if (debug_ > 0)
0722 edm::LogVerbatim("GammaJetAnalysis") << "nPhotons=" << nPhotons_ << ", nPFJets=" << nPFJets_;
0723
0724 HERE(Form("nPhotons_=%d, nPFJets_=%d", nPhotons_, nPFJets_));
0725
0726
0727 photonTrigFired_.clear();
0728 photonTrigPrescale_.clear();
0729 jetTrigFired_.clear();
0730 jetTrigPrescale_.clear();
0731
0732 HERE("trigger");
0733
0734
0735
0736 bool photonTrigFlag = (photonTrigNamesV_.size() == 0) ? true : false;
0737 bool jetTrigFlag = (jetTrigNamesV_.size() == 0) ? true : false;
0738 if ((photonTrigNamesV_.size() == 1) && (photonTrigNamesV_[0].length() == 0))
0739 photonTrigFlag = true;
0740 if ((jetTrigNamesV_.size() == 1) && (jetTrigNamesV_[0].length() == 0))
0741 jetTrigFlag = true;
0742
0743
0744 if (!photonTrigFlag || !jetTrigFlag) {
0745
0746 edm::Handle<edm::TriggerResults> triggerResults;
0747 if (!iEvent.getByToken(tok_TrigRes_, triggerResults)) {
0748 edm::LogWarning("GammaJetAnalysis") << "Could not find TriggerResults::HLT";
0749 return;
0750 }
0751 const edm::TriggerNames& evTrigNames = iEvent.triggerNames(*triggerResults);
0752
0753 if (debugHLTTrigNames > 0) {
0754 if (debug_ > 1)
0755 edm::LogVerbatim("GammaJetAnalysis") << "debugHLTTrigNames is on";
0756 const std::vector<std::string>* trNames = &evTrigNames.triggerNames();
0757 for (size_t i = 0; i < trNames->size(); ++i) {
0758 if (trNames->at(i).find("_Photon") != std::string::npos) {
0759 if (debug_ > 1)
0760 edm::LogVerbatim("GammaJetAnalysis") << " - " << trNames->at(i);
0761 }
0762 }
0763 if (debug_ > 1)
0764 edm::LogVerbatim("GammaJetAnalysis") << " ";
0765 debugHLTTrigNames--;
0766 }
0767
0768 size_t id = 0;
0769 for (size_t i = 0; i < photonTrigNamesV_.size(); ++i) {
0770 const std::string trigName = photonTrigNamesV_.at(i);
0771 id = helper_findTrigger(evTrigNames.triggerNames(), trigName);
0772 if (id == evTrigNames.size()) {
0773 photonTrigFired_.push_back(0);
0774 photonTrigPrescale_.push_back(-1);
0775 continue;
0776 }
0777 int fired = triggerResults->accept(id);
0778 if (fired)
0779 photonTrigFlag = true;
0780 photonTrigFired_.push_back(fired);
0781 if (!writeTriggerPrescale_)
0782 photonTrigPrescale_.push_back(-1);
0783 else {
0784
0785 auto const prescaleVals =
0786 hltPrescaleProvider_.prescaleValues<double>(iEvent, evSetup, evTrigNames.triggerName(id));
0787 photonTrigPrescale_.push_back(prescaleVals.first * prescaleVals.second);
0788 }
0789 }
0790 for (size_t i = 0; i < jetTrigNamesV_.size(); ++i) {
0791 const std::string trigName = jetTrigNamesV_.at(i);
0792 id = helper_findTrigger(evTrigNames.triggerNames(), trigName);
0793 if (id == evTrigNames.size()) {
0794 jetTrigFired_.push_back(0);
0795 jetTrigPrescale_.push_back(-1);
0796 continue;
0797 }
0798 int fired = triggerResults->accept(id);
0799 if (fired)
0800 jetTrigFlag = true;
0801 jetTrigFired_.push_back(fired);
0802 auto const prescaleVals =
0803 hltPrescaleProvider_.prescaleValues<double>(iEvent, evSetup, evTrigNames.triggerName(id));
0804 jetTrigPrescale_.push_back(prescaleVals.first * prescaleVals.second);
0805 }
0806 }
0807
0808 if (!photonTrigFlag && !jetTrigFlag) {
0809 if (debug_ > 0)
0810 edm::LogVerbatim("GammaJetAnalysis") << "no trigger fired";
0811 return;
0812 }
0813
0814 HERE("start isolation");
0815
0816 tagPho_pfiso_mycharged03.clear();
0817
0818 edm::Handle<std::vector<reco::GenJet>> genjets;
0819 edm::Handle<std::vector<reco::GenParticle>> genparticles;
0820 const edm::Handle<reco::PFCandidateCollection> pfHandle = iEvent.getHandle(tok_PFCand_);
0821 const edm::Handle<reco::VertexCollection> vtxHandle = iEvent.getHandle(tok_Vertex_);
0822 const edm::Handle<reco::GsfElectronCollection> gsfElectronHandle = iEvent.getHandle(tok_GsfElec_);
0823 const edm::Handle<double> rhoHandle_2012 = iEvent.getHandle(tok_Rho_);
0824 rho2012_ = *(rhoHandle_2012.product());
0825
0826
0827 const edm::Handle<reco::ConversionCollection> convH = iEvent.getHandle(tok_Conv_);
0828 const edm::Handle<reco::BeamSpot> beamSpotHandle = iEvent.getHandle(tok_BS_);
0829
0830 HERE("doGenJets");
0831
0832 if (doGenJets_) {
0833
0834 iEvent.getByToken(tok_GenJet_, genjets);
0835 if (!genjets.isValid()) {
0836 edm::LogWarning("GammaJetAnalysis") << "Could not find GenJet vector named " << genJetCollName_;
0837 return;
0838 }
0839 nGenJets_ = genjets->size();
0840
0841
0842 iEvent.getByToken(tok_GenPart_, genparticles);
0843 if (!genparticles.isValid()) {
0844 edm::LogWarning("GammaJetAnalysis") << "Could not find GenParticle vector named " << genParticleCollName_;
0845 return;
0846 }
0847
0848
0849 const edm::Handle<GenEventInfoProduct> genEventInfoProduct = iEvent.getHandle(tok_GenEvInfo_);
0850 if (!genEventInfoProduct.isValid()) {
0851 edm::LogWarning("GammaJetAnalysis") << "Could not find GenEventInfoProduct named " << genEventInfoName_;
0852 return;
0853 }
0854 eventWeight_ = genEventInfoProduct->weight();
0855 eventPtHat_ = 0.;
0856 if (genEventInfoProduct->hasBinningValues()) {
0857 eventPtHat_ = genEventInfoProduct->binningValues()[0];
0858 }
0859 }
0860
0861 runNumber_ = iEvent.id().run();
0862 lumiBlock_ = iEvent.id().luminosityBlock();
0863 eventNumber_ = iEvent.id().event();
0864
0865 HERE(Form("runNumber=%d, eventNumber=%d", runNumber_, eventNumber_));
0866
0867
0868 if (!photon_tag.isValid()) {
0869 tagPho_pt_ = -1;
0870 pho_2nd_pt_ = -1;
0871 tagPho_energy_ = -1;
0872 tagPho_eta_ = 0;
0873 tagPho_phi_ = 0;
0874 tagPho_sieie_ = 0;
0875 tagPho_HoE_ = 0;
0876 tagPho_r9_ = 0;
0877 tagPho_EcalIsoDR04_ = 0;
0878 tagPho_HcalIsoDR04_ = 0;
0879 tagPho_HcalIsoDR0412_ = 0;
0880 tagPho_TrkIsoHollowDR04_ = 0;
0881 tagPho_pfiso_myphoton03_ = 0;
0882 tagPho_pfiso_myneutral03_ = 0;
0883 tagPho_pfiso_mycharged03.clear();
0884 tagPho_pixelSeed_ = 0;
0885 tagPho_ConvSafeEleVeto_ = 0;
0886 tagPho_idTight_ = 0;
0887 tagPho_idLoose_ = 0;
0888 tagPho_genPt_ = 0;
0889 tagPho_genEnergy_ = 0;
0890 tagPho_genEta_ = 0;
0891 tagPho_genPhi_ = 0;
0892 tagPho_genDeltaR_ = 0;
0893 } else {
0894 HERE("bb");
0895
0896 tagPho_pt_ = photon_tag.photon()->pt();
0897 pho_2nd_pt_ = (photon_2nd.photon()) ? photon_2nd.photon()->pt() : -1.;
0898 tagPho_energy_ = photon_tag.photon()->energy();
0899 tagPho_eta_ = photon_tag.photon()->eta();
0900 tagPho_phi_ = photon_tag.photon()->phi();
0901 tagPho_sieie_ = photon_tag.photon()->sigmaIetaIeta();
0902 tagPho_HoE_ = photon_tag.photon()->hadTowOverEm();
0903 tagPho_r9_ = photon_tag.photon()->r9();
0904 tagPho_pixelSeed_ = photon_tag.photon()->hasPixelSeed();
0905 tagPho_TrkIsoHollowDR04_ = photon_tag.photon()->trkSumPtHollowConeDR04();
0906 tagPho_EcalIsoDR04_ = photon_tag.photon()->ecalRecHitSumEtConeDR04();
0907 tagPho_HcalIsoDR04_ = photon_tag.photon()->hcalTowerSumEtConeDR04();
0908 tagPho_HcalIsoDR0412_ = photon_tag.photon()->hcalTowerSumEtConeDR04() +
0909 (photon_tag.photon()->hadronicOverEm() - photon_tag.photon()->hadTowOverEm()) *
0910 (photon_tag.photon()->energy() / cosh((photon_tag.photon()->eta())));
0911
0912 HERE("tt");
0913
0914 tagPho_pfiso_myphoton03_ =
0915 pfEcalIso(photon_tag.photon(), pfHandle, 0.3, 0.0, 0.070, 0.015, 0.0, 0.0, 0.0, reco::PFCandidate::gamma);
0916 tagPho_pfiso_myneutral03_ = pfHcalIso(photon_tag.photon(), pfHandle, 0.3, 0.0, reco::PFCandidate::h0);
0917 HERE("calc charged pfiso");
0918 tagPho_pfiso_mycharged03.push_back(pfTkIsoWithVertex(
0919 photon_tag.photon(), pfHandle, vtxHandle, 0.3, 0.02, 0.02, 0.0, 0.2, 0.1, reco::PFCandidate::h));
0920
0921 HERE("got isolation");
0922
0923
0924 tagPho_ConvSafeEleVeto_ = -999;
0925
0926 HERE("get id");
0927 if (workOnAOD_ < 2) {
0928 HERE(Form("workOnAOD_<2. loose photon qual size=%d", int(loosePhotonQual->size())));
0929
0930 edm::Ref<reco::PhotonCollection> photonRef(photons, photon_tag.idx());
0931 HERE(Form("got photon ref, photon_tag.idx()=%d", photon_tag.idx()));
0932 tagPho_idLoose_ = (loosePhotonQual.isValid()) ? (*loosePhotonQual)[photonRef] : -1;
0933 tagPho_idTight_ = (tightPhotonQual.isValid()) ? (*tightPhotonQual)[photonRef] : -1;
0934 } else {
0935 tagPho_idLoose_ = (loosePhotonQual_Vec.isValid()) ? loosePhotonQual_Vec->at(photon_tag.idx()) : -1;
0936 tagPho_idTight_ = (tightPhotonQual_Vec.isValid()) ? tightPhotonQual_Vec->at(photon_tag.idx()) : -1;
0937 }
0938
0939 if (debug_ > 1)
0940 edm::LogVerbatim("GammaJetAnalysis") << "photon tag ID = " << tagPho_idLoose_ << " and " << tagPho_idTight_;
0941
0942 HERE(Form("tagPhoID= %d and %d", tagPho_idLoose_, tagPho_idTight_));
0943
0944 HERE("reset pho gen");
0945
0946 tagPho_genPt_ = 0;
0947 tagPho_genEnergy_ = 0;
0948 tagPho_genEta_ = 0;
0949 tagPho_genPhi_ = 0;
0950 tagPho_genDeltaR_ = 0;
0951 if (doGenJets_) {
0952 tagPho_genDeltaR_ = 9999.;
0953 for (std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin(); itmc != genparticles->end();
0954 itmc++) {
0955 if (itmc->status() == 1 && itmc->pdgId() == 22) {
0956 float dR = deltaR(tagPho_eta_, tagPho_phi_, itmc->eta(), itmc->phi());
0957 if (dR < tagPho_genDeltaR_) {
0958 tagPho_genPt_ = itmc->pt();
0959 tagPho_genEnergy_ = itmc->energy();
0960 tagPho_genEta_ = itmc->eta();
0961 tagPho_genPhi_ = itmc->phi();
0962 tagPho_genDeltaR_ = dR;
0963 }
0964 }
0965 }
0966 }
0967 }
0968
0969 HERE("run over PFJets");
0970
0971
0972
0973 if (doPFJets_ && (nPFJets_ > 0)) {
0974
0975 const edm::Handle<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>> hbhereco =
0976 iEvent.getHandle(tok_HBHE_);
0977 if (!hbhereco.isValid() && !workOnAOD_) {
0978 edm::LogWarning("GammaJetAnalysis") << "Could not find HBHERecHit named " << hbheRecHitName_;
0979 return;
0980 }
0981
0982
0983 const edm::Handle<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>> hfreco =
0984 iEvent.getHandle(tok_HF_);
0985 if (!hfreco.isValid() && !workOnAOD_) {
0986 edm::LogWarning("GammaJetAnalysis") << "Could not find HFRecHit named " << hfRecHitName_;
0987 return;
0988 }
0989
0990
0991 const edm::Handle<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>> horeco =
0992 iEvent.getHandle(tok_HO_);
0993 if (!horeco.isValid() && !workOnAOD_) {
0994 edm::LogWarning("GammaJetAnalysis") << "Could not find HORecHit named " << hoRecHitName_;
0995 return;
0996 }
0997
0998 HERE("get geometry");
0999
1000
1001 const CaloGeometry* geo = &evSetup.getData(tok_geom_);
1002 const HcalGeometry* HBGeom = dynamic_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, 1));
1003 const HcalGeometry* HEGeom = dynamic_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(DetId::Hcal, 2));
1004 const CaloSubdetectorGeometry* HOGeom = geo->getSubdetectorGeometry(DetId::Hcal, 3);
1005 const CaloSubdetectorGeometry* HFGeom = geo->getSubdetectorGeometry(DetId::Hcal, 4);
1006
1007 HERE("work");
1008
1009 int HBHE_n = 0;
1010 if (hbhereco.isValid()) {
1011 for (edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>::const_iterator ith =
1012 hbhereco->begin();
1013 ith != hbhereco->end();
1014 ++ith) {
1015 HBHE_n++;
1016 if (iEvent.id().event() == debugEvent) {
1017 if (debug_ > 1)
1018 edm::LogVerbatim("GammaJetAnalysis") << (*ith).id().ieta() << " " << (*ith).id().iphi();
1019 }
1020 }
1021 }
1022 int HF_n = 0;
1023 if (hfreco.isValid()) {
1024 for (edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith = hfreco->begin();
1025 ith != hfreco->end();
1026 ++ith) {
1027 HF_n++;
1028 if (iEvent.id().event() == debugEvent) {
1029 }
1030 }
1031 }
1032 int HO_n = 0;
1033 if (horeco.isValid()) {
1034 for (edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>::const_iterator ith = horeco->begin();
1035 ith != horeco->end();
1036 ++ith) {
1037 HO_n++;
1038 if (iEvent.id().event() == debugEvent) {
1039 }
1040 }
1041 }
1042
1043 HERE("Get primary vertices");
1044
1045
1046 const edm::Handle<std::vector<reco::Vertex>> pv = iEvent.getHandle(tok_PV_);
1047 if (!pv.isValid()) {
1048 edm::LogWarning("GammaJetAnalysis") << "Could not find Vertex named " << pvCollName_;
1049 return;
1050 }
1051 pf_NPV_ = 0;
1052 for (std::vector<reco::Vertex>::const_iterator it = pv->begin(); it != pv->end(); ++it) {
1053 if (!it->isFake() && it->ndof() > 4)
1054 ++pf_NPV_;
1055 }
1056
1057 HERE("get corrector");
1058
1059 reco::JetCorrector const& correctorPF = iEvent.get(jetCorrectorToken_);
1060
1061
1062 std::set<PFJetCorretPair, PFJetCorretPairComp> pfjetcorretpairset;
1063 for (reco::PFJetCollection::const_iterator it = pfjets->begin(); it != pfjets->end(); ++it) {
1064 const reco::PFJet* jet = &(*it);
1065
1066 if (deltaR(photon_tag, jet) < 0.5)
1067 continue;
1068 double jec = correctorPF.correction(*it);
1069 pfjetcorretpairset.insert(PFJetCorretPair(jet, jec));
1070 }
1071
1072 PFJetCorretPair pfjet_probe;
1073 PFJetCorretPair pf_2ndjet;
1074 PFJetCorretPair pf_3rdjet;
1075 int jet_cntr = 0;
1076 for (std::set<PFJetCorretPair, PFJetCorretPairComp>::const_iterator it = pfjetcorretpairset.begin();
1077 it != pfjetcorretpairset.end();
1078 ++it) {
1079 PFJetCorretPair jet = (*it);
1080 ++jet_cntr;
1081 if (jet_cntr == 1)
1082 pfjet_probe = jet;
1083 else if (jet_cntr == 2)
1084 pf_2ndjet = jet;
1085 else if (jet_cntr == 3)
1086 pf_3rdjet = jet;
1087
1088 }
1089
1090 HERE("reached selection");
1091
1092
1093 int failSelPF = 0;
1094
1095 if (!pfjet_probe.isValid())
1096 failSelPF |= 1;
1097 else {
1098 if (pfjet_probe.scaledEt() < jetEtMin_)
1099 failSelPF |= 2;
1100 if (calc_dPhi(photon_tag, pfjet_probe) < photonJetDPhiMin_)
1101 failSelPF |= 3;
1102 if (deltaR(photon_tag, pfjet_probe.jet()) < 0.5)
1103 failSelPF |= 4;
1104 if (pf_2ndjet.isValid() && (pf_2ndjet.scaledEt() > jet2EtMax_))
1105 failSelPF |= 5;
1106 if (pf_3rdjet.isValid() && (pf_3rdjet.scaledEt() > jet3EtMax_))
1107 failSelPF |= 6;
1108 }
1109
1110 if (!failSelPF) {
1111
1112 if (pf_3rdjet.isValid()) {
1113 pf_thirdjet_et_ = pf_3rdjet.jet()->et();
1114 pf_thirdjet_pt_ = pf_3rdjet.jet()->pt();
1115 pf_thirdjet_p_ = pf_3rdjet.jet()->p();
1116 pf_thirdjet_px_ = pf_3rdjet.jet()->px();
1117 pf_thirdjet_py_ = pf_3rdjet.jet()->py();
1118 pf_thirdjet_E_ = pf_3rdjet.jet()->energy();
1119 pf_thirdjet_eta_ = pf_3rdjet.jet()->eta();
1120 pf_thirdjet_phi_ = pf_3rdjet.jet()->phi();
1121 pf_thirdjet_scale_ = pf_3rdjet.scale();
1122 } else {
1123 pf_thirdjet_et_ = 0;
1124 pf_thirdjet_pt_ = pf_thirdjet_p_ = 0;
1125 pf_thirdjet_px_ = pf_thirdjet_py_ = 0;
1126 pf_thirdjet_E_ = pf_thirdjet_eta_ = pf_thirdjet_phi_ = 0;
1127 pf_thirdjet_scale_ = 0;
1128 }
1129
1130 HERE("fill PF jet");
1131
1132 int ntypes = 0;
1133
1134
1135
1136
1137
1138
1139
1140 PFJetCorretPair pfjet_probe_store = pfjet_probe;
1141 for (int iJet = 2; iJet > 0; iJet--) {
1142
1143 clear_leadingPfJetVars();
1144
1145 if (iJet == 2)
1146 pfjet_probe = pf_2ndjet;
1147 else
1148 pfjet_probe = pfjet_probe_store;
1149
1150 if (!pfjet_probe.jet()) {
1151 if (iJet == 2) {
1152
1153 copy_leadingPfJetVars_to_pfJet2();
1154 } else {
1155 edm::LogWarning("GammaJetAnalysis") << "error in the code: leading pf jet is null";
1156 }
1157 continue;
1158 }
1159
1160 HERE("work further");
1161
1162
1163 std::map<int, std::pair<int, std::set<float>>> ppfjet_rechits;
1164 std::map<float, int> ppfjet_clusters;
1165
1166
1167 ppfjet_pt_ = pfjet_probe.jet()->pt();
1168 ppfjet_p_ = pfjet_probe.jet()->p();
1169 ppfjet_eta_ = pfjet_probe.jet()->eta();
1170 ppfjet_area_ = pfjet_probe.jet()->jetArea();
1171 ppfjet_E_ = pfjet_probe.jet()->energy();
1172 ppfjet_E_NPVcorr_ =
1173 pfjet_probe.jet()->energy() - getNeutralPVCorr(ppfjet_eta_, pf_NPV_, ppfjet_area_, doGenJets_);
1174 ppfjet_phi_ = pfjet_probe.jet()->phi();
1175 ppfjet_NeutralHadronFrac_ = pfjet_probe.jet()->neutralHadronEnergyFraction();
1176 ppfjet_NeutralEMFrac_ = pfjet_probe.jet()->neutralEmEnergyFraction();
1177 ppfjet_nConstituents_ = pfjet_probe.jet()->nConstituents();
1178 ppfjet_ChargedHadronFrac_ = pfjet_probe.jet()->chargedHadronEnergyFraction();
1179 ppfjet_ChargedMultiplicity_ = pfjet_probe.jet()->chargedMultiplicity();
1180 ppfjet_ChargedEMFrac_ = pfjet_probe.jet()->chargedEmEnergyFraction();
1181 ppfjet_scale_ = pfjet_probe.scale();
1182 ppfjet_ntwrs_ = 0;
1183 ppfjet_cluster_n_ = 0;
1184 ppfjet_ncandtracks_ = 0;
1185
1186 HERE("Get PF constituents");
1187
1188
1189 std::vector<reco::PFCandidatePtr> probeconst = pfjet_probe.jet()->getPFConstituents();
1190 HERE(Form("probeconst.size=%d", int(probeconst.size())));
1191 int iPF = 0;
1192 for (std::vector<reco::PFCandidatePtr>::const_iterator it = probeconst.begin(); it != probeconst.end(); ++it) {
1193 bool hasTrack = false;
1194 if (!(*it))
1195 HERE("\tnull probeconst iterator value\n");
1196 reco::PFCandidate::ParticleType candidateType = (*it)->particleId();
1197 iPF++;
1198 HERE(Form("iPF=%d", iPF));
1199
1200
1201 switch (candidateType) {
1202 case reco::PFCandidate::X:
1203 ppfjet_unkown_E_ += (*it)->energy();
1204 ppfjet_unkown_px_ += (*it)->px();
1205 ppfjet_unkown_py_ += (*it)->py();
1206 ppfjet_unkown_pz_ += (*it)->pz();
1207 ppfjet_unkown_EcalE_ += (*it)->ecalEnergy();
1208 ppfjet_unkown_n_++;
1209 continue;
1210 case reco::PFCandidate::h: {
1211 ppfjet_had_E_.push_back((*it)->energy());
1212 ppfjet_had_px_.push_back((*it)->px());
1213 ppfjet_had_py_.push_back((*it)->py());
1214 ppfjet_had_pz_.push_back((*it)->pz());
1215 ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
1216 ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
1217 ppfjet_had_id_.push_back(0);
1218 ppfjet_had_ntwrs_.push_back(0);
1219 ppfjet_had_n_++;
1220
1221 if (doGenJets_) {
1222 float gendr = 99999;
1223 float genE = 0;
1224 int genpdgId = 0;
1225 for (std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin();
1226 itmc != genparticles->end();
1227 itmc++) {
1228 if (itmc->status() == 1 && itmc->pdgId() > 100) {
1229 double dr = deltaR((*it)->eta(), (*it)->phi(), itmc->eta(), itmc->phi());
1230 if (dr < gendr) {
1231 gendr = dr;
1232 genE = itmc->energy();
1233 genpdgId = itmc->pdgId();
1234 }
1235 }
1236 }
1237 ppfjet_had_E_mctruth_.push_back(genE);
1238 ppfjet_had_mcpdgId_.push_back(genpdgId);
1239 }
1240
1241 reco::TrackRef trackRef = (*it)->trackRef();
1242 if (trackRef.isNonnull()) {
1243 reco::Track track = *trackRef;
1244 ppfjet_candtrack_px_.push_back(track.px());
1245 ppfjet_candtrack_py_.push_back(track.py());
1246 ppfjet_candtrack_pz_.push_back(track.pz());
1247 ppfjet_candtrack_EcalE_.push_back((*it)->ecalEnergy());
1248 ppfjet_had_candtrackind_.push_back(ppfjet_ncandtracks_);
1249 hasTrack = true;
1250 ppfjet_ncandtracks_++;
1251 } else {
1252 ppfjet_had_candtrackind_.push_back(-2);
1253 }
1254 } break;
1255 case reco::PFCandidate::e:
1256 ppfjet_electron_E_ += (*it)->energy();
1257 ppfjet_electron_px_ += (*it)->px();
1258 ppfjet_electron_py_ += (*it)->py();
1259 ppfjet_electron_pz_ += (*it)->pz();
1260 ppfjet_electron_EcalE_ += (*it)->ecalEnergy();
1261 ppfjet_electron_n_++;
1262 continue;
1263 case reco::PFCandidate::mu:
1264 ppfjet_muon_E_ += (*it)->energy();
1265 ppfjet_muon_px_ += (*it)->px();
1266 ppfjet_muon_py_ += (*it)->py();
1267 ppfjet_muon_pz_ += (*it)->pz();
1268 ppfjet_muon_EcalE_ += (*it)->ecalEnergy();
1269 ppfjet_muon_n_++;
1270 continue;
1271 case reco::PFCandidate::gamma:
1272 ppfjet_photon_E_ += (*it)->energy();
1273 ppfjet_photon_px_ += (*it)->px();
1274 ppfjet_photon_py_ += (*it)->py();
1275 ppfjet_photon_pz_ += (*it)->pz();
1276 ppfjet_photon_EcalE_ += (*it)->ecalEnergy();
1277 ppfjet_photon_n_++;
1278 continue;
1279 case reco::PFCandidate::h0: {
1280 ppfjet_had_E_.push_back((*it)->energy());
1281 ppfjet_had_px_.push_back((*it)->px());
1282 ppfjet_had_py_.push_back((*it)->py());
1283 ppfjet_had_pz_.push_back((*it)->pz());
1284 ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
1285 ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
1286 ppfjet_had_id_.push_back(1);
1287 ppfjet_had_candtrackind_.push_back(-1);
1288 ppfjet_had_ntwrs_.push_back(0);
1289 ppfjet_had_n_++;
1290
1291 if (doGenJets_) {
1292 float gendr = 99999;
1293 float genE = 0;
1294 int genpdgId = 0;
1295 for (std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin();
1296 itmc != genparticles->end();
1297 itmc++) {
1298 if (itmc->status() == 1 && itmc->pdgId() > 100) {
1299 double dr = deltaR((*it)->eta(), (*it)->phi(), itmc->eta(), itmc->phi());
1300 if (dr < gendr) {
1301 gendr = dr;
1302 genE = itmc->energy();
1303 genpdgId = itmc->pdgId();
1304 }
1305 }
1306 }
1307 ppfjet_had_E_mctruth_.push_back(genE);
1308 ppfjet_had_mcpdgId_.push_back(genpdgId);
1309 }
1310
1311 break;
1312 }
1313 case reco::PFCandidate::h_HF: {
1314 ppfjet_had_E_.push_back((*it)->energy());
1315 ppfjet_had_px_.push_back((*it)->px());
1316 ppfjet_had_py_.push_back((*it)->py());
1317 ppfjet_had_pz_.push_back((*it)->pz());
1318 ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
1319 ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
1320 ppfjet_had_id_.push_back(2);
1321 ppfjet_had_candtrackind_.push_back(-1);
1322 ppfjet_had_ntwrs_.push_back(0);
1323 ppfjet_had_n_++;
1324
1325 if (doGenJets_) {
1326 float gendr = 99999;
1327 float genE = 0;
1328 int genpdgId = 0;
1329 for (std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin();
1330 itmc != genparticles->end();
1331 itmc++) {
1332 if (itmc->status() == 1 && itmc->pdgId() > 100) {
1333 double dr = deltaR((*it)->eta(), (*it)->phi(), itmc->eta(), itmc->phi());
1334 if (dr < gendr) {
1335 gendr = dr;
1336 genE = itmc->energy();
1337 genpdgId = itmc->pdgId();
1338 }
1339 }
1340 }
1341 ppfjet_had_E_mctruth_.push_back(genE);
1342 ppfjet_had_mcpdgId_.push_back(genpdgId);
1343 }
1344
1345 break;
1346 }
1347 case reco::PFCandidate::egamma_HF: {
1348 ppfjet_had_E_.push_back((*it)->energy());
1349 ppfjet_had_px_.push_back((*it)->px());
1350 ppfjet_had_py_.push_back((*it)->py());
1351 ppfjet_had_pz_.push_back((*it)->pz());
1352 ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
1353 ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
1354 ppfjet_had_id_.push_back(3);
1355 ppfjet_had_candtrackind_.push_back(-1);
1356 ppfjet_had_ntwrs_.push_back(0);
1357 ppfjet_had_n_++;
1358
1359 if (doGenJets_) {
1360 float gendr = 99999;
1361 float genE = 0;
1362 int genpdgId = 0;
1363 for (std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin();
1364 itmc != genparticles->end();
1365 itmc++) {
1366 if (itmc->status() == 1 && itmc->pdgId() > 100) {
1367 double dr = deltaR((*it)->eta(), (*it)->phi(), itmc->eta(), itmc->phi());
1368 if (dr < gendr) {
1369 gendr = dr;
1370 genE = itmc->energy();
1371 genpdgId = itmc->pdgId();
1372 }
1373 }
1374 }
1375 ppfjet_had_E_mctruth_.push_back(genE);
1376 ppfjet_had_mcpdgId_.push_back(genpdgId);
1377 }
1378
1379 break;
1380 }
1381 }
1382
1383 float HFHAD_E = 0;
1384 float HFEM_E = 0;
1385 int HFHAD_n_ = 0;
1386 int HFEM_n_ = 0;
1387 int maxElement = (*it)->elementsInBlocks().size();
1388 if (debug_ > 1)
1389 edm::LogVerbatim("GammaJetAnalysis") << "maxElement=" << maxElement;
1390 if (workOnAOD_ == 1) {
1391 maxElement = 0;
1392 if (debug_ > 1)
1393 edm::LogVerbatim("GammaJetAnalysis") << "forced 0";
1394 }
1395 HERE(Form("maxElement=%d", maxElement));
1396 for (int e = 0; e < maxElement; ++e) {
1397
1398 reco::PFBlockRef blockRef = (*it)->elementsInBlocks()[e].first;
1399 const edm::OwnVector<reco::PFBlockElement>& elements = blockRef->elements();
1400 for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
1401 if (elements[iEle].index() == (*it)->elementsInBlocks()[e].second) {
1402 if (elements[iEle].type() == reco::PFBlockElement::HCAL) {
1403
1404 reco::PFClusterRef clusterref = elements[iEle].clusterRef();
1405 reco::PFCluster cluster = *clusterref;
1406 double cluster_dR = deltaR(ppfjet_eta_, ppfjet_phi_, cluster.eta(), cluster.phi());
1407 if (ppfjet_clusters.count(cluster_dR) == 0) {
1408 ppfjet_clusters[cluster_dR] = ppfjet_cluster_n_;
1409 ppfjet_cluster_eta_.push_back(cluster.eta());
1410 ppfjet_cluster_phi_.push_back(cluster.phi());
1411 ppfjet_cluster_dR_.push_back(cluster_dR);
1412 ppfjet_cluster_n_++;
1413 }
1414 int cluster_ind = ppfjet_clusters[cluster_dR];
1415 std::vector<std::pair<DetId, float>> hitsAndFracs = cluster.hitsAndFractions();
1416
1417
1418 int nHits = hitsAndFracs.size();
1419 for (int iHit = 0; iHit < nHits; iHit++) {
1420 int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
1421
1422 for (edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>::const_iterator ith =
1423 hbhereco->begin();
1424 ith != hbhereco->end();
1425 ++ith) {
1426 int etaPhiRecHit = getEtaPhi((*ith).id());
1427 if (etaPhiPF == etaPhiRecHit) {
1428 ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1429 if (ppfjet_rechits.count((*ith).id()) == 0) {
1430 ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1431 ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1432 ppfjet_twr_depth_.push_back((*ith).id().depth());
1433 ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1434 ppfjet_twr_hade_.push_back((*ith).energy());
1435 ppfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
1436 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1437 ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1438 ppfjet_twr_elmttype_.push_back(0);
1439 ppfjet_twr_clusterind_.push_back(cluster_ind);
1440 if (hasTrack) {
1441 ppfjet_twr_candtrackind_.push_back(ppfjet_ncandtracks_ - 1);
1442 } else {
1443 ppfjet_twr_candtrackind_.push_back(-1);
1444 }
1445 switch ((*ith).id().subdet()) {
1446 case HcalSubdetector::HcalBarrel: {
1447 CaloCellGeometry::CornersVec cv = HBGeom->getCorners((*ith).id());
1448 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1449 float avgphi =
1450 (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1451 if ((cv[0].phi() < cv[2].phi()) && (debug_ > 1))
1452 edm::LogVerbatim("GammaJetAnalysis") << "pHB" << cv[0].phi() << " " << cv[2].phi();
1453 if (cv[0].phi() < cv[2].phi())
1454 avgphi = (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) +
1455 static_cast<double>(cv[2].phi())) /
1456 2.0;
1457 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1458 break;
1459 }
1460 case HcalSubdetector::HcalEndcap: {
1461 CaloCellGeometry::CornersVec cv = HEGeom->getCorners((*ith).id());
1462 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1463 float avgphi =
1464 (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1465 if ((cv[0].phi() < cv[2].phi()) && (debug_ > 1))
1466 edm::LogVerbatim("GammaJetAnalysis") << "pHE" << cv[0].phi() << " " << cv[2].phi();
1467 if (cv[0].phi() < cv[2].phi())
1468 avgphi = (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) +
1469 static_cast<double>(cv[2].phi())) /
1470 2.0;
1471 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1472 break;
1473 }
1474 default:
1475 ppfjet_twr_dR_.push_back(-1);
1476 break;
1477 }
1478 ppfjet_rechits[(*ith).id()].first = ppfjet_ntwrs_;
1479 ++ppfjet_ntwrs_;
1480 } else if (ppfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0) {
1481 ppfjet_twr_frac_.at(ppfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
1482 if (cluster_dR <
1483 ppfjet_cluster_dR_.at(ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first))) {
1484 ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first) = cluster_ind;
1485 }
1486 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1487 }
1488 }
1489 }
1490 }
1491 }
1492 else if (elements[iEle].type() == reco::PFBlockElement::HFHAD) {
1493 ntypes++;
1494 HFHAD_n_++;
1495
1496
1497
1498 for (edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith =
1499 hfreco->begin();
1500 ith != hfreco->end();
1501 ++ith) {
1502 if ((*ith).id().depth() == 1)
1503 continue;
1504 auto thisCell = HFGeom->getGeometry((*ith).id());
1505 const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
1506
1507 bool passMatch = false;
1508 if ((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()) {
1509 if ((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi())
1510 passMatch = true;
1511 else if (cv[0].phi() < cv[2].phi()) {
1512 if ((*it)->phi() < cv[0].phi())
1513 passMatch = true;
1514 else if ((*it)->phi() > cv[2].phi())
1515 passMatch = true;
1516 }
1517 }
1518
1519 if (passMatch) {
1520 ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1521 ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1522 ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1523 ppfjet_twr_depth_.push_back((*ith).id().depth());
1524 ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1525 ppfjet_twr_hade_.push_back((*ith).energy());
1526 ppfjet_twr_frac_.push_back(1.0);
1527 ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1528 ppfjet_twr_elmttype_.push_back(1);
1529 ppfjet_twr_clusterind_.push_back(-1);
1530 ppfjet_twr_candtrackind_.push_back(-1);
1531 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1532 float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1533 if ((cv[0].phi() < cv[2].phi()) && (debug_ > 1))
1534 edm::LogVerbatim("GammaJetAnalysis") << "pHFhad" << cv[0].phi() << " " << cv[2].phi();
1535 if (cv[0].phi() < cv[2].phi())
1536 avgphi =
1537 (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) /
1538 2.0;
1539 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1540 ++ppfjet_ntwrs_;
1541 HFHAD_E += (*ith).energy();
1542 }
1543 }
1544 } else if (elements[iEle].type() == reco::PFBlockElement::HFEM) {
1545 ntypes++;
1546 HFEM_n_++;
1547
1548 for (edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith =
1549 hfreco->begin();
1550 ith != hfreco->end();
1551 ++ith) {
1552 if ((*ith).id().depth() == 2)
1553 continue;
1554 auto thisCell = HFGeom->getGeometry((*ith).id());
1555 const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
1556
1557 bool passMatch = false;
1558 if ((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()) {
1559 if ((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi())
1560 passMatch = true;
1561 else if (cv[0].phi() < cv[2].phi()) {
1562 if ((*it)->phi() < cv[0].phi())
1563 passMatch = true;
1564 else if ((*it)->phi() > cv[2].phi())
1565 passMatch = true;
1566 }
1567 }
1568
1569 if (passMatch) {
1570 ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1571 ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1572 ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1573 ppfjet_twr_depth_.push_back((*ith).id().depth());
1574 ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1575 ppfjet_twr_hade_.push_back((*ith).energy());
1576 ppfjet_twr_frac_.push_back(1.0);
1577 ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1578 ppfjet_twr_elmttype_.push_back(2);
1579 ppfjet_twr_clusterind_.push_back(-1);
1580 ppfjet_twr_candtrackind_.push_back(-1);
1581 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1582 float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1583 if ((cv[0].phi() < cv[2].phi()) && (debug_ > 1))
1584 edm::LogVerbatim("GammaJetAnalysis") << "pHFem" << cv[0].phi() << " " << cv[2].phi();
1585 if (cv[0].phi() < cv[2].phi())
1586 avgphi =
1587 (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) /
1588 2.0;
1589 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1590 ++ppfjet_ntwrs_;
1591 HFEM_E += (*ith).energy();
1592 }
1593 }
1594 } else if (elements[iEle].type() == reco::PFBlockElement::HO) {
1595 ntypes++;
1596 reco::PFClusterRef clusterref = elements[iEle].clusterRef();
1597 reco::PFCluster cluster = *clusterref;
1598 double cluster_dR = deltaR(ppfjet_eta_, ppfjet_phi_, cluster.eta(), cluster.phi());
1599 if (ppfjet_clusters.count(cluster_dR) == 0) {
1600 ppfjet_clusters[cluster_dR] = ppfjet_cluster_n_;
1601 ppfjet_cluster_eta_.push_back(cluster.eta());
1602 ppfjet_cluster_phi_.push_back(cluster.phi());
1603 ppfjet_cluster_dR_.push_back(cluster_dR);
1604 ppfjet_cluster_n_++;
1605 }
1606 int cluster_ind = ppfjet_clusters[cluster_dR];
1607
1608 std::vector<std::pair<DetId, float>> hitsAndFracs = cluster.hitsAndFractions();
1609 int nHits = hitsAndFracs.size();
1610 for (int iHit = 0; iHit < nHits; iHit++) {
1611 int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
1612
1613 for (edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>::const_iterator ith =
1614 horeco->begin();
1615 ith != horeco->end();
1616 ++ith) {
1617 int etaPhiRecHit = getEtaPhi((*ith).id());
1618 if (etaPhiPF == etaPhiRecHit) {
1619 ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1620 if (ppfjet_rechits.count((*ith).id()) == 0) {
1621 ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1622 ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1623 ppfjet_twr_depth_.push_back((*ith).id().depth());
1624 ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1625 ppfjet_twr_hade_.push_back((*ith).energy());
1626 ppfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
1627 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1628 ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1629 ppfjet_twr_elmttype_.push_back(3);
1630 ppfjet_twr_clusterind_.push_back(cluster_ind);
1631 if (hasTrack) {
1632 ppfjet_twr_candtrackind_.push_back(ppfjet_ncandtracks_ - 1);
1633 } else {
1634 ppfjet_twr_candtrackind_.push_back(-1);
1635 }
1636 auto thisCell = HOGeom->getGeometry((*ith).id());
1637 const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
1638 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1639 float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1640 if ((cv[0].phi() < cv[2].phi()) && (debug_ > 1))
1641 edm::LogVerbatim("GammaJetAnalysis") << "pHO" << cv[0].phi() << " " << cv[2].phi();
1642 if (cv[0].phi() < cv[2].phi())
1643 avgphi = (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) +
1644 static_cast<double>(cv[2].phi())) /
1645 2.0;
1646
1647 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1648 ppfjet_rechits[(*ith).id()].first = ppfjet_ntwrs_;
1649 ++ppfjet_ntwrs_;
1650 } else if (ppfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0) {
1651 ppfjet_twr_frac_.at(ppfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
1652 if (cluster_dR <
1653 ppfjet_cluster_dR_.at(ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first))) {
1654 ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first) = cluster_ind;
1655 }
1656 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1657 }
1658 }
1659 }
1660 }
1661 }
1662 }
1663 }
1664 }
1665 switch (candidateType) {
1666 case reco::PFCandidate::h_HF:
1667 ppfjet_had_emf_.push_back(HFEM_E / (HFEM_E + HFHAD_E));
1668 break;
1669 case reco::PFCandidate::egamma_HF:
1670 ppfjet_had_emf_.push_back(-1);
1671 break;
1672 default:
1673 ppfjet_had_emf_.push_back(-1);
1674 break;
1675 }
1676 }
1677
1678 if (doGenJets_) {
1679
1680 ppfjet_gendr_ = 99999.;
1681 ppfjet_genpt_ = 0;
1682 ppfjet_genp_ = 0;
1683 for (std::vector<reco::GenJet>::const_iterator it = genjets->begin(); it != genjets->end(); ++it) {
1684 const reco::GenJet* jet = &(*it);
1685 double dr = deltaR(jet, pfjet_probe.jet());
1686 if (dr < ppfjet_gendr_) {
1687 ppfjet_gendr_ = dr;
1688 ppfjet_genpt_ = jet->pt();
1689 ppfjet_genp_ = jet->p();
1690 ppfjet_genE_ = jet->energy();
1691 }
1692 }
1693 }
1694 if (iJet == 2) {
1695 copy_leadingPfJetVars_to_pfJet2();
1696 }
1697 }
1698
1699
1700
1701
1702
1703 const edm::Handle<reco::PFMETCollection> pfmet_h = iEvent.getHandle(tok_PFMET_);
1704 if (!pfmet_h.isValid()) {
1705 edm::LogWarning("GammaJetAnalysis") << " could not find " << pfMETColl;
1706 return;
1707 }
1708 met_value_ = pfmet_h->begin()->et();
1709 met_phi_ = pfmet_h->begin()->phi();
1710 met_sumEt_ = pfmet_h->begin()->sumEt();
1711
1712 const edm::Handle<reco::PFMETCollection> pfmetType1_h = iEvent.getHandle(tok_PFType1MET_);
1713 if (pfmetType1_h.isValid()) {
1714 metType1_value_ = pfmetType1_h->begin()->et();
1715 metType1_phi_ = pfmetType1_h->begin()->phi();
1716 metType1_sumEt_ = pfmetType1_h->begin()->sumEt();
1717 } else {
1718
1719 metType1_value_ = -999.;
1720 metType1_phi_ = -999.;
1721 metType1_sumEt_ = -999.;
1722 }
1723
1724
1725 pf_tree_->Fill();
1726 }
1727 }
1728 return;
1729 }
1730
1731
1732 void GammaJetAnalysis::beginJob() {
1733 edm::Service<TFileService> fs;
1734 if (doPFJets_) {
1735 pf_tree_ = fs->make<TTree>("pf_gammajettree", "tree for gamma+jet balancing using PFJets");
1736 }
1737
1738 for (int iJet = 0; iJet < 2; iJet++) {
1739 bool doJet = doPFJets_;
1740 if (!doJet)
1741 continue;
1742 if (iJet > 0)
1743 continue;
1744 TTree* tree = pf_tree_;
1745
1746
1747 tree->Branch("photonTrig_fired", &photonTrigFired_);
1748 tree->Branch("photonTrig_prescale", &photonTrigPrescale_);
1749 tree->Branch("jetTrig_fired", &jetTrigFired_);
1750 tree->Branch("jetTrig_prescale", &jetTrigPrescale_);
1751
1752
1753 tree->Branch("RunNumber", &runNumber_, "RunNumber/I");
1754 tree->Branch("LumiBlock", &lumiBlock_, "LumiBlock/I");
1755 tree->Branch("EventNumber", &eventNumber_, "EventNumber/I");
1756 tree->Branch("EventWeight", &eventWeight_, "EventWeight/F");
1757 tree->Branch("EventPtHat", &eventPtHat_, "EventPtHat/F");
1758
1759
1760 tree->Branch("rho2012", &rho2012_, "rho2012/F");
1761 tree->Branch("tagPho_pt", &tagPho_pt_, "tagPho_pt/F");
1762 tree->Branch("pho_2nd_pt", &pho_2nd_pt_, "pho_2nd_pt/F");
1763 tree->Branch("tagPho_energy", &tagPho_energy_, "tagPho_energy/F");
1764 tree->Branch("tagPho_eta", &tagPho_eta_, "tagPho_eta/F");
1765 tree->Branch("tagPho_phi", &tagPho_phi_, "tagPho_phi/F");
1766 tree->Branch("tagPho_sieie", &tagPho_sieie_, "tagPho_sieie/F");
1767 tree->Branch("tagPho_HoE", &tagPho_HoE_, "tagPho_HoE/F");
1768 tree->Branch("tagPho_r9", &tagPho_r9_, "tagPho_r9/F");
1769 tree->Branch("tagPho_EcalIsoDR04", &tagPho_EcalIsoDR04_, "tagPho_EcalIsoDR04/F");
1770 tree->Branch("tagPho_HcalIsoDR04", &tagPho_HcalIsoDR04_, "tagPho_HcalIsoDR04/F");
1771 tree->Branch("tagPho_HcalIsoDR0412", &tagPho_HcalIsoDR0412_, "tagPho_HcalIsoDR0412/F");
1772 tree->Branch("tagPho_TrkIsoHollowDR04", &tagPho_TrkIsoHollowDR04_, "tagPho_TrkIsoHollowDR04/F");
1773 tree->Branch("tagPho_pfiso_myphoton03", &tagPho_pfiso_myphoton03_, "tagPho_pfiso_myphoton03/F");
1774 tree->Branch("tagPho_pfiso_myneutral03", &tagPho_pfiso_myneutral03_, "tagPho_pfiso_myneutral03/F");
1775 tree->Branch("tagPho_pfiso_mycharged03", "std::vector<std::vector<float> >", &tagPho_pfiso_mycharged03);
1776 tree->Branch("tagPho_pixelSeed", &tagPho_pixelSeed_, "tagPho_pixelSeed/I");
1777 tree->Branch("tagPho_ConvSafeEleVeto", &tagPho_ConvSafeEleVeto_, "tagPho_ConvSafeEleVeto/I");
1778 tree->Branch("tagPho_idTight", &tagPho_idTight_, "tagPho_idTight/I");
1779 tree->Branch("tagPho_idLoose", &tagPho_idLoose_, "tagPho_idLoose/I");
1780
1781 if (doGenJets_) {
1782 tree->Branch("tagPho_genPt", &tagPho_genPt_, "tagPho_genPt/F");
1783 tree->Branch("tagPho_genEnergy", &tagPho_genEnergy_, "tagPho_genEnergy/F");
1784 tree->Branch("tagPho_genEta", &tagPho_genEta_, "tagPho_genEta/F");
1785 tree->Branch("tagPho_genPhi", &tagPho_genPhi_, "tagPho_genPhi/F");
1786 tree->Branch("tagPho_genDeltaR", &tagPho_genDeltaR_, "tagPho_genDeltaR/F");
1787 }
1788
1789 tree->Branch("nPhotons", &nPhotons_, "nPhotons/I");
1790 tree->Branch("nGenJets", &nGenJets_, "nGenJets/I");
1791 }
1792
1793
1794
1795 if (doPFJets_) {
1796 pf_tree_->Branch("nPFJets", &nPFJets_, "nPFJets/I");
1797
1798
1799 pf_tree_->Branch("ppfjet_pt", &ppfjet_pt_, "ppfjet_pt/F");
1800 pf_tree_->Branch("ppfjet_p", &ppfjet_p_, "ppfjet_p/F");
1801 pf_tree_->Branch("ppfjet_E", &ppfjet_E_, "ppfjet_E/F");
1802 pf_tree_->Branch("ppfjet_E_NPVcorr", &ppfjet_E_NPVcorr_, "ppfjet_E_NPVcorr/F");
1803 pf_tree_->Branch("ppfjet_area", &ppfjet_area_, "ppfjet_area/F");
1804 pf_tree_->Branch("ppfjet_eta", &ppfjet_eta_, "ppfjet_eta/F");
1805 pf_tree_->Branch("ppfjet_phi", &ppfjet_phi_, "ppfjet_phi/F");
1806 pf_tree_->Branch("ppfjet_scale", &ppfjet_scale_, "ppfjet_scale/F");
1807 pf_tree_->Branch("ppfjet_NeutralHadronFrac", &ppfjet_NeutralHadronFrac_, "ppfjet_NeutralHadronFrac/F");
1808 pf_tree_->Branch("ppfjet_NeutralEMFrac", &ppfjet_NeutralEMFrac_, "ppfjet_NeutralEMFrac/F");
1809 pf_tree_->Branch("ppfjet_nConstituents", &ppfjet_nConstituents_, "ppfjet_nConstituents/I");
1810 pf_tree_->Branch("ppfjet_ChargedHadronFrac", &ppfjet_ChargedHadronFrac_, "ppfjet_ChargedHadronFrac/F");
1811 pf_tree_->Branch("ppfjet_ChargedMultiplicity", &ppfjet_ChargedMultiplicity_, "ppfjet_ChargedMultiplicity/F");
1812 pf_tree_->Branch("ppfjet_ChargedEMFrac", &ppfjet_ChargedEMFrac_, "ppfjet_ChargedEMFrac/F");
1813 if (doGenJets_) {
1814 pf_tree_->Branch("ppfjet_genpt", &ppfjet_genpt_, "ppfjet_genpt/F");
1815 pf_tree_->Branch("ppfjet_genp", &ppfjet_genp_, "ppfjet_genp/F");
1816 pf_tree_->Branch("ppfjet_genE", &ppfjet_genE_, "ppfjet_genE/F");
1817 pf_tree_->Branch("ppfjet_gendr", &ppfjet_gendr_, "ppfjet_gendr/F");
1818 }
1819 pf_tree_->Branch("ppfjet_unkown_E", &ppfjet_unkown_E_, "ppfjet_unkown_E/F");
1820 pf_tree_->Branch("ppfjet_electron_E", &ppfjet_electron_E_, "ppfjet_electron_E/F");
1821 pf_tree_->Branch("ppfjet_muon_E", &ppfjet_muon_E_, "ppfjet_muon_E/F");
1822 pf_tree_->Branch("ppfjet_photon_E", &ppfjet_photon_E_, "ppfjet_photon_E/F");
1823 pf_tree_->Branch("ppfjet_unkown_px", &ppfjet_unkown_px_, "ppfjet_unkown_px/F");
1824 pf_tree_->Branch("ppfjet_electron_px", &ppfjet_electron_px_, "ppfjet_electron_px/F");
1825 pf_tree_->Branch("ppfjet_muon_px", &ppfjet_muon_px_, "ppfjet_muon_px/F");
1826 pf_tree_->Branch("ppfjet_photon_px", &ppfjet_photon_px_, "ppfjet_photon_px/F");
1827 pf_tree_->Branch("ppfjet_unkown_py", &ppfjet_unkown_py_, "ppfjet_unkown_py/F");
1828 pf_tree_->Branch("ppfjet_electron_py", &ppfjet_electron_py_, "ppfjet_electron_py/F");
1829 pf_tree_->Branch("ppfjet_muon_py", &ppfjet_muon_py_, "ppfjet_muon_py/F");
1830 pf_tree_->Branch("ppfjet_photon_py", &ppfjet_photon_py_, "ppfjet_photon_py/F");
1831 pf_tree_->Branch("ppfjet_unkown_pz", &ppfjet_unkown_pz_, "ppfjet_unkown_pz/F");
1832 pf_tree_->Branch("ppfjet_electron_pz", &ppfjet_electron_pz_, "ppfjet_electron_pz/F");
1833 pf_tree_->Branch("ppfjet_muon_pz", &ppfjet_muon_pz_, "ppfjet_muon_pz/F");
1834 pf_tree_->Branch("ppfjet_photon_pz", &ppfjet_photon_pz_, "ppfjet_photon_pz/F");
1835 pf_tree_->Branch("ppfjet_unkown_EcalE", &ppfjet_unkown_EcalE_, "ppfjet_unkown_EcalE/F");
1836 pf_tree_->Branch("ppfjet_electron_EcalE", &ppfjet_electron_EcalE_, "ppfjet_electron_EcalE/F");
1837 pf_tree_->Branch("ppfjet_muon_EcalE", &ppfjet_muon_EcalE_, "ppfjet_muon_EcalE/F");
1838 pf_tree_->Branch("ppfjet_photon_EcalE", &ppfjet_photon_EcalE_, "ppfjet_photon_EcalE/F");
1839 pf_tree_->Branch("ppfjet_unkown_n", &ppfjet_unkown_n_, "ppfjet_unkown_n/I");
1840 pf_tree_->Branch("ppfjet_electron_n", &ppfjet_electron_n_, "ppfjet_electron_n/I");
1841 pf_tree_->Branch("ppfjet_muon_n", &ppfjet_muon_n_, "ppfjet_muon_n/I");
1842 pf_tree_->Branch("ppfjet_photon_n", &ppfjet_photon_n_, "ppfjet_photon_n/I");
1843 pf_tree_->Branch("ppfjet_had_n", &ppfjet_had_n_, "ppfjet_had_n/I");
1844 pf_tree_->Branch("ppfjet_had_E", &ppfjet_had_E_);
1845 pf_tree_->Branch("ppfjet_had_px", &ppfjet_had_px_);
1846 pf_tree_->Branch("ppfjet_had_py", &ppfjet_had_py_);
1847 pf_tree_->Branch("ppfjet_had_pz", &ppfjet_had_pz_);
1848 pf_tree_->Branch("ppfjet_had_EcalE", &ppfjet_had_EcalE_);
1849 pf_tree_->Branch("ppfjet_had_rawHcalE", &ppfjet_had_rawHcalE_);
1850 pf_tree_->Branch("ppfjet_had_emf", &ppfjet_had_emf_);
1851 pf_tree_->Branch("ppfjet_had_id", &ppfjet_had_id_);
1852 pf_tree_->Branch("ppfjet_had_candtrackind", &ppfjet_had_candtrackind_);
1853 if (doGenJets_) {
1854 pf_tree_->Branch("ppfjet_had_E_mctruth", &ppfjet_had_E_mctruth_);
1855 pf_tree_->Branch("ppfjet_had_mcpdgId", &ppfjet_had_mcpdgId_);
1856 }
1857 pf_tree_->Branch("ppfjet_had_ntwrs", &ppfjet_had_ntwrs_);
1858 pf_tree_->Branch("ppfjet_ntwrs", &ppfjet_ntwrs_, "ppfjet_ntwrs/I");
1859 pf_tree_->Branch("ppfjet_twr_ieta", &ppfjet_twr_ieta_);
1860 pf_tree_->Branch("ppfjet_twr_iphi", &ppfjet_twr_iphi_);
1861 pf_tree_->Branch("ppfjet_twr_depth", &ppfjet_twr_depth_);
1862 pf_tree_->Branch("ppfjet_twr_subdet", &ppfjet_twr_subdet_);
1863 pf_tree_->Branch("ppfjet_twr_hade", &ppfjet_twr_hade_);
1864 pf_tree_->Branch("ppfjet_twr_frac", &ppfjet_twr_frac_);
1865 pf_tree_->Branch("ppfjet_twr_candtrackind", &ppfjet_twr_candtrackind_);
1866 pf_tree_->Branch("ppfjet_twr_hadind", &ppfjet_twr_hadind_);
1867 pf_tree_->Branch("ppfjet_twr_elmttype", &ppfjet_twr_elmttype_);
1868 pf_tree_->Branch("ppfjet_twr_dR", &ppfjet_twr_dR_);
1869 pf_tree_->Branch("ppfjet_twr_clusterind", &ppfjet_twr_clusterind_);
1870 pf_tree_->Branch("ppfjet_cluster_n", &ppfjet_cluster_n_, "ppfjet_cluster_n/I");
1871 pf_tree_->Branch("ppfjet_cluster_eta", &ppfjet_cluster_eta_);
1872 pf_tree_->Branch("ppfjet_cluster_phi", &ppfjet_cluster_phi_);
1873 pf_tree_->Branch("ppfjet_cluster_dR", &ppfjet_cluster_dR_);
1874 pf_tree_->Branch("ppfjet_ncandtracks", &ppfjet_ncandtracks_, "ppfjet_ncandtracks/I");
1875 pf_tree_->Branch("ppfjet_candtrack_px", &ppfjet_candtrack_px_);
1876 pf_tree_->Branch("ppfjet_candtrack_py", &ppfjet_candtrack_py_);
1877 pf_tree_->Branch("ppfjet_candtrack_pz", &ppfjet_candtrack_pz_);
1878 pf_tree_->Branch("ppfjet_candtrack_EcalE", &ppfjet_candtrack_EcalE_);
1879
1880
1881 pf_tree_->Branch("pfjet2_pt", &pfjet2_pt_, "pfjet2_pt/F");
1882 pf_tree_->Branch("pfjet2_p", &pfjet2_p_, "pfjet2_p/F");
1883 pf_tree_->Branch("pfjet2_E", &pfjet2_E_, "pfjet2_E/F");
1884 pf_tree_->Branch("pfjet2_E_NPVcorr", &pfjet2_E_NPVcorr_, "pfjet2_E_NPVcorr/F");
1885 pf_tree_->Branch("pfjet2_area", &pfjet2_area_, "pfjet2_area/F");
1886 pf_tree_->Branch("pfjet2_eta", &pfjet2_eta_, "pfjet2_eta/F");
1887 pf_tree_->Branch("pfjet2_phi", &pfjet2_phi_, "pfjet2_phi/F");
1888 pf_tree_->Branch("pfjet2_scale", &pfjet2_scale_, "pfjet2_scale/F");
1889 pf_tree_->Branch("pfjet2_NeutralHadronFrac", &pfjet2_NeutralHadronFrac_, "pfjet2_NeutralHadronFrac/F");
1890 pf_tree_->Branch("pfjet2_NeutralEMFrac", &pfjet2_NeutralEMFrac_, "pfjet2_NeutralEMFrac/F");
1891 pf_tree_->Branch("pfjet2_nConstituents", &pfjet2_nConstituents_, "pfjet2_nConstituents/I");
1892 pf_tree_->Branch("pfjet2_ChargedHadronFrac", &pfjet2_ChargedHadronFrac_, "pfjet2_ChargedHadronFrac/F");
1893 pf_tree_->Branch("pfjet2_ChargedMultiplicity", &pfjet2_ChargedMultiplicity_, "pfjet2_ChargedMultiplicity/F");
1894 pf_tree_->Branch("pfjet2_ChargedEMFrac", &pfjet2_ChargedEMFrac_, "pfjet2_ChargedEMFrac/F");
1895 if (doGenJets_) {
1896 pf_tree_->Branch("pfjet2_genpt", &pfjet2_genpt_, "pfjet2_genpt/F");
1897 pf_tree_->Branch("pfjet2_genp", &pfjet2_genp_, "pfjet2_genp/F");
1898 pf_tree_->Branch("pfjet2_genE", &pfjet2_genE_, "pfjet2_genE/F");
1899 pf_tree_->Branch("pfjet2_gendr", &pfjet2_gendr_, "pfjet2_gendr/F");
1900 }
1901 pf_tree_->Branch("pfjet2_unkown_E", &pfjet2_unkown_E_, "pfjet2_unkown_E/F");
1902 pf_tree_->Branch("pfjet2_electron_E", &pfjet2_electron_E_, "pfjet2_electron_E/F");
1903 pf_tree_->Branch("pfjet2_muon_E", &pfjet2_muon_E_, "pfjet2_muon_E/F");
1904 pf_tree_->Branch("pfjet2_photon_E", &pfjet2_photon_E_, "pfjet2_photon_E/F");
1905 pf_tree_->Branch("pfjet2_unkown_px", &pfjet2_unkown_px_, "pfjet2_unkown_px/F");
1906 pf_tree_->Branch("pfjet2_electron_px", &pfjet2_electron_px_, "pfjet2_electron_px/F");
1907 pf_tree_->Branch("pfjet2_muon_px", &pfjet2_muon_px_, "pfjet2_muon_px/F");
1908 pf_tree_->Branch("pfjet2_photon_px", &pfjet2_photon_px_, "pfjet2_photon_px/F");
1909 pf_tree_->Branch("pfjet2_unkown_py", &pfjet2_unkown_py_, "pfjet2_unkown_py/F");
1910 pf_tree_->Branch("pfjet2_electron_py", &pfjet2_electron_py_, "pfjet2_electron_py/F");
1911 pf_tree_->Branch("pfjet2_muon_py", &pfjet2_muon_py_, "pfjet2_muon_py/F");
1912 pf_tree_->Branch("pfjet2_photon_py", &pfjet2_photon_py_, "pfjet2_photon_py/F");
1913 pf_tree_->Branch("pfjet2_unkown_pz", &pfjet2_unkown_pz_, "pfjet2_unkown_pz/F");
1914 pf_tree_->Branch("pfjet2_electron_pz", &pfjet2_electron_pz_, "pfjet2_electron_pz/F");
1915 pf_tree_->Branch("pfjet2_muon_pz", &pfjet2_muon_pz_, "pfjet2_muon_pz/F");
1916 pf_tree_->Branch("pfjet2_photon_pz", &pfjet2_photon_pz_, "pfjet2_photon_pz/F");
1917 pf_tree_->Branch("pfjet2_unkown_EcalE", &pfjet2_unkown_EcalE_, "pfjet2_unkown_EcalE/F");
1918 pf_tree_->Branch("pfjet2_electron_EcalE", &pfjet2_electron_EcalE_, "pfjet2_electron_EcalE/F");
1919 pf_tree_->Branch("pfjet2_muon_EcalE", &pfjet2_muon_EcalE_, "pfjet2_muon_EcalE/F");
1920 pf_tree_->Branch("pfjet2_photon_EcalE", &pfjet2_photon_EcalE_, "pfjet2_photon_EcalE/F");
1921 pf_tree_->Branch("pfjet2_unkown_n", &pfjet2_unkown_n_, "pfjet2_unkown_n/I");
1922 pf_tree_->Branch("pfjet2_electron_n", &pfjet2_electron_n_, "pfjet2_electron_n/I");
1923 pf_tree_->Branch("pfjet2_muon_n", &pfjet2_muon_n_, "pfjet2_muon_n/I");
1924 pf_tree_->Branch("pfjet2_photon_n", &pfjet2_photon_n_, "pfjet2_photon_n/I");
1925 pf_tree_->Branch("pfjet2_had_n", &pfjet2_had_n_, "pfjet2_had_n/I");
1926 pf_tree_->Branch("pfjet2_had_E", &pfjet2_had_E_);
1927 pf_tree_->Branch("pfjet2_had_px", &pfjet2_had_px_);
1928 pf_tree_->Branch("pfjet2_had_py", &pfjet2_had_py_);
1929 pf_tree_->Branch("pfjet2_had_pz", &pfjet2_had_pz_);
1930 pf_tree_->Branch("pfjet2_had_EcalE", &pfjet2_had_EcalE_);
1931 pf_tree_->Branch("pfjet2_had_rawHcalE", &pfjet2_had_rawHcalE_);
1932 pf_tree_->Branch("pfjet2_had_emf", &pfjet2_had_emf_);
1933 pf_tree_->Branch("pfjet2_had_id", &pfjet2_had_id_);
1934 pf_tree_->Branch("pfjet2_had_candtrackind", &pfjet2_had_candtrackind_);
1935 if (doGenJets_) {
1936 pf_tree_->Branch("pfjet2_had_E_mctruth", &pfjet2_had_E_mctruth_);
1937 pf_tree_->Branch("pfjet2_had_mcpdgId", &pfjet2_had_mcpdgId_);
1938 }
1939 pf_tree_->Branch("pfjet2_had_ntwrs", &pfjet2_had_ntwrs_);
1940 pf_tree_->Branch("pfjet2_ntwrs", &pfjet2_ntwrs_, "pfjet2_ntwrs/I");
1941 pf_tree_->Branch("pfjet2_twr_ieta", &pfjet2_twr_ieta_);
1942 pf_tree_->Branch("pfjet2_twr_iphi", &pfjet2_twr_iphi_);
1943 pf_tree_->Branch("pfjet2_twr_depth", &pfjet2_twr_depth_);
1944 pf_tree_->Branch("pfjet2_twr_subdet", &pfjet2_twr_subdet_);
1945 pf_tree_->Branch("pfjet2_twr_hade", &pfjet2_twr_hade_);
1946 pf_tree_->Branch("pfjet2_twr_frac", &pfjet2_twr_frac_);
1947 pf_tree_->Branch("pfjet2_twr_candtrackind", &pfjet2_twr_candtrackind_);
1948 pf_tree_->Branch("pfjet2_twr_hadind", &pfjet2_twr_hadind_);
1949 pf_tree_->Branch("pfjet2_twr_elmttype", &pfjet2_twr_elmttype_);
1950 pf_tree_->Branch("pfjet2_twr_dR", &pfjet2_twr_dR_);
1951 pf_tree_->Branch("pfjet2_twr_clusterind", &pfjet2_twr_clusterind_);
1952 pf_tree_->Branch("pfjet2_cluster_n", &pfjet2_cluster_n_, "pfjet2_cluster_n/I");
1953 pf_tree_->Branch("pfjet2_cluster_eta", &pfjet2_cluster_eta_);
1954 pf_tree_->Branch("pfjet2_cluster_phi", &pfjet2_cluster_phi_);
1955 pf_tree_->Branch("pfjet2_cluster_dR", &pfjet2_cluster_dR_);
1956 pf_tree_->Branch("pfjet2_ncandtracks", &pfjet2_ncandtracks_, "pfjet2_ncandtracks/I");
1957 pf_tree_->Branch("pfjet2_candtrack_px", &pfjet2_candtrack_px_);
1958 pf_tree_->Branch("pfjet2_candtrack_py", &pfjet2_candtrack_py_);
1959 pf_tree_->Branch("pfjet2_candtrack_pz", &pfjet2_candtrack_pz_);
1960 pf_tree_->Branch("pfjet2_candtrack_EcalE", &pfjet2_candtrack_EcalE_);
1961
1962
1963 pf_tree_->Branch("pf_thirdjet_et", &pf_thirdjet_et_, "pf_thirdjet_et/F");
1964 pf_tree_->Branch("pf_thirdjet_pt", &pf_thirdjet_pt_, "pf_thirdjet_pt/F");
1965 pf_tree_->Branch("pf_thirdjet_p", &pf_thirdjet_p_, "pf_thirdjet_p/F");
1966 pf_tree_->Branch("pf_thirdjet_px", &pf_thirdjet_px_, "pf_thirdjet_px/F");
1967 pf_tree_->Branch("pf_thirdjet_py", &pf_thirdjet_py_, "pf_thirdjet_py/F");
1968 pf_tree_->Branch("pf_thirdjet_E", &pf_thirdjet_E_, "pf_thirdjet_E/F");
1969 pf_tree_->Branch("pf_thirdjet_eta", &pf_thirdjet_eta_, "pf_thirdjet_eta/F");
1970 pf_tree_->Branch("pf_thirdjet_phi", &pf_thirdjet_phi_, "pf_thirdjet_phi/F");
1971 pf_tree_->Branch("pf_thirdjet_scale", &pf_thirdjet_scale_, "pf_thirdjet_scale/F");
1972
1973 pf_tree_->Branch("met_value", &met_value_, "met_value/F");
1974 pf_tree_->Branch("met_phi", &met_phi_, "met_phi/F");
1975 pf_tree_->Branch("met_sumEt", &met_sumEt_, "met_sumEt/F");
1976 pf_tree_->Branch("metType1_value", &metType1_value_, "metType1_value/F");
1977 pf_tree_->Branch("metType1_phi", &metType1_phi_, "metType1_phi/F");
1978 pf_tree_->Branch("metType1_sumEt", &metType1_sumEt_, "metType1_sumEt/F");
1979 pf_tree_->Branch("pf_NPV", &pf_NPV_, "pf_NPV/I");
1980 }
1981
1982 return;
1983 }
1984
1985
1986 void GammaJetAnalysis::endJob() {
1987 if (doPFJets_) {
1988 pf_tree_->Write();
1989 }
1990
1991
1992 {
1993 edm::Service<TFileService> fs;
1994 misc_tree_ = fs->make<TTree>("misc_tree", "tree for misc.info");
1995 misc_tree_->Branch("ignoreHLT", &ignoreHLT_, "ignoreHLT/O");
1996 misc_tree_->Branch("doPFJets", &doPFJets_, "doPFJets/O");
1997 misc_tree_->Branch("doGenJets", &doGenJets_, "doGenJets/O");
1998 misc_tree_->Branch("workOnAOD", &workOnAOD_, "workOnAOD/O");
1999 misc_tree_->Branch("photonTriggerNames", &photonTrigNamesV_);
2000 misc_tree_->Branch("jetTriggerNames", &jetTrigNamesV_);
2001 misc_tree_->Branch("nProcessed", &nProcessed_, "nProcessed/l");
2002
2003 time_t ltime;
2004 ltime = time(NULL);
2005 TString str = TString(asctime(localtime(<ime)));
2006 if (str[str.Length() - 1] == '\n')
2007 str.Remove(str.Length() - 1, 1);
2008 TObjString date(str);
2009 date.Write(str.Data());
2010 misc_tree_->Fill();
2011 misc_tree_->Write();
2012 }
2013 }
2014
2015
2016
2017 void GammaJetAnalysis::beginRun(const edm::Run& iRun, const edm::EventSetup& setup) {
2018 if (debug_ > 1)
2019 edm::LogVerbatim("GammaJetAnalysis") << "beginRun()";
2020
2021 if (!ignoreHLT_) {
2022 int noPhotonTrigger = (photonTrigNamesV_.size() == 0) ? 1 : 0;
2023 int noJetTrigger = (jetTrigNamesV_.size() == 0) ? 1 : 0;
2024 if (!noPhotonTrigger && (photonTrigNamesV_.size() == 1) && (photonTrigNamesV_[0].length() == 0))
2025 noPhotonTrigger = 1;
2026 if (!noJetTrigger && (jetTrigNamesV_.size() == 1) && (jetTrigNamesV_[0].length() == 0))
2027 noJetTrigger = 1;
2028 if (noPhotonTrigger && noJetTrigger) {
2029 ignoreHLT_ = true;
2030 if (debug_ > 1)
2031 edm::LogVerbatim("GammaJetAnalysis") << "HLT trigger ignored: no trigger requested";
2032 }
2033 } else {
2034
2035 photonTrigNamesV_.clear();
2036 jetTrigNamesV_.clear();
2037 }
2038
2039 if (!ignoreHLT_) {
2040 if (debug_ > 0)
2041 edm::LogVerbatim("GammaJetAnalysis") << "Initializing trigger information for individual run";
2042 bool changed(true);
2043 std::string processName = "HLT";
2044 if (hltPrescaleProvider_.init(iRun, setup, processName, changed)) {
2045
2046 if (changed) {
2047
2048
2049 }
2050 } else {
2051
2052
2053 throw edm::Exception(edm::errors::ProductNotFound)
2054 << " HLT config extraction failure with process name " << processName;
2055
2056 }
2057 }
2058 }
2059
2060
2061
2062
2063
2064 float GammaJetAnalysis::pfEcalIso(const reco::Photon* localPho1,
2065 edm::Handle<reco::PFCandidateCollection> pfHandle,
2066 float dRmax,
2067 float dRVetoBarrel,
2068 float dRVetoEndcap,
2069 float etaStripBarrel,
2070 float etaStripEndcap,
2071 float energyBarrel,
2072 float energyEndcap,
2073 reco::PFCandidate::ParticleType pfToUse) {
2074 if (debug_ > 1)
2075 edm::LogVerbatim("GammaJetAnalysis") << "Inside pfEcalIso";
2076 reco::Photon* localPho = localPho1->clone();
2077 float dRVeto;
2078 float etaStrip;
2079
2080 if (localPho->isEB()) {
2081 dRVeto = dRVetoBarrel;
2082 etaStrip = etaStripBarrel;
2083 } else {
2084 dRVeto = dRVetoEndcap;
2085 etaStrip = etaStripEndcap;
2086 }
2087 const reco::PFCandidateCollection* forIsolation = pfHandle.product();
2088 int nsize = forIsolation->size();
2089 float sum = 0;
2090 for (int i = 0; i < nsize; i++) {
2091 const reco::PFCandidate& pfc = (*forIsolation)[i];
2092 if (pfc.particleId() == pfToUse) {
2093
2094 if (pfc.superClusterRef().isNonnull() && localPho->superCluster().isNonnull()) {
2095 if (pfc.superClusterRef() == localPho->superCluster())
2096 continue;
2097 }
2098
2099 if (localPho->isEB()) {
2100 if (fabs(pfc.pt()) < energyBarrel)
2101 continue;
2102 } else {
2103 if (fabs(pfc.energy()) < energyEndcap)
2104 continue;
2105 }
2106
2107 math::XYZPoint pfvtx = pfc.vertex();
2108 math::XYZVector photon_directionWrtVtx(localPho->superCluster()->x() - pfvtx.x(),
2109 localPho->superCluster()->y() - pfvtx.y(),
2110 localPho->superCluster()->z() - pfvtx.z());
2111
2112 float dEta = fabs(photon_directionWrtVtx.Eta() - pfc.momentum().Eta());
2113 float dR = deltaR(
2114 photon_directionWrtVtx.Eta(), photon_directionWrtVtx.Phi(), pfc.momentum().Eta(), pfc.momentum().Phi());
2115
2116 if (dEta < etaStrip)
2117 continue;
2118
2119 if (dR > dRmax || dR < dRVeto)
2120 continue;
2121
2122 sum += pfc.pt();
2123 }
2124 }
2125 return sum;
2126 }
2127
2128
2129
2130 float GammaJetAnalysis::pfHcalIso(const reco::Photon* localPho,
2131 edm::Handle<reco::PFCandidateCollection> pfHandle,
2132 float dRmax,
2133 float dRveto,
2134 reco::PFCandidate::ParticleType pfToUse) {
2135 if (debug_ > 1)
2136 edm::LogVerbatim("GammaJetAnalysis") << "Inside pfHcalIso";
2137 return pfEcalIso(localPho, pfHandle, dRmax, dRveto, dRveto, 0.0, 0.0, 0.0, 0.0, pfToUse);
2138 }
2139
2140
2141
2142 std::vector<float> GammaJetAnalysis::pfTkIsoWithVertex(const reco::Photon* localPho1,
2143 edm::Handle<reco::PFCandidateCollection> pfHandle,
2144 edm::Handle<reco::VertexCollection> vtxHandle,
2145 float dRmax,
2146 float dRvetoBarrel,
2147 float dRvetoEndcap,
2148 float ptMin,
2149 float dzMax,
2150 float dxyMax,
2151 reco::PFCandidate::ParticleType pfToUse) {
2152 if (debug_ > 1)
2153 edm::LogVerbatim("GammaJetAnalysis") << "Inside pfTkIsoWithVertex()";
2154 reco::Photon* localPho = localPho1->clone();
2155
2156 float dRveto;
2157 if (localPho->isEB())
2158 dRveto = dRvetoBarrel;
2159 else
2160 dRveto = dRvetoEndcap;
2161
2162 std::vector<float> result;
2163 const reco::PFCandidateCollection* forIsolation = pfHandle.product();
2164
2165
2166 if (debug_ > 1)
2167 edm::LogVerbatim("GammaJetAnalysis") << "vtxHandle->size() = " << vtxHandle->size();
2168 for (unsigned int ivtx = 0; ivtx < (vtxHandle->size()); ++ivtx) {
2169 if (debug_ > 1)
2170 edm::LogVerbatim("GammaJetAnalysis") << "Vtx " << ivtx;
2171
2172 reco::VertexRef vtx(vtxHandle, ivtx);
2173 math::XYZVector photon_directionWrtVtx(localPho->superCluster()->x() - vtx->x(),
2174 localPho->superCluster()->y() - vtx->y(),
2175 localPho->superCluster()->z() - vtx->z());
2176 if (debug_ > 1)
2177 edm::LogVerbatim("GammaJetAnalysis") << "pfTkIsoWithVertex :: Will Loop over the PFCandidates";
2178 float sum = 0;
2179
2180 for (unsigned i = 0; i < forIsolation->size(); i++) {
2181 if (debug_ > 1)
2182 edm::LogVerbatim("GammaJetAnalysis") << "inside loop";
2183 const reco::PFCandidate& pfc = (*forIsolation)[i];
2184
2185
2186 if (debug_ > 1) {
2187 edm::LogVerbatim("GammaJetAnalysis") << "pfToUse=" << pfToUse;
2188 edm::LogVerbatim("GammaJetAnalysis") << "pfc.particleId()=" << pfc.particleId();
2189 }
2190
2191 if (pfc.particleId() == pfToUse) {
2192 if (debug_ > 1) {
2193 edm::LogVerbatim("GammaJetAnalysis") << "\n ***** HERE pfc.particleId() == pfToUse ";
2194 edm::LogVerbatim("GammaJetAnalysis") << "pfc.pt()=" << pfc.pt();
2195 }
2196 if (pfc.pt() < ptMin)
2197 continue;
2198
2199 float dz = fabs(pfc.trackRef()->dz(vtx->position()));
2200 if (dz > dzMax)
2201 continue;
2202
2203 float dxy = fabs(pfc.trackRef()->dxy(vtx->position()));
2204 if (fabs(dxy) > dxyMax)
2205 continue;
2206 float dR = deltaR(
2207 photon_directionWrtVtx.Eta(), photon_directionWrtVtx.Phi(), pfc.momentum().Eta(), pfc.momentum().Phi());
2208 if (dR > dRmax || dR < dRveto)
2209 continue;
2210 sum += pfc.pt();
2211 if (debug_ > 1)
2212 edm::LogVerbatim("GammaJetAnalysis") << "pt=" << pfc.pt();
2213 }
2214 }
2215 if (debug_ > 1)
2216 edm::LogVerbatim("GammaJetAnalysis") << "sum=" << sum;
2217 sum = sum * 1.0;
2218 result.push_back(sum);
2219 }
2220 if (debug_ > 1) {
2221 edm::LogVerbatim("GammaJetAnalysis") << "Will return result";
2222 edm::LogVerbatim("GammaJetAnalysis") << "result" << &result;
2223 edm::LogVerbatim("GammaJetAnalysis") << "Result returned";
2224 }
2225 return result;
2226 }
2227
2228
2229
2230 void GammaJetAnalysis::clear_leadingPfJetVars() {
2231 ppfjet_pt_ = ppfjet_p_ = ppfjet_E_ = 0;
2232 ppfjet_eta_ = ppfjet_phi_ = ppfjet_scale_ = 0.;
2233 ppfjet_area_ = ppfjet_E_NPVcorr_ = 0.;
2234 ppfjet_NeutralHadronFrac_ = ppfjet_NeutralEMFrac_ = 0.;
2235 ppfjet_nConstituents_ = 0;
2236 ppfjet_ChargedHadronFrac_ = ppfjet_ChargedMultiplicity_ = 0;
2237 ppfjet_ChargedEMFrac_ = 0.;
2238 ppfjet_gendr_ = ppfjet_genpt_ = ppfjet_genp_ = ppfjet_genE_ = 0.;
2239
2240 ppfjet_unkown_E_ = ppfjet_unkown_px_ = ppfjet_unkown_py_ = ppfjet_unkown_pz_ = ppfjet_unkown_EcalE_ = 0.0;
2241 ppfjet_electron_E_ = ppfjet_electron_px_ = ppfjet_electron_py_ = ppfjet_electron_pz_ = ppfjet_electron_EcalE_ = 0.0;
2242 ppfjet_muon_E_ = ppfjet_muon_px_ = ppfjet_muon_py_ = ppfjet_muon_pz_ = ppfjet_muon_EcalE_ = 0.0;
2243 ppfjet_photon_E_ = ppfjet_photon_px_ = ppfjet_photon_py_ = ppfjet_photon_pz_ = ppfjet_photon_EcalE_ = 0.0;
2244 ppfjet_unkown_n_ = ppfjet_electron_n_ = ppfjet_muon_n_ = ppfjet_photon_n_ = 0;
2245 ppfjet_had_n_ = 0;
2246 ppfjet_ntwrs_ = 0;
2247 ppfjet_cluster_n_ = 0;
2248 ppfjet_ncandtracks_ = 0;
2249
2250 ppfjet_had_E_.clear();
2251 ppfjet_had_px_.clear();
2252 ppfjet_had_py_.clear();
2253 ppfjet_had_pz_.clear();
2254 ppfjet_had_EcalE_.clear();
2255 ppfjet_had_rawHcalE_.clear();
2256 ppfjet_had_emf_.clear();
2257 ppfjet_had_E_mctruth_.clear();
2258 ppfjet_had_id_.clear();
2259 ppfjet_had_candtrackind_.clear();
2260 ppfjet_had_mcpdgId_.clear();
2261 ppfjet_had_ntwrs_.clear();
2262 ppfjet_twr_ieta_.clear();
2263 ppfjet_twr_iphi_.clear();
2264 ppfjet_twr_depth_.clear();
2265 ppfjet_twr_subdet_.clear();
2266 ppfjet_twr_candtrackind_.clear();
2267 ppfjet_twr_hadind_.clear();
2268 ppfjet_twr_elmttype_.clear();
2269 ppfjet_twr_hade_.clear();
2270 ppfjet_twr_frac_.clear();
2271 ppfjet_twr_dR_.clear();
2272 ppfjet_twr_clusterind_.clear();
2273 ppfjet_cluster_eta_.clear();
2274 ppfjet_cluster_phi_.clear();
2275 ppfjet_cluster_dR_.clear();
2276 ppfjet_candtrack_px_.clear();
2277 ppfjet_candtrack_py_.clear();
2278 ppfjet_candtrack_pz_.clear();
2279 ppfjet_candtrack_EcalE_.clear();
2280 }
2281
2282
2283
2284 void GammaJetAnalysis::copy_leadingPfJetVars_to_pfJet2() {
2285 pfjet2_pt_ = ppfjet_pt_;
2286 pfjet2_p_ = ppfjet_p_;
2287 pfjet2_E_ = ppfjet_E_;
2288 pfjet2_eta_ = ppfjet_eta_;
2289 pfjet2_phi_ = ppfjet_phi_;
2290 pfjet2_scale_ = ppfjet_scale_;
2291 pfjet2_area_ = ppfjet_area_;
2292 pfjet2_E_NPVcorr_ = ppfjet_E_NPVcorr_;
2293 pfjet2_NeutralHadronFrac_ = ppfjet_NeutralHadronFrac_;
2294 pfjet2_NeutralEMFrac_ = ppfjet_NeutralEMFrac_;
2295 pfjet2_nConstituents_ = ppfjet_nConstituents_;
2296 pfjet2_ChargedHadronFrac_ = ppfjet_ChargedHadronFrac_;
2297 pfjet2_ChargedMultiplicity_ = ppfjet_ChargedMultiplicity_;
2298 pfjet2_ChargedEMFrac_ = ppfjet_ChargedEMFrac_;
2299
2300 pfjet2_gendr_ = ppfjet_gendr_;
2301 pfjet2_genpt_ = ppfjet_genpt_;
2302 pfjet2_genp_ = ppfjet_genp_;
2303 pfjet2_genE_ = ppfjet_genE_;
2304
2305 pfjet2_unkown_E_ = ppfjet_unkown_E_;
2306 pfjet2_unkown_px_ = ppfjet_unkown_px_;
2307 pfjet2_unkown_py_ = ppfjet_unkown_py_;
2308 pfjet2_unkown_pz_ = ppfjet_unkown_pz_;
2309 pfjet2_unkown_EcalE_ = ppfjet_unkown_EcalE_;
2310
2311 pfjet2_electron_E_ = ppfjet_electron_E_;
2312 pfjet2_electron_px_ = ppfjet_electron_px_;
2313 pfjet2_electron_py_ = ppfjet_electron_py_;
2314 pfjet2_electron_pz_ = ppfjet_electron_pz_;
2315 pfjet2_electron_EcalE_ = ppfjet_electron_EcalE_;
2316
2317 pfjet2_muon_E_ = ppfjet_muon_E_;
2318 pfjet2_muon_px_ = ppfjet_muon_px_;
2319 pfjet2_muon_py_ = ppfjet_muon_py_;
2320 pfjet2_muon_pz_ = ppfjet_muon_pz_;
2321 pfjet2_muon_EcalE_ = ppfjet_muon_EcalE_;
2322
2323 pfjet2_photon_E_ = ppfjet_photon_E_;
2324 pfjet2_photon_px_ = ppfjet_photon_px_;
2325 pfjet2_photon_py_ = ppfjet_photon_py_;
2326 pfjet2_photon_pz_ = ppfjet_photon_pz_;
2327 pfjet2_photon_EcalE_ = ppfjet_photon_EcalE_;
2328
2329 pfjet2_unkown_n_ = ppfjet_unkown_n_;
2330 pfjet2_electron_n_ = ppfjet_electron_n_;
2331 pfjet2_muon_n_ = ppfjet_muon_n_;
2332 pfjet2_photon_n_ = ppfjet_photon_n_;
2333 pfjet2_had_n_ = ppfjet_had_n_;
2334
2335 pfjet2_had_E_ = ppfjet_had_E_;
2336 pfjet2_had_px_ = ppfjet_had_px_;
2337 pfjet2_had_py_ = ppfjet_had_py_;
2338 pfjet2_had_pz_ = ppfjet_had_pz_;
2339 pfjet2_had_EcalE_ = ppfjet_had_EcalE_;
2340 pfjet2_had_rawHcalE_ = ppfjet_had_rawHcalE_;
2341 pfjet2_had_emf_ = ppfjet_had_emf_;
2342 pfjet2_had_E_mctruth_ = ppfjet_had_E_mctruth_;
2343
2344 pfjet2_had_id_ = ppfjet_had_id_;
2345 pfjet2_had_candtrackind_ = ppfjet_had_candtrackind_;
2346 pfjet2_had_mcpdgId_ = ppfjet_had_mcpdgId_;
2347 pfjet2_had_ntwrs_ = ppfjet_had_ntwrs_;
2348
2349 pfjet2_ntwrs_ = ppfjet_ntwrs_;
2350 pfjet2_twr_ieta_ = ppfjet_twr_ieta_;
2351 pfjet2_twr_iphi_ = ppfjet_twr_iphi_;
2352 pfjet2_twr_depth_ = ppfjet_twr_depth_;
2353 pfjet2_twr_subdet_ = ppfjet_twr_subdet_;
2354 pfjet2_twr_candtrackind_ = ppfjet_twr_candtrackind_;
2355 pfjet2_twr_hadind_ = ppfjet_twr_hadind_;
2356 pfjet2_twr_elmttype_ = ppfjet_twr_elmttype_;
2357 pfjet2_twr_clusterind_ = ppfjet_twr_clusterind_;
2358
2359 pfjet2_twr_hade_ = ppfjet_twr_hade_;
2360 pfjet2_twr_frac_ = ppfjet_twr_frac_;
2361 pfjet2_twr_dR_ = ppfjet_twr_dR_;
2362
2363 pfjet2_cluster_n_ = ppfjet_cluster_n_;
2364 pfjet2_cluster_eta_ = ppfjet_cluster_eta_;
2365 pfjet2_cluster_phi_ = ppfjet_cluster_phi_;
2366 pfjet2_cluster_dR_ = ppfjet_cluster_dR_;
2367
2368 pfjet2_ncandtracks_ = ppfjet_ncandtracks_;
2369 pfjet2_candtrack_px_ = ppfjet_candtrack_px_;
2370 pfjet2_candtrack_py_ = ppfjet_candtrack_py_;
2371 pfjet2_candtrack_pz_ = ppfjet_candtrack_pz_;
2372 pfjet2_candtrack_EcalE_ = ppfjet_candtrack_EcalE_;
2373 }
2374
2375
2376
2377 double GammaJetAnalysis::deltaR(const reco::Jet* j1, const reco::Jet* j2) {
2378 double deta = j1->eta() - j2->eta();
2379 double dphi = std::fabs(j1->phi() - j2->phi());
2380 if (dphi > 3.1415927)
2381 dphi = 2 * 3.1415927 - dphi;
2382 return std::sqrt(deta * deta + dphi * dphi);
2383 }
2384
2385
2386
2387 double GammaJetAnalysis::deltaR(const double eta1, const double phi1, const double eta2, const double phi2) {
2388 double deta = eta1 - eta2;
2389 double dphi = std::fabs(phi1 - phi2);
2390 if (dphi > 3.1415927)
2391 dphi = 2 * 3.1415927 - dphi;
2392 return std::sqrt(deta * deta + dphi * dphi);
2393 }
2394
2395
2396
2397
2398
2399
2400
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410 int GammaJetAnalysis::getEtaPhi(const DetId id) {
2411 return id.rawId() & 0x3FFF;
2412 }
2413
2414
2415
2416 int GammaJetAnalysis::getEtaPhi(const HcalDetId id) {
2417 return id.rawId() & 0x3FFF;
2418 }
2419
2420
2421
2422
2423
2424 DEFINE_FWK_MODULE(GammaJetAnalysis);