AlCaIsoTracksProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//#define EDM_ML_DEBUG

// system include files
#include <atomic>
#include <memory>
#include <string>
#include <cmath>
#include <iostream>
#include <sstream>
#include <fstream>
#include <vector>
#include <boost/regex.hpp>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/Common/interface/TriggerNames.h"
#include "FWCore/Utilities/interface/Exception.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Common/interface/Ref.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/HcalIsolatedTrack/interface/HcalIsolatedTrackCandidate.h"
#include "DataFormats/HLTReco/interface/TriggerEvent.h"
#include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
#include "DataFormats/L1Trigger/interface/L1JetParticleFwd.h"
#include "DataFormats/L1Trigger/interface/L1JetParticle.h"
#include "DataFormats/Provenance/interface/ProductID.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/TrackExtra.h"
#include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "DataFormats/TrackReco/interface/TrackBase.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
#include "DataFormats/Common/interface/TriggerResults.h"

#include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"

#include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
#include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
#include "Calibration/IsolatedParticles/interface/eCone.h"
#include "Calibration/IsolatedParticles/interface/TrackSelection.h"

#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/Records/interface/CaloTopologyRecord.h"
#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"
#include "Geometry/CaloTopology/interface/CaloTopology.h"
#include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
#include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
#include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"

//#define EDM_ML_DEBUG
//
// class declaration
//

namespace AlCaIsoTracks {
  struct Counters {
    Counters() : nAll_(0), nGood_(0), nRange_(0) {}
    mutable std::atomic<unsigned int> nAll_, nGood_, nRange_;
  };
}  // namespace AlCaIsoTracks

class AlCaIsoTracksProducer : public edm::stream::EDProducer<edm::GlobalCache<AlCaIsoTracks::Counters> > {
public:
  explicit AlCaIsoTracksProducer(edm::ParameterSet const&, const AlCaIsoTracks::Counters* count);
  ~AlCaIsoTracksProducer() override;

  static std::unique_ptr<AlCaIsoTracks::Counters> initializeGlobalCache(edm::ParameterSet const&) {
    return std::make_unique<AlCaIsoTracks::Counters>();
  }

  void produce(edm::Event&, edm::EventSetup const&) override;
  void endStream() override;
  static void globalEndJob(const AlCaIsoTracks::Counters* counters);
  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
  void beginRun(edm::Run const&, edm::EventSetup const&) override;
  void endRun(edm::Run const&, edm::EventSetup const&) override;
  reco::HcalIsolatedTrackCandidateCollection* select(edm::Handle<edm::TriggerResults> const& triggerResults,
                                                     const std::vector<std::string>& triggerNames_,
                                                     edm::Handle<reco::TrackCollection>& trkCollection,
                                                     math::XYZPoint& leadPV,
                                                     edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
                                                     edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
                                                     edm::Handle<HBHERecHitCollection>& hbhe,
                                                     double ptL1,
                                                     double etaL1,
                                                     double phiL1);
  void setPtEtaPhi(std::vector<edm::Ref<l1extra::L1JetParticleCollection> >& objref,
                   double& ptL1,
                   double& etaL1,
                   double& phiL1);

  // ----------member data ---------------------------
  HLTConfigProvider hltConfig_;
  unsigned int nRun_, nAll_, nGood_, nRange_;
  spr::trackSelectionParameters selectionParameter_;
  const std::vector<std::string> trigNames_;
  const std::string theTrackQuality_, processName_;
  const double a_coneR_, a_mipR_;
  const double maxRestrictionP_, slopeRestrictionP_;
  const double pTrackMin_, eEcalMax_, eIsolate_;
  const double pTrackLow_, pTrackHigh_;
  const int preScale_;
  const edm::InputTag labelGenTrack_, labelRecVtx_, labelBS_;
  const edm::InputTag labelEB_, labelEE_, labelHBHE_, labelHltGT_;
  const edm::InputTag labelTriggerEvent_, labelTriggerResults_;
  const std::string labelIsoTk_;
  double a_charIsoR_;
  const MagneticField* bField;
  const CaloGeometry* geo;

  edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_hltGT_;
  edm::EDGetTokenT<trigger::TriggerEvent> tok_trigEvt_;
  edm::EDGetTokenT<edm::TriggerResults> tok_trigRes_;
  edm::EDGetTokenT<reco::TrackCollection> tok_genTrack_;
  edm::EDGetTokenT<reco::VertexCollection> tok_recVtx_;
  edm::EDGetTokenT<reco::BeamSpot> tok_bs_;
  edm::EDGetTokenT<EcalRecHitCollection> tok_EB_;
  edm::EDGetTokenT<EcalRecHitCollection> tok_EE_;
  edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
  edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
  edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_magField_;
};

AlCaIsoTracksProducer::AlCaIsoTracksProducer(edm::ParameterSet const& iConfig, const AlCaIsoTracks::Counters* counters)
    : nRun_(0),
      nAll_(0),
      nGood_(0),
      nRange_(0),
      trigNames_(iConfig.getParameter<std::vector<std::string> >("triggers")),
      theTrackQuality_(iConfig.getParameter<std::string>("trackQuality")),
      processName_(iConfig.getParameter<std::string>("processName")),
      a_coneR_(iConfig.getParameter<double>("coneRadius")),
      a_mipR_(iConfig.getParameter<double>("coneRadiusMIP")),
      maxRestrictionP_(iConfig.getParameter<double>("maxTrackP")),
      slopeRestrictionP_(iConfig.getParameter<double>("slopeTrackP")),
      pTrackMin_(iConfig.getParameter<double>("minimumTrackP")),
      eEcalMax_(iConfig.getParameter<double>("maximumEcalEnergy")),
      eIsolate_(iConfig.getParameter<double>("isolationEnergy")),
      pTrackLow_(iConfig.getParameter<double>("momentumRangeLow")),
      pTrackHigh_(iConfig.getParameter<double>("momentumRangeHigh")),
      preScale_(iConfig.getParameter<int>("preScaleFactor")),
      labelGenTrack_(iConfig.getParameter<edm::InputTag>("TrackLabel")),
      labelRecVtx_(iConfig.getParameter<edm::InputTag>("VertexLabel")),
      labelBS_(iConfig.getParameter<edm::InputTag>("BeamSpotLabel")),
      labelEB_(iConfig.getParameter<edm::InputTag>("EBRecHitLabel")),
      labelEE_(iConfig.getParameter<edm::InputTag>("EERecHitLabel")),
      labelHBHE_(iConfig.getParameter<edm::InputTag>("HBHERecHitLabel")),
      labelHltGT_(iConfig.getParameter<edm::InputTag>("L1GTSeedLabel")),
      labelTriggerEvent_(iConfig.getParameter<edm::InputTag>("TriggerEventLabel")),
      labelTriggerResults_(iConfig.getParameter<edm::InputTag>("TriggerResultLabel")),
      labelIsoTk_(iConfig.getParameter<std::string>("IsoTrackLabel")) {
  // Get the run parameters
  // Different isolation cuts are described in DN-2016/029
  // Tight cut uses 2 GeV; Loose cut uses 10 GeV
  // Eta dependent cut uses (maxRestrictionP_ * exp(|ieta|*log(2.5)/18))
  // with the factor for exponential slopeRestrictionP_ = log(2.5)/18
  // maxRestrictionP_ = 8 GeV as came from a study
  const double isolationRadius(28.9);
  selectionParameter_.minPt = iConfig.getParameter<double>("minTrackPt");
  ;
  selectionParameter_.minQuality = reco::TrackBase::qualityByName(theTrackQuality_);
  selectionParameter_.maxDxyPV = iConfig.getParameter<double>("maxDxyPV");
  selectionParameter_.maxDzPV = iConfig.getParameter<double>("maxDzPV");
  selectionParameter_.maxChi2 = iConfig.getParameter<double>("maxChi2");
  selectionParameter_.maxDpOverP = iConfig.getParameter<double>("maxDpOverP");
  selectionParameter_.minOuterHit = iConfig.getParameter<int>("minOuterHit");
  selectionParameter_.minLayerCrossed = iConfig.getParameter<int>("minLayerCrossed");
  selectionParameter_.maxInMiss = iConfig.getParameter<int>("maxInMiss");
  selectionParameter_.maxOutMiss = iConfig.getParameter<int>("maxOutMiss");
  a_charIsoR_ = a_coneR_ + isolationRadius;

  // define tokens for access
  tok_hltGT_ = consumes<trigger::TriggerFilterObjectWithRefs>(labelHltGT_);
  tok_trigEvt_ = consumes<trigger::TriggerEvent>(labelTriggerEvent_);
  tok_trigRes_ = consumes<edm::TriggerResults>(labelTriggerResults_);
  tok_genTrack_ = consumes<reco::TrackCollection>(labelGenTrack_);
  tok_recVtx_ = consumes<reco::VertexCollection>(labelRecVtx_);
  tok_bs_ = consumes<reco::BeamSpot>(labelBS_);
  tok_EB_ = consumes<EcalRecHitCollection>(labelEB_);
  tok_EE_ = consumes<EcalRecHitCollection>(labelEE_);
  tok_hbhe_ = consumes<HBHERecHitCollection>(labelHBHE_);
  // for event setup
  tok_geom_ = esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>();
  tok_magField_ = esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>();

  edm::LogVerbatim("HcalIsoTrack") << "Parameters read from config file \n"
                                   << "\t minPt " << selectionParameter_.minPt << "\t theTrackQuality "
                                   << theTrackQuality_ << "\t minQuality " << selectionParameter_.minQuality
                                   << "\t maxDxyPV " << selectionParameter_.maxDxyPV << "\t maxDzPV "
                                   << selectionParameter_.maxDzPV << "\t maxChi2 " << selectionParameter_.maxChi2
                                   << "\t maxDpOverP " << selectionParameter_.maxDpOverP << "\t minOuterHit "
                                   << selectionParameter_.minOuterHit << "\t minLayerCrossed "
                                   << selectionParameter_.minLayerCrossed << "\t maxInMiss "
                                   << selectionParameter_.maxInMiss << "\t maxOutMiss "
                                   << selectionParameter_.maxOutMiss << "\n"
                                   << "\t a_coneR " << a_coneR_ << "\t a_charIsoR " << a_charIsoR_ << "\t a_mipR "
                                   << a_mipR_ << "\t pTrackMin " << pTrackMin_ << "\t eEcalMax " << eEcalMax_
                                   << "\t maxRestrictionP_ " << maxRestrictionP_ << "\t slopeRestrictionP_ "
                                   << slopeRestrictionP_ << "\t eIsolate_ " << eIsolate_ << "\t Process "
                                   << processName_ << "\n"
                                   << "\t Precale factor " << preScale_ << "\t in momentum range " << pTrackLow_ << ":"
                                   << pTrackHigh_;
  for (unsigned int k = 0; k < trigNames_.size(); ++k)
    edm::LogVerbatim("HcalIsoTrack") << "Trigger[" << k << "] " << trigNames_[k];

  //create also IsolatedPixelTrackCandidateCollection which contains isolation info and reference to primary track
  produces<reco::HcalIsolatedTrackCandidateCollection>(labelIsoTk_);
  produces<reco::VertexCollection>(labelRecVtx_.label());
  produces<EcalRecHitCollection>(labelEB_.instance());
  produces<EcalRecHitCollection>(labelEE_.instance());
  produces<HBHERecHitCollection>(labelHBHE_.label());

  edm::LogVerbatim("HcalIsoTrack") << " Expected to produce the collections:\n"
                                   << "reco::HcalIsolatedTrackCandidateCollection "
                                   << " with label HcalIsolatedTrackCollection\n"
                                   << "reco::VertexCollection with label " << labelRecVtx_.label() << "\n"
                                   << "EcalRecHitCollection with label EcalRecHitsEB\n"
                                   << "EcalRecHitCollection with label EcalRecHitsEE\n"
                                   << "HBHERecHitCollection with label " << labelHBHE_.label();
}

