GoodSeedProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    PFTracking
// Class:      GoodSeedProducer
//
// Original Author:  Michele Pioppi
// March 2010. F. Beaudette. Produce PreId information

#include "RecoParticleFlow/PFTracking/plugins/GoodSeedProducer.h"
#include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
#include "RecoParticleFlow/PFClusterTools/interface/PFResolutionMap.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "DataFormats/EgammaReco/interface/ElectronSeed.h"
#include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
#include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
#include "DataFormats/ParticleFlowReco/interface/PreId.h"
#include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
#include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "CommonTools/MVAUtils/interface/GBRForestTools.h"

#include "DataFormats/Math/interface/deltaR.h"

#include <fstream>
#include <string>
#include "TMath.h"
#include "Math/VectorUtil.h"

using namespace edm;
using namespace std;
using namespace reco;

GoodSeedProducer::GoodSeedProducer(const ParameterSet& iConfig, const goodseedhelpers::HeavyObjectCache*)
    : pfTransformer_(nullptr), conf_(iConfig), resMapEtaECAL_(nullptr), resMapPhiECAL_(nullptr) {
  LogInfo("GoodSeedProducer") << "Electron PreIdentification started  ";

  //now do what ever initialization is needed
  std::vector<edm::InputTag> tags = iConfig.getParameter<vector<InputTag>>("TkColList");
  for (unsigned int i = 0; i < tags.size(); ++i) {
    trajContainers_.push_back(consumes<vector<Trajectory>>(tags[i]));
    tracksContainers_.push_back(consumes<reco::TrackCollection>(tags[i]));
  }

  minPt_ = iConfig.getParameter<double>("MinPt");
  maxPt_ = iConfig.getParameter<double>("MaxPt");
  maxEta_ = iConfig.getParameter<double>("MaxEta");

  HcalIsolWindow_ = iConfig.getParameter<double>("HcalWindow");
  EcalStripSumE_minClusEnergy_ = iConfig.getParameter<double>("EcalStripSumE_minClusEnergy");
  EcalStripSumE_deltaEta_ = iConfig.getParameter<double>("EcalStripSumE_deltaEta");
  EcalStripSumE_deltaPhiOverQ_minValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_minValue");
  EcalStripSumE_deltaPhiOverQ_maxValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_maxValue");
  minEoverP_ = iConfig.getParameter<double>("EOverPLead_minValue");
  maxHoverP_ = iConfig.getParameter<double>("HOverPLead_maxValue");

  pfCLusTagECLabel_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFEcalClusterLabel"));

  pfCLusTagHCLabel_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFHcalClusterLabel"));

  pfCLusTagPSLabel_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFPSClusterLabel"));

  preidgsf_ = iConfig.getParameter<string>("PreGsfLabel");
  preidckf_ = iConfig.getParameter<string>("PreCkfLabel");
  preidname_ = iConfig.getParameter<string>("PreIdLabel");

  fitterName_ = iConfig.getParameter<string>("Fitter");
  smootherName_ = iConfig.getParameter<string>("Smoother");

  nHitsInSeed_ = iConfig.getParameter<int>("NHitsInSeed");

  clusThreshold_ = iConfig.getParameter<double>("ClusterThreshold");

  minEp_ = iConfig.getParameter<double>("MinEOverP");
  maxEp_ = iConfig.getParameter<double>("MaxEOverP");

  //collection to produce
  produceCkfseed_ = iConfig.getUntrackedParameter<bool>("ProduceCkfSeed", false);

  // to disable the electron part (for HI collisions for examples)
  disablePreId_ = iConfig.getUntrackedParameter<bool>("DisablePreId", false);

  producePreId_ = iConfig.getUntrackedParameter<bool>("ProducePreId", true);
  //  if no electron, cannot produce the preid
  if (disablePreId_)
    producePreId_ = false;
  PtThresholdSavePredId_ = iConfig.getUntrackedParameter<double>("PtThresholdSavePreId", 1.);

  LogDebug("GoodSeedProducer") << "Seeds for GSF will be produced ";

  // no disablePreId_ switch here. The collection will be empty if it is true
  produces<ElectronSeedCollection>(preidgsf_);

  if (produceCkfseed_) {
    LogDebug("GoodSeedProducer") << "Seeds for CKF will be produced ";
    produces<TrajectorySeedCollection>(preidckf_);
  }

  if (producePreId_) {
    LogDebug("GoodSeedProducer") << "PreId debugging information will be produced ";

    produces<PreIdCollection>(preidname_);
    if (tracksContainers_.size() == 1)  // do not make a value map if more than one input track collection
      produces<edm::ValueMap<reco::PreIdRef>>(preidname_);
  }

  useQuality_ = iConfig.getParameter<bool>("UseQuality");
  trackQuality_ = TrackBase::qualityByName(iConfig.getParameter<std::string>("TrackQuality"));

  useTmva_ = iConfig.getUntrackedParameter<bool>("UseTMVA", false);

  Min_dr_ = iConfig.getParameter<double>("Min_dr");

  trackerRecHitBuilderName_ = iConfig.getParameter<std::string>("TTRHBuilder");
}

