File indexing completed on 2022-02-28 01:32:06
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/EDMException.h"
0028 #include "DataFormats/Common/interface/TriggerResults.h"
0029 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0030 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0031 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0032 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0033 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0034 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0035 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
0036 #include "DataFormats/HcalRecHit/interface/HFRecHit.h"
0037 #include "DataFormats/HcalRecHit/interface/HORecHit.h"
0038 #include "DataFormats/METReco/interface/METCollection.h"
0039 #include "DataFormats/METReco/interface/PFMET.h"
0040 #include "DataFormats/METReco/interface/PFMETCollection.h"
0041 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0042 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0043 #include "DataFormats/ParticleFlowReco/interface/PFBlockFwd.h"
0044 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0045 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0046 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0047 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
0048 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0049 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0050 #include "DataFormats/TrackReco/interface/Track.h"
0051 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0052 #include "DataFormats/Common/interface/TriggerResults.h"
0053 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0054 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0055 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0056 #include "DataFormats/VertexReco/interface/Vertex.h"
0057 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0058
0059 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0060 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0061 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0062 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0063
0064 #include "HLTrigger/HLTcore/interface/HLTPrescaleProvider.h"
0065 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0066 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
0067 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.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 pfJetCorrName_;
0190 const std::string genJetCollName_;
0191 const std::string genParticleCollName_;
0192 const std::string genEventInfoName_;
0193 const std::string hbheRecHitName_;
0194 const std::string hfRecHitName_;
0195 const std::string hoRecHitName_;
0196 const std::string rootHistFilename_;
0197 const std::string pvCollName_;
0198 const std::string prodProcess_;
0199
0200 const bool allowNoPhoton_;
0201 const double photonPtMin_;
0202 const double photonJetDPhiMin_;
0203 const double jetEtMin_;
0204 const double jet2EtMax_;
0205 const double jet3EtMax_;
0206 std::vector<std::string> photonTrigNamesV_;
0207 std::vector<std::string> jetTrigNamesV_;
0208 const bool writeTriggerPrescale_;
0209
0210 bool doPFJets_;
0211 bool doGenJets_;
0212 int workOnAOD_;
0213 bool ignoreHLT_;
0214
0215
0216 edm::EDGetTokenT<reco::PhotonCollection> tok_Photon_;
0217 edm::EDGetTokenT<reco::PFJetCollection> tok_PFJet_;
0218 edm::EDGetTokenT<std::vector<reco::GenJet>> tok_GenJet_;
0219 edm::EDGetTokenT<std::vector<reco::GenParticle>> tok_GenPart_;
0220 edm::EDGetTokenT<GenEventInfoProduct> tok_GenEvInfo_;
0221 edm::EDGetTokenT<edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>> tok_HBHE_;
0222 edm::EDGetTokenT<edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>> tok_HF_;
0223 edm::EDGetTokenT<edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>> tok_HO_;
0224 edm::EDGetTokenT<edm::ValueMap<Bool_t>> tok_loosePhoton_;
0225 edm::EDGetTokenT<edm::ValueMap<Bool_t>> tok_tightPhoton_;
0226 edm::EDGetTokenT<std::vector<Bool_t>> tok_loosePhotonV_;
0227 edm::EDGetTokenT<std::vector<Bool_t>> tok_tightPhotonV_;
0228 edm::EDGetTokenT<reco::PFCandidateCollection> tok_PFCand_;
0229 edm::EDGetTokenT<reco::VertexCollection> tok_Vertex_;
0230 edm::EDGetTokenT<reco::GsfElectronCollection> tok_GsfElec_;
0231 edm::EDGetTokenT<double> tok_Rho_;
0232 edm::EDGetTokenT<reco::ConversionCollection> tok_Conv_;
0233 edm::EDGetTokenT<reco::BeamSpot> tok_BS_;
0234 edm::EDGetTokenT<std::vector<reco::Vertex>> tok_PV_;
0235 edm::EDGetTokenT<reco::PFMETCollection> tok_PFMET_;
0236 edm::EDGetTokenT<reco::PFMETCollection> tok_PFType1MET_;
0237 edm::EDGetTokenT<edm::TriggerResults> tok_TrigRes_;
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<int> photonTrigPrescale_;
0249 std::vector<int> jetTrigFired_;
0250 std::vector<int> 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 pfJetCorrName_(iConfig.getParameter<std::string>("pfJetCorrName")),
0498 genJetCollName_(iConfig.getParameter<std::string>("genJetCollName")),
0499 genParticleCollName_(iConfig.getParameter<std::string>("genParticleCollName")),
0500 genEventInfoName_(iConfig.getParameter<std::string>("genEventInfoName")),
0501 hbheRecHitName_(iConfig.getParameter<std::string>("hbheRecHitName")),
0502 hfRecHitName_(iConfig.getParameter<std::string>("hfRecHitName")),
0503 hoRecHitName_(iConfig.getParameter<std::string>("hoRecHitName")),
0504 rootHistFilename_(iConfig.getParameter<std::string>("rootHistFilename")),
0505 pvCollName_(iConfig.getParameter<std::string>("pvCollName")),
0506 prodProcess_((iConfig.exists("prodProcess")) ? iConfig.getUntrackedParameter<std::string>("prodProcess")
0507 : "MYGAMMA"),
0508 allowNoPhoton_(iConfig.getParameter<bool>("allowNoPhoton")),
0509 photonPtMin_(iConfig.getParameter<double>("photonPtMin")),
0510 photonJetDPhiMin_(iConfig.getParameter<double>("photonJetDPhiMin")),
0511 jetEtMin_(iConfig.getParameter<double>("jetEtMin")),
0512 jet2EtMax_(iConfig.getParameter<double>("jet2EtMax")),
0513 jet3EtMax_(iConfig.getParameter<double>("jet3EtMax")),
0514 photonTrigNamesV_(iConfig.getParameter<std::vector<std::string>>("photonTriggers")),
0515 jetTrigNamesV_(iConfig.getParameter<std::vector<std::string>>("jetTriggers")),
0516 writeTriggerPrescale_(iConfig.getParameter<bool>("writeTriggerPrescale")),
0517 doPFJets_(iConfig.getParameter<bool>("doPFJets")),
0518 doGenJets_(iConfig.getParameter<bool>("doGenJets")),
0519 workOnAOD_(iConfig.getParameter<int>("workOnAOD")),
0520 ignoreHLT_(iConfig.getUntrackedParameter<bool>("ignoreHLT", false)),
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 std::pair<int, int> prescaleVals =
0786 hltPrescaleProvider_.prescaleValues(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 std::pair<int, int> prescaleVals =
0803 hltPrescaleProvider_.prescaleValues(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
1060 const JetCorrector* correctorPF = JetCorrector::getJetCorrector(pfJetCorrName_, evSetup);
1061
1062
1063 std::set<PFJetCorretPair, PFJetCorretPairComp> pfjetcorretpairset;
1064 for (reco::PFJetCollection::const_iterator it = pfjets->begin(); it != pfjets->end(); ++it) {
1065 const reco::PFJet* jet = &(*it);
1066
1067 if (deltaR(photon_tag, jet) < 0.5)
1068 continue;
1069 double jec = correctorPF->correction(*it, iEvent, evSetup);
1070 pfjetcorretpairset.insert(PFJetCorretPair(jet, jec));
1071 }
1072
1073 PFJetCorretPair pfjet_probe;
1074 PFJetCorretPair pf_2ndjet;
1075 PFJetCorretPair pf_3rdjet;
1076 int jet_cntr = 0;
1077 for (std::set<PFJetCorretPair, PFJetCorretPairComp>::const_iterator it = pfjetcorretpairset.begin();
1078 it != pfjetcorretpairset.end();
1079 ++it) {
1080 PFJetCorretPair jet = (*it);
1081 ++jet_cntr;
1082 if (jet_cntr == 1)
1083 pfjet_probe = jet;
1084 else if (jet_cntr == 2)
1085 pf_2ndjet = jet;
1086 else if (jet_cntr == 3)
1087 pf_3rdjet = jet;
1088
1089 }
1090
1091 HERE("reached selection");
1092
1093
1094 int failSelPF = 0;
1095
1096 if (!pfjet_probe.isValid())
1097 failSelPF |= 1;
1098 else {
1099 if (pfjet_probe.scaledEt() < jetEtMin_)
1100 failSelPF |= 2;
1101 if (calc_dPhi(photon_tag, pfjet_probe) < photonJetDPhiMin_)
1102 failSelPF |= 3;
1103 if (deltaR(photon_tag, pfjet_probe.jet()) < 0.5)
1104 failSelPF |= 4;
1105 if (pf_2ndjet.isValid() && (pf_2ndjet.scaledEt() > jet2EtMax_))
1106 failSelPF |= 5;
1107 if (pf_3rdjet.isValid() && (pf_3rdjet.scaledEt() > jet3EtMax_))
1108 failSelPF |= 6;
1109 }
1110
1111 if (!failSelPF) {
1112
1113 if (pf_3rdjet.isValid()) {
1114 pf_thirdjet_et_ = pf_3rdjet.jet()->et();
1115 pf_thirdjet_pt_ = pf_3rdjet.jet()->pt();
1116 pf_thirdjet_p_ = pf_3rdjet.jet()->p();
1117 pf_thirdjet_px_ = pf_3rdjet.jet()->px();
1118 pf_thirdjet_py_ = pf_3rdjet.jet()->py();
1119 pf_thirdjet_E_ = pf_3rdjet.jet()->energy();
1120 pf_thirdjet_eta_ = pf_3rdjet.jet()->eta();
1121 pf_thirdjet_phi_ = pf_3rdjet.jet()->phi();
1122 pf_thirdjet_scale_ = pf_3rdjet.scale();
1123 } else {
1124 pf_thirdjet_et_ = 0;
1125 pf_thirdjet_pt_ = pf_thirdjet_p_ = 0;
1126 pf_thirdjet_px_ = pf_thirdjet_py_ = 0;
1127 pf_thirdjet_E_ = pf_thirdjet_eta_ = pf_thirdjet_phi_ = 0;
1128 pf_thirdjet_scale_ = 0;
1129 }
1130
1131 HERE("fill PF jet");
1132
1133 int types = 0;
1134 int ntypes = 0;
1135
1136
1137
1138
1139
1140
1141
1142 PFJetCorretPair pfjet_probe_store = pfjet_probe;
1143 for (int iJet = 2; iJet > 0; iJet--) {
1144
1145 clear_leadingPfJetVars();
1146
1147 if (iJet == 2)
1148 pfjet_probe = pf_2ndjet;
1149 else
1150 pfjet_probe = pfjet_probe_store;
1151
1152 if (!pfjet_probe.jet()) {
1153 if (iJet == 2) {
1154
1155 copy_leadingPfJetVars_to_pfJet2();
1156 } else {
1157 edm::LogWarning("GammaJetAnalysis") << "error in the code: leading pf jet is null";
1158 }
1159 continue;
1160 }
1161
1162 HERE("work further");
1163
1164
1165 std::map<int, std::pair<int, std::set<float>>> ppfjet_rechits;
1166 std::map<float, int> ppfjet_clusters;
1167
1168
1169 ppfjet_pt_ = pfjet_probe.jet()->pt();
1170 ppfjet_p_ = pfjet_probe.jet()->p();
1171 ppfjet_eta_ = pfjet_probe.jet()->eta();
1172 ppfjet_area_ = pfjet_probe.jet()->jetArea();
1173 ppfjet_E_ = pfjet_probe.jet()->energy();
1174 ppfjet_E_NPVcorr_ =
1175 pfjet_probe.jet()->energy() - getNeutralPVCorr(ppfjet_eta_, pf_NPV_, ppfjet_area_, doGenJets_);
1176 ppfjet_phi_ = pfjet_probe.jet()->phi();
1177 ppfjet_NeutralHadronFrac_ = pfjet_probe.jet()->neutralHadronEnergyFraction();
1178 ppfjet_NeutralEMFrac_ = pfjet_probe.jet()->neutralEmEnergyFraction();
1179 ppfjet_nConstituents_ = pfjet_probe.jet()->nConstituents();
1180 ppfjet_ChargedHadronFrac_ = pfjet_probe.jet()->chargedHadronEnergyFraction();
1181 ppfjet_ChargedMultiplicity_ = pfjet_probe.jet()->chargedMultiplicity();
1182 ppfjet_ChargedEMFrac_ = pfjet_probe.jet()->chargedEmEnergyFraction();
1183 ppfjet_scale_ = pfjet_probe.scale();
1184 ppfjet_ntwrs_ = 0;
1185 ppfjet_cluster_n_ = 0;
1186 ppfjet_ncandtracks_ = 0;
1187
1188 HERE("Get PF constituents");
1189
1190
1191 std::vector<reco::PFCandidatePtr> probeconst = pfjet_probe.jet()->getPFConstituents();
1192 HERE(Form("probeconst.size=%d", int(probeconst.size())));
1193 int iPF = 0;
1194 for (std::vector<reco::PFCandidatePtr>::const_iterator it = probeconst.begin(); it != probeconst.end(); ++it) {
1195 bool hasTrack = false;
1196 if (!(*it))
1197 HERE("\tnull probeconst iterator value\n");
1198 reco::PFCandidate::ParticleType candidateType = (*it)->particleId();
1199 iPF++;
1200 HERE(Form("iPF=%d", iPF));
1201
1202
1203 switch (candidateType) {
1204 case reco::PFCandidate::X:
1205 ppfjet_unkown_E_ += (*it)->energy();
1206 ppfjet_unkown_px_ += (*it)->px();
1207 ppfjet_unkown_py_ += (*it)->py();
1208 ppfjet_unkown_pz_ += (*it)->pz();
1209 ppfjet_unkown_EcalE_ += (*it)->ecalEnergy();
1210 ppfjet_unkown_n_++;
1211 continue;
1212 case reco::PFCandidate::h: {
1213 ppfjet_had_E_.push_back((*it)->energy());
1214 ppfjet_had_px_.push_back((*it)->px());
1215 ppfjet_had_py_.push_back((*it)->py());
1216 ppfjet_had_pz_.push_back((*it)->pz());
1217 ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
1218 ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
1219 ppfjet_had_id_.push_back(0);
1220 ppfjet_had_ntwrs_.push_back(0);
1221 ppfjet_had_n_++;
1222
1223 if (doGenJets_) {
1224 float gendr = 99999;
1225 float genE = 0;
1226 int genpdgId = 0;
1227 for (std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin();
1228 itmc != genparticles->end();
1229 itmc++) {
1230 if (itmc->status() == 1 && itmc->pdgId() > 100) {
1231 double dr = deltaR((*it)->eta(), (*it)->phi(), itmc->eta(), itmc->phi());
1232 if (dr < gendr) {
1233 gendr = dr;
1234 genE = itmc->energy();
1235 genpdgId = itmc->pdgId();
1236 }
1237 }
1238 }
1239 ppfjet_had_E_mctruth_.push_back(genE);
1240 ppfjet_had_mcpdgId_.push_back(genpdgId);
1241 }
1242
1243 reco::TrackRef trackRef = (*it)->trackRef();
1244 if (trackRef.isNonnull()) {
1245 reco::Track track = *trackRef;
1246 ppfjet_candtrack_px_.push_back(track.px());
1247 ppfjet_candtrack_py_.push_back(track.py());
1248 ppfjet_candtrack_pz_.push_back(track.pz());
1249 ppfjet_candtrack_EcalE_.push_back((*it)->ecalEnergy());
1250 ppfjet_had_candtrackind_.push_back(ppfjet_ncandtracks_);
1251 hasTrack = true;
1252 ppfjet_ncandtracks_++;
1253 } else {
1254 ppfjet_had_candtrackind_.push_back(-2);
1255 }
1256 } break;
1257 case reco::PFCandidate::e:
1258 ppfjet_electron_E_ += (*it)->energy();
1259 ppfjet_electron_px_ += (*it)->px();
1260 ppfjet_electron_py_ += (*it)->py();
1261 ppfjet_electron_pz_ += (*it)->pz();
1262 ppfjet_electron_EcalE_ += (*it)->ecalEnergy();
1263 ppfjet_electron_n_++;
1264 continue;
1265 case reco::PFCandidate::mu:
1266 ppfjet_muon_E_ += (*it)->energy();
1267 ppfjet_muon_px_ += (*it)->px();
1268 ppfjet_muon_py_ += (*it)->py();
1269 ppfjet_muon_pz_ += (*it)->pz();
1270 ppfjet_muon_EcalE_ += (*it)->ecalEnergy();
1271 ppfjet_muon_n_++;
1272 continue;
1273 case reco::PFCandidate::gamma:
1274 ppfjet_photon_E_ += (*it)->energy();
1275 ppfjet_photon_px_ += (*it)->px();
1276 ppfjet_photon_py_ += (*it)->py();
1277 ppfjet_photon_pz_ += (*it)->pz();
1278 ppfjet_photon_EcalE_ += (*it)->ecalEnergy();
1279 ppfjet_photon_n_++;
1280 continue;
1281 case reco::PFCandidate::h0: {
1282 ppfjet_had_E_.push_back((*it)->energy());
1283 ppfjet_had_px_.push_back((*it)->px());
1284 ppfjet_had_py_.push_back((*it)->py());
1285 ppfjet_had_pz_.push_back((*it)->pz());
1286 ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
1287 ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
1288 ppfjet_had_id_.push_back(1);
1289 ppfjet_had_candtrackind_.push_back(-1);
1290 ppfjet_had_ntwrs_.push_back(0);
1291 ppfjet_had_n_++;
1292
1293 if (doGenJets_) {
1294 float gendr = 99999;
1295 float genE = 0;
1296 int genpdgId = 0;
1297 for (std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin();
1298 itmc != genparticles->end();
1299 itmc++) {
1300 if (itmc->status() == 1 && itmc->pdgId() > 100) {
1301 double dr = deltaR((*it)->eta(), (*it)->phi(), itmc->eta(), itmc->phi());
1302 if (dr < gendr) {
1303 gendr = dr;
1304 genE = itmc->energy();
1305 genpdgId = itmc->pdgId();
1306 }
1307 }
1308 }
1309 ppfjet_had_E_mctruth_.push_back(genE);
1310 ppfjet_had_mcpdgId_.push_back(genpdgId);
1311 }
1312
1313 break;
1314 }
1315 case reco::PFCandidate::h_HF: {
1316 ppfjet_had_E_.push_back((*it)->energy());
1317 ppfjet_had_px_.push_back((*it)->px());
1318 ppfjet_had_py_.push_back((*it)->py());
1319 ppfjet_had_pz_.push_back((*it)->pz());
1320 ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
1321 ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
1322 ppfjet_had_id_.push_back(2);
1323 ppfjet_had_candtrackind_.push_back(-1);
1324 ppfjet_had_ntwrs_.push_back(0);
1325 ppfjet_had_n_++;
1326
1327 if (doGenJets_) {
1328 float gendr = 99999;
1329 float genE = 0;
1330 int genpdgId = 0;
1331 for (std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin();
1332 itmc != genparticles->end();
1333 itmc++) {
1334 if (itmc->status() == 1 && itmc->pdgId() > 100) {
1335 double dr = deltaR((*it)->eta(), (*it)->phi(), itmc->eta(), itmc->phi());
1336 if (dr < gendr) {
1337 gendr = dr;
1338 genE = itmc->energy();
1339 genpdgId = itmc->pdgId();
1340 }
1341 }
1342 }
1343 ppfjet_had_E_mctruth_.push_back(genE);
1344 ppfjet_had_mcpdgId_.push_back(genpdgId);
1345 }
1346
1347 break;
1348 }
1349 case reco::PFCandidate::egamma_HF: {
1350 ppfjet_had_E_.push_back((*it)->energy());
1351 ppfjet_had_px_.push_back((*it)->px());
1352 ppfjet_had_py_.push_back((*it)->py());
1353 ppfjet_had_pz_.push_back((*it)->pz());
1354 ppfjet_had_EcalE_.push_back((*it)->ecalEnergy());
1355 ppfjet_had_rawHcalE_.push_back((*it)->rawHcalEnergy());
1356 ppfjet_had_id_.push_back(3);
1357 ppfjet_had_candtrackind_.push_back(-1);
1358 ppfjet_had_ntwrs_.push_back(0);
1359 ppfjet_had_n_++;
1360
1361 if (doGenJets_) {
1362 float gendr = 99999;
1363 float genE = 0;
1364 int genpdgId = 0;
1365 for (std::vector<reco::GenParticle>::const_iterator itmc = genparticles->begin();
1366 itmc != genparticles->end();
1367 itmc++) {
1368 if (itmc->status() == 1 && itmc->pdgId() > 100) {
1369 double dr = deltaR((*it)->eta(), (*it)->phi(), itmc->eta(), itmc->phi());
1370 if (dr < gendr) {
1371 gendr = dr;
1372 genE = itmc->energy();
1373 genpdgId = itmc->pdgId();
1374 }
1375 }
1376 }
1377 ppfjet_had_E_mctruth_.push_back(genE);
1378 ppfjet_had_mcpdgId_.push_back(genpdgId);
1379 }
1380
1381 break;
1382 }
1383 }
1384
1385 float HFHAD_E = 0;
1386 float HFEM_E = 0;
1387 int HFHAD_n_ = 0;
1388 int HFEM_n_ = 0;
1389 int HF_type_ = 0;
1390 int maxElement = (*it)->elementsInBlocks().size();
1391 if (debug_ > 1)
1392 edm::LogVerbatim("GammaJetAnalysis") << "maxElement=" << maxElement;
1393 if (workOnAOD_ == 1) {
1394 maxElement = 0;
1395 if (debug_ > 1)
1396 edm::LogVerbatim("GammaJetAnalysis") << "forced 0";
1397 }
1398 HERE(Form("maxElement=%d", maxElement));
1399 for (int e = 0; e < maxElement; ++e) {
1400
1401 reco::PFBlockRef blockRef = (*it)->elementsInBlocks()[e].first;
1402 const edm::OwnVector<reco::PFBlockElement>& elements = blockRef->elements();
1403 for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
1404 if (elements[iEle].index() == (*it)->elementsInBlocks()[e].second) {
1405 if (elements[iEle].type() == reco::PFBlockElement::HCAL) {
1406 HF_type_ |= 0x1;
1407
1408 reco::PFClusterRef clusterref = elements[iEle].clusterRef();
1409 reco::PFCluster cluster = *clusterref;
1410 double cluster_dR = deltaR(ppfjet_eta_, ppfjet_phi_, cluster.eta(), cluster.phi());
1411 if (ppfjet_clusters.count(cluster_dR) == 0) {
1412 ppfjet_clusters[cluster_dR] = ppfjet_cluster_n_;
1413 ppfjet_cluster_eta_.push_back(cluster.eta());
1414 ppfjet_cluster_phi_.push_back(cluster.phi());
1415 ppfjet_cluster_dR_.push_back(cluster_dR);
1416 ppfjet_cluster_n_++;
1417 }
1418 int cluster_ind = ppfjet_clusters[cluster_dR];
1419 std::vector<std::pair<DetId, float>> hitsAndFracs = cluster.hitsAndFractions();
1420
1421
1422 int nHits = hitsAndFracs.size();
1423 for (int iHit = 0; iHit < nHits; iHit++) {
1424 int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
1425
1426 for (edm::SortedCollection<HBHERecHit, edm::StrictWeakOrdering<HBHERecHit>>::const_iterator ith =
1427 hbhereco->begin();
1428 ith != hbhereco->end();
1429 ++ith) {
1430 int etaPhiRecHit = getEtaPhi((*ith).id());
1431 if (etaPhiPF == etaPhiRecHit) {
1432 ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1433 if (ppfjet_rechits.count((*ith).id()) == 0) {
1434 ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1435 ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1436 ppfjet_twr_depth_.push_back((*ith).id().depth());
1437 ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1438 ppfjet_twr_hade_.push_back((*ith).energy());
1439 ppfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
1440 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1441 ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1442 ppfjet_twr_elmttype_.push_back(0);
1443 ppfjet_twr_clusterind_.push_back(cluster_ind);
1444 if (hasTrack) {
1445 ppfjet_twr_candtrackind_.push_back(ppfjet_ncandtracks_ - 1);
1446 } else {
1447 ppfjet_twr_candtrackind_.push_back(-1);
1448 }
1449 switch ((*ith).id().subdet()) {
1450 case HcalSubdetector::HcalBarrel: {
1451 CaloCellGeometry::CornersVec cv = HBGeom->getCorners((*ith).id());
1452 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1453 float avgphi =
1454 (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1455 if ((cv[0].phi() < cv[2].phi()) && (debug_ > 1))
1456 edm::LogVerbatim("GammaJetAnalysis") << "pHB" << cv[0].phi() << " " << cv[2].phi();
1457 if (cv[0].phi() < cv[2].phi())
1458 avgphi = (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) +
1459 static_cast<double>(cv[2].phi())) /
1460 2.0;
1461 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1462 break;
1463 }
1464 case HcalSubdetector::HcalEndcap: {
1465 CaloCellGeometry::CornersVec cv = HEGeom->getCorners((*ith).id());
1466 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1467 float avgphi =
1468 (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1469 if ((cv[0].phi() < cv[2].phi()) && (debug_ > 1))
1470 edm::LogVerbatim("GammaJetAnalysis") << "pHE" << cv[0].phi() << " " << cv[2].phi();
1471 if (cv[0].phi() < cv[2].phi())
1472 avgphi = (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) +
1473 static_cast<double>(cv[2].phi())) /
1474 2.0;
1475 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1476 break;
1477 }
1478 default:
1479 ppfjet_twr_dR_.push_back(-1);
1480 break;
1481 }
1482 ppfjet_rechits[(*ith).id()].first = ppfjet_ntwrs_;
1483 ++ppfjet_ntwrs_;
1484 } else if (ppfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0) {
1485 ppfjet_twr_frac_.at(ppfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
1486 if (cluster_dR <
1487 ppfjet_cluster_dR_.at(ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first))) {
1488 ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first) = cluster_ind;
1489 }
1490 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1491 }
1492 }
1493 }
1494 }
1495 }
1496 else if (elements[iEle].type() == reco::PFBlockElement::HFHAD) {
1497 types |= 0x2;
1498 ntypes++;
1499 HFHAD_n_++;
1500 HF_type_ |= 0x2;
1501
1502
1503
1504 for (edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith =
1505 hfreco->begin();
1506 ith != hfreco->end();
1507 ++ith) {
1508 if ((*ith).id().depth() == 1)
1509 continue;
1510 auto thisCell = HFGeom->getGeometry((*ith).id());
1511 const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
1512
1513 bool passMatch = false;
1514 if ((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()) {
1515 if ((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi())
1516 passMatch = true;
1517 else if (cv[0].phi() < cv[2].phi()) {
1518 if ((*it)->phi() < cv[0].phi())
1519 passMatch = true;
1520 else if ((*it)->phi() > cv[2].phi())
1521 passMatch = true;
1522 }
1523 }
1524
1525 if (passMatch) {
1526 ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1527 ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1528 ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1529 ppfjet_twr_depth_.push_back((*ith).id().depth());
1530 ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1531 ppfjet_twr_hade_.push_back((*ith).energy());
1532 ppfjet_twr_frac_.push_back(1.0);
1533 ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1534 ppfjet_twr_elmttype_.push_back(1);
1535 ppfjet_twr_clusterind_.push_back(-1);
1536 ppfjet_twr_candtrackind_.push_back(-1);
1537 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1538 float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1539 if ((cv[0].phi() < cv[2].phi()) && (debug_ > 1))
1540 edm::LogVerbatim("GammaJetAnalysis") << "pHFhad" << cv[0].phi() << " " << cv[2].phi();
1541 if (cv[0].phi() < cv[2].phi())
1542 avgphi =
1543 (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) /
1544 2.0;
1545 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1546 ++ppfjet_ntwrs_;
1547 HFHAD_E += (*ith).energy();
1548 }
1549 }
1550 } else if (elements[iEle].type() == reco::PFBlockElement::HFEM) {
1551 types |= 0x4;
1552 ntypes++;
1553 HFEM_n_++;
1554 HF_type_ |= 0x4;
1555
1556 for (edm::SortedCollection<HFRecHit, edm::StrictWeakOrdering<HFRecHit>>::const_iterator ith =
1557 hfreco->begin();
1558 ith != hfreco->end();
1559 ++ith) {
1560 if ((*ith).id().depth() == 2)
1561 continue;
1562 auto thisCell = HFGeom->getGeometry((*ith).id());
1563 const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
1564
1565 bool passMatch = false;
1566 if ((*it)->eta() < cv[0].eta() && (*it)->eta() > cv[2].eta()) {
1567 if ((*it)->phi() < cv[0].phi() && (*it)->phi() > cv[2].phi())
1568 passMatch = true;
1569 else if (cv[0].phi() < cv[2].phi()) {
1570 if ((*it)->phi() < cv[0].phi())
1571 passMatch = true;
1572 else if ((*it)->phi() > cv[2].phi())
1573 passMatch = true;
1574 }
1575 }
1576
1577 if (passMatch) {
1578 ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1579 ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1580 ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1581 ppfjet_twr_depth_.push_back((*ith).id().depth());
1582 ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1583 ppfjet_twr_hade_.push_back((*ith).energy());
1584 ppfjet_twr_frac_.push_back(1.0);
1585 ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1586 ppfjet_twr_elmttype_.push_back(2);
1587 ppfjet_twr_clusterind_.push_back(-1);
1588 ppfjet_twr_candtrackind_.push_back(-1);
1589 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1590 float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1591 if ((cv[0].phi() < cv[2].phi()) && (debug_ > 1))
1592 edm::LogVerbatim("GammaJetAnalysis") << "pHFem" << cv[0].phi() << " " << cv[2].phi();
1593 if (cv[0].phi() < cv[2].phi())
1594 avgphi =
1595 (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) /
1596 2.0;
1597 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1598 ++ppfjet_ntwrs_;
1599 HFEM_E += (*ith).energy();
1600 }
1601 }
1602 } else if (elements[iEle].type() == reco::PFBlockElement::HO) {
1603 types |= 0x8;
1604 ntypes++;
1605 HF_type_ |= 0x8;
1606 reco::PFClusterRef clusterref = elements[iEle].clusterRef();
1607 reco::PFCluster cluster = *clusterref;
1608 double cluster_dR = deltaR(ppfjet_eta_, ppfjet_phi_, cluster.eta(), cluster.phi());
1609 if (ppfjet_clusters.count(cluster_dR) == 0) {
1610 ppfjet_clusters[cluster_dR] = ppfjet_cluster_n_;
1611 ppfjet_cluster_eta_.push_back(cluster.eta());
1612 ppfjet_cluster_phi_.push_back(cluster.phi());
1613 ppfjet_cluster_dR_.push_back(cluster_dR);
1614 ppfjet_cluster_n_++;
1615 }
1616 int cluster_ind = ppfjet_clusters[cluster_dR];
1617
1618 std::vector<std::pair<DetId, float>> hitsAndFracs = cluster.hitsAndFractions();
1619 int nHits = hitsAndFracs.size();
1620 for (int iHit = 0; iHit < nHits; iHit++) {
1621 int etaPhiPF = getEtaPhi(hitsAndFracs[iHit].first);
1622
1623 for (edm::SortedCollection<HORecHit, edm::StrictWeakOrdering<HORecHit>>::const_iterator ith =
1624 horeco->begin();
1625 ith != horeco->end();
1626 ++ith) {
1627 int etaPhiRecHit = getEtaPhi((*ith).id());
1628 if (etaPhiPF == etaPhiRecHit) {
1629 ppfjet_had_ntwrs_.at(ppfjet_had_n_ - 1)++;
1630 if (ppfjet_rechits.count((*ith).id()) == 0) {
1631 ppfjet_twr_ieta_.push_back((*ith).id().ieta());
1632 ppfjet_twr_iphi_.push_back((*ith).id().iphi());
1633 ppfjet_twr_depth_.push_back((*ith).id().depth());
1634 ppfjet_twr_subdet_.push_back((*ith).id().subdet());
1635 ppfjet_twr_hade_.push_back((*ith).energy());
1636 ppfjet_twr_frac_.push_back(hitsAndFracs[iHit].second);
1637 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1638 ppfjet_twr_hadind_.push_back(ppfjet_had_n_ - 1);
1639 ppfjet_twr_elmttype_.push_back(3);
1640 ppfjet_twr_clusterind_.push_back(cluster_ind);
1641 if (hasTrack) {
1642 ppfjet_twr_candtrackind_.push_back(ppfjet_ncandtracks_ - 1);
1643 } else {
1644 ppfjet_twr_candtrackind_.push_back(-1);
1645 }
1646 auto thisCell = HOGeom->getGeometry((*ith).id());
1647 const CaloCellGeometry::CornersVec& cv = thisCell->getCorners();
1648 float avgeta = (cv[0].eta() + cv[2].eta()) / 2.0;
1649 float avgphi = (static_cast<double>(cv[0].phi()) + static_cast<double>(cv[2].phi())) / 2.0;
1650 if ((cv[0].phi() < cv[2].phi()) && (debug_ > 1))
1651 edm::LogVerbatim("GammaJetAnalysis") << "pHO" << cv[0].phi() << " " << cv[2].phi();
1652 if (cv[0].phi() < cv[2].phi())
1653 avgphi = (2.0 * 3.141592653 + static_cast<double>(cv[0].phi()) +
1654 static_cast<double>(cv[2].phi())) /
1655 2.0;
1656
1657 ppfjet_twr_dR_.push_back(deltaR(ppfjet_eta_, ppfjet_phi_, avgeta, avgphi));
1658 ppfjet_rechits[(*ith).id()].first = ppfjet_ntwrs_;
1659 ++ppfjet_ntwrs_;
1660 } else if (ppfjet_rechits[(*ith).id()].second.count(hitsAndFracs[iHit].second) == 0) {
1661 ppfjet_twr_frac_.at(ppfjet_rechits[(*ith).id()].first) += hitsAndFracs[iHit].second;
1662 if (cluster_dR <
1663 ppfjet_cluster_dR_.at(ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first))) {
1664 ppfjet_twr_clusterind_.at(ppfjet_rechits[(*ith).id()].first) = cluster_ind;
1665 }
1666 ppfjet_rechits[(*ith).id()].second.insert(hitsAndFracs[iHit].second);
1667 }
1668 }
1669 }
1670 }
1671 }
1672 }
1673 }
1674 }
1675 switch (candidateType) {
1676 case reco::PFCandidate::h_HF:
1677 ppfjet_had_emf_.push_back(HFEM_E / (HFEM_E + HFHAD_E));
1678 break;
1679 case reco::PFCandidate::egamma_HF:
1680 ppfjet_had_emf_.push_back(-1);
1681 break;
1682 default:
1683 ppfjet_had_emf_.push_back(-1);
1684 break;
1685 }
1686 }
1687
1688 if (doGenJets_) {
1689
1690 ppfjet_gendr_ = 99999.;
1691 ppfjet_genpt_ = 0;
1692 ppfjet_genp_ = 0;
1693 for (std::vector<reco::GenJet>::const_iterator it = genjets->begin(); it != genjets->end(); ++it) {
1694 const reco::GenJet* jet = &(*it);
1695 double dr = deltaR(jet, pfjet_probe.jet());
1696 if (dr < ppfjet_gendr_) {
1697 ppfjet_gendr_ = dr;
1698 ppfjet_genpt_ = jet->pt();
1699 ppfjet_genp_ = jet->p();
1700 ppfjet_genE_ = jet->energy();
1701 }
1702 }
1703 }
1704 if (iJet == 2) {
1705 copy_leadingPfJetVars_to_pfJet2();
1706 }
1707 }
1708
1709
1710
1711
1712
1713 const edm::Handle<reco::PFMETCollection> pfmet_h = iEvent.getHandle(tok_PFMET_);
1714 if (!pfmet_h.isValid()) {
1715 edm::LogWarning("GammaJetAnalysis") << " could not find " << pfMETColl;
1716 return;
1717 }
1718 met_value_ = pfmet_h->begin()->et();
1719 met_phi_ = pfmet_h->begin()->phi();
1720 met_sumEt_ = pfmet_h->begin()->sumEt();
1721
1722 const edm::Handle<reco::PFMETCollection> pfmetType1_h = iEvent.getHandle(tok_PFType1MET_);
1723 if (pfmetType1_h.isValid()) {
1724 metType1_value_ = pfmetType1_h->begin()->et();
1725 metType1_phi_ = pfmetType1_h->begin()->phi();
1726 metType1_sumEt_ = pfmetType1_h->begin()->sumEt();
1727 } else {
1728
1729 metType1_value_ = -999.;
1730 metType1_phi_ = -999.;
1731 metType1_sumEt_ = -999.;
1732 }
1733
1734
1735 pf_tree_->Fill();
1736 }
1737 }
1738 return;
1739 }
1740
1741
1742 void GammaJetAnalysis::beginJob() {
1743 edm::Service<TFileService> fs;
1744 if (doPFJets_) {
1745 pf_tree_ = fs->make<TTree>("pf_gammajettree", "tree for gamma+jet balancing using PFJets");
1746 }
1747
1748 for (int iJet = 0; iJet < 2; iJet++) {
1749 bool doJet = doPFJets_;
1750 if (!doJet)
1751 continue;
1752 if (iJet > 0)
1753 continue;
1754 TTree* tree = pf_tree_;
1755
1756
1757 tree->Branch("photonTrig_fired", &photonTrigFired_);
1758 tree->Branch("photonTrig_prescale", &photonTrigPrescale_);
1759 tree->Branch("jetTrig_fired", &jetTrigFired_);
1760 tree->Branch("jetTrig_prescale", &jetTrigPrescale_);
1761
1762
1763 tree->Branch("RunNumber", &runNumber_, "RunNumber/I");
1764 tree->Branch("LumiBlock", &lumiBlock_, "LumiBlock/I");
1765 tree->Branch("EventNumber", &eventNumber_, "EventNumber/I");
1766 tree->Branch("EventWeight", &eventWeight_, "EventWeight/F");
1767 tree->Branch("EventPtHat", &eventPtHat_, "EventPtHat/F");
1768
1769
1770 tree->Branch("rho2012", &rho2012_, "rho2012/F");
1771 tree->Branch("tagPho_pt", &tagPho_pt_, "tagPho_pt/F");
1772 tree->Branch("pho_2nd_pt", &pho_2nd_pt_, "pho_2nd_pt/F");
1773 tree->Branch("tagPho_energy", &tagPho_energy_, "tagPho_energy/F");
1774 tree->Branch("tagPho_eta", &tagPho_eta_, "tagPho_eta/F");
1775 tree->Branch("tagPho_phi", &tagPho_phi_, "tagPho_phi/F");
1776 tree->Branch("tagPho_sieie", &tagPho_sieie_, "tagPho_sieie/F");
1777 tree->Branch("tagPho_HoE", &tagPho_HoE_, "tagPho_HoE/F");
1778 tree->Branch("tagPho_r9", &tagPho_r9_, "tagPho_r9/F");
1779 tree->Branch("tagPho_EcalIsoDR04", &tagPho_EcalIsoDR04_, "tagPho_EcalIsoDR04/F");
1780 tree->Branch("tagPho_HcalIsoDR04", &tagPho_HcalIsoDR04_, "tagPho_HcalIsoDR04/F");
1781 tree->Branch("tagPho_HcalIsoDR0412", &tagPho_HcalIsoDR0412_, "tagPho_HcalIsoDR0412/F");
1782 tree->Branch("tagPho_TrkIsoHollowDR04", &tagPho_TrkIsoHollowDR04_, "tagPho_TrkIsoHollowDR04/F");
1783 tree->Branch("tagPho_pfiso_myphoton03", &tagPho_pfiso_myphoton03_, "tagPho_pfiso_myphoton03/F");
1784 tree->Branch("tagPho_pfiso_myneutral03", &tagPho_pfiso_myneutral03_, "tagPho_pfiso_myneutral03/F");
1785 tree->Branch("tagPho_pfiso_mycharged03", "std::vector<std::vector<float> >", &tagPho_pfiso_mycharged03);
1786 tree->Branch("tagPho_pixelSeed", &tagPho_pixelSeed_, "tagPho_pixelSeed/I");
1787 tree->Branch("tagPho_ConvSafeEleVeto", &tagPho_ConvSafeEleVeto_, "tagPho_ConvSafeEleVeto/I");
1788 tree->Branch("tagPho_idTight", &tagPho_idTight_, "tagPho_idTight/I");
1789 tree->Branch("tagPho_idLoose", &tagPho_idLoose_, "tagPho_idLoose/I");
1790
1791 if (doGenJets_) {
1792 tree->Branch("tagPho_genPt", &tagPho_genPt_, "tagPho_genPt/F");
1793 tree->Branch("tagPho_genEnergy", &tagPho_genEnergy_, "tagPho_genEnergy/F");
1794 tree->Branch("tagPho_genEta", &tagPho_genEta_, "tagPho_genEta/F");
1795 tree->Branch("tagPho_genPhi", &tagPho_genPhi_, "tagPho_genPhi/F");
1796 tree->Branch("tagPho_genDeltaR", &tagPho_genDeltaR_, "tagPho_genDeltaR/F");
1797 }
1798
1799 tree->Branch("nPhotons", &nPhotons_, "nPhotons/I");
1800 tree->Branch("nGenJets", &nGenJets_, "nGenJets/I");
1801 }
1802
1803
1804
1805 if (doPFJets_) {
1806 pf_tree_->Branch("nPFJets", &nPFJets_, "nPFJets/I");
1807
1808
1809 pf_tree_->Branch("ppfjet_pt", &ppfjet_pt_, "ppfjet_pt/F");
1810 pf_tree_->Branch("ppfjet_p", &ppfjet_p_, "ppfjet_p/F");
1811 pf_tree_->Branch("ppfjet_E", &ppfjet_E_, "ppfjet_E/F");
1812 pf_tree_->Branch("ppfjet_E_NPVcorr", &ppfjet_E_NPVcorr_, "ppfjet_E_NPVcorr/F");
1813 pf_tree_->Branch("ppfjet_area", &ppfjet_area_, "ppfjet_area/F");
1814 pf_tree_->Branch("ppfjet_eta", &ppfjet_eta_, "ppfjet_eta/F");
1815 pf_tree_->Branch("ppfjet_phi", &ppfjet_phi_, "ppfjet_phi/F");
1816 pf_tree_->Branch("ppfjet_scale", &ppfjet_scale_, "ppfjet_scale/F");
1817 pf_tree_->Branch("ppfjet_NeutralHadronFrac", &ppfjet_NeutralHadronFrac_, "ppfjet_NeutralHadronFrac/F");
1818 pf_tree_->Branch("ppfjet_NeutralEMFrac", &ppfjet_NeutralEMFrac_, "ppfjet_NeutralEMFrac/F");
1819 pf_tree_->Branch("ppfjet_nConstituents", &ppfjet_nConstituents_, "ppfjet_nConstituents/I");
1820 pf_tree_->Branch("ppfjet_ChargedHadronFrac", &ppfjet_ChargedHadronFrac_, "ppfjet_ChargedHadronFrac/F");
1821 pf_tree_->Branch("ppfjet_ChargedMultiplicity", &ppfjet_ChargedMultiplicity_, "ppfjet_ChargedMultiplicity/F");
1822 pf_tree_->Branch("ppfjet_ChargedEMFrac", &ppfjet_ChargedEMFrac_, "ppfjet_ChargedEMFrac/F");
1823 if (doGenJets_) {
1824 pf_tree_->Branch("ppfjet_genpt", &ppfjet_genpt_, "ppfjet_genpt/F");
1825 pf_tree_->Branch("ppfjet_genp", &ppfjet_genp_, "ppfjet_genp/F");
1826 pf_tree_->Branch("ppfjet_genE", &ppfjet_genE_, "ppfjet_genE/F");
1827 pf_tree_->Branch("ppfjet_gendr", &ppfjet_gendr_, "ppfjet_gendr/F");
1828 }
1829 pf_tree_->Branch("ppfjet_unkown_E", &ppfjet_unkown_E_, "ppfjet_unkown_E/F");
1830 pf_tree_->Branch("ppfjet_electron_E", &ppfjet_electron_E_, "ppfjet_electron_E/F");
1831 pf_tree_->Branch("ppfjet_muon_E", &ppfjet_muon_E_, "ppfjet_muon_E/F");
1832 pf_tree_->Branch("ppfjet_photon_E", &ppfjet_photon_E_, "ppfjet_photon_E/F");
1833 pf_tree_->Branch("ppfjet_unkown_px", &ppfjet_unkown_px_, "ppfjet_unkown_px/F");
1834 pf_tree_->Branch("ppfjet_electron_px", &ppfjet_electron_px_, "ppfjet_electron_px/F");
1835 pf_tree_->Branch("ppfjet_muon_px", &ppfjet_muon_px_, "ppfjet_muon_px/F");
1836 pf_tree_->Branch("ppfjet_photon_px", &ppfjet_photon_px_, "ppfjet_photon_px/F");
1837 pf_tree_->Branch("ppfjet_unkown_py", &ppfjet_unkown_py_, "ppfjet_unkown_py/F");
1838 pf_tree_->Branch("ppfjet_electron_py", &ppfjet_electron_py_, "ppfjet_electron_py/F");
1839 pf_tree_->Branch("ppfjet_muon_py", &ppfjet_muon_py_, "ppfjet_muon_py/F");
1840 pf_tree_->Branch("ppfjet_photon_py", &ppfjet_photon_py_, "ppfjet_photon_py/F");
1841 pf_tree_->Branch("ppfjet_unkown_pz", &ppfjet_unkown_pz_, "ppfjet_unkown_pz/F");
1842 pf_tree_->Branch("ppfjet_electron_pz", &ppfjet_electron_pz_, "ppfjet_electron_pz/F");
1843 pf_tree_->Branch("ppfjet_muon_pz", &ppfjet_muon_pz_, "ppfjet_muon_pz/F");
1844 pf_tree_->Branch("ppfjet_photon_pz", &ppfjet_photon_pz_, "ppfjet_photon_pz/F");
1845 pf_tree_->Branch("ppfjet_unkown_EcalE", &ppfjet_unkown_EcalE_, "ppfjet_unkown_EcalE/F");
1846 pf_tree_->Branch("ppfjet_electron_EcalE", &ppfjet_electron_EcalE_, "ppfjet_electron_EcalE/F");
1847 pf_tree_->Branch("ppfjet_muon_EcalE", &ppfjet_muon_EcalE_, "ppfjet_muon_EcalE/F");
1848 pf_tree_->Branch("ppfjet_photon_EcalE", &ppfjet_photon_EcalE_, "ppfjet_photon_EcalE/F");
1849 pf_tree_->Branch("ppfjet_unkown_n", &ppfjet_unkown_n_, "ppfjet_unkown_n/I");
1850 pf_tree_->Branch("ppfjet_electron_n", &ppfjet_electron_n_, "ppfjet_electron_n/I");
1851 pf_tree_->Branch("ppfjet_muon_n", &ppfjet_muon_n_, "ppfjet_muon_n/I");
1852 pf_tree_->Branch("ppfjet_photon_n", &ppfjet_photon_n_, "ppfjet_photon_n/I");
1853 pf_tree_->Branch("ppfjet_had_n", &ppfjet_had_n_, "ppfjet_had_n/I");
1854 pf_tree_->Branch("ppfjet_had_E", &ppfjet_had_E_);
1855 pf_tree_->Branch("ppfjet_had_px", &ppfjet_had_px_);
1856 pf_tree_->Branch("ppfjet_had_py", &ppfjet_had_py_);
1857 pf_tree_->Branch("ppfjet_had_pz", &ppfjet_had_pz_);
1858 pf_tree_->Branch("ppfjet_had_EcalE", &ppfjet_had_EcalE_);
1859 pf_tree_->Branch("ppfjet_had_rawHcalE", &ppfjet_had_rawHcalE_);
1860 pf_tree_->Branch("ppfjet_had_emf", &ppfjet_had_emf_);
1861 pf_tree_->Branch("ppfjet_had_id", &ppfjet_had_id_);
1862 pf_tree_->Branch("ppfjet_had_candtrackind", &ppfjet_had_candtrackind_);
1863 if (doGenJets_) {
1864 pf_tree_->Branch("ppfjet_had_E_mctruth", &ppfjet_had_E_mctruth_);
1865 pf_tree_->Branch("ppfjet_had_mcpdgId", &ppfjet_had_mcpdgId_);
1866 }
1867 pf_tree_->Branch("ppfjet_had_ntwrs", &ppfjet_had_ntwrs_);
1868 pf_tree_->Branch("ppfjet_ntwrs", &ppfjet_ntwrs_, "ppfjet_ntwrs/I");
1869 pf_tree_->Branch("ppfjet_twr_ieta", &ppfjet_twr_ieta_);
1870 pf_tree_->Branch("ppfjet_twr_iphi", &ppfjet_twr_iphi_);
1871 pf_tree_->Branch("ppfjet_twr_depth", &ppfjet_twr_depth_);
1872 pf_tree_->Branch("ppfjet_twr_subdet", &ppfjet_twr_subdet_);
1873 pf_tree_->Branch("ppfjet_twr_hade", &ppfjet_twr_hade_);
1874 pf_tree_->Branch("ppfjet_twr_frac", &ppfjet_twr_frac_);
1875 pf_tree_->Branch("ppfjet_twr_candtrackind", &ppfjet_twr_candtrackind_);
1876 pf_tree_->Branch("ppfjet_twr_hadind", &ppfjet_twr_hadind_);
1877 pf_tree_->Branch("ppfjet_twr_elmttype", &ppfjet_twr_elmttype_);
1878 pf_tree_->Branch("ppfjet_twr_dR", &ppfjet_twr_dR_);
1879 pf_tree_->Branch("ppfjet_twr_clusterind", &ppfjet_twr_clusterind_);
1880 pf_tree_->Branch("ppfjet_cluster_n", &ppfjet_cluster_n_, "ppfjet_cluster_n/I");
1881 pf_tree_->Branch("ppfjet_cluster_eta", &ppfjet_cluster_eta_);
1882 pf_tree_->Branch("ppfjet_cluster_phi", &ppfjet_cluster_phi_);
1883 pf_tree_->Branch("ppfjet_cluster_dR", &ppfjet_cluster_dR_);
1884 pf_tree_->Branch("ppfjet_ncandtracks", &ppfjet_ncandtracks_, "ppfjet_ncandtracks/I");
1885 pf_tree_->Branch("ppfjet_candtrack_px", &ppfjet_candtrack_px_);
1886 pf_tree_->Branch("ppfjet_candtrack_py", &ppfjet_candtrack_py_);
1887 pf_tree_->Branch("ppfjet_candtrack_pz", &ppfjet_candtrack_pz_);
1888 pf_tree_->Branch("ppfjet_candtrack_EcalE", &ppfjet_candtrack_EcalE_);
1889
1890
1891 pf_tree_->Branch("pfjet2_pt", &pfjet2_pt_, "pfjet2_pt/F");
1892 pf_tree_->Branch("pfjet2_p", &pfjet2_p_, "pfjet2_p/F");
1893 pf_tree_->Branch("pfjet2_E", &pfjet2_E_, "pfjet2_E/F");
1894 pf_tree_->Branch("pfjet2_E_NPVcorr", &pfjet2_E_NPVcorr_, "pfjet2_E_NPVcorr/F");
1895 pf_tree_->Branch("pfjet2_area", &pfjet2_area_, "pfjet2_area/F");
1896 pf_tree_->Branch("pfjet2_eta", &pfjet2_eta_, "pfjet2_eta/F");
1897 pf_tree_->Branch("pfjet2_phi", &pfjet2_phi_, "pfjet2_phi/F");
1898 pf_tree_->Branch("pfjet2_scale", &pfjet2_scale_, "pfjet2_scale/F");
1899 pf_tree_->Branch("pfjet2_NeutralHadronFrac", &pfjet2_NeutralHadronFrac_, "pfjet2_NeutralHadronFrac/F");
1900 pf_tree_->Branch("pfjet2_NeutralEMFrac", &pfjet2_NeutralEMFrac_, "pfjet2_NeutralEMFrac/F");
1901 pf_tree_->Branch("pfjet2_nConstituents", &pfjet2_nConstituents_, "pfjet2_nConstituents/I");
1902 pf_tree_->Branch("pfjet2_ChargedHadronFrac", &pfjet2_ChargedHadronFrac_, "pfjet2_ChargedHadronFrac/F");
1903 pf_tree_->Branch("pfjet2_ChargedMultiplicity", &pfjet2_ChargedMultiplicity_, "pfjet2_ChargedMultiplicity/F");
1904 pf_tree_->Branch("pfjet2_ChargedEMFrac", &pfjet2_ChargedEMFrac_, "pfjet2_ChargedEMFrac/F");
1905 if (doGenJets_) {
1906 pf_tree_->Branch("pfjet2_genpt", &pfjet2_genpt_, "pfjet2_genpt/F");
1907 pf_tree_->Branch("pfjet2_genp", &pfjet2_genp_, "pfjet2_genp/F");
1908 pf_tree_->Branch("pfjet2_genE", &pfjet2_genE_, "pfjet2_genE/F");
1909 pf_tree_->Branch("pfjet2_gendr", &pfjet2_gendr_, "pfjet2_gendr/F");
1910 }
1911 pf_tree_->Branch("pfjet2_unkown_E", &pfjet2_unkown_E_, "pfjet2_unkown_E/F");
1912 pf_tree_->Branch("pfjet2_electron_E", &pfjet2_electron_E_, "pfjet2_electron_E/F");
1913 pf_tree_->Branch("pfjet2_muon_E", &pfjet2_muon_E_, "pfjet2_muon_E/F");
1914 pf_tree_->Branch("pfjet2_photon_E", &pfjet2_photon_E_, "pfjet2_photon_E/F");
1915 pf_tree_->Branch("pfjet2_unkown_px", &pfjet2_unkown_px_, "pfjet2_unkown_px/F");
1916 pf_tree_->Branch("pfjet2_electron_px", &pfjet2_electron_px_, "pfjet2_electron_px/F");
1917 pf_tree_->Branch("pfjet2_muon_px", &pfjet2_muon_px_, "pfjet2_muon_px/F");
1918 pf_tree_->Branch("pfjet2_photon_px", &pfjet2_photon_px_, "pfjet2_photon_px/F");
1919 pf_tree_->Branch("pfjet2_unkown_py", &pfjet2_unkown_py_, "pfjet2_unkown_py/F");
1920 pf_tree_->Branch("pfjet2_electron_py", &pfjet2_electron_py_, "pfjet2_electron_py/F");
1921 pf_tree_->Branch("pfjet2_muon_py", &pfjet2_muon_py_, "pfjet2_muon_py/F");
1922 pf_tree_->Branch("pfjet2_photon_py", &pfjet2_photon_py_, "pfjet2_photon_py/F");
1923 pf_tree_->Branch("pfjet2_unkown_pz", &pfjet2_unkown_pz_, "pfjet2_unkown_pz/F");
1924 pf_tree_->Branch("pfjet2_electron_pz", &pfjet2_electron_pz_, "pfjet2_electron_pz/F");
1925 pf_tree_->Branch("pfjet2_muon_pz", &pfjet2_muon_pz_, "pfjet2_muon_pz/F");
1926 pf_tree_->Branch("pfjet2_photon_pz", &pfjet2_photon_pz_, "pfjet2_photon_pz/F");
1927 pf_tree_->Branch("pfjet2_unkown_EcalE", &pfjet2_unkown_EcalE_, "pfjet2_unkown_EcalE/F");
1928 pf_tree_->Branch("pfjet2_electron_EcalE", &pfjet2_electron_EcalE_, "pfjet2_electron_EcalE/F");
1929 pf_tree_->Branch("pfjet2_muon_EcalE", &pfjet2_muon_EcalE_, "pfjet2_muon_EcalE/F");
1930 pf_tree_->Branch("pfjet2_photon_EcalE", &pfjet2_photon_EcalE_, "pfjet2_photon_EcalE/F");
1931 pf_tree_->Branch("pfjet2_unkown_n", &pfjet2_unkown_n_, "pfjet2_unkown_n/I");
1932 pf_tree_->Branch("pfjet2_electron_n", &pfjet2_electron_n_, "pfjet2_electron_n/I");
1933 pf_tree_->Branch("pfjet2_muon_n", &pfjet2_muon_n_, "pfjet2_muon_n/I");
1934 pf_tree_->Branch("pfjet2_photon_n", &pfjet2_photon_n_, "pfjet2_photon_n/I");
1935 pf_tree_->Branch("pfjet2_had_n", &pfjet2_had_n_, "pfjet2_had_n/I");
1936 pf_tree_->Branch("pfjet2_had_E", &pfjet2_had_E_);
1937 pf_tree_->Branch("pfjet2_had_px", &pfjet2_had_px_);
1938 pf_tree_->Branch("pfjet2_had_py", &pfjet2_had_py_);
1939 pf_tree_->Branch("pfjet2_had_pz", &pfjet2_had_pz_);
1940 pf_tree_->Branch("pfjet2_had_EcalE", &pfjet2_had_EcalE_);
1941 pf_tree_->Branch("pfjet2_had_rawHcalE", &pfjet2_had_rawHcalE_);
1942 pf_tree_->Branch("pfjet2_had_emf", &pfjet2_had_emf_);
1943 pf_tree_->Branch("pfjet2_had_id", &pfjet2_had_id_);
1944 pf_tree_->Branch("pfjet2_had_candtrackind", &pfjet2_had_candtrackind_);
1945 if (doGenJets_) {
1946 pf_tree_->Branch("pfjet2_had_E_mctruth", &pfjet2_had_E_mctruth_);
1947 pf_tree_->Branch("pfjet2_had_mcpdgId", &pfjet2_had_mcpdgId_);
1948 }
1949 pf_tree_->Branch("pfjet2_had_ntwrs", &pfjet2_had_ntwrs_);
1950 pf_tree_->Branch("pfjet2_ntwrs", &pfjet2_ntwrs_, "pfjet2_ntwrs/I");
1951 pf_tree_->Branch("pfjet2_twr_ieta", &pfjet2_twr_ieta_);
1952 pf_tree_->Branch("pfjet2_twr_iphi", &pfjet2_twr_iphi_);
1953 pf_tree_->Branch("pfjet2_twr_depth", &pfjet2_twr_depth_);
1954 pf_tree_->Branch("pfjet2_twr_subdet", &pfjet2_twr_subdet_);
1955 pf_tree_->Branch("pfjet2_twr_hade", &pfjet2_twr_hade_);
1956 pf_tree_->Branch("pfjet2_twr_frac", &pfjet2_twr_frac_);
1957 pf_tree_->Branch("pfjet2_twr_candtrackind", &pfjet2_twr_candtrackind_);
1958 pf_tree_->Branch("pfjet2_twr_hadind", &pfjet2_twr_hadind_);
1959 pf_tree_->Branch("pfjet2_twr_elmttype", &pfjet2_twr_elmttype_);
1960 pf_tree_->Branch("pfjet2_twr_dR", &pfjet2_twr_dR_);
1961 pf_tree_->Branch("pfjet2_twr_clusterind", &pfjet2_twr_clusterind_);
1962 pf_tree_->Branch("pfjet2_cluster_n", &pfjet2_cluster_n_, "pfjet2_cluster_n/I");
1963 pf_tree_->Branch("pfjet2_cluster_eta", &pfjet2_cluster_eta_);
1964 pf_tree_->Branch("pfjet2_cluster_phi", &pfjet2_cluster_phi_);
1965 pf_tree_->Branch("pfjet2_cluster_dR", &pfjet2_cluster_dR_);
1966 pf_tree_->Branch("pfjet2_ncandtracks", &pfjet2_ncandtracks_, "pfjet2_ncandtracks/I");
1967 pf_tree_->Branch("pfjet2_candtrack_px", &pfjet2_candtrack_px_);
1968 pf_tree_->Branch("pfjet2_candtrack_py", &pfjet2_candtrack_py_);
1969 pf_tree_->Branch("pfjet2_candtrack_pz", &pfjet2_candtrack_pz_);
1970 pf_tree_->Branch("pfjet2_candtrack_EcalE", &pfjet2_candtrack_EcalE_);
1971
1972
1973 pf_tree_->Branch("pf_thirdjet_et", &pf_thirdjet_et_, "pf_thirdjet_et/F");
1974 pf_tree_->Branch("pf_thirdjet_pt", &pf_thirdjet_pt_, "pf_thirdjet_pt/F");
1975 pf_tree_->Branch("pf_thirdjet_p", &pf_thirdjet_p_, "pf_thirdjet_p/F");
1976 pf_tree_->Branch("pf_thirdjet_px", &pf_thirdjet_px_, "pf_thirdjet_px/F");
1977 pf_tree_->Branch("pf_thirdjet_py", &pf_thirdjet_py_, "pf_thirdjet_py/F");
1978 pf_tree_->Branch("pf_thirdjet_E", &pf_thirdjet_E_, "pf_thirdjet_E/F");
1979 pf_tree_->Branch("pf_thirdjet_eta", &pf_thirdjet_eta_, "pf_thirdjet_eta/F");
1980 pf_tree_->Branch("pf_thirdjet_phi", &pf_thirdjet_phi_, "pf_thirdjet_phi/F");
1981 pf_tree_->Branch("pf_thirdjet_scale", &pf_thirdjet_scale_, "pf_thirdjet_scale/F");
1982
1983 pf_tree_->Branch("met_value", &met_value_, "met_value/F");
1984 pf_tree_->Branch("met_phi", &met_phi_, "met_phi/F");
1985 pf_tree_->Branch("met_sumEt", &met_sumEt_, "met_sumEt/F");
1986 pf_tree_->Branch("metType1_value", &metType1_value_, "metType1_value/F");
1987 pf_tree_->Branch("metType1_phi", &metType1_phi_, "metType1_phi/F");
1988 pf_tree_->Branch("metType1_sumEt", &metType1_sumEt_, "metType1_sumEt/F");
1989 pf_tree_->Branch("pf_NPV", &pf_NPV_, "pf_NPV/I");
1990 }
1991
1992 return;
1993 }
1994
1995
1996 void GammaJetAnalysis::endJob() {
1997 if (doPFJets_) {
1998 pf_tree_->Write();
1999 }
2000
2001
2002 {
2003 edm::Service<TFileService> fs;
2004 misc_tree_ = fs->make<TTree>("misc_tree", "tree for misc.info");
2005 misc_tree_->Branch("ignoreHLT", &ignoreHLT_, "ignoreHLT/O");
2006 misc_tree_->Branch("doPFJets", &doPFJets_, "doPFJets/O");
2007 misc_tree_->Branch("doGenJets", &doGenJets_, "doGenJets/O");
2008 misc_tree_->Branch("workOnAOD", &workOnAOD_, "workOnAOD/O");
2009 misc_tree_->Branch("photonTriggerNames", &photonTrigNamesV_);
2010 misc_tree_->Branch("jetTriggerNames", &jetTrigNamesV_);
2011 misc_tree_->Branch("nProcessed", &nProcessed_, "nProcessed/l");
2012
2013 time_t ltime;
2014 ltime = time(NULL);
2015 TString str = TString(asctime(localtime(<ime)));
2016 if (str[str.Length() - 1] == '\n')
2017 str.Remove(str.Length() - 1, 1);
2018 TObjString date(str);
2019 date.Write(str.Data());
2020 misc_tree_->Fill();
2021 misc_tree_->Write();
2022 }
2023 }
2024
2025
2026
2027 void GammaJetAnalysis::beginRun(const edm::Run& iRun, const edm::EventSetup& setup) {
2028 if (debug_ > 1)
2029 edm::LogVerbatim("GammaJetAnalysis") << "beginRun()";
2030
2031 if (!ignoreHLT_) {
2032 int noPhotonTrigger = (photonTrigNamesV_.size() == 0) ? 1 : 0;
2033 int noJetTrigger = (jetTrigNamesV_.size() == 0) ? 1 : 0;
2034 if (!noPhotonTrigger && (photonTrigNamesV_.size() == 1) && (photonTrigNamesV_[0].length() == 0))
2035 noPhotonTrigger = 1;
2036 if (!noJetTrigger && (jetTrigNamesV_.size() == 1) && (jetTrigNamesV_[0].length() == 0))
2037 noJetTrigger = 1;
2038 if (noPhotonTrigger && noJetTrigger) {
2039 ignoreHLT_ = true;
2040 if (debug_ > 1)
2041 edm::LogVerbatim("GammaJetAnalysis") << "HLT trigger ignored: no trigger requested";
2042 }
2043 } else {
2044
2045 photonTrigNamesV_.clear();
2046 jetTrigNamesV_.clear();
2047 }
2048
2049 if (!ignoreHLT_) {
2050 if (debug_ > 0)
2051 edm::LogVerbatim("GammaJetAnalysis") << "Initializing trigger information for individual run";
2052 bool changed(true);
2053 std::string processName = "HLT";
2054 if (hltPrescaleProvider_.init(iRun, setup, processName, changed)) {
2055
2056 if (changed) {
2057
2058
2059 }
2060 } else {
2061
2062
2063 throw edm::Exception(edm::errors::ProductNotFound)
2064 << " HLT config extraction failure with process name " << processName;
2065
2066 }
2067 }
2068 }
2069
2070
2071
2072
2073
2074 float GammaJetAnalysis::pfEcalIso(const reco::Photon* localPho1,
2075 edm::Handle<reco::PFCandidateCollection> pfHandle,
2076 float dRmax,
2077 float dRVetoBarrel,
2078 float dRVetoEndcap,
2079 float etaStripBarrel,
2080 float etaStripEndcap,
2081 float energyBarrel,
2082 float energyEndcap,
2083 reco::PFCandidate::ParticleType pfToUse) {
2084 if (debug_ > 1)
2085 edm::LogVerbatim("GammaJetAnalysis") << "Inside pfEcalIso";
2086 reco::Photon* localPho = localPho1->clone();
2087 float dRVeto;
2088 float etaStrip;
2089
2090 if (localPho->isEB()) {
2091 dRVeto = dRVetoBarrel;
2092 etaStrip = etaStripBarrel;
2093 } else {
2094 dRVeto = dRVetoEndcap;
2095 etaStrip = etaStripEndcap;
2096 }
2097 const reco::PFCandidateCollection* forIsolation = pfHandle.product();
2098 int nsize = forIsolation->size();
2099 float sum = 0;
2100 for (int i = 0; i < nsize; i++) {
2101 const reco::PFCandidate& pfc = (*forIsolation)[i];
2102 if (pfc.particleId() == pfToUse) {
2103
2104 if (pfc.superClusterRef().isNonnull() && localPho->superCluster().isNonnull()) {
2105 if (pfc.superClusterRef() == localPho->superCluster())
2106 continue;
2107 }
2108
2109 if (localPho->isEB()) {
2110 if (fabs(pfc.pt()) < energyBarrel)
2111 continue;
2112 } else {
2113 if (fabs(pfc.energy()) < energyEndcap)
2114 continue;
2115 }
2116
2117 math::XYZPoint pfvtx = pfc.vertex();
2118 math::XYZVector photon_directionWrtVtx(localPho->superCluster()->x() - pfvtx.x(),
2119 localPho->superCluster()->y() - pfvtx.y(),
2120 localPho->superCluster()->z() - pfvtx.z());
2121
2122 float dEta = fabs(photon_directionWrtVtx.Eta() - pfc.momentum().Eta());
2123 float dR = deltaR(
2124 photon_directionWrtVtx.Eta(), photon_directionWrtVtx.Phi(), pfc.momentum().Eta(), pfc.momentum().Phi());
2125
2126 if (dEta < etaStrip)
2127 continue;
2128
2129 if (dR > dRmax || dR < dRVeto)
2130 continue;
2131
2132 sum += pfc.pt();
2133 }
2134 }
2135 return sum;
2136 }
2137
2138
2139
2140 float GammaJetAnalysis::pfHcalIso(const reco::Photon* localPho,
2141 edm::Handle<reco::PFCandidateCollection> pfHandle,
2142 float dRmax,
2143 float dRveto,
2144 reco::PFCandidate::ParticleType pfToUse) {
2145 if (debug_ > 1)
2146 edm::LogVerbatim("GammaJetAnalysis") << "Inside pfHcalIso";
2147 return pfEcalIso(localPho, pfHandle, dRmax, dRveto, dRveto, 0.0, 0.0, 0.0, 0.0, pfToUse);
2148 }
2149
2150
2151
2152 std::vector<float> GammaJetAnalysis::pfTkIsoWithVertex(const reco::Photon* localPho1,
2153 edm::Handle<reco::PFCandidateCollection> pfHandle,
2154 edm::Handle<reco::VertexCollection> vtxHandle,
2155 float dRmax,
2156 float dRvetoBarrel,
2157 float dRvetoEndcap,
2158 float ptMin,
2159 float dzMax,
2160 float dxyMax,
2161 reco::PFCandidate::ParticleType pfToUse) {
2162 if (debug_ > 1)
2163 edm::LogVerbatim("GammaJetAnalysis") << "Inside pfTkIsoWithVertex()";
2164 reco::Photon* localPho = localPho1->clone();
2165
2166 float dRveto;
2167 if (localPho->isEB())
2168 dRveto = dRvetoBarrel;
2169 else
2170 dRveto = dRvetoEndcap;
2171
2172 std::vector<float> result;
2173 const reco::PFCandidateCollection* forIsolation = pfHandle.product();
2174
2175
2176 if (debug_ > 1)
2177 edm::LogVerbatim("GammaJetAnalysis") << "vtxHandle->size() = " << vtxHandle->size();
2178 for (unsigned int ivtx = 0; ivtx < (vtxHandle->size()); ++ivtx) {
2179 if (debug_ > 1)
2180 edm::LogVerbatim("GammaJetAnalysis") << "Vtx " << ivtx;
2181
2182 reco::VertexRef vtx(vtxHandle, ivtx);
2183 math::XYZVector photon_directionWrtVtx(localPho->superCluster()->x() - vtx->x(),
2184 localPho->superCluster()->y() - vtx->y(),
2185 localPho->superCluster()->z() - vtx->z());
2186 if (debug_ > 1)
2187 edm::LogVerbatim("GammaJetAnalysis") << "pfTkIsoWithVertex :: Will Loop over the PFCandidates";
2188 float sum = 0;
2189
2190 for (unsigned i = 0; i < forIsolation->size(); i++) {
2191 if (debug_ > 1)
2192 edm::LogVerbatim("GammaJetAnalysis") << "inside loop";
2193 const reco::PFCandidate& pfc = (*forIsolation)[i];
2194
2195
2196 if (debug_ > 1) {
2197 edm::LogVerbatim("GammaJetAnalysis") << "pfToUse=" << pfToUse;
2198 edm::LogVerbatim("GammaJetAnalysis") << "pfc.particleId()=" << pfc.particleId();
2199 }
2200
2201 if (pfc.particleId() == pfToUse) {
2202 if (debug_ > 1) {
2203 edm::LogVerbatim("GammaJetAnalysis") << "\n ***** HERE pfc.particleId() == pfToUse ";
2204 edm::LogVerbatim("GammaJetAnalysis") << "pfc.pt()=" << pfc.pt();
2205 }
2206 if (pfc.pt() < ptMin)
2207 continue;
2208
2209 float dz = fabs(pfc.trackRef()->dz(vtx->position()));
2210 if (dz > dzMax)
2211 continue;
2212
2213 float dxy = fabs(pfc.trackRef()->dxy(vtx->position()));
2214 if (fabs(dxy) > dxyMax)
2215 continue;
2216 float dR = deltaR(
2217 photon_directionWrtVtx.Eta(), photon_directionWrtVtx.Phi(), pfc.momentum().Eta(), pfc.momentum().Phi());
2218 if (dR > dRmax || dR < dRveto)
2219 continue;
2220 sum += pfc.pt();
2221 if (debug_ > 1)
2222 edm::LogVerbatim("GammaJetAnalysis") << "pt=" << pfc.pt();
2223 }
2224 }
2225 if (debug_ > 1)
2226 edm::LogVerbatim("GammaJetAnalysis") << "sum=" << sum;
2227 sum = sum * 1.0;
2228 result.push_back(sum);
2229 }
2230 if (debug_ > 1) {
2231 edm::LogVerbatim("GammaJetAnalysis") << "Will return result";
2232 edm::LogVerbatim("GammaJetAnalysis") << "result" << &result;
2233 edm::LogVerbatim("GammaJetAnalysis") << "Result returned";
2234 }
2235 return result;
2236 }
2237
2238
2239
2240 void GammaJetAnalysis::clear_leadingPfJetVars() {
2241 ppfjet_pt_ = ppfjet_p_ = ppfjet_E_ = 0;
2242 ppfjet_eta_ = ppfjet_phi_ = ppfjet_scale_ = 0.;
2243 ppfjet_area_ = ppfjet_E_NPVcorr_ = 0.;
2244 ppfjet_NeutralHadronFrac_ = ppfjet_NeutralEMFrac_ = 0.;
2245 ppfjet_nConstituents_ = 0;
2246 ppfjet_ChargedHadronFrac_ = ppfjet_ChargedMultiplicity_ = 0;
2247 ppfjet_ChargedEMFrac_ = 0.;
2248 ppfjet_gendr_ = ppfjet_genpt_ = ppfjet_genp_ = ppfjet_genE_ = 0.;
2249
2250 ppfjet_unkown_E_ = ppfjet_unkown_px_ = ppfjet_unkown_py_ = ppfjet_unkown_pz_ = ppfjet_unkown_EcalE_ = 0.0;
2251 ppfjet_electron_E_ = ppfjet_electron_px_ = ppfjet_electron_py_ = ppfjet_electron_pz_ = ppfjet_electron_EcalE_ = 0.0;
2252 ppfjet_muon_E_ = ppfjet_muon_px_ = ppfjet_muon_py_ = ppfjet_muon_pz_ = ppfjet_muon_EcalE_ = 0.0;
2253 ppfjet_photon_E_ = ppfjet_photon_px_ = ppfjet_photon_py_ = ppfjet_photon_pz_ = ppfjet_photon_EcalE_ = 0.0;
2254 ppfjet_unkown_n_ = ppfjet_electron_n_ = ppfjet_muon_n_ = ppfjet_photon_n_ = 0;
2255 ppfjet_had_n_ = 0;
2256 ppfjet_ntwrs_ = 0;
2257 ppfjet_cluster_n_ = 0;
2258 ppfjet_ncandtracks_ = 0;
2259
2260 ppfjet_had_E_.clear();
2261 ppfjet_had_px_.clear();
2262 ppfjet_had_py_.clear();
2263 ppfjet_had_pz_.clear();
2264 ppfjet_had_EcalE_.clear();
2265 ppfjet_had_rawHcalE_.clear();
2266 ppfjet_had_emf_.clear();
2267 ppfjet_had_E_mctruth_.clear();
2268 ppfjet_had_id_.clear();
2269 ppfjet_had_candtrackind_.clear();
2270 ppfjet_had_mcpdgId_.clear();
2271 ppfjet_had_ntwrs_.clear();
2272 ppfjet_twr_ieta_.clear();
2273 ppfjet_twr_iphi_.clear();
2274 ppfjet_twr_depth_.clear();
2275 ppfjet_twr_subdet_.clear();
2276 ppfjet_twr_candtrackind_.clear();
2277 ppfjet_twr_hadind_.clear();
2278 ppfjet_twr_elmttype_.clear();
2279 ppfjet_twr_hade_.clear();
2280 ppfjet_twr_frac_.clear();
2281 ppfjet_twr_dR_.clear();
2282 ppfjet_twr_clusterind_.clear();
2283 ppfjet_cluster_eta_.clear();
2284 ppfjet_cluster_phi_.clear();
2285 ppfjet_cluster_dR_.clear();
2286 ppfjet_candtrack_px_.clear();
2287 ppfjet_candtrack_py_.clear();
2288 ppfjet_candtrack_pz_.clear();
2289 ppfjet_candtrack_EcalE_.clear();
2290 }
2291
2292
2293
2294 void GammaJetAnalysis::copy_leadingPfJetVars_to_pfJet2() {
2295 pfjet2_pt_ = ppfjet_pt_;
2296 pfjet2_p_ = ppfjet_p_;
2297 pfjet2_E_ = ppfjet_E_;
2298 pfjet2_eta_ = ppfjet_eta_;
2299 pfjet2_phi_ = ppfjet_phi_;
2300 pfjet2_scale_ = ppfjet_scale_;
2301 pfjet2_area_ = ppfjet_area_;
2302 pfjet2_E_NPVcorr_ = ppfjet_E_NPVcorr_;
2303 pfjet2_NeutralHadronFrac_ = ppfjet_NeutralHadronFrac_;
2304 pfjet2_NeutralEMFrac_ = ppfjet_NeutralEMFrac_;
2305 pfjet2_nConstituents_ = ppfjet_nConstituents_;
2306 pfjet2_ChargedHadronFrac_ = ppfjet_ChargedHadronFrac_;
2307 pfjet2_ChargedMultiplicity_ = ppfjet_ChargedMultiplicity_;
2308 pfjet2_ChargedEMFrac_ = ppfjet_ChargedEMFrac_;
2309
2310 pfjet2_gendr_ = ppfjet_gendr_;
2311 pfjet2_genpt_ = ppfjet_genpt_;
2312 pfjet2_genp_ = ppfjet_genp_;
2313 pfjet2_genE_ = ppfjet_genE_;
2314
2315 pfjet2_unkown_E_ = ppfjet_unkown_E_;
2316 pfjet2_unkown_px_ = ppfjet_unkown_px_;
2317 pfjet2_unkown_py_ = ppfjet_unkown_py_;
2318 pfjet2_unkown_pz_ = ppfjet_unkown_pz_;
2319 pfjet2_unkown_EcalE_ = ppfjet_unkown_EcalE_;
2320
2321 pfjet2_electron_E_ = ppfjet_electron_E_;
2322 pfjet2_electron_px_ = ppfjet_electron_px_;
2323 pfjet2_electron_py_ = ppfjet_electron_py_;
2324 pfjet2_electron_pz_ = ppfjet_electron_pz_;
2325 pfjet2_electron_EcalE_ = ppfjet_electron_EcalE_;
2326
2327 pfjet2_muon_E_ = ppfjet_muon_E_;
2328 pfjet2_muon_px_ = ppfjet_muon_px_;
2329 pfjet2_muon_py_ = ppfjet_muon_py_;
2330 pfjet2_muon_pz_ = ppfjet_muon_pz_;
2331 pfjet2_muon_EcalE_ = ppfjet_muon_EcalE_;
2332
2333 pfjet2_photon_E_ = ppfjet_photon_E_;
2334 pfjet2_photon_px_ = ppfjet_photon_px_;
2335 pfjet2_photon_py_ = ppfjet_photon_py_;
2336 pfjet2_photon_pz_ = ppfjet_photon_pz_;
2337 pfjet2_photon_EcalE_ = ppfjet_photon_EcalE_;
2338
2339 pfjet2_unkown_n_ = ppfjet_unkown_n_;
2340 pfjet2_electron_n_ = ppfjet_electron_n_;
2341 pfjet2_muon_n_ = ppfjet_muon_n_;
2342 pfjet2_photon_n_ = ppfjet_photon_n_;
2343 pfjet2_had_n_ = ppfjet_had_n_;
2344
2345 pfjet2_had_E_ = ppfjet_had_E_;
2346 pfjet2_had_px_ = ppfjet_had_px_;
2347 pfjet2_had_py_ = ppfjet_had_py_;
2348 pfjet2_had_pz_ = ppfjet_had_pz_;
2349 pfjet2_had_EcalE_ = ppfjet_had_EcalE_;
2350 pfjet2_had_rawHcalE_ = ppfjet_had_rawHcalE_;
2351 pfjet2_had_emf_ = ppfjet_had_emf_;
2352 pfjet2_had_E_mctruth_ = ppfjet_had_E_mctruth_;
2353
2354 pfjet2_had_id_ = ppfjet_had_id_;
2355 pfjet2_had_candtrackind_ = ppfjet_had_candtrackind_;
2356 pfjet2_had_mcpdgId_ = ppfjet_had_mcpdgId_;
2357 pfjet2_had_ntwrs_ = ppfjet_had_ntwrs_;
2358
2359 pfjet2_ntwrs_ = ppfjet_ntwrs_;
2360 pfjet2_twr_ieta_ = ppfjet_twr_ieta_;
2361 pfjet2_twr_iphi_ = ppfjet_twr_iphi_;
2362 pfjet2_twr_depth_ = ppfjet_twr_depth_;
2363 pfjet2_twr_subdet_ = ppfjet_twr_subdet_;
2364 pfjet2_twr_candtrackind_ = ppfjet_twr_candtrackind_;
2365 pfjet2_twr_hadind_ = ppfjet_twr_hadind_;
2366 pfjet2_twr_elmttype_ = ppfjet_twr_elmttype_;
2367 pfjet2_twr_clusterind_ = ppfjet_twr_clusterind_;
2368
2369 pfjet2_twr_hade_ = ppfjet_twr_hade_;
2370 pfjet2_twr_frac_ = ppfjet_twr_frac_;
2371 pfjet2_twr_dR_ = ppfjet_twr_dR_;
2372
2373 pfjet2_cluster_n_ = ppfjet_cluster_n_;
2374 pfjet2_cluster_eta_ = ppfjet_cluster_eta_;
2375 pfjet2_cluster_phi_ = ppfjet_cluster_phi_;
2376 pfjet2_cluster_dR_ = ppfjet_cluster_dR_;
2377
2378 pfjet2_ncandtracks_ = ppfjet_ncandtracks_;
2379 pfjet2_candtrack_px_ = ppfjet_candtrack_px_;
2380 pfjet2_candtrack_py_ = ppfjet_candtrack_py_;
2381 pfjet2_candtrack_pz_ = ppfjet_candtrack_pz_;
2382 pfjet2_candtrack_EcalE_ = ppfjet_candtrack_EcalE_;
2383 }
2384
2385
2386
2387 double GammaJetAnalysis::deltaR(const reco::Jet* j1, const reco::Jet* j2) {
2388 double deta = j1->eta() - j2->eta();
2389 double dphi = std::fabs(j1->phi() - j2->phi());
2390 if (dphi > 3.1415927)
2391 dphi = 2 * 3.1415927 - dphi;
2392 return std::sqrt(deta * deta + dphi * dphi);
2393 }
2394
2395
2396
2397 double GammaJetAnalysis::deltaR(const double eta1, const double phi1, const double eta2, const double phi2) {
2398 double deta = eta1 - eta2;
2399 double dphi = std::fabs(phi1 - phi2);
2400 if (dphi > 3.1415927)
2401 dphi = 2 * 3.1415927 - dphi;
2402 return std::sqrt(deta * deta + dphi * dphi);
2403 }
2404
2405
2406
2407
2408
2409
2410
2411
2412
2413
2414
2415
2416
2417
2418
2419
2420 int GammaJetAnalysis::getEtaPhi(const DetId id) {
2421 return id.rawId() & 0x3FFF;
2422 }
2423
2424
2425
2426 int GammaJetAnalysis::getEtaPhi(const HcalDetId id) {
2427 return id.rawId() & 0x3FFF;
2428 }
2429
2430
2431
2432
2433
2434 DEFINE_FWK_MODULE(GammaJetAnalysis);