AlCaIsoTracksProducer::~AlCaIsoTracksProducer() {}

void AlCaIsoTracksProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  // producer for  (HCAL isolated tracks)
  desc.add<edm::InputTag>("TrackLabel", edm::InputTag("generalTracks"));
  desc.add<edm::InputTag>("VertexLabel", edm::InputTag("offlinePrimaryVertices"));
  desc.add<edm::InputTag>("BeamSpotLabel", edm::InputTag("offlineBeamSpot"));
  desc.add<edm::InputTag>("EBRecHitLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
  desc.add<edm::InputTag>("EERecHitLabel", edm::InputTag("ecalRecHit", "EcalRecHitsEE"));
  desc.add<edm::InputTag>("HBHERecHitLabel", edm::InputTag("hbhereco"));
  desc.add<edm::InputTag>("L1GTSeedLabel", edm::InputTag("hltL1sV0SingleJet60"));
  desc.add<edm::InputTag>("TriggerEventLabel", edm::InputTag("hltTriggerSummaryAOD", "", "HLT"));
  desc.add<edm::InputTag>("TriggerResultLabel", edm::InputTag("TriggerResults", "", "HLT"));
  desc.add<std::string>("IsoTrackLabel", "HcalIsolatedTrackCollection");
  std::vector<std::string> triggers = {"HLT_IsoTrackHB", "HLT_IsoTrackHE"};
  desc.add<std::vector<std::string> >("triggers", triggers);
  desc.add<std::string>("processName", "HLT");
  // following 10 parameters are parameters to select good tracks
  desc.add<std::string>("trackQuality", "highPurity");
  desc.add<double>("minTrackPt", 1.0);
  desc.add<double>("maxDxyPV", 10.0);
  desc.add<double>("maxDzPV", 100.0);
  desc.add<double>("maxChi2", 5.0);
  desc.add<double>("maxDpOverP", 0.1);
  desc.add<int>("minOuterHit", 4);
  desc.add<int>("minLayerCrossed", 8);
  desc.add<int>("maxInMiss", 2);
  desc.add<int>("maxOutMiss", 2);
  // Minimum momentum of selected isolated track and signal zone
  desc.add<double>("coneRadius", 34.98);
  desc.add<double>("minimumTrackP", 20.0);
  // signal zone in ECAL and MIP energy cutoff
  desc.add<double>("coneRadiusMIP", 14.0);
  desc.add<double>("maximumEcalEnergy", 2.0);
  // following 3 parameters are for isolation cuts and described in the code
  desc.add<double>("maxTrackP", 8.0);
  desc.add<double>("slopeTrackP", 0.05090504066);
  desc.add<double>("isolationEnergy", 10.0);
  // Prescale events only containing isolated tracks in the range
  desc.add<double>("momentumRangeLow", 20.0);
  desc.add<double>("momentumRangeHigh", 40.0);
  desc.add<int>("preScaleFactor", 10);
  descriptions.add("alcaisotrk", desc);
}

