Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-05-09 22:37:14 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!