DeDxHitInfoProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    DeDxHitInfoProducer
// Class:      DeDxHitInfoProducer
//
//
// Original Author:  loic Quertenmont (querten)
//         Created:  Mon Nov 21 14:09:02 CEST 2014
//

#include "RecoTracker/DeDx/plugins/DeDxHitInfoProducer.h"

// system include files

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

DeDxHitInfoProducer::DeDxHitInfoProducer(const edm::ParameterSet& iConfig)
    : usePixel(iConfig.getParameter<bool>("usePixel")),
      useStrip(iConfig.getParameter<bool>("useStrip")),
      MeVperADCPixel(iConfig.getParameter<double>("MeVperADCPixel")),
      MeVperADCStrip(iConfig.getParameter<double>("MeVperADCStrip")),
      minTrackHits(iConfig.getParameter<unsigned>("minTrackHits")),
      minTrackPt(iConfig.getParameter<double>("minTrackPt")),
      minTrackPtPrescale(iConfig.getParameter<double>("minTrackPtPrescale")),
      maxTrackEta(iConfig.getParameter<double>("maxTrackEta")),
      m_calibrationPath(iConfig.getParameter<string>("calibrationPath")),
      useCalibration(iConfig.getParameter<bool>("useCalibration")),
      shapetest(iConfig.getParameter<bool>("shapeTest")),
      lowPtTracksPrescalePass(iConfig.getParameter<uint32_t>("lowPtTracksPrescalePass")),
      lowPtTracksPrescaleFail(iConfig.getParameter<uint32_t>("lowPtTracksPrescaleFail")),
      lowPtTracksEstimator(iConfig.getParameter<edm::ParameterSet>("lowPtTracksEstimatorParameters")),
      lowPtTracksDeDxThreshold(iConfig.getParameter<double>("lowPtTracksDeDxThreshold")) {
  produces<reco::DeDxHitInfoCollection>();
  produces<reco::DeDxHitInfoAss>();
  produces<edm::ValueMap<int>>("prescale");

  m_tracksTag = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));

  if (!usePixel && !useStrip)
    edm::LogError("DeDxHitsProducer") << "No Pixel Hits NOR Strip Hits will be saved.  Running this module is useless";
}

DeDxHitInfoProducer::~DeDxHitInfoProducer() {}

// ------------ method called once each job just before starting event loop  ------------
void DeDxHitInfoProducer::beginRun(edm::Run const& run, const edm::EventSetup& iSetup) {
  iSetup.get<TrackerDigiGeometryRecord>().get(tkGeom);
  if (useCalibration && calibGains.empty()) {
    m_off = tkGeom->offsetDU(GeomDetEnumerators::PixelBarrel);  //index start at the first pixel

    DeDxTools::makeCalibrationMap(m_calibrationPath, *tkGeom, calibGains, m_off);
  }
}

void DeDxHitInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  edm::Handle<reco::TrackCollection> trackCollectionHandle;
  iEvent.getByToken(m_tracksTag, trackCollectionHandle);
  const TrackCollection& trackCollection(*trackCollectionHandle.product());

  // creates the output collection
  auto resultdedxHitColl = std::make_unique<reco::DeDxHitInfoCollection>();

  std::vector<int> indices;
  std::vector<int> prescales;
  uint64_t state[2] = {iEvent.id().event(), iEvent.id().luminosityBlock()};
  for (unsigned int j = 0; j < trackCollection.size(); j++) {
    const reco::Track& track = trackCollection[j];

    //track selection
    bool passPt = (track.pt() >= minTrackPt), passLowDeDx = false, passHighDeDx = false, pass = passPt;
    if (!pass && (track.pt() >= minTrackPtPrescale)) {
      if (lowPtTracksPrescalePass > 0) {
        passHighDeDx = ((xorshift128p(state) % lowPtTracksPrescalePass) == 0);
      }
      if (lowPtTracksPrescaleFail > 0) {
        passLowDeDx = ((xorshift128p(state) % lowPtTracksPrescaleFail) == 0);
      }
      pass = passHighDeDx || passLowDeDx;
    }
    if (!pass || std::abs(track.eta()) > maxTrackEta || track.numberOfValidHits() < minTrackHits) {
      indices.push_back(-1);
      continue;
    }

    reco::DeDxHitInfo hitDeDxInfo;
    auto const& trajParams = track.extra()->trajParams();
    auto hb = track.recHitsBegin();
    for (unsigned int h = 0; h < track.recHitsSize(); h++) {
      auto recHit = *(hb + h);
      if (!trackerHitRTTI::isFromDet(*recHit))
        continue;

      auto trackDirection = trajParams[h].direction();
      float cosine = trackDirection.z() / trackDirection.mag();

      processHit(recHit, track.p(), cosine, hitDeDxInfo, trajParams[h].position());
    }

    if (!passPt) {
      std::vector<DeDxHit> hits;
      hits.reserve(hitDeDxInfo.size());
      for (unsigned int i = 0, n = hitDeDxInfo.size(); i < n; ++i) {
        if (hitDeDxInfo.detId(i).subdetId() <= 2) {
          hits.push_back(DeDxHit(hitDeDxInfo.charge(i) / hitDeDxInfo.pathlength(i) * MeVperADCPixel, 0, 0, 0));
        } else {
          if (shapetest && !DeDxTools::shapeSelection(*hitDeDxInfo.stripCluster(i)))
            continue;
          hits.push_back(DeDxHit(hitDeDxInfo.charge(i) / hitDeDxInfo.pathlength(i) * MeVperADCStrip, 0, 0, 0));
        }
      }
      std::sort(hits.begin(), hits.end(), std::less<DeDxHit>());
      if (lowPtTracksEstimator.dedx(hits).first < lowPtTracksDeDxThreshold) {
        if (passLowDeDx) {
          prescales.push_back(lowPtTracksPrescaleFail);
        } else {
          indices.push_back(-1);
          continue;
        }
      } else {
        if (passHighDeDx) {
          prescales.push_back(lowPtTracksPrescalePass);
        } else {
          indices.push_back(-1);
          continue;
        }
      }
    } else {
      prescales.push_back(1);
    }
    indices.push_back(resultdedxHitColl->size());
    resultdedxHitColl->push_back(hitDeDxInfo);
  }

  edm::OrphanHandle<reco::DeDxHitInfoCollection> dedxHitCollHandle = iEvent.put(std::move(resultdedxHitColl));

  //create map passing the handle to the matched collection
  auto dedxMatch = std::make_unique<reco::DeDxHitInfoAss>(dedxHitCollHandle);
  reco::DeDxHitInfoAss::Filler filler(*dedxMatch);
  filler.insert(trackCollectionHandle, indices.begin(), indices.end());
  filler.fill();
  iEvent.put(std::move(dedxMatch));

  auto dedxPrescale = std::make_unique<edm::ValueMap<int>>();
  edm::ValueMap<int>::Filler pfiller(*dedxPrescale);
  pfiller.insert(dedxHitCollHandle, prescales.begin(), prescales.end());
  pfiller.fill();
  iEvent.put(std::move(dedxPrescale), "prescale");
}

