Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-16 06:12:12

0001 // system include files
0002 #include <atomic>
0003 #include <cmath>
0004 #include <memory>
0005 #include <string>
0006 #include <vector>
0007 
0008 // Root objects
0009 #include "TROOT.h"
0010 #include "TSystem.h"
0011 #include "TFile.h"
0012 #include "TProfile.h"
0013 #include "TDirectory.h"
0014 #include "TTree.h"
0015 #include "TLorentzVector.h"
0016 #include "TInterpreter.h"
0017 
0018 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0019 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0020 #include "Calibration/IsolatedParticles/interface/eCone.h"
0021 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0022 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0023 #include "Calibration/IsolatedParticles/interface/TrackSelection.h"
0024 
0025 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0026 
0027 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0028 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0029 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0030 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0031 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0032 
0033 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0034 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0035 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0036 //Tracks
0037 #include "DataFormats/TrackReco/interface/Track.h"
0038 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0039 #include "DataFormats/TrackReco/interface/HitPattern.h"
0040 #include "DataFormats/TrackReco/interface/TrackBase.h"
0041 // Vertices
0042 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0043 #include "DataFormats/VertexReco/interface/Vertex.h"
0044 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0045 //Triggers
0046 #include "DataFormats/Common/interface/TriggerResults.h"
0047 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0048 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0049 #include "DataFormats/Luminosity/interface/LumiDetails.h"
0050 #include "DataFormats/Math/interface/LorentzVector.h"
0051 #include "DataFormats/Math/interface/deltaR.h"
0052 #include "DataFormats/L1TGlobal/interface/GlobalAlgBlk.h"
0053 #include "DataFormats/L1TGlobal/interface/GlobalExtBlk.h"
0054 #include "DataFormats/L1Trigger/interface/BXVector.h"
0055 //Generator information
0056 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0057 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
0058 
0059 #include "FWCore/Framework/interface/Frameworkfwd.h"
0060 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0061 #include "FWCore/Framework/interface/Event.h"
0062 #include "FWCore/Framework/interface/MakerMacros.h"
0063 #include "FWCore/ServiceRegistry/interface/Service.h"
0064 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0065 #include "FWCore/Common/interface/TriggerNames.h"
0066 
0067 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0068 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0069 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0070 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0071 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0072 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0073 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0074 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0075 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0076 
0077 #include "L1Trigger/L1TGlobal/interface/L1TGlobalUtil.h"
0078 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0079 
0080 #include "MagneticField/Engine/interface/MagneticField.h"
0081 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0082 
0083 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0084 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0085 
0086 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0087 
0088 //#define EDM_ML_DEBUG
0089 
0090 class HcalIsoTrackStudy : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0091 public:
0092   explicit HcalIsoTrackStudy(edm::ParameterSet const&);
0093   ~HcalIsoTrackStudy() override {}
0094 
0095   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0096 
0097 private:
0098   void analyze(edm::Event const&, edm::EventSetup const&) override;
0099   void beginJob() override;
0100   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0101   void endRun(edm::Run const&, edm::EventSetup const&) override;
0102 
0103   std::array<int, 3> fillTree(std::vector<math::XYZTLorentzVector>& vecL1,
0104                               std::vector<math::XYZTLorentzVector>& vecL3,
0105                               math::XYZPoint& leadPV,
0106                               std::vector<spr::propagatedTrackDirection>& trkCaloDirections,
0107                               std::vector<spr::propagatedTrackID>& trkCaloDets,
0108                               const CaloGeometry* geo,
0109                               const CaloTopology* topo,
0110                               const HcalTopology* theHBHETopology,
0111                               const EcalChannelStatus* theEcalChStatus,
0112                               const EcalSeverityLevelAlgo* theEcalSevlv,
0113                               edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
0114                               edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
0115                               edm::Handle<HBHERecHitCollection>& hbhe,
0116                               edm::Handle<CaloTowerCollection>& towerHandle,
0117                               const edm::Handle<reco::GenParticleCollection>& genParticles,
0118                               const HcalRespCorrs* respCorrs);
0119   double dR(math::XYZTLorentzVector&, math::XYZTLorentzVector&);
0120   double trackP(const reco::Track*, const edm::Handle<reco::GenParticleCollection>&);
0121   double trackE(const reco::Track*, const edm::Handle<reco::GenParticleCollection>&);
0122   double rhoh(const edm::Handle<CaloTowerCollection>&);
0123   void storeEnergy(int indx,
0124                    const HcalRespCorrs* respCorrs,
0125                    const std::vector<DetId>& ids,
0126                    std::vector<double>& edet,
0127                    double& eHcal,
0128                    std::vector<unsigned int>* detIds,
0129                    std::vector<double>* hitEnergies);
0130   void fillECALmatrix(const DetId& detId,
0131                       edm::Handle<EcalRecHitCollection>& hitsEB,
0132                       edm::Handle<EcalRecHitCollection>& hitsEE,
0133                       const CaloGeometry* geo,
0134                       const CaloTopology* caloTopology,
0135                       const reco::Track* pTrack,
0136                       int ieta,
0137                       int iphi);
0138   void fillHCALmatrix(const HcalTopology* topology,
0139                       const DetId& detId,
0140                       edm::Handle<HBHERecHitCollection>& hbhe,
0141                       const CaloGeometry* geo,
0142                       const reco::Track* pTrack,
0143                       int ieta,
0144                       int iphi);
0145   void TrackMap(unsigned int trkIndex, std::vector<spr::propagatedTrackDirection>& trkDirs, double dR);
0146   double eThreshold(const DetId& id, double eta) const;
0147 
0148   l1t::L1TGlobalUtil* l1GtUtils_;
0149   HLTConfigProvider hltConfig_;
0150   const std::vector<std::string> trigNames_;
0151   spr::trackSelectionParameters selectionParameter_;
0152   const std::string theTrackQuality_;
0153   const std::string processName_, l1Filter_;
0154   const std::string l2Filter_, l3Filter_;
0155   const double a_coneR_, a_mipR_, pTrackMin_, eEcalMax_;
0156   const double maxRestrictionP_, slopeRestrictionP_;
0157   const double hcalScale_, eIsolate1_, eIsolate2_;
0158   const double pTrackLow_, pTrackHigh_;
0159   const int prescaleLow_, prescaleHigh_;
0160   const int useRaw_, dataType_, mode_;
0161   const bool ignoreTrigger_, useL1Trigger_;
0162   const bool unCorrect_, collapseDepth_;
0163   const double hitEthrEB_, hitEthrEE0_, hitEthrEE1_;
0164   const double hitEthrEE2_, hitEthrEE3_;
0165   const double hitEthrEELo_, hitEthrEEHi_;
0166   const edm::InputTag triggerEvent_, theTriggerResultsLabel_;
0167   const std::string labelGenTrack_, labelRecVtx_, labelEB_;
0168   const std::string labelEE_, labelHBHE_, labelTower_, l1TrigName_;
0169   const int matrixECAL_, matrixHCAL_;
0170   const double mapR_;
0171   const bool get2Ddist_;
0172   const bool usePFThresh_;
0173   const std::string labelBS_, modnam_, prdnam_;
0174   const edm::InputTag algTag_, extTag_;
0175 
0176   const edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
0177   const edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
0178   const edm::EDGetTokenT<reco::GenParticleCollection> tok_parts_;
0179   const edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
0180   const edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
0181   const edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
0182   const edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
0183   const edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
0184   const edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0185   const edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;
0186   const edm::EDGetTokenT<CaloTowerCollection> tok_cala_;
0187   const edm::EDGetTokenT<BXVector<GlobalAlgBlk>> tok_alg_;
0188 
0189   const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_ddrec_;
0190   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_bFieldH_;
0191   const edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> tok_ecalChStatus_;
0192   const edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> tok_sevlv_;
0193   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0194   const edm::ESGetToken<CaloTopology, CaloTopologyRecord> tok_caloTopology_;
0195   const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
0196   const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_resp_;
0197   const edm::ESGetToken<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd> tok_ecalPFRecHitThresholds_;
0198 
0199   unsigned int nRun_, nLow_, nHigh_;
0200   double a_charIsoR_, a_coneR1_, a_coneR2_;
0201   const HcalDDDRecConstants* hdc_;
0202   const EcalPFRecHitThresholds* eThresholds_;
0203   std::vector<double> etabins_, phibins_;
0204   double etadist_, phidist_, etahalfdist_, phihalfdist_;
0205 
0206   TTree *tree_, *tree2_;
0207   unsigned int t_RunNo, t_EventNo;
0208   int t_Run, t_Event, t_DataType, t_ieta, t_iphi;
0209   int t_goodPV, t_nVtx, t_nTrk;
0210   double t_EventWeight, t_p, t_pt, t_phi;
0211   std::vector<double>*t_mapP, *t_mapPt, *t_mapEta, *t_mapPhi;
0212   double t_l1pt, t_l1eta, t_l1phi;
0213   double t_l3pt, t_l3eta, t_l3phi;
0214   double t_mindR1, t_mindR2;
0215   double t_eMipDR, t_hmaxNearP, t_gentrackP, t_gentrackE;
0216   double t_emaxNearP, t_eAnnular, t_hAnnular;
0217   double t_eHcal, t_eHcal10, t_eHcal30, t_rhoh;
0218   bool t_selectTk, t_qltyFlag, t_qltyMissFlag;
0219   bool t_qltyPVFlag, t_TrigPass, t_TrigPassSel;
0220   std::vector<unsigned int>*t_DetIds, *t_DetIds1, *t_DetIds3;
0221   std::vector<double>*t_HitEnergies, *t_HitEnergies1, *t_HitEnergies3;
0222   std::vector<bool>*t_trgbits, *t_hltbits;
0223   bool t_L1Bit;
0224   int t_Tracks, t_TracksProp, t_TracksSaved;
0225   int t_TracksLoose, t_TracksTight, t_allvertex;
0226   std::vector<int>*t_ietaAll, *t_ietaGood, *t_trackType;
0227   std::vector<unsigned int>*t_DetIdEC, *t_DetIdHC;
0228   std::vector<double>*t_HitEnergyEC, *t_HitDistEC;
0229   std::vector<double>*t_HitEnergyHC, *t_HitDistHC;
0230 };
0231 
0232 HcalIsoTrackStudy::HcalIsoTrackStudy(const edm::ParameterSet& iConfig)
0233     : trigNames_(iConfig.getParameter<std::vector<std::string>>("triggers")),
0234       theTrackQuality_(iConfig.getParameter<std::string>("trackQuality")),
0235       processName_(iConfig.getParameter<std::string>("processName")),
0236       l1Filter_(iConfig.getParameter<std::string>("l1Filter")),
0237       l2Filter_(iConfig.getParameter<std::string>("l2Filter")),
0238       l3Filter_(iConfig.getParameter<std::string>("l3Filter")),
0239       a_coneR_(iConfig.getParameter<double>("coneRadius")),
0240       a_mipR_(iConfig.getParameter<double>("coneRadiusMIP")),
0241       pTrackMin_(iConfig.getParameter<double>("minimumTrackP")),
0242       eEcalMax_(iConfig.getParameter<double>("maximumEcalEnergy")),
0243       maxRestrictionP_(iConfig.getParameter<double>("maxTrackP")),
0244       slopeRestrictionP_(iConfig.getParameter<double>("slopeTrackP")),
0245       hcalScale_(iConfig.getUntrackedParameter<double>("hHcalScale", 1.0)),
0246       eIsolate1_(iConfig.getParameter<double>("isolationEnergyTight")),
0247       eIsolate2_(iConfig.getParameter<double>("isolationEnergyLoose")),
0248       pTrackLow_(iConfig.getParameter<double>("momentumLow")),
0249       pTrackHigh_(iConfig.getParameter<double>("momentumHigh")),
0250       prescaleLow_(iConfig.getParameter<int>("prescaleLow")),
0251       prescaleHigh_(iConfig.getParameter<int>("prescaleHigh")),
0252       useRaw_(iConfig.getUntrackedParameter<int>("useRaw", 0)),
0253       dataType_(iConfig.getUntrackedParameter<int>("dataType", 0)),
0254       mode_(iConfig.getUntrackedParameter<int>("outMode", 11)),
0255       ignoreTrigger_(iConfig.getUntrackedParameter<bool>("ignoreTriggers", false)),
0256       useL1Trigger_(iConfig.getUntrackedParameter<bool>("useL1Trigger", false)),
0257       unCorrect_(iConfig.getUntrackedParameter<bool>("unCorrect", false)),
0258       collapseDepth_(iConfig.getUntrackedParameter<bool>("collapseDepth", false)),
0259       hitEthrEB_(iConfig.getParameter<double>("EBHitEnergyThreshold")),
0260       hitEthrEE0_(iConfig.getParameter<double>("EEHitEnergyThreshold0")),
0261       hitEthrEE1_(iConfig.getParameter<double>("EEHitEnergyThreshold1")),
0262       hitEthrEE2_(iConfig.getParameter<double>("EEHitEnergyThreshold2")),
0263       hitEthrEE3_(iConfig.getParameter<double>("EEHitEnergyThreshold3")),
0264       hitEthrEELo_(iConfig.getParameter<double>("EEHitEnergyThresholdLow")),
0265       hitEthrEEHi_(iConfig.getParameter<double>("EEHitEnergyThresholdHigh")),
0266       triggerEvent_(iConfig.getParameter<edm::InputTag>("labelTriggerEvent")),
0267       theTriggerResultsLabel_(iConfig.getParameter<edm::InputTag>("labelTriggerResult")),
0268       labelGenTrack_(iConfig.getParameter<std::string>("labelTrack")),
0269       labelRecVtx_(iConfig.getParameter<std::string>("labelVertex")),
0270       labelEB_(iConfig.getParameter<std::string>("labelEBRecHit")),
0271       labelEE_(iConfig.getParameter<std::string>("labelEERecHit")),
0272       labelHBHE_(iConfig.getParameter<std::string>("labelHBHERecHit")),
0273       labelTower_(iConfig.getParameter<std::string>("labelCaloTower")),
0274       l1TrigName_(iConfig.getUntrackedParameter<std::string>("l1TrigName", "L1_SingleJet60")),
0275       matrixECAL_(iConfig.getUntrackedParameter<int>("matrixECAL", 5)),
0276       matrixHCAL_(iConfig.getUntrackedParameter<int>("matrixHCAL", 3)),
0277       mapR_(iConfig.getUntrackedParameter<double>("mapRadius", 34.98)),
0278       get2Ddist_(iConfig.getUntrackedParameter<bool>("get2Ddist", false)),
0279       usePFThresh_(iConfig.getParameter<bool>("usePFThreshold")),
0280       labelBS_(iConfig.getParameter<std::string>("labelBeamSpot")),
0281       modnam_(iConfig.getUntrackedParameter<std::string>("moduleName", "")),
0282       prdnam_(iConfig.getUntrackedParameter<std::string>("producerName", "")),
0283       algTag_(iConfig.getParameter<edm::InputTag>("algInputTag")),
0284       extTag_(iConfig.getParameter<edm::InputTag>("extInputTag")),
0285       tok_trigEvt_(consumes<trigger::TriggerEvent>(triggerEvent_)),
0286       tok_trigRes_(consumes<edm::TriggerResults>(theTriggerResultsLabel_)),
0287       tok_parts_(consumes<reco::GenParticleCollection>(edm::InputTag("genParticles"))),
0288       tok_bs_(consumes<reco::BeamSpot>(labelBS_)),
0289       tok_genTrack_(consumes<reco::TrackCollection>(labelGenTrack_)),
0290       tok_recVtx_((modnam_.empty()) ? consumes<reco::VertexCollection>(labelRecVtx_)
0291                                     : consumes<reco::VertexCollection>(edm::InputTag(modnam_, labelRecVtx_, prdnam_))),
0292       tok_EB_((modnam_.empty()) ? consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", labelEB_))
0293                                 : consumes<EcalRecHitCollection>(edm::InputTag(modnam_, labelEB_, prdnam_))),
0294       tok_EE_((modnam_.empty()) ? consumes<EcalRecHitCollection>(edm::InputTag("ecalRecHit", labelEE_))
0295                                 : consumes<EcalRecHitCollection>(edm::InputTag(modnam_, labelEE_, prdnam_))),
0296       tok_hbhe_((modnam_.empty()) ? consumes<HBHERecHitCollection>(labelHBHE_)
0297                                   : consumes<HBHERecHitCollection>(edm::InputTag(modnam_, labelHBHE_, prdnam_))),
0298       tok_ew_(consumes<GenEventInfoProduct>(edm::InputTag("generator"))),
0299       tok_cala_(consumes<CaloTowerCollection>(labelTower_)),
0300       tok_alg_(consumes<BXVector<GlobalAlgBlk>>(algTag_)),
0301       tok_ddrec_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()),
0302       tok_bFieldH_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0303       tok_ecalChStatus_(esConsumes<EcalChannelStatus, EcalChannelStatusRcd>()),
0304       tok_sevlv_(esConsumes<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd>()),
0305       tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0306       tok_caloTopology_(esConsumes<CaloTopology, CaloTopologyRecord>()),
0307       tok_htopo_(esConsumes<HcalTopology, HcalRecNumberingRecord>()),
0308       tok_resp_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd>()),
0309       tok_ecalPFRecHitThresholds_(esConsumes<EcalPFRecHitThresholds, EcalPFRecHitThresholdsRcd>()),
0310       nRun_(0),
0311       nLow_(0),
0312       nHigh_(0),
0313       hdc_(nullptr) {
0314   usesResource(TFileService::kSharedResource);
0315 
0316   //now do whatever initialization is needed
0317   const double isolationRadius(28.9), innerR(10.0), outerR(30.0);
0318   reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality_);
0319   selectionParameter_.minPt = iConfig.getParameter<double>("minTrackPt");
0320   ;
0321   selectionParameter_.minQuality = trackQuality_;
0322   selectionParameter_.maxDxyPV = iConfig.getParameter<double>("maxDxyPV");
0323   selectionParameter_.maxDzPV = iConfig.getParameter<double>("maxDzPV");
0324   selectionParameter_.maxChi2 = iConfig.getParameter<double>("maxChi2");
0325   selectionParameter_.maxDpOverP = iConfig.getParameter<double>("maxDpOverP");
0326   selectionParameter_.minOuterHit = iConfig.getParameter<int>("minOuterHit");
0327   selectionParameter_.minLayerCrossed = iConfig.getParameter<int>("minLayerCrossed");
0328   selectionParameter_.maxInMiss = iConfig.getParameter<int>("maxInMiss");
0329   selectionParameter_.maxOutMiss = iConfig.getParameter<int>("maxOutMiss");
0330   a_charIsoR_ = a_coneR_ + isolationRadius;
0331   a_coneR1_ = a_coneR_ + innerR;
0332   a_coneR2_ = a_coneR_ + outerR;
0333   // Different isolation cuts are described in DN-2016/029
0334   // Tight cut uses 2 GeV; Loose cut uses 10 GeV
0335   // Eta dependent cut uses (maxRestrictionP_ * exp(|ieta|*log(2.5)/18))
0336   // with the factor for exponential slopeRestrictionP_ = log(2.5)/18
0337   // maxRestrictionP_ = 8 GeV as came from a study
0338   l1GtUtils_ = new l1t::L1TGlobalUtil(iConfig, consumesCollector(), *this, algTag_, extTag_);
0339   // define tokens for access
0340 
0341   if (modnam_.empty()) {
0342     edm::LogVerbatim("HcalIsoTrack") << "Labels used " << triggerEvent_ << " " << theTriggerResultsLabel_ << " "
0343                                      << labelBS_ << " " << labelRecVtx_ << " " << labelGenTrack_ << " "
0344                                      << edm::InputTag("ecalRecHit", labelEB_) << " "
0345                                      << edm::InputTag("ecalRecHit", labelEE_) << " " << labelHBHE_ << " "
0346                                      << labelTower_;
0347   } else {
0348     edm::LogVerbatim("HcalIsoTrack") << "Labels used " << triggerEvent_ << " " << theTriggerResultsLabel_ << " "
0349                                      << labelBS_ << " " << edm::InputTag(modnam_, labelRecVtx_, prdnam_) << " "
0350                                      << labelGenTrack_ << " " << edm::InputTag(modnam_, labelEB_, prdnam_) << " "
0351                                      << edm::InputTag(modnam_, labelEE_, prdnam_) << " "
0352                                      << edm::InputTag(modnam_, labelHBHE_, prdnam_) << " " << labelTower_;
0353   }
0354 
0355   edm::LogVerbatim("HcalIsoTrack") << "Parameters read from config file \n"
0356                                    << "\t minPt " << selectionParameter_.minPt << "\t theTrackQuality "
0357                                    << theTrackQuality_ << "\t minQuality " << selectionParameter_.minQuality
0358                                    << "\t maxDxyPV " << selectionParameter_.maxDxyPV << "\t maxDzPV "
0359                                    << selectionParameter_.maxDzPV << "\t maxChi2 " << selectionParameter_.maxChi2
0360                                    << "\t maxDpOverP " << selectionParameter_.maxDpOverP << "\t minOuterHit "
0361                                    << selectionParameter_.minOuterHit << "\t minLayerCrossed "
0362                                    << selectionParameter_.minLayerCrossed << "\t maxInMiss "
0363                                    << selectionParameter_.maxInMiss << "\t maxOutMiss "
0364                                    << selectionParameter_.maxOutMiss << "\t a_coneR " << a_coneR_ << ":" << a_coneR1_
0365                                    << ":" << a_coneR2_ << "\t a_charIsoR " << a_charIsoR_ << "\t a_mipR " << a_mipR_
0366                                    << "\n pTrackMin_ " << pTrackMin_ << "\t eEcalMax_ " << eEcalMax_
0367                                    << "\t maxRestrictionP_ " << maxRestrictionP_ << "\t slopeRestrictionP_ "
0368                                    << slopeRestrictionP_ << "\t eIsolateStrong_ " << eIsolate1_ << "\t eIsolateSoft_ "
0369                                    << eIsolate2_ << "\t hcalScale_ " << hcalScale_ << "\n\t momentumLow_ " << pTrackLow_
0370                                    << "\t prescaleLow_ " << prescaleLow_ << "\t momentumHigh_ " << pTrackHigh_
0371                                    << "\t prescaleHigh_ " << prescaleHigh_ << "\n\t useRaw_ " << useRaw_
0372                                    << "\t ignoreTrigger_ " << ignoreTrigger_ << "\n\t useL1Trigegr_ " << useL1Trigger_
0373                                    << "\t dataType_      " << dataType_ << "\t mode_          " << mode_
0374                                    << "\t unCorrect_     " << unCorrect_ << "\t collapseDepth_ " << collapseDepth_
0375                                    << "\t L1TrigName_    " << l1TrigName_;
0376   edm::LogVerbatim("HcalIsoTrack") << "Process " << processName_ << " L1Filter:" << l1Filter_
0377                                    << " L2Filter:" << l2Filter_ << " L3Filter:" << l3Filter_;
0378   for (unsigned int k = 0; k < trigNames_.size(); ++k) {
0379     edm::LogVerbatim("HcalIsoTrack") << "Trigger[" << k << "] " << trigNames_[k];
0380   }
0381 
0382   for (int i = 0; i < 10; i++)
0383     phibins_.push_back(-M_PI + 0.1 * (2 * i + 1) * M_PI);
0384   for (int i = 0; i < 8; ++i)
0385     etabins_.push_back(-2.1 + 0.6 * i);
0386   etadist_ = etabins_[1] - etabins_[0];
0387   phidist_ = phibins_[1] - phibins_[0];
0388   etahalfdist_ = 0.5 * etadist_;
0389   phihalfdist_ = 0.5 * phidist_;
0390   edm::LogVerbatim("HcalIsoTrack") << "EtaDist " << etadist_ << " " << etahalfdist_ << " PhiDist " << phidist_ << " "
0391                                    << phihalfdist_;
0392   unsigned int k1(0), k2(0);
0393   for (auto phi : phibins_) {
0394     edm::LogVerbatim("HcalIsoTrack") << "phibin_[" << k1 << "] " << phi;
0395     ++k1;
0396   }
0397   for (auto eta : etabins_) {
0398     edm::LogVerbatim("HcalIsoTrack") << "etabin_[" << k2 << "] " << eta;
0399     ++k2;
0400   }
0401 }
0402 
0403 void HcalIsoTrackStudy::analyze(edm::Event const& iEvent, edm::EventSetup const& iSetup) {
0404   t_Run = iEvent.id().run();
0405   t_Event = iEvent.id().event();
0406   t_DataType = dataType_;
0407 #ifdef EDM_ML_DEBUG
0408   edm::LogVerbatim("HcalIsoTrack") << "Run " << t_Run << " Event " << t_Event << " type " << t_DataType
0409                                    << " Luminosity " << iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing();
0410 #endif
0411   //Get magnetic field and ECAL channel status
0412   const MagneticField* bField = &iSetup.getData(tok_bFieldH_);
0413   const EcalChannelStatus* theEcalChStatus = &iSetup.getData(tok_ecalChStatus_);
0414   const EcalSeverityLevelAlgo* theEcalSevlv = &iSetup.getData(tok_sevlv_);
0415 
0416   // get handles to calogeometry and calotopology
0417   const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0418   const CaloTopology* caloTopology = &iSetup.getData(tok_caloTopology_);
0419   const HcalTopology* theHBHETopology = &iSetup.getData(tok_htopo_);
0420   const HcalRespCorrs* respCorrs = &iSetup.getData(tok_resp_);
0421 
0422   // get ECAL thresholds
0423   eThresholds_ = &iSetup.getData(tok_ecalPFRecHitThresholds_);
0424 
0425   //=== genParticle information
0426   const edm::Handle<reco::GenParticleCollection> genParticles = iEvent.getHandle(tok_parts_);
0427 
0428   bool okC(true);
0429   //Get track collection
0430   edm::Handle<reco::TrackCollection> trkCollection = iEvent.getHandle(tok_genTrack_);
0431   if (!trkCollection.isValid()) {
0432     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
0433     okC = false;
0434   }
0435 
0436   //event weight for FLAT sample
0437   t_EventWeight = 1.0;
0438   const edm::Handle<GenEventInfoProduct> genEventInfo = iEvent.getHandle(tok_ew_);
0439   if (genEventInfo.isValid())
0440     t_EventWeight = genEventInfo->weight();
0441 
0442   //Define the best vertex and the beamspot
0443   const edm::Handle<reco::VertexCollection> recVtxs = iEvent.getHandle(tok_recVtx_);
0444   const edm::Handle<reco::BeamSpot> beamSpotH = iEvent.getHandle(tok_bs_);
0445   math::XYZPoint leadPV(0, 0, 0);
0446   t_goodPV = t_nVtx = 0;
0447   if (recVtxs.isValid() && !(recVtxs->empty())) {
0448     t_nVtx = recVtxs->size();
0449     for (unsigned int k = 0; k < recVtxs->size(); ++k) {
0450       if (!((*recVtxs)[k].isFake()) && ((*recVtxs)[k].ndof() > 4)) {
0451         if (t_goodPV == 0)
0452           leadPV = math::XYZPoint((*recVtxs)[k].x(), (*recVtxs)[k].y(), (*recVtxs)[k].z());
0453         t_goodPV++;
0454       }
0455     }
0456   }
0457   if (t_goodPV == 0 && beamSpotH.isValid()) {
0458     leadPV = beamSpotH->position();
0459   }
0460   t_allvertex = t_goodPV;
0461 #ifdef EDM_ML_DEBUG
0462   edm::LogVerbatim("HcalIsoTrack") << "Primary Vertex " << leadPV << " out of " << t_goodPV << " vertex";
0463   if (beamSpotH.isValid()) {
0464     edm::LogVerbatim("HcalIsoTrack") << " Beam Spot " << beamSpotH->position();
0465   }
0466 #endif
0467   // RecHits
0468   edm::Handle<EcalRecHitCollection> barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
0469   if (!barrelRecHitsHandle.isValid()) {
0470     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEB_;
0471     okC = false;
0472   }
0473   edm::Handle<EcalRecHitCollection> endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
0474   if (!endcapRecHitsHandle.isValid()) {
0475     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEE_;
0476     okC = false;
0477   }
0478   edm::Handle<HBHERecHitCollection> hbhe = iEvent.getHandle(tok_hbhe_);
0479   if (!hbhe.isValid()) {
0480     edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelHBHE_;
0481     okC = false;
0482   }
0483   edm::Handle<CaloTowerCollection> caloTower = iEvent.getHandle(tok_cala_);
0484 
0485   //Propagate tracks to calorimeter surface)
0486   std::vector<spr::propagatedTrackDirection> trkCaloDirections;
0487   spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDirections, false);
0488   std::vector<spr::propagatedTrackID> trkCaloDets;
0489   spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDets, false);
0490   std::vector<math::XYZTLorentzVector> vecL1, vecL3;
0491   t_RunNo = iEvent.id().run();
0492   t_EventNo = iEvent.id().event();
0493   t_Tracks = trkCollection->size();
0494   t_TracksProp = trkCaloDirections.size();
0495   t_ietaAll->clear();
0496   t_ietaGood->clear();
0497   t_trackType->clear();
0498   t_trgbits->clear();
0499   t_hltbits->clear();
0500 #ifdef EDM_ML_DEBUG
0501   edm::LogVerbatim("HcalIsoTrack") << "# of propagated tracks " << t_TracksProp << " out of " << t_Tracks
0502                                    << " with Trigger " << ignoreTrigger_;
0503 #endif
0504 
0505   //Trigger
0506   t_trgbits->assign(trigNames_.size(), false);
0507   t_hltbits->assign(trigNames_.size(), false);
0508   t_TracksSaved = t_TracksLoose = t_TracksTight = 0;
0509   t_L1Bit = true;
0510   t_TrigPass = false;
0511 
0512   //L1
0513   l1GtUtils_->retrieveL1(iEvent, iSetup, tok_alg_);
0514   const std::vector<std::pair<std::string, bool>>& finalDecisions = l1GtUtils_->decisionsFinal();
0515   for (const auto& decision : finalDecisions) {
0516     if (decision.first.find(l1TrigName_) != std::string::npos) {
0517       t_L1Bit = decision.second;
0518       break;
0519     }
0520   }
0521 #ifdef EDM_ML_DEBUG
0522   edm::LogVerbatim("HcalIsoTrack") << "Trigger Information for " << l1TrigName_ << " is " << t_L1Bit
0523                                    << " from a list of " << finalDecisions.size() << " decisions";
0524 #endif
0525 
0526   //HLT
0527   const edm::Handle<edm::TriggerResults> triggerResults = iEvent.getHandle(tok_trigRes_);
0528   if (triggerResults.isValid()) {
0529     const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0530     const std::vector<std::string>& names = triggerNames.triggerNames();
0531     if (!trigNames_.empty()) {
0532       for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0533         int hlt = triggerResults->accept(iHLT);
0534         for (unsigned int i = 0; i < trigNames_.size(); ++i) {
0535           if (names[iHLT].find(trigNames_[i]) != std::string::npos) {
0536             t_trgbits->at(i) = (hlt > 0);
0537             t_hltbits->at(i) = (hlt > 0);
0538             if (hlt > 0)
0539               t_TrigPass = true;
0540 #ifdef EDM_ML_DEBUG
0541             edm::LogVerbatim("HcalIsoTrack")
0542                 << "This trigger " << names[iHLT] << " Flag " << hlt << ":" << t_trgbits->at(i);
0543 #endif
0544           }
0545         }
0546       }
0547     }
0548   }
0549 #ifdef EDM_ML_DEBUG
0550   edm::LogVerbatim("HcalIsoTrack") << "HLT Information shows " << t_TrigPass << ":" << trigNames_.empty() << ":" << okC;
0551 #endif
0552 
0553   std::array<int, 3> ntksave{{0, 0, 0}};
0554   if (ignoreTrigger_ || useL1Trigger_) {
0555     t_l1pt = t_l1eta = t_l1phi = 0;
0556     t_l3pt = t_l3eta = t_l3phi = 0;
0557     if (ignoreTrigger_ || t_L1Bit)
0558       ntksave = fillTree(vecL1,
0559                          vecL3,
0560                          leadPV,
0561                          trkCaloDirections,
0562                          trkCaloDets,
0563                          geo,
0564                          caloTopology,
0565                          theHBHETopology,
0566                          theEcalChStatus,
0567                          theEcalSevlv,
0568                          barrelRecHitsHandle,
0569                          endcapRecHitsHandle,
0570                          hbhe,
0571                          caloTower,
0572                          genParticles,
0573                          respCorrs);
0574     t_TracksSaved = ntksave[0];
0575     t_TracksLoose = ntksave[1];
0576     t_TracksTight = ntksave[2];
0577   } else {
0578     trigger::TriggerEvent triggerEvent;
0579     const edm::Handle<trigger::TriggerEvent> triggerEventHandle = iEvent.getHandle(tok_trigEvt_);
0580     if (!triggerEventHandle.isValid()) {
0581       edm::LogWarning("HcalIsoTrack") << "Error! Can't get the product " << triggerEvent_.label();
0582     } else if (okC) {
0583       triggerEvent = *(triggerEventHandle.product());
0584       const trigger::TriggerObjectCollection& TOC(triggerEvent.getObjects());
0585       bool done(false);
0586       if (triggerResults.isValid()) {
0587         std::vector<std::string> modules;
0588         const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
0589         const std::vector<std::string>& names = triggerNames.triggerNames();
0590         for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
0591           bool ok = (t_TrigPass) || (trigNames_.empty());
0592           if (ok) {
0593             unsigned int triggerindx = hltConfig_.triggerIndex(names[iHLT]);
0594             const std::vector<std::string>& moduleLabels(hltConfig_.moduleLabels(triggerindx));
0595             std::vector<math::XYZTLorentzVector> vecL2;
0596             vecL1.clear();
0597             vecL3.clear();
0598             //loop over all trigger filters in event (i.e. filters passed)
0599             for (unsigned int ifilter = 0; ifilter < triggerEvent.sizeFilters(); ++ifilter) {
0600               std::vector<int> Keys;
0601               std::string label = triggerEvent.filterTag(ifilter).label();
0602               //loop over keys to objects passing this filter
0603               for (unsigned int imodule = 0; imodule < moduleLabels.size(); imodule++) {
0604                 if (label.find(moduleLabels[imodule]) != std::string::npos) {
0605 #ifdef EDM_ML_DEBUG
0606                   edm::LogVerbatim("HcalIsoTrack") << "FilterName " << label;
0607 #endif
0608                   for (unsigned int ifiltrKey = 0; ifiltrKey < triggerEvent.filterKeys(ifilter).size(); ++ifiltrKey) {
0609                     Keys.push_back(triggerEvent.filterKeys(ifilter)[ifiltrKey]);
0610                     const trigger::TriggerObject& TO(TOC[Keys[ifiltrKey]]);
0611                     math::XYZTLorentzVector v4(TO.px(), TO.py(), TO.pz(), TO.energy());
0612                     if (label.find(l2Filter_) != std::string::npos) {
0613                       vecL2.push_back(v4);
0614                     } else if (label.find(l3Filter_) != std::string::npos) {
0615                       vecL3.push_back(v4);
0616                     } else if ((label.find(l1Filter_) != std::string::npos) || (l1Filter_.empty())) {
0617                       vecL1.push_back(v4);
0618                     }
0619 #ifdef EDM_ML_DEBUG
0620                     edm::LogVerbatim("HcalIsoTrack")
0621                         << "key " << ifiltrKey << " : pt " << TO.pt() << " eta " << TO.eta() << " phi " << TO.phi()
0622                         << " mass " << TO.mass() << " Id " << TO.id();
0623 #endif
0624                   }
0625 #ifdef EDM_ML_DEBUG
0626                   edm::LogVerbatim("HcalIsoTrack")
0627                       << "sizes " << vecL1.size() << ":" << vecL2.size() << ":" << vecL3.size();
0628 #endif
0629                 }
0630               }
0631             }
0632             //// deta, dphi and dR for leading L1 object with L2 objects
0633             math::XYZTLorentzVector mindRvec1;
0634             double mindR1(999);
0635             for (unsigned int i = 0; i < vecL2.size(); i++) {
0636               double dr = dR(vecL1[0], vecL2[i]);
0637 #ifdef EDM_ML_DEBUG
0638               edm::LogVerbatim("HcalIsoTrack") << "lvl2[" << i << "] dR " << dr;
0639 #endif
0640               if (dr < mindR1) {
0641                 mindR1 = dr;
0642                 mindRvec1 = vecL2[i];
0643               }
0644             }
0645 #ifdef EDM_ML_DEBUG
0646             edm::LogVerbatim("HcalIsoTrack") << "L2 object closest to L1 " << mindRvec1 << " at Dr " << mindR1;
0647 #endif
0648 
0649             if (!vecL1.empty()) {
0650               t_l1pt = vecL1[0].pt();
0651               t_l1eta = vecL1[0].eta();
0652               t_l1phi = vecL1[0].phi();
0653             } else {
0654               t_l1pt = t_l1eta = t_l1phi = 0;
0655             }
0656             if (!vecL3.empty()) {
0657               t_l3pt = vecL3[0].pt();
0658               t_l3eta = vecL3[0].eta();
0659               t_l3phi = vecL3[0].phi();
0660             } else {
0661               t_l3pt = t_l3eta = t_l3phi = 0;
0662             }
0663             // Now fill in the tree for each selected track
0664             if (!done) {
0665               ntksave = fillTree(vecL1,
0666                                  vecL3,
0667                                  leadPV,
0668                                  trkCaloDirections,
0669                                  trkCaloDets,
0670                                  geo,
0671                                  caloTopology,
0672                                  theHBHETopology,
0673                                  theEcalChStatus,
0674                                  theEcalSevlv,
0675                                  barrelRecHitsHandle,
0676                                  endcapRecHitsHandle,
0677                                  hbhe,
0678                                  caloTower,
0679                                  genParticles,
0680                                  respCorrs);
0681               t_TracksSaved += ntksave[0];
0682               t_TracksLoose += ntksave[1];
0683               t_TracksTight += ntksave[2];
0684               done = true;
0685             }
0686           }
0687         }
0688       }
0689     }
0690   }
0691 #ifdef EDM_ML_DEBUG
0692   edm::LogVerbatim("HcalIsoTrack") << "Final results on selected tracks " << t_TracksSaved << ":" << t_TracksLoose
0693                                    << ":" << t_TracksTight;
0694 #endif
0695   t_TrigPassSel = (t_TracksSaved > 0);
0696   tree2_->Fill();
0697 }
0698 
0699 void HcalIsoTrackStudy::beginJob() {
0700   edm::Service<TFileService> fs;
0701   tree_ = fs->make<TTree>("CalibTree", "CalibTree");
0702 
0703   tree_->Branch("t_Run", &t_Run, "t_Run/I");
0704   tree_->Branch("t_Event", &t_Event, "t_Event/I");
0705   tree_->Branch("t_DataType", &t_DataType, "t_DataType/I");
0706   tree_->Branch("t_ieta", &t_ieta, "t_ieta/I");
0707   tree_->Branch("t_iphi", &t_iphi, "t_iphi/I");
0708   tree_->Branch("t_EventWeight", &t_EventWeight, "t_EventWeight/D");
0709   tree_->Branch("t_nVtx", &t_nVtx, "t_nVtx/I");
0710   tree_->Branch("t_nTrk", &t_nTrk, "t_nTrk/I");
0711   tree_->Branch("t_goodPV", &t_goodPV, "t_goodPV/I");
0712   if (((mode_ / 10) % 10) == 1) {
0713     tree_->Branch("t_l1pt", &t_l1pt, "t_l1pt/D");
0714     tree_->Branch("t_l1eta", &t_l1eta, "t_l1eta/D");
0715     tree_->Branch("t_l1phi", &t_l1phi, "t_l1phi/D");
0716     tree_->Branch("t_l3pt", &t_l3pt, "t_l3pt/D");
0717     tree_->Branch("t_l3eta", &t_l3eta, "t_l3eta/D");
0718     tree_->Branch("t_l3phi", &t_l3phi, "t_l3phi/D");
0719   }
0720   tree_->Branch("t_p", &t_p, "t_p/D");
0721   tree_->Branch("t_pt", &t_pt, "t_pt/D");
0722   tree_->Branch("t_phi", &t_phi, "t_phi/D");
0723 
0724   t_mapP = new std::vector<double>();
0725   t_mapPt = new std::vector<double>();
0726   t_mapEta = new std::vector<double>();
0727   t_mapPhi = new std::vector<double>();
0728   tree_->Branch("t_mapP", "std::vector<double>", &t_mapP);
0729   tree_->Branch("t_mapPt", "std::vector<double>", &t_mapPt);
0730   tree_->Branch("t_mapEta", "std::vector<double>", &t_mapEta);
0731   tree_->Branch("t_mapPhi", "std::vector<double>", &t_mapPhi);
0732 
0733   tree_->Branch("t_mindR1", &t_mindR1, "t_mindR1/D");
0734   tree_->Branch("t_mindR2", &t_mindR2, "t_mindR2/D");
0735   tree_->Branch("t_eMipDR", &t_eMipDR, "t_eMipDR/D");
0736   tree_->Branch("t_eHcal", &t_eHcal, "t_eHcal/D");
0737   tree_->Branch("t_eHcal10", &t_eHcal10, "t_eHcal10/D");
0738   tree_->Branch("t_eHcal30", &t_eHcal30, "t_eHcal30/D");
0739   tree_->Branch("t_hmaxNearP", &t_hmaxNearP, "t_hmaxNearP/D");
0740   tree_->Branch("t_emaxNearP", &t_emaxNearP, "t_emaxNearP/D");
0741   tree_->Branch("t_eAnnular", &t_eAnnular, "t_eAnnular/D");
0742   tree_->Branch("t_hAnnular", &t_hAnnular, "t_hAnnular/D");
0743   tree_->Branch("t_rhoh", &t_rhoh, "t_rhoh/D");
0744   tree_->Branch("t_selectTk", &t_selectTk, "t_selectTk/O");
0745   tree_->Branch("t_qltyFlag", &t_qltyFlag, "t_qltyFlag/O");
0746   tree_->Branch("t_qltyMissFlag", &t_qltyMissFlag, "t_qltyMissFlag/O");
0747   tree_->Branch("t_qltyPVFlag", &t_qltyPVFlag, "t_qltyPVFlag/O");
0748   tree_->Branch("t_gentrackP", &t_gentrackP, "t_gentrackP/D");
0749   tree_->Branch("t_gentrackE", &t_gentrackE, "t_gentrackE/D");
0750 
0751   t_DetIds = new std::vector<unsigned int>();
0752   t_DetIds1 = new std::vector<unsigned int>();
0753   t_DetIds3 = new std::vector<unsigned int>();
0754   t_DetIdEC = new std::vector<unsigned int>();
0755   t_DetIdHC = new std::vector<unsigned int>();
0756   t_HitEnergies = new std::vector<double>();
0757   t_HitEnergies1 = new std::vector<double>();
0758   t_HitEnergies3 = new std::vector<double>();
0759   t_HitEnergyEC = new std::vector<double>();
0760   t_HitDistEC = new std::vector<double>();
0761   t_HitEnergyHC = new std::vector<double>();
0762   t_HitDistHC = new std::vector<double>();
0763   t_trgbits = new std::vector<bool>();
0764   tree_->Branch("t_DetIds", "std::vector<unsigned int>", &t_DetIds);
0765   tree_->Branch("t_HitEnergies", "std::vector<double>", &t_HitEnergies);
0766   if (((mode_ / 10) % 10) == 1) {
0767     tree_->Branch("t_trgbits", "std::vector<bool>", &t_trgbits);
0768   }
0769   if ((mode_ % 10) == 1) {
0770     tree_->Branch("t_DetIds1", "std::vector<unsigned int>", &t_DetIds1);
0771     tree_->Branch("t_HitEnergies1", "std::vector<double>", &t_HitEnergies1);
0772     tree_->Branch("t_DetIds3", "std::vector<unsigned int>", &t_DetIds3);
0773     tree_->Branch("t_HitEnergies3", "std::vector<double>", &t_HitEnergies3);
0774     tree_->Branch("t_DetIdEC", "std::vector<unsigned int>", &t_DetIdEC);
0775     tree_->Branch("t_HitEnergyEC", "std::vector<double>", &t_HitEnergyEC);
0776     tree_->Branch("t_HitDistEC", "std::vector<double>", &t_HitDistEC);
0777     tree_->Branch("t_DetIdHC", "std::vector<unsigned int>", &t_DetIdHC);
0778     tree_->Branch("t_HitEnergyHC", "std::vector<double>", &t_HitEnergyHC);
0779     tree_->Branch("t_HitDistHC", "std::vector<double>", &t_HitDistHC);
0780   }
0781   tree2_ = fs->make<TTree>("EventInfo", "Event Information");
0782 
0783   tree2_->Branch("t_RunNo", &t_RunNo, "t_RunNo/i");
0784   tree2_->Branch("t_EventNo", &t_EventNo, "t_EventNo/i");
0785   tree2_->Branch("t_Tracks", &t_Tracks, "t_Tracks/I");
0786   tree2_->Branch("t_TracksProp", &t_TracksProp, "t_TracksProp/I");
0787   tree2_->Branch("t_TracksSaved", &t_TracksSaved, "t_TracksSaved/I");
0788   tree2_->Branch("t_TracksLoose", &t_TracksLoose, "t_TracksLoose/I");
0789   tree2_->Branch("t_TracksTight", &t_TracksTight, "t_TracksTight/I");
0790   tree2_->Branch("t_TrigPass", &t_TrigPass, "t_TrigPass/O");
0791   tree2_->Branch("t_TrigPassSel", &t_TrigPassSel, "t_TrigPassSel/O");
0792   tree2_->Branch("t_L1Bit", &t_L1Bit, "t_L1Bit/O");
0793   tree2_->Branch("t_allvertex", &t_allvertex, "t_allvertex/I");
0794   t_hltbits = new std::vector<bool>();
0795   t_ietaAll = new std::vector<int>();
0796   t_ietaGood = new std::vector<int>();
0797   t_trackType = new std::vector<int>();
0798   tree2_->Branch("t_ietaAll", "std::vector<int>", &t_ietaAll);
0799   tree2_->Branch("t_ietaGood", "std::vector<int>", &t_ietaGood);
0800   tree2_->Branch("t_trackType", "std::vector<int>", &t_trackType);
0801   tree2_->Branch("t_hltbits", "std::vector<bool>", &t_hltbits);
0802 }
0803 
0804 // ------------ method called when starting to processes a run  ------------
0805 void HcalIsoTrackStudy::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
0806   hdc_ = &iSetup.getData(tok_ddrec_);
0807 
0808   bool changed_(true);
0809   bool flag = hltConfig_.init(iRun, iSetup, processName_, changed_);
0810   edm::LogVerbatim("HcalIsoTrack") << "Run[" << nRun_ << "] " << iRun.run() << " process " << processName_
0811                                    << " init flag " << flag << " change flag " << changed_;
0812   // check if trigger names in (new) config
0813   if (changed_) {
0814     changed_ = false;
0815 #ifdef EDM_ML_DEBUG
0816     edm::LogVerbatim("HcalIsoTrack") << "New trigger menu found !!!";
0817 #endif
0818     const unsigned int n(hltConfig_.size());
0819     for (unsigned itrig = 0; itrig < trigNames_.size(); itrig++) {
0820       unsigned int triggerindx = hltConfig_.triggerIndex(trigNames_[itrig]);
0821       if (triggerindx >= n) {
0822         edm::LogWarning("HcalIsoTrack") << trigNames_[itrig] << " " << triggerindx << " does not exist in "
0823                                         << "the current menu";
0824 #ifdef EDM_ML_DEBUG
0825       } else {
0826         edm::LogVerbatim("HcalIsoTrack") << trigNames_[itrig] << " " << triggerindx << " exists";
0827 #endif
0828       }
0829     }
0830   }
0831 }
0832 
0833 // ------------ method called when ending the processing of a run  ------------
0834 void HcalIsoTrackStudy::endRun(edm::Run const& iRun, edm::EventSetup const&) {
0835   nRun_++;
0836   edm::LogVerbatim("HcalIsoTrack") << "endRun[" << nRun_ << "] " << iRun.run();
0837 }
0838 
0839 void HcalIsoTrackStudy::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0840   edm::ParameterSetDescription desc;
0841   std::vector<std::string> trig = {"HLT_PFJet40",
0842                                    "HLT_PFJet60",
0843                                    "HLT_PFJet80",
0844                                    "HLT_PFJet140",
0845                                    "HLT_PFJet200",
0846                                    "HLT_PFJet260",
0847                                    "HLT_PFJet320",
0848                                    "HLT_PFJet400",
0849                                    "HLT_PFJet450",
0850                                    "HLT_PFJet500"};
0851   desc.add<std::vector<std::string>>("triggers", trig);
0852   desc.add<std::string>("processName", "HLT");
0853   desc.add<std::string>("l1Filter", "");
0854   desc.add<std::string>("l2Filter", "L2Filter");
0855   desc.add<std::string>("l3Filter", "Filter");
0856   // following 10 parameters are parameters to select good tracks
0857   desc.add<std::string>("trackQuality", "highPurity");
0858   desc.add<double>("minTrackPt", 1.0);
0859   desc.add<double>("maxDxyPV", 0.02);
0860   desc.add<double>("maxDzPV", 0.02);
0861   desc.add<double>("maxChi2", 5.0);
0862   desc.add<double>("maxDpOverP", 0.1);
0863   desc.add<int>("minOuterHit", 4);
0864   desc.add<int>("minLayerCrossed", 8);
0865   desc.add<int>("maxInMiss", 0);
0866   desc.add<int>("maxOutMiss", 0);
0867   // Minimum momentum of selected isolated track and signal zone
0868   desc.add<double>("minimumTrackP", 10.0);
0869   desc.add<double>("coneRadius", 34.98);
0870   // signal zone in ECAL and MIP energy cutoff
0871   desc.add<double>("coneRadiusMIP", 14.0);
0872   desc.add<double>("maximumEcalEnergy", 2.0);
0873   // following 4 parameters are for isolation cuts and described in the code
0874   desc.add<double>("maxTrackP", 8.0);
0875   desc.add<double>("slopeTrackP", 0.05090504066);
0876   desc.add<double>("isolationEnergyTight", 2.0);
0877   desc.add<double>("isolationEnergyLoose", 10.0);
0878   // energy thershold for ECAL (from Egamma group)
0879   desc.add<double>("EBHitEnergyThreshold", 0.10);
0880   desc.add<double>("EEHitEnergyThreshold0", -41.0664);
0881   desc.add<double>("EEHitEnergyThreshold1", 68.7950);
0882   desc.add<double>("EEHitEnergyThreshold2", -38.1483);
0883   desc.add<double>("EEHitEnergyThreshold3", 7.04303);
0884   desc.add<double>("EEHitEnergyThresholdLow", 0.08);
0885   desc.add<double>("EEHitEnergyThresholdHigh", 0.30);
0886   // prescale factors
0887   desc.add<double>("momentumLow", 40.0);
0888   desc.add<double>("momentumHigh", 60.0);
0889   desc.add<int>("prescaleLow", 1);
0890   desc.add<int>("prescaleHigh", 1);
0891   // various labels for collections used in the code
0892   desc.add<edm::InputTag>("labelTriggerEvent", edm::InputTag("hltTriggerSummaryAOD", "", "HLT"));
0893   desc.add<edm::InputTag>("labelTriggerResult", edm::InputTag("TriggerResults", "", "HLT"));
0894   desc.add<std::string>("labelTrack", "generalTracks");
0895   desc.add<std::string>("labelVertex", "offlinePrimaryVertices");
0896   desc.add<std::string>("labelEBRecHit", "EcalRecHitsEB");
0897   desc.add<std::string>("labelEERecHit", "EcalRecHitsEE");
0898   desc.add<std::string>("labelHBHERecHit", "hbhereco");
0899   desc.add<std::string>("labelBeamSpot", "offlineBeamSpot");
0900   desc.add<std::string>("labelCaloTower", "towerMaker");
0901   desc.add<edm::InputTag>("algInputTag", edm::InputTag("gtStage2Digis"));
0902   desc.add<edm::InputTag>("extInputTag", edm::InputTag("gtStage2Digis"));
0903   desc.addUntracked<std::string>("moduleName", "");
0904   desc.addUntracked<std::string>("producerName", "");
0905   //  Various flags used for selecting tracks, choice of energy Method2/0
0906   //  Data type 0/1 for single jet trigger or others
0907   desc.addUntracked<int>("useRaw", 0);
0908   desc.addUntracked<bool>("ignoreTriggers", false);
0909   desc.addUntracked<bool>("useL1Trigger", false);
0910   desc.addUntracked<double>("hcalScale", 1.0);
0911   desc.addUntracked<int>("dataType", 0);
0912   desc.addUntracked<int>("outMode", 11);
0913   desc.addUntracked<bool>("unCorrect", false);
0914   desc.addUntracked<bool>("collapseDepth", false);
0915   desc.addUntracked<std::string>("l1TrigName", "L1_SingleJet60");
0916   desc.addUntracked<int>("matrixECAL", 5);
0917   desc.addUntracked<int>("matrixHCAL", 3);
0918   desc.addUntracked<double>("mapRadius", 34.98);
0919   desc.addUntracked<bool>("get2Ddist", false);
0920   desc.add<bool>("usePFThreshold", true);
0921   descriptions.add("hcalIsoTrackStudy", desc);
0922 }
0923 
0924 std::array<int, 3> HcalIsoTrackStudy::fillTree(std::vector<math::XYZTLorentzVector>& vecL1,
0925                                                std::vector<math::XYZTLorentzVector>& vecL3,
0926                                                math::XYZPoint& leadPV,
0927                                                std::vector<spr::propagatedTrackDirection>& trkCaloDirections,
0928                                                std::vector<spr::propagatedTrackID>& trkCaloDets,
0929                                                const CaloGeometry* geo,
0930                                                const CaloTopology* caloTopology,
0931                                                const HcalTopology* theHBHETopology,
0932                                                const EcalChannelStatus* theEcalChStatus,
0933                                                const EcalSeverityLevelAlgo* theEcalSevlv,
0934                                                edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
0935                                                edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
0936                                                edm::Handle<HBHERecHitCollection>& hbhe,
0937                                                edm::Handle<CaloTowerCollection>& tower,
0938                                                const edm::Handle<reco::GenParticleCollection>& genParticles,
0939                                                const HcalRespCorrs* respCorrs) {
0940   int nSave(0), nLoose(0), nTight(0);
0941   //Loop over tracks
0942   std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
0943   unsigned int nTracks(0), nselTracks(0);
0944   t_nTrk = trkCaloDirections.size();
0945   t_rhoh = (tower.isValid()) ? rhoh(tower) : 0;
0946   for (trkDetItr = trkCaloDirections.begin(), nTracks = 0; trkDetItr != trkCaloDirections.end();
0947        trkDetItr++, nTracks++) {
0948     const reco::Track* pTrack = &(*(trkDetItr->trkItr));
0949     math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
0950 #ifdef EDM_ML_DEBUG
0951     edm::LogVerbatim("HcalIsoTrack") << "This track : " << nTracks << " (pt|eta|phi|p) :" << pTrack->pt() << "|"
0952                                      << pTrack->eta() << "|" << pTrack->phi() << "|" << pTrack->p();
0953 #endif
0954     t_mindR2 = 999;
0955     for (unsigned int k = 0; k < vecL3.size(); ++k) {
0956       double dr = dR(vecL3[k], v4);
0957       if (dr < t_mindR2) {
0958         t_mindR2 = dr;
0959       }
0960     }
0961     t_mindR1 = (!vecL1.empty()) ? dR(vecL1[0], v4) : 999;
0962 #ifdef EDM_ML_DEBUG
0963     edm::LogVerbatim("HcalIsoTrack") << "Closest L3 object at dr :" << t_mindR2 << " and from L1 " << t_mindR1;
0964 #endif
0965     t_ieta = t_iphi = 0;
0966     if (trkDetItr->okHCAL) {
0967       HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
0968       t_ieta = detId.ieta();
0969       t_iphi = detId.iphi();
0970       if (t_p > 40.0 && t_p <= 60.0)
0971         t_ietaAll->emplace_back(t_ieta);
0972     }
0973     //Selection of good track
0974     t_selectTk = spr::goodTrack(pTrack, leadPV, selectionParameter_, false);
0975     spr::trackSelectionParameters oneCutParameters = selectionParameter_;
0976     oneCutParameters.maxDxyPV = 10;
0977     oneCutParameters.maxDzPV = 100;
0978     oneCutParameters.maxInMiss = 2;
0979     oneCutParameters.maxOutMiss = 2;
0980     bool qltyFlag = spr::goodTrack(pTrack, leadPV, oneCutParameters, false);
0981     oneCutParameters = selectionParameter_;
0982     oneCutParameters.maxDxyPV = 10;
0983     oneCutParameters.maxDzPV = 100;
0984     t_qltyMissFlag = spr::goodTrack(pTrack, leadPV, oneCutParameters, false);
0985     oneCutParameters = selectionParameter_;
0986     oneCutParameters.maxInMiss = 2;
0987     oneCutParameters.maxOutMiss = 2;
0988     t_qltyPVFlag = spr::goodTrack(pTrack, leadPV, oneCutParameters, false);
0989     double eIsolation = maxRestrictionP_ * exp(slopeRestrictionP_ * std::abs((double)t_ieta));
0990     if (eIsolation < eIsolate1_)
0991       eIsolation = eIsolate1_;
0992     if (eIsolation < eIsolate2_)
0993       eIsolation = eIsolate2_;
0994 #ifdef EDM_ML_DEBUG
0995     edm::LogVerbatim("HcalIsoTrack") << "qltyFlag|okECAL|okHCAL : " << qltyFlag << "|" << trkDetItr->okECAL << "|"
0996                                      << trkDetItr->okHCAL << " eIsolation " << eIsolation;
0997 #endif
0998     t_qltyFlag = (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL);
0999     if (t_qltyFlag) {
1000       nselTracks++;
1001       int nNearTRKs(0);
1002       std::vector<DetId> eIds;
1003       std::vector<double> eHit;
1004       t_eMipDR = spr::eCone_ecal(geo,
1005                                  barrelRecHitsHandle,
1006                                  endcapRecHitsHandle,
1007                                  trkDetItr->pointHCAL,
1008                                  trkDetItr->pointECAL,
1009                                  a_mipR_,
1010                                  trkDetItr->directionECAL,
1011                                  eIds,
1012                                  eHit);
1013       double eEcal(0);
1014       for (unsigned int k = 0; k < eIds.size(); ++k) {
1015         const GlobalPoint& pos = geo->getPosition(eIds[k]);
1016         double eta = std::abs(pos.eta());
1017         if (eHit[k] > eThreshold(eIds[k], eta))
1018           eEcal += eHit[k];
1019       }
1020 #ifdef EDM_ML_DEBUG
1021       edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << t_eMipDR << ":" << eEcal;
1022 #endif
1023       t_eMipDR = eEcal;
1024       t_emaxNearP = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 15, 15);
1025       const DetId cellE(trkDetItr->detIdECAL);
1026       std::pair<double, bool> e11x11P = spr::eECALmatrix(cellE,
1027                                                          barrelRecHitsHandle,
1028                                                          endcapRecHitsHandle,
1029                                                          *theEcalChStatus,
1030                                                          geo,
1031                                                          caloTopology,
1032                                                          theEcalSevlv,
1033                                                          5,
1034                                                          5,
1035                                                          -100.0,
1036                                                          -100.0,
1037                                                          -100.0,
1038                                                          100.0);
1039       std::pair<double, bool> e15x15P = spr::eECALmatrix(cellE,
1040                                                          barrelRecHitsHandle,
1041                                                          endcapRecHitsHandle,
1042                                                          *theEcalChStatus,
1043                                                          geo,
1044                                                          caloTopology,
1045                                                          theEcalSevlv,
1046                                                          7,
1047                                                          7,
1048                                                          -100.0,
1049                                                          -100.0,
1050                                                          -100.0,
1051                                                          100.0);
1052       if (e11x11P.second && e15x15P.second) {
1053         t_eAnnular = (e15x15P.first - e11x11P.first);
1054       } else {
1055         t_eAnnular = -(e15x15P.first - e11x11P.first);
1056       }
1057       t_hmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, false);
1058       const DetId cellH(trkDetItr->detIdHCAL);
1059       double h5x5 = spr::eHCALmatrix(
1060           theHBHETopology, cellH, hbhe, 2, 2, false, true, -100.0, -100.0, -100.0, -100.0, -100.0, 100.0);
1061       double h7x7 = spr::eHCALmatrix(
1062           theHBHETopology, cellH, hbhe, 3, 3, false, true, -100.0, -100.0, -100.0, -100.0, -100.0, 100.0);
1063       t_hAnnular = h7x7 - h5x5;
1064 #ifdef EDM_ML_DEBUG
1065       edm::LogVerbatim("HcalIsoTrack") << "max p Near (Ecal) " << t_emaxNearP << " (Hcal) " << t_hmaxNearP
1066                                        << " Annular E (Ecal) " << e11x11P.first << ":" << e15x15P.first << ":"
1067                                        << t_eAnnular << " (Hcal) " << h5x5 << ":" << h7x7 << ":" << t_hAnnular;
1068 #endif
1069       t_gentrackP = trackP(pTrack, genParticles);
1070       t_gentrackE = trackE(pTrack, genParticles);
1071       if (t_eMipDR < eEcalMax_ && t_hmaxNearP < eIsolation) {
1072         t_DetIds->clear();
1073         t_HitEnergies->clear();
1074         t_DetIds1->clear();
1075         t_HitEnergies1->clear();
1076         t_DetIds3->clear();
1077         t_HitEnergies3->clear();
1078         t_DetIdEC->clear();
1079         t_HitEnergyEC->clear();
1080         t_HitDistEC->clear();
1081         t_DetIdHC->clear();
1082         t_HitEnergyHC->clear();
1083         t_HitDistHC->clear();
1084         t_mapP->clear();
1085         t_mapPt->clear();
1086         t_mapEta->clear();
1087         t_mapPhi->clear();
1088 
1089         int nRecHits(-999), nRecHits1(-999), nRecHits3(-999);
1090         std::vector<DetId> ids, ids1, ids3;
1091         std::vector<double> edet0, edet1, edet3;
1092         t_eHcal = spr::eCone_hcal(geo,
1093                                   hbhe,
1094                                   trkDetItr->pointHCAL,
1095                                   trkDetItr->pointECAL,
1096                                   a_coneR_,
1097                                   trkDetItr->directionHCAL,
1098                                   nRecHits,
1099                                   ids,
1100                                   edet0,
1101                                   useRaw_);
1102         storeEnergy(0, respCorrs, ids, edet0, t_eHcal, t_DetIds, t_HitEnergies);
1103 
1104         //----- hcal energy in the extended cone 1 (a_coneR+10) --------------
1105         t_eHcal10 = spr::eCone_hcal(geo,
1106                                     hbhe,
1107                                     trkDetItr->pointHCAL,
1108                                     trkDetItr->pointECAL,
1109                                     a_coneR1_,
1110                                     trkDetItr->directionHCAL,
1111                                     nRecHits1,
1112                                     ids1,
1113                                     edet1,
1114                                     useRaw_);
1115         storeEnergy(1, respCorrs, ids1, edet1, t_eHcal10, t_DetIds1, t_HitEnergies1);
1116 
1117         //----- hcal energy in the extended cone 3 (a_coneR+30) --------------
1118         t_eHcal30 = spr::eCone_hcal(geo,
1119                                     hbhe,
1120                                     trkDetItr->pointHCAL,
1121                                     trkDetItr->pointECAL,
1122                                     a_coneR2_,
1123                                     trkDetItr->directionHCAL,
1124                                     nRecHits3,
1125                                     ids3,
1126                                     edet3,
1127                                     useRaw_);
1128         storeEnergy(3, respCorrs, ids3, edet3, t_eHcal30, t_DetIds3, t_HitEnergies3);
1129 
1130         fillECALmatrix(
1131             cellE, barrelRecHitsHandle, endcapRecHitsHandle, geo, caloTopology, pTrack, matrixECAL_, matrixECAL_);
1132         fillHCALmatrix(theHBHETopology, cellH, hbhe, geo, pTrack, matrixHCAL_, matrixHCAL_);
1133 
1134         TrackMap(nTracks, trkCaloDirections, mapR_);
1135 
1136         t_p = pTrack->p();
1137         t_pt = pTrack->pt();
1138         t_phi = pTrack->phi();
1139 
1140 #ifdef EDM_ML_DEBUG
1141         edm::LogVerbatim("HcalIsoTrack") << "This track : " << nTracks << " (pt|eta|phi|p) :" << t_pt << "|"
1142                                          << pTrack->eta() << "|" << t_phi << "|" << t_p << " Generator Level p "
1143                                          << t_gentrackP;
1144         edm::LogVerbatim("HcalIsoTrack") << "e_MIP " << t_eMipDR << " Chg Isolation " << t_hmaxNearP << " eHcal"
1145                                          << t_eHcal << " ieta " << t_ieta << " Quality " << t_qltyMissFlag << ":"
1146                                          << t_qltyPVFlag << ":" << t_selectTk;
1147         for (unsigned int ll = 0; ll < t_DetIds->size(); ll++) {
1148           edm::LogVerbatim("HcalIsoTrack")
1149               << "det id is = " << t_DetIds->at(ll) << "   hit enery is  = " << t_HitEnergies->at(ll);
1150         }
1151         for (unsigned int ll = 0; ll < t_DetIds1->size(); ll++) {
1152           edm::LogVerbatim("HcalIsoTrack")
1153               << "det id is = " << t_DetIds1->at(ll) << "   hit enery is  = " << t_HitEnergies1->at(ll);
1154         }
1155         for (unsigned int ll = 0; ll < t_DetIds3->size(); ll++)
1156           edm::LogVerbatim("HcalIsoTrack")
1157               << "det id is = " << t_DetIds3->at(ll) << "   hit enery is  = " << t_HitEnergies3->at(ll);
1158         edm::LogVerbatim("HcalIsoTrack") << t_DetIdEC->size() << ":" << t_DetIdHC->size() << " hits from "
1159                                          << "ECAL and HCAL";
1160         for (unsigned int ll = 0; ll < t_DetIdEC->size(); ++ll)
1161           edm::LogVerbatim("HcalIsoTrack") << "EC[" << ll << "] ID " << std::hex << t_DetIdEC->at(ll) << std::dec
1162                                            << " D " << t_HitDistEC->at(ll) << " E " << t_HitEnergyEC->at(ll);
1163         for (unsigned int ll = 0; ll < t_DetIdHC->size(); ++ll)
1164           edm::LogVerbatim("HcalIsoTrack") << "HC[" << ll << "] ID " << std::hex << t_DetIdHC->at(ll) << std::dec
1165                                            << " D " << t_HitDistHC->at(ll) << " E " << t_HitEnergyHC->at(ll);
1166 #endif
1167         bool accept(false);
1168         if (t_p > pTrackMin_) {
1169           if (t_p < pTrackLow_) {
1170             ++nLow_;
1171             if (prescaleLow_ <= 1)
1172               accept = true;
1173             else if (nLow_ % prescaleLow_ == 1)
1174               accept = true;
1175           } else if (t_p > pTrackHigh_) {
1176             ++nHigh_;
1177             if (prescaleHigh_ <= 1)
1178               accept = true;
1179             else if (nHigh_ % prescaleHigh_ == 1)
1180               accept = true;
1181           } else {
1182             accept = true;
1183           }
1184         }
1185         if (accept) {
1186           tree_->Fill();
1187           nSave++;
1188           int type(0);
1189           if (t_eMipDR < 1.0) {
1190             if (t_hmaxNearP < eIsolate2_) {
1191               ++nLoose;
1192               type = 1;
1193             }
1194             if (t_hmaxNearP < eIsolate1_) {
1195               ++nTight;
1196               type = 2;
1197             }
1198           }
1199           if (t_p > 40.0 && t_p <= 60.0 && t_selectTk) {
1200             t_ietaGood->emplace_back(t_ieta);
1201             t_trackType->emplace_back(type);
1202           }
1203 #ifdef EDM_ML_DEBUG
1204           for (unsigned int k = 0; k < t_trgbits->size(); k++) {
1205             edm::LogVerbatim("HcalIsoTrack") << "trigger bit is  = " << t_trgbits->at(k);
1206           }
1207 #endif
1208         }
1209       }
1210     }
1211   }
1212   std::array<int, 3> i3{{nSave, nLoose, nTight}};
1213   return i3;
1214 }
1215 
1216 double HcalIsoTrackStudy::dR(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
1217   return reco::deltaR(vec1.eta(), vec1.phi(), vec2.eta(), vec2.phi());
1218 }
1219 
1220 double HcalIsoTrackStudy::trackP(const reco::Track* pTrack,
1221                                  const edm::Handle<reco::GenParticleCollection>& genParticles) {
1222   double pmom = -1.0;
1223   if (genParticles.isValid()) {
1224     double mindR(999.9);
1225     for (const auto& p : (*genParticles)) {
1226       double dR = reco::deltaR(pTrack->eta(), pTrack->phi(), p.momentum().Eta(), p.momentum().Phi());
1227       if (dR < mindR) {
1228         mindR = dR;
1229         pmom = p.momentum().R();
1230 #ifdef EDM_ML_DEBUG
1231         edm::LogVerbatim("HcalIsoTrack") << "p.E() :" << p.energy() << "   p.p() :" << p.momentum().R()
1232                                          << " p.M  :" << p.mass();
1233 #endif
1234       }
1235     }
1236   }
1237   return pmom;
1238 }
1239 
1240 double HcalIsoTrackStudy::trackE(const reco::Track* pTrack,
1241                                  const edm::Handle<reco::GenParticleCollection>& genParticles) {
1242   double pE = -1.0;
1243   if (genParticles.isValid()) {
1244     double mindR(999.9);
1245     for (const auto& p : (*genParticles)) {
1246       double dR = reco::deltaR(pTrack->eta(), pTrack->phi(), p.momentum().Eta(), p.momentum().Phi());
1247       if (dR < mindR) {
1248         mindR = dR;
1249         pE = p.energy();
1250 #ifdef EDM_ML_DEBUG
1251         edm::LogVerbatim("HcalIsoTrack") << "p.E() :" << p.energy() << "   p.p() :" << p.momentum().R()
1252                                          << " p.M  :" << p.mass();
1253 #endif
1254       }
1255     }
1256   }
1257   return pE;
1258 }
1259 
1260 double HcalIsoTrackStudy::rhoh(const edm::Handle<CaloTowerCollection>& tower) {
1261   std::vector<double> sumPFNallSMDQH2;
1262   sumPFNallSMDQH2.reserve(phibins_.size() * etabins_.size());
1263 
1264   for (auto eta : etabins_) {
1265     for (auto phi : phibins_) {
1266       double hadder = 0;
1267       for (const auto& pf_it : (*tower)) {
1268         if (fabs(eta - pf_it.eta()) > etahalfdist_)
1269           continue;
1270         if (fabs(reco::deltaPhi(phi, pf_it.phi())) > phihalfdist_)
1271           continue;
1272         hadder += pf_it.hadEt();
1273       }
1274       sumPFNallSMDQH2.emplace_back(hadder);
1275     }
1276   }
1277 
1278   double evt_smdq(0);
1279   std::sort(sumPFNallSMDQH2.begin(), sumPFNallSMDQH2.end());
1280   if (sumPFNallSMDQH2.size() % 2)
1281     evt_smdq = sumPFNallSMDQH2[(sumPFNallSMDQH2.size() - 1) / 2];
1282   else
1283     evt_smdq = (sumPFNallSMDQH2[sumPFNallSMDQH2.size() / 2] + sumPFNallSMDQH2[(sumPFNallSMDQH2.size() - 2) / 2]) / 2.;
1284   double rhoh = evt_smdq / (etadist_ * phidist_);
1285 #ifdef EDM_ML_DEBUG
1286   edm::LogVerbatim("HcalIsoTrack") << "Rho " << evt_smdq << ":" << rhoh;
1287 #endif
1288   return rhoh;
1289 }
1290 
1291 void HcalIsoTrackStudy::storeEnergy(int indx,
1292                                     const HcalRespCorrs* respCorrs,
1293                                     const std::vector<DetId>& ids,
1294                                     std::vector<double>& edet,
1295                                     double& eHcal,
1296                                     std::vector<unsigned int>* detIds,
1297                                     std::vector<double>* hitEnergies) {
1298   double ehcal(0);
1299   if (unCorrect_) {
1300     for (unsigned int k = 0; k < ids.size(); ++k) {
1301       double corr = (respCorrs->getValues(ids[k]))->getValue();
1302       if (corr != 0)
1303         edet[k] /= corr;
1304       ehcal += edet[k];
1305     }
1306   } else {
1307     for (const auto& en : edet)
1308       ehcal += en;
1309   }
1310   if (std::abs(ehcal - eHcal) > 0.001)
1311     edm::LogWarning("HcalIsoTrack") << "Check inconsistent energies: " << indx << " " << eHcal << ":" << ehcal
1312                                     << " from " << ids.size() << " cells";
1313   eHcal = hcalScale_ * ehcal;
1314 
1315   if (collapseDepth_) {
1316     std::map<HcalDetId, double> hitMap;
1317     for (unsigned int k = 0; k < ids.size(); ++k) {
1318       HcalDetId id = hdc_->mergedDepthDetId(HcalDetId(ids[k]));
1319       auto itr = hitMap.find(id);
1320       if (itr == hitMap.end()) {
1321         hitMap[id] = edet[k];
1322       } else {
1323         (itr->second) += edet[k];
1324       }
1325     }
1326     detIds->reserve(hitMap.size());
1327     hitEnergies->reserve(hitMap.size());
1328     for (const auto& hit : hitMap) {
1329       detIds->emplace_back(hit.first.rawId());
1330       hitEnergies->emplace_back(hit.second);
1331     }
1332   } else {
1333     detIds->reserve(ids.size());
1334     hitEnergies->reserve(ids.size());
1335     for (unsigned int k = 0; k < ids.size(); ++k) {
1336       detIds->emplace_back(ids[k].rawId());
1337       hitEnergies->emplace_back(edet[k]);
1338     }
1339   }
1340 #ifdef EDM_ML_DEBUG
1341   edm::LogVerbatim("HcalIsoTrack") << "Input to storeEnergy with " << ids.size() << " cells";
1342   for (unsigned int k = 0; k < ids.size(); ++k)
1343     edm::LogVerbatim("HcalIsoTrack") << "Hit [" << k << "] " << HcalDetId(ids[k]) << " E " << edet[k];
1344   edm::LogVerbatim("HcalIsoTrack") << "Output of storeEnergy with " << detIds->size() << " cells and Etot " << eHcal;
1345   for (unsigned int k = 0; k < detIds->size(); ++k)
1346     edm::LogVerbatim("HcalIsoTrack") << "Hit [" << k << "] " << HcalDetId((*detIds)[k]) << " E " << (*hitEnergies)[k];
1347 #endif
1348 }
1349 
1350 void HcalIsoTrackStudy::fillECALmatrix(const DetId& detId,
1351                                        edm::Handle<EcalRecHitCollection>& hitsEB,
1352                                        edm::Handle<EcalRecHitCollection>& hitsEE,
1353                                        const CaloGeometry* geo,
1354                                        const CaloTopology* caloTopology,
1355                                        const reco::Track* pTrack,
1356                                        int ieta,
1357                                        int iphi) {
1358   //Find the position of the central hit
1359   std::vector<EcalRecHitCollection::const_iterator> hits;
1360 
1361   GlobalPoint center;
1362   if (detId.subdetId() == EcalEndcap) {
1363     EEDetId EEid = EEDetId(detId);
1364     center = geo->getPosition(EEid);
1365   } else if (detId.subdetId() == EcalBarrel) {
1366     EBDetId EBid = EBDetId(detId);
1367     center = geo->getPosition(EBid);
1368   }
1369   std::vector<DetId> vdets;
1370   spr::matrixECALIds(detId, ieta, iphi, geo, caloTopology, vdets, false);
1371 
1372 #ifdef EDM_ML_DEBUG
1373   edm::LogVerbatim("HcalIsoTrack") << "Inside the fill func |||||||||||||||||"
1374                                    << "center.x(): " << center.x() << "| center.y(): " << center.y()
1375                                    << "| center.z(): " << center.z() << "\n"
1376                                    << "Size of detid vector: " << vdets.size();
1377 #endif
1378   int i = 0;
1379   for (auto const& id : vdets) {
1380     i++;
1381 #ifdef EDM_ML_DEBUG
1382     edm::LogVerbatim("HcalIsoTrack") << "counter inside vdet:" << i;
1383 #endif
1384     hits.clear();
1385     if (id.subdetId() == EcalBarrel) {
1386       spr::findHit(hitsEB, id, hits, false);
1387     } else if (id.subdetId() == EcalEndcap) {
1388       spr::findHit(hitsEE, id, hits, false);
1389     }
1390     if (i == 1 && hits.empty()) {
1391       t_DetIdEC->emplace_back(id);
1392       t_HitEnergyEC->emplace_back(0.0);
1393       t_HitDistEC->emplace_back(0.0);
1394     }
1395     if (!hits.empty()) {
1396       for (auto const& hit : hits) {
1397         GlobalPoint pos = spr::getGpos(geo, hit);
1398         GlobalVector v1(pos.x() - center.x(), pos.y() - center.y(), pos.z() - center.z());
1399         double dist = v1.mag();
1400         if (get2Ddist_ && (dist > 0)) {
1401           double cth = ((v1.x() * pTrack->px() + v1.y() * pTrack->py() + v1.z() * pTrack->pz()) / (dist * pTrack->p()));
1402           dist *= std::sqrt(1.0 - cth * cth);
1403         }
1404         double ener = hit->energy();
1405         t_DetIdEC->emplace_back(hit->id());
1406         t_HitEnergyEC->emplace_back(ener);
1407         t_HitDistEC->emplace_back(dist);
1408       }
1409     }
1410   }
1411 }
1412 
1413 void HcalIsoTrackStudy::fillHCALmatrix(const HcalTopology* topology,
1414                                        const DetId& detId,
1415                                        edm::Handle<HBHERecHitCollection>& hbhe,
1416                                        const CaloGeometry* geo,
1417                                        const reco::Track* pTrack,
1418                                        int ieta,
1419                                        int iphi) {
1420   //Find the position of the central hit
1421   std::vector<HBHERecHitCollection::const_iterator> hits;
1422   spr::hitHCALmatrix(topology, detId, hbhe, ieta, iphi, hits, false, false);
1423   GlobalPoint center = (static_cast<const HcalGeometry*>(geo->getSubdetectorGeometry(detId)))->getPosition(detId);
1424   for (auto const& hit : hits) {
1425     GlobalPoint pos = spr::getGpos(geo, hit);
1426     GlobalVector v1(pos.x() - center.x(), pos.y() - center.y(), pos.z() - center.z());
1427     double dist = v1.mag();
1428     if (get2Ddist_ && (dist > 0)) {
1429       double cth = ((v1.x() * pTrack->px() + v1.y() * pTrack->py() + v1.z() * pTrack->pz()) / (dist * pTrack->p()));
1430       dist *= std::sqrt(1.0 - cth * cth);
1431     }
1432     double ener = spr::getEnergy(hit, useRaw_);
1433     t_DetIdHC->emplace_back(hit->id());
1434     t_HitEnergyHC->emplace_back(ener);
1435     t_HitDistHC->emplace_back(dist);
1436   }
1437 }
1438 
1439 void HcalIsoTrackStudy::TrackMap(unsigned int trkIndex,
1440                                  std::vector<spr::propagatedTrackDirection>& trkDirs,
1441                                  double dR) {
1442   if (trkDirs[trkIndex].okHCAL) {
1443     for (unsigned int indx = 0; indx < trkDirs.size(); ++indx) {
1444       if (trkDirs[indx].ok && trkDirs[indx].okHCAL) {
1445         int isConeChargedIso = spr::coneChargeIsolation(
1446             trkDirs[trkIndex].pointHCAL, trkDirs[indx].pointHCAL, trkDirs[trkIndex].directionHCAL, dR);
1447         if (isConeChargedIso == 0) {
1448           const reco::Track* pTrack = &(*(trkDirs[indx].trkItr));
1449           TVector3 point(trkDirs[indx].pointHCAL.x(), trkDirs[indx].pointHCAL.y(), trkDirs[indx].pointHCAL.z());
1450 
1451           t_mapP->emplace_back(pTrack->p());
1452           t_mapPt->emplace_back(pTrack->pt());
1453           t_mapEta->emplace_back(point.Eta());
1454           t_mapPhi->emplace_back(point.Phi());
1455         }
1456       }
1457     }
1458   }
1459 }
1460 
1461 double HcalIsoTrackStudy::eThreshold(const DetId& id, double eta) const {
1462   double eThr(hitEthrEB_);
1463   if (usePFThresh_) {
1464     eThr = static_cast<double>((*eThresholds_)[id]);
1465   } else {
1466     if (id.subdetId() != EcalBarrel) {
1467       eThr = (((eta * hitEthrEE3_ + hitEthrEE2_) * eta + hitEthrEE1_) * eta + hitEthrEE0_);
1468       if (eThr < hitEthrEELo_)
1469         eThr = hitEthrEELo_;
1470       else if (eThr > hitEthrEEHi_)
1471         eThr = hitEthrEEHi_;
1472     }
1473   }
1474   return eThr;
1475 }
1476 
1477 //define this as a plug-in
1478 DEFINE_FWK_MODULE(HcalIsoTrackStudy);