Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-03-18 02:30:00

0001 // -*- C++ -*-
0002 //
0003 // Package:    FastTrackDeDxProducer
0004 // Class:      FastTrackDeDxProducer
0005 //
0006 /**\class FastTrackDeDxProducer FastTrackDeDxProducer.cc RecoTracker/FastTrackDeDxProducer/src/FastTrackDeDxProducer.cc
0007 
0008    Description: <one line class summary>
0009 
0010    Implementation:
0011    <Notes on implementation>
0012 */
0013 // Original author:  Sam Bein
0014 //         Created:  Wednesday Dec 26 14:17:19 CEST 2018
0015 // Author of derivative code:  andrea
0016 //         Created:  Thu May 31 14:09:02 CEST 2007
0017 // Code Updates:  loic Quertenmont (querten)
0018 //         Created:  Thu May 10 14:09:02 CEST 2008
0019 //
0020 //
0021 
0022 #include <memory>
0023 
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/stream/EDProducer.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "FWCore/Framework/interface/ConsumesCollector.h"
0030 
0031 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0032 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0033 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0034 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0035 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0036 
0037 #include "DataFormats/Common/interface/ValueMap.h"
0038 #include "DataFormats/TrackReco/interface/DeDxData.h"
0039 #include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
0040 #include "DataFormats/TrackReco/interface/DeDxHit.h"
0041 #include "DataFormats/TrackReco/interface/Track.h"
0042 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0043 
0044 #include "RecoTracker/DeDx/interface/BaseDeDxEstimator.h"
0045 #include "RecoTracker/DeDx/interface/GenericAverageDeDxEstimator.h"
0046 #include "RecoTracker/DeDx/interface/TruncatedAverageDeDxEstimator.h"
0047 #include "RecoTracker/DeDx/interface/MedianDeDxEstimator.h"
0048 #include "RecoTracker/DeDx/interface/UnbinnedFitDeDxEstimator.h"
0049 #include "RecoTracker/DeDx/interface/ProductDeDxDiscriminator.h"
0050 #include "RecoTracker/DeDx/interface/SmirnovDeDxDiscriminator.h"
0051 #include "RecoTracker/DeDx/interface/ASmirnovDeDxDiscriminator.h"
0052 #include "RecoTracker/DeDx/interface/BTagLikeDeDxDiscriminator.h"
0053 
0054 #include "RecoTracker/DeDx/interface/DeDxTools.h"
0055 
0056 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0057 
0058 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0059 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0060 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
0061 #include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHitCollection.h"
0062 
0063 //
0064 // class declaration
0065 //
0066 
0067 class FastTrackDeDxProducer : public edm::stream::EDProducer<> {
0068 public:
0069   explicit FastTrackDeDxProducer(const edm::ParameterSet&);
0070   ~FastTrackDeDxProducer() override = default;
0071   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0072 
0073 private:
0074   void beginRun(edm::Run const&, const edm::EventSetup&) override;
0075   void produce(edm::Event&, const edm::EventSetup&) override;
0076 
0077   void makeCalibrationMap(const TrackerGeometry& tkGeom);
0078   void processHit(const FastTrackerRecHit& recHit,
0079                   float trackMomentum,
0080                   float& cosine,
0081                   reco::DeDxHitCollection& dedxHits,
0082                   int& NClusterSaturating);
0083 
0084   // ----------member data ---------------------------
0085 
0086   std::unique_ptr<BaseDeDxEstimator> m_estimator;
0087 
0088   edm::EDGetTokenT<reco::TrackCollection> m_tracksTag;
0089 
0090   float meVperADCPixel;
0091   float meVperADCStrip;
0092 
0093   unsigned int MaxNrStrips;
0094 
0095   std::string m_calibrationPath;
0096 
0097   std::vector<std::vector<float>> calibGains;
0098   unsigned int offsetDU_;
0099 
0100   edm::EDGetTokenT<edm::PSimHitContainer> simHitsToken;
0101   edm::EDGetTokenT<FastTrackerRecHitRefCollection> simHit2RecHitMapToken;
0102   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkGeomToken;
0103 
0104   bool usePixel;
0105   bool useStrip;
0106   bool useCalibration;
0107   bool shapetest;
0108   bool convertFromGeV2MeV;
0109   bool nothick;
0110 };
0111 
0112 using namespace reco;
0113 using namespace std;
0114 using namespace edm;
0115 
0116 void FastTrackDeDxProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0117   edm::ParameterSetDescription desc;
0118   desc.add<string>("estimator", "generic");
0119   desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks"));
0120   desc.add<bool>("UsePixel", false);
0121   desc.add<bool>("UseStrip", true);
0122   desc.add<double>("MeVperADCStrip", 3.61e-06 * 265);
0123   desc.add<double>("MeVperADCPixel", 3.61e-06);
0124   desc.add<bool>("ShapeTest", true);
0125   desc.add<bool>("UseCalibration", false);
0126   desc.add<string>("calibrationPath", "");
0127   desc.add<string>("Record", "SiStripDeDxMip_3D_Rcd");
0128   desc.add<string>("ProbabilityMode", "Accumulation");
0129   desc.add<double>("fraction", 0.4);
0130   desc.add<double>("exponent", -2.0);
0131   desc.add<bool>("convertFromGeV2MeV", true);
0132   desc.add<bool>("nothick", false);
0133   desc.add<edm::InputTag>("simHits");
0134   desc.add<edm::InputTag>("simHit2RecHitMap");
0135   descriptions.add("FastTrackDeDxProducer", desc);
0136 }
0137 
0138 FastTrackDeDxProducer::FastTrackDeDxProducer(const edm::ParameterSet& iConfig)
0139     : simHitsToken(consumes<edm::PSimHitContainer>(iConfig.getParameter<edm::InputTag>("simHits"))),
0140       simHit2RecHitMapToken(
0141           consumes<FastTrackerRecHitRefCollection>(iConfig.getParameter<edm::InputTag>("simHit2RecHitMap"))) {
0142   produces<ValueMap<DeDxData>>();
0143 
0144   auto cCollector = consumesCollector();
0145   string estimatorName = iConfig.getParameter<string>("estimator");
0146   if (estimatorName == "median")
0147     m_estimator = std::unique_ptr<BaseDeDxEstimator>(new MedianDeDxEstimator(iConfig));
0148   else if (estimatorName == "generic")
0149     m_estimator = std::unique_ptr<BaseDeDxEstimator>(new GenericAverageDeDxEstimator(iConfig));
0150   else if (estimatorName == "truncated")
0151     m_estimator = std::unique_ptr<BaseDeDxEstimator>(new TruncatedAverageDeDxEstimator(iConfig));
0152   //else if(estimatorName == "unbinnedFit")         m_estimator = std::unique_ptr<BaseDeDxEstimator> (new UnbinnedFitDeDxEstimator(iConfig));//estimator used in FullSimVersion
0153   else if (estimatorName == "productDiscrim")
0154     m_estimator = std::unique_ptr<BaseDeDxEstimator>(new ProductDeDxDiscriminator(iConfig, cCollector));
0155   else if (estimatorName == "btagDiscrim")
0156     m_estimator = std::unique_ptr<BaseDeDxEstimator>(new BTagLikeDeDxDiscriminator(iConfig, cCollector));
0157   else if (estimatorName == "smirnovDiscrim")
0158     m_estimator = std::unique_ptr<BaseDeDxEstimator>(new SmirnovDeDxDiscriminator(iConfig, cCollector));
0159   else if (estimatorName == "asmirnovDiscrim")
0160     m_estimator = std::unique_ptr<BaseDeDxEstimator>(new ASmirnovDeDxDiscriminator(iConfig, cCollector));
0161   else
0162     throw cms::Exception("fastsim::SimplifiedGeometry::FastTrackDeDxProducer.cc") << " estimator name does not exist";
0163 
0164   //Commented for now, might be used in the future
0165   //   MaxNrStrips         = iConfig.getUntrackedParameter<unsigned>("maxNrStrips"        ,  255);
0166 
0167   m_tracksTag = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks"));
0168 
0169   usePixel = iConfig.getParameter<bool>("UsePixel");
0170   useStrip = iConfig.getParameter<bool>("UseStrip");
0171   meVperADCPixel = iConfig.getParameter<double>("MeVperADCPixel");
0172   meVperADCStrip = iConfig.getParameter<double>("MeVperADCStrip");
0173 
0174   shapetest = iConfig.getParameter<bool>("ShapeTest");
0175   useCalibration = iConfig.getParameter<bool>("UseCalibration");
0176   m_calibrationPath = iConfig.getParameter<string>("calibrationPath");
0177 
0178   convertFromGeV2MeV = iConfig.getParameter<bool>("convertFromGeV2MeV");
0179   nothick = iConfig.getParameter<bool>("nothick");
0180 
0181   if (!usePixel && !useStrip)
0182     throw cms::Exception("fastsim::SimplifiedGeometry::FastTrackDeDxProducer.cc")
0183         << " neither pixel hits nor strips hits will be used to compute de/dx";
0184 
0185   if (useCalibration) {
0186     tkGeomToken = esConsumes<edm::Transition::BeginRun>();
0187   }
0188 }
0189 
0190 // ------------ method called once each job just before starting event loop  ------------
0191 void FastTrackDeDxProducer::beginRun(edm::Run const& run, const edm::EventSetup& iSetup) {
0192   if (useCalibration && calibGains.empty()) {
0193     auto const& tkGeom = iSetup.getData(tkGeomToken);
0194     offsetDU_ = tkGeom.offsetDU(GeomDetEnumerators::PixelBarrel);  //index start at the first pixel
0195 
0196     deDxTools::makeCalibrationMap(m_calibrationPath, tkGeom, calibGains, offsetDU_);
0197   }
0198 
0199   m_estimator->beginRun(run, iSetup);
0200 }
0201 
0202 void FastTrackDeDxProducer::produce(edm::Event& iEvent, const edm::EventSetup&) {
0203   auto trackDeDxEstimateAssociation = std::make_unique<ValueMap<DeDxData>>();
0204   ValueMap<DeDxData>::Filler filler(*trackDeDxEstimateAssociation);
0205 
0206   edm::Handle<reco::TrackCollection> trackCollectionHandle;
0207   iEvent.getByToken(m_tracksTag, trackCollectionHandle);
0208   const auto& trackCollection = *trackCollectionHandle;
0209   std::vector<DeDxData> dedxEstimate(trackCollection.size());
0210 
0211   for (unsigned int j = 0; j < trackCollection.size(); j++) {
0212     const reco::TrackRef track = reco::TrackRef(trackCollectionHandle.product(), j);
0213 
0214     int NClusterSaturating = 0;
0215     DeDxHitCollection dedxHits;
0216 
0217     auto const& trajParams = track->extra()->trajParams();
0218     assert(trajParams.size() == track->recHitsSize());
0219 
0220     auto hb = track->recHitsBegin();
0221     dedxHits.reserve(track->recHitsSize() / 2);
0222 
0223     for (unsigned int h = 0; h < track->recHitsSize(); h++) {
0224       const FastTrackerRecHit& recHit = static_cast<const FastTrackerRecHit&>(*(*(hb + h)));
0225       if (!recHit.isValid())
0226         continue;  //FastTrackerRecHit recHit = *(hb+h);
0227       auto trackDirection = trajParams[h].direction();
0228       float cosine = trackDirection.z() / trackDirection.mag();
0229       processHit(recHit, track->p(), cosine, dedxHits, NClusterSaturating);
0230     }
0231 
0232     sort(dedxHits.begin(), dedxHits.end(), less<DeDxHit>());
0233     std::pair<float, float> val_and_error = m_estimator->dedx(dedxHits);
0234     //WARNING: Since the dEdX Error is not properly computed for the moment
0235     //It was decided to store the number of saturating cluster in that dataformat
0236     val_and_error.second = NClusterSaturating;
0237     dedxEstimate[j] = DeDxData(val_and_error.first, val_and_error.second, dedxHits.size());
0238   }
0239 
0240   filler.insert(trackCollectionHandle, dedxEstimate.begin(), dedxEstimate.end());
0241   // fill the association map and put it into the event
0242   filler.fill();
0243   iEvent.put(std::move(trackDeDxEstimateAssociation));
0244 }
0245 
0246 void FastTrackDeDxProducer::processHit(const FastTrackerRecHit& recHit,
0247                                        float trackMomentum,
0248                                        float& cosine,
0249                                        reco::DeDxHitCollection& dedxHits,
0250                                        int& NClusterSaturating) {
0251   if (!recHit.isValid())
0252     return;
0253 
0254   auto const& thit = static_cast<BaseTrackerRecHit const&>(recHit);
0255   if (!thit.isValid())
0256     return;
0257   if (!thit.hasPositionAndError())
0258     return;
0259 
0260   if (recHit.isPixel()) {
0261     if (!usePixel)
0262       return;
0263 
0264     auto& detUnit = *(recHit.detUnit());
0265     float pathLen = detUnit.surface().bounds().thickness() / fabs(cosine);
0266     if (nothick)
0267       pathLen = 1.0;
0268     float charge = recHit.energyLoss() / pathLen;
0269     if (convertFromGeV2MeV)
0270       charge *= 1000;
0271     dedxHits.push_back(DeDxHit(charge, trackMomentum, pathLen, recHit.geographicalId()));
0272   } else if (!recHit.isPixel()) {  // && !recHit.isMatched()){//check what recHit.isMatched is doing
0273     if (!useStrip)
0274       return;
0275     auto& detUnit = *(recHit.detUnit());
0276     float pathLen = detUnit.surface().bounds().thickness() / fabs(cosine);
0277     if (nothick)
0278       pathLen = 1.0;
0279     float dedxOfRecHit = recHit.energyLoss() / pathLen;
0280     if (convertFromGeV2MeV)
0281       dedxOfRecHit *= 1000;
0282     if (!shapetest) {
0283       dedxHits.push_back(DeDxHit(dedxOfRecHit, trackMomentum, pathLen, recHit.geographicalId()));
0284     }
0285   }
0286 }
0287 
0288 //define this as a plug-in
0289 DEFINE_FWK_MODULE(FastTrackDeDxProducer);