IsolatedParticles/plugins/IsolatedTracksCone.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    IsolatedParticles
0004 // Class:      IsolatedTracksCone
0005 //
0006 /**\class IsolatedTracksCone IsolatedTracksCone.cc Calibration/IsolatedParticles/plugins/IsolatedTracksCone.cc
0007
0008  Description: Studies properties of isolated particles in the context of
0009               cone algorithm
0010
0011  Implementation:
0012      <Notes on implementation>
0013 */
0014 //
0015 // Original Author: Jim Hirschauer (adaptation of Seema Sharma's
0016 //                                  IsolatedTracksNew)
0017 //         Created:  Thu Nov  6 15:30:40 CST 2008
0018 //
0019 //
0020 // system include files
0021 #include <cmath>
0022 #include <map>
0023 #include <memory>
0024 #include <sstream>
0025 #include <string>
0026 #include <vector>
0027
0028 #include <Math/GenVector/VectorUtil.h>
0029
0030 // user include files
0031 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0032 #include "Calibration/IsolatedParticles/interface/FindCaloHitCone.h"
0033 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
0034 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0035 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0036 #include "Calibration/IsolatedParticles/interface/eCone.h"
0037 #include "Calibration/IsolatedParticles/interface/MatchingSimTrack.h"
0038 #include "Calibration/IsolatedParticles/interface/CaloSimInfo.h"
0039
0040 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0041
0042 #include "DataFormats/Common/interface/Ref.h"
0043 #include "DataFormats/Common/interface/Handle.h"
0044 #include "DataFormats/Candidate/interface/Candidate.h"
0045 #include "DataFormats/DetId/interface/DetId.h"
0046 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0047 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0048 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0049 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0050 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0051 // muons and tracks
0052 #include "DataFormats/TrackReco/interface/Track.h"
0053 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0054 #include "DataFormats/TrackReco/interface/HitPattern.h"
0055 // Vertices
0056 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0057 #include "DataFormats/VertexReco/interface/Vertex.h"
0058 //L1 objects
0059 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
0060 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
0061 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0062 #include "DataFormats/Common/interface/TriggerResults.h"
0063 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0064 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0065 #include "DataFormats/Math/interface/deltaPhi.h"
0066
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/EventSetup.h"
0071 #include "FWCore/Framework/interface/MakerMacros.h"
0072 #include "FWCore/Common/interface/TriggerNames.h"
0073 #include "FWCore/Framework/interface/LuminosityBlock.h"
0074 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0075 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0076
0077 //TFile Service
0078 #include "FWCore/ServiceRegistry/interface/Service.h"
0079 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0080
0081 // ecal / hcal
0082 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0083 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0084 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0085 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0086 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0087 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0088
0089 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0090 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0091
0092 // SimHit, SimTrack, ...
0093 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0094 #include "SimDataFormats/Track/interface/SimTrack.h"
0095 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0096 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0097 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0098
0099 // track associator
0100 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0101 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0102
0103 // tracker hit associator
0104 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0105 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
0106
0107 // root objects
0108 #include "TROOT.h"
0109 #include "TSystem.h"
0110 #include "TFile.h"
0111 #include "TH1F.h"
0112 #include "TH2F.h"
0113 #include "TProfile.h"
0114 #include "TDirectory.h"
0115 #include "TTree.h"
0116
0117 class IsolatedTracksCone : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0118 public:
0119   explicit IsolatedTracksCone(const edm::ParameterSet&);
0120   ~IsolatedTracksCone() override;
0121
0122   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0123
0124 private:
0125   void beginJob() override;
0126   void analyze(const edm::Event&, const edm::EventSetup&) override;
0127   void endJob() override;
0128
0129   void printTrack(const reco::Track* pTrack);
0130
0131   //   void BookHistograms();
0132   void buildTree();
0133   void clearTrackVectors();
0134
0135   static constexpr int nEtaBins_ = 4;
0136   static constexpr int nPBins_ = 21;
0137   std::array<double, nPBins_ + 1> genPartPBins_;
0138   std::array<double, nEtaBins_ + 1> genPartEtaBins;
0139
0140   const bool doMC_;
0141   const int myverbose_;
0142   const bool useJetTrigger_;
0143   const double drLeadJetVeto_, ptMinLeadJet_;
0144   const int debugTrks_;
0145   const bool printTrkHitPattern_;
0146
0147   const TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0148
0149   const edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_L1extTauJet_;
0150   const edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_L1extCenJet_;
0151   const edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_L1extFwdJet_;
0152
0153   const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0154   const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0155   const edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0156
0157   const edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0158   const edm::EDGetTokenT<edm::SimTrackContainer> tok_simTk_;
0159   const edm::EDGetTokenT<edm::SimVertexContainer> tok_simVtx_;
0160
0161   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloEB_;
0162   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloEE_;
0163   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloHH_;
0164   const edm::EDGetTokenT<edm::TriggerResults> tok_trigger_;
0165
0166   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0167   const edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_caloTopology_;
0168   const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_topo_;
0169   const edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_ecalChStatus_;
0170   const edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0171
0172   const double minTrackP_, maxTrackEta_, maxNearTrackP_;
0173   const int debugEcalSimInfo_;
0174   const bool applyEcalIsolation_;
0175
0176   // track associator to detector parameters
0177   TrackAssociatorParameters parameters_;
0178   std::unique_ptr<TrackDetectorAssociator> trackAssociator_;
0179
0180   TTree* ntp_;
0181
0182   TH1F* h_RawPt;
0183   TH1F* h_RawP;
0184   TH1F* h_RawEta;
0185   TH1F* h_RawPhi;
0186
0187   int nRawTRK;
0188   int nFailHighPurityQaul;
0189   int nFailPt;
0190   int nFailEta;
0191   int nMissEcal;
0192   int nMissHcal;
0193   int nEVT;
0194   int nEVT_failL1;
0195   int nTRK;
0196   double leadL1JetPT;
0197   double leadL1JetEta;
0198   double leadL1JetPhi;
0199
0200   std::unique_ptr<std::vector<std::vector<int> > > t_v_hnNearTRKs;
0201   std::unique_ptr<std::vector<std::vector<int> > > t_v_hnLayers_maxNearP;
0202   std::unique_ptr<std::vector<std::vector<int> > > t_v_htrkQual_maxNearP;
0203   std::unique_ptr<std::vector<std::vector<double> > > t_v_hmaxNearP_goodTrk;
0204   std::unique_ptr<std::vector<std::vector<double> > > t_v_hmaxNearP;
0205
0206   std::unique_ptr<std::vector<std::vector<int> > > t_v_cone_hnNearTRKs;
0207   std::unique_ptr<std::vector<std::vector<int> > > t_v_cone_hnLayers_maxNearP;
0208   std::unique_ptr<std::vector<std::vector<int> > > t_v_cone_htrkQual_maxNearP;
0209   std::unique_ptr<std::vector<std::vector<double> > > t_v_cone_hmaxNearP_goodTrk;
0210   std::unique_ptr<std::vector<std::vector<double> > > t_v_cone_hmaxNearP;
0211
0212   std::unique_ptr<std::vector<double> > t_trkNOuterHits;
0213   std::unique_ptr<std::vector<double> > t_trkNLayersCrossed;
0214   std::unique_ptr<std::vector<double> > t_dtFromLeadJet;
0215   std::unique_ptr<std::vector<double> > t_trkP;
0216   std::unique_ptr<std::vector<double> > t_simP;
0217   std::unique_ptr<std::vector<double> > t_trkPt;
0218   std::unique_ptr<std::vector<double> > t_trkEta;
0219   std::unique_ptr<std::vector<double> > t_trkPhi;
0220   std::unique_ptr<std::vector<double> > t_e3x3;
0221
0222   std::unique_ptr<std::vector<std::vector<double> > > t_v_eDR;
0223   std::unique_ptr<std::vector<std::vector<double> > > t_v_eMipDR;
0224
0225   std::unique_ptr<std::vector<double> > t_h3x3;
0226   std::unique_ptr<std::vector<double> > t_h5x5;
0227   std::unique_ptr<std::vector<double> > t_hsim3x3;
0228   std::unique_ptr<std::vector<double> > t_hsim5x5;
0229
0230   std::unique_ptr<std::vector<double> > t_nRH_h3x3;
0231   std::unique_ptr<std::vector<double> > t_nRH_h5x5;
0232
0233   std::unique_ptr<std::vector<double> > t_hsim3x3Matched;
0234   std::unique_ptr<std::vector<double> > t_hsim5x5Matched;
0235   std::unique_ptr<std::vector<double> > t_hsim3x3Rest;
0236   std::unique_ptr<std::vector<double> > t_hsim5x5Rest;
0237   std::unique_ptr<std::vector<double> > t_hsim3x3Photon;
0238   std::unique_ptr<std::vector<double> > t_hsim5x5Photon;
0239   std::unique_ptr<std::vector<double> > t_hsim3x3NeutHad;
0240   std::unique_ptr<std::vector<double> > t_hsim5x5NeutHad;
0241   std::unique_ptr<std::vector<double> > t_hsim3x3CharHad;
0242   std::unique_ptr<std::vector<double> > t_hsim5x5CharHad;
0243   std::unique_ptr<std::vector<double> > t_hsim3x3PdgMatched;
0244   std::unique_ptr<std::vector<double> > t_hsim5x5PdgMatched;
0245   std::unique_ptr<std::vector<double> > t_hsim3x3Total;
0246   std::unique_ptr<std::vector<double> > t_hsim5x5Total;
0247
0248   std::unique_ptr<std::vector<int> > t_hsim3x3NMatched;
0249   std::unique_ptr<std::vector<int> > t_hsim3x3NRest;
0250   std::unique_ptr<std::vector<int> > t_hsim3x3NPhoton;
0251   std::unique_ptr<std::vector<int> > t_hsim3x3NNeutHad;
0252   std::unique_ptr<std::vector<int> > t_hsim3x3NCharHad;
0253   std::unique_ptr<std::vector<int> > t_hsim3x3NTotal;
0254
0255   std::unique_ptr<std::vector<int> > t_hsim5x5NMatched;
0256   std::unique_ptr<std::vector<int> > t_hsim5x5NRest;
0257   std::unique_ptr<std::vector<int> > t_hsim5x5NPhoton;
0258   std::unique_ptr<std::vector<int> > t_hsim5x5NNeutHad;
0259   std::unique_ptr<std::vector<int> > t_hsim5x5NCharHad;
0260   std::unique_ptr<std::vector<int> > t_hsim5x5NTotal;
0261
0262   std::unique_ptr<std::vector<double> > t_distFromHotCell_h3x3;
0263   std::unique_ptr<std::vector<int> > t_ietaFromHotCell_h3x3;
0264   std::unique_ptr<std::vector<int> > t_iphiFromHotCell_h3x3;
0265   std::unique_ptr<std::vector<double> > t_distFromHotCell_h5x5;
0266   std::unique_ptr<std::vector<int> > t_ietaFromHotCell_h5x5;
0267   std::unique_ptr<std::vector<int> > t_iphiFromHotCell_h5x5;
0268
0269   std::unique_ptr<std::vector<double> > t_trkHcalEne;
0270   std::unique_ptr<std::vector<double> > t_trkEcalEne;
0271
0272   std::unique_ptr<std::vector<std::vector<double> > > t_v_hsimInfoConeMatched;
0273   std::unique_ptr<std::vector<std::vector<double> > > t_v_hsimInfoConeRest;
0274   std::unique_ptr<std::vector<std::vector<double> > > t_v_hsimInfoConePhoton;
0275   std::unique_ptr<std::vector<std::vector<double> > > t_v_hsimInfoConeNeutHad;
0276   std::unique_ptr<std::vector<std::vector<double> > > t_v_hsimInfoConeCharHad;
0277   std::unique_ptr<std::vector<std::vector<double> > > t_v_hsimInfoConePdgMatched;
0278   std::unique_ptr<std::vector<std::vector<double> > > t_v_hsimInfoConeTotal;
0279
0280   std::unique_ptr<std::vector<std::vector<int> > > t_v_hsimInfoConeNMatched;
0281   std::unique_ptr<std::vector<std::vector<int> > > t_v_hsimInfoConeNRest;
0282   std::unique_ptr<std::vector<std::vector<int> > > t_v_hsimInfoConeNPhoton;
0283   std::unique_ptr<std::vector<std::vector<int> > > t_v_hsimInfoConeNNeutHad;
0284   std::unique_ptr<std::vector<std::vector<int> > > t_v_hsimInfoConeNCharHad;
0285   std::unique_ptr<std::vector<std::vector<int> > > t_v_hsimInfoConeNTotal;
0286
0287   std::unique_ptr<std::vector<std::vector<double> > > t_v_hsimCone;
0288   std::unique_ptr<std::vector<std::vector<double> > > t_v_hCone;
0289   std::unique_ptr<std::vector<std::vector<int> > > t_v_nRecHitsCone;
0290   std::unique_ptr<std::vector<std::vector<int> > > t_v_nSimHitsCone;
0291
0292   std::unique_ptr<std::vector<std::vector<double> > > t_v_distFromHotCell;
0293   std::unique_ptr<std::vector<std::vector<int> > > t_v_ietaFromHotCell;
0294   std::unique_ptr<std::vector<std::vector<int> > > t_v_iphiFromHotCell;
0295
0296   std::unique_ptr<std::vector<std::vector<int> > > t_v_hlTriggers;
0297   std::unique_ptr<std::vector<int> > t_hltHB;
0298   std::unique_ptr<std::vector<int> > t_hltHE;
0299   std::unique_ptr<std::vector<int> > t_hltL1Jet15;
0300   std::unique_ptr<std::vector<int> > t_hltJet30;
0301   std::unique_ptr<std::vector<int> > t_hltJet50;
0302   std::unique_ptr<std::vector<int> > t_hltJet80;
0303   std::unique_ptr<std::vector<int> > t_hltJet110;
0304   std::unique_ptr<std::vector<int> > t_hltJet140;
0305   std::unique_ptr<std::vector<int> > t_hltJet180;
0306   std::unique_ptr<std::vector<int> > t_hltL1SingleEG5;
0307   std::unique_ptr<std::vector<int> > t_hltZeroBias;
0308   std::unique_ptr<std::vector<int> > t_hltMinBiasHcal;
0309   std::unique_ptr<std::vector<int> > t_hltMinBiasEcal;
0310   std::unique_ptr<std::vector<int> > t_hltMinBiasPixel;
0311   std::unique_ptr<std::vector<int> > t_hltSingleIsoTau30_Trk5;
0312   std::unique_ptr<std::vector<int> > t_hltDoubleLooseIsoTau15_Trk5;
0313
0314   std::unique_ptr<std::vector<std::vector<int> > > t_v_RH_h3x3_ieta;
0315   std::unique_ptr<std::vector<std::vector<int> > > t_v_RH_h3x3_iphi;
0316   std::unique_ptr<std::vector<std::vector<double> > > t_v_RH_h3x3_ene;
0317   std::unique_ptr<std::vector<std::vector<int> > > t_v_RH_h5x5_ieta;
0318   std::unique_ptr<std::vector<std::vector<int> > > t_v_RH_h5x5_iphi;
0319   std::unique_ptr<std::vector<std::vector<double> > > t_v_RH_h5x5_ene;
0320   std::unique_ptr<std::vector<std::vector<int> > > t_v_RH_r26_ieta;
0321   std::unique_ptr<std::vector<std::vector<int> > > t_v_RH_r26_iphi;
0322   std::unique_ptr<std::vector<std::vector<double> > > t_v_RH_r26_ene;
0323   std::unique_ptr<std::vector<std::vector<int> > > t_v_RH_r44_ieta;
0324   std::unique_ptr<std::vector<std::vector<int> > > t_v_RH_r44_iphi;
0325   std::unique_ptr<std::vector<std::vector<double> > > t_v_RH_r44_ene;
0326
0327   std::unique_ptr<std::vector<unsigned int> > t_irun;
0328   std::unique_ptr<std::vector<unsigned int> > t_ievt;
0329   std::unique_ptr<std::vector<unsigned int> > t_ilum;
0330 };
0331
0332 IsolatedTracksCone::IsolatedTracksCone(const edm::ParameterSet& iConfig)
0333     : doMC_(iConfig.getUntrackedParameter<bool>("doMC", false)),
0334       myverbose_(iConfig.getUntrackedParameter<int>("verbosity", 5)),
0335       useJetTrigger_(iConfig.getUntrackedParameter<bool>("useJetTrigger", false)),
0336       drLeadJetVeto_(iConfig.getUntrackedParameter<double>("drLeadJetVeto", 1.2)),
0337       ptMinLeadJet_(iConfig.getUntrackedParameter<double>("ptMinLeadJet", 15.0)),
0338       debugTrks_(iConfig.getUntrackedParameter<int>("debugTracks")),
0339       printTrkHitPattern_(iConfig.getUntrackedParameter<bool>("printTrkHitPattern")),
0340       trackerHitAssociatorConfig_(consumesCollector()),
0341       tok_L1extTauJet_(
0342           consumes<l1extra::L1JetParticleCollection>(iConfig.getParameter<edm::InputTag>("L1extraTauJetSource"))),
0343       tok_L1extCenJet_(
0344           consumes<l1extra::L1JetParticleCollection>(iConfig.getParameter<edm::InputTag>("L1extraCenJetSource"))),
0345       tok_L1extFwdJet_(
0346           consumes<l1extra::L1JetParticleCollection>(iConfig.getParameter<edm::InputTag>("L1extraFwdJetSource"))),
0347       tok_EB_(consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"))),
0348       tok_EE_(consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"))),
0349       tok_hbhe_(consumes<HBHERecHitCollection>(edm::InputTag("hbhereco"))),
0350       tok_genTrack_(consumes<reco::TrackCollection>(edm::InputTag("generalTracks"))),
0351       tok_simTk_(consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"))),
0352       tok_simVtx_(consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"))),
0353       tok_caloEB_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEB"))),
0354       tok_caloEE_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEE"))),
0355       tok_caloHH_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"))),
0356       tok_trigger_(consumes<edm::TriggerResults>(edm::InputTag("TriggerResults", "", "HLT"))),
0357       tok_geom_(esConsumes()),
0358       tok_caloTopology_(esConsumes()),
0359       tok_topo_(esConsumes()),
0360       tok_ecalChStatus_(esConsumes()),
0361       tok_sevlv_(esConsumes()),
0362       minTrackP_(iConfig.getUntrackedParameter<double>("minTrackP", 10.0)),
0363       maxTrackEta_(iConfig.getUntrackedParameter<double>("maxTrackEta", 5.0)),
0364       maxNearTrackP_(iConfig.getUntrackedParameter<double>("maxNearTrackP", 1.0)),
0365       debugEcalSimInfo_(iConfig.getUntrackedParameter<int>("debugEcalSimInfo")),
0366       applyEcalIsolation_(iConfig.getUntrackedParameter<bool>("applyEcalIsolation")) {
0367   //now do what ever initialization is needed
0368
0369   edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
0370   edm::ConsumesCollector iC = consumesCollector();
0371   parameters_.loadParameters(parameters, iC);
0372   trackAssociator_ = std::make_unique<TrackDetectorAssociator>();
0373   trackAssociator_->useDefaultPropagator();
0374
0375   if (myverbose_ >= 0) {
0376     edm::LogVerbatim("IsoTrack") << "Parameters read from config file \n"
0377                                  << "myverbose_ " << myverbose_ << "\n"
0378                                  << "useJetTrigger_ " << useJetTrigger_ << "\n"
0379                                  << "drLeadJetVeto_ " << drLeadJetVeto_ << "\n"
0380                                  << "minTrackP_ " << minTrackP_ << "\n"
0381                                  << "maxTrackEta_ " << maxTrackEta_ << "\n"
0382                                  << "maxNearTrackP_ " << maxNearTrackP_;
0383   }
0384 }
0385
0386 IsolatedTracksCone::~IsolatedTracksCone() {}
0387
0388 void IsolatedTracksCone::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0389   edm::ParameterSetDescription desc;
0390   desc.addUntracked<bool>("doMC", false);
0391   desc.addUntracked<int>("verbosity", 1);
0392   desc.addUntracked<bool>("useJetTrigger", false);
0393   desc.addUntracked<double>("drLeadJetVeto", 1.2);
0394   desc.addUntracked<double>("ptMinLeadJet", 15.0);
0395   desc.addUntracked<int>("debugTracks", 0);
0396   desc.addUntracked<bool>("printTrkHitPattern", true);
0397   desc.addUntracked<double>("minTrackP", 1.0);
0398   desc.addUntracked<double>("maxTrackEta", 2.6);
0399   desc.addUntracked<double>("maxNearTrackP", 1.0);
0400   desc.addUntracked<bool>("debugEcalSimInfo", false);
0401   desc.addUntracked<bool>("applyEcalIsolation", true);
0402   desc.addUntracked<bool>("debugL1Info", false);
0403   desc.add<edm::InputTag>("L1extraTauJetSource", edm::InputTag("l1extraParticles", "Tau"));
0404   desc.add<edm::InputTag>("L1extraCenJetSource", edm::InputTag("l1extraParticles", "Central"));
0405   desc.add<edm::InputTag>("L1extraFwdJetSource", edm::InputTag("l1extraParticles", "Forward"));
0406   //The following parameters are used as default for track association
0407   edm::ParameterSetDescription desc_TrackAssoc;
0408   desc_TrackAssoc.add<double>("muonMaxDistanceSigmaX", 0.0);
0409   desc_TrackAssoc.add<double>("muonMaxDistanceSigmaY", 0.0);
0410   desc_TrackAssoc.add<edm::InputTag>("CSCSegmentCollectionLabel", edm::InputTag("cscSegments"));
0411   desc_TrackAssoc.add<double>("dRHcal", 9999.0);
0412   desc_TrackAssoc.add<double>("dREcal", 9999.0);
0413   desc_TrackAssoc.add<edm::InputTag>("CaloTowerCollectionLabel", edm::InputTag("towerMaker"));
0414   desc_TrackAssoc.add<bool>("useEcal", true);
0415   desc_TrackAssoc.add<double>("dREcalPreselection", 0.05);
0416   desc_TrackAssoc.add<edm::InputTag>("HORecHitCollectionLabel", edm::InputTag("horeco"));
0417   desc_TrackAssoc.add<double>("dRMuon", 9999.0);
0418   desc_TrackAssoc.add<std::string>("crossedEnergyType", "SinglePointAlongTrajectory");
0419   desc_TrackAssoc.add<double>("muonMaxDistanceX", 5.0);
0420   desc_TrackAssoc.add<double>("muonMaxDistanceY", 5.0);
0421   desc_TrackAssoc.add<bool>("useHO", false);
0422   desc_TrackAssoc.add<bool>("accountForTrajectoryChangeCalo", false);
0423   desc_TrackAssoc.add<edm::InputTag>("DTRecSegment4DCollectionLabel", edm::InputTag("dt4DSegments"));
0424   desc_TrackAssoc.add<edm::InputTag>("EERecHitCollectionLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0425   desc_TrackAssoc.add<double>("dRHcalPreselection", 0.2);
0426   desc_TrackAssoc.add<bool>("useMuon", false);
0427   desc_TrackAssoc.add<bool>("useCalo", true);
0428   desc_TrackAssoc.add<edm::InputTag>("EBRecHitCollectionLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0429   desc_TrackAssoc.add<double>("dRMuonPreselection", 0.2);
0430   desc_TrackAssoc.add<bool>("truthMatch", false);
0431   desc_TrackAssoc.add<edm::InputTag>("HBHERecHitCollectionLabel", edm::InputTag("hbhereco"));
0432   desc_TrackAssoc.add<bool>("useHcal", true);
0433   desc_TrackAssoc.add<bool>("usePreshower", false);
0434   desc_TrackAssoc.add<double>("dRPreshowerPreselection", 0.2);
0435   desc_TrackAssoc.add<double>("trajectoryUncertaintyTolerance", 1.0);
0436   desc.add<edm::ParameterSetDescription>("TrackAssociatorParameters", desc_TrackAssoc);
0437   descriptions.add("isolatedTracksCone", desc);
0438 }
0439
0440 void IsolatedTracksCone::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0441   unsigned int irun = (unsigned int)iEvent.id().run();
0442   unsigned int ilum = (unsigned int)iEvent.getLuminosityBlock().luminosityBlock();
0443   unsigned int ievt = (unsigned int)iEvent.id().event();
0444
0445   clearTrackVectors();
0446
0447   // check the L1 objects
0448   bool L1Pass = false;
0449   leadL1JetPT = -999, leadL1JetEta = -999, leadL1JetPhi = -999;
0450   if (!useJetTrigger_) {
0451     L1Pass = true;
0452   } else {
0453     edm::Handle<l1extra::L1JetParticleCollection> l1TauHandle;
0454     iEvent.getByToken(tok_L1extTauJet_, l1TauHandle);
0455     l1extra::L1JetParticleCollection::const_iterator itr;
0456     for (itr = l1TauHandle->begin(); itr != l1TauHandle->end(); ++itr) {
0457       if (itr->pt() > leadL1JetPT) {
0458         leadL1JetPT = itr->pt();
0459         leadL1JetEta = itr->eta();
0460         leadL1JetPhi = itr->phi();
0461       }
0462     }
0463     edm::Handle<l1extra::L1JetParticleCollection> l1CenJetHandle;
0464     iEvent.getByToken(tok_L1extCenJet_, l1CenJetHandle);
0465     for (itr = l1CenJetHandle->begin(); itr != l1CenJetHandle->end(); ++itr) {
0466       if (itr->pt() > leadL1JetPT) {
0467         leadL1JetPT = itr->pt();
0468         leadL1JetEta = itr->eta();
0469         leadL1JetPhi = itr->phi();
0470       }
0471     }
0472     edm::Handle<l1extra::L1JetParticleCollection> l1FwdJetHandle;
0473     iEvent.getByToken(tok_L1extFwdJet_, l1FwdJetHandle);
0474     for (itr = l1FwdJetHandle->begin(); itr != l1FwdJetHandle->end(); ++itr) {
0475       if (itr->pt() > leadL1JetPT) {
0476         leadL1JetPT = itr->pt();
0477         leadL1JetEta = itr->eta();
0478         leadL1JetPhi = itr->phi();
0479       }
0480     }
0481     if (leadL1JetPT > ptMinLeadJet_) {
0482       L1Pass = true;
0483     }
0484   }
0485
0486   ////////////////////////////
0487   // Break now if L1Pass is false
0488   ////////////////////////////
0489   if (!L1Pass) {
0490     nEVT_failL1++;
0491     //    edm::LogVerbatim("IsoTrack") << "L1Pass is false : " << L1Pass;
0492     //     return;
0493   }
0494
0495   ///////////////////////////////////////////////
0496   // Get the collection handles
0497   ///////////////////////////////////////////////
0498
0499   const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0500   const CaloSubdetectorGeometry* gEB = (geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel));
0501   const CaloSubdetectorGeometry* gEE = (geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap));
0502   const CaloSubdetectorGeometry* gHB = (geo->getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0503   const CaloSubdetectorGeometry* gHE = (geo->getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0504
0505   const CaloTopology* caloTopology = &iSetup.getData(tok_caloTopology_);
0506   const HcalTopology* theHBHETopology = &iSetup.getData(tok_topo_);
0507
0508   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
0509   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
0510   iEvent.getByToken(tok_EB_, barrelRecHitsHandle);
0511   iEvent.getByToken(tok_EE_, endcapRecHitsHandle);
0512
0513   // Retrieve the good/bad ECAL channels from the DB
0514   const EcalChannelStatus* theEcalChStatus = &iSetup.getData(tok_ecalChStatus_);
0515
0516   edm::Handle<HBHERecHitCollection> hbhe;
0517   iEvent.getByToken(tok_hbhe_, hbhe);
0518   const HBHERecHitCollection Hithbhe = *(hbhe.product());
0519
0520   edm::Handle<reco::TrackCollection> trkCollection;
0521   iEvent.getByToken(tok_genTrack_, trkCollection);
0522   reco::TrackCollection::const_iterator trkItr;
0523   if (debugTrks_ > 1) {
0524     edm::LogVerbatim("IsoTrack") << "Track Collection: ";
0525     edm::LogVerbatim("IsoTrack") << "Number of Tracks " << trkCollection->size();
0526   }
0527   std::string theTrackQuality = "highPurity";
0528   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0529
0530   //get Handles to SimTracks and SimHits
0531   edm::Handle<edm::SimTrackContainer> SimTk;
0532   if (doMC_)
0533     iEvent.getByToken(tok_simTk_, SimTk);
0534
0535   edm::Handle<edm::SimVertexContainer> SimVtx;
0536   if (doMC_)
0537     iEvent.getByToken(tok_simVtx_, SimVtx);
0538
0539   //get Handles to PCaloHitContainers of eb/ee/hbhe
0540   edm::Handle<edm::PCaloHitContainer> pcaloeb;
0541   if (doMC_)
0542     iEvent.getByToken(tok_caloEB_, pcaloeb);
0543
0544   edm::Handle<edm::PCaloHitContainer> pcaloee;
0545   if (doMC_)
0546     iEvent.getByToken(tok_caloEE_, pcaloee);
0547
0548   edm::Handle<edm::PCaloHitContainer> pcalohh;
0549   if (doMC_)
0550     iEvent.getByToken(tok_caloHH_, pcalohh);
0551
0552   /////////////////////////////////////////////////////////
0553   // Get HLT_IsoTrackHB/HE Information
0554   /////////////////////////////////////////////////////////
0555
0556   edm::Handle<edm::TriggerResults> triggerResults;
0557   iEvent.getByToken(tok_trigger_, triggerResults);
0558
0559   std::vector<int> v_hlTriggers;
0560   int hltHB(-99);
0561   int hltHE(-99);
0562   int hltL1Jet15(-99);
0563   int hltJet30(-99);
0564   int hltJet50(-99);
0565   int hltJet80(-99);
0566   int hltJet110(-99);
0567   int hltJet140(-99);
0568   int hltJet180(-99);
0569   int hltL1SingleEG5(-99);
0570   int hltZeroBias(-99);
0571   int hltMinBiasHcal(-99);
0572   int hltMinBiasEcal(-99);
0573   int hltMinBiasPixel(-99);
0574   int hltSingleIsoTau30_Trk5(-99);
0575   int hltDoubleLooseIsoTau15_Trk5(-99);
0576
0577   if (triggerResults.isValid()) {
0578     const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0579     // TriggerNames class  triggerNames.init(*triggerResults);
0580
0581     for (unsigned int i = 0; i < triggerResults->size(); i++) {
0582       //      edm::LogVerbatim("IsoTrack") << "triggerNames.triggerName(" << i << ") = " << triggerNames.triggerName(i);
0583       if (triggerNames.triggerName(i) == "HLT_IsoTrackHE_1E31")
0584         hltHE = triggerResults->accept(i);
0585       if (triggerNames.triggerName(i) == "HLT_IsoTrackHB_1E31")
0586         hltHB = triggerResults->accept(i);
0587       if (triggerNames.triggerName(i) == "HLT_L1Jet15")
0588         hltL1Jet15 = triggerResults->accept(i);
0589       if (triggerNames.triggerName(i) == "HLT_Jet30")
0590         hltJet30 = triggerResults->accept(i);
0591       if (triggerNames.triggerName(i) == "HLT_Jet50")
0592         hltJet50 = triggerResults->accept(i);
0593       if (triggerNames.triggerName(i) == "HLT_Jet80")
0594         hltJet80 = triggerResults->accept(i);
0595       if (triggerNames.triggerName(i) == "HLT_Jet110")
0596         hltJet110 = triggerResults->accept(i);
0597       if (triggerNames.triggerName(i) == "HLT_Jet140")
0598         hltJet140 = triggerResults->accept(i);
0599       if (triggerNames.triggerName(i) == "HLT_Jet180")
0600         hltJet180 = triggerResults->accept(i);
0601       if (triggerNames.triggerName(i) == "HLT_L1SingleEG5")
0602         hltL1SingleEG5 = triggerResults->accept(i);
0603       if (triggerNames.triggerName(i) == "HLT_ZeroBias")
0604         hltZeroBias = triggerResults->accept(i);
0605       if (triggerNames.triggerName(i) == "HLT_MinBiasHcal")
0606         hltMinBiasHcal = triggerResults->accept(i);
0607       if (triggerNames.triggerName(i) == "HLT_MinBiasEcal")
0608         hltMinBiasEcal = triggerResults->accept(i);
0609       if (triggerNames.triggerName(i) == "HLT_MinBiasPixel")
0610         hltMinBiasPixel = triggerResults->accept(i);
0611       if (triggerNames.triggerName(i) == "HLT_SingleIsoTau30_Trk5")
0612         hltSingleIsoTau30_Trk5 = triggerResults->accept(i);
0613       if (triggerNames.triggerName(i) == "HLT_DoubleLooseIsoTau15_Trk5")
0614         hltDoubleLooseIsoTau15_Trk5 = triggerResults->accept(i);
0615     }
0616   }
0617
0618   ////////////////////////////
0619   // Primary loop over tracks
0620   ////////////////////////////
0621   std::unique_ptr<TrackerHitAssociator> associate;
0622   if (doMC_)
0623     associate = std::make_unique<TrackerHitAssociator>(iEvent, trackerHitAssociatorConfig_);
0624
0625   nTRK = 0;
0626   nRawTRK = 0;
0627   nFailPt = 0;
0628   nFailEta = 0;
0629   nFailHighPurityQaul = 0;
0630   nMissEcal = 0;
0631   nMissHcal = 0;
0632
0633   const EcalSeverityLevelAlgo* sevlv = &iSetup.getData(tok_sevlv_);
0634   for (trkItr = trkCollection->begin(); trkItr != trkCollection->end(); ++trkItr) {
0635     nRawTRK++;
0636
0637     const reco::Track* pTrack = &(*trkItr);
0638
0639     /////////////////////////////////////////
0640     // Check for min Pt and max Eta P
0641     /////////////////////////////////////////
0642
0643     bool trkQual = pTrack->quality(trackQuality_);
0644     bool goodPt = pTrack->p() > minTrackP_;
0645     bool goodEta = std::abs(pTrack->momentum().eta()) < maxTrackEta_;
0646
0647     double eta1 = pTrack->momentum().eta();
0648     double phi1 = pTrack->momentum().phi();
0649     double pt1 = pTrack->pt();
0650     double p1 = pTrack->p();
0651
0652     if (!goodEta) {
0653       nFailEta++;
0654     }
0655     if (!goodPt) {
0656       nFailPt++;
0657     }
0658     if (!trkQual) {
0659       nFailHighPurityQaul++;
0660     }
0661
0662     h_RawPt->Fill(pt1);
0663     h_RawP->Fill(p1);
0664     h_RawEta->Fill(eta1);
0665     h_RawPhi->Fill(phi1);
0666
0667     if (!goodEta || !goodPt || !trkQual)
0668       continue;  // Skip to next track
0669
0670     ////////////////////////////////////////////
0671     // Find track trajectory
0672     ////////////////////////////////////////////
0673
0674     const FreeTrajectoryState fts1 =
0675         trackAssociator_->getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *pTrack);
0676
0677     TrackDetMatchInfo info1 = trackAssociator_->associate(iEvent, iSetup, fts1, parameters_);
0678
0679     ////////////////////////////////////////////
0680     // First confirm track makes it to Hcal
0681     ////////////////////////////////////////////
0682
0683     if (info1.trkGlobPosAtHcal.x() == 0 && info1.trkGlobPosAtHcal.y() == 0 && info1.trkGlobPosAtHcal.z() == 0) {
0684       nMissHcal++;
0685       continue;
0686     }
0687
0688     const GlobalPoint hpoint1(info1.trkGlobPosAtHcal.x(), info1.trkGlobPosAtHcal.y(), info1.trkGlobPosAtHcal.z());
0689
0690     ////////////////////////////
0691     // Get basic quantities
0692     ////////////////////////////
0693
0694     const reco::HitPattern& hitp = pTrack->hitPattern();
0695     int nLayersCrossed = hitp.trackerLayersWithMeasurement();
0696     int nOuterHits = hitp.stripTOBLayersWithMeasurement() + hitp.stripTECLayersWithMeasurement();
0697
0698     double simP = 0;
0699     if (doMC_) {
0700       edm::SimTrackContainer::const_iterator matchedSimTrk =
0701           spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false);
0702       simP = matchedSimTrk->momentum().P();
0703     }
0704     ////////////////////////////////////////////
0705     // Get Ecal Point
0706     ////////////////////////////////////////////
0707
0708     const GlobalPoint point1(info1.trkGlobPosAtEcal.x(), info1.trkGlobPosAtEcal.y(), info1.trkGlobPosAtEcal.z());
0709
0710     // Sanity check that track hits Ecal
0711
0712     if (info1.trkGlobPosAtEcal.x() == 0 && info1.trkGlobPosAtEcal.y() == 0 && info1.trkGlobPosAtEcal.z() == 0) {
0713       edm::LogVerbatim("IsoTrack") << "Track doesn't reach Ecal.";
0714       nMissEcal++;
0715       continue;
0716     }
0717
0718     // Get Track Momentum - make sure you have latest version of TrackDetMatchInfo
0719
0720     GlobalVector trackMomAtEcal = info1.trkMomAtEcal;
0721     GlobalVector trackMomAtHcal = info1.trkMomAtHcal;
0722
0723     /////////////////////////////////////////////////////////
0724     // If using Jet trigger, get distance from leading jet
0725     /////////////////////////////////////////////////////////
0726
0727     double drFromLeadJet = 999.0;
0728     if (useJetTrigger_) {
0729       double dphi = reco::deltaPhi(phi1, leadL1JetPhi);
0730       double deta = eta1 - leadL1JetEta;
0731       drFromLeadJet = sqrt(dphi * dphi + deta * deta);
0732     }
0733
0734     ///////////////////////////////////////////////////////
0735     // Define Arrays for sizes of Charge, Neut Iso Radii and
0736     // Clustering Cone Radius.
0737     ///////////////////////////////////////////////////////
0738
0739     const int a_size = 7;
0740     double a_coneR[a_size];
0741     double a_charIsoR[a_size];
0742     double a_neutIsoR[a_size];
0743
0744     a_coneR[0] = 17.49;  // = area of 2x2
0745     a_coneR[1] = 26.23;  // = area of 3x3
0746     a_coneR[2] = 30.61;
0747     a_coneR[3] = 34.98;  // = area of 4x4
0748     a_coneR[4] = 39.35;
0749     a_coneR[5] = 43.72;  // = area of 5x5
0750     a_coneR[6] = 52.46;  // = area of 6x6
0751
0752     for (int i = 0; i < a_size; i++) {
0753       a_charIsoR[i] = a_coneR[i] + 28.9;      // 28.9 gives 55.1 for 3x3 benchmark
0754       a_neutIsoR[i] = a_charIsoR[i] * 0.726;  // Ecal radius = 0.726*Hcal radius
0755     }
0756
0757     ///////////////////////////////////////////////////////
0758     // Do Neutral Iso in radius on Ecal surface.
0759     ///////////////////////////////////////////////////////
0760
0761     // NxN cluster
0762     double e3x3 = -999.0;
0763     double trkEcalEne = -999.0;
0764
0765     if (std::abs(point1.eta()) < 1.479) {
0766       const DetId isoCell = gEB->getClosestCell(point1);
0767       e3x3 = spr::eECALmatrix(
0768                  isoCell, barrelRecHitsHandle, endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology, sevlv, 1, 1)
0769                  .first;
0770       trkEcalEne = spr::eCaloSimInfo(iEvent, geo, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate);
0771     } else {
0772       const DetId isoCell = gEE->getClosestCell(point1);
0773       e3x3 = spr::eECALmatrix(
0774                  isoCell, barrelRecHitsHandle, endcapRecHitsHandle, *theEcalChStatus, geo, caloTopology, sevlv, 1, 1)
0775                  .first;
0776       trkEcalEne = spr::eCaloSimInfo(iEvent, geo, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate);
0777     }
0778
0779     // Cone cluster
0780
0781     // Set up array of cone sizes for MIP cut
0782     const int a_mip_size = 5;
0783     double a_mipR[a_mip_size];
0784     a_mipR[0] = 3.84;  // = area of 3x3 ecal
0785     a_mipR[1] = 14.0;
0786     a_mipR[2] = 19.0;
0787     a_mipR[3] = 24.0;
0788     a_mipR[4] = 9.0;  // = from standard analyzer
0789
0790     std::vector<double> v_eDR;
0791     for (int i = 0; i < a_size; i++) {
0792       int nRH_eDR = 0;
0793
0794       // Cone in ecal
0795       double eDR = spr::eCone_ecal(
0796           geo, barrelRecHitsHandle, endcapRecHitsHandle, hpoint1, point1, a_neutIsoR[i], trackMomAtEcal, nRH_eDR);
0797       v_eDR.push_back(eDR);
0798     }
0799
0800     std::vector<double> v_eMipDR;
0801     for (int i = 0; i < a_mip_size; i++) {
0802       int nRH_eMipDR = 0;
0803       double eMipDR = spr::eCone_ecal(
0804           geo, barrelRecHitsHandle, endcapRecHitsHandle, hpoint1, point1, a_mipR[i], trackMomAtEcal, nRH_eMipDR);
0805
0806       v_eMipDR.push_back(eMipDR);
0807     }
0808
0809     ////////////////////////////////////////////
0810     // Do charge isolation in radius at Hcal surface for 5 different
0811     // radii defined above in a_charIso
0812     ////////////////////////////////////////////
0813
0814     std::vector<double> v_hmaxNearP_goodTrk;
0815     std::vector<double> v_hmaxNearP;
0816     std::vector<int> v_hnNearTRKs;
0817     std::vector<int> v_hnLayers_maxNearP;
0818     std::vector<int> v_htrkQual_maxNearP;
0819
0820     std::vector<double> v_cone_hmaxNearP_goodTrk;
0821     std::vector<double> v_cone_hmaxNearP;
0822     std::vector<int> v_cone_hnNearTRKs;
0823     std::vector<int> v_cone_hnLayers_maxNearP;
0824     std::vector<int> v_cone_htrkQual_maxNearP;
0825
0826     for (int i = 0; i < a_size; i++) {
0827       double hmaxNearP = -999.0;
0828       int hnNearTRKs = 0;
0829       int hnLayers_maxNearP = 0;
0830       int htrkQual_maxNearP = -1;
0831       double hmaxNearP_goodTrk = -999.0;
0832
0833       double conehmaxNearP = -999.0;
0834       int conehnNearTRKs = 0;
0835       int conehnLayers_maxNearP = 0;
0836       int conehtrkQual_maxNearP = -1;
0837       double conehmaxNearP_goodTrk = -999.0;
0838
0839       conehmaxNearP = spr::coneChargeIsolation(iEvent,
0840                                                iSetup,
0841                                                trkItr,
0842                                                trkCollection,
0843                                                *trackAssociator_,
0844                                                parameters_,
0845                                                theTrackQuality,
0846                                                conehnNearTRKs,
0847                                                conehnLayers_maxNearP,
0848                                                conehtrkQual_maxNearP,
0849                                                conehmaxNearP_goodTrk,
0850                                                hpoint1,
0851                                                trackMomAtHcal,
0852                                                a_charIsoR[i]);
0853
0854       v_hmaxNearP_goodTrk.push_back(hmaxNearP_goodTrk);
0855       v_hmaxNearP.push_back(hmaxNearP);
0856       v_hnNearTRKs.push_back(hnNearTRKs);
0857       v_hnLayers_maxNearP.push_back(hnLayers_maxNearP);
0858       v_htrkQual_maxNearP.push_back(htrkQual_maxNearP);
0859
0860       v_cone_hmaxNearP_goodTrk.push_back(conehmaxNearP_goodTrk);
0861       v_cone_hmaxNearP.push_back(conehmaxNearP);
0862       v_cone_hnNearTRKs.push_back(conehnNearTRKs);
0863       v_cone_hnLayers_maxNearP.push_back(conehnLayers_maxNearP);
0864       v_cone_htrkQual_maxNearP.push_back(conehtrkQual_maxNearP);
0865     }
0866
0867     double h3x3 = -999.0, h5x5 = -999.0;
0868     double hsim3x3 = -999.0, hsim5x5 = -999.0, trkHcalEne = -999.0;
0869     std::map<std::string, double> hsimInfo3x3, hsimInfo5x5;
0870     double distFromHotCell_h3x3 = -99.;
0871     int ietaFromHotCell_h3x3 = -99;
0872     int iphiFromHotCell_h3x3 = -99;
0873     double distFromHotCell_h5x5 = -99.;
0874     int ietaFromHotCell_h5x5 = -99;
0875     int iphiFromHotCell_h5x5 = -99;
0876
0877     GlobalPoint gPosHotCell_h3x3(0., 0., 0.);
0878     GlobalPoint gPosHotCell_h5x5(0., 0., 0.);
0879
0880     int nRH_h3x3(0), nRH_h5x5(0);
0881
0882     // Hcal Energy Clustering
0883
0884     // Get closetcell for ietaFromHotCell and iphiFromHotCell
0885     DetId ClosestCell;
0886     if (std::abs(pTrack->eta()) < 1.392) {
0887       ClosestCell = gHB->getClosestCell(hpoint1);
0888     } else {
0889       ClosestCell = gHE->getClosestCell(hpoint1);
0890     }
0891     // Transform into HcalDetId so that I can get ieta, iphi later.
0892     HcalDetId ClosestCell_HcalDetId(ClosestCell.rawId());
0893
0894     // Using NxN Cluster
0895     std::vector<int> v_RH_h3x3_ieta;
0896     std::vector<int> v_RH_h3x3_iphi;
0897     std::vector<double> v_RH_h3x3_ene;
0898     std::vector<int> v_RH_h5x5_ieta;
0899     std::vector<int> v_RH_h5x5_iphi;
0900     std::vector<double> v_RH_h5x5_ene;
0901
0902     h3x3 = spr::eHCALmatrix(geo,
0903                             theHBHETopology,
0904                             ClosestCell,
0905                             hbhe,
0906                             1,
0907                             1,
0908                             nRH_h3x3,
0909                             v_RH_h3x3_ieta,
0910                             v_RH_h3x3_iphi,
0911                             v_RH_h3x3_ene,
0912                             gPosHotCell_h3x3);
0913     distFromHotCell_h3x3 = spr::getDistInPlaneTrackDir(hpoint1, trackMomAtHcal, gPosHotCell_h3x3);
0914
0915     h5x5 = spr::eHCALmatrix(geo,
0916                             theHBHETopology,
0917                             ClosestCell,
0918                             hbhe,
0919                             2,
0920                             2,
0921                             nRH_h5x5,
0922                             v_RH_h5x5_ieta,
0923                             v_RH_h5x5_iphi,
0924                             v_RH_h5x5_ene,
0925                             gPosHotCell_h5x5);
0926     distFromHotCell_h5x5 = spr::getDistInPlaneTrackDir(hpoint1, trackMomAtHcal, gPosHotCell_h5x5);
0927
0928     //     double heta = (double)hpoint1.eta();
0929     //     double hphi = (double)hpoint1.phi();
0930     std::vector<int> multiplicity3x3;
0931     std::vector<int> multiplicity5x5;
0932     if (doMC_) {
0933       hsim3x3 = spr::eHCALmatrix(theHBHETopology, ClosestCell, pcalohh, 1, 1);
0934       hsim5x5 = spr::eHCALmatrix(theHBHETopology, ClosestCell, pcalohh, 2, 2);
0935
0936       hsimInfo3x3 = spr::eHCALSimInfo(
0937           iEvent, theHBHETopology, ClosestCell, pcalohh, SimTk, SimVtx, pTrack, *associate, 1, 1, multiplicity3x3);
0938       hsimInfo5x5 = spr::eHCALSimInfo(
0939           iEvent, theHBHETopology, ClosestCell, pcalohh, SimTk, SimVtx, pTrack, *associate, 2, 2, multiplicity5x5);
0940
0941       // Get energy from all simhits in hcal associated with iso track
0942       trkHcalEne = spr::eCaloSimInfo(iEvent, geo, pcalohh, SimTk, SimVtx, pTrack, *associate);
0943     }
0944
0945     // Finally for cones of varying radii.
0946     std::vector<double> v_hsimInfoConeMatched;
0947     std::vector<double> v_hsimInfoConeRest;
0948     std::vector<double> v_hsimInfoConePhoton;
0949     std::vector<double> v_hsimInfoConeNeutHad;
0950     std::vector<double> v_hsimInfoConeCharHad;
0951     std::vector<double> v_hsimInfoConePdgMatched;
0952     std::vector<double> v_hsimInfoConeTotal;
0953
0954     std::vector<int> v_hsimInfoConeNMatched;
0955     std::vector<int> v_hsimInfoConeNTotal;
0956     std::vector<int> v_hsimInfoConeNNeutHad;
0957     std::vector<int> v_hsimInfoConeNCharHad;
0958     std::vector<int> v_hsimInfoConeNPhoton;
0959     std::vector<int> v_hsimInfoConeNRest;
0960
0961     std::vector<double> v_hsimCone;
0962     std::vector<double> v_hCone;
0963
0964     std::vector<int> v_nRecHitsCone;
0965     std::vector<int> v_nSimHitsCone;
0966
0967     std::vector<double> v_distFromHotCell;
0968     std::vector<int> v_ietaFromHotCell;
0969     std::vector<int> v_iphiFromHotCell;
0970     GlobalPoint gposHotCell(0., 0., 0.);
0971
0972     std::vector<int> v_RH_r26_ieta;
0973     std::vector<int> v_RH_r26_iphi;
0974     std::vector<double> v_RH_r26_ene;
0975     std::vector<int> v_RH_r44_ieta;
0976     std::vector<int> v_RH_r44_iphi;
0977     std::vector<double> v_RH_r44_ene;
0978
0979     for (int i = 0; i < a_size; i++) {
0980       std::map<std::string, double> hsimInfoCone;
0981       double hsimCone = -999.0, hCone = -999.0;
0982       double distFromHotCell = -99.0;
0983       int ietaFromHotCell = -99;
0984       int iphiFromHotCell = -99;
0985       int ietaHotCell = -99;
0986       int iphiHotCell = -99;
0987       int nRecHitsCone = -999;
0988       int nSimHitsCone = -999;
0989
0990       std::vector<int> multiplicityCone;
0991       std::vector<DetId> coneRecHitDetIds;
0992       if (doMC_)
0993         hsimCone = spr::eCone_hcal(geo, pcalohh, hpoint1, point1, a_coneR[i], trackMomAtHcal, nSimHitsCone);
0994
0995       // If needed, get ieta and iphi of rechits for cones of 23.25
0996       // and for hitmap for debugging
0997       bool makeHitmaps = false;
0998       if (a_coneR[i] == 26.23 && makeHitmaps) {
0999         hCone = spr::eCone_hcal(geo,
1000                                 hbhe,
1001                                 hpoint1,
1002                                 point1,
1003                                 a_coneR[i],
1004                                 trackMomAtHcal,
1005                                 nRecHitsCone,
1006                                 v_RH_r26_ieta,
1007                                 v_RH_r26_iphi,
1008                                 v_RH_r26_ene,
1009                                 coneRecHitDetIds,
1010                                 distFromHotCell,
1011                                 ietaHotCell,
1012                                 iphiHotCell,
1013                                 gposHotCell);
1014       } else if (a_coneR[i] == 43.72 && makeHitmaps) {
1015         hCone = spr::eCone_hcal(geo,
1016                                 hbhe,
1017                                 hpoint1,
1018                                 point1,
1019                                 a_coneR[i],
1020                                 trackMomAtHcal,
1021                                 nRecHitsCone,
1022                                 v_RH_r44_ieta,
1023                                 v_RH_r44_iphi,
1024                                 v_RH_r44_ene,
1025                                 coneRecHitDetIds,
1026                                 distFromHotCell,
1027                                 ietaHotCell,
1028                                 iphiHotCell,
1029                                 gposHotCell);
1030       } else {
1031         hCone = spr::eCone_hcal(geo,
1032                                 hbhe,
1033                                 hpoint1,
1034                                 point1,
1035                                 a_coneR[i],
1036                                 trackMomAtHcal,
1037                                 nRecHitsCone,
1038                                 coneRecHitDetIds,
1039                                 distFromHotCell,
1040                                 ietaHotCell,
1041                                 iphiHotCell,
1042                                 gposHotCell);
1043       }
1044
1045       if (ietaHotCell != 99) {
1046         ietaFromHotCell = ietaHotCell - ClosestCell_HcalDetId.ieta();
1047         iphiFromHotCell = iphiHotCell - ClosestCell_HcalDetId.iphi();
1048       }
1049
1050       // SimHits NOT matched to RecHits
1051       if (doMC_) {
1052         hsimInfoCone = spr::eHCALSimInfoCone(iEvent,
1053                                              pcalohh,
1054                                              SimTk,
1055                                              SimVtx,
1056                                              pTrack,
1057                                              *associate,
1058                                              geo,
1059                                              hpoint1,
1060                                              point1,
1061                                              a_coneR[i],
1062                                              trackMomAtHcal,
1063                                              multiplicityCone);
1064
1065         // SimHits matched to RecHits
1066         //       hsimInfoCone = spr::eHCALSimInfoCone(iEvent,pcalohh, SimTk, SimVtx,
1067         //                             pTrack, *associate,
1068         //                             geo, hpoint1, point1,
1069         //                             a_coneR[i], trackMomAtHcal,
1070         //                             multiplicityCone,
1071         //                             coneRecHitDetIds);
1072
1073         v_hsimInfoConeMatched.push_back(hsimInfoCone["eMatched"]);
1074         v_hsimInfoConeRest.push_back(hsimInfoCone["eRest"]);
1075         v_hsimInfoConePhoton.push_back(hsimInfoCone["eGamma"]);
1076         v_hsimInfoConeNeutHad.push_back(hsimInfoCone["eNeutralHad"]);
1077         v_hsimInfoConeCharHad.push_back(hsimInfoCone["eChargedHad"]);
1078         v_hsimInfoConePdgMatched.push_back(hsimInfoCone["pdgMatched"]);
1079         v_hsimInfoConeTotal.push_back(hsimInfoCone["eTotal"]);
1080
1081         v_hsimInfoConeNMatched.push_back(multiplicityCone.at(0));
1082
1083         v_hsimInfoConeNTotal.push_back(multiplicityCone.at(1));
1084         v_hsimInfoConeNNeutHad.push_back(multiplicityCone.at(2));
1085         v_hsimInfoConeNCharHad.push_back(multiplicityCone.at(3));
1086         v_hsimInfoConeNPhoton.push_back(multiplicityCone.at(4));
1087         v_hsimInfoConeNRest.push_back(multiplicityCone.at(5));
1088
1089         v_hsimCone.push_back(hsimCone);
1090         v_nSimHitsCone.push_back(nSimHitsCone);
1091       }
1092       v_hCone.push_back(hCone);
1093       v_nRecHitsCone.push_back(nRecHitsCone);
1094
1095       v_distFromHotCell.push_back(distFromHotCell);
1096       v_ietaFromHotCell.push_back(ietaFromHotCell);
1097       v_iphiFromHotCell.push_back(iphiFromHotCell);
1098     }
1099
1100     ////////////////////////////////////////////
1101     // Fill Vectors that go into root file
1102     ////////////////////////////////////////////
1103
1104     t_v_hnNearTRKs->push_back(v_hnNearTRKs);
1105     t_v_hnLayers_maxNearP->push_back(v_hnLayers_maxNearP);
1106     t_v_htrkQual_maxNearP->push_back(v_htrkQual_maxNearP);
1107     t_v_hmaxNearP_goodTrk->push_back(v_hmaxNearP_goodTrk);
1108     t_v_hmaxNearP->push_back(v_hmaxNearP);
1109
1110     t_v_cone_hnNearTRKs->push_back(v_cone_hnNearTRKs);
1111     t_v_cone_hnLayers_maxNearP->push_back(v_cone_hnLayers_maxNearP);
1112     t_v_cone_htrkQual_maxNearP->push_back(v_cone_htrkQual_maxNearP);
1113     t_v_cone_hmaxNearP_goodTrk->push_back(v_cone_hmaxNearP_goodTrk);
1114     t_v_cone_hmaxNearP->push_back(v_cone_hmaxNearP);
1115
1116     //    t_hScale            ->push_back(hScale                    );
1117     t_trkNOuterHits->push_back(nOuterHits);
1118     t_trkNLayersCrossed->push_back(nLayersCrossed);
1119     t_dtFromLeadJet->push_back(drFromLeadJet);
1120     t_trkP->push_back(p1);
1121     t_trkPt->push_back(pt1);
1122     t_trkEta->push_back(eta1);
1123     t_trkPhi->push_back(phi1);
1124
1125     t_e3x3->push_back(e3x3);
1126     t_v_eDR->push_back(v_eDR);
1127     t_v_eMipDR->push_back(v_eMipDR);
1128
1129     t_h3x3->push_back(h3x3);
1130     t_h5x5->push_back(h5x5);
1131     t_nRH_h3x3->push_back(nRH_h3x3);
1132     t_nRH_h5x5->push_back(nRH_h5x5);
1133
1134     t_v_RH_h3x3_ieta->push_back(v_RH_h3x3_ieta);
1135     t_v_RH_h3x3_iphi->push_back(v_RH_h3x3_iphi);
1136     t_v_RH_h3x3_ene->push_back(v_RH_h3x3_ene);
1137     t_v_RH_h5x5_ieta->push_back(v_RH_h5x5_ieta);
1138     t_v_RH_h5x5_iphi->push_back(v_RH_h5x5_iphi);
1139     t_v_RH_h5x5_ene->push_back(v_RH_h5x5_ene);
1140
1141     if (doMC_) {
1142       t_simP->push_back(simP);
1143       t_hsim3x3->push_back(hsim3x3);
1144       t_hsim5x5->push_back(hsim5x5);
1145
1146       t_hsim3x3Matched->push_back(hsimInfo3x3["eMatched"]);
1147       t_hsim5x5Matched->push_back(hsimInfo5x5["eMatched"]);
1148       t_hsim3x3Rest->push_back(hsimInfo3x3["eRest"]);
1149       t_hsim5x5Rest->push_back(hsimInfo5x5["eRest"]);
1150       t_hsim3x3Photon->push_back(hsimInfo3x3["eGamma"]);
1151       t_hsim5x5Photon->push_back(hsimInfo5x5["eGamma"]);
1152       t_hsim3x3NeutHad->push_back(hsimInfo3x3["eNeutralHad"]);
1153       t_hsim5x5NeutHad->push_back(hsimInfo5x5["eNeutralHad"]);
1154       t_hsim3x3CharHad->push_back(hsimInfo3x3["eChargedHad"]);
1155       t_hsim5x5CharHad->push_back(hsimInfo5x5["eChargedHad"]);
1156       t_hsim3x3Total->push_back(hsimInfo3x3["eTotal"]);
1157       t_hsim5x5Total->push_back(hsimInfo5x5["eTotal"]);
1158       t_hsim3x3PdgMatched->push_back(hsimInfo3x3["pdgMatched"]);
1159       t_hsim5x5PdgMatched->push_back(hsimInfo5x5["pdgMatched"]);
1160
1161       t_hsim3x3NMatched->push_back(multiplicity3x3.at(0));
1162       t_hsim3x3NTotal->push_back(multiplicity3x3.at(1));
1163       t_hsim3x3NNeutHad->push_back(multiplicity3x3.at(2));
1164       t_hsim3x3NCharHad->push_back(multiplicity3x3.at(3));
1165       t_hsim3x3NPhoton->push_back(multiplicity3x3.at(4));
1166       t_hsim3x3NRest->push_back(multiplicity3x3.at(5));
1167
1168       t_hsim5x5NMatched->push_back(multiplicity5x5.at(0));
1169       t_hsim5x5NTotal->push_back(multiplicity5x5.at(1));
1170       t_hsim5x5NNeutHad->push_back(multiplicity5x5.at(2));
1171       t_hsim5x5NCharHad->push_back(multiplicity5x5.at(3));
1172       t_hsim5x5NPhoton->push_back(multiplicity5x5.at(4));
1173       t_hsim5x5NRest->push_back(multiplicity5x5.at(5));
1174     }
1175
1176     t_distFromHotCell_h3x3->push_back(distFromHotCell_h3x3);
1177     t_ietaFromHotCell_h3x3->push_back(ietaFromHotCell_h3x3);
1178     t_iphiFromHotCell_h3x3->push_back(iphiFromHotCell_h3x3);
1179     t_distFromHotCell_h5x5->push_back(distFromHotCell_h5x5);
1180     t_ietaFromHotCell_h5x5->push_back(ietaFromHotCell_h5x5);
1181     t_iphiFromHotCell_h5x5->push_back(iphiFromHotCell_h5x5);
1182
1183     if (doMC_) {
1184       t_trkHcalEne->push_back(trkHcalEne);
1185       t_trkEcalEne->push_back(trkEcalEne);
1186
1187       t_v_hsimInfoConeMatched->push_back(v_hsimInfoConeMatched);
1188       t_v_hsimInfoConeRest->push_back(v_hsimInfoConeRest);
1189       t_v_hsimInfoConePhoton->push_back(v_hsimInfoConePhoton);
1190       t_v_hsimInfoConeNeutHad->push_back(v_hsimInfoConeNeutHad);
1191       t_v_hsimInfoConeCharHad->push_back(v_hsimInfoConeCharHad);
1192       t_v_hsimInfoConePdgMatched->push_back(v_hsimInfoConePdgMatched);
1193       t_v_hsimInfoConeTotal->push_back(v_hsimInfoConeTotal);
1194
1195       t_v_hsimInfoConeNMatched->push_back(v_hsimInfoConeNMatched);
1196       t_v_hsimInfoConeNTotal->push_back(v_hsimInfoConeNTotal);
1197       t_v_hsimInfoConeNNeutHad->push_back(v_hsimInfoConeNNeutHad);
1198       t_v_hsimInfoConeNCharHad->push_back(v_hsimInfoConeNCharHad);
1199       t_v_hsimInfoConeNPhoton->push_back(v_hsimInfoConeNPhoton);
1200       t_v_hsimInfoConeNRest->push_back(v_hsimInfoConeNRest);
1201
1202       t_v_hsimCone->push_back(v_hsimCone);
1203       t_v_hCone->push_back(v_hCone);
1204       t_v_nRecHitsCone->push_back(v_nRecHitsCone);
1205       t_v_nSimHitsCone->push_back(v_nSimHitsCone);
1206     }
1207
1208     t_v_distFromHotCell->push_back(v_distFromHotCell);
1209     t_v_ietaFromHotCell->push_back(v_ietaFromHotCell);
1210     t_v_iphiFromHotCell->push_back(v_iphiFromHotCell);
1211
1212     t_v_RH_r26_ieta->push_back(v_RH_r26_ieta);
1213     t_v_RH_r26_iphi->push_back(v_RH_r26_iphi);
1214     t_v_RH_r26_ene->push_back(v_RH_r26_ene);
1215     t_v_RH_r44_ieta->push_back(v_RH_r44_ieta);
1216     t_v_RH_r44_iphi->push_back(v_RH_r44_iphi);
1217     t_v_RH_r44_ene->push_back(v_RH_r44_ene);
1218
1219     t_v_hlTriggers->push_back(v_hlTriggers);
1220     t_hltHB->push_back(hltHB);
1221     t_hltHE->push_back(hltHE);
1222     t_hltL1Jet15->push_back(hltL1Jet15);
1223     t_hltJet30->push_back(hltJet30);
1224     t_hltJet50->push_back(hltJet50);
1225     t_hltJet80->push_back(hltJet80);
1226     t_hltJet110->push_back(hltJet110);
1227     t_hltJet140->push_back(hltJet140);
1228     t_hltJet180->push_back(hltJet180);
1229     t_hltL1SingleEG5->push_back(hltL1SingleEG5);
1230     t_hltZeroBias->push_back(hltZeroBias);
1231     t_hltMinBiasHcal->push_back(hltMinBiasHcal);
1232     t_hltMinBiasEcal->push_back(hltMinBiasEcal);
1233     t_hltMinBiasPixel->push_back(hltMinBiasPixel);
1234     t_hltSingleIsoTau30_Trk5->push_back(hltSingleIsoTau30_Trk5);
1235     t_hltDoubleLooseIsoTau15_Trk5->push_back(hltDoubleLooseIsoTau15_Trk5);
1236
1237     t_irun->push_back(irun);
1238     t_ievt->push_back(ievt);
1239     t_ilum->push_back(ilum);
1240
1241     nTRK++;
1242
1243   }  // Loop over track collection
1244
1245   //  edm::LogVerbatim("IsoTrack") << "nEVT= " << nEVT;
1246
1247   ntp_->Fill();
1248   nEVT++;
1249 }
1250
1251 void IsolatedTracksCone::beginJob() {
1252   //   hbScale = 120.0;
1253   //   heScale = 135.0;
1254   nEVT = 0;
1255   nEVT_failL1 = 0;
1256   nTRK = 0;
1257
1258   genPartPBins_ = {{0.0,  1.0,  2.0,  3.0,  4.0,  5.0,  6.0,  7.0,  8.0,  9.0,  10.0,
1259                     12.0, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 100.0}};
1260
1261   genPartEtaBins = {{0.0, 0.5, 1.1, 1.7, 2.0}};
1262
1263   t_v_hnNearTRKs = std::make_unique<std::vector<std::vector<int> > >();
1264   t_v_hnLayers_maxNearP = std::make_unique<std::vector<std::vector<int> > >();
1265   t_v_htrkQual_maxNearP = std::make_unique<std::vector<std::vector<int> > >();
1266   t_v_hmaxNearP_goodTrk = std::make_unique<std::vector<std::vector<double> > >();
1267   t_v_hmaxNearP = std::make_unique<std::vector<std::vector<double> > >();
1268
1269   t_v_cone_hnNearTRKs = std::make_unique<std::vector<std::vector<int> > >();
1270   t_v_cone_hnLayers_maxNearP = std::make_unique<std::vector<std::vector<int> > >();
1271   t_v_cone_htrkQual_maxNearP = std::make_unique<std::vector<std::vector<int> > >();
1272   t_v_cone_hmaxNearP_goodTrk = std::make_unique<std::vector<std::vector<double> > >();
1273   t_v_cone_hmaxNearP = std::make_unique<std::vector<std::vector<double> > >();
1274
1275   //  t_hScale             = std::make_unique<std::vector<double> >();
1276   t_trkNOuterHits = std::make_unique<std::vector<double> >();
1277   t_trkNLayersCrossed = std::make_unique<std::vector<double> >();
1278   t_dtFromLeadJet = std::make_unique<std::vector<double> >();
1279   t_trkP = std::make_unique<std::vector<double> >();
1280   t_trkPt = std::make_unique<std::vector<double> >();
1281   t_trkEta = std::make_unique<std::vector<double> >();
1282   t_trkPhi = std::make_unique<std::vector<double> >();
1283
1284   t_e3x3 = std::make_unique<std::vector<double> >();
1285   t_v_eDR = std::make_unique<std::vector<std::vector<double> > >();
1286   t_v_eMipDR = std::make_unique<std::vector<std::vector<double> > >();
1287
1288   t_h3x3 = std::make_unique<std::vector<double> >();
1289   t_h5x5 = std::make_unique<std::vector<double> >();
1290
1291   t_nRH_h3x3 = std::make_unique<std::vector<double> >();
1292   t_nRH_h5x5 = std::make_unique<std::vector<double> >();
1293
1294   if (doMC_) {
1295     t_simP = std::make_unique<std::vector<double> >();
1296     t_hsim3x3 = std::make_unique<std::vector<double> >();
1297     t_hsim5x5 = std::make_unique<std::vector<double> >();
1298
1299     t_hsim3x3Matched = std::make_unique<std::vector<double> >();
1300     t_hsim5x5Matched = std::make_unique<std::vector<double> >();
1301     t_hsim3x3Rest = std::make_unique<std::vector<double> >();
1302     t_hsim5x5Rest = std::make_unique<std::vector<double> >();
1303     t_hsim3x3Photon = std::make_unique<std::vector<double> >();
1304     t_hsim5x5Photon = std::make_unique<std::vector<double> >();
1305     t_hsim3x3NeutHad = std::make_unique<std::vector<double> >();
1306     t_hsim5x5NeutHad = std::make_unique<std::vector<double> >();
1307     t_hsim3x3CharHad = std::make_unique<std::vector<double> >();
1308     t_hsim5x5CharHad = std::make_unique<std::vector<double> >();
1309     t_hsim3x3PdgMatched = std::make_unique<std::vector<double> >();
1310     t_hsim5x5PdgMatched = std::make_unique<std::vector<double> >();
1311     t_hsim3x3Total = std::make_unique<std::vector<double> >();
1312     t_hsim5x5Total = std::make_unique<std::vector<double> >();
1313
1314     t_hsim3x3NMatched = std::make_unique<std::vector<int> >();
1315     t_hsim3x3NTotal = std::make_unique<std::vector<int> >();
1316     t_hsim3x3NNeutHad = std::make_unique<std::vector<int> >();
1317     t_hsim3x3NCharHad = std::make_unique<std::vector<int> >();
1318     t_hsim3x3NPhoton = std::make_unique<std::vector<int> >();
1319     t_hsim3x3NRest = std::make_unique<std::vector<int> >();
1320
1321     t_hsim5x5NMatched = std::make_unique<std::vector<int> >();
1322     t_hsim5x5NTotal = std::make_unique<std::vector<int> >();
1323     t_hsim5x5NNeutHad = std::make_unique<std::vector<int> >();
1324     t_hsim5x5NCharHad = std::make_unique<std::vector<int> >();
1325     t_hsim5x5NPhoton = std::make_unique<std::vector<int> >();
1326     t_hsim5x5NRest = std::make_unique<std::vector<int> >();
1327
1328     t_trkHcalEne = std::make_unique<std::vector<double> >();
1329     t_trkEcalEne = std::make_unique<std::vector<double> >();
1330   }
1331
1332   t_distFromHotCell_h3x3 = std::make_unique<std::vector<double> >();
1333   t_ietaFromHotCell_h3x3 = std::make_unique<std::vector<int> >();
1334   t_iphiFromHotCell_h3x3 = std::make_unique<std::vector<int> >();
1335   t_distFromHotCell_h5x5 = std::make_unique<std::vector<double> >();
1336   t_ietaFromHotCell_h5x5 = std::make_unique<std::vector<int> >();
1337   t_iphiFromHotCell_h5x5 = std::make_unique<std::vector<int> >();
1338
1339   if (doMC_) {
1340     t_v_hsimInfoConeMatched = std::make_unique<std::vector<std::vector<double> > >();
1341     t_v_hsimInfoConeRest = std::make_unique<std::vector<std::vector<double> > >();
1342     t_v_hsimInfoConePhoton = std::make_unique<std::vector<std::vector<double> > >();
1343     t_v_hsimInfoConeNeutHad = std::make_unique<std::vector<std::vector<double> > >();
1344     t_v_hsimInfoConeCharHad = std::make_unique<std::vector<std::vector<double> > >();
1345     t_v_hsimInfoConePdgMatched = std::make_unique<std::vector<std::vector<double> > >();
1346     t_v_hsimInfoConeTotal = std::make_unique<std::vector<std::vector<double> > >();
1347
1348     t_v_hsimInfoConeNMatched = std::make_unique<std::vector<std::vector<int> > >();
1349     t_v_hsimInfoConeNTotal = std::make_unique<std::vector<std::vector<int> > >();
1350     t_v_hsimInfoConeNNeutHad = std::make_unique<std::vector<std::vector<int> > >();
1351     t_v_hsimInfoConeNCharHad = std::make_unique<std::vector<std::vector<int> > >();
1352     t_v_hsimInfoConeNPhoton = std::make_unique<std::vector<std::vector<int> > >();
1353     t_v_hsimInfoConeNRest = std::make_unique<std::vector<std::vector<int> > >();
1354
1355     t_v_hsimCone = std::make_unique<std::vector<std::vector<double> > >();
1356   }
1357
1358   t_v_hCone = std::make_unique<std::vector<std::vector<double> > >();
1359   t_v_nRecHitsCone = std::make_unique<std::vector<std::vector<int> > >();
1360   t_v_nSimHitsCone = std::make_unique<std::vector<std::vector<int> > >();
1361
1362   t_v_distFromHotCell = std::make_unique<std::vector<std::vector<double> > >();
1363   t_v_ietaFromHotCell = std::make_unique<std::vector<std::vector<int> > >();
1364   t_v_iphiFromHotCell = std::make_unique<std::vector<std::vector<int> > >();
1365
1366   t_v_RH_h3x3_ieta = std::make_unique<std::vector<std::vector<int> > >();
1367   t_v_RH_h3x3_iphi = std::make_unique<std::vector<std::vector<int> > >();
1368   t_v_RH_h3x3_ene = std::make_unique<std::vector<std::vector<double> > >();
1369   t_v_RH_h5x5_ieta = std::make_unique<std::vector<std::vector<int> > >();
1370   t_v_RH_h5x5_iphi = std::make_unique<std::vector<std::vector<int> > >();
1371   t_v_RH_h5x5_ene = std::make_unique<std::vector<std::vector<double> > >();
1372   t_v_RH_r26_ieta = std::make_unique<std::vector<std::vector<int> > >();
1373   t_v_RH_r26_iphi = std::make_unique<std::vector<std::vector<int> > >();
1374   t_v_RH_r26_ene = std::make_unique<std::vector<std::vector<double> > >();
1375   t_v_RH_r44_ieta = std::make_unique<std::vector<std::vector<int> > >();
1376   t_v_RH_r44_iphi = std::make_unique<std::vector<std::vector<int> > >();
1377   t_v_RH_r44_ene = std::make_unique<std::vector<std::vector<double> > >();
1378
1379   t_v_hlTriggers = std::make_unique<std::vector<std::vector<int> > >();
1380
1381   t_hltHE = std::make_unique<std::vector<int> >();
1382   t_hltHB = std::make_unique<std::vector<int> >();
1383   t_hltL1Jet15 = std::make_unique<std::vector<int> >();
1384   t_hltJet30 = std::make_unique<std::vector<int> >();
1385   t_hltJet50 = std::make_unique<std::vector<int> >();
1386   t_hltJet80 = std::make_unique<std::vector<int> >();
1387   t_hltJet110 = std::make_unique<std::vector<int> >();
1388   t_hltJet140 = std::make_unique<std::vector<int> >();
1389   t_hltJet180 = std::make_unique<std::vector<int> >();
1390   t_hltL1SingleEG5 = std::make_unique<std::vector<int> >();
1391   t_hltZeroBias = std::make_unique<std::vector<int> >();
1392   t_hltMinBiasHcal = std::make_unique<std::vector<int> >();
1393   t_hltMinBiasEcal = std::make_unique<std::vector<int> >();
1394   t_hltMinBiasPixel = std::make_unique<std::vector<int> >();
1395   t_hltSingleIsoTau30_Trk5 = std::make_unique<std::vector<int> >();
1396   t_hltDoubleLooseIsoTau15_Trk5 = std::make_unique<std::vector<int> >();
1397
1398   t_irun = std::make_unique<std::vector<unsigned int> >();
1399   t_ievt = std::make_unique<std::vector<unsigned int> >();
1400   t_ilum = std::make_unique<std::vector<unsigned int> >();
1401
1402   buildTree();
1403 }
1404
1405 void IsolatedTracksCone::clearTrackVectors() {
1406   t_v_hnNearTRKs->clear();
1407   t_v_hnLayers_maxNearP->clear();
1408   t_v_htrkQual_maxNearP->clear();
1409   t_v_hmaxNearP_goodTrk->clear();
1410   t_v_hmaxNearP->clear();
1411
1412   t_v_cone_hnNearTRKs->clear();
1413   t_v_cone_hnLayers_maxNearP->clear();
1414   t_v_cone_htrkQual_maxNearP->clear();
1415   t_v_cone_hmaxNearP_goodTrk->clear();
1416   t_v_cone_hmaxNearP->clear();
1417
1418   //  t_hScale             ->clear();
1419   t_trkNOuterHits->clear();
1420   t_trkNLayersCrossed->clear();
1421   t_dtFromLeadJet->clear();
1422   t_trkP->clear();
1423   t_trkPt->clear();
1424   t_trkEta->clear();
1425   t_trkPhi->clear();
1426   t_e3x3->clear();
1427   t_v_eDR->clear();
1428   t_v_eMipDR->clear();
1429   t_h3x3->clear();
1430   t_h5x5->clear();
1431   t_nRH_h3x3->clear();
1432   t_nRH_h5x5->clear();
1433
1434   if (doMC_) {
1435     t_simP->clear();
1436     t_hsim3x3->clear();
1437     t_hsim5x5->clear();
1438     t_hsim3x3Matched->clear();
1439     t_hsim5x5Matched->clear();
1440     t_hsim3x3Rest->clear();
1441     t_hsim5x5Rest->clear();
1442     t_hsim3x3Photon->clear();
1443     t_hsim5x5Photon->clear();
1444     t_hsim3x3NeutHad->clear();
1445     t_hsim5x5NeutHad->clear();
1446     t_hsim3x3CharHad->clear();
1447     t_hsim5x5CharHad->clear();
1448     t_hsim3x3PdgMatched->clear();
1449     t_hsim5x5PdgMatched->clear();
1450     t_hsim3x3Total->clear();
1451     t_hsim5x5Total->clear();
1452
1453     t_hsim3x3NMatched->clear();
1454     t_hsim3x3NTotal->clear();
1455     t_hsim3x3NNeutHad->clear();
1456     t_hsim3x3NCharHad->clear();
1457     t_hsim3x3NPhoton->clear();
1458     t_hsim3x3NRest->clear();
1459
1460     t_hsim5x5NMatched->clear();
1461     t_hsim5x5NTotal->clear();
1462     t_hsim5x5NNeutHad->clear();
1463     t_hsim5x5NCharHad->clear();
1464     t_hsim5x5NPhoton->clear();
1465     t_hsim5x5NRest->clear();
1466
1467     t_trkHcalEne->clear();
1468     t_trkEcalEne->clear();
1469   }
1470
1471   t_distFromHotCell_h3x3->clear();
1472   t_ietaFromHotCell_h3x3->clear();
1473   t_iphiFromHotCell_h3x3->clear();
1474   t_distFromHotCell_h5x5->clear();
1475   t_ietaFromHotCell_h5x5->clear();
1476   t_iphiFromHotCell_h5x5->clear();
1477
1478   if (doMC_) {
1479     t_v_hsimInfoConeMatched->clear();
1480     t_v_hsimInfoConeRest->clear();
1481     t_v_hsimInfoConePhoton->clear();
1482     t_v_hsimInfoConeNeutHad->clear();
1483     t_v_hsimInfoConeCharHad->clear();
1484     t_v_hsimInfoConePdgMatched->clear();
1485     t_v_hsimInfoConeTotal->clear();
1486
1487     t_v_hsimInfoConeNMatched->clear();
1488     t_v_hsimInfoConeNRest->clear();
1489     t_v_hsimInfoConeNPhoton->clear();
1490     t_v_hsimInfoConeNNeutHad->clear();
1491     t_v_hsimInfoConeNCharHad->clear();
1492     t_v_hsimInfoConeNTotal->clear();
1493
1494     t_v_hsimCone->clear();
1495   }
1496
1497   t_v_hCone->clear();
1498   t_v_nRecHitsCone->clear();
1499   t_v_nSimHitsCone->clear();
1500
1501   t_v_distFromHotCell->clear();
1502   t_v_ietaFromHotCell->clear();
1503   t_v_iphiFromHotCell->clear();
1504
1505   t_v_RH_h3x3_ieta->clear();
1506   t_v_RH_h3x3_iphi->clear();
1507   t_v_RH_h3x3_ene->clear();
1508   t_v_RH_h5x5_ieta->clear();
1509   t_v_RH_h5x5_iphi->clear();
1510   t_v_RH_h5x5_ene->clear();
1511   t_v_RH_r26_ieta->clear();
1512   t_v_RH_r26_iphi->clear();
1513   t_v_RH_r26_ene->clear();
1514   t_v_RH_r44_ieta->clear();
1515   t_v_RH_r44_iphi->clear();
1516   t_v_RH_r44_ene->clear();
1517
1518   t_v_hlTriggers->clear();
1519   t_hltHB->clear();
1520   t_hltHE->clear();
1521   t_hltL1Jet15->clear();
1522   t_hltJet30->clear();
1523   t_hltJet50->clear();
1524   t_hltJet80->clear();
1525   t_hltJet110->clear();
1526   t_hltJet140->clear();
1527   t_hltJet180->clear();
1528   t_hltL1SingleEG5->clear();
1529   t_hltZeroBias->clear();
1530   t_hltMinBiasHcal->clear();
1531   t_hltMinBiasEcal->clear();
1532   t_hltMinBiasPixel->clear();
1533   t_hltSingleIsoTau30_Trk5->clear();
1534   t_hltDoubleLooseIsoTau15_Trk5->clear();
1535
1536   t_irun->clear();
1537   t_ievt->clear();
1538   t_ilum->clear();
1539 }
1540
1541 // ----- method called once each job just after ending the event loop ----
1542 void IsolatedTracksCone::endJob() {
1543   edm::LogVerbatim("IsoTrack") << "Number of Events Processed " << nEVT << " failed L1 " << nEVT_failL1;
1544 }
1545
1546 void IsolatedTracksCone::buildTree() {
1547   edm::Service<TFileService> fs;
1548   h_RawPt = fs->make<TH1F>("hRawPt", "hRawPt", 100, 0.0, 100.0);
1549   h_RawP = fs->make<TH1F>("hRawP", "hRawP", 100, 0.0, 100.0);
1550   h_RawEta = fs->make<TH1F>("hRawEta", "hRawEta", 15, 0.0, 3.0);
1551   h_RawPhi = fs->make<TH1F>("hRawPhi", "hRawPhi", 100, -3.2, 3.2);
1552
1553   ntp_ = fs->make<TTree>("ntp", "ntp");
1554
1555   // Counters
1556   ntp_->Branch("nEVT", &nEVT, "nEVT/I");
1557   ntp_->Branch("leadL1JetPT", &leadL1JetPT, "leadL1JetPT/D");
1558   ntp_->Branch("leadL1JetEta", &leadL1JetEta, "leadL1JetEta/D");
1559   ntp_->Branch("leadL1JetPhi", &leadL1JetPhi, "leadL1JetPhi/D");
1560   ntp_->Branch("nTRK", &nTRK, "nTRK/I");
1561   ntp_->Branch("nRawTRK", &nRawTRK, "nRawTRK/I");
1562   ntp_->Branch("nFailHighPurityQaul", &nFailHighPurityQaul, "nFailHighPurityQaul/I");
1563   ntp_->Branch("nFailPt", &nFailPt, "nFailPt/I");
1564   ntp_->Branch("nFailEta", &nFailEta, "nFailEta/I");
1565   ntp_->Branch("nMissEcal", &nMissEcal, "nMissEcal/I");
1566   ntp_->Branch("nMissHcal", &nMissHcal, "nMissHcal/I");
1567
1568   ntp_->Branch("hnNearTRKs", "std::vector<std::vector<int> >   ", &t_v_hnNearTRKs);
1569   ntp_->Branch("hnLayers_maxNearP", "std::vector<std::vector<int> >   ", &t_v_hnLayers_maxNearP);
1570   ntp_->Branch("htrkQual_maxNearP", "std::vector<std::vector<int> >   ", &t_v_htrkQual_maxNearP);
1571   ntp_->Branch("hmaxNearP_goodTrk", "std::vector<std::vector<double> >", &t_v_hmaxNearP_goodTrk);
1572   ntp_->Branch("hmaxNearP", "std::vector<std::vector<double> >", &t_v_hmaxNearP);
1573
1574   ntp_->Branch("cone_hnNearTRKs", "std::vector<std::vector<int> >   ", &t_v_cone_hnNearTRKs);
1575   ntp_->Branch("cone_hnLayers_maxNearP", "std::vector<std::vector<int> >   ", &t_v_cone_hnLayers_maxNearP);
1576   ntp_->Branch("cone_htrkQual_maxNearP", "std::vector<std::vector<int> >   ", &t_v_cone_htrkQual_maxNearP);
1577   ntp_->Branch("cone_hmaxNearP_goodTrk", "std::vector<std::vector<double> >", &t_v_cone_hmaxNearP_goodTrk);
1578   ntp_->Branch("cone_hmaxNearP", "std::vector<std::vector<double> >", &t_v_cone_hmaxNearP);
1579
1580   //ntp_->Branch("hScale"           , "std::vector<double>", &t_hScale          );
1581   ntp_->Branch("trkNOuterHits", "std::vector<double>", &t_trkNOuterHits);
1582   ntp_->Branch("trkNLayersCrossed", "std::vector<double>", &t_trkNLayersCrossed);
1583   ntp_->Branch("drFromLeadJet", "std::vector<double>", &t_dtFromLeadJet);
1584   ntp_->Branch("trkP", "std::vector<double>", &t_trkP);
1585   ntp_->Branch("trkPt", "std::vector<double>", &t_trkPt);
1586   ntp_->Branch("trkEta", "std::vector<double>", &t_trkEta);
1587   ntp_->Branch("trkPhi", "std::vector<double>", &t_trkPhi);
1588   ntp_->Branch("e3x3", "std::vector<double>", &t_e3x3);
1589
1590   ntp_->Branch("e3x3", "std::vector<double>", &t_e3x3);
1591   ntp_->Branch("v_eDR", "std::vector<std::vector<double> >", &t_v_eDR);
1592   ntp_->Branch("v_eMipDR", "std::vector<std::vector<double> >", &t_v_eMipDR);
1593
1594   ntp_->Branch("h3x3", "std::vector<double>", &t_h3x3);
1595   ntp_->Branch("h5x5", "std::vector<double>", &t_h5x5);
1596   ntp_->Branch("nRH_h3x3", "std::vector<double>", &t_nRH_h3x3);
1597   ntp_->Branch("nRH_h5x5", "std::vector<double>", &t_nRH_h5x5);
1598
1599   if (doMC_) {
1600     ntp_->Branch("simP", "std::vector<double>", &t_simP);
1601     ntp_->Branch("hsim3x3", "std::vector<double>", &t_hsim3x3);
1602     ntp_->Branch("hsim5x5", "std::vector<double>", &t_hsim5x5);
1603
1604     ntp_->Branch("hsim3x3Matched", "std::vector<double>", &t_hsim3x3Matched);
1605     ntp_->Branch("hsim5x5Matched", "std::vector<double>", &t_hsim5x5Matched);
1606     ntp_->Branch("hsim3x3Rest", "std::vector<double>", &t_hsim3x3Rest);
1607     ntp_->Branch("hsim5x5Rest", "std::vector<double>", &t_hsim5x5Rest);
1608     ntp_->Branch("hsim3x3Photon", "std::vector<double>", &t_hsim3x3Photon);
1609     ntp_->Branch("hsim5x5Photon", "std::vector<double>", &t_hsim5x5Photon);
1610     ntp_->Branch("hsim3x3NeutHad", "std::vector<double>", &t_hsim3x3NeutHad);
1611     ntp_->Branch("hsim5x5NeutHad", "std::vector<double>", &t_hsim5x5NeutHad);
1612     ntp_->Branch("hsim3x3CharHad", "std::vector<double>", &t_hsim3x3CharHad);
1613     ntp_->Branch("hsim5x5CharHad", "std::vector<double>", &t_hsim5x5CharHad);
1614     ntp_->Branch("hsim3x3PdgMatched", "std::vector<double>", &t_hsim3x3PdgMatched);
1615     ntp_->Branch("hsim5x5PdgMatched", "std::vector<double>", &t_hsim5x5PdgMatched);
1616     ntp_->Branch("hsim3x3Total", "std::vector<double>", &t_hsim3x3Total);
1617     ntp_->Branch("hsim5x5Total", "std::vector<double>", &t_hsim5x5Total);
1618
1619     ntp_->Branch("hsim3x3NMatched", "std::vector<int>", &t_hsim3x3NMatched);
1620     ntp_->Branch("hsim3x3NRest", "std::vector<int>", &t_hsim3x3NRest);
1621     ntp_->Branch("hsim3x3NPhoton", "std::vector<int>", &t_hsim3x3NPhoton);
1622     ntp_->Branch("hsim3x3NNeutHad", "std::vector<int>", &t_hsim3x3NNeutHad);
1623     ntp_->Branch("hsim3x3NCharHad", "std::vector<int>", &t_hsim3x3NCharHad);
1624     ntp_->Branch("hsim3x3NTotal", "std::vector<int>", &t_hsim3x3NTotal);
1625
1626     ntp_->Branch("hsim5x5NMatched", "std::vector<int>", &t_hsim5x5NMatched);
1627     ntp_->Branch("hsim5x5NRest", "std::vector<int>", &t_hsim5x5NRest);
1628     ntp_->Branch("hsim5x5NPhoton", "std::vector<int>", &t_hsim5x5NPhoton);
1629     ntp_->Branch("hsim5x5NNeutHad", "std::vector<int>", &t_hsim5x5NNeutHad);
1630     ntp_->Branch("hsim5x5NCharHad", "std::vector<int>", &t_hsim5x5NCharHad);
1631     ntp_->Branch("hsim5x5NTotal", "std::vector<int>", &t_hsim5x5NTotal);
1632
1633     ntp_->Branch("trkHcalEne", "std::vector<double>", &t_trkHcalEne);
1634     ntp_->Branch("trkEcalEne", "std::vector<double>", &t_trkEcalEne);
1635   }
1636
1637   ntp_->Branch("distFromHotCell_h3x3", "std::vector<double>", &t_distFromHotCell_h3x3);
1638   ntp_->Branch("ietaFromHotCell_h3x3", "std::vector<int>", &t_ietaFromHotCell_h3x3);
1639   ntp_->Branch("iphiFromHotCell_h3x3", "std::vector<int>", &t_iphiFromHotCell_h3x3);
1640   ntp_->Branch("distFromHotCell_h5x5", "std::vector<double>", &t_distFromHotCell_h5x5);
1641   ntp_->Branch("ietaFromHotCell_h5x5", "std::vector<int>", &t_ietaFromHotCell_h5x5);
1642   ntp_->Branch("iphiFromHotCell_h5x5", "std::vector<int>", &t_iphiFromHotCell_h5x5);
1643
1644   if (doMC_) {
1645     ntp_->Branch("v_hsimInfoConeMatched", "std::vector<std::vector<double> >", &t_v_hsimInfoConeMatched);
1646     ntp_->Branch("v_hsimInfoConeRest", "std::vector<std::vector<double> >", &t_v_hsimInfoConeRest);
1647     ntp_->Branch("v_hsimInfoConePhoton", "std::vector<std::vector<double> >", &t_v_hsimInfoConePhoton);
1648     ntp_->Branch("v_hsimInfoConeNeutHad", "std::vector<std::vector<double> >", &t_v_hsimInfoConeNeutHad);
1649     ntp_->Branch("v_hsimInfoConeCharHad", "std::vector<std::vector<double> >", &t_v_hsimInfoConeCharHad);
1650     ntp_->Branch("v_hsimInfoConePdgMatched", "std::vector<std::vector<double> >", &t_v_hsimInfoConePdgMatched);
1651     ntp_->Branch("v_hsimInfoConeTotal", "std::vector<std::vector<double> >", &t_v_hsimInfoConeTotal);
1652
1653     ntp_->Branch("v_hsimInfoConeNMatched", "std::vector<std::vector<int> >", &t_v_hsimInfoConeNMatched);
1654     ntp_->Branch("v_hsimInfoConeNRest", "std::vector<std::vector<int> >", &t_v_hsimInfoConeNRest);
1655     ntp_->Branch("v_hsimInfoConeNPhoton", "std::vector<std::vector<int> >", &t_v_hsimInfoConeNPhoton);
1656     ntp_->Branch("v_hsimInfoConeNNeutHad", "std::vector<std::vector<int> >", &t_v_hsimInfoConeNNeutHad);
1657     ntp_->Branch("v_hsimInfoConeNCharHad", "std::vector<std::vector<int> >", &t_v_hsimInfoConeNCharHad);
1658     ntp_->Branch("v_hsimInfoConeNTotal", "std::vector<std::vector<int> >", &t_v_hsimInfoConeNTotal);
1659
1660     ntp_->Branch("v_hsimCone", "std::vector<std::vector<double> >", &t_v_hsimCone);
1661   }
1662
1663   ntp_->Branch("v_hCone", "std::vector<std::vector<double> >", &t_v_hCone);
1664   ntp_->Branch("v_nRecHitsCone", "std::vector<std::vector<int> >", &t_v_nRecHitsCone);
1665   ntp_->Branch("v_nSimHitsCone", "std::vector<std::vector<int> >", &t_v_nSimHitsCone);
1666
1667   ntp_->Branch("v_distFromHotCell", "std::vector<std::vector<double> >", &t_v_distFromHotCell);
1668   ntp_->Branch("v_ietaFromHotCell", "std::vector<std::vector<int> >", &t_v_ietaFromHotCell);
1669   ntp_->Branch("v_iphiFromHotCell", "std::vector<std::vector<int> >", &t_v_iphiFromHotCell);
1670
1671   ntp_->Branch("v_RH_h3x3_ieta", "std::vector<std::vector<int> >", &t_v_RH_h3x3_ieta);
1672   ntp_->Branch("v_RH_h3x3_iphi", "std::vector<std::vector<int> >", &t_v_RH_h3x3_iphi);
1673   ntp_->Branch("v_RH_h3x3_ene", "std::vector<std::vector<double> >", &t_v_RH_h3x3_ene);
1674   ntp_->Branch("v_RH_h5x5_ieta", "std::vector<std::vector<int> >", &t_v_RH_h5x5_ieta);
1675   ntp_->Branch("v_RH_h5x5_iphi", "std::vector<std::vector<int> >", &t_v_RH_h5x5_iphi);
1676   ntp_->Branch("v_RH_h5x5_ene", "std::vector<std::vector<double> >", &t_v_RH_h5x5_ene);
1677   ntp_->Branch("v_RH_r26_ieta", "std::vector<std::vector<int> >", &t_v_RH_r26_ieta);
1678   ntp_->Branch("v_RH_r26_iphi", "std::vector<std::vector<int> >", &t_v_RH_r26_iphi);
1679   ntp_->Branch("v_RH_r26_ene", "std::vector<std::vector<double> >", &t_v_RH_r26_ene);
1680   ntp_->Branch("v_RH_r44_ieta", "std::vector<std::vector<int> >", &t_v_RH_r44_ieta);
1681   ntp_->Branch("v_RH_r44_iphi", "std::vector<std::vector<int> >", &t_v_RH_r44_iphi);
1682   ntp_->Branch("v_RH_r44_ene", "std::vector<std::vector<double> >", &t_v_RH_r44_ene);
1683
1684   ntp_->Branch("v_hlTriggers", "std::vector<std::vector<int> >", &t_v_hlTriggers);
1685   ntp_->Branch("v_hltHB", "std::vector<int>", &t_hltHB);
1686   ntp_->Branch("v_hltHE", "std::vector<int>", &t_hltHE);
1687   ntp_->Branch("v_hltL1Jet15", "std::vector<int>", &t_hltL1Jet15);
1688   ntp_->Branch("v_hltJet30", "std::vector<int>", &t_hltJet30);
1689   ntp_->Branch("v_hltJet50", "std::vector<int>", &t_hltJet50);
1690   ntp_->Branch("v_hltJet80", "std::vector<int>", &t_hltJet80);
1691   ntp_->Branch("v_hltJet110", "std::vector<int>", &t_hltJet110);
1692   ntp_->Branch("v_hltJet140", "std::vector<int>", &t_hltJet140);
1693   ntp_->Branch("v_hltJet180", "std::vector<int>", &t_hltJet180);
1694   ntp_->Branch("v_hltL1SingleEG5", "std::vector<int>", &t_hltL1SingleEG5);
1695   ntp_->Branch("v_hltZeroBias", "std::vector<int>", &t_hltZeroBias);
1696   ntp_->Branch("v_hltMinBiasHcal", "std::vector<int>", &t_hltMinBiasHcal);
1697   ntp_->Branch("v_hltMinBiasEcal", "std::vector<int>", &t_hltMinBiasEcal);
1698   ntp_->Branch("v_hltMinBiasPixel", "std::vector<int>", &t_hltMinBiasPixel);
1699   ntp_->Branch("v_hltSingleIsoTau30_Trk5", "std::vector<int>", &t_hltSingleIsoTau30_Trk5);
1700   ntp_->Branch("v_hltDoubleLooseIsoTau15_Trk5", "std::vector<int>", &t_hltDoubleLooseIsoTau15_Trk5);
1701
1702   ntp_->Branch("irun", "std::vector<unsigned int>", &t_irun);
1703   ntp_->Branch("ievt", "std::vector<unsigned int>", &t_ievt);
1704   ntp_->Branch("ilum", "std::vector<unsigned int>", &t_ilum);
1705 }
1706
1707 void IsolatedTracksCone::printTrack(const reco::Track* pTrack) {
1708   std::string theTrackQuality = "highPurity";
1709   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
1710
1711   edm::LogVerbatim("IsoTrack") << " Reference Point " << pTrack->referencePoint() << "\n"
1712                                << " TrackMmentum " << pTrack->momentum() << " (pt,eta,phi)(" << pTrack->pt() << ","
1713                                << pTrack->eta() << "," << pTrack->phi() << ")"
1714                                << " p " << pTrack->p() << "\n"
1715                                << " Normalized chi2 " << pTrack->normalizedChi2() << "  charge " << pTrack->charge()
1716                                << " qoverp() " << pTrack->qoverp() << "+-" << pTrack->qoverpError() << " d0 "
1717                                << pTrack->d0() << "\n"
1718                                << " NValidHits " << pTrack->numberOfValidHits() << "  NLostHits "
1719                                << pTrack->numberOfLostHits() << " TrackQuality " << pTrack->qualityName(trackQuality_)
1720                                << " " << pTrack->quality(trackQuality_);
1721
1722   if (printTrkHitPattern_) {
1723     const reco::HitPattern& p = pTrack->hitPattern();
1724     for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::TRACK_HITS); i++) {
1725       std::ostringstream st1;
1726       p.printHitPattern(reco::HitPattern::TRACK_HITS, i, st1);
1727       edm::LogVerbatim("IsoTrack") << st1.str();
1728     }
1729   }
1730 }
1731
1732 //define this as a plug-in
1733 DEFINE_FWK_MODULE(IsolatedTracksCone);