//
// member functions
//

// ------------ method called to produce the data  ------------
void GoodSeedProducer::produce(Event& iEvent, const EventSetup& iSetup) {
  LogDebug("GoodSeedProducer") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run();
  //Create empty output collections
  auto output_preid = std::make_unique<ElectronSeedCollection>();
  auto output_nopre = std::make_unique<TrajectorySeedCollection>();
  auto output_preidinfo = std::make_unique<PreIdCollection>();
  auto preIdMap_p = std::make_unique<edm::ValueMap<reco::PreIdRef>>();
  edm::ValueMap<reco::PreIdRef>::Filler mapFiller(*preIdMap_p);

  //Tracking Tools
  if (!disablePreId_) {
    edm::ESHandle<TrajectoryFitter> aFitter;
    edm::ESHandle<TrajectorySmoother> aSmoother;
    iSetup.get<TrajectoryFitter::Record>().get(fitterName_, aFitter);
    iSetup.get<TrajectoryFitter::Record>().get(smootherName_, aSmoother);
    smoother_.reset(aSmoother->clone());
    fitter_ = aFitter->clone();
    edm::ESHandle<TransientTrackingRecHitBuilder> theTrackerRecHitBuilder;
    iSetup.get<TransientRecHitRecord>().get(trackerRecHitBuilderName_, theTrackerRecHitBuilder);
    hitCloner = static_cast<TkTransientTrackingRecHitBuilder const*>(theTrackerRecHitBuilder.product())->cloner();
    fitter_->setHitCloner(&hitCloner);
    smoother_->setHitCloner(&hitCloner);
  }

  // clear temporary maps
  refMap_.clear();

  //Magnetic Field
  ESHandle<MagneticField> magneticField;
  iSetup.get<IdealMagneticFieldRecord>().get(magneticField);

  //Handle input collections
  //ECAL clusters
  Handle<PFClusterCollection> theECPfClustCollection;
  iEvent.getByToken(pfCLusTagECLabel_, theECPfClustCollection);

  vector<PFCluster const*> basClus;
  for (auto const& klus : *theECPfClustCollection.product()) {
    if (klus.correctedEnergy() > clusThreshold_)
      basClus.push_back(&klus);
  }

  //HCAL clusters
  Handle<PFClusterCollection> theHCPfClustCollection;
  iEvent.getByToken(pfCLusTagHCLabel_, theHCPfClustCollection);

  //PS clusters
  Handle<PFClusterCollection> thePSPfClustCollection;
  iEvent.getByToken(pfCLusTagPSLabel_, thePSPfClustCollection);

  //Vector of track collections
  for (unsigned int istr = 0; istr < tracksContainers_.size(); ++istr) {
    //Track collection
    Handle<TrackCollection> tkRefCollection;
    iEvent.getByToken(tracksContainers_[istr], tkRefCollection);
    const TrackCollection& Tk = *(tkRefCollection.product());

    LogDebug("GoodSeedProducer") << "Number of tracks in collection "
                                 << "tracksContainers_[" << istr << "] to be analyzed " << Tk.size();

    //loop over the track collection
    for (unsigned int i = 0; i < Tk.size(); ++i) {
      if (useQuality_ && (!(Tk[i].quality(trackQuality_))))
        continue;

      reco::PreId myPreId;
      bool GoodPreId = false;

      TrackRef trackRef(tkRefCollection, i);
      math::XYZVectorF tkmom(Tk[i].momentum());
      auto tketa = tkmom.eta();
      auto tkpt = std::sqrt(tkmom.perp2());
      auto const& Seed = (*trackRef->seedRef());
      if (!disablePreId_) {
        int ipteta = getBin(Tk[i].eta(), Tk[i].pt());
        int ibin = ipteta * 9;

        // FIXME the original code was buggy should be outerMomentum...
        float oPTOB = 1.f / std::sqrt(Tk[i].innerMomentum().mag2());
        //  float chikfred=Tk[i].normalizedChi2();
        float nchi = Tk[i].normalizedChi2();

        int nhitpi = Tk[i].found();
        float EP = 0;

        // set track info
        myPreId.setTrack(trackRef);
        //CLUSTERS - TRACK matching

        auto pfmass = 0.0005;
        auto pfoutenergy = sqrt((pfmass * pfmass) + Tk[i].outerMomentum().Mag2());

        XYZTLorentzVector mom = XYZTLorentzVector(
            Tk[i].outerMomentum().x(), Tk[i].outerMomentum().y(), Tk[i].outerMomentum().z(), pfoutenergy);
        XYZTLorentzVector pos =
            XYZTLorentzVector(Tk[i].outerPosition().x(), Tk[i].outerPosition().y(), Tk[i].outerPosition().z(), 0.);

        BaseParticlePropagator theOutParticle(RawParticle(mom, pos, Tk[i].charge()), 0, 0, B_.z());

        theOutParticle.propagateToEcalEntrance(false);

        float toteta = 1000.f;
        float totphi = 1000.f;
        float dr = 1000.f;
        float EE = 0.f;
        float feta = 0.f;
        GlobalPoint ElecTrkEcalPos(0, 0, 0);

        PFClusterRef clusterRef;
        math::XYZPoint meanShowerSaved;
        if (theOutParticle.getSuccess() != 0) {
          ElecTrkEcalPos = GlobalPoint(theOutParticle.particle().vertex().x(),
                                       theOutParticle.particle().vertex().y(),
                                       theOutParticle.particle().vertex().z());

          constexpr float psLim = 2.50746495928f;  // std::sinh(1.65f);
          bool isBelowPS = (ElecTrkEcalPos.z() * ElecTrkEcalPos.z()) > (psLim * psLim) * ElecTrkEcalPos.perp2();
          // bool isBelowPS=(std::abs(ElecTrkEcalPos.eta())>1.65f);

          unsigned clusCounter = 0;
          float max_ee = 0;
          for (auto aClus : basClus) {
            float tmp_ep = float(aClus->correctedEnergy()) * oPTOB;
            if ((tmp_ep < minEp_) | (tmp_ep > maxEp_)) {
              ++clusCounter;
              continue;
            }

            double ecalShowerDepth = PFCluster::getDepthCorrection(aClus->correctedEnergy(), isBelowPS, false);
            auto mom = theOutParticle.particle().momentum().Vect();
            auto meanShower = ElecTrkEcalPos + GlobalVector(mom.x(), mom.y(), mom.z()).unit() * ecalShowerDepth;

            float etarec = meanShower.eta();
            float phirec = meanShower.phi();

            float tmp_phi = std::abs(aClus->positionREP().phi() - phirec);
            if (tmp_phi > float(TMath::Pi()))
              tmp_phi -= float(TMath::TwoPi());

            float tmp_dr = std::sqrt(std::pow(tmp_phi, 2.f) + std::pow(aClus->positionREP().eta() - etarec, 2.f));

            if (tmp_dr < dr) {
              dr = tmp_dr;
              if (dr < Min_dr_) {  // find the most closest and energetic ECAL cluster
                if (aClus->correctedEnergy() > max_ee) {
                  toteta = aClus->positionREP().eta() - etarec;
                  totphi = tmp_phi;
                  EP = tmp_ep;
                  EE = aClus->correctedEnergy();
                  feta = aClus->positionREP().eta();
                  clusterRef = PFClusterRef(theECPfClustCollection, clusCounter);
                  meanShowerSaved = meanShower;
                }
              }
            }
            ++clusCounter;
          }
        }
        float trk_ecalDeta_ = fabs(toteta);
        float trk_ecalDphi_ = fabs(totphi);

        //Resolution maps
        auto ecaletares = resMapEtaECAL_->GetBinContent(resMapEtaECAL_->FindBin(feta, EE));
        auto ecalphires = resMapPhiECAL_->GetBinContent(resMapPhiECAL_->FindBin(feta, EE));

        //geomatrical compatibility
        float chieta = (toteta != 1000.f) ? toteta / ecaletares : toteta;
        float chiphi = (totphi != 1000.f) ? totphi / ecalphires : totphi;
        float chichi = sqrt(chieta * chieta + chiphi * chiphi);

        //Matching criteria
        float eta_cut = thr[ibin + 0];
        float phi_cut = thr[ibin + 1];
        float ep_cutmin = thr[ibin + 2];
        bool GoodMatching =
            ((trk_ecalDeta_ < eta_cut) && (trk_ecalDphi_ < phi_cut) && (EP > ep_cutmin) && (nhitpi > 10));

        bool EcalMatching = GoodMatching;

        if (tkpt > maxPt_)
          GoodMatching = true;
        if (tkpt < minPt_)
          GoodMatching = false;

        math::XYZPoint myPoint(ElecTrkEcalPos.x(), ElecTrkEcalPos.y(), ElecTrkEcalPos.z());
        myPreId.setECALMatchingProperties(
            clusterRef, myPoint, meanShowerSaved, std::abs(toteta), std::abs(totphi), chieta, chiphi, chichi, EP);
        myPreId.setECALMatching(EcalMatching);

        bool GoodRange = ((std::abs(tketa) < maxEta_) & (tkpt > minPt_));
        //KF FILTERING FOR UNMATCHED EVENTS
        int hit1max = int(thr[ibin + 3]);
        float chiredmin = thr[ibin + 4];
        bool GoodKFFiltering = ((nchi > chiredmin) | (nhitpi < hit1max));

        myPreId.setTrackFiltering(GoodKFFiltering);

        bool GoodTkId = false;

        if ((!GoodMatching) && (GoodKFFiltering) && (GoodRange)) {
          chired = 1000;
          chiRatio = 1000;
          dpt = 0;
          nhit = nhitpi;
          chikfred = nchi;
          trk_ecalDeta = trk_ecalDeta_;
          trk_ecalDphi = trk_ecalDphi_;

          Trajectory::ConstRecHitContainer tmp;
          for (auto const& hit : Tk[i].recHits())
            tmp.push_back(hit->cloneSH());
          auto const& theTrack = Tk[i];
          GlobalVector gv(theTrack.innerMomentum().x(), theTrack.innerMomentum().y(), theTrack.innerMomentum().z());
          GlobalPoint gp(theTrack.innerPosition().x(), theTrack.innerPosition().y(), theTrack.innerPosition().z());
          GlobalTrajectoryParameters gtps(gp, gv, theTrack.charge(), &*magneticField);
          TrajectoryStateOnSurface tsos(gtps, theTrack.innerStateCovariance(), *tmp[0]->surface());
          Trajectory&& FitTjs = fitter_->fitOne(Seed, tmp, tsos);

          if (FitTjs.isValid()) {
            Trajectory&& SmooTjs = smoother_->trajectory(FitTjs);
            if (SmooTjs.isValid()) {
              //Track refitted with electron hypothesis

              float pt_out = SmooTjs.firstMeasurement().updatedState().globalMomentum().perp();
              float pt_in = SmooTjs.lastMeasurement().updatedState().globalMomentum().perp();
              dpt = (pt_in > 0) ? fabs(pt_out - pt_in) / pt_in : 0.;
              // the following is simply the number of degrees of freedom
              chiRatio = SmooTjs.chiSquared() / Tk[i].chi2();
              chired = chiRatio * chikfred;
            }
          }

          //TMVA Analysis
          if (useTmva_) {
            eta = tketa;
            pt = tkpt;
            eP = EP;
            float vars[10] = {nhit, chikfred, dpt, eP, chiRatio, chired, trk_ecalDeta, trk_ecalDphi, pt, eta};

            float Ytmva = globalCache()->gbr[ipteta]->GetClassifier(vars);

            float BDTcut = thr[ibin + 5];
            if (Ytmva > BDTcut)
              GoodTkId = true;
            myPreId.setMVA(GoodTkId, Ytmva);
            myPreId.setTrackProperties(chired, chiRatio, dpt);
          } else {
            float chiratiocut = thr[ibin + 6];
            float gschicut = thr[ibin + 7];
            float gsptmin = thr[ibin + 8];

            GoodTkId = ((dpt > gsptmin) & (chired < gschicut) & (chiRatio < chiratiocut));
          }
        }

        GoodPreId = GoodTkId | GoodMatching;

        myPreId.setFinalDecision(GoodPreId);

#ifdef EDM_ML_DEBUG
        if (GoodPreId)
          LogDebug("GoodSeedProducer") << "Track (pt= " << Tk[i].pt() << "GeV/c, eta= " << Tk[i].eta()
                                       << ") preidentified for agreement between  track and ECAL cluster";
        if (GoodPreId && (!GoodMatching))
          LogDebug("GoodSeedProducer") << "Track (pt= " << Tk[i].pt() << "GeV/c, eta= " << Tk[i].eta()
                                       << ") preidentified only for track properties";
#endif

      }  // end of !disablePreId_

      if (GoodPreId) {
        //NEW SEED with n hits
        ElectronSeed NewSeed(Seed);
        NewSeed.setCtfTrack(trackRef);
        output_preid->push_back(NewSeed);
      } else {
        if (produceCkfseed_) {
          output_nopre->push_back(Seed);
        }
      }
      if (producePreId_ && myPreId.pt() > PtThresholdSavePredId_) {
        // save the index of the PreId object as to be able to create a Ref later
        refMap_[trackRef] = output_preidinfo->size();
        output_preidinfo->push_back(myPreId);
      }
    }  //end loop on track collection
  }    //end loop on the vector of track collections

  // no disablePreId_ switch, it is simpler to have an empty collection rather than no collection
  iEvent.put(std::move(output_preid), preidgsf_);
  if (produceCkfseed_)
    iEvent.put(std::move(output_nopre), preidckf_);
  if (producePreId_) {
    const edm::OrphanHandle<reco::PreIdCollection> preIdRefProd = iEvent.put(std::move(output_preidinfo), preidname_);
    // now make the Value Map, but only if one input collection
    if (tracksContainers_.size() == 1) {
      Handle<TrackCollection> tkRefCollection;
      iEvent.getByToken(tracksContainers_[0], tkRefCollection);
      fillPreIdRefValueMap(tkRefCollection, preIdRefProd, mapFiller);
      mapFiller.fill();
      iEvent.put(std::move(preIdMap_p), preidname_);
    }
  }
  // clear temporary maps
  refMap_.clear();
}

