Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-11-28 03:03:22

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