Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-28 01:32:06

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