GooFit  v2.1.3
simpleFit.cpp
Go to the documentation of this file.
1 #include <goofit/Application.h>
2 #include <goofit/FitManager.h>
7 #include <goofit/Variable.h>
9 
10 #include <TCanvas.h>
11 #include <TH1F.h>
12 #include <TRandom.h>
13 #include <TStyle.h>
14 
15 #include <iostream>
16 
17 using namespace std;
18 using namespace GooFit;
19 
20 // CPU-side Novosibirsk evaluation for use in generating toy MC.
21 double novosib(double x, double peak, double width, double tail) {
22  double qa = 0, qb = 0, qc = 0, qx = 0, qy = 0;
23 
24  if(fabs(tail) < 1.e-7)
25  qc = 0.5 * pow(((x - peak) / width), 2);
26  else {
27  qa = tail * sqrt(log(4.));
28  qb = sinh(qa) / qa;
29  qx = (x - peak) / width * qb;
30  qy = 1. + tail * qx;
31 
32  //---- Cutting curve from right side
33  if(qy > 1.E-7)
34  qc = 0.5 * (pow((log(qy) / tail), 2) + tail * tail);
35  else
36  qc = 15.0;
37  }
38 
39  //---- Normalize the result
40  return exp(-qc);
41 }
42 
43 TCanvas *foo = 0;
44 
45 void fitAndPlot(GooPdf *total, UnbinnedDataSet *data, TH1F &dataHist, Observable xvar, const char *fname) {
46  total->setData(data);
47  FitManager fitter(total);
48  fitter.fit();
49 
50  if(!fitter)
51  std::exit(fitter);
52 
53  TH1F pdfHist("pdfHist", "", xvar.getNumBins(), xvar.getLowerLimit(), xvar.getUpperLimit());
54  pdfHist.SetStats(false);
55 
56  UnbinnedDataSet grid = total->makeGrid();
57  total->setData(&grid);
58  std::vector<std::vector<double>> pdfVals = total->getCompProbsAtDataPoints();
59 
60  double totalPdf = 0;
61 
62  for(int i = 0; i < grid.getNumEvents(); ++i) {
63  grid.loadEvent(i);
64  pdfHist.Fill(xvar.getValue(), pdfVals[0][i]);
65  totalPdf += pdfVals[0][i];
66  }
67 
68  for(int i = 0; i < xvar.getNumBins(); ++i) {
69  double val = pdfHist.GetBinContent(i + 1);
70  val /= totalPdf;
71  val *= data->getNumEvents();
72  pdfHist.SetBinContent(i + 1, val);
73  }
74 
75  // foo->SetLogy(true);
76  dataHist.SetMarkerStyle(8);
77  dataHist.SetMarkerSize(0.5);
78  dataHist.Draw("p");
79  pdfHist.SetLineColor(kBlue);
80  pdfHist.SetLineWidth(3);
81  pdfHist.Draw("lsame");
82  foo->SaveAs(fname);
83 }
84 
85 int main(int argc, char **argv) {
86  GooFit::Application app("Simple fit example", argc, argv);
87 
88  size_t numevents = 100000;
89  app.add_option("-n,--num", numevents, "Number of events", true);
90 
91  GOOFIT_PARSE(app);
92 
93  GooFit::setROOTStyle();
94 
95  // Independent variable.
96  Observable xvar("xvar", -100, 100);
97  xvar.setNumBins(1000); // For such a large range, want more bins for better accuracy in normalisation.
98 
99  // Data sets for the three fits.
100  UnbinnedDataSet landdata(xvar);
101  UnbinnedDataSet bifgdata(xvar);
102  UnbinnedDataSet novodata(xvar);
103 
104  // Histograms for showing the fit.
105  TH1F landHist("landHist", "", xvar.getNumBins(), xvar.getLowerLimit(), xvar.getUpperLimit());
106  TH1F bifgHist("bifgHist", "", xvar.getNumBins(), xvar.getLowerLimit(), xvar.getUpperLimit());
107  TH1F novoHist("novoHist", "", xvar.getNumBins(), xvar.getLowerLimit(), xvar.getUpperLimit());
108  landHist.SetStats(false);
109  bifgHist.SetStats(false);
110  novoHist.SetStats(false);
111 
112  TRandom donram(42);
113 
114  double maxNovo = 0;
115 
116  for(double x = xvar.getLowerLimit(); x < xvar.getUpperLimit(); x += 0.01) {
117  double curr = novosib(x, 0.3, 0.5, 1.0);
118 
119  if(curr < maxNovo)
120  continue;
121 
122  maxNovo = curr;
123  }
124 
125  double leftSigma = 13;
126  double rightSigma = 29;
127  double leftIntegral = 0.5 / (leftSigma * sqrt(2 * M_PI));
128  double rightIntegral = 0.5 / (rightSigma * sqrt(2 * M_PI));
129  double totalIntegral = leftIntegral + rightIntegral;
130  double bifpoint = -10;
131 
132  // Generating three sets of toy MC.
133  while(landdata.getNumEvents() < numevents) {
134  // Landau
135  try {
136  xvar.setValue(donram.Landau(20, 1));
137  landdata.addEvent();
138  landHist.Fill(xvar.getValue());
139  } catch(const GooFit::OutOfRange &) {
140  }
141  }
142 
143  while(bifgdata.getNumEvents() < numevents) {
144  // Bifurcated Gaussian
145  double val;
146  if(donram.Uniform() < (leftIntegral / totalIntegral)) {
147  do {
148  val = donram.Gaus(bifpoint, rightSigma);
149  } while(val < bifpoint || val > xvar.getUpperLimit());
150  xvar.setValue(val);
151 
152  } else {
153  do {
154  val = donram.Gaus(bifpoint, leftSigma);
155  } while(val > bifpoint || val < xvar.getLowerLimit());
156  xvar.setValue(val);
157  }
158 
159  bifgdata.addEvent();
160  bifgHist.Fill(xvar.getValue());
161  }
162 
163  while(novodata.getNumEvents() < numevents) {
164  // And Novosibirsk.
165  while(true) {
166  xvar.setValue(donram.Uniform(xvar.getLowerLimit(), xvar.getUpperLimit()));
167  double y = donram.Uniform(0, maxNovo);
168 
169  if(y < novosib(xvar.getValue(), 0.3, 0.5, 1.0))
170  break;
171  }
172 
173  novodata.addEvent();
174  novoHist.Fill(xvar.getValue());
175  }
176 
177  foo = new TCanvas();
178 
179  Variable mpv("mpv", 40, 0, 150);
180  Variable sigma("sigma", 5, 0, 30);
181  GooPdf *landau = new LandauPdf("landau", xvar, mpv, sigma);
182  fitAndPlot(landau, &landdata, landHist, xvar, "simple_fit_cpp_landau.png");
183 
184  Variable nmean("nmean", 0.4, -10.0, 10.0);
185  Variable nsigm("nsigm", 0.6, 0.0, 1.0);
186  Variable ntail("ntail", 1.1, 0.1, 0.0, 3.0);
187  GooPdf *novo = new NovosibirskPdf("novo", xvar, nmean, nsigm, ntail);
188  fitAndPlot(novo, &novodata, novoHist, xvar, "simple_fit_cpp_novo.png");
189 
190  Variable gmean("gmean", 3.0, 1, -15, 15);
191  Variable lsigm("lsigm", 10, 1, 10, 20);
192  Variable rsigm("rsigm", 20, 1, 10, 40);
193  GooPdf *bifur = new BifurGaussPdf("bifur", xvar, gmean, lsigm, rsigm);
194  fitAndPlot(bifur, &bifgdata, bifgHist, xvar, "simple_fit_cpp_bifur.png");
195 
196  return 0;
197 }
size_t getNumEvents() const
Definition: DataSet.h:47
void setNumBins(size_t num)
Set the number of bins.
Definition: Variable.h:128
TCanvas * foo
Definition: simpleFit.cpp:43
__host__ std::vector< std::vector< fptype > > getCompProbsAtDataPoints()
Produce a list of probabilies at points.
double novosib(double x, double peak, double width, double tail)
Definition: simpleFit.cpp:21
int main(int argc, char **argv)
Definition: simpleFit.cpp:85
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
Special class for observables. Used in DataSets.
Definition: Variable.h:109
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
UnbinnedDataSet * data
void fitAndPlot(GooPdf *total, UnbinnedDataSet *data, TH1F &dataHist, Observable xvar, const char *fname)
Definition: simpleFit.cpp:45
__host__ void setData(DataSet *data)
Observable * sigma
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
fptype getValue() const
Get the value.
Definition: Variable.h:68
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
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.
__host__ UnbinnedDataSet makeGrid()
Set an equidistant grid based on the stored variable binning.