Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-12 22:32:29

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       t_L1AlgoNames(nullptr),
0324       t_L1PreScale(nullptr),
0325       t_L1CenJetPt(nullptr),
0326       t_L1CenJetEta(nullptr),
0327       t_L1CenJetPhi(nullptr),
0328       t_L1FwdJetPt(nullptr),
0329       t_L1FwdJetEta(nullptr),
0330       t_L1FwdJetPhi(nullptr),
0331       t_L1TauJetPt(nullptr),
0332       t_L1TauJetEta(nullptr),
0333       t_L1TauJetPhi(nullptr),
0334       t_L1MuonPt(nullptr),
0335       t_L1MuonEta(nullptr),
0336       t_L1MuonPhi(nullptr),
0337       t_L1IsoEMPt(nullptr),
0338       t_L1IsoEMEta(nullptr),
0339       t_L1IsoEMPhi(nullptr),
0340       t_L1NonIsoEMPt(nullptr),
0341       t_L1NonIsoEMEta(nullptr),
0342       t_L1NonIsoEMPhi(nullptr),
0343       t_L1METPt(nullptr),
0344       t_L1METEta(nullptr),
0345       t_L1METPhi(nullptr),
0346       t_PVx(nullptr),
0347       t_PVy(nullptr),
0348       t_PVz(nullptr),
0349       t_PVTracksSumPt(nullptr),
0350       t_PVTracksSumPtWt(nullptr),
0351       t_PVTracksSumPtHP(nullptr),
0352       t_PVTracksSumPtHPWt(nullptr),
0353       t_PVisValid(nullptr),
0354       t_PVNTracks(nullptr),
0355       t_PVNTracksWt(nullptr),
0356       t_PVndof(nullptr),
0357       t_PVNTracksHP(nullptr),
0358       t_PVNTracksHPWt(nullptr),
0359       t_jetPt(nullptr),
0360       t_jetEta(nullptr),
0361       t_jetPhi(nullptr),
0362       t_nTrksJetCalo(nullptr),
0363       t_nTrksJetVtx(nullptr),
0364       t_trackPAll(nullptr),
0365       t_trackEtaAll(nullptr),
0366       t_trackPhiAll(nullptr),
0367       t_trackPdgIdAll(nullptr),
0368       t_trackPtAll(nullptr),
0369       t_trackDxyAll(nullptr),
0370       t_trackDzAll(nullptr),
0371       t_trackDxyPVAll(nullptr),
0372       t_trackDzPVAll(nullptr),
0373       t_trackChiSqAll(nullptr),
0374       t_trackP(nullptr),
0375       t_trackPt(nullptr),
0376       t_trackEta(nullptr),
0377       t_trackPhi(nullptr),
0378       t_trackEcalEta(nullptr),
0379       t_trackEcalPhi(nullptr),
0380       t_trackHcalEta(nullptr),
0381       t_trackHcalPhi(nullptr),
0382       t_trackDxy(nullptr),
0383       t_trackDxyBS(nullptr),
0384       t_trackDz(nullptr),
0385       t_trackDzBS(nullptr),
0386       t_trackDxyPV(nullptr),
0387       t_trackDzPV(nullptr),
0388       t_trackChiSq(nullptr),
0389       t_trackPVIdx(nullptr),
0390       t_NLayersCrossed(nullptr),
0391       t_trackNOuterHits(nullptr),
0392       t_trackHitsTOB(nullptr),
0393       t_trackHitsTEC(nullptr),
0394       t_trackHitInMissTOB(nullptr),
0395       t_trackHitInMissTEC(nullptr),
0396       t_trackHitInMissTIB(nullptr),
0397       t_trackHitInMissTID(nullptr),
0398       t_trackHitInMissTIBTID(nullptr),
0399       t_trackHitOutMissTOB(nullptr),
0400       t_trackHitOutMissTEC(nullptr),
0401       t_trackHitOutMissTIB(nullptr),
0402       t_trackHitOutMissTID(nullptr),
0403       t_trackHitOutMissTOBTEC(nullptr),
0404       t_trackHitInMeasTOB(nullptr),
0405       t_trackHitInMeasTEC(nullptr),
0406       t_trackHitInMeasTIB(nullptr),
0407       t_trackHitInMeasTID(nullptr),
0408       t_trackHitOutMeasTOB(nullptr),
0409       t_trackHitOutMeasTEC(nullptr),
0410       t_trackHitOutMeasTIB(nullptr),
0411       t_trackHitOutMeasTID(nullptr),
0412       t_trackOutPosOutHitDr(nullptr),
0413       t_trackL(nullptr),
0414       t_maxNearP31x31(nullptr),
0415       t_maxNearP21x21(nullptr),
0416       t_ecalSpike11x11(nullptr),
0417       t_e7x7(nullptr),
0418       t_e9x9(nullptr),
0419       t_e11x11(nullptr),
0420       t_e15x15(nullptr),
0421       t_e7x7_10Sig(nullptr),
0422       t_e9x9_10Sig(nullptr),
0423       t_e11x11_10Sig(nullptr),
0424       t_e15x15_10Sig(nullptr),
0425       t_e7x7_15Sig(nullptr),
0426       t_e9x9_15Sig(nullptr),
0427       t_e11x11_15Sig(nullptr),
0428       t_e15x15_15Sig(nullptr),
0429       t_e7x7_20Sig(nullptr),
0430       t_e9x9_20Sig(nullptr),
0431       t_e11x11_20Sig(nullptr),
0432       t_e15x15_20Sig(nullptr),
0433       t_e7x7_25Sig(nullptr),
0434       t_e9x9_25Sig(nullptr),
0435       t_e11x11_25Sig(nullptr),
0436       t_e15x15_25Sig(nullptr),
0437       t_e7x7_30Sig(nullptr),
0438       t_e9x9_30Sig(nullptr),
0439       t_e11x11_30Sig(nullptr),
0440       t_e15x15_30Sig(nullptr),
0441       t_esimPdgId(nullptr),
0442       t_simTrackP(nullptr),
0443       t_trkEcalEne(nullptr),
0444       t_esim7x7(nullptr),
0445       t_esim9x9(nullptr),
0446       t_esim11x11(nullptr),
0447       t_esim15x15(nullptr),
0448       t_esim7x7Matched(nullptr),
0449       t_esim9x9Matched(nullptr),
0450       t_esim11x11Matched(nullptr),
0451       t_esim15x15Matched(nullptr),
0452       t_esim7x7Rest(nullptr),
0453       t_esim9x9Rest(nullptr),
0454       t_esim11x11Rest(nullptr),
0455       t_esim15x15Rest(nullptr),
0456       t_esim7x7Photon(nullptr),
0457       t_esim9x9Photon(nullptr),
0458       t_esim11x11Photon(nullptr),
0459       t_esim15x15Photon(nullptr),
0460       t_esim7x7NeutHad(nullptr),
0461       t_esim9x9NeutHad(nullptr),
0462       t_esim11x11NeutHad(nullptr),
0463       t_esim15x15NeutHad(nullptr),
0464       t_esim7x7CharHad(nullptr),
0465       t_esim9x9CharHad(nullptr),
0466       t_esim11x11CharHad(nullptr),
0467       t_esim15x15CharHad(nullptr),
0468       t_maxNearHcalP3x3(nullptr),
0469       t_maxNearHcalP5x5(nullptr),
0470       t_maxNearHcalP7x7(nullptr),
0471       t_h3x3(nullptr),
0472       t_h5x5(nullptr),
0473       t_h7x7(nullptr),
0474       t_h3x3Sig(nullptr),
0475       t_h5x5Sig(nullptr),
0476       t_h7x7Sig(nullptr),
0477       t_infoHcal(nullptr),
0478       t_trkHcalEne(nullptr),
0479       t_hsim3x3(nullptr),
0480       t_hsim5x5(nullptr),
0481       t_hsim7x7(nullptr),
0482       t_hsim3x3Matched(nullptr),
0483       t_hsim5x5Matched(nullptr),
0484       t_hsim7x7Matched(nullptr),
0485       t_hsim3x3Rest(nullptr),
0486       t_hsim5x5Rest(nullptr),
0487       t_hsim7x7Rest(nullptr),
0488       t_hsim3x3Photon(nullptr),
0489       t_hsim5x5Photon(nullptr),
0490       t_hsim7x7Photon(nullptr),
0491       t_hsim3x3NeutHad(nullptr),
0492       t_hsim5x5NeutHad(nullptr),
0493       t_hsim7x7NeutHad(nullptr),
0494       t_hsim3x3CharHad(nullptr),
0495       t_hsim5x5CharHad(nullptr),
0496       t_hsim7x7CharHad(nullptr) {
0497   if (L1TriggerAlgoInfo_) {
0498     m_l1GtUtils = std::make_unique<L1GtUtils>(
0499         iConfig, consumesCollector(), useL1GtTriggerMenuLite, *this, L1GtUtils::UseEventSetupIn::Event);
0500   }
0501 
0502   usesResource(TFileService::kSharedResource);
0503 
0504   //now do what ever initialization is needed
0505 
0506   edm::InputTag L1extraTauJetSource_ = iConfig.getParameter<edm::InputTag>("l1extraTauJetSource");
0507   edm::InputTag L1extraCenJetSource_ = iConfig.getParameter<edm::InputTag>("l1extraCenJetSource");
0508   edm::InputTag L1extraFwdJetSource_ = iConfig.getParameter<edm::InputTag>("l1extraFwdJetSource");
0509   edm::InputTag L1extraMuonSource_ = iConfig.getParameter<edm::InputTag>("l1extraMuonSource");
0510   edm::InputTag L1extraIsoEmSource_ = iConfig.getParameter<edm::InputTag>("l1extraIsoEmSource");
0511   edm::InputTag L1extraNonIsoEmSource_ = iConfig.getParameter<edm::InputTag>("l1extraNonIsoEmSource");
0512   edm::InputTag L1GTReadoutRcdSource_ = iConfig.getParameter<edm::InputTag>("l1GTReadoutRcdSource");
0513   edm::InputTag L1GTObjectMapRcdSource_ = iConfig.getParameter<edm::InputTag>("l1GTObjectMapRcdSource");
0514   edm::InputTag JetSrc_ = iConfig.getParameter<edm::InputTag>("jetSource");
0515   edm::InputTag JetExtender_ = iConfig.getParameter<edm::InputTag>("jetExtender");
0516   edm::InputTag HBHERecHitSource_ = iConfig.getParameter<edm::InputTag>("hbheRecHitSource");
0517 
0518   // define tokens for access
0519   tok_L1extTauJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraTauJetSource_);
0520   tok_L1extCenJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraCenJetSource_);
0521   tok_L1extFwdJet_ = consumes<l1extra::L1JetParticleCollection>(L1extraFwdJetSource_);
0522   tok_L1extMu_ = consumes<l1extra::L1MuonParticleCollection>(L1extraMuonSource_);
0523   tok_L1extIsoEm_ = consumes<l1extra::L1EmParticleCollection>(L1extraIsoEmSource_);
0524   tok_L1extNoIsoEm_ = consumes<l1extra::L1EmParticleCollection>(L1extraNonIsoEmSource_);
0525   tok_jets_ = consumes<reco::CaloJetCollection>(JetSrc_);
0526   tok_hbhe_ = consumes<HBHERecHitCollection>(HBHERecHitSource_);
0527 
0528   tok_genTrack_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
0529   tok_recVtx_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0530   tok_bs_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
0531   tok_EB_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0532   tok_EE_ = consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
0533 
0534   tok_simTk_ = consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"));
0535   tok_simVtx_ = consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"));
0536   tok_caloEB_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEB"));
0537   tok_caloEE_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "EcalHitsEE"));
0538   tok_caloHH_ = consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"));
0539 
0540   if (myverbose_ >= 0) {
0541     edm::LogVerbatim("IsoTrack") << "Parameters read from config file \n"
0542                                  << " doMC        " << doMC_ << "\t myverbose " << myverbose_ << "\t minTrackP "
0543                                  << minTrackP_ << "\t maxTrackEta " << maxTrackEta_ << "\t tMinE " << tMinE_
0544                                  << "\t tMaxE " << tMaxE_ << "\t tMinH " << tMinH_ << "\t tMaxH " << tMaxH_
0545                                  << "\n debugL1Info " << debugL1Info_ << "\t L1TriggerAlgoInfo " << L1TriggerAlgoInfo_
0546                                  << "\n";
0547   }
0548 
0549   tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord>();
0550   tok_caloTopology_ = esConsumes<CaloTopology, CaloTopologyRecord>();
0551   tok_topo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0552   tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0553   tok_ecalChStatus_ = esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0554   tok_sevlv_ = esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>();
0555   tok_htmap_ = esConsumes<EcalTrigTowerConstituentsMap, IdealGeometryRecord>();
0556 }
0557 
0558 IsolatedTracksNxN::~IsolatedTracksNxN() {
0559   delete t_PVx;
0560   delete t_PVy;
0561   delete t_PVz;
0562   delete t_PVisValid;
0563   delete t_PVndof;
0564   delete t_PVNTracks;
0565   delete t_PVNTracksWt;
0566   delete t_PVTracksSumPt;
0567   delete t_PVTracksSumPtWt;
0568   delete t_PVNTracksHP;
0569   delete t_PVNTracksHPWt;
0570   delete t_PVTracksSumPtHP;
0571   delete t_PVTracksSumPtHPWt;
0572   delete t_L1AlgoNames;
0573   delete t_L1PreScale;
0574   delete t_L1CenJetPt;
0575   delete t_L1CenJetEta;
0576   delete t_L1CenJetPhi;
0577   delete t_L1FwdJetPt;
0578   delete t_L1FwdJetEta;
0579   delete t_L1FwdJetPhi;
0580   delete t_L1TauJetPt;
0581   delete t_L1TauJetEta;
0582   delete t_L1TauJetPhi;
0583   delete t_L1MuonPt;
0584   delete t_L1MuonEta;
0585   delete t_L1MuonPhi;
0586   delete t_L1IsoEMPt;
0587   delete t_L1IsoEMEta;
0588   delete t_L1IsoEMPhi;
0589   delete t_L1NonIsoEMPt;
0590   delete t_L1NonIsoEMEta;
0591   delete t_L1NonIsoEMPhi;
0592   delete t_L1METPt;
0593   delete t_L1METEta;
0594   delete t_L1METPhi;
0595   delete t_jetPt;
0596   delete t_jetEta;
0597   delete t_jetPhi;
0598   delete t_nTrksJetCalo;
0599   delete t_nTrksJetVtx;
0600   delete t_trackPAll;
0601   delete t_trackEtaAll;
0602   delete t_trackPhiAll;
0603   delete t_trackPdgIdAll;
0604   delete t_trackPtAll;
0605   delete t_trackDxyAll;
0606   delete t_trackDzAll;
0607   delete t_trackDxyPVAll;
0608   delete t_trackDzPVAll;
0609   delete t_trackChiSqAll;
0610   delete t_trackP;
0611   delete t_trackPt;
0612   delete t_trackEta;
0613   delete t_trackPhi;
0614   delete t_trackEcalEta;
0615   delete t_trackEcalPhi;
0616   delete t_trackHcalEta;
0617   delete t_trackHcalPhi;
0618   delete t_trackNOuterHits;
0619   delete t_NLayersCrossed;
0620   delete t_trackDxy;
0621   delete t_trackDxyBS;
0622   delete t_trackDz;
0623   delete t_trackDzBS;
0624   delete t_trackDxyPV;
0625   delete t_trackDzPV;
0626   delete t_trackPVIdx;
0627   delete t_trackChiSq;
0628   delete t_trackHitsTOB;
0629   delete t_trackHitsTEC;
0630   delete t_trackHitInMissTOB;
0631   delete t_trackHitInMissTEC;
0632   delete t_trackHitInMissTIB;
0633   delete t_trackHitInMissTID;
0634   delete t_trackHitInMissTIBTID;
0635   delete t_trackHitOutMissTOB;
0636   delete t_trackHitOutMissTEC;
0637   delete t_trackHitOutMissTIB;
0638   delete t_trackHitOutMissTID;
0639   delete t_trackHitOutMissTOBTEC;
0640   delete t_trackHitInMeasTOB;
0641   delete t_trackHitInMeasTEC;
0642   delete t_trackHitInMeasTIB;
0643   delete t_trackHitInMeasTID;
0644   delete t_trackHitOutMeasTOB;
0645   delete t_trackHitOutMeasTEC;
0646   delete t_trackHitOutMeasTIB;
0647   delete t_trackHitOutMeasTID;
0648   delete t_trackOutPosOutHitDr;
0649   delete t_trackL;
0650   delete t_maxNearP31x31;
0651   delete t_maxNearP21x21;
0652   delete t_ecalSpike11x11;
0653   delete t_e7x7;
0654   delete t_e9x9;
0655   delete t_e11x11;
0656   delete t_e15x15;
0657   delete t_e7x7_10Sig;
0658   delete t_e9x9_10Sig;
0659   delete t_e11x11_10Sig;
0660   delete t_e15x15_10Sig;
0661   delete t_e7x7_15Sig;
0662   delete t_e9x9_15Sig;
0663   delete t_e11x11_15Sig;
0664   delete t_e15x15_15Sig;
0665   delete t_e7x7_20Sig;
0666   delete t_e9x9_20Sig;
0667   delete t_e11x11_20Sig;
0668   delete t_e15x15_20Sig;
0669   delete t_e7x7_25Sig;
0670   delete t_e9x9_25Sig;
0671   delete t_e11x11_25Sig;
0672   delete t_e15x15_25Sig;
0673   delete t_e7x7_30Sig;
0674   delete t_e9x9_30Sig;
0675   delete t_e11x11_30Sig;
0676   delete t_e15x15_30Sig;
0677   delete t_esim7x7;
0678   delete t_esim9x9;
0679   delete t_esim11x11;
0680   delete t_esim15x15;
0681   delete t_esim7x7Matched;
0682   delete t_esim9x9Matched;
0683   delete t_esim11x11Matched;
0684   delete t_esim15x15Matched;
0685   delete t_esim7x7Rest;
0686   delete t_esim9x9Rest;
0687   delete t_esim11x11Rest;
0688   delete t_esim15x15Rest;
0689   delete t_esim7x7Photon;
0690   delete t_esim9x9Photon;
0691   delete t_esim11x11Photon;
0692   delete t_esim15x15Photon;
0693   delete t_esim7x7NeutHad;
0694   delete t_esim9x9NeutHad;
0695   delete t_esim11x11NeutHad;
0696   delete t_esim15x15NeutHad;
0697   delete t_esim7x7CharHad;
0698   delete t_esim9x9CharHad;
0699   delete t_esim11x11CharHad;
0700   delete t_esim15x15CharHad;
0701   delete t_trkEcalEne;
0702   delete t_simTrackP;
0703   delete t_esimPdgId;
0704   delete t_maxNearHcalP3x3;
0705   delete t_maxNearHcalP5x5;
0706   delete t_maxNearHcalP7x7;
0707   delete t_h3x3;
0708   delete t_h5x5;
0709   delete t_h7x7;
0710   delete t_h3x3Sig;
0711   delete t_h5x5Sig;
0712   delete t_h7x7Sig;
0713   delete t_infoHcal;
0714   delete t_trkHcalEne;
0715   delete t_hsim3x3;
0716   delete t_hsim5x5;
0717   delete t_hsim7x7;
0718   delete t_hsim3x3Matched;
0719   delete t_hsim5x5Matched;
0720   delete t_hsim7x7Matched;
0721   delete t_hsim3x3Rest;
0722   delete t_hsim5x5Rest;
0723   delete t_hsim7x7Rest;
0724   delete t_hsim3x3Photon;
0725   delete t_hsim5x5Photon;
0726   delete t_hsim7x7Photon;
0727   delete t_hsim3x3NeutHad;
0728   delete t_hsim5x5NeutHad;
0729   delete t_hsim7x7NeutHad;
0730   delete t_hsim3x3CharHad;
0731   delete t_hsim5x5CharHad;
0732   delete t_hsim7x7CharHad;
0733 }
0734 
0735 void IsolatedTracksNxN::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0736   edm::ParameterSetDescription desc;
0737   desc.addUntracked<bool>("doMC", false);
0738   desc.addUntracked<bool>("writeAllTracks", false);
0739   desc.addUntracked<int>("verbosity", 1);
0740   desc.addUntracked<double>("pvTracksPtMin", 0.200);
0741   desc.addUntracked<int>("debugTracks", 0);
0742   desc.addUntracked<bool>("printTrkHitPattern", true);
0743   desc.addUntracked<double>("minTrackP", 1.0);
0744   desc.addUntracked<double>("maxTrackEta", 2.6);
0745   desc.addUntracked<bool>("debugL1Info", false);
0746   desc.addUntracked<bool>("l1TriggerAlgoInfo", false);
0747   desc.add<edm::InputTag>("l1extraTauJetSource", edm::InputTag("l1extraParticles", "Tau"));
0748   desc.add<edm::InputTag>("l1extraCenJetSource", edm::InputTag("l1extraParticles", "Central"));
0749   desc.add<edm::InputTag>("l1extraFwdJetSource", edm::InputTag("l1extraParticles", "Forward"));
0750   desc.add<edm::InputTag>("l1extraMuonSource", edm::InputTag("l1extraParticles"));
0751   desc.add<edm::InputTag>("l1extraIsoEmSource", edm::InputTag("l1extraParticles", "Isolated"));
0752   desc.add<edm::InputTag>("l1extraNonIsoEmSource", edm::InputTag("l1extraParticles", "NonIsolated"));
0753   desc.add<edm::InputTag>("l1GTReadoutRcdSource", edm::InputTag("gtDigis"));
0754   desc.add<edm::InputTag>("l1GTObjectMapRcdSource", edm::InputTag("hltL1GtObjectMap"));
0755   desc.add<edm::InputTag>("jetSource", edm::InputTag("iterativeCone5CaloJets"));
0756   desc.add<edm::InputTag>("jetExtender", edm::InputTag("iterativeCone5JetExtender"));
0757   desc.add<edm::InputTag>("hbheRecHitSource", edm::InputTag("hbhereco"));
0758   desc.addUntracked<double>("maxNearTrackPT", 1.0);
0759   desc.addUntracked<double>("timeMinCutECAL", -500.0);
0760   desc.addUntracked<double>("timeMaxCutECAL", 500.0);
0761   desc.addUntracked<double>("timeMinCutHCAL", -500.0);
0762   desc.addUntracked<double>("timeMaxCutHCAL", 500.0);
0763   descriptions.add("isolatedTracksNxN", desc);
0764 }
0765 
0766 void IsolatedTracksNxN::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0767   bool haveIsoTrack = false;
0768 
0769   const MagneticField *bField = &iSetup.getData(tok_magField_);
0770 
0771   clearTreeVectors();
0772 
0773   t_RunNo = iEvent.id().run();
0774   t_EvtNo = iEvent.id().event();
0775   t_Lumi = iEvent.luminosityBlock();
0776   t_Bunch = iEvent.bunchCrossing();
0777 
0778   ++nEventProc_;
0779 
0780   edm::Handle<reco::TrackCollection> trkCollection;
0781   iEvent.getByToken(tok_genTrack_, trkCollection);
0782   if (debugTrks_ > 1) {
0783     edm::LogVerbatim("IsoTrack") << "Track Collection: ";
0784     edm::LogVerbatim("IsoTrack") << "Number of Tracks " << trkCollection->size();
0785   }
0786   std::string theTrackQuality = "highPurity";
0787   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0788 
0789   //===================== save L1 Trigger information =======================
0790   if (L1TriggerAlgoInfo_) {
0791     m_l1GtUtils->getL1GtRunCache(iEvent, iSetup, useL1EventSetup, useL1GtTriggerMenuLite);
0792 
0793     int iErrorCode = -1;
0794     int l1ConfCode = -1;
0795     const bool l1Conf = m_l1GtUtils->availableL1Configuration(iErrorCode, l1ConfCode);
0796     if (!l1Conf) {
0797       edm::LogVerbatim("IsoTrack")
0798           << "\nL1 configuration code:" << l1ConfCode << "\nNo valid L1 trigger configuration available."
0799           << "\nSee text above for error code interpretation"
0800           << "\nNo return here, in order to test each method, protected against configuration error.";
0801     }
0802 
0803     const L1GtTriggerMenu *m_l1GtMenu = m_l1GtUtils->ptrL1TriggerMenuEventSetup(iErrorCode);
0804     const AlgorithmMap &algorithmMap = m_l1GtMenu->gtAlgorithmMap();
0805     const std::string &menuName = m_l1GtMenu->gtTriggerMenuName();
0806 
0807     if (!initL1_) {
0808       initL1_ = true;
0809       edm::LogVerbatim("IsoTrack") << "menuName " << menuName;
0810       for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end(); itAlgo++) {
0811         std::string algName = itAlgo->first;
0812         int algBitNumber = (itAlgo->second).algoBitNumber();
0813         l1AlgoMap_.insert(std::pair<std::pair<unsigned int, std::string>, int>(
0814             std::pair<unsigned int, std::string>(algBitNumber, algName), 0));
0815       }
0816       std::map<std::pair<unsigned int, std::string>, int>::iterator itr;
0817       for (itr = l1AlgoMap_.begin(); itr != l1AlgoMap_.end(); itr++) {
0818         edm::LogVerbatim("IsoTrack") << " ********** " << (itr->first).first << " " << (itr->first).second << " "
0819                                      << itr->second;
0820       }
0821     }
0822 
0823     std::vector<int> algbits;
0824     for (CItAlgo itAlgo = algorithmMap.begin(); itAlgo != algorithmMap.end(); itAlgo++) {
0825       std::string algName = itAlgo->first;
0826       int algBitNumber = (itAlgo->second).algoBitNumber();
0827       bool decision = m_l1GtUtils->decision(iEvent, itAlgo->first, iErrorCode);
0828       int preScale = m_l1GtUtils->prescaleFactor(iEvent, itAlgo->first, iErrorCode);
0829 
0830       // save the algo names which fired
0831       if (decision) {
0832         l1AlgoMap_[std::pair<unsigned int, std::string>(algBitNumber, algName)] += 1;
0833         h_L1AlgoNames->Fill(algBitNumber);
0834         t_L1AlgoNames->push_back(itAlgo->first);
0835         t_L1PreScale->push_back(preScale);
0836         t_L1Decision[algBitNumber] = 1;
0837         algbits.push_back(algBitNumber);
0838       }
0839     }
0840 
0841     if (debugL1Info_) {
0842       for (unsigned int ii = 0; ii < t_L1AlgoNames->size(); ii++) {
0843         edm::LogVerbatim("IsoTrack") << ii << " " << (*t_L1AlgoNames)[ii] << " " << (*t_L1PreScale)[ii] << " "
0844                                      << algbits[ii];
0845       }
0846       for (int i = 0; i < 128; ++i) {
0847         edm::LogVerbatim("IsoTrack") << "L1Decision: " << i << ":" << t_L1Decision[i];
0848       }
0849     }
0850 
0851     // L1Taus
0852     edm::Handle<l1extra::L1JetParticleCollection> l1TauHandle;
0853     iEvent.getByToken(tok_L1extTauJet_, l1TauHandle);
0854     l1extra::L1JetParticleCollection::const_iterator itr;
0855     int iL1Obj = 0;
0856     for (itr = l1TauHandle->begin(), iL1Obj = 0; itr != l1TauHandle->end(); ++itr, iL1Obj++) {
0857       if (iL1Obj < 1) {
0858         t_L1TauJetPt->push_back(itr->pt());
0859         t_L1TauJetEta->push_back(itr->eta());
0860         t_L1TauJetPhi->push_back(itr->phi());
0861       }
0862       if (debugL1Info_) {
0863         edm::LogVerbatim("IsoTrack") << "tauJ p/pt  " << itr->momentum() << " " << itr->pt() << "  eta/phi "
0864                                      << itr->eta() << " " << itr->phi();
0865       }
0866     }
0867 
0868     // L1 Central Jets
0869     edm::Handle<l1extra::L1JetParticleCollection> l1CenJetHandle;
0870     iEvent.getByToken(tok_L1extCenJet_, l1CenJetHandle);
0871     for (itr = l1CenJetHandle->begin(), iL1Obj = 0; itr != l1CenJetHandle->end(); ++itr, iL1Obj++) {
0872       if (iL1Obj < 1) {
0873         t_L1CenJetPt->push_back(itr->pt());
0874         t_L1CenJetEta->push_back(itr->eta());
0875         t_L1CenJetPhi->push_back(itr->phi());
0876       }
0877       if (debugL1Info_) {
0878         edm::LogVerbatim("IsoTrack") << "cenJ p/pt     " << itr->momentum() << " " << itr->pt() << "  eta/phi "
0879                                      << itr->eta() << " " << itr->phi();
0880       }
0881     }
0882 
0883     // L1 Forward Jets
0884     edm::Handle<l1extra::L1JetParticleCollection> l1FwdJetHandle;
0885     iEvent.getByToken(tok_L1extFwdJet_, l1FwdJetHandle);
0886     for (itr = l1FwdJetHandle->begin(), iL1Obj = 0; itr != l1FwdJetHandle->end(); ++itr, iL1Obj++) {
0887       if (iL1Obj < 1) {
0888         t_L1FwdJetPt->push_back(itr->pt());
0889         t_L1FwdJetEta->push_back(itr->eta());
0890         t_L1FwdJetPhi->push_back(itr->phi());
0891       }
0892       if (debugL1Info_) {
0893         edm::LogVerbatim("IsoTrack") << "fwdJ p/pt     " << itr->momentum() << " " << itr->pt() << "  eta/phi "
0894                                      << itr->eta() << " " << itr->phi();
0895       }
0896     }
0897 
0898     // L1 Isolated EM onjects
0899     l1extra::L1EmParticleCollection::const_iterator itrEm;
0900     edm::Handle<l1extra::L1EmParticleCollection> l1IsoEmHandle;
0901     iEvent.getByToken(tok_L1extIsoEm_, l1IsoEmHandle);
0902     for (itrEm = l1IsoEmHandle->begin(), iL1Obj = 0; itrEm != l1IsoEmHandle->end(); ++itrEm, iL1Obj++) {
0903       if (iL1Obj < 1) {
0904         t_L1IsoEMPt->push_back(itrEm->pt());
0905         t_L1IsoEMEta->push_back(itrEm->eta());
0906         t_L1IsoEMPhi->push_back(itrEm->phi());
0907       }
0908       if (debugL1Info_) {
0909         edm::LogVerbatim("IsoTrack") << "isoEm p/pt    " << itrEm->momentum() << " " << itrEm->pt() << "  eta/phi "
0910                                      << itrEm->eta() << " " << itrEm->phi();
0911       }
0912     }
0913 
0914     // L1 Non-Isolated EM onjects
0915     edm::Handle<l1extra::L1EmParticleCollection> l1NonIsoEmHandle;
0916     iEvent.getByToken(tok_L1extNoIsoEm_, l1NonIsoEmHandle);
0917     for (itrEm = l1NonIsoEmHandle->begin(), iL1Obj = 0; itrEm != l1NonIsoEmHandle->end(); ++itrEm, iL1Obj++) {
0918       if (iL1Obj < 1) {
0919         t_L1NonIsoEMPt->push_back(itrEm->pt());
0920         t_L1NonIsoEMEta->push_back(itrEm->eta());
0921         t_L1NonIsoEMPhi->push_back(itrEm->phi());
0922       }
0923       if (debugL1Info_) {
0924         edm::LogVerbatim("IsoTrack") << "nonIsoEm p/pt " << itrEm->momentum() << " " << itrEm->pt() << "  eta/phi "
0925                                      << itrEm->eta() << " " << itrEm->phi();
0926       }
0927     }
0928 
0929     // L1 Muons
0930     l1extra::L1MuonParticleCollection::const_iterator itrMu;
0931     edm::Handle<l1extra::L1MuonParticleCollection> l1MuHandle;
0932     iEvent.getByToken(tok_L1extMu_, l1MuHandle);
0933     for (itrMu = l1MuHandle->begin(), iL1Obj = 0; itrMu != l1MuHandle->end(); ++itrMu, iL1Obj++) {
0934       if (iL1Obj < 1) {
0935         t_L1MuonPt->push_back(itrMu->pt());
0936         t_L1MuonEta->push_back(itrMu->eta());
0937         t_L1MuonPhi->push_back(itrMu->phi());
0938       }
0939       if (debugL1Info_) {
0940         edm::LogVerbatim("IsoTrack") << "l1muon p/pt   " << itrMu->momentum() << " " << itrMu->pt() << "  eta/phi "
0941                                      << itrMu->eta() << " " << itrMu->phi();
0942       }
0943     }
0944   }
0945 
0946   //============== store the information about all the Non-Fake vertices ===============
0947 
0948   edm::Handle<reco::VertexCollection> recVtxs;
0949   iEvent.getByToken(tok_recVtx_, recVtxs);
0950 
0951   std::vector<reco::Track> svTracks;
0952   math::XYZPoint leadPV(0, 0, 0);
0953   double sumPtMax = -1.0;
0954   for (unsigned int ind = 0; ind < recVtxs->size(); ind++) {
0955     if (!((*recVtxs)[ind].isFake())) {
0956       double vtxTrkSumPt = 0.0, vtxTrkSumPtWt = 0.0;
0957       int vtxTrkNWt = 0;
0958       double vtxTrkSumPtHP = 0.0, vtxTrkSumPtHPWt = 0.0;
0959       int vtxTrkNHP = 0, vtxTrkNHPWt = 0;
0960 
0961       reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
0962 
0963       for (vtxTrack = (*recVtxs)[ind].tracks_begin(); vtxTrack != (*recVtxs)[ind].tracks_end(); vtxTrack++) {
0964         if ((*vtxTrack)->pt() < pvTracksPtMin_)
0965           continue;
0966 
0967         bool trkQuality = (*vtxTrack)->quality(trackQuality_);
0968 
0969         vtxTrkSumPt += (*vtxTrack)->pt();
0970 
0971         vtxTrkSumPt += (*vtxTrack)->pt();
0972         if (trkQuality) {
0973           vtxTrkSumPtHP += (*vtxTrack)->pt();
0974           vtxTrkNHP++;
0975         }
0976 
0977         double weight = (*recVtxs)[ind].trackWeight(*vtxTrack);
0978         h_PVTracksWt->Fill(weight);
0979         if (weight > 0.5) {
0980           vtxTrkSumPtWt += (*vtxTrack)->pt();
0981           vtxTrkNWt++;
0982           if (trkQuality) {
0983             vtxTrkSumPtHPWt += (*vtxTrack)->pt();
0984             vtxTrkNHPWt++;
0985           }
0986         }
0987       }
0988 
0989       if (vtxTrkSumPt > sumPtMax) {
0990         sumPtMax = vtxTrkSumPt;
0991         leadPV = math::XYZPoint((*recVtxs)[ind].x(), (*recVtxs)[ind].y(), (*recVtxs)[ind].z());
0992       }
0993 
0994       t_PVx->push_back((*recVtxs)[ind].x());
0995       t_PVy->push_back((*recVtxs)[ind].y());
0996       t_PVz->push_back((*recVtxs)[ind].z());
0997       t_PVisValid->push_back((*recVtxs)[ind].isValid());
0998       t_PVNTracks->push_back((*recVtxs)[ind].tracksSize());
0999       t_PVndof->push_back((*recVtxs)[ind].ndof());
1000       t_PVNTracksWt->push_back(vtxTrkNWt);
1001       t_PVTracksSumPt->push_back(vtxTrkSumPt);
1002       t_PVTracksSumPtWt->push_back(vtxTrkSumPtWt);
1003 
1004       t_PVNTracksHP->push_back(vtxTrkNHP);
1005       t_PVNTracksHPWt->push_back(vtxTrkNHPWt);
1006       t_PVTracksSumPtHP->push_back(vtxTrkSumPtHP);
1007       t_PVTracksSumPtHPWt->push_back(vtxTrkSumPtHPWt);
1008 
1009       if (myverbose_ == 4) {
1010         edm::LogVerbatim("IsoTrack") << "PV " << ind << " isValid " << (*recVtxs)[ind].isValid() << " isFake "
1011                                      << (*recVtxs)[ind].isFake() << " hasRefittedTracks() " << ind << " "
1012                                      << (*recVtxs)[ind].hasRefittedTracks() << " refittedTrksSize "
1013                                      << (*recVtxs)[ind].refittedTracks().size() << "  tracksSize() "
1014                                      << (*recVtxs)[ind].tracksSize() << " sumPt " << vtxTrkSumPt;
1015       }
1016     }  // if vtx is not Fake
1017   }    // loop over PVs
1018   //===================================================================================
1019 
1020   // Get the beamspot
1021   edm::Handle<reco::BeamSpot> beamSpotH;
1022   iEvent.getByToken(tok_bs_, beamSpotH);
1023   math::XYZPoint bspot;
1024   bspot = (beamSpotH.isValid()) ? beamSpotH->position() : math::XYZPoint(0, 0, 0);
1025 
1026   //=====================================================================
1027 
1028   edm::Handle<reco::CaloJetCollection> jets;
1029   iEvent.getByToken(tok_jets_, jets);
1030   //  edm::Handle<reco::JetExtendedAssociation::Container> jetExtender;
1031   //  iEvent.getByLabel(JetExtender_,jetExtender);
1032 
1033   for (unsigned int ijet = 0; ijet < (*jets).size(); ijet++) {
1034     t_jetPt->push_back((*jets)[ijet].pt());
1035     t_jetEta->push_back((*jets)[ijet].eta());
1036     t_jetPhi->push_back((*jets)[ijet].phi());
1037     t_nTrksJetVtx->push_back(-1.0);
1038     t_nTrksJetCalo->push_back(-1.0);
1039   }
1040 
1041   //===================================================================================
1042 
1043   // get handles to calogeometry and calotopology
1044   const CaloGeometry *geo = &iSetup.getData(tok_geom_);
1045   const CaloTopology *caloTopology = &iSetup.getData(tok_caloTopology_);
1046   const HcalTopology *theHBHETopology = &iSetup.getData(tok_topo_);
1047 
1048   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle;
1049   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle;
1050   iEvent.getByToken(tok_EB_, barrelRecHitsHandle);
1051   iEvent.getByToken(tok_EE_, endcapRecHitsHandle);
1052 
1053   // Retrieve the good/bad ECAL channels from the DB
1054   const EcalChannelStatus *theEcalChStatus = &iSetup.getData(tok_ecalChStatus_);
1055   const EcalSeverityLevelAlgo *sevlv = &iSetup.getData(tok_sevlv_);
1056 
1057   // Retrieve trigger tower map
1058   const EcalTrigTowerConstituentsMap &ttMap = iSetup.getData(tok_htmap_);
1059 
1060   edm::Handle<HBHERecHitCollection> hbhe;
1061   iEvent.getByToken(tok_hbhe_, hbhe);
1062   if (!hbhe.isValid()) {
1063     ++nbad_;
1064     if (nbad_ < 10)
1065       edm::LogVerbatim("IsoTrack") << "No HBHE rechit collection";
1066     return;
1067   }
1068   const HBHERecHitCollection Hithbhe = *(hbhe.product());
1069 
1070   //get Handles to SimTracks and SimHits
1071   edm::Handle<edm::SimTrackContainer> SimTk;
1072   if (doMC_)
1073     iEvent.getByToken(tok_simTk_, SimTk);
1074 
1075   edm::Handle<edm::SimVertexContainer> SimVtx;
1076   if (doMC_)
1077     iEvent.getByToken(tok_simVtx_, SimVtx);
1078 
1079   //get Handles to PCaloHitContainers of eb/ee/hbhe
1080   edm::Handle<edm::PCaloHitContainer> pcaloeb;
1081   if (doMC_)
1082     iEvent.getByToken(tok_caloEB_, pcaloeb);
1083 
1084   edm::Handle<edm::PCaloHitContainer> pcaloee;
1085   if (doMC_)
1086     iEvent.getByToken(tok_caloEE_, pcaloee);
1087 
1088   edm::Handle<edm::PCaloHitContainer> pcalohh;
1089   if (doMC_)
1090     iEvent.getByToken(tok_caloHH_, pcalohh);
1091 
1092   //associates tracker rechits/simhits to a track
1093   std::unique_ptr<TrackerHitAssociator> associate;
1094   if (doMC_)
1095     associate = std::make_unique<TrackerHitAssociator>(iEvent, trackerHitAssociatorConfig_);
1096 
1097   //===================================================================================
1098 
1099   h_nTracks->Fill(trkCollection->size());
1100 
1101   int nTracks = 0;
1102 
1103   t_nTracks = trkCollection->size();
1104 
1105   // get the list of DetIds closest to the impact point of track on surface calorimeters
1106   std::vector<spr::propagatedTrackID> trkCaloDets;
1107   spr::propagateCALO(trkCollection, geo, bField, theTrackQuality, trkCaloDets, false);
1108   std::vector<spr::propagatedTrackID>::const_iterator trkDetItr;
1109 
1110   if (myverbose_ > 2) {
1111     for (trkDetItr = trkCaloDets.begin(); trkDetItr != trkCaloDets.end(); trkDetItr++) {
1112       edm::LogVerbatim("IsoTrack") << trkDetItr->trkItr->p() << " " << trkDetItr->trkItr->eta() << " "
1113                                    << trkDetItr->okECAL << " " << trkDetItr->okHCAL;
1114       if (trkDetItr->okECAL) {
1115         if (trkDetItr->detIdECAL.subdetId() == EcalBarrel)
1116           edm::LogVerbatim("IsoTrack") << (EBDetId)trkDetItr->detIdECAL;
1117         else
1118           edm::LogVerbatim("IsoTrack") << (EEDetId)trkDetItr->detIdECAL;
1119       }
1120       if (trkDetItr->okHCAL)
1121         edm::LogVerbatim("IsoTrack") << (HcalDetId)trkDetItr->detIdHCAL;
1122     }
1123   }
1124 
1125   int nvtxTracks = 0;
1126   for (trkDetItr = trkCaloDets.begin(), nTracks = 0; trkDetItr != trkCaloDets.end(); trkDetItr++, nTracks++) {
1127     const reco::Track *pTrack = &(*(trkDetItr->trkItr));
1128 
1129     // find vertex index the track is associated with
1130     int pVtxTkId = -1;
1131     for (unsigned int ind = 0; ind < recVtxs->size(); ind++) {
1132       if (!((*recVtxs)[ind].isFake())) {
1133         reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
1134         for (vtxTrack = (*recVtxs)[ind].tracks_begin(); vtxTrack != (*recVtxs)[ind].tracks_end(); vtxTrack++) {
1135           const edm::RefToBase<reco::Track> pvtxTrack = (*vtxTrack);
1136           if (pTrack == pvtxTrack.get()) {
1137             pVtxTkId = ind;
1138             break;
1139             if (myverbose_ > 2) {
1140               if (pTrack->pt() > 1.0) {
1141                 edm::LogVerbatim("IsoTrack") << "Debug the track association with vertex ";
1142                 edm::LogVerbatim("IsoTrack") << pTrack << " " << pvtxTrack.get();
1143                 edm::LogVerbatim("IsoTrack")
1144                     << " trkVtxIndex " << nvtxTracks << " vtx " << ind << " pt " << pTrack->pt() << " eta "
1145                     << pTrack->eta() << " " << pTrack->pt() - pvtxTrack->pt() << " "
1146                     << pTrack->eta() - pvtxTrack->eta();
1147                 nvtxTracks++;
1148               }
1149             }
1150           }
1151         }
1152       }
1153     }
1154 
1155     const reco::HitPattern &hitp = pTrack->hitPattern();
1156     int nLayersCrossed = hitp.trackerLayersWithMeasurement();
1157     int nOuterHits = hitp.stripTOBLayersWithMeasurement() + hitp.stripTECLayersWithMeasurement();
1158 
1159     bool ifGood = pTrack->quality(trackQuality_);
1160     double pt1 = pTrack->pt();
1161     double p1 = pTrack->p();
1162     double eta1 = pTrack->momentum().eta();
1163     double phi1 = pTrack->momentum().phi();
1164     double etaEcal1 = trkDetItr->etaECAL;
1165     double phiEcal1 = trkDetItr->phiECAL;
1166     double etaHcal1 = trkDetItr->etaHCAL;
1167     double phiHcal1 = trkDetItr->phiHCAL;
1168     double dxy1 = pTrack->dxy();
1169     double dz1 = pTrack->dz();
1170     double chisq1 = pTrack->normalizedChi2();
1171     double dxybs1 = beamSpotH.isValid() ? pTrack->dxy(bspot) : pTrack->dxy();
1172     double dzbs1 = beamSpotH.isValid() ? pTrack->dz(bspot) : pTrack->dz();
1173     double dxypv1 = pTrack->dxy();
1174     double dzpv1 = pTrack->dz();
1175     if (pVtxTkId >= 0) {
1176       math::XYZPoint thisTkPV =
1177           math::XYZPoint((*recVtxs)[pVtxTkId].x(), (*recVtxs)[pVtxTkId].y(), (*recVtxs)[pVtxTkId].z());
1178       dxypv1 = pTrack->dxy(thisTkPV);
1179       dzpv1 = pTrack->dz(thisTkPV);
1180     }
1181 
1182     h_recEtaPt_0->Fill(eta1, pt1);
1183     h_recEtaP_0->Fill(eta1, p1);
1184     h_recPt_0->Fill(pt1);
1185     h_recP_0->Fill(p1);
1186     h_recEta_0->Fill(eta1);
1187     h_recPhi_0->Fill(phi1);
1188 
1189     if (ifGood && nLayersCrossed > 7) {
1190       h_recEtaPt_1->Fill(eta1, pt1);
1191       h_recEtaP_1->Fill(eta1, p1);
1192       h_recPt_1->Fill(pt1);
1193       h_recP_1->Fill(p1);
1194       h_recEta_1->Fill(eta1);
1195       h_recPhi_1->Fill(phi1);
1196     }
1197 
1198     if (!ifGood)
1199       continue;
1200 
1201     if (writeAllTracks_ && p1 > 2.0 && nLayersCrossed > 7) {
1202       t_trackPAll->push_back(p1);
1203       t_trackEtaAll->push_back(eta1);
1204       t_trackPhiAll->push_back(phi1);
1205       t_trackPtAll->push_back(pt1);
1206       t_trackDxyAll->push_back(dxy1);
1207       t_trackDzAll->push_back(dz1);
1208       t_trackDxyPVAll->push_back(dxypv1);
1209       t_trackDzPVAll->push_back(dzpv1);
1210       t_trackChiSqAll->push_back(chisq1);
1211     }
1212     if (doMC_) {
1213       edm::SimTrackContainer::const_iterator matchedSimTrkAll =
1214           spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false);
1215       if (writeAllTracks_ && matchedSimTrkAll != SimTk->end())
1216         t_trackPdgIdAll->push_back(matchedSimTrkAll->type());
1217     }
1218 
1219     if (pt1 > minTrackP_ && std::abs(eta1) < maxTrackEta_ && trkDetItr->okECAL) {
1220       double maxNearP31x31 = 999.0, maxNearP25x25 = 999.0, maxNearP21x21 = 999.0, maxNearP15x15 = 999.0;
1221       maxNearP31x31 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 15, 15);
1222       maxNearP25x25 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 12, 12);
1223       maxNearP21x21 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 10, 10);
1224       maxNearP15x15 = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 7, 7);
1225 
1226       int iTrkEtaBin = -1, iTrkMomBin = -1;
1227       for (unsigned int ieta = 0; ieta < NEtaBins; ieta++) {
1228         if (std::abs(eta1) > genPartEtaBins[ieta] && std::abs(eta1) < genPartEtaBins[ieta + 1])
1229           iTrkEtaBin = ieta;
1230       }
1231       for (unsigned int ipt = 0; ipt < NPBins; ipt++) {
1232         if (p1 > genPartPBins[ipt] && p1 < genPartPBins[ipt + 1])
1233           iTrkMomBin = ipt;
1234       }
1235       if (iTrkMomBin >= 0 && iTrkEtaBin >= 0) {
1236         h_maxNearP31x31[iTrkMomBin][iTrkEtaBin]->Fill(maxNearP31x31);
1237         h_maxNearP25x25[iTrkMomBin][iTrkEtaBin]->Fill(maxNearP25x25);
1238         h_maxNearP21x21[iTrkMomBin][iTrkEtaBin]->Fill(maxNearP21x21);
1239         h_maxNearP15x15[iTrkMomBin][iTrkEtaBin]->Fill(maxNearP15x15);
1240       }
1241       if (maxNearP31x31 < 0.0 && nLayersCrossed > 7 && nOuterHits > 4) {
1242         h_recEtaPt_2->Fill(eta1, pt1);
1243         h_recEtaP_2->Fill(eta1, p1);
1244         h_recPt_2->Fill(pt1);
1245         h_recP_2->Fill(p1);
1246         h_recEta_2->Fill(eta1);
1247         h_recPhi_2->Fill(phi1);
1248       }
1249 
1250       // if charge isolated in ECAL, store the further quantities
1251       if (maxNearP31x31 < 1.0) {
1252         haveIsoTrack = true;
1253 
1254         // get the matching simTrack
1255         double simTrackP = -1;
1256         if (doMC_) {
1257           edm::SimTrackContainer::const_iterator matchedSimTrk =
1258               spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, *associate, false);
1259           if (matchedSimTrk != SimTk->end())
1260             simTrackP = matchedSimTrk->momentum().P();
1261         }
1262 
1263         // get ECal Tranverse Profile
1264         std::pair<double, bool> e7x7P, e9x9P, e11x11P, e15x15P;
1265         std::pair<double, bool> e7x7_10SigP, e9x9_10SigP, e11x11_10SigP, e15x15_10SigP;
1266         std::pair<double, bool> e7x7_15SigP, e9x9_15SigP, e11x11_15SigP, e15x15_15SigP;
1267         std::pair<double, bool> e7x7_20SigP, e9x9_20SigP, e11x11_20SigP, e15x15_20SigP;
1268         std::pair<double, bool> e7x7_25SigP, e9x9_25SigP, e11x11_25SigP, e15x15_25SigP;
1269         std::pair<double, bool> e7x7_30SigP, e9x9_30SigP, e11x11_30SigP, e15x15_30SigP;
1270 
1271         spr::caloSimInfo simInfo3x3, simInfo5x5, simInfo7x7, simInfo9x9;
1272         spr::caloSimInfo simInfo11x11, simInfo13x13, simInfo15x15, simInfo21x21, simInfo25x25, simInfo31x31;
1273         double trkEcalEne = 0;
1274 
1275         const DetId isoCell = trkDetItr->detIdECAL;
1276         e7x7P = spr::eECALmatrix(isoCell,
1277                                  barrelRecHitsHandle,
1278                                  endcapRecHitsHandle,
1279                                  *theEcalChStatus,
1280                                  geo,
1281                                  caloTopology,
1282                                  sevlv,
1283                                  3,
1284                                  3,
1285                                  -100.0,
1286                                  -100.0,
1287                                  tMinE_,
1288                                  tMaxE_);
1289         e9x9P = spr::eECALmatrix(isoCell,
1290                                  barrelRecHitsHandle,
1291                                  endcapRecHitsHandle,
1292                                  *theEcalChStatus,
1293                                  geo,
1294                                  caloTopology,
1295                                  sevlv,
1296                                  4,
1297                                  4,
1298                                  -100.0,
1299                                  -100.0,
1300                                  tMinE_,
1301                                  tMaxE_);
1302         e11x11P = spr::eECALmatrix(isoCell,
1303                                    barrelRecHitsHandle,
1304                                    endcapRecHitsHandle,
1305                                    *theEcalChStatus,
1306                                    geo,
1307                                    caloTopology,
1308                                    sevlv,
1309                                    5,
1310                                    5,
1311                                    -100.0,
1312                                    -100.0,
1313                                    tMinE_,
1314                                    tMaxE_);
1315         e15x15P = spr::eECALmatrix(isoCell,
1316                                    barrelRecHitsHandle,
1317                                    endcapRecHitsHandle,
1318                                    *theEcalChStatus,
1319                                    geo,
1320                                    caloTopology,
1321                                    sevlv,
1322                                    7,
1323                                    7,
1324                                    -100.0,
1325                                    -100.0,
1326                                    tMinE_,
1327                                    tMaxE_);
1328 
1329         e7x7_10SigP = spr::eECALmatrix(isoCell,
1330                                        barrelRecHitsHandle,
1331                                        endcapRecHitsHandle,
1332                                        *theEcalChStatus,
1333                                        geo,
1334                                        caloTopology,
1335                                        sevlv,
1336                                        3,
1337                                        3,
1338                                        0.030,
1339                                        0.150,
1340                                        tMinE_,
1341                                        tMaxE_);
1342         e9x9_10SigP = spr::eECALmatrix(isoCell,
1343                                        barrelRecHitsHandle,
1344                                        endcapRecHitsHandle,
1345                                        *theEcalChStatus,
1346                                        geo,
1347                                        caloTopology,
1348                                        sevlv,
1349                                        4,
1350                                        4,
1351                                        0.030,
1352                                        0.150,
1353                                        tMinE_,
1354                                        tMaxE_);
1355         e11x11_10SigP = spr::eECALmatrix(isoCell,
1356                                          barrelRecHitsHandle,
1357                                          endcapRecHitsHandle,
1358                                          *theEcalChStatus,
1359                                          geo,
1360                                          caloTopology,
1361                                          sevlv,
1362                                          5,
1363                                          5,
1364                                          0.030,
1365                                          0.150,
1366                                          tMinE_,
1367                                          tMaxE_);
1368         e15x15_10SigP = spr::eECALmatrix(isoCell,
1369                                          barrelRecHitsHandle,
1370                                          endcapRecHitsHandle,
1371                                          *theEcalChStatus,
1372                                          geo,
1373                                          caloTopology,
1374                                          sevlv,
1375                                          7,
1376                                          7,
1377                                          0.030,
1378                                          0.150,
1379                                          tMinE_,
1380                                          tMaxE_);
1381 
1382         e7x7_15SigP = spr::eECALmatrix(isoCell,
1383                                        barrelRecHitsHandle,
1384                                        endcapRecHitsHandle,
1385                                        *theEcalChStatus,
1386                                        geo,
1387                                        caloTopology,
1388                                        sevlv,
1389                                        ttMap,
1390                                        3,
1391                                        3,
1392                                        0.20,
1393                                        0.45,
1394                                        tMinE_,
1395                                        tMaxE_);
1396         e9x9_15SigP = spr::eECALmatrix(isoCell,
1397                                        barrelRecHitsHandle,
1398                                        endcapRecHitsHandle,
1399                                        *theEcalChStatus,
1400                                        geo,
1401                                        caloTopology,
1402                                        sevlv,
1403                                        ttMap,
1404                                        4,
1405                                        4,
1406                                        0.20,
1407                                        0.45,
1408                                        tMinE_,
1409                                        tMaxE_);
1410         e11x11_15SigP = spr::eECALmatrix(isoCell,
1411                                          barrelRecHitsHandle,
1412                                          endcapRecHitsHandle,
1413                                          *theEcalChStatus,
1414                                          geo,
1415                                          caloTopology,
1416                                          sevlv,
1417                                          ttMap,
1418                                          5,
1419                                          5,
1420                                          0.20,
1421                                          0.45,
1422                                          tMinE_,
1423                                          tMaxE_);
1424         e15x15_15SigP = spr::eECALmatrix(isoCell,
1425                                          barrelRecHitsHandle,
1426                                          endcapRecHitsHandle,
1427                                          *theEcalChStatus,
1428                                          geo,
1429                                          caloTopology,
1430                                          sevlv,
1431                                          ttMap,
1432                                          7,
1433                                          7,
1434                                          0.20,
1435                                          0.45,
1436                                          tMinE_,
1437                                          tMaxE_,
1438                                          false);
1439 
1440         e7x7_20SigP = spr::eECALmatrix(isoCell,
1441                                        barrelRecHitsHandle,
1442                                        endcapRecHitsHandle,
1443                                        *theEcalChStatus,
1444                                        geo,
1445                                        caloTopology,
1446                                        sevlv,
1447                                        3,
1448                                        3,
1449                                        0.060,
1450                                        0.300,
1451                                        tMinE_,
1452                                        tMaxE_);
1453         e9x9_20SigP = spr::eECALmatrix(isoCell,
1454                                        barrelRecHitsHandle,
1455                                        endcapRecHitsHandle,
1456                                        *theEcalChStatus,
1457                                        geo,
1458                                        caloTopology,
1459                                        sevlv,
1460                                        4,
1461                                        4,
1462                                        0.060,
1463                                        0.300,
1464                                        tMinE_,
1465                                        tMaxE_);
1466         e11x11_20SigP = spr::eECALmatrix(isoCell,
1467                                          barrelRecHitsHandle,
1468                                          endcapRecHitsHandle,
1469                                          *theEcalChStatus,
1470                                          geo,
1471                                          caloTopology,
1472                                          sevlv,
1473                                          5,
1474                                          5,
1475                                          0.060,
1476                                          0.300,
1477                                          tMinE_,
1478                                          tMaxE_);
1479         e15x15_20SigP = spr::eECALmatrix(isoCell,
1480                                          barrelRecHitsHandle,
1481                                          endcapRecHitsHandle,
1482                                          *theEcalChStatus,
1483                                          geo,
1484                                          caloTopology,
1485                                          sevlv,
1486                                          7,
1487                                          7,
1488                                          0.060,
1489                                          0.300,
1490                                          tMinE_,
1491                                          tMaxE_);
1492 
1493         e7x7_25SigP = spr::eECALmatrix(isoCell,
1494                                        barrelRecHitsHandle,
1495                                        endcapRecHitsHandle,
1496                                        *theEcalChStatus,
1497                                        geo,
1498                                        caloTopology,
1499                                        sevlv,
1500                                        3,
1501                                        3,
1502                                        0.075,
1503                                        0.375,
1504                                        tMinE_,
1505                                        tMaxE_);
1506         e9x9_25SigP = spr::eECALmatrix(isoCell,
1507                                        barrelRecHitsHandle,
1508                                        endcapRecHitsHandle,
1509                                        *theEcalChStatus,
1510                                        geo,
1511                                        caloTopology,
1512                                        sevlv,
1513                                        4,
1514                                        4,
1515                                        0.075,
1516                                        0.375,
1517                                        tMinE_,
1518                                        tMaxE_);
1519         e11x11_25SigP = spr::eECALmatrix(isoCell,
1520                                          barrelRecHitsHandle,
1521                                          endcapRecHitsHandle,
1522                                          *theEcalChStatus,
1523                                          geo,
1524                                          caloTopology,
1525                                          sevlv,
1526                                          5,
1527                                          5,
1528                                          0.075,
1529                                          0.375,
1530                                          tMinE_,
1531                                          tMaxE_);
1532         e15x15_25SigP = spr::eECALmatrix(isoCell,
1533                                          barrelRecHitsHandle,
1534                                          endcapRecHitsHandle,
1535                                          *theEcalChStatus,
1536                                          geo,
1537                                          caloTopology,
1538                                          sevlv,
1539                                          7,
1540                                          7,
1541                                          0.075,
1542                                          0.375,
1543                                          tMinE_,
1544                                          tMaxE_);
1545 
1546         e7x7_30SigP = spr::eECALmatrix(isoCell,
1547                                        barrelRecHitsHandle,
1548                                        endcapRecHitsHandle,
1549                                        *theEcalChStatus,
1550                                        geo,
1551                                        caloTopology,
1552                                        sevlv,
1553                                        3,
1554                                        3,
1555                                        0.090,
1556                                        0.450,
1557                                        tMinE_,
1558                                        tMaxE_);
1559         e9x9_30SigP = spr::eECALmatrix(isoCell,
1560                                        barrelRecHitsHandle,
1561                                        endcapRecHitsHandle,
1562                                        *theEcalChStatus,
1563                                        geo,
1564                                        caloTopology,
1565                                        sevlv,
1566                                        4,
1567                                        4,
1568                                        0.090,
1569                                        0.450,
1570                                        tMinE_,
1571                                        tMaxE_);
1572         e11x11_30SigP = spr::eECALmatrix(isoCell,
1573                                          barrelRecHitsHandle,
1574                                          endcapRecHitsHandle,
1575                                          *theEcalChStatus,
1576                                          geo,
1577                                          caloTopology,
1578                                          sevlv,
1579                                          5,
1580                                          5,
1581                                          0.090,
1582                                          0.450,
1583                                          tMinE_,
1584                                          tMaxE_);
1585         e15x15_30SigP = spr::eECALmatrix(isoCell,
1586                                          barrelRecHitsHandle,
1587                                          endcapRecHitsHandle,
1588                                          *theEcalChStatus,
1589                                          geo,
1590                                          caloTopology,
1591                                          sevlv,
1592                                          7,
1593                                          7,
1594                                          0.090,
1595                                          0.450,
1596                                          tMinE_,
1597                                          tMaxE_);
1598         if (myverbose_ == 2) {
1599           edm::LogVerbatim("IsoTrack") << "clean  ecal rechit ";
1600           edm::LogVerbatim("IsoTrack") << "e7x7 " << e7x7P.first << " e9x9 " << e9x9P.first << " e11x11 "
1601                                        << e11x11P.first << " e15x15 " << e15x15P.first;
1602           edm::LogVerbatim("IsoTrack") << "e7x7_10Sig " << e7x7_10SigP.first << " e11x11_10Sig " << e11x11_10SigP.first
1603                                        << " e15x15_10Sig " << e15x15_10SigP.first;
1604         }
1605 
1606         if (doMC_) {
1607           // check the energy from SimHits
1608           spr::eECALSimInfo(
1609               iEvent, isoCell, geo, caloTopology, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate, 1, 1, simInfo3x3);
1610           spr::eECALSimInfo(
1611               iEvent, isoCell, geo, caloTopology, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate, 2, 2, simInfo5x5);
1612           spr::eECALSimInfo(
1613               iEvent, isoCell, geo, caloTopology, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate, 3, 3, simInfo7x7);
1614           spr::eECALSimInfo(
1615               iEvent, isoCell, geo, caloTopology, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate, 4, 4, simInfo9x9);
1616           spr::eECALSimInfo(iEvent,
1617                             isoCell,
1618                             geo,
1619                             caloTopology,
1620                             pcaloeb,
1621                             pcaloee,
1622                             SimTk,
1623                             SimVtx,
1624                             pTrack,
1625                             *associate,
1626                             5,
1627                             5,
1628                             simInfo11x11);
1629           spr::eECALSimInfo(iEvent,
1630                             isoCell,
1631                             geo,
1632                             caloTopology,
1633                             pcaloeb,
1634                             pcaloee,
1635                             SimTk,
1636                             SimVtx,
1637                             pTrack,
1638                             *associate,
1639                             6,
1640                             6,
1641                             simInfo13x13);
1642           spr::eECALSimInfo(iEvent,
1643                             isoCell,
1644                             geo,
1645                             caloTopology,
1646                             pcaloeb,
1647                             pcaloee,
1648                             SimTk,
1649                             SimVtx,
1650                             pTrack,
1651                             *associate,
1652                             7,
1653                             7,
1654                             simInfo15x15,
1655                             150.0,
1656                             false);
1657           spr::eECALSimInfo(iEvent,
1658                             isoCell,
1659                             geo,
1660                             caloTopology,
1661                             pcaloeb,
1662                             pcaloee,
1663                             SimTk,
1664                             SimVtx,
1665                             pTrack,
1666                             *associate,
1667                             10,
1668                             10,
1669                             simInfo21x21);
1670           spr::eECALSimInfo(iEvent,
1671                             isoCell,
1672                             geo,
1673                             caloTopology,
1674                             pcaloeb,
1675                             pcaloee,
1676                             SimTk,
1677                             SimVtx,
1678                             pTrack,
1679                             *associate,
1680                             12,
1681                             12,
1682                             simInfo25x25);
1683           spr::eECALSimInfo(iEvent,
1684                             isoCell,
1685                             geo,
1686                             caloTopology,
1687                             pcaloeb,
1688                             pcaloee,
1689                             SimTk,
1690                             SimVtx,
1691                             pTrack,
1692                             *associate,
1693                             15,
1694                             15,
1695                             simInfo31x31);
1696 
1697           trkEcalEne = spr::eCaloSimInfo(iEvent, geo, pcaloeb, pcaloee, SimTk, SimVtx, pTrack, *associate);
1698           if (myverbose_ == 1) {
1699             edm::LogVerbatim("IsoTrack") << "Track momentum " << pt1;
1700             edm::LogVerbatim("IsoTrack") << "ecal siminfo ";
1701             edm::LogVerbatim("IsoTrack") << "simInfo3x3: eTotal " << simInfo3x3.eTotal << " eMatched "
1702                                          << simInfo3x3.eMatched << " eRest " << simInfo3x3.eRest << " eGamma "
1703                                          << simInfo3x3.eGamma << " eNeutralHad " << simInfo3x3.eNeutralHad
1704                                          << " eChargedHad " << simInfo3x3.eChargedHad;
1705             edm::LogVerbatim("IsoTrack") << "simInfo5x5: eTotal " << simInfo5x5.eTotal << " eMatched "
1706                                          << simInfo5x5.eMatched << " eRest " << simInfo5x5.eRest << " eGamma "
1707                                          << simInfo5x5.eGamma << " eNeutralHad " << simInfo5x5.eNeutralHad
1708                                          << " eChargedHad " << simInfo5x5.eChargedHad;
1709             edm::LogVerbatim("IsoTrack") << "simInfo7x7: eTotal " << simInfo7x7.eTotal << " eMatched "
1710                                          << simInfo7x7.eMatched << " eRest " << simInfo7x7.eRest << " eGamma "
1711                                          << simInfo7x7.eGamma << " eNeutralHad " << simInfo7x7.eNeutralHad
1712                                          << " eChargedHad " << simInfo7x7.eChargedHad;
1713             edm::LogVerbatim("IsoTrack") << "simInfo9x9: eTotal " << simInfo9x9.eTotal << " eMatched "
1714                                          << simInfo9x9.eMatched << " eRest " << simInfo9x9.eRest << " eGamma "
1715                                          << simInfo9x9.eGamma << " eNeutralHad " << simInfo9x9.eNeutralHad
1716                                          << " eChargedHad " << simInfo9x9.eChargedHad;
1717             edm::LogVerbatim("IsoTrack") << "simInfo11x11: eTotal " << simInfo11x11.eTotal << " eMatched "
1718                                          << simInfo11x11.eMatched << " eRest " << simInfo11x11.eRest << " eGamma "
1719                                          << simInfo11x11.eGamma << " eNeutralHad " << simInfo11x11.eNeutralHad
1720                                          << " eChargedHad " << simInfo11x11.eChargedHad;
1721             edm::LogVerbatim("IsoTrack") << "simInfo15x15: eTotal " << simInfo15x15.eTotal << " eMatched "
1722                                          << simInfo15x15.eMatched << " eRest " << simInfo15x15.eRest << " eGamma "
1723                                          << simInfo15x15.eGamma << " eNeutralHad " << simInfo15x15.eNeutralHad
1724                                          << " eChargedHad " << simInfo15x15.eChargedHad;
1725             edm::LogVerbatim("IsoTrack") << "simInfo31x31: eTotal " << simInfo31x31.eTotal << " eMatched "
1726                                          << simInfo31x31.eMatched << " eRest " << simInfo31x31.eRest << " eGamma "
1727                                          << simInfo31x31.eGamma << " eNeutralHad " << simInfo31x31.eNeutralHad
1728                                          << " eChargedHad " << simInfo31x31.eChargedHad;
1729             edm::LogVerbatim("IsoTrack") << "trkEcalEne" << trkEcalEne;
1730           }
1731         }
1732 
1733         // =======  Get HCAL information
1734         double hcalScale = 1.0;
1735         if (std::abs(pTrack->eta()) < 1.4) {
1736           hcalScale = 120.0;
1737         } else {
1738           hcalScale = 135.0;
1739         }
1740 
1741         double maxNearHcalP3x3 = -1, maxNearHcalP5x5 = -1, maxNearHcalP7x7 = -1;
1742         maxNearHcalP3x3 = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 1, 1);
1743         maxNearHcalP5x5 = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 2, 2);
1744         maxNearHcalP7x7 = spr::chargeIsolationHcal(nTracks, trkCaloDets, theHBHETopology, 3, 3);
1745 
1746         double h3x3 = 0, h5x5 = 0, h7x7 = 0;
1747         double h3x3Sig = 0, h5x5Sig = 0, h7x7Sig = 0;
1748         double trkHcalEne = 0;
1749         spr::caloSimInfo hsimInfo3x3, hsimInfo5x5, hsimInfo7x7;
1750 
1751         if (trkDetItr->okHCAL) {
1752           const DetId ClosestCell(trkDetItr->detIdHCAL);
1753           // bool includeHO=false, bool algoNew=true, bool debug=false
1754           h3x3 = spr::eHCALmatrix(
1755               theHBHETopology, ClosestCell, hbhe, 1, 1, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_);
1756           h5x5 = spr::eHCALmatrix(
1757               theHBHETopology, ClosestCell, hbhe, 2, 2, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_);
1758           h7x7 = spr::eHCALmatrix(
1759               theHBHETopology, ClosestCell, hbhe, 3, 3, false, true, -100.0, -100.0, -100.0, -100.0, tMinH_, tMaxH_);
1760           h3x3Sig = spr::eHCALmatrix(
1761               theHBHETopology, ClosestCell, hbhe, 1, 1, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_, tMaxH_);
1762           h5x5Sig = spr::eHCALmatrix(
1763               theHBHETopology, ClosestCell, hbhe, 2, 2, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_, tMaxH_);
1764           h7x7Sig = spr::eHCALmatrix(
1765               theHBHETopology, ClosestCell, hbhe, 3, 3, false, true, 0.7, 0.8, -100.0, -100.0, tMinH_, tMaxH_);
1766 
1767           if (myverbose_ == 2) {
1768             edm::LogVerbatim("IsoTrack") << "HCAL 3x3 " << h3x3 << " " << h3x3Sig << " 5x5 " << h5x5 << " " << h5x5Sig
1769                                          << " 7x7 " << h7x7 << " " << h7x7Sig;
1770           }
1771 
1772           if (doMC_) {
1773             spr::eHCALSimInfo(
1774                 iEvent, theHBHETopology, ClosestCell, geo, pcalohh, SimTk, SimVtx, pTrack, *associate, 1, 1, hsimInfo3x3);
1775             spr::eHCALSimInfo(
1776                 iEvent, theHBHETopology, ClosestCell, geo, pcalohh, SimTk, SimVtx, pTrack, *associate, 2, 2, hsimInfo5x5);
1777             spr::eHCALSimInfo(iEvent,
1778                               theHBHETopology,
1779                               ClosestCell,
1780                               geo,
1781                               pcalohh,
1782                               SimTk,
1783                               SimVtx,
1784                               pTrack,
1785                               *associate,
1786                               3,
1787                               3,
1788                               hsimInfo7x7,
1789                               150.0,
1790                               false,
1791                               false);
1792             trkHcalEne = spr::eCaloSimInfo(iEvent, geo, pcalohh, SimTk, SimVtx, pTrack, *associate);
1793             if (myverbose_ == 1) {
1794               edm::LogVerbatim("IsoTrack") << "Hcal siminfo ";
1795               edm::LogVerbatim("IsoTrack")
1796                   << "hsimInfo3x3: eTotal " << hsimInfo3x3.eTotal << " eMatched " << hsimInfo3x3.eMatched << " eRest "
1797                   << hsimInfo3x3.eRest << " eGamma " << hsimInfo3x3.eGamma << " eNeutralHad " << hsimInfo3x3.eNeutralHad
1798                   << " eChargedHad " << hsimInfo3x3.eChargedHad;
1799               edm::LogVerbatim("IsoTrack")
1800                   << "hsimInfo5x5: eTotal " << hsimInfo5x5.eTotal << " eMatched " << hsimInfo5x5.eMatched << " eRest "
1801                   << hsimInfo5x5.eRest << " eGamma " << hsimInfo5x5.eGamma << " eNeutralHad " << hsimInfo5x5.eNeutralHad
1802                   << " eChargedHad " << hsimInfo5x5.eChargedHad;
1803               edm::LogVerbatim("IsoTrack")
1804                   << "hsimInfo7x7: eTotal " << hsimInfo7x7.eTotal << " eMatched " << hsimInfo7x7.eMatched << " eRest "
1805                   << hsimInfo7x7.eRest << " eGamma " << hsimInfo7x7.eGamma << " eNeutralHad " << hsimInfo7x7.eNeutralHad
1806                   << " eChargedHad " << hsimInfo7x7.eChargedHad;
1807               edm::LogVerbatim("IsoTrack") << "trkHcalEne " << trkHcalEne;
1808             }
1809           }
1810 
1811           // debug the ecal and hcal matrix
1812           if (myverbose_ == 4) {
1813             edm::LogVerbatim("IsoTrack") << "Run " << iEvent.id().run() << "  Event " << iEvent.id().event();
1814             std::vector<std::pair<DetId, double> > v7x7 =
1815                 spr::eHCALmatrixCell(theHBHETopology, ClosestCell, hbhe, 3, 3, false, false);
1816             double sumv = 0.0;
1817 
1818             for (unsigned int iv = 0; iv < v7x7.size(); iv++) {
1819               sumv += v7x7[iv].second;
1820             }
1821             edm::LogVerbatim("IsoTrack") << "h7x7 " << h7x7 << " v7x7 " << sumv << " in " << v7x7.size();
1822             for (unsigned int iv = 0; iv < v7x7.size(); iv++) {
1823               HcalDetId id = v7x7[iv].first;
1824               edm::LogVerbatim("IsoTrack") << " Cell " << iv << " 0x" << std::hex << v7x7[iv].first() << std::dec << " "
1825                                            << id << " Energy " << v7x7[iv].second;
1826             }
1827           }
1828         }
1829         if (doMC_) {
1830           trkHcalEne = spr::eCaloSimInfo(iEvent, geo, pcalohh, SimTk, SimVtx, pTrack, *associate);
1831         }
1832 
1833         // ====================================================================================================
1834         // get diff between track outermost hit position and the propagation point at outermost surface of tracker
1835         std::pair<math::XYZPoint, double> point2_TK0 = spr::propagateTrackerEnd(pTrack, bField, false);
1836         math::XYZPoint diff(pTrack->outerPosition().X() - point2_TK0.first.X(),
1837                             pTrack->outerPosition().Y() - point2_TK0.first.Y(),
1838                             pTrack->outerPosition().Z() - point2_TK0.first.Z());
1839         double trackOutPosOutHitDr = diff.R();
1840         double trackL = point2_TK0.second;
1841         if (myverbose_ == 5) {
1842           edm::LogVerbatim("IsoTrack") << " propagted " << point2_TK0.first << " " << point2_TK0.first.eta() << " "
1843                                        << point2_TK0.first.phi();
1844           edm::LogVerbatim("IsoTrack") << " outerPosition() " << pTrack->outerPosition() << " "
1845                                        << pTrack->outerPosition().eta() << " " << pTrack->outerPosition().phi();
1846           edm::LogVerbatim("IsoTrack") << "diff " << diff << " diffR " << diff.R() << " diffR/L "
1847                                        << diff.R() / point2_TK0.second;
1848         }
1849 
1850         for (unsigned int ind = 0; ind < recVtxs->size(); ind++) {
1851           if (!((*recVtxs)[ind].isFake())) {
1852             reco::Vertex::trackRef_iterator vtxTrack = (*recVtxs)[ind].tracks_begin();
1853             if (reco::deltaR(eta1, phi1, (*vtxTrack)->eta(), (*vtxTrack)->phi()) < 0.01)
1854               t_trackPVIdx->push_back(ind);
1855             else
1856               t_trackPVIdx->push_back(-1);
1857           }
1858         }
1859 
1860         // Fill the tree Branches here
1861         t_trackPVIdx->push_back(pVtxTkId);
1862         t_trackP->push_back(p1);
1863         t_trackPt->push_back(pt1);
1864         t_trackEta->push_back(eta1);
1865         t_trackPhi->push_back(phi1);
1866         t_trackEcalEta->push_back(etaEcal1);
1867         t_trackEcalPhi->push_back(phiEcal1);
1868         t_trackHcalEta->push_back(etaHcal1);
1869         t_trackHcalPhi->push_back(phiHcal1);
1870         t_trackDxy->push_back(dxy1);
1871         t_trackDz->push_back(dz1);
1872         t_trackDxyBS->push_back(dxybs1);
1873         t_trackDzBS->push_back(dzbs1);
1874         t_trackDxyPV->push_back(dxypv1);
1875         t_trackDzPV->push_back(dzpv1);
1876         t_trackChiSq->push_back(chisq1);
1877         t_trackNOuterHits->push_back(nOuterHits);
1878         t_NLayersCrossed->push_back(nLayersCrossed);
1879 
1880         t_trackHitsTOB->push_back(hitp.stripTOBLayersWithMeasurement());
1881         t_trackHitsTEC->push_back(hitp.stripTECLayersWithMeasurement());
1882         t_trackHitInMissTOB->push_back(hitp.stripTOBLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1883         t_trackHitInMissTEC->push_back(hitp.stripTECLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1884         t_trackHitInMissTIB->push_back(hitp.stripTIBLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1885         t_trackHitInMissTID->push_back(hitp.stripTIDLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1886         t_trackHitInMissTIBTID->push_back(hitp.stripTIBLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS) +
1887                                           hitp.stripTIDLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1888 
1889         t_trackHitOutMissTOB->push_back(hitp.stripTOBLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1890         t_trackHitOutMissTEC->push_back(hitp.stripTECLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1891         t_trackHitOutMissTIB->push_back(hitp.stripTIBLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1892         t_trackHitOutMissTID->push_back(hitp.stripTIDLayersWithoutMeasurement(reco::HitPattern::MISSING_INNER_HITS));
1893         t_trackHitOutMissTOBTEC->push_back(hitp.stripTOBLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS) +
1894                                            hitp.stripTECLayersWithoutMeasurement(reco::HitPattern::MISSING_OUTER_HITS));
1895 
1896         t_trackHitInMeasTOB->push_back(hitp.stripTOBLayersWithMeasurement());
1897         t_trackHitInMeasTEC->push_back(hitp.stripTECLayersWithMeasurement());
1898         t_trackHitInMeasTIB->push_back(hitp.stripTIBLayersWithMeasurement());
1899         t_trackHitInMeasTID->push_back(hitp.stripTIDLayersWithMeasurement());
1900         t_trackHitOutMeasTOB->push_back(hitp.stripTOBLayersWithMeasurement());
1901         t_trackHitOutMeasTEC->push_back(hitp.stripTECLayersWithMeasurement());
1902         t_trackHitOutMeasTIB->push_back(hitp.stripTIBLayersWithMeasurement());
1903         t_trackHitOutMeasTID->push_back(hitp.stripTIDLayersWithMeasurement());
1904         t_trackOutPosOutHitDr->push_back(trackOutPosOutHitDr);
1905         t_trackL->push_back(trackL);
1906 
1907         t_maxNearP31x31->push_back(maxNearP31x31);
1908         t_maxNearP21x21->push_back(maxNearP21x21);
1909 
1910         t_ecalSpike11x11->push_back(e11x11P.second);
1911         t_e7x7->push_back(e7x7P.first);
1912         t_e9x9->push_back(e9x9P.first);
1913         t_e11x11->push_back(e11x11P.first);
1914         t_e15x15->push_back(e15x15P.first);
1915 
1916         t_e7x7_10Sig->push_back(e7x7_10SigP.first);
1917         t_e9x9_10Sig->push_back(e9x9_10SigP.first);
1918         t_e11x11_10Sig->push_back(e11x11_10SigP.first);
1919         t_e15x15_10Sig->push_back(e15x15_10SigP.first);
1920         t_e7x7_15Sig->push_back(e7x7_15SigP.first);
1921         t_e9x9_15Sig->push_back(e9x9_15SigP.first);
1922         t_e11x11_15Sig->push_back(e11x11_15SigP.first);
1923         t_e15x15_15Sig->push_back(e15x15_15SigP.first);
1924         t_e7x7_20Sig->push_back(e7x7_20SigP.first);
1925         t_e9x9_20Sig->push_back(e9x9_20SigP.first);
1926         t_e11x11_20Sig->push_back(e11x11_20SigP.first);
1927         t_e15x15_20Sig->push_back(e15x15_20SigP.first);
1928         t_e7x7_25Sig->push_back(e7x7_25SigP.first);
1929         t_e9x9_25Sig->push_back(e9x9_25SigP.first);
1930         t_e11x11_25Sig->push_back(e11x11_25SigP.first);
1931         t_e15x15_25Sig->push_back(e15x15_25SigP.first);
1932         t_e7x7_30Sig->push_back(e7x7_30SigP.first);
1933         t_e9x9_30Sig->push_back(e9x9_30SigP.first);
1934         t_e11x11_30Sig->push_back(e11x11_30SigP.first);
1935         t_e15x15_30Sig->push_back(e15x15_30SigP.first);
1936 
1937         if (doMC_) {
1938           t_esim7x7->push_back(simInfo7x7.eTotal);
1939           t_esim9x9->push_back(simInfo9x9.eTotal);
1940           t_esim11x11->push_back(simInfo11x11.eTotal);
1941           t_esim15x15->push_back(simInfo15x15.eTotal);
1942 
1943           t_esim7x7Matched->push_back(simInfo7x7.eMatched);
1944           t_esim9x9Matched->push_back(simInfo9x9.eMatched);
1945           t_esim11x11Matched->push_back(simInfo11x11.eMatched);
1946           t_esim15x15Matched->push_back(simInfo15x15.eMatched);
1947 
1948           t_esim7x7Rest->push_back(simInfo7x7.eRest);
1949           t_esim9x9Rest->push_back(simInfo9x9.eRest);
1950           t_esim11x11Rest->push_back(simInfo11x11.eRest);
1951           t_esim15x15Rest->push_back(simInfo15x15.eRest);
1952 
1953           t_esim7x7Photon->push_back(simInfo7x7.eGamma);
1954           t_esim9x9Photon->push_back(simInfo9x9.eGamma);
1955           t_esim11x11Photon->push_back(simInfo11x11.eGamma);
1956           t_esim15x15Photon->push_back(simInfo15x15.eGamma);
1957 
1958           t_esim7x7NeutHad->push_back(simInfo7x7.eNeutralHad);
1959           t_esim9x9NeutHad->push_back(simInfo9x9.eNeutralHad);
1960           t_esim11x11NeutHad->push_back(simInfo11x11.eNeutralHad);
1961           t_esim15x15NeutHad->push_back(simInfo15x15.eNeutralHad);
1962 
1963           t_esim7x7CharHad->push_back(simInfo7x7.eChargedHad);
1964           t_esim9x9CharHad->push_back(simInfo9x9.eChargedHad);
1965           t_esim11x11CharHad->push_back(simInfo11x11.eChargedHad);
1966           t_esim15x15CharHad->push_back(simInfo15x15.eChargedHad);
1967 
1968           t_trkEcalEne->push_back(trkEcalEne);
1969           t_simTrackP->push_back(simTrackP);
1970           t_esimPdgId->push_back(simInfo11x11.pdgMatched);
1971         }
1972 
1973         t_maxNearHcalP3x3->push_back(maxNearHcalP3x3);
1974         t_maxNearHcalP5x5->push_back(maxNearHcalP5x5);
1975         t_maxNearHcalP7x7->push_back(maxNearHcalP7x7);
1976 
1977         t_h3x3->push_back(h3x3);
1978         t_h5x5->push_back(h5x5);
1979         t_h7x7->push_back(h7x7);
1980         t_h3x3Sig->push_back(h3x3Sig);
1981         t_h5x5Sig->push_back(h5x5Sig);
1982         t_h7x7Sig->push_back(h7x7Sig);
1983 
1984         t_infoHcal->push_back(trkDetItr->okHCAL);
1985         if (doMC_) {
1986           t_trkHcalEne->push_back(hcalScale * trkHcalEne);
1987 
1988           t_hsim3x3->push_back(hcalScale * hsimInfo3x3.eTotal);
1989           t_hsim5x5->push_back(hcalScale * hsimInfo5x5.eTotal);
1990           t_hsim7x7->push_back(hcalScale * hsimInfo7x7.eTotal);
1991 
1992           t_hsim3x3Matched->push_back(hcalScale * hsimInfo3x3.eMatched);
1993           t_hsim5x5Matched->push_back(hcalScale * hsimInfo5x5.eMatched);
1994           t_hsim7x7Matched->push_back(hcalScale * hsimInfo7x7.eMatched);
1995 
1996           t_hsim3x3Rest->push_back(hcalScale * hsimInfo3x3.eRest);
1997           t_hsim5x5Rest->push_back(hcalScale * hsimInfo5x5.eRest);
1998           t_hsim7x7Rest->push_back(hcalScale * hsimInfo7x7.eRest);
1999 
2000           t_hsim3x3Photon->push_back(hcalScale * hsimInfo3x3.eGamma);
2001           t_hsim5x5Photon->push_back(hcalScale * hsimInfo5x5.eGamma);
2002           t_hsim7x7Photon->push_back(hcalScale * hsimInfo7x7.eGamma);
2003 
2004           t_hsim3x3NeutHad->push_back(hcalScale * hsimInfo3x3.eNeutralHad);
2005           t_hsim5x5NeutHad->push_back(hcalScale * hsimInfo5x5.eNeutralHad);
2006           t_hsim7x7NeutHad->push_back(hcalScale * hsimInfo7x7.eNeutralHad);
2007 
2008           t_hsim3x3CharHad->push_back(hcalScale * hsimInfo3x3.eChargedHad);
2009           t_hsim5x5CharHad->push_back(hcalScale * hsimInfo5x5.eChargedHad);
2010           t_hsim7x7CharHad->push_back(hcalScale * hsimInfo7x7.eChargedHad);
2011         }
2012 
2013       }  // if loosely isolated track
2014     }    // check p1/eta
2015   }      // loop over track collection
2016 
2017   if (haveIsoTrack)
2018     tree_->Fill();
2019 }
2020 
2021 // ----- method called once each job just before starting event loop ----
2022 void IsolatedTracksNxN::beginJob() {
2023   nEventProc_ = 0;
2024   nbad_ = 0;
2025   initL1_ = false;
2026   double tempgen_TH[NPBins + 1] = {
2027       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};
2028 
2029   for (unsigned int i = 0; i < NPBins + 1; i++)
2030     genPartPBins[i] = tempgen_TH[i];
2031 
2032   double tempgen_Eta[NEtaBins + 1] = {0.0, 1.131, 1.653, 2.172};
2033 
2034   for (unsigned int i = 0; i < NEtaBins + 1; i++)
2035     genPartEtaBins[i] = tempgen_Eta[i];
2036 
2037   bookHistograms();
2038 }
2039 
2040 // ----- method called once each job just after ending the event loop ----
2041 void IsolatedTracksNxN::endJob() {
2042   if (L1TriggerAlgoInfo_) {
2043     std::map<std::pair<unsigned int, std::string>, int>::iterator itr;
2044     for (itr = l1AlgoMap_.begin(); itr != l1AlgoMap_.end(); itr++) {
2045       edm::LogVerbatim("IsoTrack") << " ****endjob**** " << (itr->first).first << " " << (itr->first).second << " "
2046                                    << itr->second;
2047       int ibin = (itr->first).first;
2048       TString name((itr->first).second);
2049       h_L1AlgoNames->GetXaxis()->SetBinLabel(ibin + 1, name);
2050     }
2051     edm::LogVerbatim("IsoTrack") << "Number of Events Processed " << nEventProc_;
2052   }
2053 }
2054 
2055 //========================================================================================================
2056 
2057 void IsolatedTracksNxN::clearTreeVectors() {
2058   t_PVx->clear();
2059   t_PVy->clear();
2060   t_PVz->clear();
2061   t_PVisValid->clear();
2062   t_PVndof->clear();
2063   t_PVNTracks->clear();
2064   t_PVNTracksWt->clear();
2065   t_PVTracksSumPt->clear();
2066   t_PVTracksSumPtWt->clear();
2067   t_PVNTracksHP->clear();
2068   t_PVNTracksHPWt->clear();
2069   t_PVTracksSumPtHP->clear();
2070   t_PVTracksSumPtHPWt->clear();
2071 
2072   for (int i = 0; i < 128; i++)
2073     t_L1Decision[i] = 0;
2074   t_L1AlgoNames->clear();
2075   t_L1PreScale->clear();
2076 
2077   t_L1CenJetPt->clear();
2078   t_L1CenJetEta->clear();
2079   t_L1CenJetPhi->clear();
2080   t_L1FwdJetPt->clear();
2081   t_L1FwdJetEta->clear();
2082   t_L1FwdJetPhi->clear();
2083   t_L1TauJetPt->clear();
2084   t_L1TauJetEta->clear();
2085   t_L1TauJetPhi->clear();
2086   t_L1MuonPt->clear();
2087   t_L1MuonEta->clear();
2088   t_L1MuonPhi->clear();
2089   t_L1IsoEMPt->clear();
2090   t_L1IsoEMEta->clear();
2091   t_L1IsoEMPhi->clear();
2092   t_L1NonIsoEMPt->clear();
2093   t_L1NonIsoEMEta->clear();
2094   t_L1NonIsoEMPhi->clear();
2095   t_L1METPt->clear();
2096   t_L1METEta->clear();
2097   t_L1METPhi->clear();
2098 
2099   t_jetPt->clear();
2100   t_jetEta->clear();
2101   t_jetPhi->clear();
2102   t_nTrksJetCalo->clear();
2103   t_nTrksJetVtx->clear();
2104 
2105   t_trackPAll->clear();
2106   t_trackEtaAll->clear();
2107   t_trackPhiAll->clear();
2108   t_trackPdgIdAll->clear();
2109   t_trackPtAll->clear();
2110   t_trackDxyAll->clear();
2111   t_trackDzAll->clear();
2112   t_trackDxyPVAll->clear();
2113   t_trackDzPVAll->clear();
2114   t_trackChiSqAll->clear();
2115 
2116   t_trackP->clear();
2117   t_trackPt->clear();
2118   t_trackEta->clear();
2119   t_trackPhi->clear();
2120   t_trackEcalEta->clear();
2121   t_trackEcalPhi->clear();
2122   t_trackHcalEta->clear();
2123   t_trackHcalPhi->clear();
2124   t_NLayersCrossed->clear();
2125   t_trackNOuterHits->clear();
2126   t_trackDxy->clear();
2127   t_trackDxyBS->clear();
2128   t_trackDz->clear();
2129   t_trackDzBS->clear();
2130   t_trackDxyPV->clear();
2131   t_trackDzPV->clear();
2132   t_trackChiSq->clear();
2133   t_trackPVIdx->clear();
2134   t_trackHitsTOB->clear();
2135   t_trackHitsTEC->clear();
2136   t_trackHitInMissTOB->clear();
2137   t_trackHitInMissTEC->clear();
2138   t_trackHitInMissTIB->clear();
2139   t_trackHitInMissTID->clear();
2140   t_trackHitInMissTIBTID->clear();
2141   t_trackHitOutMissTOB->clear();
2142   t_trackHitOutMissTEC->clear();
2143   t_trackHitOutMissTIB->clear();
2144   t_trackHitOutMissTID->clear();
2145   t_trackHitOutMissTOBTEC->clear();
2146   t_trackHitInMeasTOB->clear();
2147   t_trackHitInMeasTEC->clear();
2148   t_trackHitInMeasTIB->clear();
2149   t_trackHitInMeasTID->clear();
2150   t_trackHitOutMeasTOB->clear();
2151   t_trackHitOutMeasTEC->clear();
2152   t_trackHitOutMeasTIB->clear();
2153   t_trackHitOutMeasTID->clear();
2154   t_trackOutPosOutHitDr->clear();
2155   t_trackL->clear();
2156 
2157   t_maxNearP31x31->clear();
2158   t_maxNearP21x21->clear();
2159 
2160   t_ecalSpike11x11->clear();
2161   t_e7x7->clear();
2162   t_e9x9->clear();
2163   t_e11x11->clear();
2164   t_e15x15->clear();
2165 
2166   t_e7x7_10Sig->clear();
2167   t_e9x9_10Sig->clear();
2168   t_e11x11_10Sig->clear();
2169   t_e15x15_10Sig->clear();
2170   t_e7x7_15Sig->clear();
2171   t_e9x9_15Sig->clear();
2172   t_e11x11_15Sig->clear();
2173   t_e15x15_15Sig->clear();
2174   t_e7x7_20Sig->clear();
2175   t_e9x9_20Sig->clear();
2176   t_e11x11_20Sig->clear();
2177   t_e15x15_20Sig->clear();
2178   t_e7x7_25Sig->clear();
2179   t_e9x9_25Sig->clear();
2180   t_e11x11_25Sig->clear();
2181   t_e15x15_25Sig->clear();
2182   t_e7x7_30Sig->clear();
2183   t_e9x9_30Sig->clear();
2184   t_e11x11_30Sig->clear();
2185   t_e15x15_30Sig->clear();
2186 
2187   if (doMC_) {
2188     t_simTrackP->clear();
2189     t_esimPdgId->clear();
2190     t_trkEcalEne->clear();
2191 
2192     t_esim7x7->clear();
2193     t_esim9x9->clear();
2194     t_esim11x11->clear();
2195     t_esim15x15->clear();
2196 
2197     t_esim7x7Matched->clear();
2198     t_esim9x9Matched->clear();
2199     t_esim11x11Matched->clear();
2200     t_esim15x15Matched->clear();
2201 
2202     t_esim7x7Rest->clear();
2203     t_esim9x9Rest->clear();
2204     t_esim11x11Rest->clear();
2205     t_esim15x15Rest->clear();
2206 
2207     t_esim7x7Photon->clear();
2208     t_esim9x9Photon->clear();
2209     t_esim11x11Photon->clear();
2210     t_esim15x15Photon->clear();
2211 
2212     t_esim7x7NeutHad->clear();
2213     t_esim9x9NeutHad->clear();
2214     t_esim11x11NeutHad->clear();
2215     t_esim15x15NeutHad->clear();
2216 
2217     t_esim7x7CharHad->clear();
2218     t_esim9x9CharHad->clear();
2219     t_esim11x11CharHad->clear();
2220     t_esim15x15CharHad->clear();
2221   }
2222 
2223   t_maxNearHcalP3x3->clear();
2224   t_maxNearHcalP5x5->clear();
2225   t_maxNearHcalP7x7->clear();
2226 
2227   t_h3x3->clear();
2228   t_h5x5->clear();
2229   t_h7x7->clear();
2230   t_h3x3Sig->clear();
2231   t_h5x5Sig->clear();
2232   t_h7x7Sig->clear();
2233 
2234   t_infoHcal->clear();
2235 
2236   if (doMC_) {
2237     t_trkHcalEne->clear();
2238 
2239     t_hsim3x3->clear();
2240     t_hsim5x5->clear();
2241     t_hsim7x7->clear();
2242     t_hsim3x3Matched->clear();
2243     t_hsim5x5Matched->clear();
2244     t_hsim7x7Matched->clear();
2245     t_hsim3x3Rest->clear();
2246     t_hsim5x5Rest->clear();
2247     t_hsim7x7Rest->clear();
2248     t_hsim3x3Photon->clear();
2249     t_hsim5x5Photon->clear();
2250     t_hsim7x7Photon->clear();
2251     t_hsim3x3NeutHad->clear();
2252     t_hsim5x5NeutHad->clear();
2253     t_hsim7x7NeutHad->clear();
2254     t_hsim3x3CharHad->clear();
2255     t_hsim5x5CharHad->clear();
2256     t_hsim7x7CharHad->clear();
2257   }
2258 }
2259 
2260 void IsolatedTracksNxN::bookHistograms() {
2261   char hname[100], htit[100];
2262 
2263   edm::Service<TFileService> fs;
2264   TFileDirectory dir = fs->mkdir("nearMaxTrackP");
2265 
2266   for (unsigned int ieta = 0; ieta < NEtaBins; ieta++) {
2267     double lowEta = -5.0, highEta = 5.0;
2268     lowEta = genPartEtaBins[ieta];
2269     highEta = genPartEtaBins[ieta + 1];
2270 
2271     for (unsigned int ipt = 0; ipt < NPBins; ipt++) {
2272       double lowP = 0.0, highP = 300.0;
2273       lowP = genPartPBins[ipt];
2274       highP = genPartPBins[ipt + 1];
2275       sprintf(hname, "h_maxNearP31x31_ptBin%i_etaBin%i", ipt, ieta);
2276       sprintf(htit, "maxNearP in 31x31 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP);
2277       h_maxNearP31x31[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
2278       h_maxNearP31x31[ipt][ieta]->Sumw2();
2279       sprintf(hname, "h_maxNearP25x25_ptBin%i_etaBin%i", ipt, ieta);
2280       sprintf(htit, "maxNearP in 25x25 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP);
2281       h_maxNearP25x25[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
2282       h_maxNearP25x25[ipt][ieta]->Sumw2();
2283       sprintf(hname, "h_maxNearP21x21_ptBin%i_etaBin%i", ipt, ieta);
2284       sprintf(htit, "maxNearP in 21x21 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP);
2285       h_maxNearP21x21[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
2286       h_maxNearP21x21[ipt][ieta]->Sumw2();
2287       sprintf(hname, "h_maxNearP15x15_ptBin%i_etaBin%i", ipt, ieta);
2288       sprintf(htit, "maxNearP in 15x15 (%3.2f<|#eta|<%3.2f), (%2.0f<trkP<%3.0f)", lowEta, highEta, lowP, highP);
2289       h_maxNearP15x15[ipt][ieta] = dir.make<TH1F>(hname, htit, 220, -2.0, 20.0);
2290       h_maxNearP15x15[ipt][ieta]->Sumw2();
2291     }
2292   }
2293 
2294   h_L1AlgoNames = fs->make<TH1I>("h_L1AlgoNames", "h_L1AlgoNames:Bin Labels", 128, -0.5, 127.5);
2295 
2296   // Reconstructed Tracks
2297 
2298   h_PVTracksWt = fs->make<TH1F>("h_PVTracksWt", "h_PVTracksWt", 600, -0.1, 1.1);
2299 
2300   h_nTracks = fs->make<TH1F>("h_nTracks", "h_nTracks", 1000, -0.5, 999.5);
2301 
2302   sprintf(hname, "h_recEtaPt_0");
2303   sprintf(htit, "h_recEtaPt (all tracks Eta vs pT)");
2304   h_recEtaPt_0 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
2305 
2306   sprintf(hname, "h_recEtaP_0");
2307   sprintf(htit, "h_recEtaP (all tracks Eta vs pT)");
2308   h_recEtaP_0 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
2309 
2310   h_recPt_0 = fs->make<TH1F>("h_recPt_0", "Pt (all tracks)", 15, genPartPBins);
2311   h_recP_0 = fs->make<TH1F>("h_recP_0", "P  (all tracks)", 15, genPartPBins);
2312   h_recEta_0 = fs->make<TH1F>("h_recEta_0", "Eta (all tracks)", 60, -3.0, 3.0);
2313   h_recPhi_0 = fs->make<TH1F>("h_recPhi_0", "Phi (all tracks)", 100, -3.2, 3.2);
2314   //-------------------------
2315   sprintf(hname, "h_recEtaPt_1");
2316   sprintf(htit, "h_recEtaPt (all good tracks Eta vs pT)");
2317   h_recEtaPt_1 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
2318 
2319   sprintf(hname, "h_recEtaP_1");
2320   sprintf(htit, "h_recEtaP (all good tracks Eta vs pT)");
2321   h_recEtaP_1 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
2322 
2323   h_recPt_1 = fs->make<TH1F>("h_recPt_1", "Pt (all good tracks)", 15, genPartPBins);
2324   h_recP_1 = fs->make<TH1F>("h_recP_1", "P  (all good tracks)", 15, genPartPBins);
2325   h_recEta_1 = fs->make<TH1F>("h_recEta_1", "Eta (all good tracks)", 60, -3.0, 3.0);
2326   h_recPhi_1 = fs->make<TH1F>("h_recPhi_1", "Phi (all good tracks)", 100, -3.2, 3.2);
2327   //-------------------------
2328   sprintf(hname, "h_recEtaPt_2");
2329   sprintf(htit, "h_recEtaPt (charge isolation Eta vs pT)");
2330   h_recEtaPt_2 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
2331 
2332   sprintf(hname, "h_recEtaP_2");
2333   sprintf(htit, "h_recEtaP (charge isolation Eta vs pT)");
2334   h_recEtaP_2 = fs->make<TH2F>(hname, htit, 30, -3.0, 3.0, 15, genPartPBins);
2335 
2336   h_recPt_2 = fs->make<TH1F>("h_recPt_2", "Pt (charge isolation)", 15, genPartPBins);
2337   h_recP_2 = fs->make<TH1F>("h_recP_2", "P  (charge isolation)", 15, genPartPBins);
2338   h_recEta_2 = fs->make<TH1F>("h_recEta_2", "Eta (charge isolation)", 60, -3.0, 3.0);
2339   h_recPhi_2 = fs->make<TH1F>("h_recPhi_2", "Phi (charge isolation)", 100, -3.2, 3.2);
2340 
2341   tree_ = fs->make<TTree>("tree", "tree");
2342   tree_->SetAutoSave(10000);
2343 
2344   tree_->Branch("t_EvtNo", &t_EvtNo, "t_EvtNo/I");
2345   tree_->Branch("t_RunNo", &t_RunNo, "t_RunNo/I");
2346   tree_->Branch("t_Lumi", &t_Lumi, "t_Lumi/I");
2347   tree_->Branch("t_Bunch", &t_Bunch, "t_Bunch/I");
2348 
2349   t_PVx = new std::vector<double>();
2350   t_PVy = new std::vector<double>();
2351   t_PVz = new std::vector<double>();
2352   t_PVisValid = new std::vector<int>();
2353   t_PVndof = new std::vector<int>();
2354   t_PVNTracks = new std::vector<int>();
2355   t_PVNTracksWt = new std::vector<int>();
2356   t_PVTracksSumPt = new std::vector<double>();
2357   t_PVTracksSumPtWt = new std::vector<double>();
2358   t_PVNTracksHP = new std::vector<int>();
2359   t_PVNTracksHPWt = new std::vector<int>();
2360   t_PVTracksSumPtHP = new std::vector<double>();
2361   t_PVTracksSumPtHPWt = new std::vector<double>();
2362 
2363   tree_->Branch("PVx", "std::vector<double>", &t_PVx);
2364   tree_->Branch("PVy", "std::vector<double>", &t_PVy);
2365   tree_->Branch("PVz", "std::vector<double>", &t_PVz);
2366   tree_->Branch("PVisValid", "std::vector<int>", &t_PVisValid);
2367   tree_->Branch("PVndof", "std::vector<int>", &t_PVndof);
2368   tree_->Branch("PVNTracks", "std::vector<int>", &t_PVNTracks);
2369   tree_->Branch("PVNTracksWt", "std::vector<int>", &t_PVNTracksWt);
2370   tree_->Branch("t_PVTracksSumPt", "std::vector<double>", &t_PVTracksSumPt);
2371   tree_->Branch("t_PVTracksSumPtWt", "std::vector<double>", &t_PVTracksSumPtWt);
2372   tree_->Branch("PVNTracksHP", "std::vector<int>", &t_PVNTracksHP);
2373   tree_->Branch("PVNTracksHPWt", "std::vector<int>", &t_PVNTracksHPWt);
2374   tree_->Branch("t_PVTracksSumPtHP", "std::vector<double>", &t_PVTracksSumPtHP);
2375   tree_->Branch("t_PVTracksSumPtHPWt", "std::vector<double>", &t_PVTracksSumPtHPWt);
2376 
2377   //----- L1Trigger information
2378   for (int i = 0; i < 128; i++)
2379     t_L1Decision[i] = 0;
2380   t_L1AlgoNames = new std::vector<std::string>();
2381   t_L1PreScale = new std::vector<int>();
2382   t_L1CenJetPt = new std::vector<double>();
2383   t_L1CenJetEta = new std::vector<double>();
2384   t_L1CenJetPhi = new std::vector<double>();
2385   t_L1FwdJetPt = new std::vector<double>();
2386   t_L1FwdJetEta = new std::vector<double>();
2387   t_L1FwdJetPhi = new std::vector<double>();
2388   t_L1TauJetPt = new std::vector<double>();
2389   t_L1TauJetEta = new std::vector<double>();
2390   t_L1TauJetPhi = new std::vector<double>();
2391   t_L1MuonPt = new std::vector<double>();
2392   t_L1MuonEta = new std::vector<double>();
2393   t_L1MuonPhi = new std::vector<double>();
2394   t_L1IsoEMPt = new std::vector<double>();
2395   t_L1IsoEMEta = new std::vector<double>();
2396   t_L1IsoEMPhi = new std::vector<double>();
2397   t_L1NonIsoEMPt = new std::vector<double>();
2398   t_L1NonIsoEMEta = new std::vector<double>();
2399   t_L1NonIsoEMPhi = new std::vector<double>();
2400   t_L1METPt = new std::vector<double>();
2401   t_L1METEta = new std::vector<double>();
2402   t_L1METPhi = new std::vector<double>();
2403 
2404   tree_->Branch("t_L1Decision", t_L1Decision, "t_L1Decision[128]/I");
2405   tree_->Branch("t_L1AlgoNames", "std::vector<string>", &t_L1AlgoNames);
2406   tree_->Branch("t_L1PreScale", "std::vector<int>", &t_L1PreScale);
2407   tree_->Branch("t_L1CenJetPt", "std::vector<double>", &t_L1CenJetPt);
2408   tree_->Branch("t_L1CenJetEta", "std::vector<double>", &t_L1CenJetEta);
2409   tree_->Branch("t_L1CenJetPhi", "std::vector<double>", &t_L1CenJetPhi);
2410   tree_->Branch("t_L1FwdJetPt", "std::vector<double>", &t_L1FwdJetPt);
2411   tree_->Branch("t_L1FwdJetEta", "std::vector<double>", &t_L1FwdJetEta);
2412   tree_->Branch("t_L1FwdJetPhi", "std::vector<double>", &t_L1FwdJetPhi);
2413   tree_->Branch("t_L1TauJetPt", "std::vector<double>", &t_L1TauJetPt);
2414   tree_->Branch("t_L1TauJetEta", "std::vector<double>", &t_L1TauJetEta);
2415   tree_->Branch("t_L1TauJetPhi", "std::vector<double>", &t_L1TauJetPhi);
2416   tree_->Branch("t_L1MuonPt", "std::vector<double>", &t_L1MuonPt);
2417   tree_->Branch("t_L1MuonEta", "std::vector<double>", &t_L1MuonEta);
2418   tree_->Branch("t_L1MuonPhi", "std::vector<double>", &t_L1MuonPhi);
2419   tree_->Branch("t_L1IsoEMPt", "std::vector<double>", &t_L1IsoEMPt);
2420   tree_->Branch("t_L1IsoEMEta", "std::vector<double>", &t_L1IsoEMEta);
2421   tree_->Branch("t_L1IsoEMPhi", "std::vector<double>", &t_L1IsoEMPhi);
2422   tree_->Branch("t_L1NonIsoEMPt", "std::vector<double>", &t_L1NonIsoEMPt);
2423   tree_->Branch("t_L1NonIsoEMEta", "std::vector<double>", &t_L1NonIsoEMEta);
2424   tree_->Branch("t_L1NonIsoEMPhi", "std::vector<double>", &t_L1NonIsoEMPhi);
2425   tree_->Branch("t_L1METPt", "std::vector<double>", &t_L1METPt);
2426   tree_->Branch("t_L1METEta", "std::vector<double>", &t_L1METEta);
2427   tree_->Branch("t_L1METPhi", "std::vector<double>", &t_L1METPhi);
2428 
2429   t_jetPt = new std::vector<double>();
2430   t_jetEta = new std::vector<double>();
2431   t_jetPhi = new std::vector<double>();
2432   t_nTrksJetCalo = new std::vector<double>();
2433   t_nTrksJetVtx = new std::vector<double>();
2434   tree_->Branch("t_jetPt", "std::vector<double>", &t_jetPt);
2435   tree_->Branch("t_jetEta", "std::vector<double>", &t_jetEta);
2436   tree_->Branch("t_jetPhi", "std::vector<double>", &t_jetPhi);
2437   tree_->Branch("t_nTrksJetCalo", "std::vector<double>", &t_nTrksJetCalo);
2438   tree_->Branch("t_nTrksJetVtx", "std::vector<double>", &t_nTrksJetVtx);
2439 
2440   t_trackPAll = new std::vector<double>();
2441   t_trackEtaAll = new std::vector<double>();
2442   t_trackPhiAll = new std::vector<double>();
2443   t_trackPdgIdAll = new std::vector<double>();
2444   t_trackPtAll = new std::vector<double>();
2445   t_trackDxyAll = new std::vector<double>();
2446   t_trackDzAll = new std::vector<double>();
2447   t_trackDxyPVAll = new std::vector<double>();
2448   t_trackDzPVAll = new std::vector<double>();
2449   t_trackChiSqAll = new std::vector<double>();
2450   tree_->Branch("t_trackPAll", "std::vector<double>", &t_trackPAll);
2451   tree_->Branch("t_trackPhiAll", "std::vector<double>", &t_trackPhiAll);
2452   tree_->Branch("t_trackEtaAll", "std::vector<double>", &t_trackEtaAll);
2453   tree_->Branch("t_trackPtAll", "std::vector<double>", &t_trackPtAll);
2454   tree_->Branch("t_trackDxyAll", "std::vector<double>", &t_trackDxyAll);
2455   tree_->Branch("t_trackDzAll", "std::vector<double>", &t_trackDzAll);
2456   tree_->Branch("t_trackDxyPVAll", "std::vector<double>", &t_trackDxyPVAll);
2457   tree_->Branch("t_trackDzPVAll", "std::vector<double>", &t_trackDzPVAll);
2458   tree_->Branch("t_trackChiSqAll", "std::vector<double>", &t_trackChiSqAll);
2459   //tree_->Branch("t_trackPdgIdAll",     "std::vector<double>", &t_trackPdgIdAll);
2460 
2461   t_trackP = new std::vector<double>();
2462   t_trackPt = new std::vector<double>();
2463   t_trackEta = new std::vector<double>();
2464   t_trackPhi = new std::vector<double>();
2465   t_trackEcalEta = new std::vector<double>();
2466   t_trackEcalPhi = new std::vector<double>();
2467   t_trackHcalEta = new std::vector<double>();
2468   t_trackHcalPhi = new std::vector<double>();
2469   t_trackNOuterHits = new std::vector<int>();
2470   t_NLayersCrossed = new std::vector<int>();
2471   t_trackDxy = new std::vector<double>();
2472   t_trackDxyBS = new std::vector<double>();
2473   t_trackDz = new std::vector<double>();
2474   t_trackDzBS = new std::vector<double>();
2475   t_trackDxyPV = new std::vector<double>();
2476   t_trackDzPV = new std::vector<double>();
2477   t_trackPVIdx = new std::vector<int>();
2478   t_trackChiSq = new std::vector<double>();
2479   t_trackHitsTOB = new std::vector<int>();
2480   t_trackHitsTEC = new std::vector<int>();
2481   t_trackHitInMissTOB = new std::vector<int>();
2482   t_trackHitInMissTEC = new std::vector<int>();
2483   t_trackHitInMissTIB = new std::vector<int>();
2484   t_trackHitInMissTID = new std::vector<int>();
2485   t_trackHitOutMissTOB = new std::vector<int>();
2486   t_trackHitOutMissTEC = new std::vector<int>();
2487   t_trackHitOutMissTIB = new std::vector<int>();
2488   t_trackHitOutMissTID = new std::vector<int>();
2489   t_trackHitInMissTIBTID = new std::vector<int>();
2490   t_trackHitOutMissTOB = new std::vector<int>();
2491   t_trackHitOutMissTEC = new std::vector<int>();
2492   t_trackHitOutMissTIB = new std::vector<int>();
2493   t_trackHitOutMissTID = new std::vector<int>();
2494   t_trackHitOutMissTOBTEC = new std::vector<int>();
2495   t_trackHitInMeasTOB = new std::vector<int>();
2496   t_trackHitInMeasTEC = new std::vector<int>();
2497   t_trackHitInMeasTIB = new std::vector<int>();
2498   t_trackHitInMeasTID = new std::vector<int>();
2499   t_trackHitOutMeasTOB = new std::vector<int>();
2500   t_trackHitOutMeasTEC = new std::vector<int>();
2501   t_trackHitOutMeasTIB = new std::vector<int>();
2502   t_trackHitOutMeasTID = new std::vector<int>();
2503   t_trackOutPosOutHitDr = new std::vector<double>();
2504   t_trackL = new std::vector<double>();
2505 
2506   tree_->Branch("t_trackP", "std::vector<double>", &t_trackP);
2507   tree_->Branch("t_trackPt", "std::vector<double>", &t_trackPt);
2508   tree_->Branch("t_trackEta", "std::vector<double>", &t_trackEta);
2509   tree_->Branch("t_trackPhi", "std::vector<double>", &t_trackPhi);
2510   tree_->Branch("t_trackEcalEta", "std::vector<double>", &t_trackEcalEta);
2511   tree_->Branch("t_trackEcalPhi", "std::vector<double>", &t_trackEcalPhi);
2512   tree_->Branch("t_trackHcalEta", "std::vector<double>", &t_trackHcalEta);
2513   tree_->Branch("t_trackHcalPhi", "std::vector<double>", &t_trackHcalPhi);
2514 
2515   tree_->Branch("t_trackNOuterHits", "std::vector<int>", &t_trackNOuterHits);
2516   tree_->Branch("t_NLayersCrossed", "std::vector<int>", &t_NLayersCrossed);
2517   tree_->Branch("t_trackHitsTOB", "std::vector<int>", &t_trackHitsTOB);
2518   tree_->Branch("t_trackHitsTEC", "std::vector<int>", &t_trackHitsTEC);
2519   tree_->Branch("t_trackHitInMissTOB", "std::vector<int>", &t_trackHitInMissTOB);
2520   tree_->Branch("t_trackHitInMissTEC", "std::vector<int>", &t_trackHitInMissTEC);
2521   tree_->Branch("t_trackHitInMissTIB", "std::vector<int>", &t_trackHitInMissTIB);
2522   tree_->Branch("t_trackHitInMissTID", "std::vector<int>", &t_trackHitInMissTID);
2523   tree_->Branch("t_trackHitInMissTIBTID", "std::vector<int>", &t_trackHitInMissTIBTID);
2524   tree_->Branch("t_trackHitOutMissTOB", "std::vector<int>", &t_trackHitOutMissTOB);
2525   tree_->Branch("t_trackHitOutMissTEC", "std::vector<int>", &t_trackHitOutMissTEC);
2526   tree_->Branch("t_trackHitOutMissTIB", "std::vector<int>", &t_trackHitOutMissTIB);
2527   tree_->Branch("t_trackHitOutMissTID", "std::vector<int>", &t_trackHitOutMissTID);
2528   tree_->Branch("t_trackHitOutMissTOBTEC", "std::vector<int>", &t_trackHitOutMissTOBTEC);
2529   tree_->Branch("t_trackHitInMeasTOB", "std::vector<int>", &t_trackHitInMeasTOB);
2530   tree_->Branch("t_trackHitInMeasTEC", "std::vector<int>", &t_trackHitInMeasTEC);
2531   tree_->Branch("t_trackHitInMeasTIB", "std::vector<int>", &t_trackHitInMeasTIB);
2532   tree_->Branch("t_trackHitInMeasTID", "std::vector<int>", &t_trackHitInMeasTID);
2533   tree_->Branch("t_trackHitOutMeasTOB", "std::vector<int>", &t_trackHitOutMeasTOB);
2534   tree_->Branch("t_trackHitOutMeasTEC", "std::vector<int>", &t_trackHitOutMeasTEC);
2535   tree_->Branch("t_trackHitOutMeasTIB", "std::vector<int>", &t_trackHitOutMeasTIB);
2536   tree_->Branch("t_trackHitOutMeasTID", "std::vector<int>", &t_trackHitOutMeasTID);
2537   tree_->Branch("t_trackOutPosOutHitDr", "std::vector<double>", &t_trackOutPosOutHitDr);
2538   tree_->Branch("t_trackL", "std::vector<double>", &t_trackL);
2539 
2540   tree_->Branch("t_trackDxy", "std::vector<double>", &t_trackDxy);
2541   tree_->Branch("t_trackDxyBS", "std::vector<double>", &t_trackDxyBS);
2542   tree_->Branch("t_trackDz", "std::vector<double>", &t_trackDz);
2543   tree_->Branch("t_trackDzBS", "std::vector<double>", &t_trackDzBS);
2544   tree_->Branch("t_trackDxyPV", "std::vector<double>", &t_trackDxyPV);
2545   tree_->Branch("t_trackDzPV", "std::vector<double>", &t_trackDzPV);
2546   tree_->Branch("t_trackChiSq", "std::vector<double>", &t_trackChiSq);
2547   tree_->Branch("t_trackPVIdx", "std::vector<int>", &t_trackPVIdx);
2548 
2549   t_maxNearP31x31 = new std::vector<double>();
2550   t_maxNearP21x21 = new std::vector<double>();
2551 
2552   tree_->Branch("t_maxNearP31x31", "std::vector<double>", &t_maxNearP31x31);
2553   tree_->Branch("t_maxNearP21x21", "std::vector<double>", &t_maxNearP21x21);
2554 
2555   t_ecalSpike11x11 = new std::vector<int>();
2556   t_e7x7 = new std::vector<double>();
2557   t_e9x9 = new std::vector<double>();
2558   t_e11x11 = new std::vector<double>();
2559   t_e15x15 = new std::vector<double>();
2560 
2561   tree_->Branch("t_ecalSpike11x11", "std::vector<int>", &t_ecalSpike11x11);
2562   tree_->Branch("t_e7x7", "std::vector<double>", &t_e7x7);
2563   tree_->Branch("t_e9x9", "std::vector<double>", &t_e9x9);
2564   tree_->Branch("t_e11x11", "std::vector<double>", &t_e11x11);
2565   tree_->Branch("t_e15x15", "std::vector<double>", &t_e15x15);
2566 
2567   t_e7x7_10Sig = new std::vector<double>();
2568   t_e9x9_10Sig = new std::vector<double>();
2569   t_e11x11_10Sig = new std::vector<double>();
2570   t_e15x15_10Sig = new std::vector<double>();
2571   t_e7x7_15Sig = new std::vector<double>();
2572   t_e9x9_15Sig = new std::vector<double>();
2573   t_e11x11_15Sig = new std::vector<double>();
2574   t_e15x15_15Sig = new std::vector<double>();
2575   t_e7x7_20Sig = new std::vector<double>();
2576   t_e9x9_20Sig = new std::vector<double>();
2577   t_e11x11_20Sig = new std::vector<double>();
2578   t_e15x15_20Sig = new std::vector<double>();
2579   t_e7x7_25Sig = new std::vector<double>();
2580   t_e9x9_25Sig = new std::vector<double>();
2581   t_e11x11_25Sig = new std::vector<double>();
2582   t_e15x15_25Sig = new std::vector<double>();
2583   t_e7x7_30Sig = new std::vector<double>();
2584   t_e9x9_30Sig = new std::vector<double>();
2585   t_e11x11_30Sig = new std::vector<double>();
2586   t_e15x15_30Sig = new std::vector<double>();
2587 
2588   tree_->Branch("t_e7x7_10Sig", "std::vector<double>", &t_e7x7_10Sig);
2589   tree_->Branch("t_e9x9_10Sig", "std::vector<double>", &t_e9x9_10Sig);
2590   tree_->Branch("t_e11x11_10Sig", "std::vector<double>", &t_e11x11_10Sig);
2591   tree_->Branch("t_e15x15_10Sig", "std::vector<double>", &t_e15x15_10Sig);
2592   tree_->Branch("t_e7x7_15Sig", "std::vector<double>", &t_e7x7_15Sig);
2593   tree_->Branch("t_e9x9_15Sig", "std::vector<double>", &t_e9x9_15Sig);
2594   tree_->Branch("t_e11x11_15Sig", "std::vector<double>", &t_e11x11_15Sig);
2595   tree_->Branch("t_e15x15_15Sig", "std::vector<double>", &t_e15x15_15Sig);
2596   tree_->Branch("t_e7x7_20Sig", "std::vector<double>", &t_e7x7_20Sig);
2597   tree_->Branch("t_e9x9_20Sig", "std::vector<double>", &t_e9x9_20Sig);
2598   tree_->Branch("t_e11x11_20Sig", "std::vector<double>", &t_e11x11_20Sig);
2599   tree_->Branch("t_e15x15_20Sig", "std::vector<double>", &t_e15x15_20Sig);
2600   tree_->Branch("t_e7x7_25Sig", "std::vector<double>", &t_e7x7_25Sig);
2601   tree_->Branch("t_e9x9_25Sig", "std::vector<double>", &t_e9x9_25Sig);
2602   tree_->Branch("t_e11x11_25Sig", "std::vector<double>", &t_e11x11_25Sig);
2603   tree_->Branch("t_e15x15_25Sig", "std::vector<double>", &t_e15x15_25Sig);
2604   tree_->Branch("t_e7x7_30Sig", "std::vector<double>", &t_e7x7_30Sig);
2605   tree_->Branch("t_e9x9_30Sig", "std::vector<double>", &t_e9x9_30Sig);
2606   tree_->Branch("t_e11x11_30Sig", "std::vector<double>", &t_e11x11_30Sig);
2607   tree_->Branch("t_e15x15_30Sig", "std::vector<double>", &t_e15x15_30Sig);
2608 
2609   if (doMC_) {
2610     t_esim7x7 = new std::vector<double>();
2611     t_esim9x9 = new std::vector<double>();
2612     t_esim11x11 = new std::vector<double>();
2613     t_esim15x15 = new std::vector<double>();
2614 
2615     t_esim7x7Matched = new std::vector<double>();
2616     t_esim9x9Matched = new std::vector<double>();
2617     t_esim11x11Matched = new std::vector<double>();
2618     t_esim15x15Matched = new std::vector<double>();
2619 
2620     t_esim7x7Rest = new std::vector<double>();
2621     t_esim9x9Rest = new std::vector<double>();
2622     t_esim11x11Rest = new std::vector<double>();
2623     t_esim15x15Rest = new std::vector<double>();
2624 
2625     t_esim7x7Photon = new std::vector<double>();
2626     t_esim9x9Photon = new std::vector<double>();
2627     t_esim11x11Photon = new std::vector<double>();
2628     t_esim15x15Photon = new std::vector<double>();
2629 
2630     t_esim7x7NeutHad = new std::vector<double>();
2631     t_esim9x9NeutHad = new std::vector<double>();
2632     t_esim11x11NeutHad = new std::vector<double>();
2633     t_esim15x15NeutHad = new std::vector<double>();
2634 
2635     t_esim7x7CharHad = new std::vector<double>();
2636     t_esim9x9CharHad = new std::vector<double>();
2637     t_esim11x11CharHad = new std::vector<double>();
2638     t_esim15x15CharHad = new std::vector<double>();
2639 
2640     t_trkEcalEne = new std::vector<double>();
2641     t_simTrackP = new std::vector<double>();
2642     t_esimPdgId = new std::vector<double>();
2643 
2644     tree_->Branch("t_esim7x7", "std::vector<double>", &t_esim7x7);
2645     tree_->Branch("t_esim9x9", "std::vector<double>", &t_esim9x9);
2646     tree_->Branch("t_esim11x11", "std::vector<double>", &t_esim11x11);
2647     tree_->Branch("t_esim15x15", "std::vector<double>", &t_esim15x15);
2648 
2649     tree_->Branch("t_esim7x7Matched", "std::vector<double>", &t_esim7x7Matched);
2650     tree_->Branch("t_esim9x9Matched", "std::vector<double>", &t_esim9x9Matched);
2651     tree_->Branch("t_esim11x11Matched", "std::vector<double>", &t_esim11x11Matched);
2652     tree_->Branch("t_esim15x15Matched", "std::vector<double>", &t_esim15x15Matched);
2653 
2654     tree_->Branch("t_esim7x7Rest", "std::vector<double>", &t_esim7x7Rest);
2655     tree_->Branch("t_esim9x9Rest", "std::vector<double>", &t_esim9x9Rest);
2656     tree_->Branch("t_esim11x11Rest", "std::vector<double>", &t_esim11x11Rest);
2657     tree_->Branch("t_esim15x15Rest", "std::vector<double>", &t_esim15x15Rest);
2658 
2659     tree_->Branch("t_esim7x7Photon", "std::vector<double>", &t_esim7x7Photon);
2660     tree_->Branch("t_esim9x9Photon", "std::vector<double>", &t_esim9x9Photon);
2661     tree_->Branch("t_esim11x11Photon", "std::vector<double>", &t_esim11x11Photon);
2662     tree_->Branch("t_esim15x15Photon", "std::vector<double>", &t_esim15x15Photon);
2663 
2664     tree_->Branch("t_esim7x7NeutHad", "std::vector<double>", &t_esim7x7NeutHad);
2665     tree_->Branch("t_esim9x9NeutHad", "std::vector<double>", &t_esim9x9NeutHad);
2666     tree_->Branch("t_esim11x11NeutHad", "std::vector<double>", &t_esim11x11NeutHad);
2667     tree_->Branch("t_esim15x15NeutHad", "std::vector<double>", &t_esim15x15NeutHad);
2668 
2669     tree_->Branch("t_esim7x7CharHad", "std::vector<double>", &t_esim7x7CharHad);
2670     tree_->Branch("t_esim9x9CharHad", "std::vector<double>", &t_esim9x9CharHad);
2671     tree_->Branch("t_esim11x11CharHad", "std::vector<double>", &t_esim11x11CharHad);
2672     tree_->Branch("t_esim15x15CharHad", "std::vector<double>", &t_esim15x15CharHad);
2673 
2674     tree_->Branch("t_trkEcalEne", "std::vector<double>", &t_trkEcalEne);
2675     tree_->Branch("t_simTrackP", "std::vector<double>", &t_simTrackP);
2676     tree_->Branch("t_esimPdgId", "std::vector<double>", &t_esimPdgId);
2677   }
2678 
2679   t_maxNearHcalP3x3 = new std::vector<double>();
2680   t_maxNearHcalP5x5 = new std::vector<double>();
2681   t_maxNearHcalP7x7 = new std::vector<double>();
2682   t_h3x3 = new std::vector<double>();
2683   t_h5x5 = new std::vector<double>();
2684   t_h7x7 = new std::vector<double>();
2685   t_h3x3Sig = new std::vector<double>();
2686   t_h5x5Sig = new std::vector<double>();
2687   t_h7x7Sig = new std::vector<double>();
2688   t_infoHcal = new std::vector<int>();
2689 
2690   if (doMC_) {
2691     t_trkHcalEne = new std::vector<double>();
2692     t_hsim3x3 = new std::vector<double>();
2693     t_hsim5x5 = new std::vector<double>();
2694     t_hsim7x7 = new std::vector<double>();
2695     t_hsim3x3Matched = new std::vector<double>();
2696     t_hsim5x5Matched = new std::vector<double>();
2697     t_hsim7x7Matched = new std::vector<double>();
2698     t_hsim3x3Rest = new std::vector<double>();
2699     t_hsim5x5Rest = new std::vector<double>();
2700     t_hsim7x7Rest = new std::vector<double>();
2701     t_hsim3x3Photon = new std::vector<double>();
2702     t_hsim5x5Photon = new std::vector<double>();
2703     t_hsim7x7Photon = new std::vector<double>();
2704     t_hsim3x3NeutHad = new std::vector<double>();
2705     t_hsim5x5NeutHad = new std::vector<double>();
2706     t_hsim7x7NeutHad = new std::vector<double>();
2707     t_hsim3x3CharHad = new std::vector<double>();
2708     t_hsim5x5CharHad = new std::vector<double>();
2709     t_hsim7x7CharHad = new std::vector<double>();
2710   }
2711 
2712   tree_->Branch("t_maxNearHcalP3x3", "std::vector<double>", &t_maxNearHcalP3x3);
2713   tree_->Branch("t_maxNearHcalP5x5", "std::vector<double>", &t_maxNearHcalP5x5);
2714   tree_->Branch("t_maxNearHcalP7x7", "std::vector<double>", &t_maxNearHcalP7x7);
2715   tree_->Branch("t_h3x3", "std::vector<double>", &t_h3x3);
2716   tree_->Branch("t_h5x5", "std::vector<double>", &t_h5x5);
2717   tree_->Branch("t_h7x7", "std::vector<double>", &t_h7x7);
2718   tree_->Branch("t_h3x3Sig", "std::vector<double>", &t_h3x3Sig);
2719   tree_->Branch("t_h5x5Sig", "std::vector<double>", &t_h5x5Sig);
2720   tree_->Branch("t_h7x7Sig", "std::vector<double>", &t_h7x7Sig);
2721   tree_->Branch("t_infoHcal", "std::vector<int>", &t_infoHcal);
2722 
2723   if (doMC_) {
2724     tree_->Branch("t_trkHcalEne", "std::vector<double>", &t_trkHcalEne);
2725     tree_->Branch("t_hsim3x3", "std::vector<double>", &t_hsim3x3);
2726     tree_->Branch("t_hsim5x5", "std::vector<double>", &t_hsim5x5);
2727     tree_->Branch("t_hsim7x7", "std::vector<double>", &t_hsim7x7);
2728     tree_->Branch("t_hsim3x3Matched", "std::vector<double>", &t_hsim3x3Matched);
2729     tree_->Branch("t_hsim5x5Matched", "std::vector<double>", &t_hsim5x5Matched);
2730     tree_->Branch("t_hsim7x7Matched", "std::vector<double>", &t_hsim7x7Matched);
2731     tree_->Branch("t_hsim3x3Rest", "std::vector<double>", &t_hsim3x3Rest);
2732     tree_->Branch("t_hsim5x5Rest", "std::vector<double>", &t_hsim5x5Rest);
2733     tree_->Branch("t_hsim7x7Rest", "std::vector<double>", &t_hsim7x7Rest);
2734     tree_->Branch("t_hsim3x3Photon", "std::vector<double>", &t_hsim3x3Photon);
2735     tree_->Branch("t_hsim5x5Photon", "std::vector<double>", &t_hsim5x5Photon);
2736     tree_->Branch("t_hsim7x7Photon", "std::vector<double>", &t_hsim7x7Photon);
2737     tree_->Branch("t_hsim3x3NeutHad", "std::vector<double>", &t_hsim3x3NeutHad);
2738     tree_->Branch("t_hsim5x5NeutHad", "std::vector<double>", &t_hsim5x5NeutHad);
2739     tree_->Branch("t_hsim7x7NeutHad", "std::vector<double>", &t_hsim7x7NeutHad);
2740     tree_->Branch("t_hsim3x3CharHad", "std::vector<double>", &t_hsim3x3CharHad);
2741     tree_->Branch("t_hsim5x5CharHad", "std::vector<double>", &t_hsim5x5CharHad);
2742     tree_->Branch("t_hsim7x7CharHad", "std::vector<double>", &t_hsim7x7CharHad);
2743   }
2744   tree_->Branch("t_nTracks", &t_nTracks, "t_nTracks/I");
2745 }
2746 
2747 void IsolatedTracksNxN::printTrack(const reco::Track *pTrack) {
2748   std::string theTrackQuality = "highPurity";
2749   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
2750 
2751   edm::LogVerbatim("IsoTrack") << " Reference Point " << pTrack->referencePoint() << "\n TrackMmentum "
2752                                << pTrack->momentum() << " (pt,eta,phi)(" << pTrack->pt() << "," << pTrack->eta() << ","
2753                                << pTrack->phi() << ")"
2754                                << " p " << pTrack->p() << "\n Normalized chi2 " << pTrack->normalizedChi2()
2755                                << "  charge " << pTrack->charge() << " qoverp() " << pTrack->qoverp() << "+-"
2756                                << pTrack->qoverpError() << " d0 " << pTrack->d0() << "\n NValidHits "
2757                                << pTrack->numberOfValidHits() << "  NLostHits " << pTrack->numberOfLostHits()
2758                                << " TrackQuality " << pTrack->qualityName(trackQuality_) << " "
2759                                << pTrack->quality(trackQuality_);
2760 
2761   if (printTrkHitPattern_) {
2762     const reco::HitPattern &p = pTrack->hitPattern();
2763 
2764     std::ostringstream st1;
2765     st1 << "default ";
2766     for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::TRACK_HITS); i++) {
2767       p.printHitPattern(reco::HitPattern::TRACK_HITS, i, std::cout);
2768     }
2769     edm::LogVerbatim("IsoTrack") << st1.str();
2770     std::ostringstream st2;
2771     st2 << "trackerMissingInnerHits() ";
2772     for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::MISSING_INNER_HITS); i++) {
2773       p.printHitPattern(reco::HitPattern::MISSING_INNER_HITS, i, st2);
2774     }
2775     edm::LogVerbatim("IsoTrack") << st2.str();
2776     std::ostringstream st3;
2777     st3 << "trackerMissingOuterHits() ";
2778     for (int i = 0; i < p.numberOfAllHits(reco::HitPattern::MISSING_OUTER_HITS); i++) {
2779       p.printHitPattern(reco::HitPattern::MISSING_OUTER_HITS, i, st3);
2780     }
2781     edm::LogVerbatim("IsoTrack") << st3.str();
2782 
2783     edm::LogVerbatim("IsoTrack") << "\n \t trackerLayersWithMeasurement() " << p.trackerLayersWithMeasurement()
2784                                  << "\n \t pixelLayersWithMeasurement() " << p.pixelLayersWithMeasurement()
2785                                  << "\n \t stripLayersWithMeasurement() " << p.stripLayersWithMeasurement()
2786                                  << "\n \t pixelBarrelLayersWithMeasurement() " << p.pixelBarrelLayersWithMeasurement()
2787                                  << "\n \t pixelEndcapLayersWithMeasurement() " << p.pixelEndcapLayersWithMeasurement()
2788                                  << "\n \t stripTIBLayersWithMeasurement() " << p.stripTIBLayersWithMeasurement()
2789                                  << "\n \t stripTIDLayersWithMeasurement() " << p.stripTIDLayersWithMeasurement()
2790                                  << "\n \t stripTOBLayersWithMeasurement() " << p.stripTOBLayersWithMeasurement()
2791                                  << "\n \t stripTECLayersWithMeasurement() " << p.stripTECLayersWithMeasurement();
2792   }
2793 }
2794 
2795 //define this as a plug-in
2796 DEFINE_FWK_MODULE(IsolatedTracksNxN);