void AlCaIsoTracksProducer::produce(edm::Event& iEvent, edm::EventSetup const& iSetup) {
  nAll_++;
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalIsoTrack") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event() << " Luminosity "
                                   << iEvent.luminosityBlock() << " Bunch " << iEvent.bunchCrossing();
#endif
  bool valid(true);
  //Step1: Get all the relevant containers
  trigger::TriggerEvent triggerEvent;
  auto const& triggerEventHandle = iEvent.getHandle(tok_trigEvt_);
  if (!triggerEventHandle.isValid()) {
    edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelTriggerEvent_;
    valid = false;
  }
  auto const& triggerResults = iEvent.getHandle(tok_trigRes_);
  if (!triggerResults.isValid()) {
    edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelTriggerResults_;
    valid = false;
  }

  auto trkCollection = iEvent.getHandle(tok_genTrack_);
  if (!trkCollection.isValid()) {
    edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
    valid = false;
  }

  auto const& recVtxs = iEvent.getHandle(tok_recVtx_);
  if (!recVtxs.isValid()) {
    edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelGenTrack_;
    valid = false;
  }

  auto const& beamSpotH = iEvent.getHandle(tok_bs_);
  math::XYZPoint leadPV(0, 0, 0);
  if (valid) {
    if (!recVtxs->empty() && !((*recVtxs)[0].isFake())) {
      leadPV = math::XYZPoint((*recVtxs)[0].x(), (*recVtxs)[0].y(), (*recVtxs)[0].z());
    } else if (beamSpotH.isValid()) {
      leadPV = beamSpotH->position();
    }
  }
#ifdef EDM_ML_DEBUG
  edm::LogVerbatim("HcalIsoTrack") << "Primary Vertex " << leadPV;
#endif

  auto barrelRecHitsHandle = iEvent.getHandle(tok_EB_);
  if (!barrelRecHitsHandle.isValid()) {
    edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEB_;
    valid = false;
  }
  auto endcapRecHitsHandle = iEvent.getHandle(tok_EE_);
  if (!endcapRecHitsHandle.isValid()) {
    edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelEE_;
    valid = false;
  }
  auto hbhe = iEvent.getHandle(tok_hbhe_);
  if (!hbhe.isValid()) {
    edm::LogWarning("HcalIsoTrack") << "Cannot access the collection " << labelHBHE_;
    valid = false;
  }

  //Get L1 trigger object
  double ptL1(0), etaL1(0), phiL1(0);
  auto const& l1trigobj = iEvent.getHandle(tok_hltGT_);

  if (l1trigobj.isValid()) {
    std::vector<edm::Ref<l1extra::L1JetParticleCollection> > l1tauobjref;
    l1trigobj->getObjects(trigger::TriggerL1TauJet, l1tauobjref);
    setPtEtaPhi(l1tauobjref, ptL1, etaL1, phiL1);

    std::vector<edm::Ref<l1extra::L1JetParticleCollection> > l1jetobjref;
    l1trigobj->getObjects(trigger::TriggerL1CenJet, l1jetobjref);
    setPtEtaPhi(l1jetobjref, ptL1, etaL1, phiL1);

    std::vector<edm::Ref<l1extra::L1JetParticleCollection> > l1forjetobjref;
    l1trigobj->getObjects(trigger::TriggerL1ForJet, l1forjetobjref);
    setPtEtaPhi(l1forjetobjref, ptL1, etaL1, phiL1);
  } else {
    valid = false;
  }

  auto outputHcalIsoTrackColl = std::make_unique<reco::HcalIsolatedTrackCandidateCollection>();
  auto outputVColl = std::make_unique<reco::VertexCollection>();
  auto outputEBColl = std::make_unique<EBRecHitCollection>();
  auto outputEEColl = std::make_unique<EERecHitCollection>();
  auto outputHBHEColl = std::make_unique<HBHERecHitCollection>();

  //For valid HLT record
  if (!valid) {
    edm::LogWarning("HcalIsoTrack") << "Error! Can't get some of the products";
  } else {
    trigger::TriggerEvent triggerEvent = *(triggerEventHandle.product());
    if (triggerResults.isValid()) {
      const edm::TriggerNames& triggerNames = iEvent.triggerNames(*triggerResults);
      const std::vector<std::string>& triggerNames_ = triggerNames.triggerNames();
      reco::HcalIsolatedTrackCandidateCollection* isotk = select(triggerResults,
                                                                 triggerNames_,
                                                                 trkCollection,
                                                                 leadPV,
                                                                 barrelRecHitsHandle,
                                                                 endcapRecHitsHandle,
                                                                 hbhe,
                                                                 ptL1,
                                                                 etaL1,
                                                                 phiL1);
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HcalIsoTrack") << "AlCaIsoTracksProducer::select returns " << isotk->size()
                                       << " isolated tracks";
#endif

      if (!isotk->empty()) {
        int ntrin(0), ntrout(0);
        for (reco::HcalIsolatedTrackCandidateCollection::const_iterator itr = isotk->begin(); itr != isotk->end();
             ++itr) {
          if (itr->p() > pTrackLow_ && itr->p() < pTrackHigh_)
            ntrin++;
          else
            ntrout++;
        }
        bool selectEvent = ntrout > 0;
        if (!selectEvent && ntrin > 0) {
          ++nRange_;
          if (preScale_ <= 1)
            selectEvent = true;
          else if (nRange_ % preScale_ == 1)
            selectEvent = true;
        }
        if (selectEvent) {
          for (reco::HcalIsolatedTrackCandidateCollection::const_iterator itr = isotk->begin(); itr != isotk->end();
               ++itr)
            outputHcalIsoTrackColl->push_back(*itr);

          for (reco::VertexCollection::const_iterator vtx = recVtxs->begin(); vtx != recVtxs->end(); ++vtx)
            outputVColl->push_back(*vtx);

          for (edm::SortedCollection<EcalRecHit>::const_iterator ehit = barrelRecHitsHandle->begin();
               ehit != barrelRecHitsHandle->end();
               ++ehit)
            outputEBColl->push_back(*ehit);

          for (edm::SortedCollection<EcalRecHit>::const_iterator ehit = endcapRecHitsHandle->begin();
               ehit != endcapRecHitsHandle->end();
               ++ehit)
            outputEEColl->push_back(*ehit);

          for (std::vector<HBHERecHit>::const_iterator hhit = hbhe->begin(); hhit != hbhe->end(); ++hhit)
            outputHBHEColl->push_back(*hhit);
          ++nGood_;
        }
      }
    }
  }
  iEvent.put(std::move(outputHcalIsoTrackColl), labelIsoTk_);
  iEvent.put(std::move(outputVColl), labelRecVtx_.label());
  iEvent.put(std::move(outputEBColl), labelEB_.instance());
  iEvent.put(std::move(outputEEColl), labelEE_.instance());
  iEvent.put(std::move(outputHBHEColl), labelHBHE_.label());
}

