Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-01-18 23:28:32 

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 alcaHcalIsoTracks {
0087   struct Counters {
0088     Counters() : nAll_(0), nGood_(0), nRange_(0) {}
0089     mutable std::atomic<unsigned int> nAll_, nGood_, nRange_;
0090   };
0091 }  // namespace alcaHcalIsoTracks
0092 
0093 class AlCaHcalIsotrkProducer : public edm::stream::EDProducer<edm::GlobalCache<alcaHcalIsoTracks::Counters>> {
0094 public:
0095   explicit AlCaHcalIsotrkProducer(edm::ParameterSet const&, const alcaHcalIsoTracks::Counters*);
0096   ~AlCaHcalIsotrkProducer() override = default;
0097 
0098   static std::unique_ptr<alcaHcalIsoTracks::Counters> initializeGlobalCache(edm::ParameterSet const&) {
0099     return std::make_unique<alcaHcalIsoTracks::Counters>();
0100   }
0101 
0102   void produce(edm::Event&, edm::EventSetup const&) override;
0103   void endStream() override;
0104   static void globalEndJob(const alcaHcalIsoTracks::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 alcaHcalIsoTracks::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               std::string label = triggerEvent.filterTag(ifilter).label();
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]) != std::string::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 alcaHcalIsoTracks::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), nselTracks(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       nselTracks++;
0888       int nNearTRKs(0);
0889       ////////////////////////////////-MIP STUFF-//////////////////////////////
0890       std::vector<DetId> eIds;
0891       std::vector<double> eHit;
0892 #ifdef EDM_ML_DEBUG
0893       double eMipDR =
0894 #endif
0895           spr::eCone_ecal(geo,
0896                           barrelRecHitsHandle,
0897                           endcapRecHitsHandle,
0898                           trkDetItr->pointHCAL,
0899                           trkDetItr->pointECAL,
0900                           a_mipR_,
0901                           trkDetItr->directionECAL,
0902                           eIds,
0903                           eHit);
0904       double eEcal(0);
0905       for (unsigned int k = 0; k < eIds.size(); ++k) {
0906         if (eHit[k] > eThreshold(eIds[k], geo))
0907           eEcal += eHit[k];
0908       }
0909 #ifdef EDM_ML_DEBUG
0910       if (debug_)
0911         edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR << ":" << eEcal;
0912 #endif
0913       isoTk.eMipDR_.emplace_back(eEcal);
0914       ////////////////////////////////-MIP STUFF-///////////////////////////////
0915       ////////////////////////////////-MIP STUFF-2//////////////////////////////
0916       std::vector<DetId> eIds2;
0917       std::vector<double> eHit2;
0918 #ifdef EDM_ML_DEBUG
0919       double eMipDR2 =
0920 #endif
0921           spr::eCone_ecal(geo,
0922                           barrelRecHitsHandle,
0923                           endcapRecHitsHandle,
0924                           trkDetItr->pointHCAL,
0925                           trkDetItr->pointECAL,
0926                           a_mipR2_,
0927                           trkDetItr->directionECAL,
0928                           eIds2,
0929                           eHit2);
0930       double eEcal2(0);
0931       for (unsigned int k = 0; k < eIds2.size(); ++k) {
0932         if (eHit2[k] > eThreshold(eIds2[k], geo))
0933           eEcal2 += eHit2[k];
0934       }
0935 #ifdef EDM_ML_DEBUG
0936       if (debug_)
0937         edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR2 << ":" << eEcal2;
0938 #endif
0939       isoTk.eMipDR_.emplace_back(eEcal2);
0940       ////////////////////////////////-MIP STUFF-2/////////////////////////////
0941       ////////////////////////////////-MIP STUFF-3/////////////////////////////
0942       std::vector<DetId> eIds3;
0943       std::vector<double> eHit3;
0944 #ifdef EDM_ML_DEBUG
0945       double eMipDR3 =
0946 #endif
0947           spr::eCone_ecal(geo,
0948                           barrelRecHitsHandle,
0949                           endcapRecHitsHandle,
0950                           trkDetItr->pointHCAL,
0951                           trkDetItr->pointECAL,
0952                           a_mipR3_,
0953                           trkDetItr->directionECAL,
0954                           eIds3,
0955                           eHit3);
0956       double eEcal3(0);
0957       for (unsigned int k = 0; k < eIds3.size(); ++k) {
0958         if (eHit3[k] > eThreshold(eIds3[k], geo))
0959           eEcal3 += eHit3[k];
0960       }
0961 #ifdef EDM_ML_DEBUG
0962       if (debug_)
0963         edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR3 << ":" << eEcal3;
0964 #endif
0965       isoTk.eMipDR_.emplace_back(eEcal3);
0966       ////////////////////////////////-MIP STUFF-3/////////////////////////////
0967       ////////////////////////////////-MIP STUFF-4/////////////////////////////
0968       std::vector<DetId> eIds4;
0969       std::vector<double> eHit4;
0970 #ifdef EDM_ML_DEBUG
0971       double eMipDR4 =
0972 #endif
0973           spr::eCone_ecal(geo,
0974                           barrelRecHitsHandle,
0975                           endcapRecHitsHandle,
0976                           trkDetItr->pointHCAL,
0977                           trkDetItr->pointECAL,
0978                           a_mipR4_,
0979                           trkDetItr->directionECAL,
0980                           eIds4,
0981                           eHit4);
0982       double eEcal4(0);
0983       for (unsigned int k = 0; k < eIds4.size(); ++k) {
0984         if (eHit4[k] > eThreshold(eIds4[k], geo))
0985           eEcal4 += eHit4[k];
0986       }
0987 #ifdef EDM_ML_DEBUG
0988       if (debug_)
0989         edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR4 << ":" << eEcal4;
0990 #endif
0991       isoTk.eMipDR_.emplace_back(eEcal4);
0992       ////////////////////////////////-MIP STUFF-4/////////////////////////////
0993       ////////////////////////////////-MIP STUFF-5/////////////////////////////
0994       std::vector<DetId> eIds5;
0995       std::vector<double> eHit5;
0996 #ifdef EDM_ML_DEBUG
0997       double eMipDR5 =
0998 #endif
0999           spr::eCone_ecal(geo,
1000                           barrelRecHitsHandle,
1001                           endcapRecHitsHandle,
1002                           trkDetItr->pointHCAL,
1003                           trkDetItr->pointECAL,
1004                           a_mipR5_,
1005                           trkDetItr->directionECAL,
1006                           eIds5,
1007                           eHit5);
1008       double eEcal5(0);
1009       for (unsigned int k = 0; k < eIds5.size(); ++k) {
1010         if (eHit5[k] > eThreshold(eIds5[k], geo))
1011           eEcal5 += eHit5[k];
1012       }
1013 #ifdef EDM_ML_DEBUG
1014       if (debug_)
1015         edm::LogVerbatim("HcalIsoTrack") << "eMIP before and after: " << eMipDR5 << ":" << eEcal5;
1016 #endif
1017       isoTk.eMipDR_.emplace_back(eEcal5);
1018       ////////////////////////////////-MIP STUFF-5/////////////////////////////
1019 
1020       isoTk.emaxNearP_ = spr::chargeIsolationEcal(nTracks, trkCaloDets, geo, caloTopology, 15, 15);
1021       const DetId cellE(trkDetItr->detIdECAL);
1022       std::pair<double, bool> e11x11P = spr::eECALmatrix(cellE,
1023                                                          barrelRecHitsHandle,
1024                                                          endcapRecHitsHandle,
1025                                                          *theEcalChStatus,
1026                                                          geo,
1027                                                          caloTopology,
1028                                                          theEcalSevlv,
1029                                                          5,
1030                                                          5,
1031                                                          -100.0,
1032                                                          -100.0,
1033                                                          -100.0,
1034                                                          100.0);
1035       std::pair<double, bool> e15x15P = spr::eECALmatrix(cellE,
1036                                                          barrelRecHitsHandle,
1037                                                          endcapRecHitsHandle,
1038                                                          *theEcalChStatus,
1039                                                          geo,
1040                                                          caloTopology,
1041                                                          theEcalSevlv,
1042                                                          7,
1043                                                          7,
1044                                                          -100.0,
1045                                                          -100.0,
1046                                                          -100.0,
1047                                                          100.0);
1048       if (e11x11P.second && e15x15P.second) {
1049         isoTk.eAnnular_ = (e15x15P.first - e11x11P.first);
1050       } else {
1051         isoTk.eAnnular_ = -(e15x15P.first - e11x11P.first);
1052       }
1053       isoTk.hmaxNearP_ = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, false);
1054       const DetId cellH(trkDetItr->detIdHCAL);
1055       double h5x5 = spr::eHCALmatrix(
1056           theHBHETopology, cellH, hbhe, 2, 2, false, true, -100.0, -100.0, -100.0, -100.0, -100.0, 100.0);
1057       double h7x7 = spr::eHCALmatrix(
1058           theHBHETopology, cellH, hbhe, 3, 3, false, true, -100.0, -100.0, -100.0, -100.0, -100.0, 100.0);
1059       isoTk.hAnnular_ = h7x7 - h5x5;
1060 #ifdef EDM_ML_DEBUG
1061       if (debug_)
1062         edm::LogVerbatim("HcalIsoTrack") << "max p Near (Ecal) " << isoTk.emaxNearP_ << " (Hcal) " << isoTk.hmaxNearP_
1063                                          << " Annular E (Ecal) " << e11x11P.first << ":" << e15x15P.first << ":"
1064                                          << isoTk.eAnnular_ << " (Hcal) " << h5x5 << ":" << h7x7 << ":"
1065                                          << isoTk.hAnnular_;
1066       if (isoTk.eMipDR_[0] < eEcalMax_)
1067         flag += 16;
1068       if (isoTk.hmaxNearP_ < eIsolation)
1069         flag += 32;
1070 #endif
1071       isoTk.gentrackP_ = trackP(pTrack, genParticles);
1072       if (isoTk.eMipDR_[0] < eEcalMax_ && isoTk.hmaxNearP_ < eIsolation) {
1073         int nRecHits(-999), nRecHits1(-999), nRecHits3(-999);
1074         std::vector<DetId> ids, ids1, ids3;
1075         std::vector<double> edet0, edet1, edet3;
1076         isoTk.eHcal_ = spr::eCone_hcal(geo,
1077                                        hbhe,
1078                                        trkDetItr->pointHCAL,
1079                                        trkDetItr->pointECAL,
1080                                        a_coneR_,
1081                                        trkDetItr->directionHCAL,
1082                                        nRecHits,
1083                                        ids,
1084                                        edet0,
1085                                        0);
1086         if (!oldID_.empty()) {
1087           for (unsigned k = 0; k < ids.size(); ++k)
1088             ids[k] = newId(ids[k]);
1089         }
1090         storeEnergy(respCorrs, ids, edet0, isoTk.eHcal_, isoTk.detIds_, isoTk.hitEnergies_);
1091         std::pair<double, double> ehcal0 =
1092             storeEnergy(respCorrs, hbhe, ids, isoTk.hitEnergiesRaw_, isoTk.hitEnergiesAux_);
1093         isoTk.eHcalRaw_ = ehcal0.first;
1094         isoTk.eHcalAux_ = ehcal0.second;
1095 
1096         //----- hcal energy in the extended cone 1 (a_coneR+10) --------------
1097         isoTk.eHcal10_ = spr::eCone_hcal(geo,
1098                                          hbhe,
1099                                          trkDetItr->pointHCAL,
1100                                          trkDetItr->pointECAL,
1101                                          a_coneR1_,
1102                                          trkDetItr->directionHCAL,
1103                                          nRecHits1,
1104                                          ids1,
1105                                          edet1,
1106                                          0);
1107         if (!oldID_.empty()) {
1108           for (unsigned k = 0; k < ids1.size(); ++k)
1109             ids1[k] = newId(ids1[k]);
1110         }
1111         storeEnergy(respCorrs, ids1, edet1, isoTk.eHcal10_, isoTk.detIds1_, isoTk.hitEnergies1_);
1112         std::pair<double, double> ehcal1 =
1113             storeEnergy(respCorrs, hbhe, ids1, isoTk.hitEnergies1Raw_, isoTk.hitEnergies1Aux_);
1114         isoTk.eHcal10Raw_ = ehcal1.first;
1115         isoTk.eHcal10Aux_ = ehcal1.second;
1116 
1117         //----- hcal energy in the extended cone 3 (a_coneR+30) --------------
1118         isoTk.eHcal30_ = spr::eCone_hcal(geo,
1119                                          hbhe,
1120                                          trkDetItr->pointHCAL,
1121                                          trkDetItr->pointECAL,
1122                                          a_coneR2_,
1123                                          trkDetItr->directionHCAL,
1124                                          nRecHits3,
1125                                          ids3,
1126                                          edet3,
1127                                          0);
1128         if (!oldID_.empty()) {
1129           for (unsigned k = 0; k < ids3.size(); ++k)
1130             ids3[k] = newId(ids3[k]);
1131         }
1132         storeEnergy(respCorrs, ids3, edet3, isoTk.eHcal30_, isoTk.detIds3_, isoTk.hitEnergies3_);
1133         std::pair<double, double> ehcal3 =
1134             storeEnergy(respCorrs, hbhe, ids3, isoTk.hitEnergies3Raw_, isoTk.hitEnergies3Aux_);
1135         isoTk.eHcal30Raw_ = ehcal3.first;
1136         isoTk.eHcal30Aux_ = ehcal3.second;
1137 
1138         if (isoTk.p_ > pTrackMin_)
1139           accept = true;
1140 #ifdef EDM_ML_DEBUG
1141         if (accept)
1142           flag += 64;
1143         if (debug_) {
1144           std::string ctype = accept ? " ***** ACCEPT *****" : "";
1145           edm::LogVerbatim("HcalIsoTrack")
1146               << "This track : " << nTracks << " (pt|eta|phi|p) : " << isoTk.pt_ << "|" << pTrack->eta() << "|"
1147               << isoTk.phi_ << "|" << isoTk.p_ << " Generator Level p " << isoTk.gentrackP_;
1148           edm::LogVerbatim("HcalIsoTrack")
1149               << "e_MIP " << isoTk.eMipDR_[0] << " Chg Isolation " << isoTk.hmaxNearP_ << " eHcal" << isoTk.eHcal_
1150               << ":" << isoTk.eHcalRaw_ << ":" << isoTk.eHcalAux_ << " ieta " << isoTk.ieta_ << " Quality "
1151               << isoTk.qltyMissFlag_ << ":" << isoTk.qltyPVFlag_ << ":" << isoTk.selectTk_ << ctype;
1152           for (unsigned int ll = 0; ll < isoTk.detIds_.size(); ll++) {
1153             edm::LogVerbatim("HcalIsoTrack")
1154                 << "det id is = " << HcalDetId(isoTk.detIds_[ll]) << "   hit enery is  = " << isoTk.hitEnergies_[ll]
1155                 << " : " << isoTk.hitEnergiesRaw_[ll] << " : " << isoTk.hitEnergiesAux_[ll];
1156           }
1157           for (unsigned int ll = 0; ll < isoTk.detIds1_.size(); ll++) {
1158             edm::LogVerbatim("HcalIsoTrack")
1159                 << "det id is = " << HcalDetId(isoTk.detIds1_[ll]) << "   hit enery is  = " << isoTk.hitEnergies1_[ll]
1160                 << " : " << isoTk.hitEnergies1Raw_[ll] << " : " << isoTk.hitEnergies1Aux_[ll];
1161           }
1162           for (unsigned int ll = 0; ll < isoTk.detIds3_.size(); ll++) {
1163             edm::LogVerbatim("HcalIsoTrack")
1164                 << "det id is = " << HcalDetId(isoTk.detIds3_[ll]) << "   hit enery is  = " << isoTk.hitEnergies3_[ll]
1165                 << " : " << isoTk.hitEnergies3Raw_[ll] << " : " << isoTk.hitEnergies3Aux_[ll];
1166           }
1167         }
1168 #endif
1169         if (accept) {
1170           edm::LogVerbatim("HcalIsoTrackX")
1171               << "Run " << eventId.run() << " Event " << eventId.event() << " Track " << nTracks << " p " << isoTk.p_;
1172           hocalib.emplace_back(isoTk);
1173           nSave++;
1174           int type(0);
1175           if (isoTk.eMipDR_[0] < 1.0) {
1176             if (isoTk.hmaxNearP_ < eIsolate2_) {
1177               ++nLoose;
1178               type = 1;
1179             }
1180             if (isoTk.hmaxNearP_ < eIsolate1_) {
1181               ++nTight;
1182               type = 2;
1183             }
1184           }
1185           if (isoTk.p_ > 40.0 && isoTk.p_ <= 60.0 && isoTk.selectTk_) {
1186             hocalibEvent.ietaGood_.emplace_back(isoTk.ieta_);
1187             hocalibEvent.trackType_.emplace_back(type);
1188           }
1189 #ifdef EDM_ML_DEBUG
1190           for (unsigned int k = 0; k < isoTk.trgbits_.size(); k++) {
1191             edm::LogVerbatim("HcalIsoTrack") << "trigger bit is  = " << isoTk.trgbits_[k];
1192           }
1193 #endif
1194         }
1195       }
1196     }
1197 #ifdef EDM_ML_DEBUG
1198     if (debug_) {
1199       if (isoTk.eMipDR_.empty())
1200         edm::LogVerbatim("HcalIsoTrack") << "Track " << nTracks << " Selection Flag " << std::hex << flag << std::dec
1201                                          << " Accept " << accept << " Momentum " << isoTk.p_ << ":" << pTrackMin_;
1202       else
1203         edm::LogVerbatim("HcalIsoTrack") << "Track " << nTracks << " Selection Flag " << std::hex << flag << std::dec
1204                                          << " Accept " << accept << " Momentum " << isoTk.p_ << ":" << pTrackMin_
1205                                          << " Ecal Energy " << isoTk.eMipDR_[0] << ":" << eEcalMax_
1206                                          << " Charge Isolation " << isoTk.hmaxNearP_ << ":" << eIsolation;
1207     }
1208 #endif
1209   }
1210   std::array<int, 3> i3{{nSave, nLoose, nTight}};
1211   return i3;
1212 }
1213 
1214 double AlCaHcalIsotrkProducer::dR(math::XYZTLorentzVector& vec1, math::XYZTLorentzVector& vec2) {
1215   return reco::deltaR(vec1.eta(), vec1.phi(), vec2.eta(), vec2.phi());
1216 }
1217 
1218 double AlCaHcalIsotrkProducer::trackP(const reco::Track* pTrack,
1219                                       const edm::Handle<reco::GenParticleCollection>& genParticles) {
1220   double pmom = -1.0;
1221   if (genParticles.isValid()) {
1222     double mindR(999.9);
1223     for (const auto& p : (*genParticles)) {
1224       double dR = reco::deltaR(pTrack->eta(), pTrack->phi(), p.momentum().Eta(), p.momentum().Phi());
1225       if (dR < mindR) {
1226         mindR = dR;
1227         pmom = p.momentum().R();
1228       }
1229     }
1230   }
1231   return pmom;
1232 }
1233 
1234 double AlCaHcalIsotrkProducer::rhoh(const edm::Handle<CaloTowerCollection>& tower) {
1235   std::vector<double> sumPFNallSMDQH2;
1236   sumPFNallSMDQH2.reserve(phibins_.size() * etabins_.size());
1237   for (auto eta : etabins_) {
1238     for (auto phi : phibins_) {
1239       double hadder = 0;
1240       for (const auto& pf_it : (*tower)) {
1241         if (fabs(eta - pf_it.eta()) > etahalfdist_)
1242           continue;
1243         if (fabs(reco::deltaPhi(phi, pf_it.phi())) > phihalfdist_)
1244           continue;
1245         hadder += pf_it.hadEt();
1246       }
1247       sumPFNallSMDQH2.emplace_back(hadder);
1248     }
1249   }
1250 
1251   double evt_smdq(0);
1252   std::sort(sumPFNallSMDQH2.begin(), sumPFNallSMDQH2.end());
1253   if (sumPFNallSMDQH2.size() % 2)
1254     evt_smdq = sumPFNallSMDQH2[(sumPFNallSMDQH2.size() - 1) / 2];
1255   else
1256     evt_smdq = (sumPFNallSMDQH2[sumPFNallSMDQH2.size() / 2] + sumPFNallSMDQH2[(sumPFNallSMDQH2.size() - 2) / 2]) / 2.;
1257   double rhoh = evt_smdq / (etadist_ * phidist_);
1258 #ifdef EDM_ML_DEBUG
1259   if (debug_)
1260     edm::LogVerbatim("HcalIsoTrack") << "Rho " << evt_smdq << ":" << rhoh;
1261 #endif
1262   return rhoh;
1263 }
1264 
1265 double AlCaHcalIsotrkProducer::eThreshold(const DetId& id, const CaloGeometry* geo) const {
1266   const GlobalPoint& pos = geo->getPosition(id);
1267   double eta = std::abs(pos.eta());
1268   double eThr(hitEthrEB_);
1269   if (usePFThresh_) {
1270     eThr = static_cast<double>((*eThresholds_)[id]);
1271   } else {
1272     if (id.subdetId() != EcalBarrel) {
1273       eThr = (((eta * hitEthrEE3_ + hitEthrEE2_) * eta + hitEthrEE1_) * eta + hitEthrEE0_);
1274       if (eThr < hitEthrEELo_)
1275         eThr = hitEthrEELo_;
1276       else if (eThr > hitEthrEEHi_)
1277         eThr = hitEthrEEHi_;
1278     }
1279   }
1280   return eThr;
1281 }
1282 
1283 DetId AlCaHcalIsotrkProducer::newId(const DetId& id) {
1284   HcalDetId hid(id);
1285   if (hep17_ && ((hid.iphi() < 63) || (hid.iphi() > 66) || (hid.zside() < 0)))
1286     return id;
1287   for (unsigned int k = 0; k < oldID_.size(); ++k) {
1288     if ((hid.subdetId() == oldDet_[k]) && (hid.ietaAbs() == oldEta_[k]) && (hid.depth() == oldDepth_[k])) {
1289       return static_cast<DetId>(HcalDetId(hid.subdet(), hid.ieta(), hid.iphi(), newDepth_[k]));
1290     }
1291   }
1292   return id;
1293 }
1294 
1295 void AlCaHcalIsotrkProducer::storeEnergy(const HcalRespCorrs* respCorrs,
1296                                          const std::vector<DetId>& ids,
1297                                          std::vector<double>& edet,
1298                                          double& eHcal,
1299                                          std::vector<unsigned int>& detIds,
1300                                          std::vector<double>& hitEnergies) {
1301   double ehcal(0);
1302   if (unCorrect_) {
1303     for (unsigned int k = 0; k < ids.size(); ++k) {
1304       double corr = (respCorrs->getValues(ids[k]))->getValue();
1305       if (corr != 0)
1306         edet[k] /= corr;
1307       ehcal += edet[k];
1308     }
1309   } else {
1310     for (const auto& en : edet)
1311       ehcal += en;
1312   }
1313   if ((std::abs(ehcal - eHcal) > 0.001) && (!unCorrect_))
1314     edm::LogWarning("HcalIsoTrack") << "Check inconsistent energies: " << eHcal << ":" << ehcal << " from "
1315                                     << ids.size() << " cells";
1316   eHcal = hcalScale_ * ehcal;
1317 
1318   if (collapseDepth_) {
1319     std::map<HcalDetId, double> hitMap;
1320     for (unsigned int k = 0; k < ids.size(); ++k) {
1321       HcalDetId id = hdc_->mergedDepthDetId(HcalDetId(ids[k]));
1322       auto itr = hitMap.find(id);
1323       if (itr == hitMap.end()) {
1324         hitMap[id] = edet[k];
1325       } else {
1326         (itr->second) += edet[k];
1327       }
1328     }
1329     detIds.reserve(hitMap.size());
1330     hitEnergies.reserve(hitMap.size());
1331     for (const auto& hit : hitMap) {
1332       detIds.emplace_back(hit.first.rawId());
1333       hitEnergies.emplace_back(hit.second);
1334     }
1335   } else {
1336     detIds.reserve(ids.size());
1337     hitEnergies.reserve(ids.size());
1338     for (unsigned int k = 0; k < ids.size(); ++k) {
1339       detIds.emplace_back(ids[k].rawId());
1340       hitEnergies.emplace_back(edet[k]);
1341     }
1342   }
1343 #ifdef EDM_ML_DEBUG
1344   if (debug_) {
1345     edm::LogVerbatim("HcalIsoTrack") << "StoreEnergy1::Input to storeEnergy with " << ids.size() << " cells";
1346     for (unsigned int k = 0; k < ids.size(); ++k)
1347       edm::LogVerbatim("HcalIsoTrack") << "Hit [" << k << "] " << HcalDetId(ids[k]) << " E " << edet[k];
1348     edm::LogVerbatim("HcalIsoTrack") << "StoreEnergy1::Output of storeEnergy with " << detIds.size()
1349                                      << " cells and Etot " << eHcal;
1350     for (unsigned int k = 0; k < detIds.size(); ++k)
1351       edm::LogVerbatim("HcalIsoTrack") << "Hit [" << k << "] " << HcalDetId(detIds[k]) << " E " << hitEnergies[k];
1352   }
1353 #endif
1354 }
1355 
1356 std::pair<double, double> AlCaHcalIsotrkProducer::storeEnergy(const HcalRespCorrs* respCorrs,
1357                                                               edm::Handle<HBHERecHitCollection>& hbhe,
1358                                                               const std::vector<DetId>& ids,
1359                                                               std::vector<double>& hitEnergy1,
1360                                                               std::vector<double>& hitEnergy2) {
1361   double ehcal1(0), ehcal2(0);
1362   std::vector<double> edet1, edet2;
1363   for (unsigned int k = 0; k < ids.size(); ++k) {
1364     double e1(0), e2(0);
1365     for (auto itr = hbhe->begin(); itr != hbhe->end(); ++itr) {
1366       if (itr->id() == ids[k]) {
1367         e1 = itr->eraw();
1368         e2 = itr->eaux();
1369         break;
1370       }
1371     }
1372     if (e1 < 1.e-10)
1373       e1 = 0;
1374     if (e2 < 1.e-10)
1375       e2 = 0;
1376     edet1.emplace_back(e1);
1377     edet2.emplace_back(e2);
1378   }
1379   if (unCorrect_) {
1380     for (unsigned int k = 0; k < ids.size(); ++k) {
1381       double corr = (respCorrs->getValues(ids[k]))->getValue();
1382       if (corr != 0) {
1383         edet1[k] /= corr;
1384         edet2[k] /= corr;
1385       }
1386     }
1387   }
1388   for (unsigned int k = 0; k < ids.size(); ++k) {
1389     ehcal1 += edet1[k];
1390     ehcal2 += edet2[k];
1391   }
1392   ehcal1 *= hcalScale_;
1393   ehcal2 *= hcalScale_;
1394 
1395   if (collapseDepth_) {
1396     std::map<HcalDetId, std::pair<double, double>> hitMap;
1397     for (unsigned int k = 0; k < ids.size(); ++k) {
1398       HcalDetId id = hdc_->mergedDepthDetId(HcalDetId(ids[k]));
1399       auto itr = hitMap.find(id);
1400       if (itr == hitMap.end()) {
1401         hitMap[id] = std::make_pair(edet1[k], edet2[k]);
1402       } else {
1403         (itr->second).first += edet1[k];
1404         (itr->second).second += edet2[k];
1405       }
1406     }
1407     hitEnergy1.reserve(hitMap.size());
1408     hitEnergy2.reserve(hitMap.size());
1409     for (const auto& hit : hitMap) {
1410       hitEnergy1.emplace_back(hit.second.first);
1411       hitEnergy2.emplace_back(hit.second.second);
1412     }
1413   } else {
1414     hitEnergy1.reserve(ids.size());
1415     hitEnergy2.reserve(ids.size());
1416     for (unsigned int k = 0; k < ids.size(); ++k) {
1417       hitEnergy1.emplace_back(edet1[k]);
1418       hitEnergy2.emplace_back(edet2[k]);
1419     }
1420   }
1421 #ifdef EDM_ML_DEBUG
1422   if (debug_) {
1423     edm::LogVerbatim("HcalIsoTrack") << "StoreEnergy2::Input to storeEnergy with " << ids.size() << " cells";
1424     edm::LogVerbatim("HcalIsoTrack") << "StoreEnergy2::Output of storeEnergy with " << hitEnergy1.size()
1425                                      << " cells and Etot " << ehcal1 << ":" << ehcal2;
1426     for (unsigned int k = 0; k < hitEnergy1.size(); ++k)
1427       edm::LogVerbatim("HcalIsoTrack") << "Hit [" << k << "] " << hitEnergy1[k] << " : " << hitEnergy2[k];
1428   }
1429 #endif
1430   return std::make_pair(ehcal1, ehcal2);
1431 }
1432 
1433 bool AlCaHcalIsotrkProducer::notaMuon(const reco::Track* pTrack0, const edm::Handle<reco::MuonCollection>& muonh) {
1434   bool flag(true);
1435   for (reco::MuonCollection::const_iterator recMuon = muonh->begin(); recMuon != muonh->end(); ++recMuon) {
1436     if (recMuon->innerTrack().isNonnull()) {
1437       const reco::Track* pTrack = (recMuon->innerTrack()).get();
1438       bool mediumMuon = (((recMuon->isPFMuon()) && (recMuon->isGlobalMuon() || recMuon->isTrackerMuon())) &&
1439                          (recMuon->innerTrack()->validFraction() > 0.49));
1440       if (mediumMuon) {
1441         double chiGlobal = ((recMuon->globalTrack().isNonnull()) ? recMuon->globalTrack()->normalizedChi2() : 999);
1442         bool goodGlob = (recMuon->isGlobalMuon() && chiGlobal < 3 &&
1443                          recMuon->combinedQuality().chi2LocalPosition < 12 && recMuon->combinedQuality().trkKink < 20);
1444         mediumMuon = muon::segmentCompatibility(*recMuon) > (goodGlob ? 0.303 : 0.451);
1445       }
1446       if (mediumMuon) {
1447         double dR = reco::deltaR(pTrack->eta(), pTrack->phi(), pTrack0->eta(), pTrack0->phi());
1448         if (dR < 0.1) {
1449           flag = false;
1450           break;
1451         }
1452       }
1453     }
1454   }
1455   return flag;
1456 }
1457 
1458 #include "FWCore/Framework/interface/MakerMacros.h"
1459 
1460 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!