KaliVeda
Toolkit for HIC analysis
bayesian_estimator.h
1 
4 #ifndef __bayesian_estimator_H
5 #define __bayesian_estimator_H
6 
7 #include "KVBase.h"
8 #include "impact_parameter_distribution.h"
9 
10 #include <vector>
11 #include <numeric>
12 #include <TF2.h>
13 #include <TH1.h>
14 #include <TCanvas.h>
15 #include <TGraph.h>
16 #include <KVValueRange.h>
17 #include <KVNameValueList.h>
18 #include <TGraphErrors.h>
19 #include <TH2.h>
20 #include "TMath.h"
21 #include "TNamed.h"
22 
23 #include "Math/PdfFuncMathCore.h"
24 
25 namespace KVImpactParameters {
41  class gamma_kernel {
42  public:
43  double operator()(double X, double mean, double reduced_variance)
44  {
46  return ROOT::Math::gamma_pdf(X, mean / reduced_variance, reduced_variance);
47  }
48 
50  };
51 
67  class BD_kernel {
68  public:
69  double operator()(double X, double mean, double reduced_variance)
70  {
72  return ROOT::Math::binomial_pdf(TMath::Nint(X), 1. - reduced_variance, mean / (1. - reduced_variance));
73  }
74 
76  };
77 
93  class NBD_kernel {
94  public:
95  double operator()(double X, double mean, double reduced_variance)
96  {
98  return ROOT::Math::negative_binomial_pdf(TMath::Nint(X), 1. / reduced_variance, mean / (reduced_variance - 1.));
99  }
100 
101  ClassDef(NBD_kernel, 0)
102  };
103 
104 
219  template
220  <class FittingFunction, class FluctuationKernel>
221  class bayesian_estimator : public KVBase {
222 
223  FittingFunction theFitter;
224  FluctuationKernel theKernel;
225 
228  std::vector<double> sel_rapp;
241 
243 
244  impact_parameter_distribution fIPDist;// used to convert centrality to impact parameter
245 
246  double mean_X_vs_cb(double* x, double* par)
247  {
252 
253  theFitter.fill_params_from_array(par);
254  return theFitter.meanX(x[0]);
255  }
256  double mean_X_vs_b(double* x, double* par)
257  {
262 
263  theFitter.fill_params_from_array(par);
264  return theFitter.meanX(fIPDist.GetCentrality().Eval(x[0]));
265  }
266  double P_X_cb(double X, double cb)
267  {
271  return fIntegerVariable ? theKernel(TMath::Nint(X), theFitter.meanX(cb), theFitter.redVar(cb))
272  : theKernel(X, theFitter.meanX(cb), theFitter.redVar(cb));
273  }
274  double P_X_cb_for_TF2_obs_vs_b(double* x, double*)
275  {
279 
280  auto X = x[1];
281  auto cb = fIPDist.GetCentrality().Eval(x[0]);
282  return x[0] * theKernel(X, theFitter.meanX(cb), theFitter.redVar(cb));
283  }
284  double P_X_cb_for_integral(double* x, double* par)
285  {
292  return P_X_cb(par[0], x[0]);
293  }
294  double P_X_cb_for_X_integral(double* x, double* par)
295  {
303  return P_X_cb(x[0], par[0]);
304  }
305  double P_X_cb_for_X_integral_with_selection(double* x, double* par)
306  {
314 
315  int bin = histo->FindBin(x[0]);
316  if (bin < 1 || bin > histo->GetNbinsX()) return 0;
317  return sel_rapp[bin - 1] * P_X_cb(x[0], par[0]);
318  }
319  double cb_integrated_P_X(double* x, double* p)
320  {
328 
329  theFitter.fill_params_from_array(p);
330  p_X_cb_integrator.SetParameter(0, x[0]); // set value of X in integrand
331  return p_X_cb_integrator.Integral(0, 1, 1.e-4);
332  }
333  double P_X_from_fit(double* x, double* par)
334  {
339  p_X_cb_integrator.SetParameter(0, x[0]); // set value of X in integrand
340  return par[0] * p_X_cb_integrator.Integral(0, 1, 1.e-4);
341  }
342  double cb_dist_for_X_selection(double* x, double* p)
343  {
353 
355  double num = p_X_X_integrator.Integral(p[0], p[1], 1.e-4);
356  double den = fitted_P_X.Integral(p[0], p[1], 1.e-4);
357  if (den > 0) return num / den;
358  return 0;
359  }
360  double cb_dist_for_arb_X_selection(double* x, double* p)
361  {
370 
372  double num = p_X_X_integrator_with_selection.Integral(p[0], p[1], 1.e-4);
374 
375  return num;
376  }
377  double b_dist_for_X_selection(double* x, double* p)
378  {
388 
390  double num = p_X_X_integrator.Integral(p[0], p[1], 1.e-4);
391  return num * fIPDist.GetDifferentialCrossSection(x[0]);
392  }
393  double b_dist_for_arb_X_selection(double* x, double* p)
394  {
403 
405  double num = p_X_X_integrator_with_selection.Integral(p[0], p[1], 1.e-4);
406  return num * fIPDist.GetDifferentialCrossSection(x[0]);
407  }
408  double b_dist_for_arb_X_selection_from_histo(double* x, double* p)
409  {
418  TH1* hh_centb = fIPDist.GetCentralityFromHisto();
419  Double_t centb = hh_centb->GetBinContent(hh_centb->FindBin(x[0]));
420 
422  double num = p_X_X_integrator_with_selection.Integral(p[0], p[1], 1.e-4);
424  }
425 
426  public:
427  bayesian_estimator(Bool_t integer_variable = false)
428  : KVBase(),
429  theFitter(),
430  p_X_cb_integrator("p_X_cb_integrator", this, &bayesian_estimator::P_X_cb_for_integral, 0, 1, 1),
431  P_X_fit_function("P_X_fit_function", this, &bayesian_estimator::cb_integrated_P_X, 0, 1, theFitter.npar()),
432  mean_X_vs_cb_function("mean_X_vs_cb", this, &bayesian_estimator::mean_X_vs_cb, 0, 1, theFitter.npar()),
433  mean_X_vs_b_function("mean_X_vs_b", this, &bayesian_estimator::mean_X_vs_b, 0, 20, theFitter.npar()),
434  p_X_X_integrator("p_X_X_integrator", this, &bayesian_estimator::P_X_cb_for_X_integral, 0, 1000, 1),
435  p_X_X_integrator_with_selection("p_X_X_integrator_with_selection", this, &bayesian_estimator::P_X_cb_for_X_integral_with_selection, 0, 1000, 1),
436  fitted_P_X("fitted_P_X", this, &bayesian_estimator::P_X_from_fit, 0, 1000, 1),
437  Cb_dist_for_X_select("Cb_dist_for_X_select", this, &bayesian_estimator::cb_dist_for_X_selection, 0, 1, 2),
438  Cb_dist_for_arb_X_select("Cb_dist_for_arb_X_select", this, &bayesian_estimator::cb_dist_for_arb_X_selection, 0, 1, 2),
439  B_dist_for_X_select("b_dist_for_X_select", this, &bayesian_estimator::b_dist_for_X_selection, 0, 20, 2),
440  B_dist_for_arb_X_select("b_dist_for_arb_X_select", this, &bayesian_estimator::b_dist_for_arb_X_selection, 0, 20, 2),
441  B_dist_for_arb_X_select_from_histo("b_dist_for_arb_X_select_from_histo", this, &bayesian_estimator::b_dist_for_arb_X_selection_from_histo, 0, 20, 2),
442  fIntegerVariable(integer_variable)
443  {
448  theFitter.set_par_names(P_X_fit_function);
449  theFitter.set_par_names(mean_X_vs_cb_function);
450  theFitter.set_par_names(mean_X_vs_b_function);
451  }
452 
453  bayesian_estimator(const FittingFunction& previous_fit, Bool_t integer_variable = false)
454  : KVBase(),
455  theFitter(previous_fit),
456  p_X_cb_integrator("p_X_cb_integrator", this, &bayesian_estimator::P_X_cb_for_integral, 0, 1, 1),
457  P_X_fit_function("P_X_fit_function", this, &bayesian_estimator::cb_integrated_P_X, 0, 1, theFitter.npar()),
458  mean_X_vs_cb_function("mean_X_vs_cb", this, &bayesian_estimator::mean_X_vs_cb, 0, 1, theFitter.npar()),
459  mean_X_vs_b_function("mean_X_vs_b", this, &bayesian_estimator::mean_X_vs_b, 0, 20, theFitter.npar()),
460  p_X_X_integrator("p_X_X_integrator", this, &bayesian_estimator::P_X_cb_for_X_integral, 0, 1000, 1),
461  p_X_X_integrator_with_selection("p_X_X_integrator_with_selection", this, &bayesian_estimator::P_X_cb_for_X_integral_with_selection, 0, 1000, 1),
462  fitted_P_X("fitted_P_X", this, &bayesian_estimator::P_X_from_fit, 0, 1000, 1),
463  Cb_dist_for_X_select("Cb_dist_for_X_select", this, &bayesian_estimator::cb_dist_for_X_selection, 0, 1, 2),
464  Cb_dist_for_arb_X_select("Cb_dist_for_arb_X_select", this, &bayesian_estimator::cb_dist_for_arb_X_selection, 0, 1, 2),
465  B_dist_for_X_select("b_dist_for_X_select", this, &bayesian_estimator::b_dist_for_X_selection, 0, 20, 2),
466  B_dist_for_arb_X_select("b_dist_for_arb_X_select", this, &bayesian_estimator::b_dist_for_arb_X_selection, 0, 20, 2),
467  B_dist_for_arb_X_select_from_histo("b_dist_for_arb_X_select_from_histo", this, &bayesian_estimator::b_dist_for_arb_X_selection_from_histo, 0, 20, 2),
468  fIntegerVariable(integer_variable)
469  {
475  theFitter.set_par_names(P_X_fit_function);
476  theFitter.set_par_names(mean_X_vs_cb_function);
477  theFitter.set_par_names(mean_X_vs_b_function);
478  }
479 
480  virtual ~bayesian_estimator() {}
481 
483  {
489 
490  histo = h;
491  histo->Scale(1. / h->Integral("width"));
492  }
493 
494  void FitHisto(TH1* h = nullptr)
495  {
506 
507  if (h) histo = h;
508  auto xmin = histo->GetXaxis()->GetBinLowEdge(1);
511  theFitter.set_initial_parameters(histo, P_X_fit_function);
512  histo->Fit(&P_X_fit_function);//,"EMG");
513  }
514 
515  void SetIPDistParams(double sigmaR, double deltab)
516  {
533 
535  }
536 
537  void SetIPDistFromHisto(TH1* ip_histo)
538  {
552 
553  fIPDist.FitIPDist(ip_histo);
554  }
555 
557  {
564  return fIPDist;
565  }
566 
567  void DrawMeanXvsCb(const TString& title = "", Color_t color = -1, Option_t* opt = "")
568  {
577 
578  GetMeanXvsCb().Draw();
579  }
581  {
588 
589  Double_t par[theFitter.npar()];
590  theFitter.fill_array_from_params(par);
592  return mean_X_vs_cb_function;
593  }
594  TGraph* GraphMeanXvsb(int npts = 500)
595  {
605 
606  Double_t par[theFitter.npar()];
607  theFitter.fill_array_from_params(par);
608  KVValueRange<double> brange(0, GetIPDist().GetB0() + 4 * GetIPDist().GetDeltaB());
610  auto gr = new TGraph;
611  for (int i = 0; i <= npts; ++i) gr->SetPoint(i, brange.ValueIofN(i, npts), mean_X_vs_b_function.Eval(brange.ValueIofN(i, npts)));
612  return gr;
613  }
614  TH2* TransformXaxisToBaxis(const TH2* obs_vs_X, int baxis_nbins = 100, double baxis_bmin = 0, double baxis_bmax = 1)
615  {
620 
621  auto Nx = obs_vs_X->GetNbinsX();
622  auto Ny = obs_vs_X->GetNbinsY();
623  auto ymin = obs_vs_X->GetYaxis()->GetBinLowEdge(1);
624  auto ymax = obs_vs_X->GetYaxis()->GetBinUpEdge(Ny);
625  auto ywidth = obs_vs_X->GetYaxis()->GetBinWidth(1);
626  auto obs_vs_b = new TH2F("obs_vs_b", obs_vs_X->GetTitle(),
627  baxis_nbins, baxis_bmin, baxis_bmax, Ny, ymin, ymax);
628  for (int binx = 1; binx <= Nx; ++binx) {
630  obs_vs_X->GetXaxis()->GetBinUpEdge(binx));
631  for (int biny = 1; biny <= Ny; ++biny) {
632  auto stat = obs_vs_X->GetBinContent(binx, biny);
633  auto yval = obs_vs_X->GetYaxis()->GetBinCenter(biny);
634  for (int i = 0; i < stat; ++i) {
635  obs_vs_b->Fill(B_dist_for_X_select.GetRandom(0., 1.),
636  gRandom->Uniform(yval - .5 * ywidth, yval + .5 * ywidth));
637  }
638  }
639  }
640  return obs_vs_b;
641  }
643  {
645  if (!histo) {
646  Error("GraphMeanbvsX", "Need to give histogram of X observable to NormalizeHisto() first");
647  return nullptr;
648  }
649 
650  auto g = new TGraphErrors;
651  auto N = histo->GetNbinsX();
652  auto W = histo->GetBinWidth(1);
653 
654  std::cout << N << std::endl;
655 
656  for (int i = 1; i <= N; ++i) {
657  double mean, sigma;
658  if (histo->GetBinContent(i) > 0) {
659  auto X = histo->GetBinCenter(i);
660  std::cout << i << " " << X;
661  GetMeanAndSigmaBDistForXSelection({X - W / 2, X + W / 2}, mean, sigma);
662 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,24,0)
663  g->AddPoint(X, mean);
664 #else
665  g->SetPoint(g->GetN(), X, mean);
666 #endif
667  g->SetPointError(g->GetN() - 1, 0, sigma);
668  std::cout << " " << mean << " " << sigma << std::endl;
669  }
670  }
671  return g;
672  }
674  {
679 
680  if (!histo) {
681  Warning("update_fit_params", "no histogram set with FitHisto(TH1*)");
682  return;
683  }
684  TF1* fit = (TF1*)histo->FindObject("P_X_fit_function");
685  if (!fit) {
686  Warning("update_fit_params", "no fit function found in histogram");
687  return;
688  }
689  theFitter.fill_params_from_array(fit->GetParameters());
690  }
691  double DrawCbDistForXSelection(double X1, double X2, Option_t* opt = "", Color_t color = kRed, const TString& title = "")
692  {
707  f->SetNpx(500);
708  f->SetLineColor(color);
709  f->SetMarkerColor(color);
710  f->SetLineWidth(2);
711  f->SetTitle(title);
712  return f->GetMaximum();
713  }
714  void DrawCbDistForSelection(TH1* sel, TH1* incl, double& mean_cb, double& sigma_cb, Option_t* opt = "", Color_t color = kRed, const TString& title = "")
715  {
730 
731  assert(sel->GetNbinsX() == incl->GetNbinsX());
732 
735  sel_rapp.assign((std::vector<double>::size_type)incl->GetNbinsX(), 0.0);
736  int first_bin(0), last_bin(0);
737  for (int i = 1; i <= incl->GetNbinsX(); ++i) {
738  if (incl->GetBinContent(i) > 0) {
739  sel_rapp[i - 1] = sel->GetBinContent(i) / incl->GetBinContent(i);
740  if (sel->GetBinContent(i) > 0) {
741  if (!first_bin) first_bin = i;
742  last_bin = i;
743  }
744  }
745  }
746  double Xmin = incl->GetBinLowEdge(first_bin);
747  double Xmax = incl->GetXaxis()->GetBinUpEdge(last_bin);
748 
749  histo = incl;
750  h_selection = sel;
751 
756 
758  double cb_mean(0), cb_sqrmean(0), sum_pcb(0);
759  TGraph* f = new TGraph;
760  for (int i = 0; i < 500; ++i) {//graph with 500 points
761  double cb = i / 499.;
762  double p_cb = Cb_dist_for_arb_X_select.Eval(cb);
763 
764  cb_mean += p_cb * cb;
765  cb_sqrmean += p_cb * cb * cb;
766  sum_pcb += p_cb;
767 
768  f->SetPoint(i, cb, p_cb);
769  }
770 
772  cb_sqrmean /= sum_pcb;
773  cb_mean /= sum_pcb;
774  mean_cb = cb_mean;
775  sigma_cb = TMath::Sqrt(cb_sqrmean - cb_mean * cb_mean);
776 
778  f->SetLineColor(color);
779  f->SetMarkerColor(color);
780  f->SetLineWidth(2);
781  f->SetTitle(title);
782  if (TString(opt) == "same") f->Draw("l");
783  else f->Draw("al");
784  }
785  double DrawBDistForXSelection(KVValueRange<double> Xrange, Option_t* opt = "", Color_t color = kRed, const TString& title = "")
786  {
803 
804  B_dist_for_X_select.SetParameters(Xrange.Min(), Xrange.Max());
805  TGraph* f = new TGraph;
806  double maxS = 0;
807  for (int i = 0; i < 500; ++i) {
808  double b = 2 * i * GetIPDist().GetB0() / 499.;
809  double sig = B_dist_for_X_select.Eval(b);
810  if (sig > maxS) maxS = sig;
811  f->SetPoint(i, b, sig);
812  }
813  f->SetLineColor(color);
814  f->SetMarkerColor(color);
815  f->SetLineWidth(2);
816  f->SetTitle(title);
817  if (TString(opt) == "same") f->Draw("l");
818  else f->Draw("al");
819  return maxS;
820  }
822  {
823  return B_dist_for_X_select;
824  }
826  {
839 
840 
841  B_dist_for_X_select.SetParameters(Xrange.Min(), Xrange.Max());
842  TGraph* f = new TGraph;
843  for (int i = 0; i < npts; ++i) {
844  double b = 2 * i * GetIPDist().GetB0() / 499.;
845  double sig = B_dist_for_X_select.Eval(b);
846  f->SetPoint(i, b, sig);
847  }
848  f->SetLineWidth(2);
849  return f;
850  }
851  void DrawBDistForSelection(TH1* sel, TH1* incl, double& mean, double& sigma, Option_t* opt = "", Color_t color = kRed, const TString& title = "")
852  {
872 
873 
874  assert(sel->GetNbinsX() == incl->GetNbinsX());
875 
876  sel_rapp.assign((std::vector<double>::size_type)incl->GetNbinsX(), 0.0);
877  int first_bin(0), last_bin(0);
878  for (int i = 1; i <= incl->GetNbinsX(); ++i) {
879  if (incl->GetBinContent(i) > 0) {
880  sel_rapp[i - 1] = sel->GetBinContent(i) / incl->GetBinContent(i);
881  if (sel->GetBinContent(i) > 0) {
882  if (!first_bin) first_bin = i;
883  last_bin = i;
884  }
885  }
886  }
887  double Xmin = incl->GetBinLowEdge(first_bin);
888  double Xmax = incl->GetXaxis()->GetBinUpEdge(last_bin);
889 
890  histo = incl;
891  h_selection = sel;
892 
894  B_dist_for_arb_X_select.SetRange(0, 2 * GetIPDist().GetB0());
896 
897  double bmean(0), bsqrmean(0), sigtot(0);
898  TGraph* f = new TGraph;
899  for (int i = 0; i < 500; ++i) {
900  double b = 2 * i * GetIPDist().GetB0() / 499.;
901  double sig = B_dist_for_arb_X_select.Eval(b);
902  bmean += sig * b;
903  bsqrmean += sig * b * b;
904  sigtot += sig;
905  f->SetPoint(i, b, sig);
906  }
907  f->SetLineColor(color);
908  f->SetMarkerColor(color);
909  f->SetLineWidth(2);
910  f->SetTitle(title);
911  mean = bmean / sigtot;
912  bsqrmean /= sigtot;
913  sigma = TMath::Sqrt(bsqrmean - mean * mean);
914  if (TString(opt) == "same") f->Draw("l");
915  else f->Draw("al");
916  }
917 
918  void DrawBDistForSelectionFromIPHisto(TH1* sel, TH1* incl, double& mean, double& sigma, Option_t* opt = "", Color_t color = kRed, const TString& title = "")
919  {
937 
938  assert(sel->GetNbinsX() == incl->GetNbinsX());
939 
940  sel_rapp.assign((std::vector<double>::size_type)incl->GetNbinsX(), 0.0);
941  int first_bin(0), last_bin(0);
942  for (int i = 1; i <= incl->GetNbinsX(); ++i) {
943  if (incl->GetBinContent(i) > 0) {
944  sel_rapp[i - 1] = sel->GetBinContent(i) / incl->GetBinContent(i);
945  if (sel->GetBinContent(i) > 0) {
946  if (!first_bin) first_bin = i;
947  last_bin = i;
948  }
949  }
950  }
951  double Xmin = incl->GetBinLowEdge(first_bin);
952  double Xmax = incl->GetXaxis()->GetBinUpEdge(last_bin);
953 
954  histo = incl;
955  h_selection = sel;
956 
958  B_dist_for_arb_X_select_from_histo.SetRange(0, 2 * GetIPDist().GetB0()); //if SetIPDistFromHisto() was called, a fit have been applied thus B0 was estimated
960 
961  double bmean(0), bsqrmean(0), sigtot(0);
962  TGraph* f = new TGraph;
963  for (int i = 0; i < 500; ++i) {
964  double b = 2 * i * GetIPDist().GetB0() / 499.;
966  bmean += sig * b;
967  bsqrmean += sig * b * b;
968  sigtot += sig;
969  f->SetPoint(i, b, sig);
970  }
971  f->SetLineColor(color);
972  f->SetMarkerColor(color);
973  f->SetLineWidth(2);
974  f->SetTitle(title);
975  mean = bmean / sigtot;
976  bsqrmean /= sigtot;
977  sigma = TMath::Sqrt(bsqrmean - mean * mean);
978  if (TString(opt) == "same") f->Draw("l");
979  else f->Draw("al");
980  }
981 
982  void GetMeanAndSigmaBDistForSelection(TH1* sel, TH1* incl, double& mean, double& sigma)
983  {
1000 
1001 
1002  assert(sel->GetNbinsX() == incl->GetNbinsX());
1003 
1004  sel_rapp.assign((std::vector<double>::size_type)incl->GetNbinsX(), 0.0);
1005  int first_bin(0), last_bin(0);
1006  for (int i = 1; i <= incl->GetNbinsX(); ++i) {
1007  if (incl->GetBinContent(i) > 0) {
1008  sel_rapp[i - 1] = sel->GetBinContent(i) / incl->GetBinContent(i);
1009  if (sel->GetBinContent(i) > 0) {
1010  if (!first_bin) first_bin = i;
1011  last_bin = i;
1012  }
1013  }
1014  }
1015  double Xmin = incl->GetBinLowEdge(first_bin);
1016  double Xmax = incl->GetXaxis()->GetBinUpEdge(last_bin);
1017 
1018  histo = incl;
1019  h_selection = sel;
1020 
1022  mean = B_dist_for_arb_X_select.Mean(0, 20);
1023  double var = B_dist_for_arb_X_select.Variance(0, 20);
1024  sigma = TMath::Sqrt(var);
1025  }
1026 
1027  void GetMeanAndSigmaBDistForXSelection(KVValueRange<double> Xrange, double& mean, double& sigma)
1028  {
1043 
1044  B_dist_for_X_select.SetParameters(Xrange.Min(), Xrange.Max());
1045  mean = sigma = 0;
1046  double sig_tot = 0;
1047  for (int i = 0; i < 500; ++i) {
1048  double b = 2 * i * GetIPDist().GetB0() / 499.;
1049  double sig = B_dist_for_X_select.Eval(b);
1050  mean += b * sig;
1051  sigma += pow(b, 2) * sig;
1052  sig_tot += sig;
1053  }
1054  if (sig_tot <= 0) {
1055  Warning("GetMeanAndSigmaBDistForXSelection", "Total cross-section = 0!");
1056  return;
1057  }
1058  mean /= sig_tot;
1059  sigma /= sig_tot;
1060  sigma = pow(sigma - pow(mean, 2), 0.5);
1061  }
1062  void DrawFittedP_X(double norm = 1.0, Option_t* opt = "", Color_t color = kRed, const TString& title = "")
1063  {
1070 
1071  TF1* f = GetFittedP_X(norm)->DrawCopy(opt);
1072  f->SetNpx(500);
1073  f->SetLineColor(color);
1074  f->SetLineWidth(2);
1075  f->SetTitle(title);
1076  }
1077  TF1* GetFittedP_X(double norm = 1.0)
1078  {
1083  fitted_P_X.SetParameter(0, norm);
1084  return &fitted_P_X;
1085  }
1086  TGraph* GraphP_XForGivenB(double b, KVValueRange<double> Xrange, int npts = 500)
1087  {
1094  auto graph = new TGraph;
1095  for (int i = 0; i < npts; ++i) graph->SetPoint(i, Xrange.ValueIofN(i, npts), p_X_X_integrator.Eval(Xrange.ValueIofN(i, npts)));
1096  return graph;
1097  }
1098  void Print(Option_t* = "") const
1099  {
1101  theFitter.print_fit_params();
1102  }
1104  {
1106  return theFitter.get_params();
1107  }
1108  void DrawNormalisedMeanXvsb(const TString& title, Color_t color, Option_t* opt)
1109  {
1125 
1126 
1127  theFitter.backup_params();
1128  theFitter.normalise_shape_function();
1129  Double_t par[5];
1130  theFitter.fill_array_from_params(par);
1131  mean_X_vs_b_function.SetRange(0, GetIPDist().GetB0() + 2 * GetIPDist().GetDeltaB());
1133  TF1* copy = mean_X_vs_b_function.DrawCopy(opt);
1134  if (!title.IsNull()) copy->SetTitle(title);
1135  if (color >= 0) copy->SetLineColor(color);
1136  theFitter.restore_params();
1137  }
1138 
1140  {
1143 
1144  auto jpd = new TF2("joint_proba_dist", this, &bayesian_estimator::P_X_cb_for_TF2_obs_vs_b,
1145  b_range.Min(), b_range.Max(), X_range.Min(), X_range.Max(), 0);
1146  return jpd;
1147  }
1148 
1149  ClassDef(bayesian_estimator, 1) //Estimate impact parameter distribution by fits to data
1150  };
1151 
1152 }
1153 
1154 #endif
#define f(i)
bool Bool_t
short Color_t
double Double_t
const char Option_t
#define ClassDef(name, id)
#define X(type, name)
#define N
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t sel
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t g
float xmin
float ymin
float xmax
float ymax
R__EXTERN TRandom * gRandom
Base class for KaliVeda framework.
Definition: KVBase.h:139
void Error(const char *method, const char *msgfmt,...) const override
Definition: KVBase.cpp:1650
void Warning(const char *method, const char *msgfmt,...) const override
Definition: KVBase.cpp:1637
Fluctuation kernel using binomial distribution for use with bayesian_estimator.
double operator()(double X, double mean, double reduced_variance)
Fluctuation kernel using negative binomial distribution for use with bayesian_estimator.
double operator()(double X, double mean, double reduced_variance)
Impact parameter distribution reconstruction from experimental data.
double P_X_cb_for_integral(double *x, double *par)
impact_parameter_distribution fIPDist
void DrawBDistForSelection(TH1 *sel, TH1 *incl, double &mean, double &sigma, Option_t *opt="", Color_t color=kRed, const TString &title="")
void DrawCbDistForSelection(TH1 *sel, TH1 *incl, double &mean_cb, double &sigma_cb, Option_t *opt="", Color_t color=kRed, const TString &title="")
double mean_X_vs_b(double *x, double *par)
double DrawBDistForXSelection(KVValueRange< double > Xrange, Option_t *opt="", Color_t color=kRed, const TString &title="")
void GetMeanAndSigmaBDistForSelection(TH1 *sel, TH1 *incl, double &mean, double &sigma)
bayesian_estimator(Bool_t integer_variable=false)
double cb_dist_for_X_selection(double *x, double *p)
void SetIPDistParams(double sigmaR, double deltab)
double b_dist_for_arb_X_selection_from_histo(double *x, double *p)
double P_X_cb_for_TF2_obs_vs_b(double *x, double *)
double P_X_from_fit(double *x, double *par)
double b_dist_for_arb_X_selection(double *x, double *p)
double DrawCbDistForXSelection(double X1, double X2, Option_t *opt="", Color_t color=kRed, const TString &title="")
bayesian_estimator(const FittingFunction &previous_fit, Bool_t integer_variable=false)
double b_dist_for_X_selection(double *x, double *p)
impact_parameter_distribution & GetIPDist()
void GetMeanAndSigmaBDistForXSelection(KVValueRange< double > Xrange, double &mean, double &sigma)
double cb_integrated_P_X(double *x, double *p)
double P_X_cb_for_X_integral(double *x, double *par)
void DrawNormalisedMeanXvsb(const TString &title, Color_t color, Option_t *opt)
double cb_dist_for_arb_X_selection(double *x, double *p)
double mean_X_vs_cb(double *x, double *par)
TGraph * GraphBDistForXSelection(KVValueRange< double > Xrange, int npts=500)
void DrawBDistForSelectionFromIPHisto(TH1 *sel, TH1 *incl, double &mean, double &sigma, Option_t *opt="", Color_t color=kRed, const TString &title="")
double P_X_cb_for_X_integral_with_selection(double *x, double *par)
TGraph * GraphP_XForGivenB(double b, KVValueRange< double > Xrange, int npts=500)
TF2 * GetJointProbabilityDistribution(KVValueRange< double > b_range, KVValueRange< double > X_range)
void DrawMeanXvsCb(const TString &title="", Color_t color=-1, Option_t *opt="")
TH2 * TransformXaxisToBaxis(const TH2 *obs_vs_X, int baxis_nbins=100, double baxis_bmin=0, double baxis_bmax=1)
void DrawFittedP_X(double norm=1.0, Option_t *opt="", Color_t color=kRed, const TString &title="")
Fluctuation kernel using gamma distribution for use with bayesian_estimator.
double operator()(double X, double mean, double reduced_variance)
Class implementing parametrizable impact parameter distributions.
void SetDeltaB_WithConstantCrossSection(Double_t deltab, Double_t sigmaR=0)
Handles lists of named parameters with different types, a list of KVNamedParameter objects.
Range of values specified by minimum, maximum.
Definition: KVValueRange.h:19
ValueType Max() const
Definition: KVValueRange.h:48
ValueType ValueIofN(Int_t i, Int_t n) const
Definition: KVValueRange.h:60
ValueType Min() const
Definition: KVValueRange.h:44
virtual void SetLineColor(Color_t lcolor)
virtual Double_t GetBinCenter(Int_t bin) const
virtual Double_t GetBinLowEdge(Int_t bin) const
virtual Double_t GetBinWidth(Int_t bin) const
virtual Double_t GetBinUpEdge(Int_t bin) const
virtual Double_t Mean(Double_t a, Double_t b, const Double_t *params=nullptr, Double_t epsilon=0.000001)
virtual TH1 * GetHistogram() const
virtual void SetRange(Double_t xmin, Double_t xmax)
virtual Double_t Integral(Double_t a, Double_t b, Double_t epsrel=1.e-12)
void SetTitle(const char *title="") override
void Draw(Option_t *option="") override
virtual Double_t Variance(Double_t a, Double_t b, const Double_t *params=nullptr, Double_t epsilon=0.000001)
virtual TF1 * DrawCopy(Option_t *option="") const
virtual Double_t GetRandom(Double_t xmin, Double_t xmax, TRandom *rng=nullptr, Option_t *opt=nullptr)
virtual void SetParameters(const Double_t *params)
virtual void SetParNames(const char *name0="p0", const char *name1="p1", const char *name2="p2", const char *name3="p3", const char *name4="p4", const char *name5="p5", const char *name6="p6", const char *name7="p7", const char *name8="p8", const char *name9="p9", const char *name10="p10")
virtual Double_t Eval(Double_t x, Double_t y=0, Double_t z=0, Double_t t=0) const
virtual void SetParameter(const TString &name, Double_t value)
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
virtual Double_t GetBinCenter(Int_t bin) const
virtual Int_t GetNbinsY() const
TAxis * GetXaxis()
TObject * FindObject(const char *name) const override
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
virtual Int_t GetNbinsX() const
TAxis * GetYaxis()
virtual Double_t GetBinLowEdge(Int_t bin) const
virtual Double_t Integral(Int_t binx1, Int_t binx2, Option_t *option="") const
virtual Double_t GetBinContent(Int_t bin) const
virtual Double_t GetBinWidth(Int_t bin) const
virtual void Scale(Double_t c1=1, Option_t *option="")
virtual Int_t FindBin(Double_t x, Double_t y=0, Double_t z=0)
virtual Double_t GetBinContent(Int_t bin) const
const char * GetTitle() const override
virtual Double_t Uniform(Double_t x1, Double_t x2)
Bool_t IsNull() const
double binomial_pdf(unsigned int k, double p, unsigned int n)
double negative_binomial_pdf(unsigned int k, double p, double n)
double gamma_pdf(double x, double alpha, double theta, double x0=0)
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
const Double_t sigma
Double_t x[n]
TGraphErrors * gr
TH1 * h
fit(model, train_loader, val_loader, num_epochs, batch_size, optimizer, criterion, save_best, scheduler)
Int_t Nint(T x)
Double_t Sqrt(Double_t x)