// intialize the cross-thread cache to hold gbr trees and resolution maps
namespace goodseedhelpers {
  HeavyObjectCache::HeavyObjectCache(const edm::ParameterSet& conf) {
    // mvas
    const bool useTmva = conf.getUntrackedParameter<bool>("UseTMVA", false);

    if (useTmva) {
      std::array<edm::FileInPath, kMaxWeights> weights = {{edm::FileInPath(conf.getParameter<string>("Weights1")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights2")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights3")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights4")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights5")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights6")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights7")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights8")),
                                                           edm::FileInPath(conf.getParameter<string>("Weights9"))}};

      for (UInt_t j = 0; j < gbr.size(); ++j) {
        gbr[j] = createGBRForest(weights[j]);
      }
    }
  }
}  // namespace goodseedhelpers

// ------------ method called once each job just before starting event loop  ------------
void GoodSeedProducer::beginRun(const edm::Run& run, const EventSetup& es) {
  //Magnetic Field
  ESHandle<MagneticField> magneticField;
  es.get<IdealMagneticFieldRecord>().get(magneticField);
  B_ = magneticField->inTesla(GlobalPoint(0, 0, 0));

  pfTransformer_.reset(new PFTrackTransformer(B_));
  pfTransformer_->OnlyProp();

  //Resolution maps
  FileInPath ecalEtaMap(conf_.getParameter<string>("EtaMap"));
  FileInPath ecalPhiMap(conf_.getParameter<string>("PhiMap"));
  resMapEtaECAL_.reset(new PFResolutionMap("ECAL_eta", ecalEtaMap.fullPath().c_str()));
  resMapPhiECAL_.reset(new PFResolutionMap("ECAL_phi", ecalPhiMap.fullPath().c_str()));

  //read threshold
  FileInPath parFile(conf_.getParameter<string>("ThresholdFile"));
  ifstream ifs(parFile.fullPath().c_str());
  for (int iy = 0; iy < 81; ++iy)
    ifs >> thr[iy];
}

