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