DeDx/plugins/DeDxEstimatorProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    DeDxEstimatorProducer
0004 // Class:      DeDxEstimatorProducer
0005 //
0006 /**\class DeDxEstimatorProducer DeDxEstimatorProducer.cc RecoTracker/DeDxEstimatorProducer/src/DeDxEstimatorProducer.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  andrea
0015 //         Created:  Thu May 31 14:09:02 CEST 2007
0016 //    Code Updates:  loic Quertenmont (querten)
0017 //         Created:  Thu May 10 14:09:02 CEST 2008
0018 //
0019 //
0020
0021 // system include files
0022 #include "RecoTracker/DeDx/plugins/DeDxEstimatorProducer.h"
0023 #include "FWCore/Framework/interface/ConsumesCollector.h"
0024
0025 using namespace reco;
0026 using namespace std;
0027 using namespace edm;
0028
0029 void DeDxEstimatorProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0030   edm::ParameterSetDescription desc;
0031   desc.add<string>("estimator", "generic");
0032   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0033   desc.add<bool>("UseDeDxHits", false);
0034   desc.add<edm::InputTag>("pixelDeDxHits", {});
0035   desc.add<edm::InputTag>("stripDeDxHits", {});
0036   desc.add<bool>("UsePixel", false);
0037   desc.add<bool>("UseStrip", true);
0038   desc.add<double>("MeVperADCPixel", 3.61e-06);
0039   desc.add<double>("MeVperADCStrip", 3.61e-06 * 265);
0040   desc.add<bool>("ShapeTest", true);
0041   desc.add<bool>("UseCalibration", false);
0042   desc.add<string>("calibrationPath", "");
0043   desc.add<string>("Record", "SiStripDeDxMip_3D_Rcd");
0044   desc.add<string>("ProbabilityMode", "Accumulation");
0045   desc.add<double>("fraction", 0.4);
0046   desc.add<double>("exponent", -2.0);
0047   desc.add<bool>("truncate", true);
0048
0049   descriptions.add("DeDxEstimatorProducer", desc);
0050 }
0051
0052 DeDxEstimatorProducer::DeDxEstimatorProducer(const edm::ParameterSet& iConfig)
0053     : tkGeomToken(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()) {
0054   produces<ValueMap<DeDxData>>();
0055
0056   auto cCollector = consumesCollector();
0057   string estimatorName = iConfig.getParameter<string>("estimator");
0058   if (estimatorName == "median")
0059     m_estimator = std::make_unique<MedianDeDxEstimator>(iConfig);
0060   else if (estimatorName == "generic")
0061     m_estimator = std::make_unique<GenericAverageDeDxEstimator>(iConfig);
0062   else if (estimatorName == "truncated")
0063     m_estimator = std::make_unique<TruncatedAverageDeDxEstimator>(iConfig);
0064   else if (estimatorName == "genericTruncated")
0065     m_estimator = std::make_unique<GenericTruncatedAverageDeDxEstimator>(iConfig);
0066   else if (estimatorName == "unbinnedFit")
0067     m_estimator = std::make_unique<UnbinnedFitDeDxEstimator>(iConfig);
0068   else if (estimatorName == "likelihoodFit")
0069     m_estimator = std::make_unique<LikelihoodFitDeDxEstimator>(iConfig);
0070   else if (estimatorName == "productDiscrim")
0071     m_estimator = std::make_unique<ProductDeDxDiscriminator>(iConfig, cCollector);
0072   else if (estimatorName == "btagDiscrim")
0073     m_estimator = std::make_unique<BTagLikeDeDxDiscriminator>(iConfig, cCollector);
0074   else if (estimatorName == "smirnovDiscrim")
0075     m_estimator = std::make_unique<SmirnovDeDxDiscriminator>(iConfig, cCollector);
0076   else if (estimatorName == "asmirnovDiscrim")
0077     m_estimator = std::make_unique<ASmirnovDeDxDiscriminator>(iConfig, cCollector);
0078
0079   //Commented for now, might be used in the future
0080   //   MaxNrStrips         = iConfig.getUntrackedParameter<unsigned>("maxNrStrips"        ,  255);
0081
0082   m_tracksTag = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0083
0084   useDeDxHits = iConfig.getParameter<bool>("UseDeDxHits");
0085   m_pixelDeDxHitsTag = consumes<reco::TrackDeDxHitsCollection>(iConfig.getParameter<edm::InputTag>("pixelDeDxHits"));
0086   m_stripDeDxHitsTag = consumes<reco::TrackDeDxHitsCollection>(iConfig.getParameter<edm::InputTag>("stripDeDxHits"));
0087
0088   usePixel = iConfig.getParameter<bool>("UsePixel");
0089   useStrip = iConfig.getParameter<bool>("UseStrip");
0090   meVperADCPixel = iConfig.getParameter<double>("MeVperADCPixel");
0091   meVperADCStrip = iConfig.getParameter<double>("MeVperADCStrip");
0092
0093   shapetest = iConfig.getParameter<bool>("ShapeTest");
0094   useCalibration = iConfig.getParameter<bool>("UseCalibration");
0095   m_calibrationPath = iConfig.getParameter<string>("calibrationPath");
0096
0097   if (!usePixel && !useStrip)
0098     edm::LogWarning("DeDxHitsProducer")
0099         << "Pixel Hits AND Strip Hits will not be used to estimate dEdx --> BUG, Please Update the config file";
0100 }
0101
0102 DeDxEstimatorProducer::~DeDxEstimatorProducer() {}
0103
0104 // ------------ method called once each job just before starting event loop  ------------
0105 void DeDxEstimatorProducer::beginRun(edm::Run const& run, const edm::EventSetup& iSetup) {
0106   tkGeom = &iSetup.getData(tkGeomToken);
0107   if (useCalibration && calibGains.empty()) {
0108     m_off = tkGeom->offsetDU(GeomDetEnumerators::PixelBarrel);  //index start at the first pixel
0109
0110     deDxTools::makeCalibrationMap(m_calibrationPath, *tkGeom, calibGains, m_off);
0111   }
0112
0113   m_estimator->beginRun(run, iSetup);
0114 }
0115
0116 void DeDxEstimatorProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0117   auto trackDeDxEstimateAssociation = std::make_unique<ValueMap<DeDxData>>();
0118   ValueMap<DeDxData>::Filler filler(*trackDeDxEstimateAssociation);
0119
0120   edm::Handle<reco::TrackCollection> trackCollectionHandle;
0121   iEvent.getByToken(m_tracksTag, trackCollectionHandle);
0122   const auto& pixelDeDxHits = iEvent.getHandle(m_pixelDeDxHitsTag);
0123   const auto& stripDeDxHits = iEvent.getHandle(m_stripDeDxHitsTag);
0124
0125   std::vector<DeDxData> dedxEstimate(trackCollectionHandle->size());
0126
0127   for (unsigned int j = 0; j < trackCollectionHandle->size(); j++) {
0128     const reco::TrackRef track = reco::TrackRef(trackCollectionHandle, j);
0129
0130     int NClusterSaturating = 0;
0131     DeDxHitCollection dedxHits;
0132
0133     if (useDeDxHits) {
0134       if (usePixel)
0135         dedxHits = (*pixelDeDxHits)[track];
0136       if (useStrip) {
0137         const auto& hits = (*stripDeDxHits)[track];
0138         dedxHits.insert(dedxHits.end(), hits.begin(), hits.end());
0139       }
0140     } else {
0141       auto const& trajParams = track->extra()->trajParams();
0142       assert(trajParams.size() == track->recHitsSize());
0143       auto hb = track->recHitsBegin();
0144       dedxHits.reserve(track->recHitsSize() / 2);
0145       for (unsigned int h = 0; h < track->recHitsSize(); h++) {
0146         auto recHit = *(hb + h);
0147         if (!trackerHitRTTI::isFromDet(*recHit))
0148           continue;
0149
0150         auto trackDirection = trajParams[h].direction();
0151         float cosine = trackDirection.z() / trackDirection.mag();
0152         processHit(recHit, track->p(), cosine, dedxHits, NClusterSaturating);
0153       }
0154     }
0155
0156     sort(dedxHits.begin(), dedxHits.end(), less<DeDxHit>());
0157     std::pair<float, float> val_and_error = m_estimator->dedx(dedxHits);
0158
0159     //WARNING: Since the dEdX Error is not properly computed for the moment
0160     //It was decided to store the number of saturating cluster in that dataformat
0161     val_and_error.second = NClusterSaturating;
0162     dedxEstimate[j] = DeDxData(val_and_error.first, val_and_error.second, dedxHits.size());
0163   }
0164   ///////////////////////////////////////
0165
0166   filler.insert(trackCollectionHandle, dedxEstimate.begin(), dedxEstimate.end());
0167
0168   // fill the association map and put it into the event
0169   filler.fill();
0170   iEvent.put(std::move(trackDeDxEstimateAssociation));
0171 }
0172
0173 void DeDxEstimatorProducer::processHit(const TrackingRecHit* recHit,
0174                                        float trackMomentum,
0175                                        float& cosine,
0176                                        reco::DeDxHitCollection& dedxHits,
0177                                        int& NClusterSaturating) {
0178   auto const& thit = static_cast<BaseTrackerRecHit const&>(*recHit);
0179   if (!thit.isValid())
0180     return;
0181
0182   auto const& clus = thit.firstClusterRef();
0183   if (!clus.isValid())
0184     return;
0185
0186   if (clus.isPixel()) {
0187     if (!usePixel)
0188       return;
0189
0190     const auto* detUnit = recHit->detUnit();
0191     if (detUnit == nullptr)
0192       detUnit = tkGeom->idToDet(thit.geographicalId());
0193     float pathLen = detUnit->surface().bounds().thickness() / fabs(cosine);
0194     float chargeAbs = clus.pixelCluster().charge();
0195     float charge = meVperADCPixel * chargeAbs / pathLen;
0196     dedxHits.push_back(DeDxHit(charge, trackMomentum, pathLen, thit.geographicalId()));
0197   } else if (clus.isStrip() && !thit.isMatched()) {
0198     if (!useStrip)
0199       return;
0200
0201     const auto* detUnit = recHit->detUnit();
0202     if (detUnit == nullptr)
0203       detUnit = tkGeom->idToDet(thit.geographicalId());
0204     int NSaturating = 0;
0205     float pathLen = detUnit->surface().bounds().thickness() / fabs(cosine);
0206     float chargeAbs = deDxTools::getCharge(&(clus.stripCluster()), NSaturating, *detUnit, calibGains, m_off);
0207     float charge = meVperADCStrip * chargeAbs / pathLen;
0208     if (!shapetest || (shapetest && deDxTools::shapeSelection(clus.stripCluster()))) {
0209       dedxHits.push_back(DeDxHit(charge, trackMomentum, pathLen, thit.geographicalId()));
0210       if (NSaturating > 0)
0211         NClusterSaturating++;
0212     }
0213   } else if (clus.isStrip() && thit.isMatched()) {
0214     if (!useStrip)
0215       return;
0216     const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(recHit);
0217     if (!matchedHit)
0218       return;
0219     const GluedGeomDet* gdet = static_cast<const GluedGeomDet*>(matchedHit->det());
0220     if (gdet == nullptr)
0221       gdet = static_cast<const GluedGeomDet*>(tkGeom->idToDet(thit.geographicalId()));
0222
0223     auto& detUnitM = *(gdet->monoDet());
0224     int NSaturating = 0;
0225     float pathLen = detUnitM.surface().bounds().thickness() / fabs(cosine);
0226     float chargeAbs = deDxTools::getCharge(&(matchedHit->monoCluster()), NSaturating, detUnitM, calibGains, m_off);
0227     float charge = meVperADCStrip * chargeAbs / pathLen;
0228     if (!shapetest || (shapetest && deDxTools::shapeSelection(matchedHit->monoCluster()))) {
0229       dedxHits.push_back(DeDxHit(charge, trackMomentum, pathLen, matchedHit->monoId()));
0230       if (NSaturating > 0)
0231         NClusterSaturating++;
0232     }
0233
0234     auto& detUnitS = *(gdet->stereoDet());
0235     NSaturating = 0;
0236     pathLen = detUnitS.surface().bounds().thickness() / fabs(cosine);
0237     chargeAbs = deDxTools::getCharge(&(matchedHit->stereoCluster()), NSaturating, detUnitS, calibGains, m_off);
0238     charge = meVperADCStrip * chargeAbs / pathLen;
0239     if (!shapetest || (shapetest && deDxTools::shapeSelection(matchedHit->stereoCluster()))) {
0240       dedxHits.push_back(DeDxHit(charge, trackMomentum, pathLen, matchedHit->stereoId()));
0241       if (NSaturating > 0)
0242         NClusterSaturating++;
0243     }
0244   }
0245 }
0246
0247 //define this as a plug-in
0248 DEFINE_FWK_MODULE(DeDxEstimatorProducer);