int GoodSeedProducer::getBin(float eta, float pt) {
  int ie = 0;
  int ip = 0;
  if (fabs(eta) < 0.8)
    ie = 0;
  else {
    if (fabs(eta) < 1.479)
      ie = 1;
    else
      ie = 2;
  }
  if (pt < 6)
    ip = 0;
  else {
    if (pt < 12)
      ip = 1;
    else
      ip = 2;
  }
  int iep = ie * 3 + ip;
  LogDebug("GoodSeedProducer") << "Track pt =" << pt << " eta=" << eta << " bin=" << iep;
  return iep;
}

void GoodSeedProducer::fillPreIdRefValueMap(Handle<TrackCollection> tracks,
                                            const edm::OrphanHandle<reco::PreIdCollection>& preidhandle,
                                            edm::ValueMap<reco::PreIdRef>::Filler& filler) {
  std::vector<reco::PreIdRef> values;

  unsigned ntracks = tracks->size();
  for (unsigned itrack = 0; itrack < ntracks; ++itrack) {
    reco::TrackRef theTrackRef(tracks, itrack);
    std::map<reco::TrackRef, unsigned>::const_iterator itcheck = refMap_.find(theTrackRef);
    if (itcheck == refMap_.end()) {
      // the track has been early discarded
      values.push_back(reco::PreIdRef());
    } else {
      edm::Ref<reco::PreIdCollection> preIdRef(preidhandle, itcheck->second);
      values.push_back(preIdRef);
      //     std::cout << " Checking Refs " << (theTrackRef==preIdRef->trackRef()) << std::endl;
    }
  }
  filler.insert(tracks, values.begin(), values.end());
}

DEFINE_FWK_MODULE(GoodSeedProducer);