GooFit  v2.1.3
2d_plot.cpp
Go to the documentation of this file.
1 #include <goofit/Application.h>
2 #include <goofit/FitManager.h>
6 #include <goofit/Variable.h>
8 
9 #include <RVersion.h>
10 #include <TCanvas.h>
11 #include <TH1F.h>
12 #include <TH2F.h>
13 
14 #include <iostream>
15 #include <random>
16 #include <sys/time.h>
17 #include <sys/times.h>
18 
19 using namespace std;
20 using namespace GooFit;
21 
22 int main(int argc, char **argv) {
23  GooFit::Application app("2D plot example", argc, argv);
24 
25  GOOFIT_PARSE(app);
26 
27  // In real code, use a random device here
28  std::mt19937 gen(137);
29  std::normal_distribution<> dx(0.2, 1.1);
30  std::normal_distribution<> dy(0.5, 0.3);
31 
32  GooFit::setROOTStyle();
33 
34  Observable xvar{"xvar", -5, 5};
35  Observable yvar{"yvar", -5, 5};
36  UnbinnedDataSet data({xvar, yvar});
37 
38  TH2F dataHist("dataHist",
39  "",
40  xvar.getNumBins(),
41  xvar.getLowerLimit(),
42  xvar.getUpperLimit(),
43  yvar.getNumBins(),
44  yvar.getLowerLimit(),
45  yvar.getUpperLimit());
46  TH1F xvarHist("xvarHist", "", xvar.getNumBins(), xvar.getLowerLimit(), xvar.getUpperLimit());
47  TH1F yvarHist("yvarHist", "", yvar.getNumBins(), yvar.getLowerLimit(), yvar.getUpperLimit());
48 
49  dataHist.SetStats(false);
50  xvarHist.SetStats(false);
51  yvarHist.SetStats(false);
52 
53  size_t totalData = 0;
54 
55  while(totalData < 100000) {
56  xvar.setValue(dx(gen));
57  yvar.setValue(dy(gen));
58 
59  if(fabs(xvar.getValue()) < 5 && fabs(yvar.getValue()) < 5) {
60  data.addEvent();
61  dataHist.Fill(xvar.getValue(), yvar.getValue());
62  xvarHist.Fill(xvar.getValue());
63  yvarHist.Fill(yvar.getValue());
64  totalData++;
65  }
66  }
67 
68  Variable xmean{"xmean", 0, 1, -10, 10};
69  Variable xsigm{"xsigm", 1, 0.5, 1.5};
70  GaussianPdf xgauss{"xgauss", xvar, xmean, xsigm};
71 
72  Variable ymean{"ymean", 0, 1, -10, 10};
73  Variable ysigm{"ysigm", 0.4, 0.1, 0.6};
74  GaussianPdf ygauss{"ygauss", yvar, ymean, ysigm};
75 
76  ProdPdf total("total", {&xgauss, &ygauss});
77  total.setData(&data);
78 
79  FitManager fitter(&total);
80  fitter.fit();
81 
82  TH2F pdfHist("pdfHist",
83  "",
84  xvar.getNumBins(),
85  xvar.getLowerLimit(),
86  xvar.getUpperLimit(),
87  yvar.getNumBins(),
88  yvar.getLowerLimit(),
89  yvar.getUpperLimit());
90  TH1F xpdfHist("xpdfHist", "", xvar.getNumBins(), xvar.getLowerLimit(), xvar.getUpperLimit());
91  TH1F ypdfHist("ypdfHist", "", yvar.getNumBins(), yvar.getLowerLimit(), yvar.getUpperLimit());
92 
93  pdfHist.SetStats(false);
94  xpdfHist.SetStats(false);
95  ypdfHist.SetStats(false);
96 
97  UnbinnedDataSet grid = total.makeGrid();
98  total.setData(&grid);
99  std::vector<std::vector<double>> pdfVals = total.getCompProbsAtDataPoints();
100 
101  TCanvas foo;
102  dataHist.Draw("colz");
103  foo.SaveAs("data.png");
104 
105  double totalPdf = 0;
106 
107  for(int i = 0; i < grid.getNumEvents(); ++i) {
108  grid.loadEvent(i);
109  pdfHist.Fill(xvar.getValue(), yvar.getValue(), pdfVals[0][i]);
110  xpdfHist.Fill(xvar.getValue(), pdfVals[0][i]);
111  ypdfHist.Fill(yvar.getValue(), pdfVals[0][i]);
112  totalPdf += pdfVals[0][i];
113  }
114 
115  for(int i = 0; i < xvar.getNumBins(); ++i) {
116  double val = xpdfHist.GetBinContent(i + 1);
117  val /= totalPdf;
118  val *= totalData;
119  xpdfHist.SetBinContent(i + 1, val);
120  }
121 
122  for(int i = 0; i < yvar.getNumBins(); ++i) {
123  double val = ypdfHist.GetBinContent(i + 1);
124  val /= totalPdf;
125  val *= totalData;
126  ypdfHist.SetBinContent(i + 1, val);
127 
128  for(int j = 0; j < xvar.getNumBins(); ++j) {
129  val = pdfHist.GetBinContent(j + 1, i + 1);
130  val /= totalPdf;
131  val *= totalData;
132  pdfHist.SetBinContent(j + 1, i + 1, val);
133  }
134  }
135 
136  pdfHist.Draw("colz");
137  foo.SaveAs("pdf.png");
138 
139  for(int i = 0; i < yvar.getNumBins(); ++i) {
140  for(int j = 0; j < xvar.getNumBins(); ++j) {
141  double pval = pdfHist.GetBinContent(j + 1, i + 1);
142  double dval = dataHist.GetBinContent(j + 1, i + 1);
143  pval -= dval;
144  pval /= std::max(1.0, sqrt(dval));
145  pdfHist.SetBinContent(j + 1, i + 1, pval);
146  }
147  }
148 
149  pdfHist.GetZaxis()->SetRangeUser(-5, 5);
150  pdfHist.Draw("colz");
151  foo.SaveAs("pull.png");
152 
153  xvarHist.SetMarkerStyle(8);
154  xvarHist.SetMarkerSize(0.5);
155  xvarHist.Draw("p");
156  xpdfHist.SetLineColor(kBlue);
157  xpdfHist.SetLineWidth(3);
158  xpdfHist.Draw("lsame");
159  foo.SaveAs("xhist.png");
160 
161  yvarHist.SetMarkerStyle(8);
162  yvarHist.SetMarkerSize(0.5);
163  yvarHist.Draw("p");
164  ypdfHist.SetLineColor(kBlue);
165  ypdfHist.SetLineWidth(3);
166  ypdfHist.Draw("lsame");
167  foo.SaveAs("yhist.png");
168 
169  return fitter;
170 }
size_t getNumEvents() const
Definition: DataSet.h:47
int main(int argc, char **argv)
Definition: 2d_plot.cpp:22
Special class for observables. Used in DataSets.
Definition: Variable.h:109
UnbinnedDataSet * data
__host__ void setData(DataSet *data)
TCanvas foo
Definition: chisquare.cpp:22
ROOT::Minuit2::FunctionMinimum fit()
This runs the fit.
#define GOOFIT_PARSE(app,...)
Definition: Application.h:11
void loadEvent(size_t idx)
Set all the variables to the current event values.