void AlCaIsoTracksProducer::endStream() {
  globalCache()->nAll_ += nAll_;
  globalCache()->nGood_ += nGood_;
  globalCache()->nRange_ += nRange_;
}

void AlCaIsoTracksProducer::globalEndJob(const AlCaIsoTracks::Counters* count) {
  edm::LogVerbatim("HcalIsoTrack") << "Finds " << count->nGood_ << " good tracks in " << count->nAll_ << " events and "
                                   << count->nRange_ << " events in the momentum raange";
}

void AlCaIsoTracksProducer::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
  bool changed(false);
  edm::LogVerbatim("HcalIsoTrack") << "Run[" << nRun_ << "] " << iRun.run() << " hltconfig.init "
                                   << hltConfig_.init(iRun, iSetup, processName_, changed);

  bField = &(iSetup.getData(tok_magField_));
  geo = &(iSetup.getData(tok_geom_));
}

void AlCaIsoTracksProducer::endRun(edm::Run const& iRun, edm::EventSetup const&) {
  edm::LogVerbatim("HcalIsoTrack") << "endRun [" << nRun_ << "] " << iRun.run();
  ++nRun_;
}

reco::HcalIsolatedTrackCandidateCollection* AlCaIsoTracksProducer::select(
    edm::Handle<edm::TriggerResults> const& triggerResults,
    const std::vector<std::string>& triggerNames_,
    edm::Handle<reco::TrackCollection>& trkCollection,
    math::XYZPoint& leadPV,
    edm::Handle<EcalRecHitCollection>& barrelRecHitsHandle,
    edm::Handle<EcalRecHitCollection>& endcapRecHitsHandle,
    edm::Handle<HBHERecHitCollection>& hbhe,
    double ptL1,
    double etaL1,
    double phiL1) {
  reco::HcalIsolatedTrackCandidateCollection* trackCollection = new reco::HcalIsolatedTrackCandidateCollection;
  bool ok(false);

  // Find a good HLT trigger
  for (unsigned int iHLT = 0; iHLT < triggerResults->size(); iHLT++) {
    int hlt = triggerResults->accept(iHLT);
    for (unsigned int i = 0; i < trigNames_.size(); ++i) {
      if (triggerNames_[iHLT].find(trigNames_[i]) != std::string::npos) {
        if (hlt > 0) {
          ok = true;
        }
        edm::LogVerbatim("HcalIsoTrack") << "The trigger we are looking for " << triggerNames_[iHLT] << " Flag " << hlt
                                         << ":" << ok;
      }
    }
  }

  //Propagate tracks to calorimeter surface)
  std::vector<spr::propagatedTrackDirection> trkCaloDirections;
  spr::propagateCALO(trkCollection, geo, bField, theTrackQuality_, trkCaloDirections, false);

  std::vector<spr::propagatedTrackDirection>::const_iterator trkDetItr;
  unsigned int nTracks(0), nselTracks(0);
  for (trkDetItr = trkCaloDirections.begin(), nTracks = 0; trkDetItr != trkCaloDirections.end();
       trkDetItr++, nTracks++) {
    const reco::Track* pTrack = &(*(trkDetItr->trkItr));
    math::XYZTLorentzVector v4(pTrack->px(), pTrack->py(), pTrack->pz(), pTrack->p());
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalIsoTrack") << "This track : " << nTracks << " (pt|eta|phi|p) :" << pTrack->pt() << "|"
                                     << pTrack->eta() << "|" << pTrack->phi() << "|" << pTrack->p();
#endif
    //Selection of good track
    bool qltyFlag = spr::goodTrack(pTrack, leadPV, selectionParameter_, false);
#ifdef EDM_ML_DEBUG
    edm::LogVerbatim("HcalIsoTrack") << "qltyFlag|okECAL|okHCAL : " << qltyFlag << "|" << trkDetItr->okECAL << "|"
                                     << trkDetItr->okHCAL;
#endif
    if (qltyFlag && trkDetItr->okECAL && trkDetItr->okHCAL) {
      double t_p = pTrack->p();
      nselTracks++;
      int nRH_eMipDR(0), nNearTRKs(0);
      double eMipDR = spr::eCone_ecal(geo,
                                      barrelRecHitsHandle,
                                      endcapRecHitsHandle,
                                      trkDetItr->pointHCAL,
                                      trkDetItr->pointECAL,
                                      a_mipR_,
                                      trkDetItr->directionECAL,
                                      nRH_eMipDR);
      double hmaxNearP = spr::chargeIsolationCone(nTracks, trkCaloDirections, a_charIsoR_, nNearTRKs, false);
      HcalDetId detId = (HcalDetId)(trkDetItr->detIdHCAL);
      int ieta = detId.ietaAbs();
      double eIsolation = (maxRestrictionP_ * exp(slopeRestrictionP_ * ((double)(ieta))));
      if (eIsolation < eIsolate_)
        eIsolation = eIsolate_;
#ifdef EDM_ML_DEBUG
      edm::LogVerbatim("HcalIsoTrack") << "This track : " << nTracks << " (pt|eta|phi|p) :" << pTrack->pt() << "|"
                                       << pTrack->eta() << "|" << pTrack->phi() << "|" << t_p << " e_MIP " << eMipDR
                                       << " Chg Isolation " << hmaxNearP << ":" << eIsolation;
#endif
      if (t_p > pTrackMin_ && eMipDR < eEcalMax_ && hmaxNearP < eIsolation) {
        reco::HcalIsolatedTrackCandidate newCandidate(v4);
        newCandidate.SetMaxP(hmaxNearP);
        newCandidate.SetEnergyEcal(eMipDR);
        newCandidate.setL1(ptL1, etaL1, phiL1);
        newCandidate.SetEtaPhiEcal((trkDetItr->pointECAL).eta(), (trkDetItr->pointECAL).phi());
        HcalDetId detId = HcalDetId(trkDetItr->detIdHCAL);
        newCandidate.SetEtaPhiHcal(
            (trkDetItr->pointHCAL).eta(), (trkDetItr->pointHCAL).phi(), detId.ieta(), detId.iphi());
        int indx(0);
        for (reco::TrackCollection::const_iterator trkItr1 = trkCollection->begin(); trkItr1 != trkCollection->end();
             ++trkItr1, ++indx) {
          const reco::Track* pTrack1 = &(*trkItr1);
          if (pTrack1 == pTrack) {
            reco::TrackRef tRef = reco::TrackRef(trkCollection, indx);
            newCandidate.setTrack(tRef);
            break;
          }
        }
        trackCollection->push_back(newCandidate);
      }
    }
  }
  return trackCollection;
}

void AlCaIsoTracksProducer::setPtEtaPhi(std::vector<edm::Ref<l1extra::L1JetParticleCollection> >& objref,
                                        double& ptL1,
                                        double& etaL1,
                                        double& phiL1) {
  for (unsigned int p = 0; p < objref.size(); p++) {
    if (objref[p]->pt() > ptL1) {
      ptL1 = objref[p]->pt();
      phiL1 = objref[p]->phi();
      etaL1 = objref[p]->eta();
    }
  }
}

#include "FWCore/Framework/interface/MakerMacros.h"

DEFINE_FWK_MODULE(AlCaIsoTracksProducer);