GooFit  v2.1.3
chisquare.cpp
Go to the documentation of this file.
1 #include <goofit/Application.h>
2 #include <goofit/BinnedDataSet.h>
3 #include <goofit/FitControl.h>
4 #include <goofit/FitManager.h>
6 #include <goofit/Variable.h>
7 
8 #include <CLI/Timer.hpp>
9 
10 #include <TCanvas.h>
11 #include <TH1F.h>
12 #include <TLatex.h>
13 #include <TMinuit.h>
14 #include <TRandom.h>
16 
17 #include <iostream>
18 #include <string>
19 #include <tuple>
20 #include <vector>
21 
22 TCanvas foo;
23 
24 using namespace std;
25 using namespace GooFit;
26 
27 vector<double> ratios;
28 vector<double> errors;
29 
30 // Global needed to put it in the fit function for ROOT
31 Observable decayTime{"decayTime", 0, 10};
32 
33 double integralExpCon(double lo, double hi) { return (exp(-lo) - exp(-hi)); }
34 
35 double integralExpLin(double lo, double hi) { return ((lo + 1) * exp(-lo) - (hi + 1) * exp(-hi)); }
36 
37 double integralExpSqu(double lo, double hi) {
38  return ((lo * lo + 2 * lo + 2) * exp(-lo) - (hi * hi + 2 * hi + 2) * exp(-hi));
39 }
40 
42  vector<int> &rsEvtVec,
43  vector<int> &wsEvtVec,
44  double conCoef,
45  double linCoef,
46  double squCoef,
47  int eventsToGenerate) {
48  static TRandom donram(24);
49  double totalRSintegral = integralExpCon(0, 100);
50  double step = (decayTime.getUpperLimit() - decayTime.getLowerLimit()) / decayTime.getNumBins();
51 
52  for(int i = 0; i < decayTime.getNumBins(); ++i) {
53  double binStart = i * step;
54  binStart += decayTime.getLowerLimit();
55  double binFinal = binStart + step;
56 
57  double rsIntegral = integralExpCon(binStart, binFinal);
58  double wsIntegral = conCoef * integralExpCon(binStart, binFinal);
59  wsIntegral += linCoef * integralExpLin(binStart, binFinal);
60  wsIntegral += squCoef * integralExpSqu(binStart, binFinal);
61 
62  double expectedRSevts = eventsToGenerate * rsIntegral / totalRSintegral;
63  double expectedWSevts = eventsToGenerate * wsIntegral / totalRSintegral;
64 
65  int rsEvts = donram.Poisson(expectedRSevts);
66  int wsEvts = donram.Poisson(expectedWSevts);
67  rsEvtVec[i] = rsEvts;
68  wsEvtVec[i] = wsEvts;
69 
70  if(0 == (i % 10))
71  std::cout << "Events in bin " << i << " : " << rsEvts << " (" << expectedRSevts << ") " << wsEvts << " ("
72  << expectedWSevts << ")\n";
73  }
74 }
75 
76 std::tuple<int, std::string> fitRatio(Observable decayTime,
77  vector<Variable> weights,
78  vector<int> &rsEvts,
79  vector<int> &wsEvts,
80  std::string plotName = "") {
81  TH1D ratioHist("ratioHist", "", decayTime.getNumBins(), decayTime.getLowerLimit(), decayTime.getUpperLimit());
82 
83  BinnedDataSet ratioData(decayTime);
84 
85  for(unsigned int i = 0; i < wsEvts.size(); ++i) {
86  double ratio = wsEvts[i];
87 
88  if(0 == rsEvts[i])
89  rsEvts[i] = 1; // Cheating to avoid div by zero.
90 
91  ratio /= rsEvts[i];
92 
93  if(0 == wsEvts[i])
94  wsEvts[i] = 1; // Avoid zero errors
95 
96  double error = wsEvts[i] / pow(rsEvts[i], 2);
97  error += pow(wsEvts[i], 2) / pow(rsEvts[i], 3);
98  error = sqrt(error);
99 
100  ratioData.setBinContent(i, ratio);
101  ratioData.setBinError(i, error);
102  ratioHist.SetBinContent(i + 1, ratio);
103  ratioHist.SetBinError(i + 1, error);
104  }
105 
106  PolynomialPdf poly("poly", decayTime, weights);
107  poly.setFitControl(std::make_shared<BinnedErrorFit>());
108  poly.setData(&ratioData);
109  FitManager fitter{&poly};
110 
111  CLI::Timer timer_cpu{"GPU"};
112  fitter.fit();
113  std::string timer_str = timer_cpu.to_string();
114 
115  if(!plotName.empty()) {
116  vector<fptype> values = poly.evaluateAtPoints(decayTime);
117  TH1D pdfHist("pdfHist", "", decayTime.getNumBins(), decayTime.getLowerLimit(), decayTime.getUpperLimit());
118 
119  for(int i = 0; i < values.size(); ++i) {
120  pdfHist.SetBinContent(i + 1, values[i]);
121  }
122 
123  ratioHist.SetMarkerStyle(8);
124  ratioHist.SetMarkerSize(0.5);
125  ratioHist.SetStats(false);
126  ratioHist.Draw("p");
127 
128  TString str1 = fmt::format(
129  "Constant [10^{{-2}}] : {:.3} #pm {:.3}", weights[0].getValue() * 1e2, weights[0].getError() * 1e2);
130  TLatex res1(0.14, 0.83, str1);
131  res1.SetNDC(true);
132 
133  TString str2 = fmt::format(
134  "Linear [10^{{-4}}] : {:.3} #pm {:.3}", weights[1].getValue() * 1e4, weights[1].getError() * 1e4);
135  TLatex res2(0.14, 0.73, str2);
136  res2.SetNDC(true);
137 
138  TString str3 = fmt::format(
139  "Quadratic [10^{{-6}}]: {:.3} #pm {:.3}", weights[2].getValue() * 1e6, weights[2].getError() * 1e6);
140  TLatex res3(0.14, 0.63, str3);
141  res3.SetNDC(true);
142 
143  res1.Draw();
144  res2.Draw();
145  res3.Draw();
146 
147  pdfHist.SetLineColor(kBlue);
148  pdfHist.SetLineWidth(3);
149  pdfHist.SetStats(false);
150  pdfHist.Draw("lsame");
151  foo.SaveAs(plotName.c_str());
152  }
153 
154  std::cout << "Polynomial function: " << poly.getCoefficient(2) << " * t^2 + " << poly.getCoefficient(1) << " * t + "
155  << poly.getCoefficient(0) << std::endl;
156 
157  return make_tuple(int(fitter), timer_str);
158 }
159 
160 void cpvFitFcn(int &npar, double *gin, double &fun, double *fp, int iflag) {
161  double conCoef = fp[0];
162  double linCoef = fp[1];
163  double squCoef = fp[2];
164 
165  double chisq = 0;
167 
168  for(unsigned int i = 0; i < ratios.size(); ++i) {
169  double currDTime = decayTime.getLowerLimit() + (i + 0.5) * step;
170  double pdfval = conCoef + linCoef * currDTime + squCoef * currDTime * currDTime;
171  chisq += pow((pdfval - ratios[i]) / errors[i], 2);
172  }
173 
174  fun = chisq;
175 }
176 
177 void fitRatioCPU(Observable decayTime, vector<int> &rsEvts, vector<int> &wsEvts) {
178  TH1D *ratioHist
179  = new TH1D("ratioHist", "", decayTime.getNumBins(), decayTime.getLowerLimit(), decayTime.getUpperLimit());
180 
181  ratios.resize(wsEvts.size());
182  errors.resize(wsEvts.size());
183 
184  for(unsigned int i = 0; i < wsEvts.size(); ++i) {
185  if(0 == rsEvts[i])
186  rsEvts[i] = 1; // Cheating to avoid div by zero.
187 
188  fptype ratio = wsEvts[i] / rsEvts[i];
189 
190  if(0 == wsEvts[i])
191  wsEvts[i] = 1; // Avoid zero errors
192 
193  double error = wsEvts[i] / pow(rsEvts[i], 2);
194  error += pow(wsEvts[i], 2) / pow(rsEvts[i], 3);
195  error = sqrt(error);
196 
197  ratios[i] = ratio;
198  errors[i] = error;
199  ratioHist->SetBinContent(i + 1, ratio);
200  ratioHist->SetBinError(i + 1, error);
201  }
202 
203  TMinuit *minuit = new TMinuit(3);
204  minuit->DefineParameter(0, "constaCoef", 0.03, 0.01, -1, 1);
205  minuit->DefineParameter(1, "linearCoef", 0, 0.01, -1, 1);
206  minuit->DefineParameter(2, "secondCoef", 0, 0.01, -1, 1);
207  minuit->SetFCN(cpvFitFcn);
208 
209  minuit->Migrad();
210 }
211 
212 int main(int argc, char **argv) {
213  GooFit::Application app("Chi-square example", argc, argv);
214 
215  int numbins = 100;
216  app.add_option("-n,--numbins", numbins, "Number of bins", true);
217 
218  int eventsToGenerate = 10000000;
219  app.add_option("-e,--events", eventsToGenerate, "Events to generate", true);
220 
221  GOOFIT_PARSE(app);
222 
223  // Time is in units of lifetime
224 
225  decayTime.setValue(100.);
226  decayTime.setNumBins(numbins);
227 
228  double rSubD = 0.03;
229  double rBarD = 0.03;
230  double delta = 0;
231  double wpPhi = 0;
232  double x_mix = 0.0016;
233  double y_mix = 0.0055;
234  double magPQ = 1.0;
235  double magQP = 1.0 / magPQ;
236 
237  vector<int> dZeroEvtsWS(decayTime.getNumBins());
238  vector<int> dZeroEvtsRS(decayTime.getNumBins());
239  vector<int> d0barEvtsWS(decayTime.getNumBins());
240  vector<int> d0barEvtsRS(decayTime.getNumBins());
241 
242  double dZeroLinearCoef = magPQ * sqrt(rSubD) * (y_mix * cos(delta + wpPhi) - x_mix * sin(delta + wpPhi));
243  double d0barLinearCoef = magQP * sqrt(rBarD) * (y_mix * cos(delta - wpPhi) - x_mix * sin(delta - wpPhi));
244 
245  double dZeroSecondCoef = 0.25 * magPQ * magPQ * (x_mix * x_mix + y_mix * y_mix);
246  double d0barSecondCoef = 0.25 * magQP * magQP * (x_mix * x_mix + y_mix * y_mix);
247 
248  generateEvents(decayTime, dZeroEvtsRS, dZeroEvtsWS, rSubD, dZeroLinearCoef, dZeroSecondCoef, eventsToGenerate);
249  generateEvents(decayTime, d0barEvtsRS, d0barEvtsWS, rBarD, d0barLinearCoef, d0barSecondCoef, eventsToGenerate);
250 
251  Variable constaCoef("constaCoef", 0.03, 0.01, -1, 1);
252  Variable linearCoef("linearCoef", 0, 0.01, -1, 1);
253  Variable secondCoef("secondCoef", 0, 0.01, -1, 1);
254 
255  vector<Variable> weights = {constaCoef, linearCoef, secondCoef};
256 
257  int retval1, retval2;
258  std::string fit1, fit2;
259  std::tie(retval1, fit1)
260  = fitRatio(decayTime, weights, dZeroEvtsRS, dZeroEvtsWS, "chisquare_dzeroEvtRatio_goo_cpp.png");
261  std::tie(retval2, fit2)
262  = fitRatio(decayTime, weights, d0barEvtsRS, d0barEvtsWS, "chisquare_dzbarEvtRatio_goo_cpp.png");
263 
264  CLI::Timer timer_cpu{"Total CPU (2x fits)"};
265  fitRatioCPU(decayTime, dZeroEvtsRS, dZeroEvtsWS);
266  fitRatioCPU(decayTime, d0barEvtsRS, d0barEvtsWS);
267  std::string cpu_string = timer_cpu.to_string();
268 
269  std::cout << fit1 << "\n" << fit2 << "\n" << cpu_string << std::endl;
270 
271  fmt::print("Exit codes (should be 0): {} and {}\n", retval1, retval2);
272 
273  return retval1 + retval2;
274 }
double fptype
void setNumBins(size_t num)
Set the number of bins.
Definition: Variable.h:128
std::tuple< int, std::string > fitRatio(Observable decayTime, vector< Variable > weights, vector< int > &rsEvts, vector< int > &wsEvts, std::string plotName="")
Definition: chisquare.cpp:76
double integralExpSqu(double lo, double hi)
Definition: chisquare.cpp:37
std::vector< Variable > weights
void generateEvents(Observable decayTime, vector< int > &rsEvtVec, vector< int > &wsEvtVec, double conCoef, double linCoef, double squCoef, int eventsToGenerate)
Definition: chisquare.cpp:41
__host__ fptype getCoefficient(int coef) const
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
void fitRatioCPU(Observable decayTime, vector< int > &rsEvts, vector< int > &wsEvts)
Definition: chisquare.cpp:177
Special class for observables. Used in DataSets.
Definition: Variable.h:109
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
int main(int argc, char **argv)
Definition: chisquare.cpp:212
__host__ void setFitControl(std::shared_ptr< FitControl > fc) override
__host__ void setData(DataSet *data)
double integralExpCon(double lo, double hi)
Definition: chisquare.cpp:33
TCanvas foo
Definition: chisquare.cpp:22
void cpvFitFcn(int &npar, double *gin, double &fun, double *fp, int iflag)
Definition: chisquare.cpp:160
__host__ std::vector< fptype > evaluateAtPoints(Observable var)
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
vector< double > ratios
Definition: chisquare.cpp:27
#define GOOFIT_PARSE(app,...)
Definition: Application.h:11
Observable decayTime
Definition: chisquare.cpp:31
vector< double > errors
Definition: chisquare.cpp:28
double integralExpLin(double lo, double hi)
Definition: chisquare.cpp:35