DeDx/interface/BTagLikeDeDxDiscriminator.h

0001 #ifndef RecoTrackerDeDx_BTagLikeDeDxDiscriminator_h
0002 #define RecoTrackerDeDx_BTagLikeDeDxDiscriminator_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 BTagLikeDeDxDiscriminator : public BaseDeDxEstimator {
0009 public:
0010   BTagLikeDeDxDiscriminator(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& run, 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     std::sort(vect_probs.begin(), vect_probs.end(), std::less<float>());
0044     float SumJet = 0.;
0045     for (size_t i = 0; i < size; i++) {
0046       if (vect_probs[i] <= 0.0001)
0047         vect_probs[i] = 0.0001;
0048       SumJet += log(vect_probs[i]);
0049     }
0050
0051     float Loginvlog = log(-SumJet);
0052     float Prob = 1;
0053     float lfact = 1;
0054
0055     for (size_t l = 1; l != size; l++) {
0056       lfact *= l;
0057       Prob += exp(l * Loginvlog - log(lfact));
0058     }
0059
0060     float LogProb = log(Prob);
0061     float ProbJet = std::min((float)exp(std::max(LogProb + SumJet, -30.0f)), 1.0f);
0062     float TotalProb = -log10(ProbJet) / 4;
0063     TotalProb = 1 - TotalProb;
0064
0065     return std::make_pair(TotalProb, -1);
0066   }
0067
0068 private:
0069   float meVperADCStrip;
0070   deDxTools::ESGetTokenH3DDVariant token_;
0071   std::string ProbabilityMode;
0072   TH3F* Prob_ChargePath;
0073 };
0074
0075 #endif