void DeDxHitInfoProducer::processHit(const TrackingRecHit* recHit,
                                     const float trackMomentum,
                                     const float cosine,
                                     reco::DeDxHitInfo& hitDeDxInfo,
                                     const LocalPoint& hitLocalPos) {
  auto const& thit = static_cast<BaseTrackerRecHit const&>(*recHit);
  if (!thit.isValid())
    return;

  float cosineAbs = std::max(0.00000001f, std::abs(cosine));  //make sure cosine is not 0

  auto const& clus = thit.firstClusterRef();
  if (!clus.isValid())
    return;

  if (clus.isPixel()) {
    if (!usePixel)
      return;

    const auto* detUnit = recHit->detUnit();
    if (detUnit == nullptr)
      detUnit = tkGeom->idToDet(thit.geographicalId());
    float pathLen = detUnit->surface().bounds().thickness() / cosineAbs;
    float chargeAbs = clus.pixelCluster().charge();
    hitDeDxInfo.addHit(chargeAbs, pathLen, thit.geographicalId(), hitLocalPos, clus.pixelCluster());
  } else if (clus.isStrip() && !thit.isMatched()) {
    if (!useStrip)
      return;

    const auto* detUnit = recHit->detUnit();
    if (detUnit == nullptr)
      detUnit = tkGeom->idToDet(thit.geographicalId());
    int NSaturating = 0;
    float pathLen = detUnit->surface().bounds().thickness() / cosineAbs;
    float chargeAbs = DeDxTools::getCharge(&(clus.stripCluster()), NSaturating, *detUnit, calibGains, m_off);
    hitDeDxInfo.addHit(chargeAbs, pathLen, thit.geographicalId(), hitLocalPos, clus.stripCluster());
  } else if (clus.isStrip() && thit.isMatched()) {
    if (!useStrip)
      return;
    const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
    if (!matchedHit)
      return;

    const auto* detUnitM = matchedHit->monoHit().detUnit();
    if (detUnitM == nullptr)
      detUnitM = tkGeom->idToDet(matchedHit->monoHit().geographicalId());
    int NSaturating = 0;
    float pathLen = detUnitM->surface().bounds().thickness() / cosineAbs;
    float chargeAbs =
        DeDxTools::getCharge(&(matchedHit->monoHit().stripCluster()), NSaturating, *detUnitM, calibGains, m_off);
    hitDeDxInfo.addHit(chargeAbs, pathLen, thit.geographicalId(), hitLocalPos, matchedHit->monoHit().stripCluster());

    const auto* detUnitS = matchedHit->stereoHit().detUnit();
    if (detUnitS == nullptr)
      detUnitS = tkGeom->idToDet(matchedHit->stereoHit().geographicalId());
    NSaturating = 0;
    pathLen = detUnitS->surface().bounds().thickness() / cosineAbs;
    chargeAbs =
        DeDxTools::getCharge(&(matchedHit->stereoHit().stripCluster()), NSaturating, *detUnitS, calibGains, m_off);
    hitDeDxInfo.addHit(chargeAbs, pathLen, thit.geographicalId(), hitLocalPos, matchedHit->stereoHit().stripCluster());
  }
}

//define this as a plug-in
DEFINE_FWK_MODULE(DeDxHitInfoProducer);