base classes

bb/CMakeLists.txt

-include_directories(${CMAKE_CURRENT_SOURCE_DIR})
+include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_SOURCE_DIR})
 ROOT_GENERATE_DICTIONARY(G__PidBB
         BetheBlochHelper.h Shine.h
         LINKDEF BBLinkDef.h

src/CMakeLists.txt

 
 set(SOURCES
-#    ParticleFit.cpp
     Getter.cpp
     CutGGetter.cpp
     ParticlePdf1D.cpp

src/Constants.h

@@ -38,6 +38,8 @@ const std::unordered_map<int, float> masses =
     };
 }
 
+typedef int PdgCode_t;
+
 namespace PidFunction {
 enum eNames {
   kA = 0,

src/ParticlePdf1D.h

@@ -2,6 +2,7 @@
 #define PID_SRC_PARTICLEPDF1D_H_
 
 #include "TF1.h"
+#include "Constants.h"
 
 namespace Pid {
 
@@ -10,7 +11,7 @@ class ParticlePdf1DBase {
   ParticlePdf1DBase() = default;
   virtual ~ParticlePdf1DBase() = default;
 
-  virtual float GetPdfValue(float y) = 0;
+  virtual float GetPdfValue(float y) const = 0;
 
  protected:
   bool is_init_{false};
@@ -19,14 +20,22 @@ class ParticlePdf1DBase {
 class ParticlePdf1D : public ParticlePdf1DBase {
  public:
   ParticlePdf1D() = default;
-  ParticlePdf1D(const TF1& pdf) : pdf_(pdf), is_init_(true) {}
+  ParticlePdf1D(const TF1& pdf, const std::vector<double>& par, PdgCode_t pdg) : pdf_(pdf), pdg_(pdg), is_init_(true) {
+    assert(pdf_.GetNpar() == par.size());
+    pdf_.SetParameters(&par[0]);
+  }
 
-  float GetPdfValue(float y) final {
+  float GetPdfValue(float y) const final {
     return pdf_.Eval(y);
   }
 
+  PdgCode_t GetPdg() const {
+    return pdg_;
+  }
+
  protected:
   TF1 pdf_;
+  PdgCode_t pdg_{-1};
   bool is_init_{false};
 };
 }

src/ParticlePdf2D.h

@@ -13,7 +13,7 @@ class ParticlePdf2DBase{
   virtual ~ParticlePdf2DBase() = default;
 
   virtual float GetPdfValue(float x, float y) const = 0;
-  virtual ParticlePdf1DBase* GetPdf1D(float x) const = 0;
+  virtual const ParticlePdf1DBase* GetPdf1D(float x) const = 0;
 
  protected:
   float min_x_{0.f};
@@ -31,11 +31,11 @@ class ParticlePdf2DBinned : public ParticlePdf2DBase {
   ~ParticlePdf2DBinned() override = default;
 
   float GetPdfValue(float x, float y) const override {
-    const int bin = FindBin(x);
-    return bin>=0 && bin<pdfs_.size() ? pdfs_.at(FindBin(x))->GetPdfValue(y) : 0.f;
+    const auto* pdf = GetPdf1D(x);
+    return pdf != nullptr ? pdf->GetPdfValue(y) : 0.f;
   }
 
-  ParticlePdf1DBase* GetPdf1D(float x) const override {
+  const ParticlePdf1DBase* GetPdf1D(float x) const override {
     const int bin = FindBin(x);
     return bin>=0 && bin<pdfs_.size() ? pdfs_.at(FindBin(x)) : nullptr;
   }

src/PidLinkDef.h

@@ -4,10 +4,15 @@
 #pragma link off all classes;
 #pragma link off all functions;
 
-#pragma link C++ class Pid::BaseFitter+;
-#pragma link C++ class Pid::ParticleFit+;
+#pragma add namespace Pid;
+
 #pragma link C++ class Pid::BaseGetter--;
 #pragma link C++ class Pid::Getter+;
 #pragma link C++ class Pid::CutGGetter+;
 
+#pragma link C++ class Pid::ParticlePdf1DBase+;
+#pragma link C++ class Pid::ParticlePdf1D+;
+#pragma link C++ class Pid::ParticlePdf2DBase+;
+#pragma link C++ class Pid::ParticlePdf2DBinned+;
+
 #endif

tasks/fit_tof.cpp

@@ -14,7 +14,7 @@
 
 using namespace RooFit;
 
-std::vector<Pid::ParticlePdf1D*> Fit1D(TH1* input, float momentum);
+std::vector<Pid::ParticlePdf1DBase*> Fit1D(TH1* input, float momentum);
 
 int main(int argc, char* argv[]){
 
@@ -34,26 +34,44 @@ int main(int argc, char* argv[]){
 
   auto* out_file{TFile::Open("result.root", "recreate")};
 
-  for(int i=0; i<m2p->GetXaxis()->GetNbins(); ++i)
+  std::vector<Pid::ParticlePdf1DBase*> pdfs_pion;
+  std::vector<Pid::ParticlePdf1DBase*> pdfs_kaon;
+  std::vector<Pid::ParticlePdf1DBase*> pdfs_proton;
+
+  std::vector<double> bins{ m2p->GetXaxis()->GetBinLowEdge(1) };
+  for(int i=1; i<m2p->GetXaxis()->GetNbins(); ++i)
   {
     TH1* slice = m2p->ProjectionY(Form("px_%d",i), i, i);
 
     if(slice->GetEntries() < 1000) continue;
+    bins.emplace_back(m2p->GetXaxis()->GetBinUpEdge(i));
+    std::cout << "bins " << bins.back() << std::endl;
+
     const float momentum = m2p->GetXaxis()->GetBinCenter(i);
 
     auto result = Fit1D(slice, momentum);
+    pdfs_pion.emplace_back( result[0] );
+    pdfs_kaon.emplace_back( result[1] );
+    pdfs_proton.emplace_back( result[2] );
+  }
+  TAxis a(static_cast<Int_t>(bins.size()-1), &(bins[0]));
 
+  Pid::ParticlePdf2DBinned* protons = new Pid::ParticlePdf2DBinned(pdfs_proton, {static_cast<Int_t>(bins.size()-1), &(bins[0])});
+  Pid::ParticlePdf2DBinned* kaons = new Pid::ParticlePdf2DBinned(pdfs_kaon, {static_cast<Int_t>(bins.size()-1), &(bins[0])});
+  Pid::ParticlePdf2DBinned* pions = new Pid::ParticlePdf2DBinned(pdfs_pion, {static_cast<Int_t>(bins.size()-1), &(bins[0])});
 
-  }
+  Pid::Getter tof;
+  tof.AddParticle(protons, 2212);
+  tof.AddParticle(kaons, 321);
+  tof.AddParticle(pions, 211);
 
-  Pid::ParticlePdf2DBinned p;
-  Pid::Getter tof();
+  tof.Write("tof_getter");
 
   out_file->Close();
 
 }
 
-std::vector<Pid::ParticlePdf1D*> Fit1D(TH1* input, float momentum){
+std::vector<Pid::ParticlePdf1DBase*> Fit1D(TH1* input, float momentum){
 
   const float p2 = momentum*momentum;
 
@@ -96,5 +114,12 @@ std::vector<Pid::ParticlePdf1D*> Fit1D(TH1* input, float momentum){
 
   m2frame->Write(Form("bin_%f", momentum));
 
-  return std::vector<Pid::ParticlePdf1D*>();
+  TF1 gaus("g", "gaus", -0.5, 1.5);
+  std::vector<Pid::ParticlePdf1DBase*> ret = {
+    new Pid::ParticlePdf1D(gaus, {pi_frac.getVal(), pi_mean.getVal(), pi_width.getVal()}, 211),
+    new Pid::ParticlePdf1D(gaus, {K_frac.getVal(), K_mean.getVal(), K_width.getVal()}, 321),
+    new Pid::ParticlePdf1D(gaus, { 1-pi_frac.getVal()-K_frac.getVal(), p_mean.getVal(), p_width.getVal()}, 2212),
+  };
+
+  return ret;
 }
\ No newline at end of file