Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 09:34:22 

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