Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-09-21 02:12:12

0001 // -*- C++ -*-//
0002 // Package:    IsoTrig
0003 // Class:      IsoTrig
0004 //
0005 /**\class IsoTrig IsoTrig.cc IsoTrig/IsoTrig/src/IsoTrig.cc
0006 
0007  Description: [one line class summary]
0008 
0009  Implementation:
0010      [Notes on implementation]
0011 */
0012 //
0013 // Original Author:  Ruchi Gupta
0014 //         Created:  Fri May 25 12:02:48 CDT 2012
0015 // $Id$
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <unordered_map>
0022 
0023 // Root objects
0024 #include "TROOT.h"
0025 #include "TH1.h"
0026 #include "TH2.h"
0027 #include "TSystem.h"
0028 #include "TFile.h"
0029 #include "TProfile.h"
0030 #include "TDirectory.h"
0031 #include "TTree.h"
0032 #include "TMath.h"
0033 
0034 // user include files
0035 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0036 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0037 #include "Calibration/IsolatedParticles/interface/eCone.h"
0038 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0039 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0040 
0041 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0042 
0043 #include "DataFormats/Common/interface/RefToBase.h"
0044 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0045 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0046 #include "DataFormats/Math/interface/LorentzVector.h"
0047 #include "DataFormats/Math/interface/Point3D.h"
0048 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0049 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0050 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0051 #include "DataFormats/TrackReco/interface/Track.h"
0052 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0053 #include "DataFormats/TrackReco/interface/TrackBase.h"
0054 //Tracks
0055 #include "DataFormats/TrackReco/interface/HitPattern.h"
0056 // Vertices
0057 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0058 #include "DataFormats/VertexReco/interface/Vertex.h"
0059 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0060 //Triggers
0061 #include "DataFormats/Common/interface/TriggerResults.h"
0062 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0063 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0064 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0065 #include "DataFormats/Luminosity/interface/LumiDetails.h"
0066 
0067 #include "FWCore/Common/interface/TriggerNames.h"
0068 #include "FWCore/Framework/interface/Frameworkfwd.h"
0069 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0070 #include "FWCore/Framework/interface/Event.h"
0071 #include "FWCore/Framework/interface/MakerMacros.h"
0072 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0073 #include "FWCore/ServiceRegistry/interface/Service.h"
0074 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0075 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0076 
0077 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0078 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0079 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0080 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0081 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0082 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0083 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0084 
0085 #include "HLTrigger/HLTcore/interface/HLTPrescaleProvider.h"
0086 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0087 
0088 #include "MagneticField/Engine/interface/MagneticField.h"
0089 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0090 
0091 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0092 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0093 
0094 #include "TrackingTools/TransientTrackingRecHit/interface/SeedingLayerSetsHits.h"
0095 
0096 class IsoTrig : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0097 public:
0098   explicit IsoTrig(const edm::ParameterSet &);
0099   ~IsoTrig() override;
0100 
0101   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0102 
0103 private:
0104   void analyze(const edm::Event &, const edm::EventSetup &) override;
0105   void beginJob() override;
0106   void endJob() override;
0107   void beginRun(edm::Run const &, edm::EventSetup const &) override;
0108   void endRun(edm::Run const &, edm::EventSetup const &) override {}
0109 
0110   void clearMipCutTreeVectors();
0111   void clearChgIsolnTreeVectors();
0112   void pushChgIsolnTreeVecs(math::XYZTLorentzVector &Pixcand,
0113                             math::XYZTLorentzVector &Trkcand,
0114                             std::vector<double> &PixMaxP,
0115                             double &TrkMaxP,
0116                             bool &selTk);
0117   void pushMipCutTreeVecs(math::XYZTLorentzVector &NFcand,
0118                           math::XYZTLorentzVector &Trkcand,
0119                           double &EmipNFcand,
0120                           double &EmipTrkcand,
0121                           double &mindR,
0122                           double &mindP1,
0123                           std::vector<bool> &Flags,
0124                           double hCone);
0125   void StudyTrkEbyP(edm::Handle<reco::TrackCollection> &trkCollection);
0126   void studyTiming(const edm::Event &theEvent);
0127   void studyMipCut(edm::Handle<reco::TrackCollection> &trkCollection,
0128                    edm::Handle<reco::IsolatedPixelTrackCandidateCollection> &L2cands);
0129   void studyTrigger(edm::Handle<reco::TrackCollection> &, std::vector<reco::TrackCollection::const_iterator> &);
0130   void studyIsolation(edm::Handle<reco::TrackCollection> &, std::vector<reco::TrackCollection::const_iterator> &);
0131   void chgIsolation(double &etaTriggered,
0132                     double &phiTriggered,
0133                     edm::Handle<reco::TrackCollection> &trkCollection,
0134                     const edm::Event &theEvent);
0135   void getGoodTracks(const edm::Event &, edm::Handle<reco::TrackCollection> &);
0136   void fillHist(int, math::XYZTLorentzVector &);
0137   void fillDifferences(int, math::XYZTLorentzVector &, math::XYZTLorentzVector &, bool);
0138   void fillCuts(int, double, double, double, math::XYZTLorentzVector &, int, bool);
0139   void fillEnergy(int, int, double, double, math::XYZTLorentzVector &);
0140   double dEta(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0141   double dPhi(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0142   double dR(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0143   double dPt(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0144   double dP(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0145   double dinvPt(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0146   std::pair<double, double> etaPhiTrigger();
0147   std::pair<double, double> GetEtaPhiAtEcal(double etaIP, double phiIP, double pT, int charge, double vtxZ);
0148   double getDistInCM(double eta1, double phi1, double eta2, double phi2);
0149 
0150   // ----------member data ---------------------------
0151   HLTPrescaleProvider hltPrescaleProvider_;
0152   const std::vector<std::string> trigNames_;
0153   const edm::InputTag pixCandTag_, l1CandTag_, l2CandTag_;
0154   const std::vector<edm::InputTag> pixelTracksSources_;
0155   const bool doL2L3_, doTiming_, doMipCutTree_;
0156   const bool doTrkResTree_, doChgIsolTree_, doStudyIsol_;
0157   const int verbosity_;
0158   const std::vector<double> pixelIsolationConeSizeAtEC_;
0159   const double minPTrackValue_, vtxCutSeed_, vtxCutIsol_;
0160   const double tauUnbiasCone_, prelimCone_;
0161   std::string theTrackQuality_;
0162   const std::string processName_;
0163   double rEB_, zEE_, bfVal_;
0164   spr::trackSelectionParameters selectionParameters_;
0165   const double dr_L1_, a_coneR_, a_charIsoR_, a_neutIsoR_;
0166   const double a_mipR_, a_neutR1_, a_neutR2_, cutMip_;
0167   const double cutCharge_, cutNeutral_;
0168   const int minRunNo_, maxRunNo_;
0169   edm::EDGetTokenT<LumiDetails> tok_lumi_;
0170   edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
0171   edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0172   edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_hlt_;
0173   edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0174   edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0175   edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0176   edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0177   edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0178   edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0179   edm::EDGetTokenT<reco::VertexCollection> tok_verthb_, tok_verthe_;
0180   edm::EDGetTokenT<SeedingLayerSetsHits> tok_SeedingLayerHB_;
0181   edm::EDGetTokenT<SeedingLayerSetsHits> tok_SeedingLayerHE_;
0182   edm::EDGetTokenT<SiPixelRecHitCollection> tok_SiPixelRecHits_;
0183   edm::EDGetTokenT<reco::IsolatedPixelTrackCandidateCollection> tok_pixtk_;
0184   edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_l1cand_;
0185   edm::EDGetTokenT<reco::IsolatedPixelTrackCandidateCollection> tok_l2cand_;
0186   std::vector<edm::EDGetTokenT<reco::TrackCollection>> tok_pixtks_;
0187 
0188   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0189   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0190 
0191   std::vector<reco::TrackRef> pixelTrackRefsHB_, pixelTrackRefsHE_;
0192   edm::Handle<HBHERecHitCollection> hbhe_;
0193   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle_;
0194   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle_;
0195   edm::Handle<reco::BeamSpot> beamSpotH_;
0196   edm::Handle<reco::VertexCollection> recVtxs_;
0197 
0198   const MagneticField *bField_;
0199   const CaloGeometry *geo_;
0200   math::XYZPoint leadPV_;
0201 
0202   std::unordered_map<unsigned int, unsigned int> trigList_;
0203   std::unordered_map<unsigned int, const std::pair<double, double>> trigPreList_;
0204   bool changed_;
0205   double pLimits_[6];
0206   edm::Service<TFileService> fs_;
0207   TTree *MipCutTree_, *ChgIsolnTree_, *TrkResTree_, *TimingTree_;
0208   std::vector<double> *t_timeL2Prod;
0209   std::vector<int> *t_nPixCand;
0210   std::vector<int> *t_nPixSeed;
0211   std::vector<int> *t_nGoodTk;
0212 
0213   std::vector<double> *t_TrkhCone;
0214   std::vector<double> *t_TrkP;
0215   std::vector<bool> *t_TrkselTkFlag;
0216   std::vector<bool> *t_TrkqltyFlag;
0217   std::vector<bool> *t_TrkMissFlag;
0218   std::vector<bool> *t_TrkPVFlag;
0219   std::vector<bool> *t_TrkNuIsolFlag;
0220 
0221   std::vector<double> *t_PixcandP;
0222   std::vector<double> *t_PixcandPt;
0223   std::vector<double> *t_PixcandEta;
0224   std::vector<double> *t_PixcandPhi;
0225   std::vector<std::vector<double>> *t_PixcandMaxP;
0226   std::vector<double> *t_PixTrkcandP;
0227   std::vector<double> *t_PixTrkcandPt;
0228   std::vector<double> *t_PixTrkcandEta;
0229   std::vector<double> *t_PixTrkcandPhi;
0230   std::vector<double> *t_PixTrkcandMaxP;
0231   std::vector<bool> *t_PixTrkcandselTk;
0232 
0233   std::vector<double> *t_NFcandP;
0234   std::vector<double> *t_NFcandPt;
0235   std::vector<double> *t_NFcandEta;
0236   std::vector<double> *t_NFcandPhi;
0237   std::vector<double> *t_NFcandEmip;
0238   std::vector<double> *t_NFTrkcandP;
0239   std::vector<double> *t_NFTrkcandPt;
0240   std::vector<double> *t_NFTrkcandEta;
0241   std::vector<double> *t_NFTrkcandPhi;
0242   std::vector<double> *t_NFTrkcandEmip;
0243   std::vector<double> *t_NFTrkMinDR;
0244   std::vector<double> *t_NFTrkMinDP1;
0245   std::vector<bool> *t_NFTrkselTkFlag;
0246   std::vector<bool> *t_NFTrkqltyFlag;
0247   std::vector<bool> *t_NFTrkMissFlag;
0248   std::vector<bool> *t_NFTrkPVFlag;
0249   std::vector<bool> *t_NFTrkPropFlag;
0250   std::vector<bool> *t_NFTrkChgIsoFlag;
0251   std::vector<bool> *t_NFTrkNeuIsoFlag;
0252   std::vector<bool> *t_NFTrkMipFlag;
0253   std::vector<double> *t_ECone;
0254 
0255   TH1D *h_EnIn, *h_EnOut;
0256   TH2D *h_MipEnMatch, *h_MipEnNoMatch;
0257   TH1I *h_nHLT, *h_HLT, *h_PreL1, *h_PreHLT;
0258   TH1I *h_Pre, *h_nL3Objs, *h_Filters;
0259   TH1D *h_PreL1wt, *h_PreHLTwt, *h_L1ObjEnergy;
0260   TH1D *h_p[20], *h_pt[20], *h_eta[20], *h_phi[20];
0261   TH1D *h_dEtaL1[2], *h_dPhiL1[2], *h_dRL1[2];
0262   TH1D *h_dEta[9], *h_dPhi[9], *h_dPt[9], *h_dP[9];
0263   TH1D *h_dinvPt[9], *h_mindR[9], *h_eMip[2];
0264   TH1D *h_eMaxNearP[2], *h_eNeutIso[2];
0265   TH1D *h_etaCalibTracks[5][2][2], *h_etaMipTracks[5][2][2];
0266   TH1D *h_eHcal[5][6][48], *h_eCalo[5][6][48];
0267   TH1D *g_Pre, *g_PreL1, *g_PreHLT;
0268   TH1I *g_Accepts;
0269   std::vector<math::XYZTLorentzVector> vec_[3];
0270 };
0271 
0272 IsoTrig::IsoTrig(const edm::ParameterSet &iConfig)
0273     : hltPrescaleProvider_(iConfig, consumesCollector(), *this),
0274       trigNames_(iConfig.getUntrackedParameter<std::vector<std::string>>("Triggers")),
0275       pixCandTag_(iConfig.getUntrackedParameter<edm::InputTag>("pixCandTag")),
0276       l1CandTag_(iConfig.getUntrackedParameter<edm::InputTag>("l1CandTag")),
0277       l2CandTag_(iConfig.getUntrackedParameter<edm::InputTag>("l2CandTag")),
0278       pixelTracksSources_(iConfig.getUntrackedParameter<std::vector<edm::InputTag>>("pixelTracksSources")),
0279       doL2L3_(iConfig.getUntrackedParameter<bool>("doL2L3", true)),
0280       doTiming_(iConfig.getUntrackedParameter<bool>("doTimingTree", true)),
0281       doMipCutTree_(iConfig.getUntrackedParameter<bool>("doMipCutTree", true)),
0282       doTrkResTree_(iConfig.getUntrackedParameter<bool>("doTrkResTree", true)),
0283       doChgIsolTree_(iConfig.getUntrackedParameter<bool>("doChgIsolTree", true)),
0284       doStudyIsol_(iConfig.getUntrackedParameter<bool>("doStudyIsol", true)),
0285       verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0286       pixelIsolationConeSizeAtEC_(iConfig.getUntrackedParameter<std::vector<double>>("pixelIsolationConeSizeAtEC")),
0287       minPTrackValue_(iConfig.getUntrackedParameter<double>("minPTrackValue")),
0288       vtxCutSeed_(iConfig.getUntrackedParameter<double>("vertexCutSeed")),
0289       vtxCutIsol_(iConfig.getUntrackedParameter<double>("vertexCutIsol")),
0290       tauUnbiasCone_(iConfig.getUntrackedParameter<double>("tauUnbiasCone")),
0291       prelimCone_(iConfig.getUntrackedParameter<double>("prelimCone")),
0292       theTrackQuality_(iConfig.getUntrackedParameter<std::string>("trackQuality", "highPurity")),
0293       processName_(iConfig.getUntrackedParameter<std::string>("processName", "HLT")),
0294       dr_L1_(iConfig.getUntrackedParameter<double>("isolationL1", 1.0)),
0295       a_coneR_(iConfig.getUntrackedParameter<double>("coneRadius", 34.98)),
0296       a_charIsoR_(a_coneR_ + 28.9),
0297       a_neutIsoR_(a_charIsoR_ * 0.726),
0298       a_mipR_(iConfig.getUntrackedParameter<double>("coneRadiusMIP", 14.0)),
0299       a_neutR1_(iConfig.getUntrackedParameter<double>("coneRadiusNeut1", 21.0)),
0300       a_neutR2_(iConfig.getUntrackedParameter<double>("coneRadiusNeut2", 29.0)),
0301       cutMip_(iConfig.getUntrackedParameter<double>("cutMIP", 1.0)),
0302       cutCharge_(iConfig.getUntrackedParameter<double>("chargeIsolation", 2.0)),
0303       cutNeutral_(iConfig.getUntrackedParameter<double>("neutralIsolation", 2.0)),
0304       minRunNo_(iConfig.getUntrackedParameter<int>("minRun")),
0305       maxRunNo_(iConfig.getUntrackedParameter<int>("maxRun")),
0306       changed_(false),
0307       t_timeL2Prod(nullptr),
0308       t_nPixCand(nullptr),
0309       t_nPixSeed(nullptr),
0310       t_nGoodTk(nullptr),
0311       t_TrkhCone(nullptr),
0312       t_TrkP(nullptr),
0313       t_TrkselTkFlag(nullptr),
0314       t_TrkqltyFlag(nullptr),
0315       t_TrkMissFlag(nullptr),
0316       t_TrkPVFlag(nullptr),
0317       t_TrkNuIsolFlag(nullptr),
0318       t_PixcandP(nullptr),
0319       t_PixcandPt(nullptr),
0320       t_PixcandEta(nullptr),
0321       t_PixcandPhi(nullptr),
0322       t_PixcandMaxP(nullptr),
0323       t_PixTrkcandP(nullptr),
0324       t_PixTrkcandPt(nullptr),
0325       t_PixTrkcandEta(nullptr),
0326       t_PixTrkcandPhi(nullptr),
0327       t_PixTrkcandMaxP(nullptr),
0328       t_PixTrkcandselTk(nullptr),
0329       t_NFcandP(nullptr),
0330       t_NFcandPt(nullptr),
0331       t_NFcandEta(nullptr),
0332       t_NFcandPhi(nullptr),
0333       t_NFcandEmip(nullptr),
0334       t_NFTrkcandP(nullptr),
0335       t_NFTrkcandPt(nullptr),
0336       t_NFTrkcandEta(nullptr),
0337       t_NFTrkcandPhi(nullptr),
0338       t_NFTrkcandEmip(nullptr),
0339       t_NFTrkMinDR(nullptr),
0340       t_NFTrkMinDP1(nullptr),
0341       t_NFTrkselTkFlag(nullptr),
0342       t_NFTrkqltyFlag(nullptr),
0343       t_NFTrkMissFlag(nullptr),
0344       t_NFTrkPVFlag(nullptr),
0345       t_NFTrkPropFlag(nullptr),
0346       t_NFTrkChgIsoFlag(nullptr),
0347       t_NFTrkNeuIsoFlag(nullptr),
0348       t_NFTrkMipFlag(nullptr),
0349       t_ECone(nullptr) {
0350   usesResource(TFileService::kSharedResource);
0351 
0352   //now do whatever initialization is needed
0353   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality_);
0354   selectionParameters_.minPt = iConfig.getUntrackedParameter<double>("minTrackPt", 10.0);
0355   selectionParameters_.minQuality = trackQuality_;
0356   selectionParameters_.maxDxyPV = iConfig.getUntrackedParameter<double>("maxDxyPV", 0.2);
0357   selectionParameters_.maxDzPV = iConfig.getUntrackedParameter<double>("maxDzPV", 5.0);
0358   selectionParameters_.maxChi2 = iConfig.getUntrackedParameter<double>("maxChi2", 5.0);
0359   selectionParameters_.maxDpOverP = iConfig.getUntrackedParameter<double>("maxDpOverP", 0.1);
0360   selectionParameters_.minOuterHit = iConfig.getUntrackedParameter<int>("minOuterHit", 4);
0361   selectionParameters_.minLayerCrossed = iConfig.getUntrackedParameter<int>("minLayerCrossed", 8);
0362   selectionParameters_.maxInMiss = iConfig.getUntrackedParameter<int>("maxInMiss", 0);
0363   selectionParameters_.maxOutMiss = iConfig.getUntrackedParameter<int>("maxOutMiss", 0);
0364 
0365   // define tokens for access
0366   tok_lumi_ = consumes<LumiDetails, edm::InLumi>(edm::InputTag("lumiProducer"));
0367   edm::InputTag triggerEvent_("hltTriggerSummaryAOD", "", processName_);
0368   tok_trigEvt_ = consumes<trigger::TriggerEvent>(triggerEvent_);
0369   edm::InputTag theTriggerResultsLabel("TriggerResults", "", processName_);
0370   tok_trigRes_ = consumes<edm::TriggerResults>(theTriggerResultsLabel);
0371   tok_genTrack_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
0372   tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0373   tok_bs_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
0374   tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0375   tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0376   tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
0377   tok_pixtk_ = consumes<reco::IsolatedPixelTrackCandidateCollection>(pixCandTag_);
0378   tok_l1cand_ = consumes<trigger::TriggerFilterObjectWithRefs>(l1CandTag_);
0379   tok_l2cand_ = consumes<reco::IsolatedPixelTrackCandidateCollection>(l2CandTag_);
0380   if (doTiming_) {
0381     tok_verthb_ = consumes<reco::VertexCollection>(edm::InputTag("hltHITPixelVerticesHB"));
0382     tok_verthe_ = consumes<reco::VertexCollection>(edm::InputTag("hltHITPixelVerticesHB"));
0383     tok_hlt_ = consumes<trigger::TriggerFilterObjectWithRefs>(edm::InputTag("hltL1sL1SingleJet68"));
0384     tok_SeedingLayerHB_ = consumes<SeedingLayerSetsHits>(edm::InputTag("hltPixelLayerTripletsHITHB"));
0385     tok_SeedingLayerHE_ = consumes<SeedingLayerSetsHits>(edm::InputTag("hltPixelLayerTripletsHITHE"));
0386     tok_SiPixelRecHits_ = consumes<SiPixelRecHitCollection>(edm::InputTag("hltSiPixelRecHits"));
0387   }
0388   if (doChgIsolTree_) {
0389     for (unsigned int k = 0; k < pixelTracksSources_.size(); ++k) {
0390       //      edm::InputTag  pix (pixelTracksSources_[k],"",processName_);
0391       //      tok_pixtks_.push_back(consumes<reco::TrackCollection>(pix));
0392       tok_pixtks_.push_back(consumes<reco::TrackCollection>(pixelTracksSources_[k]));
0393     }
0394   }
0395   if (verbosity_ >= 0) {
0396     edm::LogVerbatim("IsoTrack") << "Parameters read from config file \n"
0397                                  << "\t minPt " << selectionParameters_.minPt << "\t theTrackQuality "
0398                                  << theTrackQuality_ << "\t minQuality " << selectionParameters_.minQuality
0399                                  << "\t maxDxyPV " << selectionParameters_.maxDxyPV << "\t maxDzPV "
0400                                  << selectionParameters_.maxDzPV << "\t maxChi2 " << selectionParameters_.maxChi2
0401                                  << "\t maxDpOverP " << selectionParameters_.maxDpOverP << "\t minOuterHit "
0402                                  << selectionParameters_.minOuterHit << "\t minLayerCrossed "
0403                                  << selectionParameters_.minLayerCrossed << "\t maxInMiss "
0404                                  << selectionParameters_.maxInMiss << "\t maxOutMiss "
0405                                  << selectionParameters_.maxOutMiss << "\t a_coneR " << a_coneR_ << "\t a_charIsoR "
0406                                  << a_charIsoR_ << "\t a_neutIsoR " << a_neutIsoR_ << "\t a_mipR " << a_mipR_
0407                                  << "\t a_neutR " << a_neutR1_ << ":" << a_neutR2_ << "\t cuts (MIP " << cutMip_
0408                                  << " : Charge " << cutCharge_ << " : Neutral " << cutNeutral_ << ")";
0409     edm::LogVerbatim("IsoTrack") << "Charge Isolation parameters:"
0410                                  << "\t minPTrackValue " << minPTrackValue_ << "\t vtxCutSeed " << vtxCutSeed_
0411                                  << "\t vtxCutIsol " << vtxCutIsol_ << "\t tauUnbiasCone " << tauUnbiasCone_
0412                                  << "\t prelimCone " << prelimCone_ << "\t pixelIsolationConeSizeAtEC";
0413     for (unsigned int k = 0; k < pixelIsolationConeSizeAtEC_.size(); ++k)
0414       edm::LogVerbatim("IsoTrack") << "[" << k << "] " << pixelIsolationConeSizeAtEC_[k];
0415   }
0416   double pl[] = {20, 30, 40, 60, 80, 120};
0417   for (int i = 0; i < 6; ++i)
0418     pLimits_[i] = pl[i];
0419   rEB_ = 123.8;
0420   zEE_ = 317.0;
0421 
0422   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0423   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0424 }
0425 
0426 IsoTrig::~IsoTrig() {
0427   // do anything here that needs to be done at desctruction time
0428   // (e.g. close files, deallocate resources etc.)
0429   if (t_timeL2Prod)
0430     delete t_timeL2Prod;
0431   if (t_nPixCand)
0432     delete t_nPixCand;
0433   if (t_nPixSeed)
0434     delete t_nPixSeed;
0435   if (t_nGoodTk)
0436     delete t_nGoodTk;
0437   if (t_TrkhCone)
0438     delete t_TrkhCone;
0439   if (t_TrkP)
0440     delete t_TrkP;
0441   if (t_TrkselTkFlag)
0442     delete t_TrkselTkFlag;
0443   if (t_TrkqltyFlag)
0444     delete t_TrkqltyFlag;
0445   if (t_TrkMissFlag)
0446     delete t_TrkMissFlag;
0447   if (t_TrkPVFlag)
0448     delete t_TrkPVFlag;
0449   if (t_TrkNuIsolFlag)
0450     delete t_TrkNuIsolFlag;
0451   if (t_PixcandP)
0452     delete t_PixcandP;
0453   if (t_PixcandPt)
0454     delete t_PixcandPt;
0455   if (t_PixcandEta)
0456     delete t_PixcandEta;
0457   if (t_PixcandPhi)
0458     delete t_PixcandPhi;
0459   if (t_PixcandMaxP)
0460     delete t_PixcandMaxP;
0461   if (t_PixTrkcandP)
0462     delete t_PixTrkcandP;
0463   if (t_PixTrkcandPt)
0464     delete t_PixTrkcandPt;
0465   if (t_PixTrkcandEta)
0466     delete t_PixTrkcandEta;
0467   if (t_PixTrkcandPhi)
0468     delete t_PixTrkcandPhi;
0469   if (t_PixTrkcandMaxP)
0470     delete t_PixTrkcandMaxP;
0471   if (t_PixTrkcandselTk)
0472     delete t_PixTrkcandselTk;
0473   if (t_NFcandP)
0474     delete t_NFcandP;
0475   if (t_NFcandPt)
0476     delete t_NFcandPt;
0477   if (t_NFcandEta)
0478     delete t_NFcandEta;
0479   if (t_NFcandPhi)
0480     delete t_NFcandPhi;
0481   if (t_NFcandEmip)
0482     delete t_NFcandEmip;
0483   if (t_NFTrkcandP)
0484     delete t_NFTrkcandP;
0485   if (t_NFTrkcandPt)
0486     delete t_NFTrkcandPt;
0487   if (t_NFTrkcandEta)
0488     delete t_NFTrkcandEta;
0489   if (t_NFTrkcandPhi)
0490     delete t_NFTrkcandPhi;
0491   if (t_NFTrkcandEmip)
0492     delete t_NFTrkcandEmip;
0493   if (t_NFTrkMinDR)
0494     delete t_NFTrkMinDR;
0495   if (t_NFTrkMinDP1)
0496     delete t_NFTrkMinDP1;
0497   if (t_NFTrkselTkFlag)
0498     delete t_NFTrkselTkFlag;
0499   if (t_NFTrkqltyFlag)
0500     delete t_NFTrkqltyFlag;
0501   if (t_NFTrkMissFlag)
0502     delete t_NFTrkMissFlag;
0503   if (t_NFTrkPVFlag)
0504     delete t_NFTrkPVFlag;
0505   if (t_NFTrkPropFlag)
0506     delete t_NFTrkPropFlag;
0507   if (t_NFTrkChgIsoFlag)
0508     delete t_NFTrkChgIsoFlag;
0509   if (t_NFTrkNeuIsoFlag)
0510     delete t_NFTrkNeuIsoFlag;
0511   if (t_NFTrkMipFlag)
0512     delete t_NFTrkMipFlag;
0513   if (t_ECone)
0514     delete t_ECone;
0515 }
0516 
0517 void IsoTrig::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0518   std::vector<std::string> triggers = {"HLT_IsoTrackHB"};
0519   std::vector<edm::InputTag> tags = {edm::InputTag("hltHITPixelTracksHB"), edm::InputTag("hltHITPixelTracksHE")};
0520   std::vector<double> cones = {35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 63.9, 70.0};
0521   edm::ParameterSetDescription desc;
0522   desc.addUntracked<std::vector<std::string>>("Triggers", triggers);
0523   desc.addUntracked<edm::InputTag>("pixCandTag", edm::InputTag(" "));
0524   desc.addUntracked<edm::InputTag>("l1CandTag", edm::InputTag("hltL1sV0SingleJet60"));
0525   desc.addUntracked<edm::InputTag>("l2CandTag", edm::InputTag("isolEcalPixelTrackProd"));
0526   desc.addUntracked<bool>("doL2L3", false);
0527   desc.addUntracked<bool>("doTimingTree", false);
0528   desc.addUntracked<bool>("doMipCutTree", false);
0529   desc.addUntracked<bool>("doTrkResTree", true);
0530   desc.addUntracked<bool>("doChgIsolTree", false);
0531   desc.addUntracked<bool>("doStudyIsol", false);
0532   desc.addUntracked<int>("verbosity", 0);
0533   desc.addUntracked<std::string>("processName", "HLT");
0534   desc.addUntracked<std::string>("trackQuality", "highPurity");
0535   desc.addUntracked<double>("minTrackPt", 10.0);
0536   desc.addUntracked<double>("maxDxyPV", 0.02);
0537   desc.addUntracked<double>("maxDzPV", 0.02);
0538   desc.addUntracked<double>("maxChi2", 5.0);
0539   desc.addUntracked<double>("maxDpOverP", 0.1);
0540   desc.addUntracked<int>("minOuterHit", 4);
0541   desc.addUntracked<int>("minLayerCrossed", 8);
0542   desc.addUntracked<int>("maxInMiss", 0);
0543   desc.addUntracked<int>("maxOutMiss", 0);
0544   desc.addUntracked<double>("isolationL1", 1.0);
0545   desc.addUntracked<double>("coneRadius", 34.98);
0546   desc.addUntracked<double>("coneRadiusMIP", 14.0);
0547   desc.addUntracked<double>("coneRadiusNeut1", 21.0);
0548   desc.addUntracked<double>("coneRadiusNeut2", 29.0);
0549   desc.addUntracked<double>("cutMIP", 1.0);
0550   desc.addUntracked<double>("chargeIsolation", 2.0);
0551   desc.addUntracked<double>("neutralIsolation", 2.0);
0552   desc.addUntracked<int>("minRun", 190456);
0553   desc.addUntracked<int>("maxRun", 203002);
0554   desc.addUntracked<std::vector<edm::InputTag>>("pixelTracksSources", tags);
0555   desc.addUntracked<std::vector<double>>("pixelIsolationConeSizeAtEC", cones);
0556   desc.addUntracked<double>("minPTrackValue", 0.0);
0557   desc.addUntracked<double>("vertexCutSeed", 101.0);
0558   desc.addUntracked<double>("vertexCutIsol", 101.0);
0559   desc.addUntracked<double>("tauUnbiasCone", 1.2);
0560   desc.addUntracked<double>("prelimCone", 1.0);
0561   desc.add<unsigned int>("stageL1Trigger", 1);
0562   descriptions.add("isoTrigDefault", desc);
0563 }
0564 
0565 void IsoTrig::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0566   if (verbosity_ % 10 > 0)
0567     edm::LogVerbatim("IsoTrack") << "Event starts====================================";
0568 
0569   int RunNo = iEvent.id().run();
0570 
0571   HLTConfigProvider const &hltConfig = hltPrescaleProvider_.hltConfigProvider();
0572 
0573   bField_ = &iSetup.getData(tok_magField_);
0574   geo_ = &iSetup.getData(tok_geom_);
0575   GlobalVector BField = bField_->inTesla(GlobalPoint(0, 0, 0));
0576   bfVal_ = BField.mag();
0577 
0578   trigger::TriggerEvent triggerEvent;
0579   edm::Handle<trigger::TriggerEvent> triggerEventHandle;
0580   iEvent.getByToken(tok_trigEvt_, triggerEventHandle);
0581   if (!triggerEventHandle.isValid()) {
0582     edm::LogWarning("IsoTrack") << "Error! Can't get the product hltTriggerSummaryAOD";
0583 
0584   } else {
0585     triggerEvent = *(triggerEventHandle.product());
0586   }
0587   const trigger::TriggerObjectCollection &TOC(triggerEvent.getObjects());
0588   /////////////////////////////TriggerResults
0589   edm::Handle<edm::TriggerResults> triggerResults;
0590   iEvent.getByToken(tok_trigRes_, triggerResults);
0591 
0592   edm::Handle<reco::TrackCollection> trkCollection;
0593   iEvent.getByToken(tok_genTrack_, trkCollection);
0594 
0595   iEvent.getByToken(tok_EB_, barrelRecHitsHandle_);
0596   iEvent.getByToken(tok_EE_, endcapRecHitsHandle_);
0597 
0598   iEvent.getByToken(tok_hbhe_, hbhe_);
0599 
0600   iEvent.getByToken(tok_recVtx_, recVtxs_);
0601   iEvent.getByToken(tok_bs_, beamSpotH_);
0602   if (!recVtxs_->empty() && !((*recVtxs_)[0].isFake())) {
0603     leadPV_ = math::XYZPoint((*recVtxs_)[0].x(), (*recVtxs_)[0].y(), (*recVtxs_)[0].z());
0604   } else if (beamSpotH_.isValid()) {
0605     leadPV_ = beamSpotH_->position();
0606   }
0607 
0608   if ((verbosity_ / 100) % 10 > 0) {
0609     edm::LogVerbatim("IsoTrack") << "Primary Vertex " << leadPV_;
0610     if (beamSpotH_.isValid())
0611       edm::LogVerbatim("IsoTrack") << "Beam Spot " << beamSpotH_->position();
0612   }
0613   pixelTrackRefsHE_.clear();
0614   pixelTrackRefsHB_.clear();
0615   for (unsigned int iPix = 0; iPix < pixelTracksSources_.size(); iPix++) {
0616     edm::Handle<reco::TrackCollection> iPixCol;
0617     iEvent.getByToken(tok_pixtks_[iPix], iPixCol);
0618     if (iPixCol.isValid()) {
0619       for (reco::TrackCollection::const_iterator pit = iPixCol->begin(); pit != iPixCol->end(); pit++) {
0620         if (iPix == 0)
0621           pixelTrackRefsHB_.push_back(reco::TrackRef(iPixCol, pit - iPixCol->begin()));
0622         pixelTrackRefsHE_.push_back(reco::TrackRef(iPixCol, pit - iPixCol->begin()));
0623       }
0624     }
0625   }
0626   if (doTiming_)
0627     getGoodTracks(iEvent, trkCollection);
0628 
0629   for (unsigned int ifilter = 0; ifilter < triggerEvent.sizeFilters(); ++ifilter) {
0630     std::string FilterNames[7] = {"hltL1sL1SingleJet68",
0631                                   "hltIsolPixelTrackL2FilterHE",
0632                                   "ecalIsolPixelTrackFilterHE",
0633                                   "hltIsolPixelTrackL3FilterHE",
0634                                   "hltIsolPixelTrackL2FilterHB",
0635                                   "ecalIsolPixelTrackFilterHB",
0636                                   "hltIsolPixelTrackL3FilterHB"};
0637     std::string label = triggerEvent.filterTag(ifilter).label();
0638     for (int i = 0; i < 7; i++) {
0639       if (label == FilterNames[i])
0640         h_Filters->Fill(i);
0641     }
0642   }
0643   edm::InputTag lumiProducer("LumiProducer", "", "RECO");
0644   edm::Handle<LumiDetails> Lumid;
0645   iEvent.getLuminosityBlock().getByToken(tok_lumi_, Lumid);
0646   float mybxlumi = -1;
0647   if (Lumid.isValid())
0648     mybxlumi = Lumid->lumiValue(LumiDetails::kOCC1, iEvent.bunchCrossing()) * 6.37;
0649   if (verbosity_ % 10 > 0)
0650     edm::LogVerbatim("IsoTrack") << "RunNo " << RunNo << " EvtNo " << iEvent.id().event() << " Lumi "
0651                                  << iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing() << " mybxlumi "
0652                                  << mybxlumi;
0653   if (!triggerResults.isValid()) {
0654     edm::LogWarning("IsoTrack") << "Error! Can't get the product triggerResults";
0655     //      std::shared_ptr<cms::Exception> const & error = triggerResults.whyFailed();
0656     //      edm::LogWarning(error->category()) << error->what();
0657   } else {
0658     std::vector<std::string> modules;
0659     h_nHLT->Fill(triggerResults->size());
0660     const edm::TriggerNames &triggerNames = iEvent.triggerNames(*triggerResults);
0661 
0662     const std::vector<std::string> &triggerNames_ = triggerNames.triggerNames();
0663     if (verbosity_ % 10 > 1)
0664       edm::LogVerbatim("IsoTrack") << "number of HLTs " << triggerNames_.size();
0665     double preL1(-1), preHLT(-1), prescale(-1);
0666     int hlt(-1);
0667     for (unsigned int i = 0; i < triggerResults->size(); i++) {
0668       unsigned int triggerindx = hltConfig.triggerIndex(triggerNames_[i]);
0669       const std::vector<std::string> &moduleLabels(hltConfig.moduleLabels(triggerindx));
0670 
0671       for (unsigned int in = 0; in < trigNames_.size(); ++in) {
0672         //    if (triggerNames_[i].find(trigNames_[in].c_str())!=std::string::npos || triggerNames_[i]==" ") {
0673         if (triggerNames_[i].find(trigNames_[in]) != std::string::npos) {
0674           if (verbosity_ % 10 > 0)
0675             edm::LogVerbatim("IsoTrack") << "trigger that i want " << triggerNames_[i] << " accept "
0676                                          << triggerResults->accept(i);
0677           hlt = triggerResults->accept(i);
0678           h_HLT->Fill(hlt);
0679           //        if (hlt>0 || triggerNames_[i]==" ") {
0680           if (hlt > 0) {
0681             edm::Handle<reco::IsolatedPixelTrackCandidateCollection> Pixcands;
0682             iEvent.getByToken(tok_pixtk_, Pixcands);
0683             edm::Handle<trigger::TriggerFilterObjectWithRefs> L1cands;
0684             iEvent.getByToken(tok_l1cand_, L1cands);
0685 
0686             auto const prescales = hltPrescaleProvider_.prescaleValues<double>(iEvent, iSetup, triggerNames_[i]);
0687             preL1 = prescales.first;
0688             preHLT = prescales.second;
0689             prescale = preL1 * preHLT;
0690             if (verbosity_ % 10 > 0)
0691               edm::LogVerbatim("IsoTrack")
0692                   << triggerNames_[i] << " accept " << hlt << " preL1 " << preL1 << " preHLT " << preHLT;
0693             for (int iv = 0; iv < 3; ++iv)
0694               vec_[iv].clear();
0695             if (trigList_.find(RunNo) != trigList_.end()) {
0696               trigList_[RunNo] += 1;
0697             } else {
0698               trigList_.insert({RunNo, 1});
0699               trigPreList_.insert({RunNo, prescales});
0700             }
0701             //loop over all trigger filters in event (i.e. filters passed)
0702             for (unsigned int ifilter = 0; ifilter < triggerEvent.sizeFilters(); ++ifilter) {
0703               std::vector<int> Keys;
0704               std::string label = triggerEvent.filterTag(ifilter).label();
0705               //loop over keys to objects passing this filter
0706               for (unsigned int imodule = 0; imodule < moduleLabels.size(); imodule++) {
0707                 if (label.find(moduleLabels[imodule]) != std::string::npos) {
0708                   if (verbosity_ % 10 > 0)
0709                     edm::LogVerbatim("IsoTrack") << "FILTERNAME " << label;
0710                   for (unsigned int ifiltrKey = 0; ifiltrKey < triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
0711                     Keys.push_back(triggerEvent.filterKeys(ifilter)[ifiltrKey]);
0712                     const trigger::TriggerObject &TO(TOC[Keys[ifiltrKey]]);
0713                     math::XYZTLorentzVector v4(TO.px(), TO.py(), TO.pz(), TO.energy());
0714                     if (label.find("L2Filter") != std::string::npos) {
0715                       vec_[1].push_back(v4);
0716                     } else if (label.find("L3Filter") != std::string::npos) {
0717                       vec_[2].push_back(v4);
0718                     } else {
0719                       vec_[0].push_back(v4);
0720                       h_L1ObjEnergy->Fill(TO.energy());
0721                     }
0722                     if (verbosity_ % 10 > 0)
0723                       edm::LogVerbatim("IsoTrack") << "key " << ifiltrKey << " : pt " << TO.pt() << " eta " << TO.eta()
0724                                                    << " phi " << TO.phi() << " mass " << TO.mass() << " Id " << TO.id();
0725                   }
0726                 }
0727               }
0728             }
0729             std::vector<reco::TrackCollection::const_iterator> goodTks;
0730             if (doL2L3_) {
0731               h_nL3Objs->Fill(vec_[2].size());
0732               studyTrigger(trkCollection, goodTks);
0733             } else {
0734               if (trkCollection.isValid()) {
0735                 reco::TrackCollection::const_iterator trkItr;
0736                 for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++)
0737                   goodTks.push_back(trkItr);
0738               }
0739             }
0740             // Now study isolation etc
0741             if (doStudyIsol_)
0742               studyIsolation(trkCollection, goodTks);
0743             if (doTrkResTree_)
0744               StudyTrkEbyP(trkCollection);
0745 
0746             std::pair<double, double> etaphi = etaPhiTrigger();
0747             edm::Handle<reco::IsolatedPixelTrackCandidateCollection> L2cands;
0748             iEvent.getByToken(tok_l2cand_, L2cands);
0749             if (!L2cands.isValid()) {
0750               if (verbosity_ % 10 > 0)
0751                 edm::LogVerbatim("IsoTrack") << " trigCand is not valid ";
0752             } else {
0753               if (doMipCutTree_)
0754                 studyMipCut(trkCollection, L2cands);
0755             }
0756             if (!pixelTracksSources_.empty())
0757               if (doChgIsolTree_ && !pixelTrackRefsHE_.empty())
0758                 chgIsolation(etaphi.first, etaphi.second, trkCollection, iEvent);
0759           }
0760           break;
0761         }
0762       }
0763     }
0764     h_PreL1->Fill(preL1);
0765     h_PreHLT->Fill(preHLT);
0766     h_Pre->Fill(prescale);
0767     h_PreL1wt->Fill(preL1, mybxlumi);
0768     h_PreHLTwt->Fill(preHLT, mybxlumi);
0769 
0770     // check if trigger names in (new) config
0771     //      edm::LogVerbatim("IsoTrack") << "changed " << changed_;
0772     if (changed_) {
0773       changed_ = false;
0774       if ((verbosity_ / 10) % 10 > 1) {
0775         edm::LogVerbatim("IsoTrack") << "New trigger menu found !!!";
0776         const unsigned int n(hltConfig.size());
0777         for (unsigned itrig = 0; itrig < triggerNames_.size(); itrig++) {
0778           unsigned int triggerindx = hltConfig.triggerIndex(triggerNames_[itrig]);
0779           if (triggerindx >= n)
0780             edm::LogVerbatim("IsoTrack") << triggerNames_[itrig] << " " << triggerindx << " does not exist in"
0781                                          << " the current menu";
0782           else
0783             edm::LogVerbatim("IsoTrack") << triggerNames_[itrig] << " " << triggerindx << " exists";
0784         }
0785       }
0786     }
0787   }
0788   if (doTiming_)
0789     studyTiming(iEvent);
0790 }
0791 
0792 void IsoTrig::clearChgIsolnTreeVectors() {
0793   t_PixcandP->clear();
0794   t_PixcandPt->clear();
0795   t_PixcandEta->clear();
0796   t_PixcandPhi->clear();
0797   for (unsigned int i = 0; i < t_PixcandMaxP->size(); i++)
0798     t_PixcandMaxP[i].clear();
0799   t_PixcandMaxP->clear();
0800   t_PixTrkcandP->clear();
0801   t_PixTrkcandPt->clear();
0802   t_PixTrkcandEta->clear();
0803   t_PixTrkcandPhi->clear();
0804   t_PixTrkcandMaxP->clear();
0805   t_PixTrkcandselTk->clear();
0806 }
0807 
0808 void IsoTrig::clearMipCutTreeVectors() {
0809   t_NFcandP->clear();
0810   t_NFcandPt->clear();
0811   t_NFcandEta->clear();
0812   t_NFcandPhi->clear();
0813   t_NFcandEmip->clear();
0814   t_NFTrkcandP->clear();
0815   t_NFTrkcandPt->clear();
0816   t_NFTrkcandEta->clear();
0817   t_NFTrkcandPhi->clear();
0818   t_NFTrkcandEmip->clear();
0819   t_NFTrkMinDR->clear();
0820   t_NFTrkMinDP1->clear();
0821   t_NFTrkselTkFlag->clear();
0822   t_NFTrkqltyFlag->clear();
0823   t_NFTrkMissFlag->clear();
0824   t_NFTrkPVFlag->clear();
0825   t_NFTrkPropFlag->clear();
0826   t_NFTrkChgIsoFlag->clear();
0827   t_NFTrkNeuIsoFlag->clear();
0828   t_NFTrkMipFlag->clear();
0829   t_ECone->clear();
0830 }
0831 
0832 void IsoTrig::pushChgIsolnTreeVecs(math::XYZTLorentzVector &Pixcand,
0833                                    math::XYZTLorentzVector &Trkcand,
0834                                    std::vector<double> &PixMaxP,
0835                                    double &TrkMaxP,
0836                                    bool &selTk) {
0837   t_PixcandP->push_back(Pixcand.r());
0838   t_PixcandPt->push_back(Pixcand.pt());
0839   t_PixcandEta->push_back(Pixcand.eta());
0840   t_PixcandPhi->push_back(Pixcand.phi());
0841   t_PixcandMaxP->push_back(PixMaxP);
0842   t_PixTrkcandP->push_back(Trkcand.r());
0843   t_PixTrkcandPt->push_back(Trkcand.pt());
0844   t_PixTrkcandEta->push_back(Trkcand.eta());
0845   t_PixTrkcandPhi->push_back(Trkcand.phi());
0846   t_PixTrkcandMaxP->push_back(TrkMaxP);
0847   t_PixTrkcandselTk->push_back(selTk);
0848 }
0849 
0850 void IsoTrig::pushMipCutTreeVecs(math::XYZTLorentzVector &NFcand,
0851                                  math::XYZTLorentzVector &Trkcand,
0852                                  double &EmipNFcand,
0853                                  double &EmipTrkcand,
0854                                  double &mindR,
0855                                  double &mindP1,
0856                                  std::vector<bool> &Flags,
0857                                  double hCone) {
0858   t_NFcandP->push_back(NFcand.r());
0859   t_NFcandPt->push_back(NFcand.pt());
0860   t_NFcandEta->push_back(NFcand.eta());
0861   t_NFcandPhi->push_back(NFcand.phi());
0862   t_NFcandEmip->push_back(EmipNFcand);
0863   t_NFTrkcandP->push_back(Trkcand.r());
0864   t_NFTrkcandPt->push_back(Trkcand.pt());
0865   t_NFTrkcandEta->push_back(Trkcand.eta());
0866   t_NFTrkcandPhi->push_back(Trkcand.phi());
0867   t_NFTrkcandEmip->push_back(EmipTrkcand);
0868   t_NFTrkMinDR->push_back(mindR);
0869   t_NFTrkMinDP1->push_back(mindP1);
0870   t_NFTrkselTkFlag->push_back(Flags[0]);
0871   t_NFTrkqltyFlag->push_back(Flags[1]);
0872   t_NFTrkMissFlag->push_back(Flags[2]);
0873   t_NFTrkPVFlag->push_back(Flags[3]);
0874   t_NFTrkPropFlag->push_back(Flags[4]);
0875   t_NFTrkChgIsoFlag->push_back(Flags[5]);
0876   t_NFTrkNeuIsoFlag->push_back(Flags[6]);
0877   t_NFTrkMipFlag->push_back(Flags[7]);
0878   t_ECone->push_back(hCone);
0879 }
0880 
0881 void IsoTrig::beginJob() {
0882   char hname[100], htit[100];
0883   std::string levels[20] = {"L1",        "L2",          "L3",        "Reco",      "RecoMatch",     "RecoNoMatch",
0884                             "L2Match",   "L2NoMatch",   "L3Match",   "L3NoMatch", "HLTTrk",        "HLTGoodTrk",
0885                             "HLTIsoTrk", "HLTMip",      "HLTSelect", "nonHLTTrk", "nonHLTGoodTrk", "nonHLTIsoTrk",
0886                             "nonHLTMip", "nonHLTSelect"};
0887   if (doTiming_) {
0888     TimingTree_ = fs_->make<TTree>("TimingTree", "TimingTree");
0889     t_timeL2Prod = new std::vector<double>();
0890     t_nPixCand = new std::vector<int>();
0891     t_nPixSeed = new std::vector<int>();
0892     t_nGoodTk = new std::vector<int>();
0893 
0894     TimingTree_->Branch("t_timeL2Prod", "std::vector<double>", &t_timeL2Prod);
0895     TimingTree_->Branch("t_nPixCand", "std::vector<int>", &t_nPixCand);
0896     TimingTree_->Branch("t_nPixSeed", "std::vector<int>", &t_nPixSeed);
0897     TimingTree_->Branch("t_nGoodTk", "std::vector<int>", &t_nGoodTk);
0898   }
0899   if (doTrkResTree_) {
0900     TrkResTree_ = fs_->make<TTree>("TrkRestree", "TrkResTree");
0901     t_TrkhCone = new std::vector<double>();
0902     t_TrkP = new std::vector<double>();
0903     t_TrkselTkFlag = new std::vector<bool>();
0904     t_TrkqltyFlag = new std::vector<bool>();
0905     t_TrkMissFlag = new std::vector<bool>();
0906     t_TrkPVFlag = new std::vector<bool>();
0907     t_TrkNuIsolFlag = new std::vector<bool>();
0908 
0909     TrkResTree_->Branch("t_TrkhCone", "std::vector<double>", &t_TrkhCone);
0910     TrkResTree_->Branch("t_TrkP", "std::vector<double>", &t_TrkP);
0911     TrkResTree_->Branch("t_TrkselTkFlag", "std::vector<bool>", &t_TrkselTkFlag);
0912     TrkResTree_->Branch("t_TrkqltyFlag", "std::vector<bool>", &t_TrkqltyFlag);
0913     TrkResTree_->Branch("t_TrkMissFlag", "std::vector<bool>", &t_TrkMissFlag);
0914     TrkResTree_->Branch("t_TrkPVFlag", "std::vector<bool>", &t_TrkPVFlag);
0915     TrkResTree_->Branch("t_TrkNuIsolFlag", "std::vector<bool>", &t_TrkNuIsolFlag);
0916   }
0917   if (doChgIsolTree_) {
0918     ChgIsolnTree_ = fs_->make<TTree>("ChgIsolnTree", "ChgIsolnTree");
0919     t_PixcandP = new std::vector<double>();
0920     t_PixcandPt = new std::vector<double>();
0921     t_PixcandEta = new std::vector<double>();
0922     t_PixcandPhi = new std::vector<double>();
0923     t_PixcandMaxP = new std::vector<std::vector<double>>();
0924     t_PixTrkcandP = new std::vector<double>();
0925     t_PixTrkcandPt = new std::vector<double>();
0926     t_PixTrkcandEta = new std::vector<double>();
0927     t_PixTrkcandPhi = new std::vector<double>();
0928     t_PixTrkcandMaxP = new std::vector<double>();
0929     t_PixTrkcandselTk = new std::vector<bool>();
0930 
0931     ChgIsolnTree_->Branch("t_PixcandP", "std::vector<double>", &t_PixcandP);
0932     ChgIsolnTree_->Branch("t_PixcandPt", "std::vector<double>", &t_PixcandPt);
0933     ChgIsolnTree_->Branch("t_PixcandEta", "std::vector<double>", &t_PixcandEta);
0934     ChgIsolnTree_->Branch("t_PixcandPhi", "std::vector<double>", &t_PixcandPhi);
0935     ChgIsolnTree_->Branch("t_PixcandMaxP", "std::vector<std::vector<double> >", &t_PixcandMaxP);
0936     ChgIsolnTree_->Branch("t_PixTrkcandP", "std::vector<double>", &t_PixTrkcandP);
0937     ChgIsolnTree_->Branch("t_PixTrkcandPt", "std::vector<double>", &t_PixTrkcandPt);
0938     ChgIsolnTree_->Branch("t_PixTrkcandEta", "std::vector<double>", &t_PixTrkcandEta);
0939     ChgIsolnTree_->Branch("t_PixTrkcandPhi", "std::vector<double>", &t_PixTrkcandPhi);
0940     ChgIsolnTree_->Branch("t_PixTrkcandMaxP", "std::vector<double>", &t_PixTrkcandMaxP);
0941     ChgIsolnTree_->Branch("t_PixTrkcandselTk", "std::vector<bool>", &t_PixTrkcandselTk);
0942   }
0943   if (doMipCutTree_) {
0944     MipCutTree_ = fs_->make<TTree>("MipCutTree", "MipCutTree");
0945     t_NFcandP = new std::vector<double>();
0946     t_NFcandPt = new std::vector<double>();
0947     t_NFcandEta = new std::vector<double>();
0948     t_NFcandPhi = new std::vector<double>();
0949     t_NFcandEmip = new std::vector<double>();
0950     t_NFTrkcandP = new std::vector<double>();
0951     t_NFTrkcandPt = new std::vector<double>();
0952     t_NFTrkcandEta = new std::vector<double>();
0953     t_NFTrkcandPhi = new std::vector<double>();
0954     t_NFTrkcandEmip = new std::vector<double>();
0955     t_NFTrkMinDR = new std::vector<double>();
0956     t_NFTrkMinDP1 = new std::vector<double>();
0957     t_NFTrkselTkFlag = new std::vector<bool>();
0958     t_NFTrkqltyFlag = new std::vector<bool>();
0959     t_NFTrkMissFlag = new std::vector<bool>();
0960     t_NFTrkPVFlag = new std::vector<bool>();
0961     t_NFTrkPropFlag = new std::vector<bool>();
0962     t_NFTrkChgIsoFlag = new std::vector<bool>();
0963     t_NFTrkNeuIsoFlag = new std::vector<bool>();
0964     t_NFTrkMipFlag = new std::vector<bool>();
0965     t_ECone = new std::vector<double>();
0966 
0967     MipCutTree_->Branch("t_NFcandP", "std::vector<double>", &t_NFcandP);
0968     MipCutTree_->Branch("t_NFcandPt", "std::vector<double>", &t_NFcandPt);
0969     MipCutTree_->Branch("t_NFcandEta", "std::vector<double>", &t_NFcandEta);
0970     MipCutTree_->Branch("t_NFcandPhi", "std::vector<double>", &t_NFcandPhi);
0971     MipCutTree_->Branch("t_NFcandEmip", "std::vector<double>", &t_NFcandEmip);
0972     MipCutTree_->Branch("t_NFTrkcandP", "std::vector<double>", &t_NFTrkcandP);
0973     MipCutTree_->Branch("t_NFTrkcandPt", "std::vector<double>", &t_NFTrkcandPt);
0974     MipCutTree_->Branch("t_NFTrkcandEta", "std::vector<double>", &t_NFTrkcandEta);
0975     MipCutTree_->Branch("t_NFTrkcandPhi", "std::vector<double>", &t_NFTrkcandPhi);
0976     MipCutTree_->Branch("t_NFTrkcandEmip", "std::vector<double>", &t_NFTrkcandEmip);
0977     MipCutTree_->Branch("t_NFTrkMinDR", "std::vector<double>", &t_NFTrkMinDR);
0978     MipCutTree_->Branch("t_NFTrkMinDP1", "std::vector<double>", &t_NFTrkMinDP1);
0979     MipCutTree_->Branch("t_NFTrkselTkFlag", "std::vector<bool>", &t_NFTrkselTkFlag);
0980     MipCutTree_->Branch("t_NFTrkqltyFlag", "std::vector<bool>", &t_NFTrkqltyFlag);
0981     MipCutTree_->Branch("t_NFTrkMissFlag", "std::vector<bool>", &t_NFTrkMissFlag);
0982     MipCutTree_->Branch("t_NFTrkPVFlag", "std::vector<bool>", &t_NFTrkPVFlag);
0983     MipCutTree_->Branch("t_NFTrkPropFlag", "std::vector<bool>", &t_NFTrkPropFlag);
0984     MipCutTree_->Branch("t_NFTrkChgIsoFlag", "std::vector<bool>", &t_NFTrkChgIsoFlag);
0985     MipCutTree_->Branch("t_NFTrkNeuIsoFlag", "std::vector<bool>", &t_NFTrkNeuIsoFlag);
0986     MipCutTree_->Branch("t_NFTrkMipFlag", "std::vector<bool>", &t_NFTrkMipFlag);
0987     MipCutTree_->Branch("t_ECone", "std::vector<double>", &t_ECone);
0988   }
0989   h_Filters = fs_->make<TH1I>("h_Filters", "Filter Accepts", 10, 0, 10);
0990   std::string FilterNames[7] = {"hltL1sL1SingleJet68",
0991                                 "hltIsolPixelTrackL2FilterHE",
0992                                 "ecalIsolPixelTrackFilterHE",
0993                                 "hltIsolPixelTrackL3FilterHE",
0994                                 "hltIsolPixelTrackL2FilterHB",
0995                                 "ecalIsolPixelTrackFilterHB",
0996                                 "hltIsolPixelTrackL3FilterHB"};
0997   for (int i = 0; i < 7; i++) {
0998     h_Filters->GetXaxis()->SetBinLabel(i + 1, FilterNames[i].c_str());
0999   }
1000 
1001   h_nHLT = fs_->make<TH1I>("h_nHLT", "Size of trigger Names", 1000, 1, 1000);
1002   h_HLT = fs_->make<TH1I>("h_HLT", "HLT accept", 3, -1, 2);
1003   h_PreL1 = fs_->make<TH1I>("h_PreL1", "L1 Prescale", 500, 0, 500);
1004   h_PreHLT = fs_->make<TH1I>("h_PreHLT", "HLT Prescale", 50, 0, 50);
1005   h_Pre = fs_->make<TH1I>("h_Pre", "Prescale", 3000, 0, 3000);
1006 
1007   h_PreL1wt = fs_->make<TH1D>("h_PreL1wt", "Weighted L1 Prescale", 500, 0, 500);
1008   h_PreHLTwt = fs_->make<TH1D>("h_PreHLTwt", "Weighted HLT Prescale", 500, 0, 100);
1009   h_L1ObjEnergy = fs_->make<TH1D>("h_L1ObjEnergy", "Energy of L1Object", 500, 0.0, 500.0);
1010 
1011   h_EnIn = fs_->make<TH1D>("h_EnInEcal", "EnergyIn Ecal", 200, 0.0, 20.0);
1012   h_EnOut = fs_->make<TH1D>("h_EnOutEcal", "EnergyOut Ecal", 200, 0.0, 20.0);
1013   h_MipEnMatch =
1014       fs_->make<TH2D>("h_MipEnMatch", "MipEn: HLT level vs Reco Level (Matched)", 200, 0.0, 20.0, 200, 0.0, 20.0);
1015   h_MipEnNoMatch = fs_->make<TH2D>(
1016       "h_MipEnNoMatch", "MipEn: HLT level vs Reco Level (No Match Found)", 200, 0.0, 20.0, 200, 0.0, 20.0);
1017 
1018   if (doL2L3_) {
1019     h_nL3Objs = fs_->make<TH1I>("h_nL3Objs", "Number of L3 objects", 10, 0, 10);
1020 
1021     std::string pairs[9] = {"L2L3",
1022                             "L2L3Match",
1023                             "L2L3NoMatch",
1024                             "L3Reco",
1025                             "L3RecoMatch",
1026                             "L3RecoNoMatch",
1027                             "NewFilterReco",
1028                             "NewFilterRecoMatch",
1029                             "NewFilterRecoNoMatch"};
1030     for (int ipair = 0; ipair < 9; ipair++) {
1031       sprintf(hname, "h_dEta%s", pairs[ipair].c_str());
1032       sprintf(htit, "#Delta#eta for %s", pairs[ipair].c_str());
1033       h_dEta[ipair] = fs_->make<TH1D>(hname, htit, 200, -10.0, 10.0);
1034       h_dEta[ipair]->GetXaxis()->SetTitle("d#eta");
1035 
1036       sprintf(hname, "h_dPhi%s", pairs[ipair].c_str());
1037       sprintf(htit, "#Delta#phi for %s", pairs[ipair].c_str());
1038       h_dPhi[ipair] = fs_->make<TH1D>(hname, htit, 140, -7.0, 7.0);
1039       h_dPhi[ipair]->GetXaxis()->SetTitle("d#phi");
1040 
1041       sprintf(hname, "h_dPt%s", pairs[ipair].c_str());
1042       sprintf(htit, "#Delta dp_{T} for %s objects", pairs[ipair].c_str());
1043       h_dPt[ipair] = fs_->make<TH1D>(hname, htit, 400, -200.0, 200.0);
1044       h_dPt[ipair]->GetXaxis()->SetTitle("dp_{T} (GeV)");
1045 
1046       sprintf(hname, "h_dP%s", pairs[ipair].c_str());
1047       sprintf(htit, "#Delta p for %s objects", pairs[ipair].c_str());
1048       h_dP[ipair] = fs_->make<TH1D>(hname, htit, 400, -200.0, 200.0);
1049       h_dP[ipair]->GetXaxis()->SetTitle("dP (GeV)");
1050 
1051       sprintf(hname, "h_dinvPt%s", pairs[ipair].c_str());
1052       sprintf(htit, "#Delta (1/p_{T}) for %s objects", pairs[ipair].c_str());
1053       h_dinvPt[ipair] = fs_->make<TH1D>(hname, htit, 500, -0.4, 0.1);
1054       h_dinvPt[ipair]->GetXaxis()->SetTitle("d(1/p_{T})");
1055       sprintf(hname, "h_mindR%s", pairs[ipair].c_str());
1056       sprintf(htit, "min(#Delta R) for %s objects", pairs[ipair].c_str());
1057       h_mindR[ipair] = fs_->make<TH1D>(hname, htit, 500, 0.0, 1.0);
1058       h_mindR[ipair]->GetXaxis()->SetTitle("dR");
1059     }
1060 
1061     for (int lvl = 0; lvl < 2; lvl++) {
1062       sprintf(hname, "h_dEtaL1%s", levels[lvl + 1].c_str());
1063       sprintf(htit, "#Delta#eta for L1 and %s objects", levels[lvl + 1].c_str());
1064       h_dEtaL1[lvl] = fs_->make<TH1D>(hname, htit, 400, -10.0, 10.0);
1065 
1066       sprintf(hname, "h_dPhiL1%s", levels[lvl + 1].c_str());
1067       sprintf(htit, "#Delta#phi for L1 and %s objects", levels[lvl + 1].c_str());
1068       h_dPhiL1[lvl] = fs_->make<TH1D>(hname, htit, 280, -7.0, 7.0);
1069 
1070       sprintf(hname, "h_dRL1%s", levels[lvl + 1].c_str());
1071       sprintf(htit, "#Delta R for L1 and %s objects", levels[lvl + 1].c_str());
1072       h_dRL1[lvl] = fs_->make<TH1D>(hname, htit, 100, 0.0, 10.0);
1073     }
1074   }
1075 
1076   int levmin = (doL2L3_ ? 0 : 10);
1077   for (int ilevel = levmin; ilevel < 20; ilevel++) {
1078     sprintf(hname, "h_p%s", levels[ilevel].c_str());
1079     sprintf(htit, "p for %s objects", levels[ilevel].c_str());
1080     h_p[ilevel] = fs_->make<TH1D>(hname, htit, 100, 0.0, 500.0);
1081     h_p[ilevel]->GetXaxis()->SetTitle("p (GeV)");
1082 
1083     sprintf(hname, "h_pt%s", levels[ilevel].c_str());
1084     sprintf(htit, "p_{T} for %s objects", levels[ilevel].c_str());
1085     h_pt[ilevel] = fs_->make<TH1D>(hname, htit, 100, 0.0, 500.0);
1086     h_pt[ilevel]->GetXaxis()->SetTitle("p_{T} (GeV)");
1087 
1088     sprintf(hname, "h_eta%s", levels[ilevel].c_str());
1089     sprintf(htit, "#eta for %s objects", levels[ilevel].c_str());
1090     h_eta[ilevel] = fs_->make<TH1D>(hname, htit, 100, -5.0, 5.0);
1091     h_eta[ilevel]->GetXaxis()->SetTitle("#eta");
1092 
1093     sprintf(hname, "h_phi%s", levels[ilevel].c_str());
1094     sprintf(htit, "#phi for %s objects", levels[ilevel].c_str());
1095     h_phi[ilevel] = fs_->make<TH1D>(hname, htit, 70, -3.5, 3.50);
1096     h_phi[ilevel]->GetXaxis()->SetTitle("#phi");
1097   }
1098 
1099   std::string cuts[2] = {"HLTMatched", "HLTNotMatched"};
1100   std::string cuts2[2] = {"All", "Away from L1"};
1101   for (int icut = 0; icut < 2; icut++) {
1102     sprintf(hname, "h_eMip%s", cuts[icut].c_str());
1103     sprintf(htit, "eMip for %s tracks", cuts[icut].c_str());
1104     h_eMip[icut] = fs_->make<TH1D>(hname, htit, 200, 0.0, 10.0);
1105     h_eMip[icut]->GetXaxis()->SetTitle("E_{Mip} (GeV)");
1106 
1107     sprintf(hname, "h_eMaxNearP%s", cuts[icut].c_str());
1108     sprintf(htit, "eMaxNearP for %s tracks", cuts[icut].c_str());
1109     h_eMaxNearP[icut] = fs_->make<TH1D>(hname, htit, 240, -2.0, 10.0);
1110     h_eMaxNearP[icut]->GetXaxis()->SetTitle("E_{MaxNearP} (GeV)");
1111 
1112     sprintf(hname, "h_eNeutIso%s", cuts[icut].c_str());
1113     sprintf(htit, "eNeutIso for %s ", cuts[icut].c_str());
1114     h_eNeutIso[icut] = fs_->make<TH1D>(hname, htit, 200, 0.0, 10.0);
1115     h_eNeutIso[icut]->GetXaxis()->SetTitle("E_{NeutIso} (GeV)");
1116 
1117     for (int kcut = 0; kcut < 2; ++kcut) {
1118       for (int lim = 0; lim < 5; ++lim) {
1119         sprintf(hname, "h_etaCalibTracks%sCut%dLim%d", cuts[icut].c_str(), kcut, lim);
1120         sprintf(htit,
1121                 "#eta for %s isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)",
1122                 cuts[icut].c_str(),
1123                 pLimits_[lim],
1124                 pLimits_[lim + 1],
1125                 cuts2[kcut].c_str());
1126         h_etaCalibTracks[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 60, -30.0, 30.0);
1127         h_etaCalibTracks[lim][icut][kcut]->GetXaxis()->SetTitle("i#eta");
1128 
1129         sprintf(hname, "h_etaMipTracks%sCut%dLim%d", cuts[icut].c_str(), kcut, lim);
1130         sprintf(htit,
1131                 "#eta for %s charge isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)",
1132                 cuts[icut].c_str(),
1133                 pLimits_[lim],
1134                 pLimits_[lim + 1],
1135                 cuts2[kcut].c_str());
1136         h_etaMipTracks[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 60, -30.0, 30.0);
1137         h_etaMipTracks[lim][icut][kcut]->GetXaxis()->SetTitle("i#eta");
1138       }
1139     }
1140   }
1141 
1142   std::string ecut1[3] = {"all", "HLTMatched", "HLTNotMatched"};
1143   std::string ecut2[2] = {"without", "with"};
1144   int etac[48] = {-1,  -2,  -3,  -4,  -5,  -6,  -7,  -8,  -9, -10, -11, -12, -13, -14, -15, -16,
1145                   -17, -18, -19, -20, -21, -22, -23, -24, 1,  2,   3,   4,   5,   6,   7,   8,
1146                   9,   10,  11,  12,  13,  14,  15,  16,  17, 18,  19,  20,  21,  22,  23,  24};
1147   for (int icut = 0; icut < 6; icut++) {
1148     //    int i1 = (icut>3 ? 1 : 0);
1149     int i1 = (icut > 2 ? 1 : 0);
1150     int i2 = icut - i1 * 3;
1151     for (int kcut = 0; kcut < 48; kcut++) {
1152       for (int lim = 0; lim < 5; ++lim) {
1153         sprintf(hname, "h_eta%dEnHcal%s%s%d", etac[kcut], ecut1[i2].c_str(), ecut2[i1].c_str(), lim);
1154         sprintf(htit,
1155                 "HCAL energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation",
1156                 etac[kcut],
1157                 ecut1[i2].c_str(),
1158                 pLimits_[lim],
1159                 pLimits_[lim + 1],
1160                 ecut2[i1].c_str());
1161         h_eHcal[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 750, 0.0, 150.0);
1162         h_eHcal[lim][icut][kcut]->GetXaxis()->SetTitle("Energy (GeV)");
1163         sprintf(hname, "h_eta%dEnCalo%s%s%d", etac[kcut], ecut1[i2].c_str(), ecut2[i1].c_str(), lim);
1164         sprintf(htit,
1165                 "Calorimter energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation",
1166                 etac[kcut],
1167                 ecut1[i2].c_str(),
1168                 pLimits_[lim],
1169                 pLimits_[lim + 1],
1170                 ecut2[i1].c_str());
1171         h_eCalo[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 750, 0.0, 150.0);
1172         h_eCalo[lim][icut][kcut]->GetXaxis()->SetTitle("Energy (GeV)");
1173       }
1174     }
1175   }
1176 }
1177 
1178 // ------------ method called once each job just after ending the event loop  ------------
1179 void IsoTrig::endJob() {
1180   edm::LogWarning("IsoTrack") << trigNames_.size() << "Triggers were run. RunNo vs HLT accepts for";
1181   for (unsigned int i = 0; i < trigNames_.size(); ++i)
1182     edm::LogWarning("IsoTrack") << "[" << i << "]: " << trigNames_[i];
1183   unsigned int n = maxRunNo_ - minRunNo_ + 1;
1184   g_Pre = fs_->make<TH1D>("h_PrevsRN", "PreScale Vs Run Number", n, minRunNo_, maxRunNo_);
1185   g_PreL1 = fs_->make<TH1D>("h_PreL1vsRN", "L1 PreScale Vs Run Number", n, minRunNo_, maxRunNo_);
1186   g_PreHLT = fs_->make<TH1D>("h_PreHLTvsRN", "HLT PreScale Vs Run Number", n, minRunNo_, maxRunNo_);
1187   g_Accepts = fs_->make<TH1I>("h_HLTAcceptsvsRN", "HLT Accepts Vs Run Number", n, minRunNo_, maxRunNo_);
1188 
1189   for (auto const &[runNum, nAccept] : trigList_) {
1190     auto const &triggerPrescales = trigPreList_[runNum];
1191     auto const preL1 = triggerPrescales.first;
1192     auto const preHLT = triggerPrescales.second;
1193     edm::LogVerbatim("IsoTrack") << runNum << " " << nAccept << " " << preL1 << " " << preHLT;
1194     g_Accepts->Fill(runNum, nAccept);
1195     g_PreL1->Fill(runNum, preL1);
1196     g_PreHLT->Fill(runNum, preHLT);
1197     g_Pre->Fill(runNum, preL1 * preHLT);
1198   }
1199 }
1200 
1201 void IsoTrig::beginRun(edm::Run const &iRun, edm::EventSetup const &iSetup) {
1202   edm::LogWarning("IsoTrack") << "Run " << iRun.run() << " hltconfig.init "
1203                               << hltPrescaleProvider_.init(iRun, iSetup, processName_, changed_);
1204 }
1205 
1206 void IsoTrig::StudyTrkEbyP(edm::Handle<reco::TrackCollection> &trkCollection) {
1207   t_TrkselTkFlag->clear();
1208   t_TrkqltyFlag->clear();
1209   t_TrkMissFlag->clear();
1210   t_TrkPVFlag->clear();
1211   t_TrkNuIsolFlag->clear();
1212   t_TrkhCone->clear();
1213   t_TrkP->clear();
1214 
1215   if (!trkCollection.isValid()) {
1216     edm::LogVerbatim("IsoTrack") << "trkCollection.isValid is false";
1217   } else {
1218     std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1219     std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1220     spr::propagateCALO(
1221         trkCollection, geo_, bField_, theTrackQuality_, trkCaloDirections1, ((verbosity_ / 100) % 10 > 2));
1222     unsigned int nTracks = 0;
1223     int nRH_eMipDR = 0, nNearTRKs = 0;
1224     std::vector<bool> selFlags;
1225     for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++, nTracks++) {
1226       double conehmaxNearP = 0, hCone = 0, eMipDR = 0.0;
1227       const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1228       if (verbosity_ % 10 > 0)
1229         edm::LogVerbatim("IsoTrack") << "track no. " << nTracks << " p(): " << pTrack->p();
1230       if (pTrack->p() > 20) {
1231         math::XYZTLorentzVector v2(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
1232         eMipDR = spr::eCone_ecal(geo_,
1233                                  barrelRecHitsHandle_,
1234                                  endcapRecHitsHandle_,
1235                                  trkDetItr->pointHCAL,
1236                                  trkDetItr->pointECAL,
1237                                  a_mipR_,
1238                                  trkDetItr->directionECAL,
1239                                  nRH_eMipDR);
1240         bool selectTk = spr::goodTrack(pTrack, leadPV_, selectionParameters_, ((verbosity_ / 100) % 10 > 1));
1241         spr::trackSelectionParameters oneCutParameters = selectionParameters_;
1242         oneCutParameters.maxDxyPV = 10;
1243         oneCutParameters.maxDzPV = 100;
1244         oneCutParameters.maxInMiss = 2;
1245         oneCutParameters.maxOutMiss = 2;
1246         bool qltyFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1247         oneCutParameters = selectionParameters_;
1248         oneCutParameters.maxDxyPV = 10;
1249         oneCutParameters.maxDzPV = 100;
1250         bool qltyMissFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1251         oneCutParameters = selectionParameters_;
1252         oneCutParameters.maxInMiss = 2;
1253         oneCutParameters.maxOutMiss = 2;
1254         bool qltyPVFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1255         if ((verbosity_ / 1000) % 10 > 1)
1256           edm::LogVerbatim("IsoTrack") << "sel " << selectTk << "ntracks " << nTracks << " a_charIsoR " << a_charIsoR_
1257                                        << " nNearTRKs " << nNearTRKs;
1258         conehmaxNearP = spr::chargeIsolationCone(
1259             nTracks, trkCaloDirections1, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 1));
1260         if ((verbosity_ / 1000) % 10 > 1)
1261           edm::LogVerbatim("IsoTrack") << "coneh " << conehmaxNearP << "ok " << trkDetItr->okECAL << " "
1262                                        << trkDetItr->okHCAL;
1263         double e1 = spr::eCone_ecal(geo_,
1264                                     barrelRecHitsHandle_,
1265                                     endcapRecHitsHandle_,
1266                                     trkDetItr->pointHCAL,
1267                                     trkDetItr->pointECAL,
1268                                     a_neutR1_,
1269                                     trkDetItr->directionECAL,
1270                                     nRH_eMipDR);
1271         double e2 = spr::eCone_ecal(geo_,
1272                                     barrelRecHitsHandle_,
1273                                     endcapRecHitsHandle_,
1274                                     trkDetItr->pointHCAL,
1275                                     trkDetItr->pointECAL,
1276                                     a_neutR2_,
1277                                     trkDetItr->directionECAL,
1278                                     nRH_eMipDR);
1279         double e_inCone = e2 - e1;
1280         bool chgIsolFlag = (conehmaxNearP < cutCharge_);
1281         bool mipFlag = (eMipDR < cutMip_);
1282         bool neuIsolFlag = (e_inCone < cutNeutral_);
1283         bool trkpropFlag = ((trkDetItr->okECAL) && (trkDetItr->okHCAL));
1284         selFlags.clear();
1285         selFlags.push_back(selectTk);
1286         selFlags.push_back(qltyFlag);
1287         selFlags.push_back(qltyMissFlag);
1288         selFlags.push_back(qltyPVFlag);
1289         if (verbosity_ % 10 > 0)
1290           edm::LogVerbatim("IsoTrack") << "emip: " << eMipDR << "<" << cutMip_ << "(" << mipFlag << ")"
1291                                        << " ; ok: " << trkDetItr->okECAL << "/" << trkDetItr->okHCAL
1292                                        << " ; chgiso: " << conehmaxNearP << "<" << cutCharge_ << "(" << chgIsolFlag
1293                                        << ")";
1294 
1295         if (chgIsolFlag && mipFlag && trkpropFlag) {
1296           double distFromHotCell = -99.0;
1297           int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1298           GlobalPoint gposHotCell(0., 0., 0.);
1299           std::vector<DetId> coneRecHitDetIds;
1300           hCone = spr::eCone_hcal(geo_,
1301                                   hbhe_,
1302                                   trkDetItr->pointHCAL,
1303                                   trkDetItr->pointECAL,
1304                                   a_coneR_,
1305                                   trkDetItr->directionHCAL,
1306                                   nRecHitsCone,
1307                                   coneRecHitDetIds,
1308                                   distFromHotCell,
1309                                   ietaHotCell,
1310                                   iphiHotCell,
1311                                   gposHotCell,
1312                                   -1);
1313           // push vectors into the Tree
1314           t_TrkselTkFlag->push_back(selFlags[0]);
1315           t_TrkqltyFlag->push_back(selFlags[1]);
1316           t_TrkMissFlag->push_back(selFlags[2]);
1317           t_TrkPVFlag->push_back(selFlags[3]);
1318           t_TrkNuIsolFlag->push_back(neuIsolFlag);
1319           t_TrkhCone->push_back(hCone);
1320           t_TrkP->push_back(pTrack->p());
1321         }
1322       }
1323     }
1324     if (verbosity_ % 10 > 0)
1325       edm::LogVerbatim("IsoTrack") << "Filling " << t_TrkP->size() << " tracks in TrkRestree out of " << nTracks;
1326   }
1327   TrkResTree_->Fill();
1328 }
1329 
1330 void IsoTrig::studyTiming(const edm::Event &theEvent) {
1331   t_timeL2Prod->clear();
1332   t_nPixCand->clear();
1333   t_nPixSeed->clear();
1334 
1335   if (verbosity_ % 10 > 0) {
1336     edm::Handle<SeedingLayerSetsHits> hblayers, helayers;
1337     theEvent.getByToken(tok_SeedingLayerHB_, hblayers);
1338     theEvent.getByToken(tok_SeedingLayerHE_, helayers);
1339     const SeedingLayerSetsHits *layershb = hblayers.product();
1340     const SeedingLayerSetsHits *layershe = helayers.product();
1341     edm::LogVerbatim("IsoTrack") << "size of Seeding TripletLayers hb/he " << layershb->size() << "/"
1342                                  << layershe->size();
1343     edm::Handle<SiPixelRecHitCollection> rchts;
1344     theEvent.getByToken(tok_SiPixelRecHits_, rchts);
1345     const SiPixelRecHitCollection *rechits = rchts.product();
1346     edm::LogVerbatim("IsoTrack") << "size of SiPixelRechits " << rechits->size();
1347   }
1348   int nCandHB = pixelTrackRefsHB_.size(), nCandHE = pixelTrackRefsHE_.size();
1349   int nSeedHB = 0, nSeedHE = 0;
1350 
1351   if (nCandHE > 0) {
1352     edm::Handle<reco::VertexCollection> pVertHB, pVertHE;
1353     theEvent.getByToken(tok_verthb_, pVertHB);
1354     theEvent.getByToken(tok_verthe_, pVertHE);
1355     edm::Handle<trigger::TriggerFilterObjectWithRefs> l1trigobj;
1356     theEvent.getByToken(tok_l1cand_, l1trigobj);
1357 
1358     std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1tauobjref;
1359     std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1jetobjref;
1360     std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1forjetobjref;
1361 
1362     l1trigobj->getObjects(trigger::TriggerL1TauJet, l1tauobjref);
1363     l1trigobj->getObjects(trigger::TriggerL1CenJet, l1jetobjref);
1364     l1trigobj->getObjects(trigger::TriggerL1ForJet, l1forjetobjref);
1365 
1366     double ptTriggered = -10;
1367     double etaTriggered = -100;
1368     double phiTriggered = -100;
1369     for (unsigned int p = 0; p < l1tauobjref.size(); p++) {
1370       if (l1tauobjref[p]->pt() > ptTriggered) {
1371         ptTriggered = l1tauobjref[p]->pt();
1372         phiTriggered = l1tauobjref[p]->phi();
1373         etaTriggered = l1tauobjref[p]->eta();
1374       }
1375     }
1376     for (unsigned int p = 0; p < l1jetobjref.size(); p++) {
1377       if (l1jetobjref[p]->pt() > ptTriggered) {
1378         ptTriggered = l1jetobjref[p]->pt();
1379         phiTriggered = l1jetobjref[p]->phi();
1380         etaTriggered = l1jetobjref[p]->eta();
1381       }
1382     }
1383     for (unsigned int p = 0; p < l1forjetobjref.size(); p++) {
1384       if (l1forjetobjref[p]->pt() > ptTriggered) {
1385         ptTriggered = l1forjetobjref[p]->pt();
1386         phiTriggered = l1forjetobjref[p]->phi();
1387         etaTriggered = l1forjetobjref[p]->eta();
1388       }
1389     }
1390     for (unsigned iS = 0; iS < pixelTrackRefsHE_.size(); iS++) {
1391       reco::VertexCollection::const_iterator vitSel;
1392       double minDZ = 100;
1393       bool vtxMatch;
1394       for (reco::VertexCollection::const_iterator vit = pVertHE->begin(); vit != pVertHE->end(); vit++) {
1395         if (fabs(pixelTrackRefsHE_[iS]->dz(vit->position())) < minDZ) {
1396           minDZ = fabs(pixelTrackRefsHE_[iS]->dz(vit->position()));
1397           vitSel = vit;
1398         }
1399       }
1400       //cut on dYX:
1401       if ((fabs(pixelTrackRefsHE_[iS]->dxy(vitSel->position())) < vtxCutSeed_) || (minDZ == 100))
1402         vtxMatch = true;
1403 
1404       //select tracks not matched to triggered L1 jet
1405       double R = deltaR(etaTriggered, phiTriggered, pixelTrackRefsHE_[iS]->eta(), pixelTrackRefsHE_[iS]->phi());
1406       if (R > tauUnbiasCone_ && vtxMatch)
1407         nSeedHE++;
1408     }
1409     for (unsigned iS = 0; iS < pixelTrackRefsHB_.size(); iS++) {
1410       reco::VertexCollection::const_iterator vitSel;
1411       double minDZ = 100;
1412       bool vtxMatch(false);
1413       for (reco::VertexCollection::const_iterator vit = pVertHB->begin(); vit != pVertHB->end(); vit++) {
1414         if (fabs(pixelTrackRefsHB_[iS]->dz(vit->position())) < minDZ) {
1415           minDZ = fabs(pixelTrackRefsHB_[iS]->dz(vit->position()));
1416           vitSel = vit;
1417         }
1418       }
1419       //cut on dYX:
1420       if ((fabs(pixelTrackRefsHB_[iS]->dxy(vitSel->position())) < 101.0) || (minDZ == 100))
1421         vtxMatch = true;
1422 
1423       //select tracks not matched to triggered L1 jet
1424       double R = deltaR(etaTriggered, phiTriggered, pixelTrackRefsHB_[iS]->eta(), pixelTrackRefsHB_[iS]->phi());
1425       if (R > 1.2 && vtxMatch)
1426         nSeedHB++;
1427     }
1428 
1429     if (verbosity_ % 10 > 0) {
1430       edm::LogVerbatim("IsoTrack") << "(HB/HE) nCand: " << nCandHB << "/" << nCandHE << "nSeed: " << nSeedHB << "/"
1431                                    << nSeedHE;
1432     }
1433   }
1434   t_nPixSeed->push_back(nSeedHB);
1435   t_nPixSeed->push_back(nSeedHE);
1436   t_nPixCand->push_back(nCandHB);
1437   t_nPixCand->push_back(nCandHE);
1438 
1439   TimingTree_->Fill();
1440 }
1441 void IsoTrig::studyMipCut(edm::Handle<reco::TrackCollection> &trkCollection,
1442                           edm::Handle<reco::IsolatedPixelTrackCandidateCollection> &L2cands) {
1443   clearMipCutTreeVectors();
1444   if (verbosity_ % 10 > 0)
1445     edm::LogVerbatim("IsoTrack") << "inside studymipcut";
1446   if (!trkCollection.isValid()) {
1447     edm::LogWarning("IsoTrack") << "trkCollection.isValid is false";
1448   } else {
1449     std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1450     std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1451     spr::propagateCALO(
1452         trkCollection, geo_, bField_, theTrackQuality_, trkCaloDirections1, ((verbosity_ / 100) % 10 > 2));
1453     if (verbosity_ % 10 > 0)
1454       edm::LogVerbatim("IsoTrack") << "Number of L2cands:" << L2cands->size() << " to be matched to something out of "
1455                                    << trkCaloDirections1.size() << " reco tracks";
1456     for (unsigned int i = 0; i < L2cands->size(); i++) {
1457       edm::Ref<reco::IsolatedPixelTrackCandidateCollection> candref =
1458           edm::Ref<reco::IsolatedPixelTrackCandidateCollection>(L2cands, i);
1459       double enIn = candref->energyIn();
1460       h_EnIn->Fill(candref->energyIn());
1461       h_EnOut->Fill(candref->energyOut());
1462       math::XYZTLorentzVector v1(
1463           candref->track()->px(), candref->track()->py(), candref->track()->pz(), candref->track()->p());
1464       if (verbosity_ % 10 > 1)
1465         edm::LogVerbatim("IsoTrack") << "HLT Level Candidate eta/phi/pt/enIn:" << candref->track()->eta() << "/"
1466                                      << candref->track()->phi() << "/" << candref->track()->pt() << "/"
1467                                      << candref->energyIn();
1468       math::XYZTLorentzVector mindRvec;
1469       double mindR = 999.9, mindP1 = 999.9, eMipDR = 0.0;
1470       std::vector<bool> selFlags;
1471       unsigned int nTracks = 0;
1472       double conehmaxNearP = 0, hCone = 0;
1473       for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++, nTracks++) {
1474         const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1475         math::XYZTLorentzVector v2(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
1476         double dr = dR(v1, v2);
1477         double dp1 = std::abs(1. / v1.r() - 1. / v2.r());
1478         if (verbosity_ % 1000 > 0)
1479           edm::LogVerbatim("IsoTrack") << "This recotrack(eta/phi/pt) " << pTrack->eta() << "/" << pTrack->phi() << "/"
1480                                        << pTrack->pt() << " has dr/dp= " << dr << "/" << dp1;
1481         if (dr < mindR) {
1482           mindR = dr;
1483           mindP1 = dp1;
1484           mindRvec = v2;
1485           int nRH_eMipDR = 0, nNearTRKs = 0;
1486           eMipDR = spr::eCone_ecal(geo_,
1487                                    barrelRecHitsHandle_,
1488                                    endcapRecHitsHandle_,
1489                                    trkDetItr->pointHCAL,
1490                                    trkDetItr->pointECAL,
1491                                    a_mipR_,
1492                                    trkDetItr->directionECAL,
1493                                    nRH_eMipDR);
1494           bool selectTk = spr::goodTrack(pTrack, leadPV_, selectionParameters_, ((verbosity_ / 100) % 10 > 1));
1495           spr::trackSelectionParameters oneCutParameters = selectionParameters_;
1496           oneCutParameters.maxDxyPV = 10;
1497           oneCutParameters.maxDzPV = 100;
1498           oneCutParameters.maxInMiss = 2;
1499           oneCutParameters.maxOutMiss = 2;
1500           bool qltyFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1501           oneCutParameters = selectionParameters_;
1502           oneCutParameters.maxDxyPV = 10;
1503           oneCutParameters.maxDzPV = 100;
1504           bool qltyMissFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1505           oneCutParameters = selectionParameters_;
1506           oneCutParameters.maxInMiss = 2;
1507           oneCutParameters.maxOutMiss = 2;
1508           bool qltyPVFlag = spr::goodTrack(pTrack, leadPV_, oneCutParameters, ((verbosity_ / 100) % 10 > 1));
1509           conehmaxNearP = spr::chargeIsolationCone(
1510               nTracks, trkCaloDirections1, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 1));
1511           double e1 = spr::eCone_ecal(geo_,
1512                                       barrelRecHitsHandle_,
1513                                       endcapRecHitsHandle_,
1514                                       trkDetItr->pointHCAL,
1515                                       trkDetItr->pointECAL,
1516                                       a_neutR1_,
1517                                       trkDetItr->directionECAL,
1518                                       nRH_eMipDR);
1519           double e2 = spr::eCone_ecal(geo_,
1520                                       barrelRecHitsHandle_,
1521                                       endcapRecHitsHandle_,
1522                                       trkDetItr->pointHCAL,
1523                                       trkDetItr->pointECAL,
1524                                       a_neutR2_,
1525                                       trkDetItr->directionECAL,
1526                                       nRH_eMipDR);
1527           double e_inCone = e2 - e1;
1528           bool chgIsolFlag = (conehmaxNearP < cutCharge_);
1529           bool mipFlag = (eMipDR < cutMip_);
1530           bool neuIsolFlag = (e_inCone < cutNeutral_);
1531           bool trkpropFlag = ((trkDetItr->okECAL) && (trkDetItr->okHCAL));
1532           selFlags.clear();
1533           selFlags.push_back(selectTk);
1534           selFlags.push_back(qltyFlag);
1535           selFlags.push_back(qltyMissFlag);
1536           selFlags.push_back(qltyPVFlag);
1537           selFlags.push_back(trkpropFlag);
1538           selFlags.push_back(chgIsolFlag);
1539           selFlags.push_back(neuIsolFlag);
1540           selFlags.push_back(mipFlag);
1541           double distFromHotCell = -99.0;
1542           int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1543           GlobalPoint gposHotCell(0., 0., 0.);
1544           std::vector<DetId> coneRecHitDetIds;
1545           hCone = spr::eCone_hcal(geo_,
1546                                   hbhe_,
1547                                   trkDetItr->pointHCAL,
1548                                   trkDetItr->pointECAL,
1549                                   a_coneR_,
1550                                   trkDetItr->directionHCAL,
1551                                   nRecHitsCone,
1552                                   coneRecHitDetIds,
1553                                   distFromHotCell,
1554                                   ietaHotCell,
1555                                   iphiHotCell,
1556                                   gposHotCell,
1557                                   -1);
1558         }
1559       }
1560       pushMipCutTreeVecs(v1, mindRvec, enIn, eMipDR, mindR, mindP1, selFlags, hCone);
1561       fillDifferences(6, v1, mindRvec, (verbosity_ % 10 > 0));
1562       if (mindR >= 0.05) {
1563         fillDifferences(8, v1, mindRvec, (verbosity_ % 10 > 0));
1564         h_MipEnNoMatch->Fill(candref->energyIn(), eMipDR);
1565       } else {
1566         fillDifferences(7, v1, mindRvec, (verbosity_ % 10 > 0));
1567         h_MipEnMatch->Fill(candref->energyIn(), eMipDR);
1568       }
1569     }
1570   }
1571   MipCutTree_->Fill();
1572 }
1573 
1574 void IsoTrig::studyTrigger(edm::Handle<reco::TrackCollection> &trkCollection,
1575                            std::vector<reco::TrackCollection::const_iterator> &goodTks) {
1576   if (verbosity_ % 10 > 0)
1577     edm::LogVerbatim("IsoTrack") << "Inside StudyTrigger";
1578   //// Filling Pt, eta, phi of L1, L2 and L3 objects
1579   for (int j = 0; j < 3; j++) {
1580     for (unsigned int k = 0; k < vec_[j].size(); k++) {
1581       if (verbosity_ % 10 > 0)
1582         edm::LogVerbatim("IsoTrack") << "vec[" << j << "][" << k << "] pt " << vec_[j][k].pt() << " eta "
1583                                      << vec_[j][k].eta() << " phi " << vec_[j][k].phi();
1584       fillHist(j, vec_[j][k]);
1585     }
1586   }
1587 
1588   double deta, dphi, dr;
1589   //// deta, dphi and dR for leading L1 object with L2 and L3 objects
1590   for (int lvl = 1; lvl < 3; lvl++) {
1591     for (unsigned int i = 0; i < vec_[lvl].size(); i++) {
1592       deta = dEta(vec_[0][0], vec_[lvl][i]);
1593       dphi = dPhi(vec_[0][0], vec_[lvl][i]);
1594       dr = dR(vec_[0][0], vec_[lvl][i]);
1595       if (verbosity_ % 10 > 1)
1596         edm::LogVerbatim("IsoTrack") << "lvl " << lvl << " i " << i << " deta " << deta << " dphi " << dphi << " dR "
1597                                      << dr;
1598       h_dEtaL1[lvl - 1]->Fill(deta);
1599       h_dPhiL1[lvl - 1]->Fill(dphi);
1600       h_dRL1[lvl - 1]->Fill(std::sqrt(dr));
1601     }
1602   }
1603 
1604   math::XYZTLorentzVector mindRvec;
1605   double mindR;
1606   for (unsigned int k = 0; k < vec_[2].size(); ++k) {
1607     //// Find min of deta/dphi/dR for each of L3 objects with L2 objects
1608     mindR = 999.9;
1609     if (verbosity_ % 10 > 1)
1610       edm::LogVerbatim("IsoTrack") << "L3obj: pt " << vec_[2][k].pt() << " eta " << vec_[2][k].eta() << " phi "
1611                                    << vec_[2][k].phi();
1612     for (unsigned int j = 0; j < vec_[1].size(); j++) {
1613       dr = dR(vec_[2][k], vec_[1][j]);
1614       if (dr < mindR) {
1615         mindR = dr;
1616         mindRvec = vec_[1][j];
1617       }
1618     }
1619     fillDifferences(0, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1620     if (mindR < 0.03) {
1621       fillDifferences(1, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1622       fillHist(6, mindRvec);
1623       fillHist(8, vec_[2][k]);
1624     } else {
1625       fillDifferences(2, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1626       fillHist(7, mindRvec);
1627       fillHist(9, vec_[2][k]);
1628     }
1629 
1630     ////// Minimum deltaR for each of L3 objs with Reco::tracks
1631     mindR = 999.9;
1632     if (verbosity_ % 10 > 0)
1633       edm::LogVerbatim("IsoTrack") << "Now Matching L3 track with reco: L3 Track (eta, phi) " << vec_[2][k].eta() << ":"
1634                                    << vec_[2][k].phi() << " L2 Track " << mindRvec.eta() << ":" << mindRvec.phi()
1635                                    << " dR " << mindR;
1636     reco::TrackCollection::const_iterator goodTk = trkCollection->end();
1637     if (trkCollection.isValid()) {
1638       double mindP(9999.9);
1639       reco::TrackCollection::const_iterator trkItr;
1640       for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++) {
1641         math::XYZTLorentzVector v4(trkItr->px(), trkItr->py(), trkItr->pz(), trkItr->p());
1642         double deltaR = dR(v4, vec_[2][k]);
1643         double dp = std::abs(v4.r() / vec_[2][k].r() - 1.0);
1644         if (deltaR < mindR) {
1645           mindR = deltaR;
1646           mindP = dp;
1647           mindRvec = v4;
1648           goodTk = trkItr;
1649         }
1650         if ((verbosity_ / 10) % 10 > 1 && deltaR < 1.0)
1651           edm::LogVerbatim("IsoTrack") << "reco track: pt " << v4.pt() << " eta " << v4.eta() << " phi " << v4.phi()
1652                                        << " DR " << deltaR;
1653       }
1654       if (verbosity_ % 10 > 0)
1655         edm::LogVerbatim("IsoTrack") << "Now Matching at Reco level in step 1 DR: " << mindR << ":" << mindP
1656                                      << " eta:phi " << mindRvec.eta() << ":" << mindRvec.phi();
1657       if (mindR < 0.03 && mindP > 0.1) {
1658         for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++) {
1659           math::XYZTLorentzVector v4(trkItr->px(), trkItr->py(), trkItr->pz(), trkItr->p());
1660           double deltaR = dR(v4, vec_[2][k]);
1661           double dp = std::abs(v4.r() / vec_[2][k].r() - 1.0);
1662           if (dp < mindP && deltaR < 0.03) {
1663             mindR = deltaR;
1664             mindP = dp;
1665             mindRvec = v4;
1666             goodTk = trkItr;
1667           }
1668         }
1669         if (verbosity_ % 10 > 0)
1670           edm::LogVerbatim("IsoTrack") << "Now Matching at Reco level in step 2 DR: " << mindR << ":" << mindP
1671                                        << " eta:phi " << mindRvec.eta() << ":" << mindRvec.phi();
1672       }
1673       fillDifferences(3, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1674       fillHist(3, mindRvec);
1675       if (mindR < 0.03) {
1676         fillDifferences(4, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1677         fillHist(4, mindRvec);
1678       } else {
1679         fillDifferences(5, vec_[2][k], mindRvec, (verbosity_ % 10 > 0));
1680         fillHist(5, mindRvec);
1681       }
1682       if (goodTk != trkCollection->end())
1683         goodTks.push_back(goodTk);
1684     }
1685   }
1686 }
1687 
1688 void IsoTrig::studyIsolation(edm::Handle<reco::TrackCollection> &trkCollection,
1689                              std::vector<reco::TrackCollection::const_iterator> &goodTks) {
1690   if (trkCollection.isValid()) {
1691     std::vector<spr::propagatedTrackDirection> trkCaloDirections;
1692     spr::propagateCALO(
1693         trkCollection, geo_, bField_, theTrackQuality_, trkCaloDirections, ((verbosity_ / 100) % 10 > 2));
1694 
1695     std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1696     if ((verbosity_ / 1000) % 10 > 1) {
1697       edm::LogVerbatim("IsoTrack") << "n of barrelRecHitsHandle " << barrelRecHitsHandle_->size();
1698       for (EcalRecHitCollection::const_iterator hit = barrelRecHitsHandle_->begin(); hit != barrelRecHitsHandle_->end();
1699            ++hit) {
1700         edm::LogVerbatim("IsoTrack") << "hit : detid(ieta,iphi) " << (EBDetId)(hit->id()) << " time " << hit->time()
1701                                      << " energy " << hit->energy();
1702       }
1703       edm::LogVerbatim("IsoTrack") << "n of endcapRecHitsHandle " << endcapRecHitsHandle_->size();
1704       for (EcalRecHitCollection::const_iterator hit = endcapRecHitsHandle_->begin(); hit != endcapRecHitsHandle_->end();
1705            ++hit) {
1706         edm::LogVerbatim("IsoTrack") << "hit : detid(ieta,iphi) " << (EEDetId)(hit->id()) << " time " << hit->time()
1707                                      << " energy " << hit->energy();
1708       }
1709       edm::LogVerbatim("IsoTrack") << "n of hbhe " << hbhe_->size();
1710       for (HBHERecHitCollection::const_iterator hit = hbhe_->begin(); hit != hbhe_->end(); ++hit) {
1711         edm::LogVerbatim("IsoTrack") << "hit : detid(ieta,iphi) " << hit->id() << " time " << hit->time() << " energy "
1712                                      << hit->energy();
1713       }
1714     }
1715     unsigned int nTracks = 0, ngoodTk = 0, nselTk = 0;
1716     int ieta = 999;
1717     for (trkDetItr = trkCaloDirections.begin(); trkDetItr != trkCaloDirections.end(); trkDetItr++, nTracks++) {
1718       bool l3Track = (std::find(goodTks.begin(), goodTks.end(), trkDetItr->trkItr) != goodTks.end());
1719       const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1720       math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
1721       bool selectTk = spr::goodTrack(pTrack, leadPV_, selectionParameters_, ((verbosity_ / 100) % 10 > 1));
1722       double eMipDR = 9999., e_inCone = 0, conehmaxNearP = 0, mindR = 999.9, hCone = 0;
1723       if (trkDetItr->okHCAL) {
1724         HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
1725         ieta = detId.ieta();
1726       }
1727       for (unsigned k = 0; k < vec_[0].size(); ++k) {
1728         double deltaR = dR(v4, vec_[0][k]);
1729         if (deltaR < mindR)
1730           mindR = deltaR;
1731       }
1732       if ((verbosity_ / 100) % 10 > 1)
1733         edm::LogVerbatim("IsoTrack") << "Track ECAL " << trkDetItr->okECAL << " HCAL " << trkDetItr->okHCAL << " Flag "
1734                                      << selectTk;
1735       if (selectTk && trkDetItr->okECAL && trkDetItr->okHCAL) {
1736         ngoodTk++;
1737         int nRH_eMipDR = 0, nNearTRKs = 0;
1738         double e1 = spr::eCone_ecal(geo_,
1739                                     barrelRecHitsHandle_,
1740                                     endcapRecHitsHandle_,
1741                                     trkDetItr->pointHCAL,
1742                                     trkDetItr->pointECAL,
1743                                     a_neutR1_,
1744                                     trkDetItr->directionECAL,
1745                                     nRH_eMipDR);
1746         double e2 = spr::eCone_ecal(geo_,
1747                                     barrelRecHitsHandle_,
1748                                     endcapRecHitsHandle_,
1749                                     trkDetItr->pointHCAL,
1750                                     trkDetItr->pointECAL,
1751                                     a_neutR2_,
1752                                     trkDetItr->directionECAL,
1753                                     nRH_eMipDR);
1754         eMipDR = spr::eCone_ecal(geo_,
1755                                  barrelRecHitsHandle_,
1756                                  endcapRecHitsHandle_,
1757                                  trkDetItr->pointHCAL,
1758                                  trkDetItr->pointECAL,
1759                                  a_mipR_,
1760                                  trkDetItr->directionECAL,
1761                                  nRH_eMipDR);
1762         conehmaxNearP =
1763             spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 1));
1764         e_inCone = e2 - e1;
1765         double distFromHotCell = -99.0;
1766         int nRecHitsCone = -99, ietaHotCell = -99, iphiHotCell = -99;
1767         GlobalPoint gposHotCell(0., 0., 0.);
1768         std::vector<DetId> coneRecHitDetIds;
1769         hCone = spr::eCone_hcal(geo_,
1770                                 hbhe_,
1771                                 trkDetItr->pointHCAL,
1772                                 trkDetItr->pointECAL,
1773                                 a_coneR_,
1774                                 trkDetItr->directionHCAL,
1775                                 nRecHitsCone,
1776                                 coneRecHitDetIds,
1777                                 distFromHotCell,
1778                                 ietaHotCell,
1779                                 iphiHotCell,
1780                                 gposHotCell,
1781                                 -1);
1782         if (eMipDR < 1.0)
1783           nselTk++;
1784       }
1785       if (l3Track) {
1786         fillHist(10, v4);
1787         if (selectTk) {
1788           fillHist(11, v4);
1789           fillCuts(0, eMipDR, conehmaxNearP, e_inCone, v4, ieta, (mindR > dr_L1_));
1790           if (conehmaxNearP < cutCharge_) {
1791             fillHist(12, v4);
1792             if (eMipDR < cutMip_) {
1793               fillHist(13, v4);
1794               fillEnergy(1, ieta, hCone, eMipDR, v4);
1795               fillEnergy(0, ieta, hCone, eMipDR, v4);
1796               if (e_inCone < cutNeutral_) {
1797                 fillHist(14, v4);
1798                 fillEnergy(4, ieta, hCone, eMipDR, v4);
1799                 fillEnergy(3, ieta, hCone, eMipDR, v4);
1800               }
1801             }
1802           }
1803         }
1804       } else if (doL2L3_) {
1805         fillHist(15, v4);
1806         if (selectTk) {
1807           fillHist(16, v4);
1808           fillCuts(1, eMipDR, conehmaxNearP, e_inCone, v4, ieta, (mindR > dr_L1_));
1809           if (conehmaxNearP < cutCharge_) {
1810             fillHist(17, v4);
1811             if (eMipDR < cutMip_) {
1812               fillHist(18, v4);
1813               fillEnergy(2, ieta, hCone, eMipDR, v4);
1814               fillEnergy(0, ieta, hCone, eMipDR, v4);
1815               if (e_inCone < cutNeutral_) {
1816                 fillHist(19, v4);
1817                 fillEnergy(5, ieta, hCone, eMipDR, v4);
1818                 fillEnergy(3, ieta, hCone, eMipDR, v4);
1819               }
1820             }
1821           }
1822         }
1823       }
1824     }
1825     //   edm::LogVerbatim("IsoTrack") << "Number of tracks selected offline " << nselTk;
1826   }
1827 }
1828 
1829 void IsoTrig::chgIsolation(double &etaTriggered,
1830                            double &phiTriggered,
1831                            edm::Handle<reco::TrackCollection> &trkCollection,
1832                            const edm::Event &theEvent) {
1833   clearChgIsolnTreeVectors();
1834   if (verbosity_ % 10 > 0)
1835     edm::LogVerbatim("IsoTrack") << "Inside chgIsolation() with eta/phi Triggered: " << etaTriggered << "/"
1836                                  << phiTriggered;
1837   std::vector<double> maxP;
1838 
1839   std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
1840   std::vector<spr::propagatedTrackDirection> trkCaloDirections1;
1841   spr::propagateCALO(trkCollection, geo_, bField_, theTrackQuality_, trkCaloDirections1, ((verbosity_ / 100) % 10 > 2));
1842   if (verbosity_ % 10 > 0)
1843     edm::LogVerbatim("IsoTrack") << "Propagated TrkCollection";
1844   for (unsigned int k = 0; k < pixelIsolationConeSizeAtEC_.size(); ++k)
1845     maxP.push_back(0);
1846   unsigned i = pixelTrackRefsHE_.size();
1847   std::vector<std::pair<unsigned int, std::pair<double, double>>> VecSeedsatEC;
1848   //loop to select isolated tracks
1849   for (unsigned iS = 0; iS < pixelTrackRefsHE_.size(); iS++) {
1850     if (pixelTrackRefsHE_[iS]->p() > minPTrackValue_) {
1851       bool vtxMatch = false;
1852       //associate to vertex (in Z)
1853       unsigned int ivSel = recVtxs_->size();
1854       double minDZ = 100;
1855       for (unsigned int iv = 0; iv < recVtxs_->size(); ++iv) {
1856         if (fabs(pixelTrackRefsHE_[iS]->dz((*recVtxs_)[iv].position())) < minDZ) {
1857           minDZ = fabs(pixelTrackRefsHE_[iS]->dz((*recVtxs_)[iv].position()));
1858           ivSel = iv;
1859         }
1860       }
1861       //cut on dYX:
1862       if (ivSel == recVtxs_->size()) {
1863         vtxMatch = true;
1864       } else if (fabs(pixelTrackRefsHE_[iS]->dxy((*recVtxs_)[ivSel].position())) < vtxCutSeed_) {
1865         vtxMatch = true;
1866       }
1867       //select tracks not matched to triggered L1 jet
1868       double R = deltaR(etaTriggered, phiTriggered, pixelTrackRefsHE_[iS]->eta(), pixelTrackRefsHE_[iS]->phi());
1869       if (R > tauUnbiasCone_ && vtxMatch) {
1870         //propagate seed track to ECAL surface:
1871         std::pair<double, double> seedCooAtEC;
1872         // in case vertex is found:
1873         if (minDZ != 100)
1874           seedCooAtEC = GetEtaPhiAtEcal(pixelTrackRefsHE_[iS]->eta(),
1875                                         pixelTrackRefsHE_[iS]->phi(),
1876                                         pixelTrackRefsHE_[iS]->pt(),
1877                                         pixelTrackRefsHE_[iS]->charge(),
1878                                         (*recVtxs_)[ivSel].z());
1879         //in case vertex is not found:
1880         else
1881           seedCooAtEC = GetEtaPhiAtEcal(pixelTrackRefsHE_[iS]->eta(),
1882                                         pixelTrackRefsHE_[iS]->phi(),
1883                                         pixelTrackRefsHE_[iS]->pt(),
1884                                         pixelTrackRefsHE_[iS]->charge(),
1885                                         0);
1886         VecSeedsatEC.push_back(std::make_pair(iS, seedCooAtEC));
1887       }
1888     }
1889   }
1890   for (unsigned int l = 0; l < VecSeedsatEC.size(); l++) {
1891     unsigned int iSeed = VecSeedsatEC[l].first;
1892     math::XYZTLorentzVector v1(pixelTrackRefsHE_[iSeed]->px(),
1893                                pixelTrackRefsHE_[iSeed]->py(),
1894                                pixelTrackRefsHE_[iSeed]->pz(),
1895                                pixelTrackRefsHE_[iSeed]->p());
1896 
1897     for (unsigned int j = 0; j < VecSeedsatEC.size(); j++) {
1898       unsigned int iSurr = VecSeedsatEC[j].first;
1899       if (iSeed != iSurr) {
1900         //define preliminary cone around seed track impact point from which tracks will be extrapolated:
1901         //  edm::Ref<reco::IsolatedPixelTrackCandidateCollection> cand2ref =
1902         //    edm::Ref<reco::IsolatedPixelTrackCandidateCollection>(L2cands, iSurr);
1903         if (deltaR(pixelTrackRefsHE_[iSeed]->eta(),
1904                    pixelTrackRefsHE_[iSeed]->phi(),
1905                    pixelTrackRefsHE_[iSurr]->eta(),
1906                    pixelTrackRefsHE_[iSurr]->phi()) < prelimCone_) {
1907           unsigned int ivSel = recVtxs_->size();
1908           double minDZ2 = 100;
1909           for (unsigned int iv = 0; iv < recVtxs_->size(); ++iv) {
1910             if (fabs(pixelTrackRefsHE_[iSurr]->dz((*recVtxs_)[iv].position())) < minDZ2) {
1911               minDZ2 = fabs(pixelTrackRefsHE_[iSurr]->dz((*recVtxs_)[iv].position()));
1912               ivSel = iv;
1913             }
1914           }
1915           //cut ot dXY:
1916           if (minDZ2 == 100 || fabs(pixelTrackRefsHE_[iSurr]->dxy((*recVtxs_)[ivSel].position())) < vtxCutIsol_) {
1917             //calculate distance at ECAL surface and update isolation:
1918             double dist = getDistInCM(VecSeedsatEC[i].second.first,
1919                                       VecSeedsatEC[i].second.second,
1920                                       VecSeedsatEC[j].second.first,
1921                                       VecSeedsatEC[j].second.second);
1922             for (unsigned int k = 0; k < pixelIsolationConeSizeAtEC_.size(); ++k) {
1923               if (dist < pixelIsolationConeSizeAtEC_[k]) {
1924                 if (pixelTrackRefsHE_[iSurr]->p() > maxP[k])
1925                   maxP[k] = pixelTrackRefsHE_[iSurr]->p();
1926               }
1927             }
1928           }
1929         }
1930       }
1931     }
1932 
1933     double conehmaxNearP = -1;
1934     bool selectTk = false;
1935     double mindR = 999.9;
1936     int nTracks = 0;
1937     math::XYZTLorentzVector mindRvec;
1938     for (trkDetItr = trkCaloDirections1.begin(); trkDetItr != trkCaloDirections1.end(); trkDetItr++, nTracks++) {
1939       int nNearTRKs = 0;
1940       const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1941       math::XYZTLorentzVector v2(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
1942       double dr = dR(v1, v2);
1943       if (dr < mindR) {
1944         selectTk = spr::goodTrack(pTrack, leadPV_, selectionParameters_, ((verbosity_ / 100) % 10 > 1));
1945         conehmaxNearP = spr::chargeIsolationCone(
1946             nTracks, trkCaloDirections1, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 1));
1947         mindR = dr;
1948         mindRvec = v2;
1949       }
1950     }
1951     pushChgIsolnTreeVecs(v1, mindRvec, maxP, conehmaxNearP, selectTk);
1952   }
1953   ChgIsolnTree_->Fill();
1954 }
1955 
1956 void IsoTrig::getGoodTracks(const edm::Event &iEvent, edm::Handle<reco::TrackCollection> &trkCollection) {
1957   t_nGoodTk->clear();
1958   std::vector<int> nGood(4, 0);
1959   if (trkCollection.isValid()) {
1960     std::vector<spr::propagatedTrackDirection> trkCaloDirections;
1961     spr::propagateCALO(
1962         trkCollection, geo_, bField_, theTrackQuality_, trkCaloDirections, ((verbosity_ / 100) % 10 > 2));
1963 
1964     // get the trigger jet
1965     edm::Handle<trigger::TriggerFilterObjectWithRefs> l1trigobj;
1966     iEvent.getByToken(tok_l1cand_, l1trigobj);
1967 
1968     std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1tauobjref;
1969     l1trigobj->getObjects(trigger::TriggerL1TauJet, l1tauobjref);
1970     std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1jetobjref;
1971     l1trigobj->getObjects(trigger::TriggerL1CenJet, l1jetobjref);
1972     std::vector<edm::Ref<l1extra::L1JetParticleCollection>> l1forjetobjref;
1973     l1trigobj->getObjects(trigger::TriggerL1ForJet, l1forjetobjref);
1974 
1975     double ptTriggered(-10), etaTriggered(-100), phiTriggered(-100);
1976     for (unsigned int p = 0; p < l1tauobjref.size(); p++) {
1977       if (l1tauobjref[p]->pt() > ptTriggered) {
1978         ptTriggered = l1tauobjref[p]->pt();
1979         phiTriggered = l1tauobjref[p]->phi();
1980         etaTriggered = l1tauobjref[p]->eta();
1981       }
1982     }
1983     for (unsigned int p = 0; p < l1jetobjref.size(); p++) {
1984       if (l1jetobjref[p]->pt() > ptTriggered) {
1985         ptTriggered = l1jetobjref[p]->pt();
1986         phiTriggered = l1jetobjref[p]->phi();
1987         etaTriggered = l1jetobjref[p]->eta();
1988       }
1989     }
1990     for (unsigned int p = 0; p < l1forjetobjref.size(); p++) {
1991       if (l1forjetobjref[p]->pt() > ptTriggered) {
1992         ptTriggered = l1forjetobjref[p]->pt();
1993         phiTriggered = l1forjetobjref[p]->phi();
1994         etaTriggered = l1forjetobjref[p]->eta();
1995       }
1996     }
1997     double pTriggered = ptTriggered * cosh(etaTriggered);
1998     math::XYZTLorentzVector pTrigger(
1999         ptTriggered * cos(phiTriggered), ptTriggered * sin(phiTriggered), pTriggered * tanh(etaTriggered), pTriggered);
2000 
2001     std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
2002     unsigned int nTracks(0);
2003     for (trkDetItr = trkCaloDirections.begin(); trkDetItr != trkCaloDirections.end(); trkDetItr++, nTracks++) {
2004       const reco::Track *pTrack = &(*(trkDetItr->trkItr));
2005       math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
2006       bool selectTk = spr::goodTrack(pTrack, leadPV_, selectionParameters_, ((verbosity_ / 100) % 10 > 1));
2007       double mindR = dR(v4, pTrigger);
2008       if ((verbosity_ / 100) % 10 > 1)
2009         edm::LogVerbatim("IsoTrack") << "Track ECAL " << trkDetItr->okECAL << " HCAL " << trkDetItr->okHCAL << " Flag "
2010                                      << selectTk;
2011       if (selectTk && trkDetItr->okECAL && trkDetItr->okHCAL && mindR > 1.0) {
2012         int nRH_eMipDR(0), nNearTRKs(0);
2013         double eMipDR = spr::eCone_ecal(geo_,
2014                                         barrelRecHitsHandle_,
2015                                         endcapRecHitsHandle_,
2016                                         trkDetItr->pointHCAL,
2017                                         trkDetItr->pointECAL,
2018                                         a_mipR_,
2019                                         trkDetItr->directionECAL,
2020                                         nRH_eMipDR);
2021         double conehmaxNearP =
2022             spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, ((verbosity_ / 100) % 10 > 1));
2023         if (conehmaxNearP < 2.0 && eMipDR < 1.0) {
2024           if (pTrack->p() >= 20 && pTrack->p() < 30) {
2025             ++nGood[0];
2026           } else if (pTrack->p() >= 30 && pTrack->p() < 40) {
2027             ++nGood[1];
2028           } else if (pTrack->p() >= 40 && pTrack->p() < 60) {
2029             ++nGood[2];
2030           } else if (pTrack->p() >= 60 && pTrack->p() < 100) {
2031             ++nGood[3];
2032           }
2033         }
2034       }
2035     }
2036   }
2037 
2038   for (unsigned int ii = 0; ii < nGood.size(); ++ii)
2039     t_nGoodTk->push_back(nGood[ii]);
2040 }
2041 
2042 void IsoTrig::fillHist(int indx, math::XYZTLorentzVector &vec) {
2043   h_p[indx]->Fill(vec.r());
2044   h_pt[indx]->Fill(vec.pt());
2045   h_eta[indx]->Fill(vec.eta());
2046   h_phi[indx]->Fill(vec.phi());
2047 }
2048 
2049 void IsoTrig::fillDifferences(int indx, math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2, bool debug) {
2050   double dr = dR(vec1, vec2);
2051   double deta = dEta(vec1, vec2);
2052   double dphi = dPhi(vec1, vec2);
2053   double dpt = dPt(vec1, vec2);
2054   double dp = dP(vec1, vec2);
2055   double dinvpt = dinvPt(vec1, vec2);
2056   h_dEta[indx]->Fill(deta);
2057   h_dPhi[indx]->Fill(dphi);
2058   h_dPt[indx]->Fill(dpt);
2059   h_dP[indx]->Fill(dp);
2060   h_dinvPt[indx]->Fill(dinvpt);
2061   h_mindR[indx]->Fill(dr);
2062   if (debug)
2063     edm::LogVerbatim("IsoTrack") << "mindR for index " << indx << " is " << dr << " deta " << deta << " dphi " << dphi
2064                                  << " dpt " << dpt << " dinvpt " << dinvpt;
2065 }
2066 
2067 void IsoTrig::fillCuts(
2068     int indx, double eMipDR, double conehmaxNearP, double e_inCone, math::XYZTLorentzVector &vec, int ieta, bool cut) {
2069   h_eMip[indx]->Fill(eMipDR);
2070   h_eMaxNearP[indx]->Fill(conehmaxNearP);
2071   h_eNeutIso[indx]->Fill(e_inCone);
2072   if ((conehmaxNearP < cutCharge_) && (eMipDR < cutMip_)) {
2073     for (int lim = 0; lim < 5; ++lim) {
2074       if ((vec.r() > pLimits_[lim]) && (vec.r() <= pLimits_[lim + 1])) {
2075         h_etaMipTracks[lim][indx][0]->Fill((double)(ieta));
2076         if (cut)
2077           h_etaMipTracks[lim][indx][1]->Fill((double)(ieta));
2078         if (e_inCone < cutNeutral_) {
2079           h_etaCalibTracks[lim][indx][0]->Fill((double)(ieta));
2080           if (cut)
2081             h_etaCalibTracks[lim][indx][1]->Fill((double)(ieta));
2082         }
2083       }
2084     }
2085   }
2086 }
2087 
2088 void IsoTrig::fillEnergy(int indx, int ieta, double hCone, double eMipDR, math::XYZTLorentzVector &vec) {
2089   int kk(-1);
2090   if (ieta > 0 && ieta < 25)
2091     kk = 23 + ieta;
2092   else if (ieta > -25 && ieta < 0)
2093     kk = -(ieta + 1);
2094   if (kk >= 0 && eMipDR > 0.01 && hCone > 1.0) {
2095     for (int lim = 0; lim < 5; ++lim) {
2096       if ((vec.r() > pLimits_[lim]) && (vec.r() <= pLimits_[lim + 1])) {
2097         h_eHcal[lim][indx][kk]->Fill(hCone);
2098         h_eCalo[lim][indx][kk]->Fill(hCone + eMipDR);
2099       }
2100     }
2101   }
2102 }
2103 
2104 double IsoTrig::dEta(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) { return (vec1.eta() - vec2.eta()); }
2105 
2106 double IsoTrig::dPhi(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) {
2107   double phi1 = vec1.phi();
2108   if (phi1 < 0)
2109     phi1 += 2.0 * M_PI;
2110   double phi2 = vec2.phi();
2111   if (phi2 < 0)
2112     phi2 += 2.0 * M_PI;
2113   double dphi = phi1 - phi2;
2114   if (dphi > M_PI)
2115     dphi -= 2. * M_PI;
2116   else if (dphi < -M_PI)
2117     dphi += 2. * M_PI;
2118   return dphi;
2119 }
2120 
2121 double IsoTrig::dR(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) {
2122   double deta = dEta(vec1, vec2);
2123   double dphi = dPhi(vec1, vec2);
2124   return std::sqrt(deta * deta + dphi * dphi);
2125 }
2126 
2127 double IsoTrig::dPt(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) { return (vec1.pt() - vec2.pt()); }
2128 
2129 double IsoTrig::dP(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) {
2130   return (std::abs(vec1.r() - vec2.r()));
2131 }
2132 
2133 double IsoTrig::dinvPt(math::XYZTLorentzVector &vec1, math::XYZTLorentzVector &vec2) {
2134   return ((1 / vec1.pt()) - (1 / vec2.pt()));
2135 }
2136 
2137 std::pair<double, double> IsoTrig::etaPhiTrigger() {
2138   double eta(0), phi(0), ptmax(0);
2139   for (unsigned int k = 0; k < vec_[0].size(); ++k) {
2140     if (k == 0) {
2141       eta = vec_[0][k].eta();
2142       phi = vec_[0][k].phi();
2143       ptmax = vec_[0][k].pt();
2144     } else if (vec_[0][k].pt() > ptmax) {
2145       eta = vec_[0][k].eta();
2146       phi = vec_[0][k].phi();
2147       ptmax = vec_[0][k].pt();
2148     }
2149   }
2150   return std::pair<double, double>(eta, phi);
2151 }
2152 
2153 std::pair<double, double> IsoTrig::GetEtaPhiAtEcal(double etaIP, double phiIP, double pT, int charge, double vtxZ) {
2154   double deltaPhi = 0;
2155   double etaEC = 100;
2156   double phiEC = 100;
2157 
2158   double Rcurv = 9999999;
2159   if (bfVal_ != 0)
2160     Rcurv = pT * 33.3 * 100 / (bfVal_ * 10);  //r(m)=pT(GeV)*33.3/B(kG)
2161 
2162   double ecDist = zEE_;
2163   double ecRad = rEB_;  //radius of ECAL barrel (cm)
2164   double theta = 2 * atan(exp(-etaIP));
2165   double zNew = 0;
2166   if (theta > 0.5 * M_PI)
2167     theta = M_PI - theta;
2168   if (fabs(etaIP) < 1.479) {
2169     if ((0.5 * ecRad / Rcurv) > 1) {
2170       etaEC = 10000;
2171       deltaPhi = 0;
2172     } else {
2173       deltaPhi = -charge * asin(0.5 * ecRad / Rcurv);
2174       double alpha1 = 2 * asin(0.5 * ecRad / Rcurv);
2175       double z = ecRad / tan(theta);
2176       if (etaIP > 0)
2177         zNew = z * (Rcurv * alpha1) / ecRad + vtxZ;  //new z-coordinate of track
2178       else
2179         zNew = -z * (Rcurv * alpha1) / ecRad + vtxZ;  //new z-coordinate of track
2180       double zAbs = fabs(zNew);
2181       if (zAbs < ecDist) {
2182         etaEC = -log(tan(0.5 * atan(ecRad / zAbs)));
2183         deltaPhi = -charge * asin(0.5 * ecRad / Rcurv);
2184       }
2185       if (zAbs > ecDist) {
2186         zAbs = (fabs(etaIP) / etaIP) * ecDist;
2187         double Zflight = fabs(zAbs - vtxZ);
2188         double alpha = (Zflight * ecRad) / (z * Rcurv);
2189         double Rec = 2 * Rcurv * sin(alpha / 2);
2190         deltaPhi = -charge * alpha / 2;
2191         etaEC = -log(tan(0.5 * atan(Rec / ecDist)));
2192       }
2193     }
2194   } else {
2195     zNew = (fabs(etaIP) / etaIP) * ecDist;
2196     double Zflight = fabs(zNew - vtxZ);
2197     double Rvirt = fabs(Zflight * tan(theta));
2198     double Rec = 2 * Rcurv * sin(Rvirt / (2 * Rcurv));
2199     deltaPhi = -(charge) * (Rvirt / (2 * Rcurv));
2200     etaEC = -log(tan(0.5 * atan(Rec / ecDist)));
2201   }
2202 
2203   if (zNew < 0)
2204     etaEC = -etaEC;
2205   phiEC = phiIP + deltaPhi;
2206 
2207   if (phiEC < -M_PI)
2208     phiEC += 2 * M_PI;
2209   if (phiEC > M_PI)
2210     phiEC -= 2 * M_PI;
2211 
2212   std::pair<double, double> retVal(etaEC, phiEC);
2213   return retVal;
2214 }
2215 
2216 double IsoTrig::getDistInCM(double eta1, double phi1, double eta2, double phi2) {
2217   double Rec;
2218   double theta1 = 2 * atan(exp(-eta1));
2219   double theta2 = 2 * atan(exp(-eta2));
2220   if (fabs(eta1) < 1.479)
2221     Rec = rEB_;  //radius of ECAL barrel
2222   else if (fabs(eta1) > 1.479 && fabs(eta1) < 7.0)
2223     Rec = tan(theta1) * zEE_;  //distance from IP to ECAL endcap
2224   else
2225     return 1000;
2226 
2227   //|vect| times tg of acos(scalar product)
2228   double angle =
2229       acos((sin(theta1) * sin(theta2) * (sin(phi1) * sin(phi2) + cos(phi1) * cos(phi2)) + cos(theta1) * cos(theta2)));
2230   if (angle < 0.5 * M_PI)
2231     return fabs((Rec / sin(theta1)) * tan(angle));
2232   else
2233     return 1000;
2234 }
2235 
2236 //define this as a plug-in
2237 DEFINE_FWK_MODULE(IsoTrig);