Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-15 23:49:19 

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