HSCPDeDxInfoProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    HSCPDeDxInfoProducer
// Class:      HSCPDeDxInfoProducer
//
//
// Original Author:  andrea
//         Created:  Thu May 31 14:09:02 CEST 2007
//    Code Updates:  loic Quertenmont (querten)
//         Created:  Thu May 10 14:09:02 CEST 2008
//
//

#include "SUSYBSMAnalysis/HSCP/plugins/HSCPDeDxInfoProducer.h"

// system include files

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

HSCPDeDxInfoProducer::HSCPDeDxInfoProducer(const edm::ParameterSet& iConfig) {
  produces<ValueMap<susybsm::HSCPDeDxInfo> >();

  MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips", 255);
  MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits", 4);
  MinTrackMomentum = iConfig.getUntrackedParameter<double>("minTrackMomentum", 0.0);
  MaxTrackMomentum = iConfig.getUntrackedParameter<double>("maxTrackMomentum", 99999.0);
  MinTrackEta = iConfig.getUntrackedParameter<double>("minTrackEta", -5.0);
  MaxTrackEta = iConfig.getUntrackedParameter<double>("maxTrackEta", 5.0);

  m_tracksTag = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
  m_trajTrackAssociationTag =
      consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("trajectoryTrackAssociation"));
  useTrajectory = iConfig.getUntrackedParameter<bool>("UseTrajectory", true);

  usePixel = iConfig.getParameter<bool>("UsePixel");
  useStrip = iConfig.getParameter<bool>("UseStrip");
  meVperADCPixel = iConfig.getParameter<double>("MeVperADCPixel");
  meVperADCStrip = iConfig.getParameter<double>("MeVperADCStrip");

  shapetest = iConfig.getParameter<bool>("ShapeTest");
  useCalibration = iConfig.getParameter<bool>("UseCalibration");
  m_calibrationPath = iConfig.getParameter<string>("calibrationPath");

  //   Reccord             = iConfig.getUntrackedParameter<std::string>  ("Reccord"            , "SiStripDeDxMip_3D_Rcd");
  //   ProbabilityMode     = iConfig.getUntrackedParameter<std::string>  ("ProbabilityMode"    , "Accumulation");
  //   Prob_ChargePath     = NULL;

  if (!usePixel && !useStrip)
    edm::LogWarning("DeDxHitsProducer")
        << "Pixel Hits AND Strip Hits will not be used to estimate dEdx --> BUG, Please Update the config file";
}

HSCPDeDxInfoProducer::~HSCPDeDxInfoProducer() {}

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

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

  //   DeDxTools::buildDiscrimMap(run, iSetup, Reccord,  ProbabilityMode, Prob_ChargePath);
}

void HSCPDeDxInfoProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  unique_ptr<ValueMap<susybsm::HSCPDeDxInfo> > trackDeDxAssociation(new ValueMap<susybsm::HSCPDeDxInfo>);
  ValueMap<susybsm::HSCPDeDxInfo>::Filler filler(*trackDeDxAssociation);

  edm::Handle<reco::TrackCollection> trackCollectionHandle;
  iEvent.getByToken(m_tracksTag, trackCollectionHandle);

  Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
  if (useTrajectory)
    iEvent.getByToken(m_trajTrackAssociationTag, trajTrackAssociationHandle);

  std::vector<susybsm::HSCPDeDxInfo> dEdxInfos(trackCollectionHandle->size());

  TrajTrackAssociationCollection::const_iterator cit;
  if (useTrajectory)
    cit = trajTrackAssociationHandle->begin();
  for (unsigned int j = 0; j < trackCollectionHandle->size(); j++) {
    const reco::TrackRef track = reco::TrackRef(trackCollectionHandle.product(), j);

    susybsm::HSCPDeDxInfo hscpDeDxInfo;

    if (useTrajectory) {  //trajectory allows to take into account the local direction of the particle on the module sensor --> muc much better 'dx' measurement
      const edm::Ref<std::vector<Trajectory> > traj = cit->key;
      cit++;
      const vector<TrajectoryMeasurement>& measurements = traj->measurements();
      for (vector<TrajectoryMeasurement>::const_iterator it = measurements.begin(); it != measurements.end(); it++) {
        TrajectoryStateOnSurface trajState = it->updatedState();
        if (!trajState.isValid())
          continue;

        const TrackingRecHit* recHit = (*it->recHit()).hit();
        if (!recHit)
          continue;
        LocalVector trackDirection = trajState.localDirection();
        float cosine = trackDirection.z() / trackDirection.mag();

        processHit(recHit, trajState.localMomentum().mag(), cosine, hscpDeDxInfo, trajState.localPosition());
      }

    } else {  //assume that the particles trajectory is a straight line originating from the center of the detector  (can be improved)
      for (unsigned int h = 0; h < track->recHitsSize(); h++) {
        const TrackingRecHit* recHit = &(*(track->recHit(h)));
        auto const& thit = static_cast<BaseTrackerRecHit const&>(*recHit);
        if (!thit.isValid())
          continue;  //make sure it's a tracker hit

        const GlobalVector& ModuleNormal = recHit->detUnit()->surface().normalVector();
        float cosine =
            (track->px() * ModuleNormal.x() + track->py() * ModuleNormal.y() + track->pz() * ModuleNormal.z()) /
            track->p();

        processHit(recHit, track->p(), cosine, hscpDeDxInfo, LocalPoint(0.0, 0.0));
      }
    }

    dEdxInfos[j] = hscpDeDxInfo;
  }

  filler.insert(trackCollectionHandle, dEdxInfos.begin(), dEdxInfos.end());
  filler.fill();
  iEvent.put(std::move(trackDeDxAssociation));
}

