IsolatedParticles/plugins/IsolatedTracksNxN.cc

0001 // -*- C++ -*
0002 //
0003 // Package:    IsolatedParticles
0004 // Class:      IsolatedTracksNxN
0005 //
0006 /**\class IsolatedTracksNxN IsolatedTracksNxN.cc Calibration/IsolatedParticles/plugins/IsolatedTracksNxN.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Seema Sharma
0015 //         Created:  Mon Aug 10 15:30:40 CST 2009
0016 //
0017 //
0018
0019 // system include files
0020 #include <cmath>
0021 #include <map>
0022 #include <memory>
0023 #include <sstream>
0024 #include <string>
0025 #include <vector>
0026
0027 // user include files
0028 #include <Math/GenVector/VectorUtil.h>
0029
0030 // root objects
0031 #include "TROOT.h"
0032 #include "TSystem.h"
0033 #include "TFile.h"
0034 #include "TH1F.h"
0035 #include "TH2F.h"
0036 #include "TProfile.h"
0037 #include "TDirectory.h"
0038 #include "TTree.h"
0039
0040 #include "Calibration/IsolatedParticles/interface/CaloSimInfo.h"
0041 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0042 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0043 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0044 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0045 #include "Calibration/IsolatedParticles/interface/eHCALMatrixExtra.h"
0046 #include "Calibration/IsolatedParticles/interface/FindCaloHit.h"
0047 #include "Calibration/IsolatedParticles/interface/MatchingSimTrack.h"
0048
0049 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0050 //L1 trigger Menus etc
0051 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
0052 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
0053 #include "CondFormats/DataRecord/interface/L1GtTriggerMaskAlgoTrigRcd.h"
0054 #include "CondFormats/DataRecord/interface/L1GtTriggerMaskTechTrigRcd.h"
0055 #include "CondFormats/L1TObjects/interface/L1GtTriggerMask.h"
0056 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
0057 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenuFwd.h"
0058 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
0059
0060 #include "DataFormats/Common/interface/Ref.h"
0061 #include "DataFormats/Common/interface/Handle.h"
0062 #include "DataFormats/Math/interface/Point3D.h"
0063 #include "DataFormats/Candidate/interface/Candidate.h"
0064 // muons and tracks
0065 #include "DataFormats/TrackReco/interface/Track.h"
0066 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0067 #include "DataFormats/TrackReco/interface/TrackBase.h"
0068 #include "DataFormats/TrackReco/interface/HitPattern.h"
0069 // Vertices
0070 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0071 #include "DataFormats/VertexReco/interface/Vertex.h"
0072 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0073 // Calorimeters
0074 #include "DataFormats/DetId/interface/DetId.h"
0075 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0076 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0077 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0078 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0079 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0080 // Trigger
0081 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
0082 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
0083 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0084 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0085 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapFwd.h"
0086 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
0087 #include "L1Trigger/GlobalTriggerAnalyzer/interface/L1GtUtils.h"
0088 //L1 objects
0089 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
0090 #include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
0091 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
0092 #include "DataFormats/L1Trigger/interface/L1EmParticleFwd.h"
0093 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
0094 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
0095 // Jets in the event
0096 #include "DataFormats/JetReco/interface/CaloJet.h"
0097 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0098 #include "DataFormats/JetReco/interface/JetExtendedAssociation.h"
0099
0100 #include "DataFormats/Math/interface/deltaPhi.h"
0101 #include "DataFormats/Math/interface/deltaR.h"
0102
0103 #include "FWCore/Framework/interface/Frameworkfwd.h"
0104 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0105 #include "FWCore/Framework/interface/Event.h"
0106 #include "FWCore/Framework/interface/EventSetup.h"
0107 #include "FWCore/Framework/interface/MakerMacros.h"
0108 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0109 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0110 // TFile Service
0111 #include "FWCore/ServiceRegistry/interface/Service.h"
0112 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0113
0114 // ecal / hcal
0115 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0116 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0117 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0118 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0119 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0120 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0121 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0122
0123 #include "MagneticField/Engine/interface/MagneticField.h"
0124 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0125
0126 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
0127 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0128 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0129
0130 // SimHit + SimTrack
0131 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0132 #include "SimDataFormats/Track/interface/SimTrack.h"
0133 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0134 #include "SimDataFormats/Vertex/interface/SimVertex.h"
0135 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0136
0137 // track associator
0138 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0139 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0140 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0141
0142 class IsolatedTracksNxN : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0143 public:
0144   explicit IsolatedTracksNxN(const edm::ParameterSet &);
0145   ~IsolatedTracksNxN() override;
0146
0147   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0148
0149 private:
0150   void beginJob() override;
0151   void analyze(const edm::Event &, const edm::EventSetup &) override;
0152   void endJob() override;
0153
0154   void printTrack(const reco::Track *pTrack);
0155
0156   void bookHistograms();
0157
0158   void clearTreeVectors();
0159
0160 private:
0161   std::unique_ptr<L1GtUtils> m_l1GtUtils;
0162   static constexpr size_t nL1BitsMax = 128;
0163   TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0164
0165   // map of trig bit, algo name and num events passed
0166   std::map<std::pair<unsigned int, std::string>, int> l1AlgoMap_;
0167   std::vector<unsigned int> m_triggerMaskAlgoTrig;
0168
0169   const bool doMC_, writeAllTracks_;
0170   const int myverbose_;
0171   const double pvTracksPtMin_;
0172   const int debugTrks_;
0173   const bool printTrkHitPattern_;
0174   const double minTrackP_, maxTrackEta_;
0175   const bool debugL1Info_, L1TriggerAlgoInfo_;
0176   const double tMinE_, tMaxE_, tMinH_, tMaxH_;
0177   bool initL1_;
0178   int nEventProc_, nbad_;
0179
0180   edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_L1extTauJet_;
0181   edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_L1extCenJet_;
0182   edm::EDGetTokenT<l1extra::L1JetParticleCollection> tok_L1extFwdJet_;
0183
0184   edm::EDGetTokenT<l1extra::L1MuonParticleCollection> tok_L1extMu_;
0185   edm::EDGetTokenT<l1extra::L1EmParticleCollection> tok_L1extIsoEm_;
0186   edm::EDGetTokenT<l1extra::L1EmParticleCollection> tok_L1extNoIsoEm_;
0187
0188   edm::EDGetTokenT<reco::CaloJetCollection> tok_jets_;
0189   edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0190
0191   edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0192   edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0193   edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0194
0195   edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0196   edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0197   edm::EDGetTokenT<edm::SimTrackContainer> tok_simTk_;
0198   edm::EDGetTokenT<edm::SimVertexContainer> tok_simVtx_;
0199   edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloEB_;
0200   edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloEE_;
0201   edm::EDGetTokenT<edm::PCaloHitContainer> tok_caloHH_;
0202
0203   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0204   edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_caloTopology_;
0205   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_topo_;
0206   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
0207   edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_ecalChStatus_;
0208   edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0209   edm::ESGetToken<EcalTrigTowerConstituentsMap, IdealGeometryRecord> tok_htmap_;
0210   static constexpr size_t NPBins = 15;
0211   static constexpr size_t NEtaBins = 3;
0212   double genPartPBins[NPBins + 1], genPartEtaBins[NEtaBins + 1];
0213
0214   TH1F *h_maxNearP15x15[NPBins][NEtaBins], *h_maxNearP21x21[NPBins][NEtaBins], *h_maxNearP25x25[NPBins][NEtaBins],
0215       *h_maxNearP31x31[NPBins][NEtaBins];
0216   TH1I *h_L1AlgoNames;
0217   TH1F *h_PVTracksWt;
0218   TH1F *h_nTracks;
0219   TH1F *h_recPt_0, *h_recP_0, *h_recEta_0, *h_recPhi_0;
0220   TH2F *h_recEtaPt_0, *h_recEtaP_0;
0221   TH1F *h_recPt_1, *h_recP_1, *h_recEta_1, *h_recPhi_1;
0222   TH2F *h_recEtaPt_1, *h_recEtaP_1;
0223   TH1F *h_recPt_2, *h_recP_2, *h_recEta_2, *h_recPhi_2;
0224   TH2F *h_recEtaPt_2, *h_recEtaP_2;
0225
0226   TTree *tree_;
0227
0228   int t_nTracks;
0229   int t_RunNo, t_EvtNo, t_Lumi, t_Bunch;
0230   std::vector<std::string> t_L1AlgoNames;
0231   std::vector<int> t_L1PreScale;
0232   int t_L1Decision[128];
0233
0234   std::vector<double> t_L1CenJetPt, t_L1CenJetEta, t_L1CenJetPhi;
0235   std::vector<double> t_L1FwdJetPt, t_L1FwdJetEta, t_L1FwdJetPhi;
0236   std::vector<double> t_L1TauJetPt, t_L1TauJetEta, t_L1TauJetPhi;
0237   std::vector<double> t_L1MuonPt, t_L1MuonEta, t_L1MuonPhi;
0238   std::vector<double> t_L1IsoEMPt, t_L1IsoEMEta, t_L1IsoEMPhi;
0239   std::vector<double> t_L1NonIsoEMPt, t_L1NonIsoEMEta, t_L1NonIsoEMPhi;
0240   std::vector<double> t_L1METPt, t_L1METEta, t_L1METPhi;
0241
0242   std::vector<double> t_PVx, t_PVy, t_PVz, t_PVTracksSumPt;
0243   std::vector<double> t_PVTracksSumPtWt, t_PVTracksSumPtHP, t_PVTracksSumPtHPWt;
0244   std::vector<int> t_PVisValid, t_PVNTracks, t_PVNTracksWt, t_PVndof;
0245   std::vector<int> t_PVNTracksHP, t_PVNTracksHPWt;
0246
0247   std::vector<double> t_jetPt, t_jetEta, t_jetPhi;
0248   std::vector<double> t_nTrksJetCalo, t_nTrksJetVtx;
0249
0250   std::vector<double> t_trackPAll, t_trackEtaAll, t_trackPhiAll, t_trackPdgIdAll;
0251   std::vector<double> t_trackPtAll;
0252   std::vector<double> t_trackDxyAll, t_trackDzAll, t_trackDxyPVAll, t_trackDzPVAll, t_trackChiSqAll;
0253
0254   std::vector<double> t_trackP, t_trackPt, t_trackEta, t_trackPhi;
0255   std::vector<double> t_trackEcalEta, t_trackEcalPhi, t_trackHcalEta, t_trackHcalPhi;
0256   std::vector<double> t_trackDxy, t_trackDxyBS, t_trackDz, t_trackDzBS;
0257   std::vector<double> t_trackDxyPV, t_trackDzPV;
0258   std::vector<double> t_trackChiSq;
0259   std::vector<int> t_trackPVIdx;
0260
0261   std::vector<int> t_NLayersCrossed, t_trackNOuterHits;
0262   std::vector<int> t_trackHitsTOB, t_trackHitsTEC;
0263   std::vector<int> t_trackHitInMissTOB, t_trackHitInMissTEC, t_trackHitInMissTIB, t_trackHitInMissTID,
0264       t_trackHitInMissTIBTID;
0265   std::vector<int> t_trackHitOutMissTOB, t_trackHitOutMissTEC, t_trackHitOutMissTIB, t_trackHitOutMissTID,
0266       t_trackHitOutMissTOBTEC;
0267   std::vector<int> t_trackHitInMeasTOB, t_trackHitInMeasTEC, t_trackHitInMeasTIB, t_trackHitInMeasTID;
0268   std::vector<int> t_trackHitOutMeasTOB, t_trackHitOutMeasTEC, t_trackHitOutMeasTIB, t_trackHitOutMeasTID;
0269   std::vector<double> t_trackOutPosOutHitDr, t_trackL;
0270
0271   std::vector<double> t_maxNearP31x31;
0272   std::vector<double> t_maxNearP21x21;
0273   std::vector<int> t_ecalSpike11x11;
0274
0275   std::vector<double> t_e7x7, t_e9x9, t_e11x11, t_e15x15;
0276   std::vector<double> t_e7x7_10Sig, t_e9x9_10Sig, t_e11x11_10Sig, t_e15x15_10Sig;
0277   std::vector<double> t_e7x7_15Sig, t_e9x9_15Sig, t_e11x11_15Sig, t_e15x15_15Sig;
0278   std::vector<double> t_e7x7_20Sig, t_e9x9_20Sig, t_e11x11_20Sig, t_e15x15_20Sig;
0279   std::vector<double> t_e7x7_25Sig, t_e9x9_25Sig, t_e11x11_25Sig, t_e15x15_25Sig;
0280   std::vector<double> t_e7x7_30Sig, t_e9x9_30Sig, t_e11x11_30Sig, t_e15x15_30Sig;
0281
0282   std::vector<double> t_esimPdgId, t_simTrackP, t_trkEcalEne;
0283   std::vector<double> t_esim7x7, t_esim9x9, t_esim11x11, t_esim15x15;
0284   std::vector<double> t_esim7x7Matched, t_esim9x9Matched, t_esim11x11Matched, t_esim15x15Matched;
0285   std::vector<double> t_esim7x7Rest, t_esim9x9Rest, t_esim11x11Rest, t_esim15x15Rest;
0286   std::vector<double> t_esim7x7Photon, t_esim9x9Photon, t_esim11x11Photon, t_esim15x15Photon;
0287   std::vector<double> t_esim7x7NeutHad, t_esim9x9NeutHad, t_esim11x11NeutHad, t_esim15x15NeutHad;
0288   std::vector<double> t_esim7x7CharHad, t_esim9x9CharHad, t_esim11x11CharHad, t_esim15x15CharHad;
0289
0290   std::vector<double> t_maxNearHcalP3x3, t_maxNearHcalP5x5, t_maxNearHcalP7x7;
0291   std::vector<double> t_h3x3, t_h5x5, t_h7x7;
0292   std::vector<double> t_h3x3Sig, t_h5x5Sig, t_h7x7Sig;
0293   std::vector<int> t_infoHcal;
0294
0295   std::vector<double> t_trkHcalEne;
0296   std::vector<double> t_hsim3x3, t_hsim5x5, t_hsim7x7;
0297   std::vector<double> t_hsim3x3Matched, t_hsim5x5Matched, t_hsim7x7Matched;
0298   std::vector<double> t_hsim3x3Rest, t_hsim5x5Rest, t_hsim7x7Rest;
0299   std::vector<double> t_hsim3x3Photon, t_hsim5x5Photon, t_hsim7x7Photon;
0300   std::vector<double> t_hsim3x3NeutHad, t_hsim5x5NeutHad, t_hsim7x7NeutHad;
0301   std::vector<double> t_hsim3x3CharHad, t_hsim5x5CharHad, t_hsim7x7CharHad;
0302 };
0303
0304 static const bool useL1EventSetup(true);
0305 static const bool useL1GtTriggerMenuLite(true);
0306
0307 IsolatedTracksNxN::IsolatedTracksNxN(const edm::ParameterSet &iConfig)
0308     : trackerHitAssociatorConfig_(consumesCollector()),
0309       doMC_(iConfig.getUntrackedParameter<bool>("doMC", false)),
0310       writeAllTracks_(iConfig.getUntrackedParameter<bool>("writeAllTracks", false)),
0311       myverbose_(iConfig.getUntrackedParameter<int>("verbosity", 5)),
0312       pvTracksPtMin_(iConfig.getUntrackedParameter<double>("pvTracksPtMin", 1.0)),
0313       debugTrks_(iConfig.getUntrackedParameter<int>("debugTracks", 0)),
0314       printTrkHitPattern_(iConfig.getUntrackedParameter<bool>("printTrkHitPattern", false)),
0315       minTrackP_(iConfig.getUntrackedParameter<double>("minTrackP", 1.0)),
0316       maxTrackEta_(iConfig.getUntrackedParameter<double>("maxTrackEta", 5.0)),
0317       debugL1Info_(iConfig.getUntrackedParameter<bool>("debugL1Info", false)),
0318       L1TriggerAlgoInfo_(iConfig.getUntrackedParameter<bool>("l1TriggerAlgoInfo", false)),
0319       tMinE_(iConfig.getUntrackedParameter<double>("timeMinCutECAL", -500.)),
0320       tMaxE_(iConfig.getUntrackedParameter<double>("timeMaxCutECAL", 500.)),
0321       tMinH_(iConfig.getUntrackedParameter<double>("timeMinCutHCAL", -500.)),
0322       tMaxH_(iConfig.getUntrackedParameter<double>("timeMaxCutHCAL", 500.)) {
0323   if (L1TriggerAlgoInfo_) {
0324     m_l1GtUtils = std::make_unique<L1GtUtils>(
0325         iConfig, consumesCollector(), useL1GtTriggerMenuLite, *this, L1GtUtils::UseEventSetupIn::Event);
0326   }
0327
0328   usesResource(TFileService::kSharedResource);
0329
0330   //now do what ever initialization is needed
0331
0332   edm::InputTag L1extraTauJetSource_ = iConfig.getParameter<edm::InputTag>("l1extraTauJetSource");
0333   edm::InputTag L1extraCenJetSource_ = iConfig.getParameter<edm::InputTag>("l1extraCenJetSource");
0334   edm::InputTag L1extraFwdJetSource_ = iConfig.getParameter<edm::InputTag>("l1extraFwdJetSource");
0335   edm::InputTag L1extraMuonSource_ = iConfig.getParameter<edm::InputTag>("l1extraMuonSource");
0336   edm::InputTag L1extraIsoEmSource_ = iConfig.getParameter<edm::InputTag>("l1extraIsoEmSource");
0337   edm::InputTag L1extraNonIsoEmSource_ = iConfig.getParameter<edm::InputTag>("l1extraNonIsoEmSource");
0338   edm::InputTag L1GTReadoutRcdSource_ = iConfig.getParameter<edm::InputTag>("l1GTReadoutRcdSource");
0339   edm::InputTag L1GTObjectMapRcdSource_ = iConfig.getParameter<edm::InputTag>("l1GTObjectMapRcdSource");
0340   edm::InputTag JetSrc_ = iConfig.getParameter<edm::InputTag>("jetSource");
0341   edm::InputTag JetExtender_ = iConfig.getParameter<edm::InputTag>("jetExtender");
0342   edm::InputTag HBHERecHitSource_ = iConfig.getParameter<edm::InputTag>("hbheRecHitSource");
0343
0344   // define tokens for access
0345   tok_L1extTauJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraTauJetSource_);
0346   tok_L1extCenJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraCenJetSource_);
0347   tok_L1extFwdJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraFwdJetSource_);
0348   tok_L1extMu_ = consumes<l1extra::L1MuonParticleCollection>(L1extraMuonSource_);
0349   tok_L1extIsoEm_ = consumes<l1extra::L1EmParticleCollection>(L1extraIsoEmSource_);
0350   tok_L1extNoIsoEm_ = consumes<l1extra::L1EmParticleCollection>(L1extraNonIsoEmSource_);
0351   tok_jets_ = consumes<reco::CaloJetCollection>(JetSrc_);
0352   tok_hbhe_ = consumes<HBHERecHitCollection>(HBHERecHitSource_);
0353
0354   tok_genTrack_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
0355   tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0356   tok_bs_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
0357   tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0358   tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0359
0360   tok_simTk_ = consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"));
0361   tok_simVtx_ = consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"));
0362   tok_caloEB_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEB"));
0363   tok_caloEE_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEE"));
0364   tok_caloHH_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"));
0365
0366   if (myverbose_ >= 0) {
0367     edm::LogVerbatim("IsoTrack") << "Parameters read from config file \n"
0368                                  << " doMC        " << doMC_ << "\t myverbose " << myverbose_ << "\t minTrackP "
0369                                  << minTrackP_ << "\t maxTrackEta " << maxTrackEta_ << "\t tMinE " << tMinE_
0370                                  << "\t tMaxE " << tMaxE_ << "\t tMinH " << tMinH_ << "\t tMaxH " << tMaxH_
0371                                  << "\n debugL1Info " << debugL1Info_ << "\t L1TriggerAlgoInfo " << L1TriggerAlgoInfo_
0372                                  << "\n";
0373   }
0374
0375   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0376   tok_caloTopology_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0377   tok_topo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0378   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0379   tok_ecalChStatus_ = esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0380   tok_sevlv_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
0381   tok_htmap_ = esConsumes<EcalTrigTowerConstituentsMap, IdealGeometryRecord>();
0382 }
0383
0384 IsolatedTracksNxN::~IsolatedTracksNxN() {}
0385
0386 void IsolatedTracksNxN::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0387   edm::ParameterSetDescription desc;
0388   desc.addUntracked<bool>("doMC", false);
0389   desc.addUntracked<bool>("writeAllTracks", false);
0390   desc.addUntracked<int>("verbosity", 1);
0391   desc.addUntracked<double>("pvTracksPtMin", 0.200);
0392   desc.addUntracked<int>("debugTracks", 0);
0393   desc.addUntracked<bool>("printTrkHitPattern", true);
0394   desc.addUntracked<double>("minTrackP", 1.0);
0395   desc.addUntracked<double>("maxTrackEta", 2.6);
0396   desc.addUntracked<bool>("debugL1Info", false);
0397   desc.addUntracked<bool>("l1TriggerAlgoInfo", false);
0398   desc.add<edm::InputTag>("l1extraTauJetSource", edm::InputTag("l1extraParticles", "Tau"));
0399   desc.add<edm::InputTag>("l1extraCenJetSource", edm::InputTag("l1extraParticles", "Central"));
0400   desc.add<edm::InputTag>("l1extraFwdJetSource", edm::InputTag("l1extraParticles", "Forward"));
0401   desc.add<edm::InputTag>("l1extraMuonSource", edm::InputTag("l1extraParticles"));
0402   desc.add<edm::InputTag>("l1extraIsoEmSource", edm::InputTag("l1extraParticles", "Isolated"));
0403   desc.add<edm::InputTag>("l1extraNonIsoEmSource", edm::InputTag("l1extraParticles", "NonIsolated"));
0404   desc.add<edm::InputTag>("l1GTReadoutRcdSource", edm::InputTag("gtDigis"));
0405   desc.add<edm::InputTag>("l1GTObjectMapRcdSource", edm::InputTag("hltL1GtObjectMap"));
0406   desc.add<edm::InputTag>("jetSource", edm::InputTag("iterativeCone5CaloJets"));
0407   desc.add<edm::InputTag>("jetExtender", edm::InputTag("iterativeCone5JetExtender"));
0408   desc.add<edm::InputTag>("hbheRecHitSource", edm::InputTag("hbhereco"));
0409   desc.addUntracked<double>("maxNearTrackPT", 1.0);
0410   desc.addUntracked<double>("timeMinCutECAL", -500.0);
0411   desc.addUntracked<double>("timeMaxCutECAL", 500.0);
0412   desc.addUntracked<double>("timeMinCutHCAL", -500.0);
0413   desc.addUntracked<double>("timeMaxCutHCAL", 500.0);
0414   descriptions.add("isolatedTracksNxN", desc);
0415 }
0416
0417 void IsolatedTracksNxN::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0418   bool haveIsoTrack = false;
0419
0420   const MagneticField *bField = &iSetup.getData(tok_magField_);
0421
0422   clearTreeVectors();
0423
0424   t_RunNo = iEvent.id().run();
0425   t_EvtNo = iEvent.id().event();
0426   t_Lumi = iEvent.luminosityBlock();
0427   t_Bunch = iEvent.bunchCrossing();
0428
0429   ++nEventProc_;
0430
0431   edm::Handle<reco::TrackCollection> trkCollection;
0432   iEvent.getByToken(tok_genTrack_, trkCollection);
0433   if (debugTrks_ > 1) {
0434     edm::LogVerbatim("IsoTrack") << "Track Collection: ";
0435     edm::LogVerbatim("IsoTrack") << "Number of Tracks " << trkCollection->size();
0436   }
0437   std::string theTrackQuality = "highPurity";
0438   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0439
0440   //===================== save L1 Trigger information =======================
0441   if (L1TriggerAlgoInfo_) {
0442     m_l1GtUtils->getL1GtRunCache(iEvent, iSetup, useL1EventSetup, useL1GtTriggerMenuLite);
0443
0444     int iErrorCode = -1;
0445     int l1ConfCode = -1;
0446     const bool l1Conf = m_l1GtUtils->availableL1Configuration(iErrorCode, l1ConfCode);
0447     if (!l1Conf) {
0448       edm::LogVerbatim("IsoTrack")
0449           << "\nL1 configuration code:" << l1ConfCode << "\nNo valid L1 trigger configuration available."
0450           << "\nSee text above for error code interpretation"
0451           << "\nNo return here, in order to test each method, protected against configuration error.";
0452     }
0453
0454     const L1GtTriggerMenu *m_l1GtMenu = m_l1GtUtils->ptrL1TriggerMenuEventSetup(iErrorCode);
0455     const AlgorithmMap &algorithmMap = m_l1GtMenu->gtAlgorithmMap();
0456     const std::string &menuName = m_l1GtMenu->gtTriggerMenuName();
0457
0458     if (!initL1_) {
0459       initL1_ = true;
0460       edm::LogVerbatim("IsoTrack") << "menuName " << menuName;
0461       for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end(); itAlgo++) {
0462         std::string algName = itAlgo->first;
0463         int algBitNumber = (itAlgo->second).algoBitNumber();
0464         l1AlgoMap_.insert(std::pair<std::pair<unsigned int, std::string>, int>(
0465             std::pair<unsigned int, std::string>(algBitNumber, algName), 0));
0466       }
0467       std::map<std::pair<unsigned int, std::string>, int>::iterator itr;
0468       for (itr = l1AlgoMap_.begin(); itr != l1AlgoMap_.end(); itr++) {
0469         edm::LogVerbatim("IsoTrack") << " ********** " << (itr->first).first << " " << (itr->first).second << " "
0470                                      << itr->second;
0471       }
0472     }
0473
0474     std::vector<int> algbits;
0475     for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end(); itAlgo++) {
0476       std::string algName = itAlgo->first;
0477       int algBitNumber = (itAlgo->second).algoBitNumber();
0478       bool decision = m_l1GtUtils->decision(iEvent, itAlgo->first, iErrorCode);
0479       int preScale = m_l1GtUtils->prescaleFactor(iEvent, itAlgo->first, iErrorCode);
0480
0481       // save the algo names which fired
0482       if (decision) {
0483         l1AlgoMap_[std::pair<unsigned int, std::string>(algBitNumber, algName)] += 1;
0484         h_L1AlgoNames->Fill(algBitNumber);
0485         t_L1AlgoNames.push_back(itAlgo->first);
0486         t_L1PreScale.push_back(preScale);
0487         t_L1Decision[algBitNumber] = 1;
0488         algbits.push_back(algBitNumber);
0489       }
0490     }
0491
0492     if (debugL1Info_) {
0493       for (unsigned int ii = 0; ii < t_L1AlgoNames.size(); ii++) {
0494         edm::LogVerbatim("IsoTrack") << ii << " " << t_L1AlgoNames[ii] << " " << t_L1PreScale[ii] << " " << algbits[ii];
0495       }
0496       for (int i = 0; i < 128; ++i) {
0497         edm::LogVerbatim("IsoTrack") << "L1Decision: " << i << ":" << t_L1Decision[i];
0498       }
0499     }
0500
0501     // L1Taus
0502     const edm::Handle<l1extra::L1JetParticleCollection> &l1TauHandle = iEvent.getHandle(tok_L1extTauJet_);
0503     l1extra::L1JetParticleCollection::const_iterator itr;
0504     int iL1Obj = 0;
0505     for (itr = l1TauHandle->begin(), iL1Obj = 0; itr != l1TauHandle->end(); ++itr, iL1Obj++) {
0506       if (iL1Obj < 1) {
0507         t_L1TauJetPt.push_back(itr->pt());
0508         t_L1TauJetEta.push_back(itr->eta());
0509         t_L1TauJetPhi.push_back(itr->phi());
0510       }
0511       if (debugL1Info_) {
0512         edm::LogVerbatim("IsoTrack") << "tauJ p/pt  " << itr->momentum() << " " << itr->pt() << "  eta/phi "
0513                                      << itr->eta() << " " << itr->phi();
0514       }
0515     }
0516
0517     // L1 Central Jets
0518     const edm::Handle<l1extra::L1JetParticleCollection> &l1CenJetHandle = iEvent.getHandle(tok_L1extCenJet_);
0519     for (itr = l1CenJetHandle->begin(), iL1Obj = 0; itr != l1CenJetHandle->end(); ++itr, iL1Obj++) {
0520       if (iL1Obj < 1) {
0521         t_L1CenJetPt.push_back(itr->pt());
0522         t_L1CenJetEta.push_back(itr->eta());
0523         t_L1CenJetPhi.push_back(itr->phi());
0524       }
0525       if (debugL1Info_) {
0526         edm::LogVerbatim("IsoTrack") << "cenJ p/pt     " << itr->momentum() << " " << itr->pt() << "  eta/phi "
0527                                      << itr->eta() << " " << itr->phi();
0528       }
0529     }
0530
0531     // L1 Forward Jets
0532     const edm::Handle<l1extra::L1JetParticleCollection> &l1FwdJetHandle = iEvent.getHandle(tok_L1extFwdJet_);
0533     for (itr = l1FwdJetHandle->begin(), iL1Obj = 0; itr != l1FwdJetHandle->end(); ++itr, iL1Obj++) {
0534       if (iL1Obj < 1) {
0535         t_L1FwdJetPt.push_back(itr->pt());
0536         t_L1FwdJetEta.push_back(itr->eta());
0537         t_L1FwdJetPhi.push_back(itr->phi());
0538       }
0539       if (debugL1Info_) {
0540         edm::LogVerbatim("IsoTrack") << "fwdJ p/pt     " << itr->momentum() << " " << itr->pt() << "  eta/phi "
0541                                      << itr->eta() << " " << itr->phi();
0542       }
0543     }
0544
0545     // L1 Isolated EM onjects
0546     l1extra::L1EmParticleCollection::const_iterator itrEm;
0547     const edm::Handle<l1extra::L1EmParticleCollection> &l1IsoEmHandle = iEvent.getHandle(tok_L1extIsoEm_);
0548     for (itrEm = l1IsoEmHandle->begin(), iL1Obj = 0; itrEm != l1IsoEmHandle->end(); ++itrEm, iL1Obj++) {
0549       if (iL1Obj < 1) {
0550         t_L1IsoEMPt.push_back(itrEm->pt());
0551         t_L1IsoEMEta.push_back(itrEm->eta());
0552         t_L1IsoEMPhi.push_back(itrEm->phi());
0553       }
0554       if (debugL1Info_) {
0555         edm::LogVerbatim("IsoTrack") << "isoEm p/pt    " << itrEm->momentum() << " " << itrEm->pt() << "  eta/phi "
0556                                      << itrEm->eta() << " " << itrEm->phi();
0557       }
0558     }
0559
0560     // L1 Non-Isolated EM onjects
0561     const edm::Handle<l1extra::L1EmParticleCollection> &l1NonIsoEmHandle = iEvent.getHandle(tok_L1extNoIsoEm_);
0562     for (itrEm = l1NonIsoEmHandle->begin(), iL1Obj = 0; itrEm != l1NonIsoEmHandle->end(); ++itrEm, iL1Obj++) {
0563       if (iL1Obj < 1) {
0564         t_L1NonIsoEMPt.push_back(itrEm->pt());
0565         t_L1NonIsoEMEta.push_back(itrEm->eta());
0566         t_L1NonIsoEMPhi.push_back(itrEm->phi());
0567       }
0568       if (debugL1Info_) {
0569         edm::LogVerbatim("IsoTrack") << "nonIsoEm p/pt " << itrEm->momentum() << " " << itrEm->pt() << "  eta/phi "
0570                                      << itrEm->eta() << " " << itrEm->phi();
0571       }
0572     }
0573
0574     // L1 Muons
0575     l1extra::L1MuonParticleCollection::const_iterator itrMu;
0576     const edm::Handle<l1extra::L1MuonParticleCollection> &l1MuHandle = iEvent.getHandle(tok_L1extMu_);
0577     for (itrMu = l1MuHandle->begin(), iL1Obj = 0; itrMu != l1MuHandle->end(); ++itrMu, iL1Obj++) {
0578       if (iL1Obj < 1) {
0579         t_L1MuonPt.push_back(itrMu->pt());
0580         t_L1MuonEta.push_back(itrMu->eta());
0581         t_L1MuonPhi.push_back(itrMu->phi());
0582       }
0583       if (debugL1Info_) {
0584         edm::LogVerbatim("IsoTrack") << "l1muon p/pt   " << itrMu->momentum() << " " << itrMu->pt() << "  eta/phi "
0585                                      << itrMu->eta() << " " << itrMu->phi();
0586       }
0587     }
0588   }
0589
0590   //============== store the information about all the Non-Fake vertices ===============
0591
0592   const edm::Handle<reco::VertexCollection> &recVtxs = iEvent.getHandle(tok_recVtx_);
0593
0594   std::vector<reco::Track> svTracks;
0595   math::XYZPoint leadPV(0, 0, 0);
0596   double sumPtMax = -1.0;
0597   for (unsigned int ind = 0; ind < recVtxs->size(); ind++) {
0598     if (!((*recVtxs)[ind].isFake())) {
0599       double vtxTrkSumPt = 0.0, vtxTrkSumPtWt = 0.0;
0600       int vtxTrkNWt = 0;
0601       double vtxTrkSumPtHP = 0.0, vtxTrkSumPtHPWt = 0.0;
0602       int vtxTrkNHP = 0, vtxTrkNHPWt = 0;
0603
0604       reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
0605
0606       for (vtxTrack = (*recVtxs)[ind].tracks_begin(); vtxTrack != (*recVtxs)[ind].tracks_end(); vtxTrack++) {
0607         if ((*vtxTrack)->pt() < pvTracksPtMin_)
0608           continue;
0609
0610         bool trkQuality = (*vtxTrack)->quality(trackQuality_);
0611
0612         vtxTrkSumPt += (*vtxTrack)->pt();
0613
0614         vtxTrkSumPt += (*vtxTrack)->pt();
0615         if (trkQuality) {
0616           vtxTrkSumPtHP += (*vtxTrack)->pt();
0617           vtxTrkNHP++;
0618         }
0619
0620         double weight = (*recVtxs)[ind].trackWeight(*vtxTrack);
0621         h_PVTracksWt->Fill(weight);
0622         if (weight > 0.5) {
0623           vtxTrkSumPtWt += (*vtxTrack)->pt();
0624           vtxTrkNWt++;
0625           if (trkQuality) {
0626             vtxTrkSumPtHPWt += (*vtxTrack)->pt();
0627             vtxTrkNHPWt++;
0628           }
0629         }
0630       }
0631
0632       if (vtxTrkSumPt > sumPtMax) {
0633         sumPtMax = vtxTrkSumPt;
0634         leadPV = math::XYZPoint((*recVtxs)[ind].x(), (*recVtxs)[ind].y(), (*recVtxs)[ind].z());
0635       }
0636
0637       t_PVx.push_back((*recVtxs)[ind].x());
0638       t_PVy.push_back((*recVtxs)[ind].y());
0639       t_PVz.push_back((*recVtxs)[ind].z());
0640       t_PVisValid.push_back((*recVtxs)[ind].isValid());
0641       t_PVNTracks.push_back((*recVtxs)[ind].tracksSize());
0642       t_PVndof.push_back((*recVtxs)[ind].ndof());
0643       t_PVNTracksWt.push_back(vtxTrkNWt);
0644       t_PVTracksSumPt.push_back(vtxTrkSumPt);
0645       t_PVTracksSumPtWt.push_back(vtxTrkSumPtWt);
0646
0647       t_PVNTracksHP.push_back(vtxTrkNHP);
0648       t_PVNTracksHPWt.push_back(vtxTrkNHPWt);
0649       t_PVTracksSumPtHP.push_back(vtxTrkSumPtHP);
0650       t_PVTracksSumPtHPWt.push_back(vtxTrkSumPtHPWt);
0651
0652       if (myverbose_ == 4) {
0653         edm::LogVerbatim("IsoTrack") << "PV " << ind << " isValid " << (*recVtxs)[ind].isValid() << " isFake "
0654                                      << (*recVtxs)[ind].isFake() << " hasRefittedTracks() " << ind << " "
0655                                      << (*recVtxs)[ind].hasRefittedTracks() << " refittedTrksSize "
0656                                      << (*recVtxs)[ind].refittedTracks().size() << "  tracksSize() "
0657                                      << (*recVtxs)[ind].tracksSize() << " sumPt " << vtxTrkSumPt;
0658       }
0659     }  // if vtx is not Fake
0660   }  // loop over PVs
0661   //===================================================================================
0662
0663   // Get the beamspot
0664   const edm::Handle<reco::BeamSpot> &beamSpotH = iEvent.getHandle(tok_bs_);
0665   math::XYZPoint bspot;
0666   bspot = (beamSpotH.isValid()) ? beamSpotH->position() : math::XYZPoint(0, 0, 0);
0667
0668   //=====================================================================
0669
0670   const edm::Handle<reco::CaloJetCollection> &jets = iEvent.getHandle(tok_jets_);
0671
0672   for (unsigned int ijet = 0; ijet < (*jets).size(); ijet++) {
0673     t_jetPt.push_back((*jets)[ijet].pt());
0674     t_jetEta.push_back((*jets)[ijet].eta());
0675     t_jetPhi.push_back((*jets)[ijet].phi());
0676     t_nTrksJetVtx.push_back(-1.0);
0677     t_nTrksJetCalo.push_back(-1.0);
0678   }
0679
0680   //===================================================================================
0681
0682   // get handles to calogeometry and calotopology
0683   const CaloGeometry *geo = &iSetup.getData(tok_geom_);
0684   const CaloTopology *caloTopology = &iSetup.getData(tok_caloTopology_);
0685   const HcalTopology *theHBHETopology = &iSetup.getData(tok_topo_);
0686
0687   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
0688   iEvent.getByToken(tok_EB_, barrelRecHitsHandle);
0689   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
0690   iEvent.getByToken(tok_EE_, endcapRecHitsHandle);
0691
0692   // Retrieve the good/bad ECAL channels from the DB
0693   const EcalChannelStatus *theEcalChStatus = &iSetup.getData(tok_ecalChStatus_);
0694   const EcalSeverityLevelAlgo *sevlv = &iSetup.getData(tok_sevlv_);
0695
0696   // Retrieve trigger tower map
0697   const EcalTrigTowerConstituentsMap &ttMap = iSetup.getData(tok_htmap_);
0698
0699   edm::Handle<HBHERecHitCollection> hbhe;
0700   iEvent.getByToken(tok_hbhe_, hbhe);
0701   if (!hbhe.isValid()) {
0702     ++nbad_;
0703     if (nbad_ < 10)
0704       edm::LogVerbatim("IsoTrack") << "No HBHE rechit collection";
0705     return;
0706   }
0707   const HBHERecHitCollection Hithbhe = *(hbhe.product());
0708
0709   //get Handles to SimTracks and SimHits
0710   edm::Handle<edm::SimTrackContainer> SimTk;
0711   if (doMC_)
0712     iEvent.getByToken(tok_simTk_, SimTk);
0713
0714   edm::Handle<edm::SimVertexContainer> SimVtx;
0715   if (doMC_)
0716     iEvent.getByToken(tok_simVtx_, SimVtx);
0717
0718   //get Handles to PCaloHitContainers of eb/ee/hbhe
0719   edm::Handle<edm::PCaloHitContainer> pcaloeb;
0720   if (doMC_)
0721     iEvent.getByToken(tok_caloEB_, pcaloeb);
0722
0723   edm::Handle<edm::PCaloHitContainer> pcaloee;
0724   if (doMC_)
0725     iEvent.getByToken(tok_caloEE_, pcaloee);
0726
0727   edm::Handle<edm::PCaloHitContainer> pcalohh;
0728   if (doMC_)
0729     iEvent.getByToken(tok_caloHH_, pcalohh);
0730
0731   //associates tracker rechits/simhits to a track
0732   std::unique_ptr<TrackerHitAssociator> associate;
0733   if (doMC_)
0734     associate = std::make_unique<TrackerHitAssociator>(iEvent, trackerHitAssociatorConfig_);
0735
0736   //===================================================================================
0737
0738   h_nTracks->Fill(trkCollection->size());
0739
0740   int nTracks = 0;
0741
0742   t_nTracks = trkCollection->size();
0743
0744   // get the list of DetIds closest to the impact point of track on surface calorimeters
0745   std::vector<spr::propagatedTrackID> trkCaloDets;
0746   spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, trkCaloDets, false);
0747   std::vector<spr::propagatedTrackID>::const_iterator trkDetItr;
0748
0749   if (myverbose_ > 2) {
0750     for (trkDetItr = trkCaloDets.begin(); trkDetItr != trkCaloDets.end(); trkDetItr++) {
0751       edm::LogVerbatim("IsoTrack") << trkDetItr->trkItr->p() << " " << trkDetItr->trkItr->eta() << " "
0752                                    << trkDetItr->okECAL << " " << trkDetItr->okHCAL;
0753       if (trkDetItr->okECAL) {
0754         if (trkDetItr->detIdECAL.subdetId() == EcalBarrel)
0755           edm::LogVerbatim("IsoTrack") << (EBDetId)trkDetItr->detIdECAL;
0756         else
0757           edm::LogVerbatim("IsoTrack") << (EEDetId)trkDetItr->detIdECAL;
0758       }
0759       if (trkDetItr->okHCAL)
0760         edm::LogVerbatim("IsoTrack") << (HcalDetId)trkDetItr->detIdHCAL;
0761     }
0762   }
0763
0764   int nvtxTracks = 0;
0765   for (trkDetItr = trkCaloDets.begin(), nTracks = 0; trkDetItr != trkCaloDets.end(); trkDetItr++, nTracks++) {
0766     const reco::Track *pTrack = &(*(trkDetItr->trkItr));
0767
0768     // find vertex index the track is associated with
0769     int pVtxTkId = -1;
0770     for (unsigned int ind = 0; ind < recVtxs->size(); ind++) {
0771       if (!((*recVtxs)[ind].isFake())) {
0772         reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
0773         for (vtxTrack = (*recVtxs)[ind].tracks_begin(); vtxTrack != (*recVtxs)[ind].tracks_end(); vtxTrack++) {
0774           const edm::RefToBase<reco::Track> &pvtxTrack = (*vtxTrack);
0775           if (pTrack == pvtxTrack.get()) {
0776             pVtxTkId = ind;
0777             break;
0778             if (myverbose_ > 2) {
0779               if (pTrack->pt() > 1.0) {
0780                 edm::LogVerbatim("IsoTrack") << "Debug the track association with vertex ";
0781                 edm::LogVerbatim("IsoTrack") << pTrack << " " << pvtxTrack.get();
0782                 edm::LogVerbatim("IsoTrack")
0783                     << " trkVtxIndex " << nvtxTracks << " vtx " << ind << " pt " << pTrack->pt() << " eta "
0784                     << pTrack->eta() << " " << pTrack->pt() - pvtxTrack->pt() << " "
0785                     << pTrack->eta() - pvtxTrack->eta();
0786                 nvtxTracks++;
0787               }
0788             }
0789           }
0790         }
0791       }
0792     }
0793
0794     const reco::HitPattern &hitp = pTrack->hitPattern();
0795     int nLayersCrossed = hitp.trackerLayersWithMeasurement();
0796     int nOuterHits = hitp.stripTOBLayersWithMeasurement() + hitp.stripTECLayersWithMeasurement();
0797
0798     bool ifGood = pTrack->quality(trackQuality_);
0799     double pt1 = pTrack->pt();
0800     double p1 = pTrack->p();
0801     double eta1 = pTrack->momentum().eta();
0802     double phi1 = pTrack->momentum().phi();
0803     double etaEcal1 = trkDetItr->etaECAL;
0804     double phiEcal1 = trkDetItr->phiECAL;
0805     double etaHcal1 = trkDetItr->etaHCAL;
0806     double phiHcal1 = trkDetItr->phiHCAL;
0807     double dxy1 = pTrack->dxy();
0808     double dz1 = pTrack->dz();
0809     double chisq1 = pTrack->normalizedChi2();
0810     double dxybs1 = beamSpotH.isValid() ? pTrack->dxy(bspot) : pTrack->dxy();
0811     double dzbs1 = beamSpotH.isValid() ? pTrack->dz(bspot) : pTrack->dz();
0812     double dxypv1 = pTrack->dxy();
0813     double dzpv1 = pTrack->dz();
0814     if (pVtxTkId >= 0) {
0815       math::XYZPoint thisTkPV =
0816           math::XYZPoint((*recVtxs)[pVtxTkId].x(), (*recVtxs)[pVtxTkId].y(), (*recVtxs)[pVtxTkId].z());
0817       dxypv1 = pTrack->dxy(thisTkPV);
0818       dzpv1 = pTrack->dz(thisTkPV);
0819     }
0820
0821     h_recEtaPt_0->Fill(eta1, pt1);
0822     h_recEtaP_0->Fill(eta1, p1);
0823     h_recPt_0->Fill(pt1);
0824     h_recP_0->Fill(p1);
0825     h_recEta_0->Fill(eta1);
0826     h_recPhi_0->Fill(phi1);
0827
0828     if (ifGood && nLayersCrossed > 7) {
0829       h_recEtaPt_1->Fill(eta1, pt1);
0830       h_recEtaP_1->Fill(eta1, p1);
0831       h_recPt_1->Fill(pt1);
0832       h_recP_1->Fill(p1);
0833       h_recEta_1->Fill(eta1);
0834       h_recPhi_1->Fill(phi1);
0835     }
0836
0837     if (!ifGood)
0838       continue;
0839
0840     if (writeAllTracks_ && p1 > 2.0 && nLayersCrossed > 7) {
0841       t_trackPAll.push_back(p1);
0842       t_trackEtaAll.push_back(eta1);
0843       t_trackPhiAll.push_back(phi1);
0844       t_trackPtAll.push_back(pt1);
0845       t_trackDxyAll.push_back(dxy1);
0846       t_trackDzAll.push_back(dz1);
0847       t_trackDxyPVAll.push_back(dxypv1);
0848       t_trackDzPVAll.push_back(dzpv1);
0849       t_trackChiSqAll.push_back(chisq1);
0850     }
0851     if (doMC_) {
0852       edm::SimTrackContainer::const_iterator matchedSimTrkAll =
0853           spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false);
0854       if (writeAllTracks_ && matchedSimTrkAll != SimTk->end())
0855         t_trackPdgIdAll.push_back(matchedSimTrkAll->type());
0856     }
0857
0858     if (pt1 > minTrackP_ && std::abs(eta1) < maxTrackEta_ && trkDetItr->okECAL) {
0859       double maxNearP31x31 = 999.0, maxNearP25x25 = 999.0, maxNearP21x21 = 999.0, maxNearP15x15 = 999.0;
0860       maxNearP31x31 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 15, 15);
0861       maxNearP25x25 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 12, 12);
0862       maxNearP21x21 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 10, 10);
0863       maxNearP15x15 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 7, 7);
0864
0865       int iTrkEtaBin = -1, iTrkMomBin = -1;
0866       for (unsigned int ieta = 0; ieta < NEtaBins; ieta++) {
0867         if (std::abs(eta1) > genPartEtaBins[ieta] && std::abs(eta1) < genPartEtaBins[ieta + 1])
0868           iTrkEtaBin = ieta;
0869       }
0870       for (unsigned int ipt = 0; ipt < NPBins; ipt++) {
0871         if (p1 > genPartPBins[ipt] && p1 < genPartPBins[ipt + 1])
0872           iTrkMomBin = ipt;
0873       }
0874       if (iTrkMomBin >= 0 && iTrkEtaBin >= 0) {
0875         h_maxNearP31x31[iTrkMomBin][iTrkEtaBin]->Fill(maxNearP31x31);
0876         h_maxNearP25x25[iTrkMomBin][iTrkEtaBin]->Fill(maxNearP25x25);
0877         h_maxNearP21x21[iTrkMomBin][iTrkEtaBin]->Fill(maxNearP21x21);
0878         h_maxNearP15x15[iTrkMomBin][iTrkEtaBin]->Fill(maxNearP15x15);
0879       }
0880       if (maxNearP31x31 < 0.0 && nLayersCrossed > 7 && nOuterHits > 4) {
0881         h_recEtaPt_2->Fill(eta1, pt1);
0882         h_recEtaP_2->Fill(eta1, p1);
0883         h_recPt_2->Fill(pt1);
0884         h_recP_2->Fill(p1);
0885         h_recEta_2->Fill(eta1);
0886         h_recPhi_2->Fill(phi1);
0887       }
0888
0889       // if charge isolated in ECAL, store the further quantities
0890       if (maxNearP31x31 < 1.0) {
0891         haveIsoTrack = true;
0892
0893         // get the matching simTrack
0894         double simTrackP = -1;
0895         if (doMC_) {
0896           edm::SimTrackContainer::const_iterator matchedSimTrk =
0897               spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false);
0898           if (matchedSimTrk != SimTk->end())
0899             simTrackP = matchedSimTrk->momentum().P();
0900         }
0901
0902         // get ECal Tranverse Profile
0903         std::pair<double, bool> e7x7P, e9x9P, e11x11P, e15x15P;
0904         std::pair<double, bool> e7x7_10SigP, e9x9_10SigP, e11x11_10SigP, e15x15_10SigP;
0905         std::pair<double, bool> e7x7_15SigP, e9x9_15SigP, e11x11_15SigP, e15x15_15SigP;
0906         std::pair<double, bool> e7x7_20SigP, e9x9_20SigP, e11x11_20SigP, e15x15_20SigP;
0907         std::pair<double, bool> e7x7_25SigP, e9x9_25SigP, e11x11_25SigP, e15x15_25SigP;
0908         std::pair<double, bool> e7x7_30SigP, e9x9_30SigP, e11x11_30SigP, e15x15_30SigP;
0909
0910         spr::caloSimInfo simInfo3x3, simInfo5x5, simInfo7x7, simInfo9x9;
0911         spr::caloSimInfo simInfo11x11, simInfo13x13, simInfo15x15, simInfo21x21, simInfo25x25, simInfo31x31;
0912         double trkEcalEne = 0;
0913
0914         const DetId isoCell = trkDetItr->detIdECAL;
0915         e7x7P = spr::eECALmatrix(isoCell,
0916                                  barrelRecHitsHandle,
0917                                  endcapRecHitsHandle,
0918                                  *theEcalChStatus,
0919                                  geo,
0920                                  caloTopology,
0921                                  sevlv,
0922                                  3,
0923                                  3,
0924                                  -100.0,
0925                                  -100.0,
0926                                  tMinE_,
0927                                  tMaxE_);
0928         e9x9P = spr::eECALmatrix(isoCell,
0929                                  barrelRecHitsHandle,
0930                                  endcapRecHitsHandle,
0931                                  *theEcalChStatus,
0932                                  geo,
0933                                  caloTopology,
0934                                  sevlv,
0935                                  4,
0936                                  4,
0937                                  -100.0,
0938                                  -100.0,
0939                                  tMinE_,
0940                                  tMaxE_);
0941         e11x11P = spr::eECALmatrix(isoCell,
0942                                    barrelRecHitsHandle,
0943                                    endcapRecHitsHandle,
0944                                    *theEcalChStatus,
0945                                    geo,
0946                                    caloTopology,
0947                                    sevlv,
0948                                    5,
0949                                    5,
0950                                    -100.0,
0951                                    -100.0,
0952                                    tMinE_,
0953                                    tMaxE_);
0954         e15x15P = spr::eECALmatrix(isoCell,
0955                                    barrelRecHitsHandle,
0956                                    endcapRecHitsHandle,
0957                                    *theEcalChStatus,
0958                                    geo,
0959                                    caloTopology,
0960                                    sevlv,
0961                                    7,
0962                                    7,
0963                                    -100.0,
0964                                    -100.0,
0965                                    tMinE_,
0966                                    tMaxE_);
0967
0968         e7x7_10SigP = spr::eECALmatrix(isoCell,
0969                                        barrelRecHitsHandle,
0970                                        endcapRecHitsHandle,
0971                                        *theEcalChStatus,
0972                                        geo,
0973                                        caloTopology,
0974                                        sevlv,
0975                                        3,
0976                                        3,
0977                                        0.030,
0978                                        0.150,
0979                                        tMinE_,
0980                                        tMaxE_);
0981         e9x9_10SigP = spr::eECALmatrix(isoCell,
0982                                        barrelRecHitsHandle,
0983                                        endcapRecHitsHandle,
0984                                        *theEcalChStatus,
0985                                        geo,
0986                                        caloTopology,
0987                                        sevlv,
0988                                        4,
0989                                        4,
0990                                        0.030,
0991                                        0.150,
0992                                        tMinE_,
0993                                        tMaxE_);
0994         e11x11_10SigP = spr::eECALmatrix(isoCell,
0995                                          barrelRecHitsHandle,
0996                                          endcapRecHitsHandle,
0997                                          *theEcalChStatus,
0998                                          geo,
0999                                          caloTopology,
1000                                          sevlv,
1001                                          5,
1002                                          5,
1003                                          0.030,
1004                                          0.150,
1005                                          tMinE_,
1006                                          tMaxE_);
1007         e15x15_10SigP = spr::eECALmatrix(isoCell,
1008                                          barrelRecHitsHandle,
1009                                          endcapRecHitsHandle,
1010                                          *theEcalChStatus,
1011                                          geo,
1012                                          caloTopology,
1013                                          sevlv,
1014                                          7,
1015                                          7,
1016                                          0.030,
1017                                          0.150,
1018                                          tMinE_,
1019                                          tMaxE_);
1020
1021         e7x7_15SigP = spr::eECALmatrix(isoCell,
1022                                        barrelRecHitsHandle,
1023                                        endcapRecHitsHandle,
1024                                        *theEcalChStatus,
1025                                        geo,
1026                                        caloTopology,
1027                                        sevlv,
1028                                        ttMap,
1029                                        3,
1030                                        3,
1031                                        0.20,
1032                                        0.45,
1033                                        tMinE_,
1034                                        tMaxE_);
1035         e9x9_15SigP = spr::eECALmatrix(isoCell,
1036                                        barrelRecHitsHandle,
1037                                        endcapRecHitsHandle,
1038                                        *theEcalChStatus,
1039                                        geo,
1040                                        caloTopology,
1041                                        sevlv,
1042                                        ttMap,
1043                                        4,
1044                                        4,
1045                                        0.20,
1046                                        0.45,
1047                                        tMinE_,
1048                                        tMaxE_);
1049         e11x11_15SigP = spr::eECALmatrix(isoCell,
1050                                          barrelRecHitsHandle,
1051                                          endcapRecHitsHandle,
1052                                          *theEcalChStatus,
1053                                          geo,
1054                                          caloTopology,
1055                                          sevlv,
1056                                          ttMap,
1057                                          5,
1058                                          5,
1059                                          0.20,
1060                                          0.45,
1061                                          tMinE_,
1062                                          tMaxE_);
1063         e15x15_15SigP = spr::eECALmatrix(isoCell,
1064                                          barrelRecHitsHandle,
1065                                          endcapRecHitsHandle,
1066                                          *theEcalChStatus,
1067                                          geo,
1068                                          caloTopology,
1069                                          sevlv,
1070                                          ttMap,
1071                                          7,
1072                                          7,
1073                                          0.20,
1074                                          0.45,
1075                                          tMinE_,
1076                                          tMaxE_,
1077                                          false);
1078
1079         e7x7_20SigP = spr::eECALmatrix(isoCell,
1080                                        barrelRecHitsHandle,
1081                                        endcapRecHitsHandle,
1082                                        *theEcalChStatus,
1083                                        geo,
1084                                        caloTopology,
1085                                        sevlv,
1086                                        3,
1087                                        3,
1088                                        0.060,
1089                                        0.300,
1090                                        tMinE_,
1091                                        tMaxE_);
1092         e9x9_20SigP = spr::eECALmatrix(isoCell,
1093                                        barrelRecHitsHandle,
1094                                        endcapRecHitsHandle,
1095                                        *theEcalChStatus,
1096                                        geo,
1097                                        caloTopology,
1098                                        sevlv,
1099                                        4,
1100                                        4,
1101                                        0.060,
1102                                        0.300,
1103                                        tMinE_,
1104                                        tMaxE_);
1105         e11x11_20SigP = spr::eECALmatrix(isoCell,
1106                                          barrelRecHitsHandle,
1107                                          endcapRecHitsHandle,
1108                                          *theEcalChStatus,
1109                                          geo,
1110                                          caloTopology,
1111                                          sevlv,
1112                                          5,
1113                                          5,
1114                                          0.060,
1115                                          0.300,
1116                                          tMinE_,
1117                                          tMaxE_);
1118         e15x15_20SigP = spr::eECALmatrix(isoCell,
1119                                          barrelRecHitsHandle,
1120                                          endcapRecHitsHandle,
1121                                          *theEcalChStatus,
1122                                          geo,
1123                                          caloTopology,
1124                                          sevlv,
1125                                          7,
1126                                          7,
1127                                          0.060,
1128                                          0.300,
1129                                          tMinE_,
1130                                          tMaxE_);
1131
1132         e7x7_25SigP = spr::eECALmatrix(isoCell,
1133                                        barrelRecHitsHandle,
1134                                        endcapRecHitsHandle,
1135                                        *theEcalChStatus,
1136                                        geo,
1137                                        caloTopology,
1138                                        sevlv,
1139                                        3,
1140                                        3,
1141                                        0.075,
1142                                        0.375,
1143                                        tMinE_,
1144                                        tMaxE_);
1145         e9x9_25SigP = spr::eECALmatrix(isoCell,
1146                                        barrelRecHitsHandle,
1147                                        endcapRecHitsHandle,
1148                                        *theEcalChStatus,
1149                                        geo,
1150                                        caloTopology,
1151                                        sevlv,
1152                                        4,
1153                                        4,
1154                                        0.075,
1155                                        0.375,
1156                                        tMinE_,
1157                                        tMaxE_);
1158         e11x11_25SigP = spr::eECALmatrix(isoCell,
1159                                          barrelRecHitsHandle,
1160                                          endcapRecHitsHandle,
1161                                          *theEcalChStatus,
1162                                          geo,
1163                                          caloTopology,
1164                                          sevlv,
1165                                          5,
1166                                          5,
1167                                          0.075,
1168                                          0.375,
1169                                          tMinE_,
1170                                          tMaxE_);
1171         e15x15_25SigP = spr::eECALmatrix(isoCell,
1172                                          barrelRecHitsHandle,
1173                                          endcapRecHitsHandle,
1174                                          *theEcalChStatus,
1175                                          geo,
1176                                          caloTopology,
1177                                          sevlv,
1178                                          7,
1179                                          7,
1180                                          0.075,
1181                                          0.375,
1182                                          tMinE_,
1183                                          tMaxE_);
1184
1185         e7x7_30SigP = spr::eECALmatrix(isoCell,
1186                                        barrelRecHitsHandle,
1187                                        endcapRecHitsHandle,
1188                                        *theEcalChStatus,
1189                                        geo,
1190                                        caloTopology,
1191                                        sevlv,
1192                                        3,
1193                                        3,
1194                                        0.090,
1195                                        0.450,
1196                                        tMinE_,
1197                                        tMaxE_);
1198         e9x9_30SigP = spr::eECALmatrix(isoCell,
1199                                        barrelRecHitsHandle,
1200                                        endcapRecHitsHandle,
1201                                        *theEcalChStatus,
1202                                        geo,
1203                                        caloTopology,
1204                                        sevlv,
1205                                        4,
1206                                        4,
1207                                        0.090,
1208                                        0.450,
1209                                        tMinE_,
1210                                        tMaxE_);
1211         e11x11_30SigP = spr::eECALmatrix(isoCell,
1212                                          barrelRecHitsHandle,
1213                                          endcapRecHitsHandle,
1214                                          *theEcalChStatus,
1215                                          geo,
1216                                          caloTopology,
1217                                          sevlv,
1218                                          5,
1219                                          5,
1220                                          0.090,
1221                                          0.450,
1222                                          tMinE_,
1223                                          tMaxE_);
1224         e15x15_30SigP = spr::eECALmatrix(isoCell,
1225                                          barrelRecHitsHandle,
1226                                          endcapRecHitsHandle,
1227                                          *theEcalChStatus,
1228                                          geo,
1229                                          caloTopology,
1230                                          sevlv,
1231                                          7,
1232                                          7,
1233                                          0.090,
1234                                          0.450,
1235                                          tMinE_,
1236                                          tMaxE_);
1237         if (myverbose_ == 2) {
1238           edm::LogVerbatim("IsoTrack") << "clean  ecal rechit ";
1239           edm::LogVerbatim("IsoTrack") << "e7x7 " << e7x7P.first << " e9x9 " << e9x9P.first << " e11x11 "
1240                                        << e11x11P.first << " e15x15 " << e15x15P.first;
1241           edm::LogVerbatim("IsoTrack") << "e7x7_10Sig " << e7x7_10SigP.first << " e11x11_10Sig " << e11x11_10SigP.first
1242                                        << " e15x15_10Sig " << e15x15_10SigP.first;
1243         }
1244
1245         if (doMC_) {
1246           // check the energy from SimHits
1247           spr::eECALSimInfo(
1248               iEvent, isoCell, geo, caloTopology, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate, 1, 1, simInfo3x3);
1249           spr::eECALSimInfo(
1250               iEvent, isoCell, geo, caloTopology, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate, 2, 2, simInfo5x5);
1251           spr::eECALSimInfo(
1252               iEvent, isoCell, geo, caloTopology, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate, 3, 3, simInfo7x7);
1253           spr::eECALSimInfo(
1254               iEvent, isoCell, geo, caloTopology, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate, 4, 4, simInfo9x9);
1255           spr::eECALSimInfo(iEvent,
1256                             isoCell,
1257                             geo,
1258                             caloTopology,
1259                             pcaloeb,
1260                             pcaloee,
1261                             SimTk,
1262                             SimVtx,
1263                             pTrack,
1264                             *associate,
1265                             5,
1266                             5,
1267                             simInfo11x11);
1268           spr::eECALSimInfo(iEvent,
1269                             isoCell,
1270                             geo,
1271                             caloTopology,
1272                             pcaloeb,
1273                             pcaloee,
1274                             SimTk,
1275                             SimVtx,
1276                             pTrack,
1277                             *associate,
1278                             6,
1279                             6,
1280                             simInfo13x13);
1281           spr::eECALSimInfo(iEvent,
1282                             isoCell,
1283                             geo,
1284                             caloTopology,
1285                             pcaloeb,
1286                             pcaloee,
1287                             SimTk,
1288                             SimVtx,
1289                             pTrack,
1290                             *associate,
1291                             7,
1292                             7,
1293                             simInfo15x15,
1294                             150.0,
1295                             false);
1296           spr::eECALSimInfo(iEvent,
1297                             isoCell,
1298                             geo,
1299                             caloTopology,
1300                             pcaloeb,
1301                             pcaloee,
1302                             SimTk,
1303                             SimVtx,
1304                             pTrack,
1305                             *associate,
1306                             10,
1307                             10,
1308                             simInfo21x21);
1309           spr::eECALSimInfo(iEvent,
1310                             isoCell,
1311                             geo,
1312                             caloTopology,
1313                             pcaloeb,
1314                             pcaloee,
1315                             SimTk,
1316                             SimVtx,
1317                             pTrack,
1318                             *associate,
1319                             12,
1320                             12,
1321                             simInfo25x25);
1322           spr::eECALSimInfo(iEvent,
1323                             isoCell,
1324                             geo,
1325                             caloTopology,
1326                             pcaloeb,
1327                             pcaloee,
1328                             SimTk,
1329                             SimVtx,
1330                             pTrack,
1331                             *associate,
1332                             15,
1333                             15,
1334                             simInfo31x31);
1335
1336           trkEcalEne = spr::eCaloSimInfo(iEvent, geo, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate);
1337           if (myverbose_ == 1) {
1338             edm::LogVerbatim("IsoTrack") << "Track momentum " << pt1;
1339             edm::LogVerbatim("IsoTrack") << "ecal siminfo ";
1340             edm::LogVerbatim("IsoTrack") << "simInfo3x3: eTotal " << simInfo3x3.eTotal << " eMatched "
1341                                          << simInfo3x3.eMatched << " eRest " << simInfo3x3.eRest << " eGamma "
1342                                          << simInfo3x3.eGamma << " eNeutralHad " << simInfo3x3.eNeutralHad
1343                                          << " eChargedHad " << simInfo3x3.eChargedHad;
1344             edm::LogVerbatim("IsoTrack") << "simInfo5x5: eTotal " << simInfo5x5.eTotal << " eMatched "
1345                                          << simInfo5x5.eMatched << " eRest " << simInfo5x5.eRest << " eGamma "
1346                                          << simInfo5x5.eGamma << " eNeutralHad " << simInfo5x5.eNeutralHad
1347                                          << " eChargedHad " << simInfo5x5.eChargedHad;
1348             edm::LogVerbatim("IsoTrack") << "simInfo7x7: eTotal " << simInfo7x7.eTotal << " eMatched "
1349                                          << simInfo7x7.eMatched << " eRest " << simInfo7x7.eRest << " eGamma "
1350                                          << simInfo7x7.eGamma << " eNeutralHad " << simInfo7x7.eNeutralHad
1351                                          << " eChargedHad " << simInfo7x7.eChargedHad;
1352             edm::LogVerbatim("IsoTrack") << "simInfo9x9: eTotal " << simInfo9x9.eTotal << " eMatched "
1353                                          << simInfo9x9.eMatched << " eRest " << simInfo9x9.eRest << " eGamma "
1354                                          << simInfo9x9.eGamma << " eNeutralHad " << simInfo9x9.eNeutralHad
1355                                          << " eChargedHad " << simInfo9x9.eChargedHad;
1356             edm::LogVerbatim("IsoTrack") << "simInfo11x11: eTotal " << simInfo11x11.eTotal << " eMatched "
1357                                          << simInfo11x11.eMatched << " eRest " << simInfo11x11.eRest << " eGamma "
1358                                          << simInfo11x11.eGamma << " eNeutralHad " << simInfo11x11.eNeutralHad
1359                                          << " eChargedHad " << simInfo11x11.eChargedHad;
1360             edm::LogVerbatim("IsoTrack") << "simInfo15x15: eTotal " << simInfo15x15.eTotal << " eMatched "
1361                                          << simInfo15x15.eMatched << " eRest " << simInfo15x15.eRest << " eGamma "
1362                                          << simInfo15x15.eGamma << " eNeutralHad " << simInfo15x15.eNeutralHad
1363                                          << " eChargedHad " << simInfo15x15.eChargedHad;
1364             edm::LogVerbatim("IsoTrack") << "simInfo31x31: eTotal " << simInfo31x31.eTotal << " eMatched "
1365                                          << simInfo31x31.eMatched << " eRest " << simInfo31x31.eRest << " eGamma "
1366                                          << simInfo31x31.eGamma << " eNeutralHad " << simInfo31x31.eNeutralHad
1367                                          << " eChargedHad " << simInfo31x31.eChargedHad;
1368             edm::LogVerbatim("IsoTrack") << "trkEcalEne" << trkEcalEne;
1369           }
1370         }
1371
1372         // =======  Get HCAL information
1373         double hcalScale = 1.0;
1374         if (std::abs(pTrack->eta()) < 1.4) {
1375           hcalScale = 120.0;
1376         } else {
1377           hcalScale = 135.0;
1378         }
1379
1380         double maxNearHcalP3x3 = -1, maxNearHcalP5x5 = -1, maxNearHcalP7x7 = -1;
1381         maxNearHcalP3x3 = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 1, 1);
1382         maxNearHcalP5x5 = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 2, 2);
1383         maxNearHcalP7x7 = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 3, 3);
1384
1385         double h3x3 = 0, h5x5 = 0, h7x7 = 0;
1386         double h3x3Sig = 0, h5x5Sig = 0, h7x7Sig = 0;
1387         double trkHcalEne = 0;
1388         spr::caloSimInfo hsimInfo3x3, hsimInfo5x5, hsimInfo7x7;
1389
1390         if (trkDetItr->okHCAL) {
1391           const DetId ClosestCell(trkDetItr->detIdHCAL);
1392           // bool includeHO=false, bool algoNew=true, bool debug=false
1393           h3x3 = spr::eHCALmatrix(
1394               theHBHETopology, ClosestCell, hbhe, 1, 1, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_);
1395           h5x5 = spr::eHCALmatrix(
1396               theHBHETopology, ClosestCell, hbhe, 2, 2, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_);
1397           h7x7 = spr::eHCALmatrix(
1398               theHBHETopology, ClosestCell, hbhe, 3, 3, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_);
1399           h3x3Sig = spr::eHCALmatrix(
1400               theHBHETopology, ClosestCell, hbhe, 1, 1, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_, tMaxH_);
1401           h5x5Sig = spr::eHCALmatrix(
1402               theHBHETopology, ClosestCell, hbhe, 2, 2, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_, tMaxH_);
1403           h7x7Sig = spr::eHCALmatrix(
1404               theHBHETopology, ClosestCell, hbhe, 3, 3, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_, tMaxH_);
1405
1406           if (myverbose_ == 2) {
1407             edm::LogVerbatim("IsoTrack") << "HCAL 3x3 " << h3x3 << " " << h3x3Sig << " 5x5 " << h5x5 << " " << h5x5Sig
1408                                          << " 7x7 " << h7x7 << " " << h7x7Sig;
1409           }
1410
1411           if (doMC_) {
1412             spr::eHCALSimInfo(
1413                 iEvent, theHBHETopology, ClosestCell, geo, pcalohh, SimTk, SimVtx, pTrack, *associate, 1, 1, hsimInfo3x3);
1414             spr::eHCALSimInfo(
1415                 iEvent, theHBHETopology, ClosestCell, geo, pcalohh, SimTk, SimVtx, pTrack, *associate, 2, 2, hsimInfo5x5);
1416             spr::eHCALSimInfo(iEvent,
1417                               theHBHETopology,
1418                               ClosestCell,
1419                               geo,
1420                               pcalohh,
1421                               SimTk,
1422                               SimVtx,
1423                               pTrack,
1424                               *associate,
1425                               3,
1426                               3,
1427                               hsimInfo7x7,
1428                               150.0,
1429                               false,
1430                               false);
1431             trkHcalEne = spr::eCaloSimInfo(iEvent, geo, pcalohh, SimTk, SimVtx, pTrack, *associate);
1432             if (myverbose_ == 1) {
1433               edm::LogVerbatim("IsoTrack") << "Hcal siminfo ";
1434               edm::LogVerbatim("IsoTrack")
1435                   << "hsimInfo3x3: eTotal " << hsimInfo3x3.eTotal << " eMatched " << hsimInfo3x3.eMatched << " eRest "
1436                   << hsimInfo3x3.eRest << " eGamma " << hsimInfo3x3.eGamma << " eNeutralHad " << hsimInfo3x3.eNeutralHad
1437                   << " eChargedHad " << hsimInfo3x3.eChargedHad;
1438               edm::LogVerbatim("IsoTrack")
1439                   << "hsimInfo5x5: eTotal " << hsimInfo5x5.eTotal << " eMatched " << hsimInfo5x5.eMatched << " eRest "
1440                   << hsimInfo5x5.eRest << " eGamma " << hsimInfo5x5.eGamma << " eNeutralHad " << hsimInfo5x5.eNeutralHad
1441                   << " eChargedHad " << hsimInfo5x5.eChargedHad;
1442               edm::LogVerbatim("IsoTrack")
1443                   << "hsimInfo7x7: eTotal " << hsimInfo7x7.eTotal << " eMatched " << hsimInfo7x7.eMatched << " eRest "
1444                   << hsimInfo7x7.eRest << " eGamma " << hsimInfo7x7.eGamma << " eNeutralHad " << hsimInfo7x7.eNeutralHad
1445                   << " eChargedHad " << hsimInfo7x7.eChargedHad;
1446               edm::LogVerbatim("IsoTrack") << "trkHcalEne " << trkHcalEne;
1447             }
1448           }
1449
1450           // debug the ecal and hcal matrix
1451           if (myverbose_ == 4) {
1452             edm::LogVerbatim("IsoTrack") << "Run " << iEvent.id().run() << "  Event " << iEvent.id().event();
1453             std::vector<std::pair<DetId, double> > v7x7 =
1454                 spr::eHCALmatrixCell(theHBHETopology, ClosestCell, hbhe, 3, 3, false, false);
1455             double sumv = 0.0;
1456
1457             for (unsigned int iv = 0; iv < v7x7.size(); iv++) {
1458               sumv += v7x7[iv].second;
1459             }
1460             edm::LogVerbatim("IsoTrack") << "h7x7 " << h7x7 << " v7x7 " << sumv << " in " << v7x7.size();
1461             for (unsigned int iv = 0; iv < v7x7.size(); iv++) {
1462               HcalDetId id = v7x7[iv].first;
1463               edm::LogVerbatim("IsoTrack") << " Cell " << iv << " 0x" << std::hex << v7x7[iv].first() << std::dec << " "
1464                                            << id << " Energy " << v7x7[iv].second;
1465             }
1466           }
1467         }
1468         if (doMC_) {
1469           trkHcalEne = spr::eCaloSimInfo(iEvent, geo, pcalohh, SimTk, SimVtx, pTrack, *associate);
1470         }
1471
1472         // ====================================================================================================
1473         // get diff between track outermost hit position and the propagation point at outermost surface of tracker
1474         std::pair<math::XYZPoint, double> point2_TK0 = spr::propagateTrackerEnd(pTrack, bField, false);
1475         math::XYZPoint diff(pTrack->outerPosition().X() - point2_TK0.first.X(),
1476                             pTrack->outerPosition().Y() - point2_TK0.first.Y(),
1477                             pTrack->outerPosition().Z() - point2_TK0.first.Z());
1478         double trackOutPosOutHitDr = diff.R();
1479         double trackL = point2_TK0.second;
1480         if (myverbose_ == 5) {
1481           edm::LogVerbatim("IsoTrack") << " propagted " << point2_TK0.first << " " << point2_TK0.first.eta() << " "
1482                                        << point2_TK0.first.phi();
1483           edm::LogVerbatim("IsoTrack") << " outerPosition() " << pTrack->outerPosition() << " "
1484                                        << pTrack->outerPosition().eta() << " " << pTrack->outerPosition().phi();
1485           edm::LogVerbatim("IsoTrack") << "diff " << diff << " diffR " << diff.R() << " diffR/L "
1486                                        << diff.R() / point2_TK0.second;
1487         }
1488
1489         for (unsigned int ind = 0; ind < recVtxs->size(); ind++) {
1490           if (!((*recVtxs)[ind].isFake())) {
1491             reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
1492             if (reco::deltaR(eta1, phi1, (*vtxTrack)->eta(), (*vtxTrack)->phi()) < 0.01)
1493               t_trackPVIdx.push_back(ind);
1494             else
1495               t_trackPVIdx.push_back(-1);
1496           }
1497         }
1498
1499         // Fill the tree Branches here
1500         t_trackPVIdx.push_back(pVtxTkId);
1501         t_trackP.push_back(p1);
1502         t_trackPt.push_back(pt1);
1503         t_trackEta.push_back(eta1);
1504         t_trackPhi.push_back(phi1);
1505         t_trackEcalEta.push_back(etaEcal1);
1506         t_trackEcalPhi.push_back(phiEcal1);
1507         t_trackHcalEta.push_back(etaHcal1);
1508         t_trackHcalPhi.push_back(phiHcal1);
1509         t_trackDxy.push_back(dxy1);
1510         t_trackDz.push_back(dz1);
1511         t_trackDxyBS.push_back(dxybs1);
1512         t_trackDzBS.push_back(dzbs1);
1513         t_trackDxyPV.push_back(dxypv1);
1514         t_trackDzPV.push_back(dzpv1);
1515         t_trackChiSq.push_back(chisq1);
1516         t_trackNOuterHits.push_back(nOuterHits);
1517         t_NLayersCrossed.push_back(nLayersCrossed);
1518
1519         t_trackHitsTOB.push_back(hitp.stripTOBLayersWithMeasurement());
1520         t_trackHitsTEC.push_back(hitp.stripTECLayersWithMeasurement());
1521         t_trackHitInMissTOB.push_back(hitp.stripTOBLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1522         t_trackHitInMissTEC.push_back(hitp.stripTECLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1523         t_trackHitInMissTIB.push_back(hitp.stripTIBLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1524         t_trackHitInMissTID.push_back(hitp.stripTIDLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1525         t_trackHitInMissTIBTID.push_back(hitp.stripTIBLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS) +
1526                                          hitp.stripTIDLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1527
1528         t_trackHitOutMissTOB.push_back(hitp.stripTOBLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1529         t_trackHitOutMissTEC.push_back(hitp.stripTECLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1530         t_trackHitOutMissTIB.push_back(hitp.stripTIBLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1531         t_trackHitOutMissTID.push_back(hitp.stripTIDLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1532         t_trackHitOutMissTOBTEC.push_back(hitp.stripTOBLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS) +
1533                                           hitp.stripTECLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1534
1535         t_trackHitInMeasTOB.push_back(hitp.stripTOBLayersWithMeasurement());
1536         t_trackHitInMeasTEC.push_back(hitp.stripTECLayersWithMeasurement());
1537         t_trackHitInMeasTIB.push_back(hitp.stripTIBLayersWithMeasurement());
1538         t_trackHitInMeasTID.push_back(hitp.stripTIDLayersWithMeasurement());
1539         t_trackHitOutMeasTOB.push_back(hitp.stripTOBLayersWithMeasurement());
1540         t_trackHitOutMeasTEC.push_back(hitp.stripTECLayersWithMeasurement());
1541         t_trackHitOutMeasTIB.push_back(hitp.stripTIBLayersWithMeasurement());
1542         t_trackHitOutMeasTID.push_back(hitp.stripTIDLayersWithMeasurement());
1543         t_trackOutPosOutHitDr.push_back(trackOutPosOutHitDr);
1544         t_trackL.push_back(trackL);
1545
1546         t_maxNearP31x31.push_back(maxNearP31x31);
1547         t_maxNearP21x21.push_back(maxNearP21x21);
1548
1549         t_ecalSpike11x11.push_back(e11x11P.second);
1550         t_e7x7.push_back(e7x7P.first);
1551         t_e9x9.push_back(e9x9P.first);
1552         t_e11x11.push_back(e11x11P.first);
1553         t_e15x15.push_back(e15x15P.first);
1554
1555         t_e7x7_10Sig.push_back(e7x7_10SigP.first);
1556         t_e9x9_10Sig.push_back(e9x9_10SigP.first);
1557         t_e11x11_10Sig.push_back(e11x11_10SigP.first);
1558         t_e15x15_10Sig.push_back(e15x15_10SigP.first);
1559         t_e7x7_15Sig.push_back(e7x7_15SigP.first);
1560         t_e9x9_15Sig.push_back(e9x9_15SigP.first);
1561         t_e11x11_15Sig.push_back(e11x11_15SigP.first);
1562         t_e15x15_15Sig.push_back(e15x15_15SigP.first);
1563         t_e7x7_20Sig.push_back(e7x7_20SigP.first);
1564         t_e9x9_20Sig.push_back(e9x9_20SigP.first);
1565         t_e11x11_20Sig.push_back(e11x11_20SigP.first);
1566         t_e15x15_20Sig.push_back(e15x15_20SigP.first);
1567         t_e7x7_25Sig.push_back(e7x7_25SigP.first);
1568         t_e9x9_25Sig.push_back(e9x9_25SigP.first);
1569         t_e11x11_25Sig.push_back(e11x11_25SigP.first);
1570         t_e15x15_25Sig.push_back(e15x15_25SigP.first);
1571         t_e7x7_30Sig.push_back(e7x7_30SigP.first);
1572         t_e9x9_30Sig.push_back(e9x9_30SigP.first);
1573         t_e11x11_30Sig.push_back(e11x11_30SigP.first);
1574         t_e15x15_30Sig.push_back(e15x15_30SigP.first);
1575
1576         if (doMC_) {
1577           t_esim7x7.push_back(simInfo7x7.eTotal);
1578           t_esim9x9.push_back(simInfo9x9.eTotal);
1579           t_esim11x11.push_back(simInfo11x11.eTotal);
1580           t_esim15x15.push_back(simInfo15x15.eTotal);
1581
1582           t_esim7x7Matched.push_back(simInfo7x7.eMatched);
1583           t_esim9x9Matched.push_back(simInfo9x9.eMatched);
1584           t_esim11x11Matched.push_back(simInfo11x11.eMatched);
1585           t_esim15x15Matched.push_back(simInfo15x15.eMatched);
1586
1587           t_esim7x7Rest.push_back(simInfo7x7.eRest);
1588           t_esim9x9Rest.push_back(simInfo9x9.eRest);
1589           t_esim11x11Rest.push_back(simInfo11x11.eRest);
1590           t_esim15x15Rest.push_back(simInfo15x15.eRest);
1591
1592           t_esim7x7Photon.push_back(simInfo7x7.eGamma);
1593           t_esim9x9Photon.push_back(simInfo9x9.eGamma);
1594           t_esim11x11Photon.push_back(simInfo11x11.eGamma);
1595           t_esim15x15Photon.push_back(simInfo15x15.eGamma);
1596
1597           t_esim7x7NeutHad.push_back(simInfo7x7.eNeutralHad);
1598           t_esim9x9NeutHad.push_back(simInfo9x9.eNeutralHad);
1599           t_esim11x11NeutHad.push_back(simInfo11x11.eNeutralHad);
1600           t_esim15x15NeutHad.push_back(simInfo15x15.eNeutralHad);
1601
1602           t_esim7x7CharHad.push_back(simInfo7x7.eChargedHad);
1603           t_esim9x9CharHad.push_back(simInfo9x9.eChargedHad);
1604           t_esim11x11CharHad.push_back(simInfo11x11.eChargedHad);
1605           t_esim15x15CharHad.push_back(simInfo15x15.eChargedHad);
1606
1607           t_trkEcalEne.push_back(trkEcalEne);
1608           t_simTrackP.push_back(simTrackP);
1609           t_esimPdgId.push_back(simInfo11x11.pdgMatched);
1610         }
1611
1612         t_maxNearHcalP3x3.push_back(maxNearHcalP3x3);
1613         t_maxNearHcalP5x5.push_back(maxNearHcalP5x5);
1614         t_maxNearHcalP7x7.push_back(maxNearHcalP7x7);
1615
1616         t_h3x3.push_back(h3x3);
1617         t_h5x5.push_back(h5x5);
1618         t_h7x7.push_back(h7x7);
1619         t_h3x3Sig.push_back(h3x3Sig);
1620         t_h5x5Sig.push_back(h5x5Sig);
1621         t_h7x7Sig.push_back(h7x7Sig);
1622
1623         t_infoHcal.push_back(trkDetItr->okHCAL);
1624         if (doMC_) {
1625           t_trkHcalEne.push_back(hcalScale * trkHcalEne);
1626
1627           t_hsim3x3.push_back(hcalScale * hsimInfo3x3.eTotal);
1628           t_hsim5x5.push_back(hcalScale * hsimInfo5x5.eTotal);
1629           t_hsim7x7.push_back(hcalScale * hsimInfo7x7.eTotal);
1630
1631           t_hsim3x3Matched.push_back(hcalScale * hsimInfo3x3.eMatched);
1632           t_hsim5x5Matched.push_back(hcalScale * hsimInfo5x5.eMatched);
1633           t_hsim7x7Matched.push_back(hcalScale * hsimInfo7x7.eMatched);
1634
1635           t_hsim3x3Rest.push_back(hcalScale * hsimInfo3x3.eRest);
1636           t_hsim5x5Rest.push_back(hcalScale * hsimInfo5x5.eRest);
1637           t_hsim7x7Rest.push_back(hcalScale * hsimInfo7x7.eRest);
1638
1639           t_hsim3x3Photon.push_back(hcalScale * hsimInfo3x3.eGamma);
1640           t_hsim5x5Photon.push_back(hcalScale * hsimInfo5x5.eGamma);
1641           t_hsim7x7Photon.push_back(hcalScale * hsimInfo7x7.eGamma);
1642
1643           t_hsim3x3NeutHad.push_back(hcalScale * hsimInfo3x3.eNeutralHad);
1644           t_hsim5x5NeutHad.push_back(hcalScale * hsimInfo5x5.eNeutralHad);
1645           t_hsim7x7NeutHad.push_back(hcalScale * hsimInfo7x7.eNeutralHad);
1646
1647           t_hsim3x3CharHad.push_back(hcalScale * hsimInfo3x3.eChargedHad);
1648           t_hsim5x5CharHad.push_back(hcalScale * hsimInfo5x5.eChargedHad);
1649           t_hsim7x7CharHad.push_back(hcalScale * hsimInfo7x7.eChargedHad);
1650         }
1651
1652       }  // if loosely isolated track
1653     }  // check p1/eta
1654   }  // loop over track collection
1655
1656   if (haveIsoTrack)
1657     tree_->Fill();
1658 }
1659
1660 // ----- method called once each job just before starting event loop ----
1661 void IsolatedTracksNxN::beginJob() {
1662   nEventProc_ = 0;
1663   nbad_ = 0;
1664   initL1_ = false;
1665   double tempgen_TH[NPBins + 1] = {
1666       0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 9.0, 11.0, 15.0, 20.0, 30.0, 50.0, 75.0, 100.0};
1667
1668   for (unsigned int i = 0; i < NPBins + 1; i++)
1669     genPartPBins[i] = tempgen_TH[i];
1670
1671   double tempgen_Eta[NEtaBins + 1] = {0.0, 1.131, 1.653, 2.172};
1672
1673   for (unsigned int i = 0; i < NEtaBins + 1; i++)
1674     genPartEtaBins[i] = tempgen_Eta[i];
1675
1676   bookHistograms();
1677 }
1678
1679 // ----- method called once each job just after ending the event loop ----
1680 void IsolatedTracksNxN::endJob() {
1681   if (L1TriggerAlgoInfo_) {
1682     std::map<std::pair<unsigned int, std::string>, int>::iterator itr;
1683     for (itr = l1AlgoMap_.begin(); itr != l1AlgoMap_.end(); itr++) {
1684       edm::LogVerbatim("IsoTrack") << " ****endjob**** " << (itr->first).first << " " << (itr->first).second << " "
1685                                    << itr->second;
1686       int ibin = (itr->first).first;
1687       TString name((itr->first).second);
1688       h_L1AlgoNames->GetXaxis()->SetBinLabel(ibin + 1, name);
1689     }
1690     edm::LogVerbatim("IsoTrack") << "Number of Events Processed " << nEventProc_;
1691   }
1692 }
1693
1694 //========================================================================================================
1695
1696 void IsolatedTracksNxN::clearTreeVectors() {
1697   t_PVx.clear();
1698   t_PVy.clear();
1699   t_PVz.clear();
1700   t_PVisValid.clear();
1701   t_PVndof.clear();
1702   t_PVNTracks.clear();
1703   t_PVNTracksWt.clear();
1704   t_PVTracksSumPt.clear();
1705   t_PVTracksSumPtWt.clear();
1706   t_PVNTracksHP.clear();
1707   t_PVNTracksHPWt.clear();
1708   t_PVTracksSumPtHP.clear();
1709   t_PVTracksSumPtHPWt.clear();
1710
1711   for (int i = 0; i < 128; i++)
1712     t_L1Decision[i] = 0;
1713   t_L1AlgoNames.clear();
1714   t_L1PreScale.clear();
1715
1716   t_L1CenJetPt.clear();
1717   t_L1CenJetEta.clear();
1718   t_L1CenJetPhi.clear();
1719   t_L1FwdJetPt.clear();
1720   t_L1FwdJetEta.clear();
1721   t_L1FwdJetPhi.clear();
1722   t_L1TauJetPt.clear();
1723   t_L1TauJetEta.clear();
1724   t_L1TauJetPhi.clear();
1725   t_L1MuonPt.clear();
1726   t_L1MuonEta.clear();
1727   t_L1MuonPhi.clear();
1728   t_L1IsoEMPt.clear();
1729   t_L1IsoEMEta.clear();
1730   t_L1IsoEMPhi.clear();
1731   t_L1NonIsoEMPt.clear();
1732   t_L1NonIsoEMEta.clear();
1733   t_L1NonIsoEMPhi.clear();
1734   t_L1METPt.clear();
1735   t_L1METEta.clear();
1736   t_L1METPhi.clear();
1737
1738   t_jetPt.clear();
1739   t_jetEta.clear();
1740   t_jetPhi.clear();
1741   t_nTrksJetCalo.clear();
1742   t_nTrksJetVtx.clear();
1743
1744   t_trackPAll.clear();
1745   t_trackEtaAll.clear();
1746   t_trackPhiAll.clear();
1747   t_trackPdgIdAll.clear();
1748   t_trackPtAll.clear();
1749   t_trackDxyAll.clear();
1750   t_trackDzAll.clear();
1751   t_trackDxyPVAll.clear();
1752   t_trackDzPVAll.clear();
1753   t_trackChiSqAll.clear();
1754
1755   t_trackP.clear();
1756   t_trackPt.clear();
1757   t_trackEta.clear();
1758   t_trackPhi.clear();
1759   t_trackEcalEta.clear();
1760   t_trackEcalPhi.clear();
1761   t_trackHcalEta.clear();
1762   t_trackHcalPhi.clear();
1763   t_NLayersCrossed.clear();
1764   t_trackNOuterHits.clear();
1765   t_trackDxy.clear();
1766   t_trackDxyBS.clear();
1767   t_trackDz.clear();
1768   t_trackDzBS.clear();
1769   t_trackDxyPV.clear();
1770   t_trackDzPV.clear();
1771   t_trackChiSq.clear();
1772   t_trackPVIdx.clear();
1773   t_trackHitsTOB.clear();
1774   t_trackHitsTEC.clear();
1775   t_trackHitInMissTOB.clear();
1776   t_trackHitInMissTEC.clear();
1777   t_trackHitInMissTIB.clear();
1778   t_trackHitInMissTID.clear();
1779   t_trackHitInMissTIBTID.clear();
1780   t_trackHitOutMissTOB.clear();
1781   t_trackHitOutMissTEC.clear();
1782   t_trackHitOutMissTIB.clear();
1783   t_trackHitOutMissTID.clear();
1784   t_trackHitOutMissTOBTEC.clear();
1785   t_trackHitInMeasTOB.clear();
1786   t_trackHitInMeasTEC.clear();
1787   t_trackHitInMeasTIB.clear();
1788   t_trackHitInMeasTID.clear();
1789   t_trackHitOutMeasTOB.clear();
1790   t_trackHitOutMeasTEC.clear();
1791   t_trackHitOutMeasTIB.clear();
1792   t_trackHitOutMeasTID.clear();
1793   t_trackOutPosOutHitDr.clear();
1794   t_trackL.clear();
1795
1796   t_maxNearP31x31.clear();
1797   t_maxNearP21x21.clear();
1798
1799   t_ecalSpike11x11.clear();
1800   t_e7x7.clear();
1801   t_e9x9.clear();
1802   t_e11x11.clear();
1803   t_e15x15.clear();
1804
1805   t_e7x7_10Sig.clear();
1806   t_e9x9_10Sig.clear();
1807   t_e11x11_10Sig.clear();
1808   t_e15x15_10Sig.clear();
1809   t_e7x7_15Sig.clear();
1810   t_e9x9_15Sig.clear();
1811   t_e11x11_15Sig.clear();
1812   t_e15x15_15Sig.clear();
1813   t_e7x7_20Sig.clear();
1814   t_e9x9_20Sig.clear();
1815   t_e11x11_20Sig.clear();
1816   t_e15x15_20Sig.clear();
1817   t_e7x7_25Sig.clear();
1818   t_e9x9_25Sig.clear();
1819   t_e11x11_25Sig.clear();
1820   t_e15x15_25Sig.clear();
1821   t_e7x7_30Sig.clear();
1822   t_e9x9_30Sig.clear();
1823   t_e11x11_30Sig.clear();
1824   t_e15x15_30Sig.clear();
1825
1826   if (doMC_) {
1827     t_simTrackP.clear();
1828     t_esimPdgId.clear();
1829     t_trkEcalEne.clear();
1830
1831     t_esim7x7.clear();
1832     t_esim9x9.clear();
1833     t_esim11x11.clear();
1834     t_esim15x15.clear();
1835
1836     t_esim7x7Matched.clear();
1837     t_esim9x9Matched.clear();
1838     t_esim11x11Matched.clear();
1839     t_esim15x15Matched.clear();
1840
1841     t_esim7x7Rest.clear();
1842     t_esim9x9Rest.clear();
1843     t_esim11x11Rest.clear();
1844     t_esim15x15Rest.clear();
1845
1846     t_esim7x7Photon.clear();
1847     t_esim9x9Photon.clear();
1848     t_esim11x11Photon.clear();
1849     t_esim15x15Photon.clear();
1850
1851     t_esim7x7NeutHad.clear();
1852     t_esim9x9NeutHad.clear();
1853     t_esim11x11NeutHad.clear();
1854     t_esim15x15NeutHad.clear();
1855
1856     t_esim7x7CharHad.clear();
1857     t_esim9x9CharHad.clear();
1858     t_esim11x11CharHad.clear();
1859     t_esim15x15CharHad.clear();
1860   }
1861
1862   t_maxNearHcalP3x3.clear();
1863   t_maxNearHcalP5x5.clear();
1864   t_maxNearHcalP7x7.clear();
1865
1866   t_h3x3.clear();
1867   t_h5x5.clear();
1868   t_h7x7.clear();
1869   t_h3x3Sig.clear();
1870   t_h5x5Sig.clear();
1871   t_h7x7Sig.clear();
1872
1873   t_infoHcal.clear();
1874
1875   if (doMC_) {
1876     t_trkHcalEne.clear();
1877
1878     t_hsim3x3.clear();
1879     t_hsim5x5.clear();
1880     t_hsim7x7.clear();
1881     t_hsim3x3Matched.clear();
1882     t_hsim5x5Matched.clear();
1883     t_hsim7x7Matched.clear();
1884     t_hsim3x3Rest.clear();
1885     t_hsim5x5Rest.clear();
1886     t_hsim7x7Rest.clear();
1887     t_hsim3x3Photon.clear();
1888     t_hsim5x5Photon.clear();
1889     t_hsim7x7Photon.clear();
1890     t_hsim3x3NeutHad.clear();
1891     t_hsim5x5NeutHad.clear();
1892     t_hsim7x7NeutHad.clear();
1893     t_hsim3x3CharHad.clear();
1894     t_hsim5x5CharHad.clear();
1895     t_hsim7x7CharHad.clear();
1896   }
1897 }
1898
1899 void IsolatedTracksNxN::bookHistograms() {
1900   char hname[100], htit[100];
1901
1902   edm::Service<TFileService> fs;
1903   TFileDirectory dir = fs->mkdir("nearMaxTrackP");
1904
1905   for (unsigned int ieta = 0; ieta < NEtaBins; ieta++) {
1906     double lowEta = -5.0, highEta = 5.0;
1907     lowEta = genPartEtaBins[ieta];
1908     highEta = genPartEtaBins[ieta + 1];
1909
1910     for (unsigned int ipt = 0; ipt < NPBins; ipt++) {
1911       double lowP = 0.0, highP = 300.0;
1912       lowP = genPartPBins[ipt];
1913       highP = genPartPBins[ipt + 1];
1914       sprintf(hname, "h_maxNearP31x31_ptBin%i_etaBin%i", ipt, ieta);
1915       sprintf(htit, "maxNearP in 31x31 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP);
1916       h_maxNearP31x31[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
1917       h_maxNearP31x31[ipt][ieta]->Sumw2();
1918       sprintf(hname, "h_maxNearP25x25_ptBin%i_etaBin%i", ipt, ieta);
1919       sprintf(htit, "maxNearP in 25x25 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP);
1920       h_maxNearP25x25[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
1921       h_maxNearP25x25[ipt][ieta]->Sumw2();
1922       sprintf(hname, "h_maxNearP21x21_ptBin%i_etaBin%i", ipt, ieta);
1923       sprintf(htit, "maxNearP in 21x21 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP);
1924       h_maxNearP21x21[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
1925       h_maxNearP21x21[ipt][ieta]->Sumw2();
1926       sprintf(hname, "h_maxNearP15x15_ptBin%i_etaBin%i", ipt, ieta);
1927       sprintf(htit, "maxNearP in 15x15 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP);
1928       h_maxNearP15x15[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
1929       h_maxNearP15x15[ipt][ieta]->Sumw2();
1930     }
1931   }
1932
1933   h_L1AlgoNames = fs->make<TH1I>("h_L1AlgoNames", "h_L1AlgoNames:Bin Labels", 128, -0.5, 127.5);
1934
1935   // Reconstructed Tracks
1936
1937   h_PVTracksWt = fs->make<TH1F>("h_PVTracksWt", "h_PVTracksWt", 600, -0.1, 1.1);
1938
1939   h_nTracks = fs->make<TH1F>("h_nTracks", "h_nTracks", 1000, -0.5, 999.5);
1940
1941   sprintf(hname, "h_recEtaPt_0");
1942   sprintf(htit, "h_recEtaPt (all tracks Eta vs pT)");
1943   h_recEtaPt_0 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
1944
1945   sprintf(hname, "h_recEtaP_0");
1946   sprintf(htit, "h_recEtaP (all tracks Eta vs pT)");
1947   h_recEtaP_0 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
1948
1949   h_recPt_0 = fs->make<TH1F>("h_recPt_0", "Pt (all tracks)", 15, genPartPBins);
1950   h_recP_0 = fs->make<TH1F>("h_recP_0", "P  (all tracks)", 15, genPartPBins);
1951   h_recEta_0 = fs->make<TH1F>("h_recEta_0", "Eta (all tracks)", 60, -3.0, 3.0);
1952   h_recPhi_0 = fs->make<TH1F>("h_recPhi_0", "Phi (all tracks)", 100, -3.2, 3.2);
1953   //-------------------------
1954   sprintf(hname, "h_recEtaPt_1");
1955   sprintf(htit, "h_recEtaPt (all good tracks Eta vs pT)");
1956   h_recEtaPt_1 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
1957
1958   sprintf(hname, "h_recEtaP_1");
1959   sprintf(htit, "h_recEtaP (all good tracks Eta vs pT)");
1960   h_recEtaP_1 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
1961
1962   h_recPt_1 = fs->make<TH1F>("h_recPt_1", "Pt (all good tracks)", 15, genPartPBins);
1963   h_recP_1 = fs->make<TH1F>("h_recP_1", "P  (all good tracks)", 15, genPartPBins);
1964   h_recEta_1 = fs->make<TH1F>("h_recEta_1", "Eta (all good tracks)", 60, -3.0, 3.0);
1965   h_recPhi_1 = fs->make<TH1F>("h_recPhi_1", "Phi (all good tracks)", 100, -3.2, 3.2);
1966   //-------------------------
1967   sprintf(hname, "h_recEtaPt_2");
1968   sprintf(htit, "h_recEtaPt (charge isolation Eta vs pT)");
1969   h_recEtaPt_2 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
1970
1971   sprintf(hname, "h_recEtaP_2");
1972   sprintf(htit, "h_recEtaP (charge isolation Eta vs pT)");
1973   h_recEtaP_2 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
1974
1975   h_recPt_2 = fs->make<TH1F>("h_recPt_2", "Pt (charge isolation)", 15, genPartPBins);
1976   h_recP_2 = fs->make<TH1F>("h_recP_2", "P  (charge isolation)", 15, genPartPBins);
1977   h_recEta_2 = fs->make<TH1F>("h_recEta_2", "Eta (charge isolation)", 60, -3.0, 3.0);
1978   h_recPhi_2 = fs->make<TH1F>("h_recPhi_2", "Phi (charge isolation)", 100, -3.2, 3.2);
1979
1980   tree_ = fs->make<TTree>("tree", "tree");
1981   tree_->SetAutoSave(10000);
1982
1983   tree_->Branch("t_EvtNo", &t_EvtNo, "t_EvtNo/I");
1984   tree_->Branch("t_RunNo", &t_RunNo, "t_RunNo/I");
1985   tree_->Branch("t_Lumi", &t_Lumi, "t_Lumi/I");
1986   tree_->Branch("t_Bunch", &t_Bunch, "t_Bunch/I");
1987
1988   //----- L1Trigger information
1989   for (int i = 0; i < 128; i++)
1990     t_L1Decision[i] = 0;
1991
1992   tree_->Branch("t_L1Decision", t_L1Decision, "t_L1Decision[128]/I");
1993   tree_->Branch("t_L1AlgoNames", &t_L1AlgoNames);
1994   tree_->Branch("t_L1PreScale", &t_L1PreScale);
1995   tree_->Branch("t_L1CenJetPt", &t_L1CenJetPt);
1996   tree_->Branch("t_L1CenJetEta", &t_L1CenJetEta);
1997   tree_->Branch("t_L1CenJetPhi", &t_L1CenJetPhi);
1998   tree_->Branch("t_L1FwdJetPt", &t_L1FwdJetPt);
1999   tree_->Branch("t_L1FwdJetEta", &t_L1FwdJetEta);
2000   tree_->Branch("t_L1FwdJetPhi", &t_L1FwdJetPhi);
2001   tree_->Branch("t_L1TauJetPt", &t_L1TauJetPt);
2002   tree_->Branch("t_L1TauJetEta", &t_L1TauJetEta);
2003   tree_->Branch("t_L1TauJetPhi", &t_L1TauJetPhi);
2004   tree_->Branch("t_L1MuonPt", &t_L1MuonPt);
2005   tree_->Branch("t_L1MuonEta", &t_L1MuonEta);
2006   tree_->Branch("t_L1MuonPhi", &t_L1MuonPhi);
2007   tree_->Branch("t_L1IsoEMPt", &t_L1IsoEMPt);
2008   tree_->Branch("t_L1IsoEMEta", &t_L1IsoEMEta);
2009   tree_->Branch("t_L1IsoEMPhi", &t_L1IsoEMPhi);
2010   tree_->Branch("t_L1NonIsoEMPt", &t_L1NonIsoEMPt);
2011   tree_->Branch("t_L1NonIsoEMEta", &t_L1NonIsoEMEta);
2012   tree_->Branch("t_L1NonIsoEMPhi", &t_L1NonIsoEMPhi);
2013   tree_->Branch("t_L1METPt", &t_L1METPt);
2014   tree_->Branch("t_L1METEta", &t_L1METEta);
2015   tree_->Branch("t_L1METPhi", &t_L1METPhi);
2016
2017   tree_->Branch("t_jetPt", &t_jetPt);
2018   tree_->Branch("t_jetEta", &t_jetEta);
2019   tree_->Branch("t_jetPhi", &t_jetPhi);
2020   tree_->Branch("t_nTrksJetCalo", &t_nTrksJetCalo);
2021   tree_->Branch("t_nTrksJetVtx", &t_nTrksJetVtx);
2022
2023   tree_->Branch("t_trackPAll", &t_trackPAll);
2024   tree_->Branch("t_trackPhiAll", &t_trackPhiAll);
2025   tree_->Branch("t_trackEtaAll", &t_trackEtaAll);
2026   tree_->Branch("t_trackPtAll", &t_trackPtAll);
2027   tree_->Branch("t_trackDxyAll", &t_trackDxyAll);
2028   tree_->Branch("t_trackDzAll", &t_trackDzAll);
2029   tree_->Branch("t_trackDxyPVAll", &t_trackDxyPVAll);
2030   tree_->Branch("t_trackDzPVAll", &t_trackDzPVAll);
2031   tree_->Branch("t_trackChiSqAll", &t_trackChiSqAll);
2032   //tree_->Branch("t_trackPdgIdAll",     &t_trackPdgIdAll);
2033
2034   tree_->Branch("t_trackP", &t_trackP);
2035   tree_->Branch("t_trackPt", &t_trackPt);
2036   tree_->Branch("t_trackEta", &t_trackEta);
2037   tree_->Branch("t_trackPhi", &t_trackPhi);
2038   tree_->Branch("t_trackEcalEta", &t_trackEcalEta);
2039   tree_->Branch("t_trackEcalPhi", &t_trackEcalPhi);
2040   tree_->Branch("t_trackHcalEta", &t_trackHcalEta);
2041   tree_->Branch("t_trackHcalPhi", &t_trackHcalPhi);
2042
2043   tree_->Branch("t_trackNOuterHits", &t_trackNOuterHits);
2044   tree_->Branch("t_NLayersCrossed", &t_NLayersCrossed);
2045   tree_->Branch("t_trackHitsTOB", &t_trackHitsTOB);
2046   tree_->Branch("t_trackHitsTEC", &t_trackHitsTEC);
2047   tree_->Branch("t_trackHitInMissTOB", &t_trackHitInMissTOB);
2048   tree_->Branch("t_trackHitInMissTEC", &t_trackHitInMissTEC);
2049   tree_->Branch("t_trackHitInMissTIB", &t_trackHitInMissTIB);
2050   tree_->Branch("t_trackHitInMissTID", &t_trackHitInMissTID);
2051   tree_->Branch("t_trackHitInMissTIBTID", &t_trackHitInMissTIBTID);
2052   tree_->Branch("t_trackHitOutMissTOB", &t_trackHitOutMissTOB);
2053   tree_->Branch("t_trackHitOutMissTEC", &t_trackHitOutMissTEC);
2054   tree_->Branch("t_trackHitOutMissTIB", &t_trackHitOutMissTIB);
2055   tree_->Branch("t_trackHitOutMissTID", &t_trackHitOutMissTID);
2056   tree_->Branch("t_trackHitOutMissTOBTEC", &t_trackHitOutMissTOBTEC);
2057   tree_->Branch("t_trackHitInMeasTOB", &t_trackHitInMeasTOB);
2058   tree_->Branch("t_trackHitInMeasTEC", &t_trackHitInMeasTEC);
2059   tree_->Branch("t_trackHitInMeasTIB", &t_trackHitInMeasTIB);
2060   tree_->Branch("t_trackHitInMeasTID", &t_trackHitInMeasTID);
2061   tree_->Branch("t_trackHitOutMeasTOB", &t_trackHitOutMeasTOB);
2062   tree_->Branch("t_trackHitOutMeasTEC", &t_trackHitOutMeasTEC);
2063   tree_->Branch("t_trackHitOutMeasTIB", &t_trackHitOutMeasTIB);
2064   tree_->Branch("t_trackHitOutMeasTID", &t_trackHitOutMeasTID);
2065   tree_->Branch("t_trackOutPosOutHitDr", &t_trackOutPosOutHitDr);
2066   tree_->Branch("t_trackL", &t_trackL);
2067
2068   tree_->Branch("t_trackDxy", &t_trackDxy);
2069   tree_->Branch("t_trackDxyBS", &t_trackDxyBS);
2070   tree_->Branch("t_trackDz", &t_trackDz);
2071   tree_->Branch("t_trackDzBS", &t_trackDzBS);
2072   tree_->Branch("t_trackDxyPV", &t_trackDxyPV);
2073   tree_->Branch("t_trackDzPV", &t_trackDzPV);
2074   tree_->Branch("t_trackChiSq", &t_trackChiSq);
2075   tree_->Branch("t_trackPVIdx", &t_trackPVIdx);
2076
2077   tree_->Branch("t_maxNearP31x31", &t_maxNearP31x31);
2078   tree_->Branch("t_maxNearP21x21", &t_maxNearP21x21);
2079
2080   tree_->Branch("t_ecalSpike11x11", &t_ecalSpike11x11);
2081   tree_->Branch("t_e7x7", &t_e7x7);
2082   tree_->Branch("t_e9x9", &t_e9x9);
2083   tree_->Branch("t_e11x11", &t_e11x11);
2084   tree_->Branch("t_e15x15", &t_e15x15);
2085
2086   tree_->Branch("t_e7x7_10Sig", &t_e7x7_10Sig);
2087   tree_->Branch("t_e9x9_10Sig", &t_e9x9_10Sig);
2088   tree_->Branch("t_e11x11_10Sig", &t_e11x11_10Sig);
2089   tree_->Branch("t_e15x15_10Sig", &t_e15x15_10Sig);
2090   tree_->Branch("t_e7x7_15Sig", &t_e7x7_15Sig);
2091   tree_->Branch("t_e9x9_15Sig", &t_e9x9_15Sig);
2092   tree_->Branch("t_e11x11_15Sig", &t_e11x11_15Sig);
2093   tree_->Branch("t_e15x15_15Sig", &t_e15x15_15Sig);
2094   tree_->Branch("t_e7x7_20Sig", &t_e7x7_20Sig);
2095   tree_->Branch("t_e9x9_20Sig", &t_e9x9_20Sig);
2096   tree_->Branch("t_e11x11_20Sig", &t_e11x11_20Sig);
2097   tree_->Branch("t_e15x15_20Sig", &t_e15x15_20Sig);
2098   tree_->Branch("t_e7x7_25Sig", &t_e7x7_25Sig);
2099   tree_->Branch("t_e9x9_25Sig", &t_e9x9_25Sig);
2100   tree_->Branch("t_e11x11_25Sig", &t_e11x11_25Sig);
2101   tree_->Branch("t_e15x15_25Sig", &t_e15x15_25Sig);
2102   tree_->Branch("t_e7x7_30Sig", &t_e7x7_30Sig);
2103   tree_->Branch("t_e9x9_30Sig", &t_e9x9_30Sig);
2104   tree_->Branch("t_e11x11_30Sig", &t_e11x11_30Sig);
2105   tree_->Branch("t_e15x15_30Sig", &t_e15x15_30Sig);
2106
2107   if (doMC_) {
2108     tree_->Branch("t_esim7x7", &t_esim7x7);
2109     tree_->Branch("t_esim9x9", &t_esim9x9);
2110     tree_->Branch("t_esim11x11", &t_esim11x11);
2111     tree_->Branch("t_esim15x15", &t_esim15x15);
2112
2113     tree_->Branch("t_esim7x7Matched", &t_esim7x7Matched);
2114     tree_->Branch("t_esim9x9Matched", &t_esim9x9Matched);
2115     tree_->Branch("t_esim11x11Matched", &t_esim11x11Matched);
2116     tree_->Branch("t_esim15x15Matched", &t_esim15x15Matched);
2117
2118     tree_->Branch("t_esim7x7Rest", &t_esim7x7Rest);
2119     tree_->Branch("t_esim9x9Rest", &t_esim9x9Rest);
2120     tree_->Branch("t_esim11x11Rest", &t_esim11x11Rest);
2121     tree_->Branch("t_esim15x15Rest", &t_esim15x15Rest);
2122
2123     tree_->Branch("t_esim7x7Photon", &t_esim7x7Photon);
2124     tree_->Branch("t_esim9x9Photon", &t_esim9x9Photon);
2125     tree_->Branch("t_esim11x11Photon", &t_esim11x11Photon);
2126     tree_->Branch("t_esim15x15Photon", &t_esim15x15Photon);
2127
2128     tree_->Branch("t_esim7x7NeutHad", &t_esim7x7NeutHad);
2129     tree_->Branch("t_esim9x9NeutHad", &t_esim9x9NeutHad);
2130     tree_->Branch("t_esim11x11NeutHad", &t_esim11x11NeutHad);
2131     tree_->Branch("t_esim15x15NeutHad", &t_esim15x15NeutHad);
2132
2133     tree_->Branch("t_esim7x7CharHad", &t_esim7x7CharHad);
2134     tree_->Branch("t_esim9x9CharHad", &t_esim9x9CharHad);
2135     tree_->Branch("t_esim11x11CharHad", &t_esim11x11CharHad);
2136     tree_->Branch("t_esim15x15CharHad", &t_esim15x15CharHad);
2137
2138     tree_->Branch("t_trkEcalEne", &t_trkEcalEne);
2139     tree_->Branch("t_simTrackP", &t_simTrackP);
2140     tree_->Branch("t_esimPdgId", &t_esimPdgId);
2141   }
2142
2143   tree_->Branch("t_maxNearHcalP3x3", &t_maxNearHcalP3x3);
2144   tree_->Branch("t_maxNearHcalP5x5", &t_maxNearHcalP5x5);
2145   tree_->Branch("t_maxNearHcalP7x7", &t_maxNearHcalP7x7);
2146   tree_->Branch("t_h3x3", &t_h3x3);
2147   tree_->Branch("t_h5x5", &t_h5x5);
2148   tree_->Branch("t_h7x7", &t_h7x7);
2149   tree_->Branch("t_h3x3Sig", &t_h3x3Sig);
2150   tree_->Branch("t_h5x5Sig", &t_h5x5Sig);
2151   tree_->Branch("t_h7x7Sig", &t_h7x7Sig);
2152   tree_->Branch("t_infoHcal", &t_infoHcal);
2153
2154   if (doMC_) {
2155     tree_->Branch("t_trkHcalEne", &t_trkHcalEne);
2156     tree_->Branch("t_hsim3x3", &t_hsim3x3);
2157     tree_->Branch("t_hsim5x5", &t_hsim5x5);
2158     tree_->Branch("t_hsim7x7", &t_hsim7x7);
2159     tree_->Branch("t_hsim3x3Matched", &t_hsim3x3Matched);
2160     tree_->Branch("t_hsim5x5Matched", &t_hsim5x5Matched);
2161     tree_->Branch("t_hsim7x7Matched", &t_hsim7x7Matched);
2162     tree_->Branch("t_hsim3x3Rest", &t_hsim3x3Rest);
2163     tree_->Branch("t_hsim5x5Rest", &t_hsim5x5Rest);
2164     tree_->Branch("t_hsim7x7Rest", &t_hsim7x7Rest);
2165     tree_->Branch("t_hsim3x3Photon", &t_hsim3x3Photon);
2166     tree_->Branch("t_hsim5x5Photon", &t_hsim5x5Photon);
2167     tree_->Branch("t_hsim7x7Photon", &t_hsim7x7Photon);
2168     tree_->Branch("t_hsim3x3NeutHad", &t_hsim3x3NeutHad);
2169     tree_->Branch("t_hsim5x5NeutHad", &t_hsim5x5NeutHad);
2170     tree_->Branch("t_hsim7x7NeutHad", &t_hsim7x7NeutHad);
2171     tree_->Branch("t_hsim3x3CharHad", &t_hsim3x3CharHad);
2172     tree_->Branch("t_hsim5x5CharHad", &t_hsim5x5CharHad);
2173     tree_->Branch("t_hsim7x7CharHad", &t_hsim7x7CharHad);
2174   }
2175   tree_->Branch("t_nTracks", &t_nTracks, "t_nTracks/I");
2176 }
2177
2178 void IsolatedTracksNxN::printTrack(const reco::Track *pTrack) {
2179   std::string theTrackQuality = "highPurity";
2180   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
2181
2182   edm::LogVerbatim("IsoTrack") << " Reference Point " << pTrack->referencePoint() << "\n TrackMmentum "
2183                                << pTrack->momentum() << " (pt,eta,phi)(" << pTrack->pt() << "," << pTrack->eta() << ","
2184                                << pTrack->phi() << ")"
2185                                << " p " << pTrack->p() << "\n Normalized chi2 " << pTrack->normalizedChi2()
2186                                << "  charge " << pTrack->charge() << " qoverp() " << pTrack->qoverp() << "+-"
2187                                << pTrack->qoverpError() << " d0 " << pTrack->d0() << "\n NValidHits "
2188                                << pTrack->numberOfValidHits() << "  NLostHits " << pTrack->numberOfLostHits()
2189                                << " TrackQuality " << pTrack->qualityName(trackQuality_) << " "
2190                                << pTrack->quality(trackQuality_);
2191
2192   if (printTrkHitPattern_) {
2193     const reco::HitPattern &p = pTrack->hitPattern();
2194
2195     std::ostringstream st1;
2196     st1 << "default ";
2197     for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::TRACK_HITS); i++) {
2198       p.printHitPattern(reco::HitPattern::TRACK_HITS, i, st1);
2199     }
2200     edm::LogVerbatim("IsoTrack") << st1.str();
2201     std::ostringstream st2;
2202     st2 << "trackerMissingInnerHits() ";
2203     for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::MISSING_INNER_HITS); i++) {
2204       p.printHitPattern(reco::HitPattern::MISSING_INNER_HITS, i, st2);
2205     }
2206     edm::LogVerbatim("IsoTrack") << st2.str();
2207     std::ostringstream st3;
2208     st3 << "trackerMissingOuterHits() ";
2209     for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::MISSING_OUTER_HITS); i++) {
2210       p.printHitPattern(reco::HitPattern::MISSING_OUTER_HITS, i, st3);
2211     }
2212     edm::LogVerbatim("IsoTrack") << st3.str();
2213
2214     edm::LogVerbatim("IsoTrack") << "\n \t trackerLayersWithMeasurement() " << p.trackerLayersWithMeasurement()
2215                                  << "\n \t pixelLayersWithMeasurement() " << p.pixelLayersWithMeasurement()
2216                                  << "\n \t stripLayersWithMeasurement() " << p.stripLayersWithMeasurement()
2217                                  << "\n \t pixelBarrelLayersWithMeasurement() " << p.pixelBarrelLayersWithMeasurement()
2218                                  << "\n \t pixelEndcapLayersWithMeasurement() " << p.pixelEndcapLayersWithMeasurement()
2219                                  << "\n \t stripTIBLayersWithMeasurement() " << p.stripTIBLayersWithMeasurement()
2220                                  << "\n \t stripTIDLayersWithMeasurement() " << p.stripTIDLayersWithMeasurement()
2221                                  << "\n \t stripTOBLayersWithMeasurement() " << p.stripTOBLayersWithMeasurement()
2222                                  << "\n \t stripTECLayersWithMeasurement() " << p.stripTECLayersWithMeasurement();
2223   }
2224 }
2225
2226 //define this as a plug-in
2227 DEFINE_FWK_MODULE(IsolatedTracksNxN);