Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:34:44 

0001 // system include files
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 // user include files
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 // class declarations
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;  // index in the photon collection
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 // Main class
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;  //(const edm::EventSetup&);
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   // parameters
0180   const int debug_;  // print debug statements
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_;       // label for the photon collection
0188   const std::string pfJetCollName_;        // label for the PF jet collection
0189   const std::string genJetCollName_;       // label for the genjet collection
0190   const std::string genParticleCollName_;  // label for the genparticle collection
0191   const std::string genEventInfoName_;     // label for the generator event info collection
0192   const std::string hbheRecHitName_;       // label for HBHERecHits collection
0193   const std::string hfRecHitName_;         // label for HFRecHit collection
0194   const std::string hoRecHitName_;         // label for HORecHit collection
0195   const std::string rootHistFilename_;     // name of the histogram file
0196   const std::string pvCollName_;           // label for primary vertex collection
0197   const std::string prodProcess_;          // the producer process for AOD=2
0198 
0199   const bool allowNoPhoton_;                   // whether module is used for dijet analysis
0200   const double photonPtMin_;                   // lowest value of the leading photon pT
0201   const double photonJetDPhiMin_;              // phi angle between the leading photon and the leading jet
0202   const double jetEtMin_;                      // lowest value of the leading jet ET
0203   const double jet2EtMax_;                     // largest value of the subleading jet ET
0204   const double jet3EtMax_;                     // largest value of the third jet ET
0205   std::vector<std::string> photonTrigNamesV_;  // photon trigger names
0206   std::vector<std::string> jetTrigNamesV_;     // jet trigger names
0207   const bool writeTriggerPrescale_;            // whether attempt to record the prescale
0208 
0209   bool doPFJets_;   // use PFJets
0210   bool doGenJets_;  // use GenJets
0211   int workOnAOD_;
0212   bool ignoreHLT_;
0213 
0214   //Tokens
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   // root file/histograms
0241   TTree* misc_tree_;  // misc.information. Will be filled only once
0242   TTree* calo_tree_;
0243   TTree* pf_tree_;
0244 
0245   // trigger info
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   // Event info
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   /// MET info
0261   float met_value_, met_phi_, met_sumEt_;
0262   float metType1_value_, metType1_phi_, metType1_sumEt_;
0263 
0264   // photon info
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   // Particle-flow jets
0278   // leading Et jet info
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   // subleading Et jet info
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   // helper functions
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   // set parameters
0525 
0526   eventWeight_ = 1.0;
0527   eventPtHat_ = 0.;
0528   nProcessed_ = 0;
0529 
0530   //Get the tokens
0531   // FAST FIX
0532   if (workOnAOD_ < 2) {  // origin data file
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     // FAST FIX
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     //tok_loosePhoton_ = consumes<edm::ValueMap<Bool_t> >(edm::InputTag("PhotonIDProdGED","PhotonCutBasedIDLoose"));
0576     //tok_tightPhoton_ = consumes<edm::ValueMap<Bool_t> >(edm::InputTag("PhotonIDProdGED:PhotonCutBasedIDTight"));
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 // member functions
0597 //
0598 
0599 // ------------ method called to for each event  ------------
0600 void GammaJetAnalysis::analyze(const edm::Event& iEvent, const edm::EventSetup& evSetup) {
0601   nProcessed_++;
0602 
0603   edm::LogVerbatim("GammaJetAnalysis") << "nProcessed=" << nProcessed_ << "\n";
0604 
0605   // 1st. Get Photons //
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   // Get photon quality flags
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   // sort photons by Et //
0641   // counter is needed later to get the reference to the ptr
0642   std::set<PhotonPair, PhotonPairComp> photonpairset;
0643   int counter = 0;
0644   for (reco::PhotonCollection::const_iterator it = photons->begin(); it != photons->end(); ++it) {
0645     //HERE(Form("photon counter=%d",counter));
0646     const reco::Photon* photon = &(*it);
0647     //if(loosePhotonQual.isValid()){
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   // TAG = Highest Et photon
0663   ///////////////////////////////
0664 
0665   // find highest Et photon //
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   // cut on photon pt
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   // 2nd. Get Jets
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   // 3rd. Check the trigger
0727   photonTrigFired_.clear();
0728   photonTrigPrescale_.clear();
0729   jetTrigFired_.clear();
0730   jetTrigPrescale_.clear();
0731 
0732   HERE("trigger");
0733 
0734   // HLT Trigger
0735   // assign "trig fired" if no triggers are specified
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   // If needed, process trigger information
0744   if (!photonTrigFlag || !jetTrigFlag) {
0745     // check the triggers
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         // for triggers with two L1 seeds this fails
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   //rho2012_ = -1e6;
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     // Get GenJets
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     // Get GenParticles
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     // Get weights
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   // fill tag photon variables
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     //tagPho_ConvSafeEleVeto_ = ((int)ConversionTools::hasMatchedPromptElectron(photon_tag.photon()->superCluster(), gsfElectronHandle, convH, beamSpotHandle->position()));
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   // Run over PFJets //
0972 
0973   if (doPFJets_ && (nPFJets_ > 0)) {
0974     // Get RecHits in HB and HE
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     // Get RecHits in HF
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     // Get RecHits in HO
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     // Get geometry
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     // Get primary vertices
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     // sort jets by corrected et
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       // do not let the jet to be close to the tag photon
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       //else break; // don't break for the statistics
1066     }
1067 
1068     HERE("reached selection");
1069 
1070     // Check selection
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       // put values into 3rd jet quantities
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       // Get PF constituents and fill HCAL towers
1112       /////////////////////////////////////////////
1113 
1114       // fill jet variables
1115       // First start from a second jet, then fill the first jet
1116       PFJetCorretPair pfjet_probe_store = pfjet_probe;
1117       for (int iJet = 2; iJet > 0; iJet--) {
1118         // prepare the container
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             // put zeros into 2nd jet quantities
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         // temporary variables
1139         std::map<int, std::pair<int, std::set<float>>> ppfjet_rechits;
1140         std::map<float, int> ppfjet_clusters;
1141 
1142         // fill the values
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         // Get PF constituents and fill HCAL towers
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           // store information
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             // Get elements from block
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) {  // Element is HB or HE
1377                   // Get cluster and hits
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                   // Run over hits and match
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                       }                                                           // Test if ieta,iphi matches
1463                     }                                                             // Loop over rechits
1464                   }                                                               // Loop over hits
1465                 }                                                                 // Test if element is from HCAL
1466                 else if (elements[iEle].type() == reco::PFBlockElement::HFHAD) {  // Element is HF
1467                   ////  h_etaHFHAD_->Fill((*it)->eta());
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;  // Remove long fibers
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) {  // Element is HF
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;  // Remove short fibers
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) {  // Element is 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                       }  // Test if ieta,iphi match
1626                     }    // Loop over rechits
1627                   }      // Loop over hits
1628                 }        // Test if element is from HO
1629               }          // Test for right element index
1630             }            // Loop over elements
1631           }              // Loop over elements in blocks
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         }  // Loop over PF constitutents
1644 
1645         if (doGenJets_) {
1646           // fill genjet variables
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         }  // doGenJets_
1661         if (iJet == 2) {
1662           copy_leadingPfJetVars_to_pfJet2();
1663         }
1664       }
1665       // double jer= (ppfjet_genpt_==double(0)) ?
1666       //  0. : pfjet_probe.jet()->et()/ppfjet_genpt_;
1667       ///h_pfrecoOgen_et->Fill(jer, eventWeight_);
1668 
1669       ///// MET /////
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         // do we need an exception here?
1686         metType1_value_ = -999.;
1687         metType1_phi_ = -999.;
1688         metType1_sumEt_ = -999.;
1689       }
1690 
1691       // fill photon+jet variables
1692       pf_tree_->Fill();
1693     }
1694   }
1695   return;
1696 }
1697 
1698 // ------------ method called once each job just before starting event loop  ------------
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;  // ! caloJet are no longer there, so store only once
1711     TTree* tree = pf_tree_;
1712 
1713     // Event triggers
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     // Event info
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     // Photon info
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     // gen.info on photon
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     // counters
1756     tree->Branch("nPhotons", &nPhotons_, "nPhotons/I");
1757     tree->Branch("nGenJets", &nGenJets_, "nGenJets/I");
1758   }
1759 
1760   //////// Particle Flow ////////
1761 
1762   if (doPFJets_) {
1763     pf_tree_->Branch("nPFJets", &nPFJets_, "nPFJets/I");
1764 
1765     // Leading jet info
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     // Subleading jet info
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     // third pf jet
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 // ------------ method called once each job just after ending the event loop  ------------
1953 void GammaJetAnalysis::endJob() {
1954   if (doPFJets_) {
1955     pf_tree_->Write();
1956   }
1957   // write miscItems
1958   // Save info about the triggers and other misc items
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     // put time stamp
1970     time_t ltime;
1971     ltime = time(NULL);
1972     TString str = TString(asctime(localtime(&ltime)));
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     // clear trigger names, if needed
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       // if init returns TRUE, initialisation has succeeded!
2013       if (changed) {
2014         // The HLT config has actually changed wrt the previous Run, hence rebook your
2015         // histograms or do anything else dependent on the revised HLT config
2016       }
2017     } else {
2018       // if init returns FALSE, initialisation has NOT succeeded, which indicates a problem
2019       // with the file and/or code and needs to be investigated!
2020       throw edm::Exception(edm::errors::ProductNotFound)
2021           << " HLT config extraction failure with process name " << processName;
2022       // In this case, all access methods will return empty values!
2023     }
2024   }
2025 }
2026 
2027 // ---------------------------------------------------------------------
2028 
2029 // helper function
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       // Do not include the PFCandidate associated by SC Ref to the reco::Photon
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       // Shift the photon direction vector according to the PF vertex
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   //Calculate isolation sum separately for each vertex
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     // Shift the photon according to the vertex
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     // Loop over the PFCandidates
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       //require that PFCandidate is a charged hadron
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   // Reset particle variables
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 // DetId rawId bits xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2366 //                  1111222      3333345555556666666
2367 //   1 = detector
2368 //   2 = subdetector
2369 //   3 = depth
2370 //   4 = zside: 0 = negative z, 1 = positive z \
2371 //   5 = abs(ieta)                              | ieta,iphi
2372 //   6 = abs(iphi)                             /
2373 */
2374 
2375 // ---------------------------------------------------------------------
2376 
2377 int GammaJetAnalysis::getEtaPhi(const DetId id) {
2378   return id.rawId() & 0x3FFF;  // Get 14 least-significant digits
2379 }
2380 
2381 // ---------------------------------------------------------------------
2382 
2383 int GammaJetAnalysis::getEtaPhi(const HcalDetId id) {
2384   return id.rawId() & 0x3FFF;  // Get 14 least-significant digits
2385 }
2386 
2387 // ---------------------------------------------------------------------
2388 
2389 //define this as a plug-in
2390 
2391 DEFINE_FWK_MODULE(GammaJetAnalysis);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!