DeDx/interface/ASmirnovDeDxDiscriminator.h

0001 #ifndef RecoTrackerDeDx_ASmirnovDeDxDiscriminator_h
0002 #define RecoTrackerDeDx_ASmirnovDeDxDiscriminator_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 ASmirnovDeDxDiscriminator : public BaseDeDxEstimator {
0009 public:
0010   ASmirnovDeDxDiscriminator(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 TotalProb = 1.0 / (12 * size);
0045     for (size_t i = 1; i <= size; i++) {
0046       TotalProb += vect_probs[i - 1] * pow(vect_probs[i - 1] - ((2.0 * i - 1.0) / (2.0 * size)), 2);
0047     }
0048     TotalProb *= (3.0 / size);
0049     return std::make_pair(TotalProb, -1);
0050   }
0051
0052 private:
0053   float meVperADCStrip;
0054   deDxTools::ESGetTokenH3DDVariant token_;
0055   std::string ProbabilityMode;
0056   TH3F* Prob_ChargePath;
0057 };
0058
0059 #endif