Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-28 01:32:03

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