void HSCPDeDxInfoProducer::processHit(const TrackingRecHit* recHit,
                                      float trackMomentum,
                                      float& cosine,
                                      susybsm::HSCPDeDxInfo& hscpDeDxInfo,
                                      LocalPoint HitLocalPos) {
  auto const& thit = static_cast<BaseTrackerRecHit const&>(*recHit);
  if (!thit.isValid())
    return;

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

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

    auto& detUnit = *(recHit->detUnit());
    float pathLen = detUnit.surface().bounds().thickness() / fabs(cosine);
    float chargeAbs = clus.pixelCluster().charge();
    hscpDeDxInfo.charges.push_back(chargeAbs);
    hscpDeDxInfo.pathlengths.push_back(pathLen);
    hscpDeDxInfo.detIds.push_back(thit.geographicalId());
    hscpDeDxInfo.localPosXs.push_back(HitLocalPos.x());
    hscpDeDxInfo.localPosYs.push_back(HitLocalPos.y());
    hscpDeDxInfo.clusterIndices.push_back(clus.key());
  } else if (clus.isStrip() && !thit.isMatched()) {
    if (!useStrip)
      return;

    auto& detUnit = *(recHit->detUnit());
    int NSaturating = 0;
    float pathLen = detUnit.surface().bounds().thickness() / fabs(cosine);
    float chargeAbs = DeDxTools::getCharge(&(clus.stripCluster()), NSaturating, detUnit, calibGains, m_off);
    hscpDeDxInfo.charges.push_back(chargeAbs);
    hscpDeDxInfo.pathlengths.push_back(pathLen);
    hscpDeDxInfo.detIds.push_back(thit.geographicalId());
    hscpDeDxInfo.localPosXs.push_back(HitLocalPos.x());
    hscpDeDxInfo.localPosYs.push_back(HitLocalPos.y());
    hscpDeDxInfo.clusterIndices.push_back(clus.key());
  } else if (clus.isStrip() && thit.isMatched()) {
    if (!useStrip)
      return;
    const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
    if (!matchedHit)
      return;

    auto& detUnitM = *(matchedHit->monoHit().detUnit());
    int NSaturating = 0;
    float pathLen = detUnitM.surface().bounds().thickness() / fabs(cosine);
    float chargeAbs =
        DeDxTools::getCharge(&(matchedHit->monoHit().stripCluster()), NSaturating, detUnitM, calibGains, m_off);
    hscpDeDxInfo.charges.push_back(chargeAbs);
    hscpDeDxInfo.pathlengths.push_back(pathLen);
    hscpDeDxInfo.detIds.push_back(thit.geographicalId());
    hscpDeDxInfo.localPosXs.push_back(HitLocalPos.x());
    hscpDeDxInfo.localPosYs.push_back(HitLocalPos.y());
    const OmniClusterRef monoClusterRef = matchedHit->monoClusterRef();
    hscpDeDxInfo.clusterIndices.push_back(monoClusterRef.key());

    auto& detUnitS = *(matchedHit->stereoHit().detUnit());
    NSaturating = 0;
    pathLen = detUnitS.surface().bounds().thickness() / fabs(cosine);
    chargeAbs =
        DeDxTools::getCharge(&(matchedHit->stereoHit().stripCluster()), NSaturating, detUnitS, calibGains, m_off);
    hscpDeDxInfo.charges.push_back(chargeAbs);
    hscpDeDxInfo.pathlengths.push_back(pathLen);
    hscpDeDxInfo.detIds.push_back(thit.geographicalId());
    hscpDeDxInfo.localPosXs.push_back(HitLocalPos.x());
    hscpDeDxInfo.localPosYs.push_back(HitLocalPos.y());
    const OmniClusterRef stereoClusterRef = matchedHit->stereoClusterRef();
    hscpDeDxInfo.clusterIndices.push_back(stereoClusterRef.key());
  }
}

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