Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-26 00:02:10

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 
0022 // Root objects
0023 #include "TROOT.h"
0024 #include "TH1.h"
0025 #include "TH2.h"
0026 #include "TSystem.h"
0027 #include "TFile.h"
0028 #include "TProfile.h"
0029 #include "TDirectory.h"
0030 #include "TTree.h"
0031 #include "TMath.h"
0032 
0033 // user include files
0034 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0035 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0036 #include "Calibration/IsolatedParticles/interface/eCone.h"
0037 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0038 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0039 
0040 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0041 
0042 #include "DataFormats/Common/interface/RefToBase.h"
0043 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0044 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0045 #include "DataFormats/Math/interface/LorentzVector.h"
0046 #include "DataFormats/Math/interface/Point3D.h"
0047 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0048 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0049 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0050 #include "DataFormats/TrackReco/interface/Track.h"
0051 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0052 #include "DataFormats/TrackReco/interface/TrackBase.h"
0053 //Tracks
0054 #include "DataFormats/TrackReco/interface/HitPattern.h"
0055 // Vertices
0056 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0057 #include "DataFormats/VertexReco/interface/Vertex.h"
0058 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0059 //Triggers
0060 #include "DataFormats/Common/interface/TriggerResults.h"
0061 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0062 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0063 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0064 #include "DataFormats/Luminosity/interface/LumiDetails.h"
0065 
0066 #include "FWCore/Common/interface/TriggerNames.h"
0067 #include "FWCore/Framework/interface/Frameworkfwd.h"
0068 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0069 #include "FWCore/Framework/interface/Event.h"
0070 #include "FWCore/Framework/interface/MakerMacros.h"
0071 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0072 #include "FWCore/ServiceRegistry/interface/Service.h"
0073 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0074 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0075 
0076 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0077 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0078 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0079 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0080 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0081 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0082 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0083 
0084 #include "HLTrigger/HLTcore/interface/HLTPrescaleProvider.h"
0085 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0086 
0087 #include "MagneticField/Engine/interface/MagneticField.h"
0088 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0089 
0090 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0091 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0092 
0093 #include "TrackingTools/TransientTrackingRecHit/interface/SeedingLayerSetsHits.h"
0094 
0095 class IsoTrig : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0096 public:
0097   explicit IsoTrig(const edm::ParameterSet &);
0098   ~IsoTrig() override;
0099 
0100   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0101 
0102 private:
0103   void analyze(const edm::Event &, const edm::EventSetup &) override;
0104   void beginJob() override;
0105   void endJob() override;
0106   void beginRun(edm::Run const &, edm::EventSetup const &) override;
0107   void endRun(edm::Run const &, edm::EventSetup const &) override {}
0108 
0109   void clearMipCutTreeVectors();
0110   void clearChgIsolnTreeVectors();
0111   void pushChgIsolnTreeVecs(math::XYZTLorentzVector &Pixcand,
0112                             math::XYZTLorentzVector &Trkcand,
0113                             std::vector<double> &PixMaxP,
0114                             double &TrkMaxP,
0115                             bool &selTk);
0116   void pushMipCutTreeVecs(math::XYZTLorentzVector &NFcand,
0117                           math::XYZTLorentzVector &Trkcand,
0118                           double &EmipNFcand,
0119                           double &EmipTrkcand,
0120                           double &mindR,
0121                           double &mindP1,
0122                           std::vector<bool> &Flags,
0123                           double hCone);
0124   void StudyTrkEbyP(edm::Handle<reco::TrackCollection> &trkCollection);
0125   void studyTiming(const edm::Event &theEvent);
0126   void studyMipCut(edm::Handle<reco::TrackCollection> &trkCollection,
0127                    edm::Handle<reco::IsolatedPixelTrackCandidateCollection> &L2cands);
0128   void studyTrigger(edm::Handle<reco::TrackCollection> &, std::vector<reco::TrackCollection::const_iterator> &);
0129   void studyIsolation(edm::Handle<reco::TrackCollection> &, std::vector<reco::TrackCollection::const_iterator> &);
0130   void chgIsolation(double &etaTriggered,
0131                     double &phiTriggered,
0132                     edm::Handle<reco::TrackCollection> &trkCollection,
0133                     const edm::Event &theEvent);
0134   void getGoodTracks(const edm::Event &, edm::Handle<reco::TrackCollection> &);
0135   void fillHist(int, math::XYZTLorentzVector &);
0136   void fillDifferences(int, math::XYZTLorentzVector &, math::XYZTLorentzVector &, bool);
0137   void fillCuts(int, double, double, double, math::XYZTLorentzVector &, int, bool);
0138   void fillEnergy(int, int, double, double, math::XYZTLorentzVector &);
0139   double dEta(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0140   double dPhi(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0141   double dR(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0142   double dPt(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0143   double dP(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0144   double dinvPt(math::XYZTLorentzVector &, math::XYZTLorentzVector &);
0145   std::pair<double, double> etaPhiTrigger();
0146   std::pair<double, double> GetEtaPhiAtEcal(double etaIP, double phiIP, double pT, int charge, double vtxZ);
0147   double getDistInCM(double eta1, double phi1, double eta2, double phi2);
0148 
0149   // ----------member data ---------------------------
0150   HLTPrescaleProvider hltPrescaleProvider_;
0151   const std::vector<std::string> trigNames_;
0152   const edm::InputTag pixCandTag_, l1CandTag_, l2CandTag_;
0153   const std::vector<edm::InputTag> pixelTracksSources_;
0154   const bool doL2L3_, doTiming_, doMipCutTree_;
0155   const bool doTrkResTree_, doChgIsolTree_, doStudyIsol_;
0156   const int verbosity_;
0157   const std::vector<double> pixelIsolationConeSizeAtEC_;
0158   const double minPTrackValue_, vtxCutSeed_, vtxCutIsol_;
0159   const double tauUnbiasCone_, prelimCone_;
0160   std::string theTrackQuality_;
0161   const std::string processName_;
0162   double rEB_, zEE_, bfVal_;
0163   spr::trackSelectionParameters selectionParameters_;
0164   const double dr_L1_, a_coneR_, a_charIsoR_, a_neutIsoR_;
0165   const double a_mipR_, a_neutR1_, a_neutR2_, cutMip_;
0166   const double cutCharge_, cutNeutral_;
0167   const int minRunNo_, maxRunNo_;
0168   edm::EDGetTokenT<LumiDetails> tok_lumi_;
0169   edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
0170   edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0171   edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_hlt_;
0172   edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0173   edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0174   edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0175   edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0176   edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0177   edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0178   edm::EDGetTokenT<reco::VertexCollection> tok_verthb_, tok_verthe_;
0179   edm::EDGetTokenT<SeedingLayerSetsHits> tok_SeedingLayerHB_;
0180   edm::EDGetTokenT<SeedingLayerSetsHits> tok_SeedingLayerHE_;
0181   edm::EDGetTokenT<SiPixelRecHitCollection> tok_SiPixelRecHits_;
0182   edm::EDGetTokenT<reco::IsolatedPixelTrackCandidateCollection> tok_pixtk_;
0183   edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_l1cand_;
0184   edm::EDGetTokenT<reco::IsolatedPixelTrackCandidateCollection> tok_l2cand_;
0185   std::vector<edm::EDGetTokenT<reco::TrackCollection>> tok_pixtks_;
0186 
0187   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0188   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0189 
0190   std::vector<reco::TrackRef> pixelTrackRefsHB_, pixelTrackRefsHE_;
0191   edm::Handle<HBHERecHitCollection> hbhe_;
0192   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle_;
0193   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle_;
0194   edm::Handle<reco::BeamSpot> beamSpotH_;
0195   edm::Handle<reco::VertexCollection> recVtxs_;
0196 
0197   const MagneticField *bField_;
0198   const CaloGeometry *geo_;
0199   math::XYZPoint leadPV_;
0200 
0201   std::map<unsigned int, unsigned int> trigList_;
0202   std::map<unsigned int, const std::pair<int, int>> trigPreList_;
0203   bool changed_;
0204   double pLimits_[6];
0205   edm::Service<TFileService> fs_;
0206   TTree *MipCutTree_, *ChgIsolnTree_, *TrkResTree_, *TimingTree_;
0207   std::vector<double> *t_timeL2Prod;
0208   std::vector<int> *t_nPixCand;
0209   std::vector<int> *t_nPixSeed;
0210   std::vector<int> *t_nGoodTk;
0211 
0212   std::vector<double> *t_TrkhCone;
0213   std::vector<double> *t_TrkP;
0214   std::vector<bool> *t_TrkselTkFlag;
0215   std::vector<bool> *t_TrkqltyFlag;
0216   std::vector<bool> *t_TrkMissFlag;
0217   std::vector<bool> *t_TrkPVFlag;
0218   std::vector<bool> *t_TrkNuIsolFlag;
0219 
0220   std::vector<double> *t_PixcandP;
0221   std::vector<double> *t_PixcandPt;
0222   std::vector<double> *t_PixcandEta;
0223   std::vector<double> *t_PixcandPhi;
0224   std::vector<std::vector<double>> *t_PixcandMaxP;
0225   std::vector<double> *t_PixTrkcandP;
0226   std::vector<double> *t_PixTrkcandPt;
0227   std::vector<double> *t_PixTrkcandEta;
0228   std::vector<double> *t_PixTrkcandPhi;
0229   std::vector<double> *t_PixTrkcandMaxP;
0230   std::vector<bool> *t_PixTrkcandselTk;
0231 
0232   std::vector<double> *t_NFcandP;
0233   std::vector<double> *t_NFcandPt;
0234   std::vector<double> *t_NFcandEta;
0235   std::vector<double> *t_NFcandPhi;
0236   std::vector<double> *t_NFcandEmip;
0237   std::vector<double> *t_NFTrkcandP;
0238   std::vector<double> *t_NFTrkcandPt;
0239   std::vector<double> *t_NFTrkcandEta;
0240   std::vector<double> *t_NFTrkcandPhi;
0241   std::vector<double> *t_NFTrkcandEmip;
0242   std::vector<double> *t_NFTrkMinDR;
0243   std::vector<double> *t_NFTrkMinDP1;
0244   std::vector<bool> *t_NFTrkselTkFlag;
0245   std::vector<bool> *t_NFTrkqltyFlag;
0246   std::vector<bool> *t_NFTrkMissFlag;
0247   std::vector<bool> *t_NFTrkPVFlag;
0248   std::vector<bool> *t_NFTrkPropFlag;
0249   std::vector<bool> *t_NFTrkChgIsoFlag;
0250   std::vector<bool> *t_NFTrkNeuIsoFlag;
0251   std::vector<bool> *t_NFTrkMipFlag;
0252   std::vector<double> *t_ECone;
0253 
0254   TH1D *h_EnIn, *h_EnOut;
0255   TH2D *h_MipEnMatch, *h_MipEnNoMatch;
0256   TH1I *h_nHLT, *h_HLT, *h_PreL1, *h_PreHLT;
0257   TH1I *h_Pre, *h_nL3Objs, *h_Filters;
0258   TH1D *h_PreL1wt, *h_PreHLTwt, *h_L1ObjEnergy;
0259   TH1D *h_p[20], *h_pt[20], *h_eta[20], *h_phi[20];
0260   TH1D *h_dEtaL1[2], *h_dPhiL1[2], *h_dRL1[2];
0261   TH1D *h_dEta[9], *h_dPhi[9], *h_dPt[9], *h_dP[9];
0262   TH1D *h_dinvPt[9], *h_mindR[9], *h_eMip[2];
0263   TH1D *h_eMaxNearP[2], *h_eNeutIso[2];
0264   TH1D *h_etaCalibTracks[5][2][2], *h_etaMipTracks[5][2][2];
0265   TH1D *h_eHcal[5][6][48], *h_eCalo[5][6][48];
0266   TH1I *g_Pre, *g_PreL1, *g_PreHLT, *g_Accepts;
0267   std::vector<math::XYZTLorentzVector> vec_[3];
0268 };
0269 
0270 IsoTrig::IsoTrig(const edm::ParameterSet &iConfig)
0271     : hltPrescaleProvider_(iConfig, consumesCollector(), *this),
0272       trigNames_(iConfig.getUntrackedParameter<std::vector<std::string>>("Triggers")),
0273       pixCandTag_(iConfig.getUntrackedParameter<edm::InputTag>("pixCandTag")),
0274       l1CandTag_(iConfig.getUntrackedParameter<edm::InputTag>("l1CandTag")),
0275       l2CandTag_(iConfig.getUntrackedParameter<edm::InputTag>("l2CandTag")),
0276       pixelTracksSources_(iConfig.getUntrackedParameter<std::vector<edm::InputTag>>("pixelTracksSources")),
0277       doL2L3_(iConfig.getUntrackedParameter<bool>("doL2L3", true)),
0278       doTiming_(iConfig.getUntrackedParameter<bool>("doTimingTree", true)),
0279       doMipCutTree_(iConfig.getUntrackedParameter<bool>("doMipCutTree", true)),
0280       doTrkResTree_(iConfig.getUntrackedParameter<bool>("doTrkResTree", true)),
0281       doChgIsolTree_(iConfig.getUntrackedParameter<bool>("doChgIsolTree", true)),
0282       doStudyIsol_(iConfig.getUntrackedParameter<bool>("doStudyIsol", true)),
0283       verbosity_(iConfig.getUntrackedParameter<int>("verbosity", 0)),
0284       pixelIsolationConeSizeAtEC_(iConfig.getUntrackedParameter<std::vector<double>>("pixelIsolationConeSizeAtEC")),
0285       minPTrackValue_(iConfig.getUntrackedParameter<double>("minPTrackValue")),
0286       vtxCutSeed_(iConfig.getUntrackedParameter<double>("vertexCutSeed")),
0287       vtxCutIsol_(iConfig.getUntrackedParameter<double>("vertexCutIsol")),
0288       tauUnbiasCone_(iConfig.getUntrackedParameter<double>("tauUnbiasCone")),
0289       prelimCone_(iConfig.getUntrackedParameter<double>("prelimCone")),
0290       theTrackQuality_(iConfig.getUntrackedParameter<std::string>("trackQuality", "highPurity")),
0291       processName_(iConfig.getUntrackedParameter<std::string>("processName", "HLT")),
0292       dr_L1_(iConfig.getUntrackedParameter<double>("isolationL1", 1.0)),
0293       a_coneR_(iConfig.getUntrackedParameter<double>("coneRadius", 34.98)),
0294       a_charIsoR_(a_coneR_ + 28.9),
0295       a_neutIsoR_(a_charIsoR_ * 0.726),
0296       a_mipR_(iConfig.getUntrackedParameter<double>("coneRadiusMIP", 14.0)),
0297       a_neutR1_(iConfig.getUntrackedParameter<double>("coneRadiusNeut1", 21.0)),
0298       a_neutR2_(iConfig.getUntrackedParameter<double>("coneRadiusNeut2", 29.0)),
0299       cutMip_(iConfig.getUntrackedParameter<double>("cutMIP", 1.0)),
0300       cutCharge_(iConfig.getUntrackedParameter<double>("chargeIsolation", 2.0)),
0301       cutNeutral_(iConfig.getUntrackedParameter<double>("neutralIsolation", 2.0)),
0302       minRunNo_(iConfig.getUntrackedParameter<int>("minRun")),
0303       maxRunNo_(iConfig.getUntrackedParameter<int>("maxRun")),
0304       changed_(false),
0305       t_timeL2Prod(nullptr),
0306       t_nPixCand(nullptr),
0307       t_nPixSeed(nullptr),
0308       t_nGoodTk(nullptr),
0309       t_TrkhCone(nullptr),
0310       t_TrkP(nullptr),
0311       t_TrkselTkFlag(nullptr),
0312       t_TrkqltyFlag(nullptr),
0313       t_TrkMissFlag(nullptr),
0314       t_TrkPVFlag(nullptr),
0315       t_TrkNuIsolFlag(nullptr),
0316       t_PixcandP(nullptr),
0317       t_PixcandPt(nullptr),
0318       t_PixcandEta(nullptr),
0319       t_PixcandPhi(nullptr),
0320       t_PixcandMaxP(nullptr),
0321       t_PixTrkcandP(nullptr),
0322       t_PixTrkcandPt(nullptr),
0323       t_PixTrkcandEta(nullptr),
0324       t_PixTrkcandPhi(nullptr),
0325       t_PixTrkcandMaxP(nullptr),
0326       t_PixTrkcandselTk(nullptr),
0327       t_NFcandP(nullptr),
0328       t_NFcandPt(nullptr),
0329       t_NFcandEta(nullptr),
0330       t_NFcandPhi(nullptr),
0331       t_NFcandEmip(nullptr),
0332       t_NFTrkcandP(nullptr),
0333       t_NFTrkcandPt(nullptr),
0334       t_NFTrkcandEta(nullptr),
0335       t_NFTrkcandPhi(nullptr),
0336       t_NFTrkcandEmip(nullptr),
0337       t_NFTrkMinDR(nullptr),
0338       t_NFTrkMinDP1(nullptr),
0339       t_NFTrkselTkFlag(nullptr),
0340       t_NFTrkqltyFlag(nullptr),
0341       t_NFTrkMissFlag(nullptr),
0342       t_NFTrkPVFlag(nullptr),
0343       t_NFTrkPropFlag(nullptr),
0344       t_NFTrkChgIsoFlag(nullptr),
0345       t_NFTrkNeuIsoFlag(nullptr),
0346       t_NFTrkMipFlag(nullptr),
0347       t_ECone(nullptr) {
0348   usesResource(TFileService::kSharedResource);
0349 
0350   //now do whatever initialization is needed
0351   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality_);
0352   selectionParameters_.minPt = iConfig.getUntrackedParameter<double>("minTrackPt", 10.0);
0353   selectionParameters_.minQuality = trackQuality_;
0354   selectionParameters_.maxDxyPV = iConfig.getUntrackedParameter<double>("maxDxyPV", 0.2);
0355   selectionParameters_.maxDzPV = iConfig.getUntrackedParameter<double>("maxDzPV", 5.0);
0356   selectionParameters_.maxChi2 = iConfig.getUntrackedParameter<double>("maxChi2", 5.0);
0357   selectionParameters_.maxDpOverP = iConfig.getUntrackedParameter<double>("maxDpOverP", 0.1);
0358   selectionParameters_.minOuterHit = iConfig.getUntrackedParameter<int>("minOuterHit", 4);
0359   selectionParameters_.minLayerCrossed = iConfig.getUntrackedParameter<int>("minLayerCrossed", 8);
0360   selectionParameters_.maxInMiss = iConfig.getUntrackedParameter<int>("maxInMiss", 0);
0361   selectionParameters_.maxOutMiss = iConfig.getUntrackedParameter<int>("maxOutMiss", 0);
0362 
0363   // define tokens for access
0364   tok_lumi_ = consumes<LumiDetails, edm::InLumi>(edm::InputTag("lumiProducer"));
0365   edm::InputTag triggerEvent_("hltTriggerSummaryAOD", "", processName_);
0366   tok_trigEvt_ = consumes<trigger::TriggerEvent>(triggerEvent_);
0367   edm::InputTag theTriggerResultsLabel("TriggerResults", "", processName_);
0368   tok_trigRes_ = consumes<edm::TriggerResults>(theTriggerResultsLabel);
0369   tok_genTrack_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
0370   tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0371   tok_bs_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
0372   tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0373   tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0374   tok_hbhe_ = consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"));
0375   tok_pixtk_ = consumes<reco::IsolatedPixelTrackCandidateCollection>(pixCandTag_);
0376   tok_l1cand_ = consumes<trigger::TriggerFilterObjectWithRefs>(l1CandTag_);
0377   tok_l2cand_ = consumes<reco::IsolatedPixelTrackCandidateCollection>(l2CandTag_);
0378   if (doTiming_) {
0379     tok_verthb_ = consumes<reco::VertexCollection>(edm::InputTag("hltHITPixelVerticesHB"));
0380     tok_verthe_ = consumes<reco::VertexCollection>(edm::InputTag("hltHITPixelVerticesHB"));
0381     tok_hlt_ = consumes<trigger::TriggerFilterObjectWithRefs>(edm::InputTag("hltL1sL1SingleJet68"));
0382     tok_SeedingLayerHB_ = consumes<SeedingLayerSetsHits>(edm::InputTag("hltPixelLayerTripletsHITHB"));
0383     tok_SeedingLayerHE_ = consumes<SeedingLayerSetsHits>(edm::InputTag("hltPixelLayerTripletsHITHE"));
0384     tok_SiPixelRecHits_ = consumes<SiPixelRecHitCollection>(edm::InputTag("hltSiPixelRecHits"));
0385   }
0386   if (doChgIsolTree_) {
0387     for (unsigned int k = 0; k < pixelTracksSources_.size(); ++k) {
0388       //      edm::InputTag  pix (pixelTracksSources_[k],"",processName_);
0389       //      tok_pixtks_.push_back(consumes<reco::TrackCollection>(pix));
0390       tok_pixtks_.push_back(consumes<reco::TrackCollection>(pixelTracksSources_[k]));
0391     }
0392   }
0393   if (verbosity_ >= 0) {
0394     edm::LogVerbatim("IsoTrack") << "Parameters read from config file \n"
0395                                  << "\t minPt " << selectionParameters_.minPt << "\t theTrackQuality "
0396                                  << theTrackQuality_ << "\t minQuality " << selectionParameters_.minQuality
0397                                  << "\t maxDxyPV " << selectionParameters_.maxDxyPV << "\t maxDzPV "
0398                                  << selectionParameters_.maxDzPV << "\t maxChi2 " << selectionParameters_.maxChi2
0399                                  << "\t maxDpOverP " << selectionParameters_.maxDpOverP << "\t minOuterHit "
0400                                  << selectionParameters_.minOuterHit << "\t minLayerCrossed "
0401                                  << selectionParameters_.minLayerCrossed << "\t maxInMiss "
0402                                  << selectionParameters_.maxInMiss << "\t maxOutMiss "
0403                                  << selectionParameters_.maxOutMiss << "\t a_coneR " << a_coneR_ << "\t a_charIsoR "
0404                                  << a_charIsoR_ << "\t a_neutIsoR " << a_neutIsoR_ << "\t a_mipR " << a_mipR_
0405                                  << "\t a_neutR " << a_neutR1_ << ":" << a_neutR2_ << "\t cuts (MIP " << cutMip_
0406                                  << " : Charge " << cutCharge_ << " : Neutral " << cutNeutral_ << ")";
0407     edm::LogVerbatim("IsoTrack") << "Charge Isolation parameters:"
0408                                  << "\t minPTrackValue " << minPTrackValue_ << "\t vtxCutSeed " << vtxCutSeed_
0409                                  << "\t vtxCutIsol " << vtxCutIsol_ << "\t tauUnbiasCone " << tauUnbiasCone_
0410                                  << "\t prelimCone " << prelimCone_ << "\t pixelIsolationConeSizeAtEC";
0411     for (unsigned int k = 0; k < pixelIsolationConeSizeAtEC_.size(); ++k)
0412       edm::LogVerbatim("IsoTrack") << "[" << k << "] " << pixelIsolationConeSizeAtEC_[k];
0413   }
0414   double pl[] = {20, 30, 40, 60, 80, 120};
0415   for (int i = 0; i < 6; ++i)
0416     pLimits_[i] = pl[i];
0417   rEB_ = 123.8;
0418   zEE_ = 317.0;
0419 
0420   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0421   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0422 }
0423 
0424 IsoTrig::~IsoTrig() {
0425   // do anything here that needs to be done at desctruction time
0426   // (e.g. close files, deallocate resources etc.)
0427   if (t_timeL2Prod)
0428     delete t_timeL2Prod;
0429   if (t_nPixCand)
0430     delete t_nPixCand;
0431   if (t_nPixSeed)
0432     delete t_nPixSeed;
0433   if (t_nGoodTk)
0434     delete t_nGoodTk;
0435   if (t_TrkhCone)
0436     delete t_TrkhCone;
0437   if (t_TrkP)
0438     delete t_TrkP;
0439   if (t_TrkselTkFlag)
0440     delete t_TrkselTkFlag;
0441   if (t_TrkqltyFlag)
0442     delete t_TrkqltyFlag;
0443   if (t_TrkMissFlag)
0444     delete t_TrkMissFlag;
0445   if (t_TrkPVFlag)
0446     delete t_TrkPVFlag;
0447   if (t_TrkNuIsolFlag)
0448     delete t_TrkNuIsolFlag;
0449   if (t_PixcandP)
0450     delete t_PixcandP;
0451   if (t_PixcandPt)
0452     delete t_PixcandPt;
0453   if (t_PixcandEta)
0454     delete t_PixcandEta;
0455   if (t_PixcandPhi)
0456     delete t_PixcandPhi;
0457   if (t_PixcandMaxP)
0458     delete t_PixcandMaxP;
0459   if (t_PixTrkcandP)
0460     delete t_PixTrkcandP;
0461   if (t_PixTrkcandPt)
0462     delete t_PixTrkcandPt;
0463   if (t_PixTrkcandEta)
0464     delete t_PixTrkcandEta;
0465   if (t_PixTrkcandPhi)
0466     delete t_PixTrkcandPhi;
0467   if (t_PixTrkcandMaxP)
0468     delete t_PixTrkcandMaxP;
0469   if (t_PixTrkcandselTk)
0470     delete t_PixTrkcandselTk;
0471   if (t_NFcandP)
0472     delete t_NFcandP;
0473   if (t_NFcandPt)
0474     delete t_NFcandPt;
0475   if (t_NFcandEta)
0476     delete t_NFcandEta;
0477   if (t_NFcandPhi)
0478     delete t_NFcandPhi;
0479   if (t_NFcandEmip)
0480     delete t_NFcandEmip;
0481   if (t_NFTrkcandP)
0482     delete t_NFTrkcandP;
0483   if (t_NFTrkcandPt)
0484     delete t_NFTrkcandPt;
0485   if (t_NFTrkcandEta)
0486     delete t_NFTrkcandEta;
0487   if (t_NFTrkcandPhi)
0488     delete t_NFTrkcandPhi;
0489   if (t_NFTrkcandEmip)
0490     delete t_NFTrkcandEmip;
0491   if (t_NFTrkMinDR)
0492     delete t_NFTrkMinDR;
0493   if (t_NFTrkMinDP1)
0494     delete t_NFTrkMinDP1;
0495   if (t_NFTrkselTkFlag)
0496     delete t_NFTrkselTkFlag;
0497   if (t_NFTrkqltyFlag)
0498     delete t_NFTrkqltyFlag;
0499   if (t_NFTrkMissFlag)
0500     delete t_NFTrkMissFlag;
0501   if (t_NFTrkPVFlag)
0502     delete t_NFTrkPVFlag;
0503   if (t_NFTrkPropFlag)
0504     delete t_NFTrkPropFlag;
0505   if (t_NFTrkChgIsoFlag)
0506     delete t_NFTrkChgIsoFlag;
0507   if (t_NFTrkNeuIsoFlag)
0508     delete t_NFTrkNeuIsoFlag;
0509   if (t_NFTrkMipFlag)
0510     delete t_NFTrkMipFlag;
0511   if (t_ECone)
0512     delete t_ECone;
0513 }
0514 
0515 void IsoTrig::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0516   std::vector<std::string> triggers = {"HLT_IsoTrackHB"};
0517   std::vector<edm::InputTag> tags = {edm::InputTag("hltHITPixelTracksHB"), edm::InputTag("hltHITPixelTracksHE")};
0518   std::vector<double> cones = {35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 63.9, 70.0};
0519   edm::ParameterSetDescription desc;
0520   desc.addUntracked<std::vector<std::string>>("Triggers", triggers);
0521   desc.addUntracked<edm::InputTag>("pixCandTag", edm::InputTag(" "));
0522   desc.addUntracked<edm::InputTag>("l1CandTag", edm::InputTag("hltL1sV0SingleJet60"));
0523   desc.addUntracked<edm::InputTag>("l2CandTag", edm::InputTag("isolEcalPixelTrackProd"));
0524   desc.addUntracked<bool>("doL2L3", false);
0525   desc.addUntracked<bool>("doTimingTree", false);
0526   desc.addUntracked<bool>("doMipCutTree", false);
0527   desc.addUntracked<bool>("doTrkResTree", true);
0528   desc.addUntracked<bool>("doChgIsolTree", false);
0529   desc.addUntracked<bool>("doStudyIsol", false);
0530   desc.addUntracked<int>("verbosity", 0);
0531   desc.addUntracked<std::string>("processName", "HLT");
0532   desc.addUntracked<std::string>("trackQuality", "highPurity");
0533   desc.addUntracked<double>("minTrackPt", 10.0);
0534   desc.addUntracked<double>("maxDxyPV", 0.02);
0535   desc.addUntracked<double>("maxDzPV", 0.02);
0536   desc.addUntracked<double>("maxChi2", 5.0);
0537   desc.addUntracked<double>("maxDpOverP", 0.1);
0538   desc.addUntracked<int>("minOuterHit", 4);
0539   desc.addUntracked<int>("minLayerCrossed", 8);
0540   desc.addUntracked<int>("maxInMiss", 0);
0541   desc.addUntracked<int>("maxOutMiss", 0);
0542   desc.addUntracked<double>("isolationL1", 1.0);
0543   desc.addUntracked<double>("coneRadius", 34.98);
0544   desc.addUntracked<double>("coneRadiusMIP", 14.0);
0545   desc.addUntracked<double>("coneRadiusNeut1", 21.0);
0546   desc.addUntracked<double>("coneRadiusNeut2", 29.0);
0547   desc.addUntracked<double>("cutMIP", 1.0);
0548   desc.addUntracked<double>("chargeIsolation", 2.0);
0549   desc.addUntracked<double>("neutralIsolation", 2.0);
0550   desc.addUntracked<int>("minRun", 190456);
0551   desc.addUntracked<int>("maxRun", 203002);
0552   desc.addUntracked<std::vector<edm::InputTag>>("pixelTracksSources", tags);
0553   desc.addUntracked<std::vector<double>>("pixelIsolationConeSizeAtEC", cones);
0554   desc.addUntracked<double>("minPTrackValue", 0.0);
0555   desc.addUntracked<double>("vertexCutSeed", 101.0);
0556   desc.addUntracked<double>("vertexCutIsol", 101.0);
0557   desc.addUntracked<double>("tauUnbiasCone", 1.2);
0558   desc.addUntracked<double>("prelimCone", 1.0);
0559   desc.add<unsigned int>("stageL1Trigger", 1);
0560   descriptions.add("isoTrigDefault", desc);
0561 }
0562 
0563 void IsoTrig::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0564   if (verbosity_ % 10 > 0)
0565     edm::LogVerbatim("IsoTrack") << "Event starts====================================";
0566 
0567   int RunNo = iEvent.id().run();
0568 
0569   HLTConfigProvider const &hltConfig = hltPrescaleProvider_.hltConfigProvider();
0570 
0571   bField_ = &iSetup.getData(tok_magField_);
0572   geo_ = &iSetup.getData(tok_geom_);
0573   GlobalVector BField = bField_->inTesla(GlobalPoint(0, 0, 0));
0574   bfVal_ = BField.mag();
0575 
0576   trigger::TriggerEvent triggerEvent;
0577   edm::Handle<trigger::TriggerEvent> triggerEventHandle;
0578   iEvent.getByToken(tok_trigEvt_, triggerEventHandle);
0579   if (!triggerEventHandle.isValid()) {
0580     edm::LogWarning("IsoTrack") << "Error! Can't get the product hltTriggerSummaryAOD";
0581 
0582   } else {
0583     triggerEvent = *(triggerEventHandle.product());
0584   }
0585   const trigger::TriggerObjectCollection &TOC(triggerEvent.getObjects());
0586   /////////////////////////////TriggerResults
0587   edm::Handle<edm::TriggerResults> triggerResults;
0588   iEvent.getByToken(tok_trigRes_, triggerResults);
0589 
0590   edm::Handle<reco::TrackCollection> trkCollection;
0591   iEvent.getByToken(tok_genTrack_, trkCollection);
0592 
0593   iEvent.getByToken(tok_EB_, barrelRecHitsHandle_);
0594   iEvent.getByToken(tok_EE_, endcapRecHitsHandle_);
0595 
0596   iEvent.getByToken(tok_hbhe_, hbhe_);
0597 
0598   iEvent.getByToken(tok_recVtx_, recVtxs_);
0599   iEvent.getByToken(tok_bs_, beamSpotH_);
0600   if (!recVtxs_->empty() && !((*recVtxs_)[0].isFake())) {
0601     leadPV_ = math::XYZPoint((*recVtxs_)[0].x(), (*recVtxs_)[0].y(), (*recVtxs_)[0].z());
0602   } else if (beamSpotH_.isValid()) {
0603     leadPV_ = beamSpotH_->position();
0604   }
0605 
0606   if ((verbosity_ / 100) % 10 > 0) {
0607     edm::LogVerbatim("IsoTrack") << "Primary Vertex " << leadPV_;
0608     if (beamSpotH_.isValid())
0609       edm::LogVerbatim("IsoTrack") << "Beam Spot " << beamSpotH_->position();
0610   }
0611   pixelTrackRefsHE_.clear();
0612   pixelTrackRefsHB_.clear();
0613   for (unsigned int iPix = 0; iPix < pixelTracksSources_.size(); iPix++) {
0614     edm::Handle<reco::TrackCollection> iPixCol;
0615     iEvent.getByToken(tok_pixtks_[iPix], iPixCol);
0616     if (iPixCol.isValid()) {
0617       for (reco::TrackCollection::const_iterator pit = iPixCol->begin(); pit != iPixCol->end(); pit++) {
0618         if (iPix == 0)
0619           pixelTrackRefsHB_.push_back(reco::TrackRef(iPixCol, pit - iPixCol->begin()));
0620         pixelTrackRefsHE_.push_back(reco::TrackRef(iPixCol, pit - iPixCol->begin()));
0621       }
0622     }
0623   }
0624   if (doTiming_)
0625     getGoodTracks(iEvent, trkCollection);
0626 
0627   for (unsigned int ifilter = 0; ifilter < triggerEvent.sizeFilters(); ++ifilter) {
0628     std::string FilterNames[7] = {"hltL1sL1SingleJet68",
0629                                   "hltIsolPixelTrackL2FilterHE",
0630                                   "ecalIsolPixelTrackFilterHE",
0631                                   "hltIsolPixelTrackL3FilterHE",
0632                                   "hltIsolPixelTrackL2FilterHB",
0633                                   "ecalIsolPixelTrackFilterHB",
0634                                   "hltIsolPixelTrackL3FilterHB"};
0635     std::string label = triggerEvent.filterTag(ifilter).label();
0636     for (int i = 0; i < 7; i++) {
0637       if (label == FilterNames[i])
0638         h_Filters->Fill(i);
0639     }
0640   }
0641   edm::InputTag lumiProducer("LumiProducer", "", "RECO");
0642   edm::Handle<LumiDetails> Lumid;
0643   iEvent.getLuminosityBlock().getByToken(tok_lumi_, Lumid);
0644   float mybxlumi = -1;
0645   if (Lumid.isValid())
0646     mybxlumi = Lumid->lumiValue(LumiDetails::kOCC1, iEvent.bunchCrossing()) * 6.37;
0647   if (verbosity_ % 10 > 0)
0648     edm::LogVerbatim("IsoTrack") << "RunNo " << RunNo << " EvtNo " << iEvent.id().event() << " Lumi "
0649                                  << iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing() << " mybxlumi "
0650                                  << mybxlumi;
0651   if (!triggerResults.isValid()) {
0652     edm::LogWarning("IsoTrack") << "Error! Can't get the product triggerResults";
0653     //      std::shared_ptr<cms::Exception> const & error = triggerResults.whyFailed();
0654     //      edm::LogWarning(error->category()) << error->what();
0655   } else {
0656     std::vector<std::string> modules;
0657     h_nHLT->Fill(triggerResults->size());
0658     const edm::TriggerNames &triggerNames = iEvent.triggerNames(*triggerResults);
0659 
0660     const std::vector<std::string> &triggerNames_ = triggerNames.triggerNames();
0661     if (verbosity_ % 10 > 1)
0662       edm::LogVerbatim("IsoTrack") << "number of HLTs " << triggerNames_.size();
0663     int hlt(-1), preL1(-1), preHLT(-1), prescale(-1);
0664     for (unsigned int i = 0; i < triggerResults->size(); i++) {
0665       unsigned int triggerindx = hltConfig.triggerIndex(triggerNames_[i]);
0666       const std::vector<std::string> &moduleLabels(hltConfig.moduleLabels(triggerindx));
0667 
0668       for (unsigned int in = 0; in < trigNames_.size(); ++in) {
0669         //    if (triggerNames_[i].find(trigNames_[in].c_str())!=std::string::npos || triggerNames_[i]==" ") {
0670         if (triggerNames_[i].find(trigNames_[in]) != std::string::npos) {
0671           if (verbosity_ % 10 > 0)
0672             edm::LogVerbatim("IsoTrack") << "trigger that i want " << triggerNames_[i] << " accept "
0673                                          << triggerResults->accept(i);
0674           hlt = triggerResults->accept(i);
0675           h_HLT->Fill(hlt);
0676           //        if (hlt>0 || triggerNames_[i]==" ") {
0677           if (hlt > 0) {
0678             edm::Handle<reco::IsolatedPixelTrackCandidateCollection> Pixcands;
0679             iEvent.getByToken(tok_pixtk_, Pixcands);
0680             edm::Handle<trigger::TriggerFilterObjectWithRefs> L1cands;
0681             iEvent.getByToken(tok_l1cand_, L1cands);
0682 
0683             const std::pair<int, int> prescales(hltPrescaleProvider_.prescaleValues(iEvent, iSetup, triggerNames_[i]));
0684             preL1 = prescales.first;
0685             preHLT = prescales.second;
0686             prescale = preL1 * preHLT;
0687             if (verbosity_ % 10 > 0)
0688               edm::LogVerbatim("IsoTrack")
0689                   << triggerNames_[i] << " accept " << hlt << " preL1 " << preL1 << " preHLT " << preHLT;
0690             for (int iv = 0; iv < 3; ++iv)
0691               vec_[iv].clear();
0692             if (trigList_.find(RunNo) != trigList_.end()) {
0693               trigList_[RunNo] += 1;
0694             } else {
0695               trigList_.insert(std::pair<unsigned int, unsigned int>(RunNo, 1));
0696               trigPreList_.insert(std::pair<unsigned int, std::pair<int, int>>(RunNo, prescales));
0697             }
0698             //loop over all trigger filters in event (i.e. filters passed)
0699             for (unsigned int ifilter = 0; ifilter < triggerEvent.sizeFilters(); ++ifilter) {
0700               std::vector<int> Keys;
0701               std::string label = triggerEvent.filterTag(ifilter).label();
0702               //loop over keys to objects passing this filter
0703               for (unsigned int imodule = 0; imodule < moduleLabels.size(); imodule++) {
0704                 if (label.find(moduleLabels[imodule]) != std::string::npos) {
0705                   if (verbosity_ % 10 > 0)
0706                     edm::LogVerbatim("IsoTrack") << "FILTERNAME " << label;
0707                   for (unsigned int ifiltrKey = 0; ifiltrKey < triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
0708                     Keys.push_back(triggerEvent.filterKeys(ifilter)[ifiltrKey]);
0709                     const trigger::TriggerObject &TO(TOC[Keys[ifiltrKey]]);
0710                     math::XYZTLorentzVector v4(TO.px(), TO.py(), TO.pz(), TO.energy());
0711                     if (label.find("L2Filter") != std::string::npos) {
0712                       vec_[1].push_back(v4);
0713                     } else if (label.find("L3Filter") != std::string::npos) {
0714                       vec_[2].push_back(v4);
0715                     } else {
0716                       vec_[0].push_back(v4);
0717                       h_L1ObjEnergy->Fill(TO.energy());
0718                     }
0719                     if (verbosity_ % 10 > 0)
0720                       edm::LogVerbatim("IsoTrack") << "key " << ifiltrKey << " : pt " << TO.pt() << " eta " << TO.eta()
0721                                                    << " phi " << TO.phi() << " mass " << TO.mass() << " Id " << TO.id();
0722                   }
0723                 }
0724               }
0725             }
0726             std::vector<reco::TrackCollection::const_iterator> goodTks;
0727             if (doL2L3_) {
0728               h_nL3Objs->Fill(vec_[2].size());
0729               studyTrigger(trkCollection, goodTks);
0730             } else {
0731               if (trkCollection.isValid()) {
0732                 reco::TrackCollection::const_iterator trkItr;
0733                 for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); trkItr++)
0734                   goodTks.push_back(trkItr);
0735               }
0736             }
0737             // Now study isolation etc
0738             if (doStudyIsol_)
0739               studyIsolation(trkCollection, goodTks);
0740             if (doTrkResTree_)
0741               StudyTrkEbyP(trkCollection);
0742 
0743             std::pair<double, double> etaphi = etaPhiTrigger();
0744             edm::Handle<reco::IsolatedPixelTrackCandidateCollection> L2cands;
0745             iEvent.getByToken(tok_l2cand_, L2cands);
0746             if (!L2cands.isValid()) {
0747               if (verbosity_ % 10 > 0)
0748                 edm::LogVerbatim("IsoTrack") << " trigCand is not valid ";
0749             } else {
0750               if (doMipCutTree_)
0751                 studyMipCut(trkCollection, L2cands);
0752             }
0753             if (!pixelTracksSources_.empty())
0754               if (doChgIsolTree_ && !pixelTrackRefsHE_.empty())
0755                 chgIsolation(etaphi.first, etaphi.second, trkCollection, iEvent);
0756           }
0757           break;
0758         }
0759       }
0760     }
0761     h_PreL1->Fill(preL1);
0762     h_PreHLT->Fill(preHLT);
0763     h_Pre->Fill(prescale);
0764     h_PreL1wt->Fill(preL1, mybxlumi);
0765     h_PreHLTwt->Fill(preHLT, mybxlumi);
0766 
0767     // check if trigger names in (new) config
0768     //      edm::LogVerbatim("IsoTrack") << "changed " << changed_;
0769     if (changed_) {
0770       changed_ = false;
0771       if ((verbosity_ / 10) % 10 > 1) {
0772         edm::LogVerbatim("IsoTrack") << "New trigger menu found !!!";
0773         const unsigned int n(hltConfig.size());
0774         for (unsigned itrig = 0; itrig < triggerNames_.size(); itrig++) {
0775           unsigned int triggerindx = hltConfig.triggerIndex(triggerNames_[itrig]);
0776           if (triggerindx >= n)
0777             edm::LogVerbatim("IsoTrack") << triggerNames_[itrig] << " " << triggerindx << " does not exist in"
0778                                          << " the current menu";
0779           else
0780             edm::LogVerbatim("IsoTrack") << triggerNames_[itrig] << " " << triggerindx << " exists";
0781         }
0782       }
0783     }
0784   }
0785   if (doTiming_)
0786     studyTiming(iEvent);
0787 }
0788 
0789 void IsoTrig::clearChgIsolnTreeVectors() {
0790   t_PixcandP->clear();
0791   t_PixcandPt->clear();
0792   t_PixcandEta->clear();
0793   t_PixcandPhi->clear();
0794   for (unsigned int i = 0; i < t_PixcandMaxP->size(); i++)
0795     t_PixcandMaxP[i].clear();
0796   t_PixcandMaxP->clear();
0797   t_PixTrkcandP->clear();
0798   t_PixTrkcandPt->clear();
0799   t_PixTrkcandEta->clear();
0800   t_PixTrkcandPhi->clear();
0801   t_PixTrkcandMaxP->clear();
0802   t_PixTrkcandselTk->clear();
0803 }
0804 
0805 void IsoTrig::clearMipCutTreeVectors() {
0806   t_NFcandP->clear();
0807   t_NFcandPt->clear();
0808   t_NFcandEta->clear();
0809   t_NFcandPhi->clear();
0810   t_NFcandEmip->clear();
0811   t_NFTrkcandP->clear();
0812   t_NFTrkcandPt->clear();
0813   t_NFTrkcandEta->clear();
0814   t_NFTrkcandPhi->clear();
0815   t_NFTrkcandEmip->clear();
0816   t_NFTrkMinDR->clear();
0817   t_NFTrkMinDP1->clear();
0818   t_NFTrkselTkFlag->clear();
0819   t_NFTrkqltyFlag->clear();
0820   t_NFTrkMissFlag->clear();
0821   t_NFTrkPVFlag->clear();
0822   t_NFTrkPropFlag->clear();
0823   t_NFTrkChgIsoFlag->clear();
0824   t_NFTrkNeuIsoFlag->clear();
0825   t_NFTrkMipFlag->clear();
0826   t_ECone->clear();
0827 }
0828 
0829 void IsoTrig::pushChgIsolnTreeVecs(math::XYZTLorentzVector &Pixcand,
0830                                    math::XYZTLorentzVector &Trkcand,
0831                                    std::vector<double> &PixMaxP,
0832                                    double &TrkMaxP,
0833                                    bool &selTk) {
0834   t_PixcandP->push_back(Pixcand.r());
0835   t_PixcandPt->push_back(Pixcand.pt());
0836   t_PixcandEta->push_back(Pixcand.eta());
0837   t_PixcandPhi->push_back(Pixcand.phi());
0838   t_PixcandMaxP->push_back(PixMaxP);
0839   t_PixTrkcandP->push_back(Trkcand.r());
0840   t_PixTrkcandPt->push_back(Trkcand.pt());
0841   t_PixTrkcandEta->push_back(Trkcand.eta());
0842   t_PixTrkcandPhi->push_back(Trkcand.phi());
0843   t_PixTrkcandMaxP->push_back(TrkMaxP);
0844   t_PixTrkcandselTk->push_back(selTk);
0845 }
0846 
0847 void IsoTrig::pushMipCutTreeVecs(math::XYZTLorentzVector &NFcand,
0848                                  math::XYZTLorentzVector &Trkcand,
0849                                  double &EmipNFcand,
0850                                  double &EmipTrkcand,
0851                                  double &mindR,
0852                                  double &mindP1,
0853                                  std::vector<bool> &Flags,
0854                                  double hCone) {
0855   t_NFcandP->push_back(NFcand.r());
0856   t_NFcandPt->push_back(NFcand.pt());
0857   t_NFcandEta->push_back(NFcand.eta());
0858   t_NFcandPhi->push_back(NFcand.phi());
0859   t_NFcandEmip->push_back(EmipNFcand);
0860   t_NFTrkcandP->push_back(Trkcand.r());
0861   t_NFTrkcandPt->push_back(Trkcand.pt());
0862   t_NFTrkcandEta->push_back(Trkcand.eta());
0863   t_NFTrkcandPhi->push_back(Trkcand.phi());
0864   t_NFTrkcandEmip->push_back(EmipTrkcand);
0865   t_NFTrkMinDR->push_back(mindR);
0866   t_NFTrkMinDP1->push_back(mindP1);
0867   t_NFTrkselTkFlag->push_back(Flags[0]);
0868   t_NFTrkqltyFlag->push_back(Flags[1]);
0869   t_NFTrkMissFlag->push_back(Flags[2]);
0870   t_NFTrkPVFlag->push_back(Flags[3]);
0871   t_NFTrkPropFlag->push_back(Flags[4]);
0872   t_NFTrkChgIsoFlag->push_back(Flags[5]);
0873   t_NFTrkNeuIsoFlag->push_back(Flags[6]);
0874   t_NFTrkMipFlag->push_back(Flags[7]);
0875   t_ECone->push_back(hCone);
0876 }
0877 
0878 void IsoTrig::beginJob() {
0879   char hname[100], htit[100];
0880   std::string levels[20] = {"L1",        "L2",          "L3",        "Reco",      "RecoMatch",     "RecoNoMatch",
0881                             "L2Match",   "L2NoMatch",   "L3Match",   "L3NoMatch", "HLTTrk",        "HLTGoodTrk",
0882                             "HLTIsoTrk", "HLTMip",      "HLTSelect", "nonHLTTrk", "nonHLTGoodTrk", "nonHLTIsoTrk",
0883                             "nonHLTMip", "nonHLTSelect"};
0884   if (doTiming_) {
0885     TimingTree_ = fs_->make<TTree>("TimingTree", "TimingTree");
0886     t_timeL2Prod = new std::vector<double>();
0887     t_nPixCand = new std::vector<int>();
0888     t_nPixSeed = new std::vector<int>();
0889     t_nGoodTk = new std::vector<int>();
0890 
0891     TimingTree_->Branch("t_timeL2Prod", "std::vector<double>", &t_timeL2Prod);
0892     TimingTree_->Branch("t_nPixCand", "std::vector<int>", &t_nPixCand);
0893     TimingTree_->Branch("t_nPixSeed", "std::vector<int>", &t_nPixSeed);
0894     TimingTree_->Branch("t_nGoodTk", "std::vector<int>", &t_nGoodTk);
0895   }
0896   if (doTrkResTree_) {
0897     TrkResTree_ = fs_->make<TTree>("TrkRestree", "TrkResTree");
0898     t_TrkhCone = new std::vector<double>();
0899     t_TrkP = new std::vector<double>();
0900     t_TrkselTkFlag = new std::vector<bool>();
0901     t_TrkqltyFlag = new std::vector<bool>();
0902     t_TrkMissFlag = new std::vector<bool>();
0903     t_TrkPVFlag = new std::vector<bool>();
0904     t_TrkNuIsolFlag = new std::vector<bool>();
0905 
0906     TrkResTree_->Branch("t_TrkhCone", "std::vector<double>", &t_TrkhCone);
0907     TrkResTree_->Branch("t_TrkP", "std::vector<double>", &t_TrkP);
0908     TrkResTree_->Branch("t_TrkselTkFlag", "std::vector<bool>", &t_TrkselTkFlag);
0909     TrkResTree_->Branch("t_TrkqltyFlag", "std::vector<bool>", &t_TrkqltyFlag);
0910     TrkResTree_->Branch("t_TrkMissFlag", "std::vector<bool>", &t_TrkMissFlag);
0911     TrkResTree_->Branch("t_TrkPVFlag", "std::vector<bool>", &t_TrkPVFlag);
0912     TrkResTree_->Branch("t_TrkNuIsolFlag", "std::vector<bool>", &t_TrkNuIsolFlag);
0913   }
0914   if (doChgIsolTree_) {
0915     ChgIsolnTree_ = fs_->make<TTree>("ChgIsolnTree", "ChgIsolnTree");
0916     t_PixcandP = new std::vector<double>();
0917     t_PixcandPt = new std::vector<double>();
0918     t_PixcandEta = new std::vector<double>();
0919     t_PixcandPhi = new std::vector<double>();
0920     t_PixcandMaxP = new std::vector<std::vector<double>>();
0921     t_PixTrkcandP = new std::vector<double>();
0922     t_PixTrkcandPt = new std::vector<double>();
0923     t_PixTrkcandEta = new std::vector<double>();
0924     t_PixTrkcandPhi = new std::vector<double>();
0925     t_PixTrkcandMaxP = new std::vector<double>();
0926     t_PixTrkcandselTk = new std::vector<bool>();
0927 
0928     ChgIsolnTree_->Branch("t_PixcandP", "std::vector<double>", &t_PixcandP);
0929     ChgIsolnTree_->Branch("t_PixcandPt", "std::vector<double>", &t_PixcandPt);
0930     ChgIsolnTree_->Branch("t_PixcandEta", "std::vector<double>", &t_PixcandEta);
0931     ChgIsolnTree_->Branch("t_PixcandPhi", "std::vector<double>", &t_PixcandPhi);
0932     ChgIsolnTree_->Branch("t_PixcandMaxP", "std::vector<std::vector<double> >", &t_PixcandMaxP);
0933     ChgIsolnTree_->Branch("t_PixTrkcandP", "std::vector<double>", &t_PixTrkcandP);
0934     ChgIsolnTree_->Branch("t_PixTrkcandPt", "std::vector<double>", &t_PixTrkcandPt);
0935     ChgIsolnTree_->Branch("t_PixTrkcandEta", "std::vector<double>", &t_PixTrkcandEta);
0936     ChgIsolnTree_->Branch("t_PixTrkcandPhi", "std::vector<double>", &t_PixTrkcandPhi);
0937     ChgIsolnTree_->Branch("t_PixTrkcandMaxP", "std::vector<double>", &t_PixTrkcandMaxP);
0938     ChgIsolnTree_->Branch("t_PixTrkcandselTk", "std::vector<bool>", &t_PixTrkcandselTk);
0939   }
0940   if (doMipCutTree_) {
0941     MipCutTree_ = fs_->make<TTree>("MipCutTree", "MipCutTree");
0942     t_NFcandP = new std::vector<double>();
0943     t_NFcandPt = new std::vector<double>();
0944     t_NFcandEta = new std::vector<double>();
0945     t_NFcandPhi = new std::vector<double>();
0946     t_NFcandEmip = new std::vector<double>();
0947     t_NFTrkcandP = new std::vector<double>();
0948     t_NFTrkcandPt = new std::vector<double>();
0949     t_NFTrkcandEta = new std::vector<double>();
0950     t_NFTrkcandPhi = new std::vector<double>();
0951     t_NFTrkcandEmip = new std::vector<double>();
0952     t_NFTrkMinDR = new std::vector<double>();
0953     t_NFTrkMinDP1 = new std::vector<double>();
0954     t_NFTrkselTkFlag = new std::vector<bool>();
0955     t_NFTrkqltyFlag = new std::vector<bool>();
0956     t_NFTrkMissFlag = new std::vector<bool>();
0957     t_NFTrkPVFlag = new std::vector<bool>();
0958     t_NFTrkPropFlag = new std::vector<bool>();
0959     t_NFTrkChgIsoFlag = new std::vector<bool>();
0960     t_NFTrkNeuIsoFlag = new std::vector<bool>();
0961     t_NFTrkMipFlag = new std::vector<bool>();
0962     t_ECone = new std::vector<double>();
0963 
0964     MipCutTree_->Branch("t_NFcandP", "std::vector<double>", &t_NFcandP);
0965     MipCutTree_->Branch("t_NFcandPt", "std::vector<double>", &t_NFcandPt);
0966     MipCutTree_->Branch("t_NFcandEta", "std::vector<double>", &t_NFcandEta);
0967     MipCutTree_->Branch("t_NFcandPhi", "std::vector<double>", &t_NFcandPhi);
0968     MipCutTree_->Branch("t_NFcandEmip", "std::vector<double>", &t_NFcandEmip);
0969     MipCutTree_->Branch("t_NFTrkcandP", "std::vector<double>", &t_NFTrkcandP);
0970     MipCutTree_->Branch("t_NFTrkcandPt", "std::vector<double>", &t_NFTrkcandPt);
0971     MipCutTree_->Branch("t_NFTrkcandEta", "std::vector<double>", &t_NFTrkcandEta);
0972     MipCutTree_->Branch("t_NFTrkcandPhi", "std::vector<double>", &t_NFTrkcandPhi);
0973     MipCutTree_->Branch("t_NFTrkcandEmip", "std::vector<double>", &t_NFTrkcandEmip);
0974     MipCutTree_->Branch("t_NFTrkMinDR", "std::vector<double>", &t_NFTrkMinDR);
0975     MipCutTree_->Branch("t_NFTrkMinDP1", "std::vector<double>", &t_NFTrkMinDP1);
0976     MipCutTree_->Branch("t_NFTrkselTkFlag", "std::vector<bool>", &t_NFTrkselTkFlag);
0977     MipCutTree_->Branch("t_NFTrkqltyFlag", "std::vector<bool>", &t_NFTrkqltyFlag);
0978     MipCutTree_->Branch("t_NFTrkMissFlag", "std::vector<bool>", &t_NFTrkMissFlag);
0979     MipCutTree_->Branch("t_NFTrkPVFlag", "std::vector<bool>", &t_NFTrkPVFlag);
0980     MipCutTree_->Branch("t_NFTrkPropFlag", "std::vector<bool>", &t_NFTrkPropFlag);
0981     MipCutTree_->Branch("t_NFTrkChgIsoFlag", "std::vector<bool>", &t_NFTrkChgIsoFlag);
0982     MipCutTree_->Branch("t_NFTrkNeuIsoFlag", "std::vector<bool>", &t_NFTrkNeuIsoFlag);
0983     MipCutTree_->Branch("t_NFTrkMipFlag", "std::vector<bool>", &t_NFTrkMipFlag);
0984     MipCutTree_->Branch("t_ECone", "std::vector<double>", &t_ECone);
0985   }
0986   h_Filters = fs_->make<TH1I>("h_Filters", "Filter Accepts", 10, 0, 10);
0987   std::string FilterNames[7] = {"hltL1sL1SingleJet68",
0988                                 "hltIsolPixelTrackL2FilterHE",
0989                                 "ecalIsolPixelTrackFilterHE",
0990                                 "hltIsolPixelTrackL3FilterHE",
0991                                 "hltIsolPixelTrackL2FilterHB",
0992                                 "ecalIsolPixelTrackFilterHB",
0993                                 "hltIsolPixelTrackL3FilterHB"};
0994   for (int i = 0; i < 7; i++) {
0995     h_Filters->GetXaxis()->SetBinLabel(i + 1, FilterNames[i].c_str());
0996   }
0997 
0998   h_nHLT = fs_->make<TH1I>("h_nHLT", "Size of trigger Names", 1000, 1, 1000);
0999   h_HLT = fs_->make<TH1I>("h_HLT", "HLT accept", 3, -1, 2);
1000   h_PreL1 = fs_->make<TH1I>("h_PreL1", "L1 Prescale", 500, 0, 500);
1001   h_PreHLT = fs_->make<TH1I>("h_PreHLT", "HLT Prescale", 50, 0, 50);
1002   h_Pre = fs_->make<TH1I>("h_Pre", "Prescale", 3000, 0, 3000);
1003 
1004   h_PreL1wt = fs_->make<TH1D>("h_PreL1wt", "Weighted L1 Prescale", 500, 0, 500);
1005   h_PreHLTwt = fs_->make<TH1D>("h_PreHLTwt", "Weighted HLT Prescale", 500, 0, 100);
1006   h_L1ObjEnergy = fs_->make<TH1D>("h_L1ObjEnergy", "Energy of L1Object", 500, 0.0, 500.0);
1007 
1008   h_EnIn = fs_->make<TH1D>("h_EnInEcal", "EnergyIn Ecal", 200, 0.0, 20.0);
1009   h_EnOut = fs_->make<TH1D>("h_EnOutEcal", "EnergyOut Ecal", 200, 0.0, 20.0);
1010   h_MipEnMatch =
1011       fs_->make<TH2D>("h_MipEnMatch", "MipEn: HLT level vs Reco Level (Matched)", 200, 0.0, 20.0, 200, 0.0, 20.0);
1012   h_MipEnNoMatch = fs_->make<TH2D>(
1013       "h_MipEnNoMatch", "MipEn: HLT level vs Reco Level (No Match Found)", 200, 0.0, 20.0, 200, 0.0, 20.0);
1014 
1015   if (doL2L3_) {
1016     h_nL3Objs = fs_->make<TH1I>("h_nL3Objs", "Number of L3 objects", 10, 0, 10);
1017 
1018     std::string pairs[9] = {"L2L3",
1019                             "L2L3Match",
1020                             "L2L3NoMatch",
1021                             "L3Reco",
1022                             "L3RecoMatch",
1023                             "L3RecoNoMatch",
1024                             "NewFilterReco",
1025                             "NewFilterRecoMatch",
1026                             "NewFilterRecoNoMatch"};
1027     for (int ipair = 0; ipair < 9; ipair++) {
1028       sprintf(hname, "h_dEta%s", pairs[ipair].c_str());
1029       sprintf(htit, "#Delta#eta for %s", pairs[ipair].c_str());
1030       h_dEta[ipair] = fs_->make<TH1D>(hname, htit, 200, -10.0, 10.0);
1031       h_dEta[ipair]->GetXaxis()->SetTitle("d#eta");
1032 
1033       sprintf(hname, "h_dPhi%s", pairs[ipair].c_str());
1034       sprintf(htit, "#Delta#phi for %s", pairs[ipair].c_str());
1035       h_dPhi[ipair] = fs_->make<TH1D>(hname, htit, 140, -7.0, 7.0);
1036       h_dPhi[ipair]->GetXaxis()->SetTitle("d#phi");
1037 
1038       sprintf(hname, "h_dPt%s", pairs[ipair].c_str());
1039       sprintf(htit, "#Delta dp_{T} for %s objects", pairs[ipair].c_str());
1040       h_dPt[ipair] = fs_->make<TH1D>(hname, htit, 400, -200.0, 200.0);
1041       h_dPt[ipair]->GetXaxis()->SetTitle("dp_{T} (GeV)");
1042 
1043       sprintf(hname, "h_dP%s", pairs[ipair].c_str());
1044       sprintf(htit, "#Delta p for %s objects", pairs[ipair].c_str());
1045       h_dP[ipair] = fs_->make<TH1D>(hname, htit, 400, -200.0, 200.0);
1046       h_dP[ipair]->GetXaxis()->SetTitle("dP (GeV)");
1047 
1048       sprintf(hname, "h_dinvPt%s", pairs[ipair].c_str());
1049       sprintf(htit, "#Delta (1/p_{T}) for %s objects", pairs[ipair].c_str());
1050       h_dinvPt[ipair] = fs_->make<TH1D>(hname, htit, 500, -0.4, 0.1);
1051       h_dinvPt[ipair]->GetXaxis()->SetTitle("d(1/p_{T})");
1052       sprintf(hname, "h_mindR%s", pairs[ipair].c_str());
1053       sprintf(htit, "min(#Delta R) for %s objects", pairs[ipair].c_str());
1054       h_mindR[ipair] = fs_->make<TH1D>(hname, htit, 500, 0.0, 1.0);
1055       h_mindR[ipair]->GetXaxis()->SetTitle("dR");
1056     }
1057 
1058     for (int lvl = 0; lvl < 2; lvl++) {
1059       sprintf(hname, "h_dEtaL1%s", levels[lvl + 1].c_str());
1060       sprintf(htit, "#Delta#eta for L1 and %s objects", levels[lvl + 1].c_str());
1061       h_dEtaL1[lvl] = fs_->make<TH1D>(hname, htit, 400, -10.0, 10.0);
1062 
1063       sprintf(hname, "h_dPhiL1%s", levels[lvl + 1].c_str());
1064       sprintf(htit, "#Delta#phi for L1 and %s objects", levels[lvl + 1].c_str());
1065       h_dPhiL1[lvl] = fs_->make<TH1D>(hname, htit, 280, -7.0, 7.0);
1066 
1067       sprintf(hname, "h_dRL1%s", levels[lvl + 1].c_str());
1068       sprintf(htit, "#Delta R for L1 and %s objects", levels[lvl + 1].c_str());
1069       h_dRL1[lvl] = fs_->make<TH1D>(hname, htit, 100, 0.0, 10.0);
1070     }
1071   }
1072 
1073   int levmin = (doL2L3_ ? 0 : 10);
1074   for (int ilevel = levmin; ilevel < 20; ilevel++) {
1075     sprintf(hname, "h_p%s", levels[ilevel].c_str());
1076     sprintf(htit, "p for %s objects", levels[ilevel].c_str());
1077     h_p[ilevel] = fs_->make<TH1D>(hname, htit, 100, 0.0, 500.0);
1078     h_p[ilevel]->GetXaxis()->SetTitle("p (GeV)");
1079 
1080     sprintf(hname, "h_pt%s", levels[ilevel].c_str());
1081     sprintf(htit, "p_{T} for %s objects", levels[ilevel].c_str());
1082     h_pt[ilevel] = fs_->make<TH1D>(hname, htit, 100, 0.0, 500.0);
1083     h_pt[ilevel]->GetXaxis()->SetTitle("p_{T} (GeV)");
1084 
1085     sprintf(hname, "h_eta%s", levels[ilevel].c_str());
1086     sprintf(htit, "#eta for %s objects", levels[ilevel].c_str());
1087     h_eta[ilevel] = fs_->make<TH1D>(hname, htit, 100, -5.0, 5.0);
1088     h_eta[ilevel]->GetXaxis()->SetTitle("#eta");
1089 
1090     sprintf(hname, "h_phi%s", levels[ilevel].c_str());
1091     sprintf(htit, "#phi for %s objects", levels[ilevel].c_str());
1092     h_phi[ilevel] = fs_->make<TH1D>(hname, htit, 70, -3.5, 3.50);
1093     h_phi[ilevel]->GetXaxis()->SetTitle("#phi");
1094   }
1095 
1096   std::string cuts[2] = {"HLTMatched", "HLTNotMatched"};
1097   std::string cuts2[2] = {"All", "Away from L1"};
1098   for (int icut = 0; icut < 2; icut++) {
1099     sprintf(hname, "h_eMip%s", cuts[icut].c_str());
1100     sprintf(htit, "eMip for %s tracks", cuts[icut].c_str());
1101     h_eMip[icut] = fs_->make<TH1D>(hname, htit, 200, 0.0, 10.0);
1102     h_eMip[icut]->GetXaxis()->SetTitle("E_{Mip} (GeV)");
1103 
1104     sprintf(hname, "h_eMaxNearP%s", cuts[icut].c_str());
1105     sprintf(htit, "eMaxNearP for %s tracks", cuts[icut].c_str());
1106     h_eMaxNearP[icut] = fs_->make<TH1D>(hname, htit, 240, -2.0, 10.0);
1107     h_eMaxNearP[icut]->GetXaxis()->SetTitle("E_{MaxNearP} (GeV)");
1108 
1109     sprintf(hname, "h_eNeutIso%s", cuts[icut].c_str());
1110     sprintf(htit, "eNeutIso for %s ", cuts[icut].c_str());
1111     h_eNeutIso[icut] = fs_->make<TH1D>(hname, htit, 200, 0.0, 10.0);
1112     h_eNeutIso[icut]->GetXaxis()->SetTitle("E_{NeutIso} (GeV)");
1113 
1114     for (int kcut = 0; kcut < 2; ++kcut) {
1115       for (int lim = 0; lim < 5; ++lim) {
1116         sprintf(hname, "h_etaCalibTracks%sCut%dLim%d", cuts[icut].c_str(), kcut, lim);
1117         sprintf(htit,
1118                 "#eta for %s isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)",
1119                 cuts[icut].c_str(),
1120                 pLimits_[lim],
1121                 pLimits_[lim + 1],
1122                 cuts2[kcut].c_str());
1123         h_etaCalibTracks[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 60, -30.0, 30.0);
1124         h_etaCalibTracks[lim][icut][kcut]->GetXaxis()->SetTitle("i#eta");
1125 
1126         sprintf(hname, "h_etaMipTracks%sCut%dLim%d", cuts[icut].c_str(), kcut, lim);
1127         sprintf(htit,
1128                 "#eta for %s charge isolated MIP tracks (%4.1f < p < %5.1f Gev/c %s)",
1129                 cuts[icut].c_str(),
1130                 pLimits_[lim],
1131                 pLimits_[lim + 1],
1132                 cuts2[kcut].c_str());
1133         h_etaMipTracks[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 60, -30.0, 30.0);
1134         h_etaMipTracks[lim][icut][kcut]->GetXaxis()->SetTitle("i#eta");
1135       }
1136     }
1137   }
1138 
1139   std::string ecut1[3] = {"all", "HLTMatched", "HLTNotMatched"};
1140   std::string ecut2[2] = {"without", "with"};
1141   int etac[48] = {-1,  -2,  -3,  -4,  -5,  -6,  -7,  -8,  -9, -10, -11, -12, -13, -14, -15, -16,
1142                   -17, -18, -19, -20, -21, -22, -23, -24, 1,  2,   3,   4,   5,   6,   7,   8,
1143                   9,   10,  11,  12,  13,  14,  15,  16,  17, 18,  19,  20,  21,  22,  23,  24};
1144   for (int icut = 0; icut < 6; icut++) {
1145     //    int i1 = (icut>3 ? 1 : 0);
1146     int i1 = (icut > 2 ? 1 : 0);
1147     int i2 = icut - i1 * 3;
1148     for (int kcut = 0; kcut < 48; kcut++) {
1149       for (int lim = 0; lim < 5; ++lim) {
1150         sprintf(hname, "h_eta%dEnHcal%s%s%d", etac[kcut], ecut1[i2].c_str(), ecut2[i1].c_str(), lim);
1151         sprintf(htit,
1152                 "HCAL energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation",
1153                 etac[kcut],
1154                 ecut1[i2].c_str(),
1155                 pLimits_[lim],
1156                 pLimits_[lim + 1],
1157                 ecut2[i1].c_str());
1158         h_eHcal[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 750, 0.0, 150.0);
1159         h_eHcal[lim][icut][kcut]->GetXaxis()->SetTitle("Energy (GeV)");
1160         sprintf(hname, "h_eta%dEnCalo%s%s%d", etac[kcut], ecut1[i2].c_str(), ecut2[i1].c_str(), lim);
1161         sprintf(htit,
1162                 "Calorimter energy for #eta=%d for %s tracks (p=%4.1f:%5.1f Gev) %s neutral isolation",
1163                 etac[kcut],
1164                 ecut1[i2].c_str(),
1165                 pLimits_[lim],
1166                 pLimits_[lim + 1],
1167                 ecut2[i1].c_str());
1168         h_eCalo[lim][icut][kcut] = fs_->make<TH1D>(hname, htit, 750, 0.0, 150.0);
1169         h_eCalo[lim][icut][kcut]->GetXaxis()->SetTitle("Energy (GeV)");
1170       }
1171     }
1172   }
1173 }
1174 
1175 // ------------ method called once each job just after ending the event loop  ------------
1176 void IsoTrig::endJob() {
1177   unsigned int preL1, preHLT;
1178   std::map<unsigned int, unsigned int>::iterator itr;
1179   std::map<unsigned int, const std::pair<int, int>>::iterator itrPre;
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<TH1I>("h_PrevsRN", "PreScale Vs Run Number", n, minRunNo_, maxRunNo_);
1185   g_PreL1 = fs_->make<TH1I>("h_PreL1vsRN", "L1 PreScale Vs Run Number", n, minRunNo_, maxRunNo_);
1186   g_PreHLT = fs_->make<TH1I>("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 (itr = trigList_.begin(), itrPre = trigPreList_.begin(); itr != trigList_.end(); itr++, itrPre++) {
1190     preL1 = (itrPre->second).first;
1191     preHLT = (itrPre->second).second;
1192     edm::LogVerbatim("IsoTrack") << itr->first << " " << itr->second << " " << itrPre->first << " " << preL1 << " "
1193                                  << preHLT;
1194     g_Accepts->Fill(itr->first, itr->second);
1195     g_PreL1->Fill(itr->first, preL1);
1196     g_PreHLT->Fill(itr->first, preHLT);
1197     g_Pre->Fill(itr->first, 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);