DeDx/interface/ProductDeDxDiscriminator.h

0001 #ifndef RecoTrackerDeDx_ProductDeDxDiscriminator_h
0002 #define RecoTrackerDeDx_ProductDeDxDiscriminator_h
0003
0004 #include "RecoTracker/DeDx/interface/BaseDeDxEstimator.h"
0005 #include "RecoTracker/DeDx/interface/DeDxTools.h"
0006 #include "DataFormats/TrackReco/interface/DeDxHit.h"
0007
0008 class ProductDeDxDiscriminator : public BaseDeDxEstimator {
0009 public:
0010   ProductDeDxDiscriminator(const edm::ParameterSet& iConfig, edm::ConsumesCollector& iCollector)
0011       : token_(deDxTools::esConsumes(iConfig.getParameter<std::string>("Record"), iCollector)) {
0012     meVperADCStrip =
0013         iConfig.getParameter<double>("MeVperADCStrip");  //currently needed until the map on the database are redone
0014     ProbabilityMode = iConfig.getParameter<std::string>("ProbabilityMode");
0015     Prob_ChargePath = nullptr;
0016   }
0017
0018   void beginRun(edm::Run const&, const edm::EventSetup& iSetup) override {
0019     auto const& histD3D = deDxTools::getHistogramD3D(iSetup, token_);
0020     deDxTools::buildDiscrimMap(histD3D, ProbabilityMode, Prob_ChargePath);
0021   }
0022
0023   std::pair<float, float> dedx(const reco::DeDxHitCollection& Hits) override {
0024     std::vector<float> vect_probs;
0025     for (size_t i = 0; i < Hits.size(); i++) {
0026       float path = Hits[i].pathLength() * 10.0;  //x10 in order to be compatible with the map content
0027       float charge =
0028           Hits[i].charge() /
0029           (10.0 *
0030            meVperADCStrip);  // 10/meVperADCStrip in order to be compatible with the map content in ADC/mm instead of MeV/cm
0031
0032       int BinX = Prob_ChargePath->GetXaxis()->FindBin(Hits[i].momentum());
0033       int BinY = Prob_ChargePath->GetYaxis()->FindBin(path);
0034       int BinZ = Prob_ChargePath->GetZaxis()->FindBin(charge);
0035       float prob = Prob_ChargePath->GetBinContent(BinX, BinY, BinZ);
0036       if (prob >= 0)
0037         vect_probs.push_back(prob);
0038     }
0039
0040     size_t size = vect_probs.size();
0041     if (size <= 0)
0042       return std::make_pair(-1, -1);
0043     float TotalProb = 1;
0044     for (size_t i = 0; i < size; i++) {
0045       if (vect_probs[i] <= 0.0001) {
0046         TotalProb *= pow(0.0001, 1.0 / size);
0047       } else {
0048         TotalProb *= pow(vect_probs[i], 1.0 / size);
0049       }
0050     }
0051     return std::make_pair(TotalProb, -1);
0052   }
0053
0054 private:
0055   float meVperADCStrip;
0056   deDxTools::ESGetTokenH3DDVariant token_;
0057   std::string ProbabilityMode;
0058   TH3F* Prob_ChargePath;
0059 };
0060
0061 #endif