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