GooFit  v2.1.3
zachFit.cpp
Go to the documentation of this file.
1 // ROOT stuff
2 #include <TCanvas.h>
3 #include <TFile.h>
4 #include <TH1F.h>
5 #include <TRandom.h>
6 #include <TStyle.h>
7 
8 // System stuff
9 #include <fstream>
10 
11 // GooFit stuff
12 #include <goofit/Application.h>
13 #include <goofit/BinnedDataSet.h>
14 #include <goofit/FitControl.h>
15 #include <goofit/FitManager.h>
16 #include <goofit/Log.h>
17 #include <goofit/UnbinnedDataSet.h>
18 #include <goofit/Variable.h>
19 
26 
27 #include <fmt/format.h>
28 
29 using namespace fmt::literals;
30 using namespace GooFit;
31 
32 TH1D *getData(DataSet *data, Observable var, std::string filename, size_t reduce = 1, std::string keyw = "data_hist") {
33  TH1D *hist = new TH1D(keyw.c_str(), "", 300, 0.1365, 0.1665);
34 
35  std::ifstream reader{filename};
36 
37  size_t val_read_in = 0;
38  while(reader >> var) {
39  if(var) {
40  if(val_read_in % reduce == 0) {
41  data->addEvent();
42  hist->Fill(var.getValue());
43  }
44  val_read_in++;
45  }
46  }
47 
48  hist->SetStats(false);
49  hist->SetMarkerStyle(8);
50  hist->SetMarkerSize(0.6);
51 
52  std::size_t pos = keyw.find("_");
53  std::string keywu = keyw.substr(0, pos);
54 
55  std::string outp1 = keywu + " events: {}";
56  std::string outp2 = "No" + keywu + " events read in!";
57 
58  GOOFIT_INFO(outp1, data->getNumEvents());
59  if(data->getNumEvents() == 0)
60  throw GooFit::GeneralError(outp2);
61  return hist;
62 }
63 
64 int main(int argc, char **argv) {
66  R"raw(This example performs a staged fit measuring the mass difference between the `D*(2010)+`
67 and `D0` using D*+ -> D0 pi+ events recorded by the BaBar detector (approximately 477
68 inverse femtobarn).
69 )raw",
70  argc,
71  argv};
72 
73  app.set_footer(R"raw(Dataset descriptions:
74 0-simple Early testing sample for GooFit before nominal dataset was released.
75  MC resolution sample and data for channel D*+ -> D0 pi+; D0 -> K- pi+
76  Samples are composed of events that pass the majority of selection criteria, but
77  fail at least one of the stricter tracking cuts. The resulting resolution is worse
78  than in the events of the nominal samples used in the official analysis/publication
79  marked below as data set options "1" and "2".
80 1-kpi Nominal MC resolution sample and data for channel D*+ -> D0 pi+; D0 -> K- pi+
81 2-k3pi Nominal MC resolution sample and data for channel D*+ -> D0 pi+; D0 -> K- pi+ pi- pi+)raw");
82 
83  int mode = 0, data = 0;
84  bool plot;
85  size_t reduce = 1;
86  app.add_set("-m,--mode,mode", mode, {0, 1, 2}, "Program mode: 0-unbinned, 1-binned, 2-binned chisq", true);
87  app.add_set("-d,--data,data", data, {0, 1, 2}, "Dataset: 0-simple, 1-kpi, 2-k3pi", true);
88  app.add_option("--reduce", reduce, "Load every X line of data", true);
89  app.add_flag("-p,--plot", plot, "Make and save plots of results");
90 
91  GOOFIT_PARSE(app);
92 
93  // Style
94  gStyle->SetCanvasBorderMode(0);
95  gStyle->SetCanvasColor(10);
96  gStyle->SetFrameFillColor(10);
97  gStyle->SetFrameBorderMode(0);
98  gStyle->SetPadColor(0);
99  gStyle->SetTitleColor(1);
100  gStyle->SetStatColor(0);
101  gStyle->SetFillColor(0);
102  gStyle->SetFuncWidth(1);
103  gStyle->SetLineWidth(1);
104  gStyle->SetLineColor(1);
105  gStyle->SetPalette(1, 0);
106 
107  TCanvas foo;
108  foo.SetLogy(true);
109 
110  // Get the name of the files to use
111  std::string mcfile, datafile;
112  if(data == 0) {
113  mcfile = app.get_filename("dataFiles/dstwidth_kpi_resMC.dat", "examples/zachFit");
114  datafile = app.get_filename("dataFiles/dstwidth_kpi_data.dat", "examples/zachFit");
115  } else if(data == 1) {
116  mcfile = app.get_filename("dataFiles/DstarWidth_D0ToKpi_deltaM_MC.dat", "examples/zachFit");
117  datafile = app.get_filename("dataFiles/DstarWidth_D0ToKpi_deltaM_Data.dat", "examples/zachFit");
118  } else {
119  mcfile = app.get_filename("dataFiles/DstarWidth_D0ToK3pi_deltaM_MC.dat", "examples/zachFit");
120  datafile = app.get_filename("dataFiles/DstarWidth_D0ToK3pi_deltaM_Data.dat", "examples/zachFit");
121  }
122 
123  Observable dm{"dm", 0.1395, 0.1665};
124  dm.setNumBins(2700);
125 
126  // This would be clearer with std::optional from C++17
127  std::unique_ptr<DataSet> mc_dataset, data_dataset;
128 
129  if(mode == 0) {
130  mc_dataset.reset(new UnbinnedDataSet{dm});
131  data_dataset.reset(new UnbinnedDataSet{dm});
132  } else {
133  mc_dataset.reset(new BinnedDataSet{dm});
134  data_dataset.reset(new BinnedDataSet{dm});
135  }
136 
137  TH1D *mc_hist = getData(mc_dataset.get(), dm, mcfile, reduce, "mc_hist");
138 
139  Variable mean1("kpi_mc_mean1", 0.145402, 0.00001, 0.143, 0.148);
140  Variable mean2("kpi_mc_mean2", 0.145465, 0.00001, 0.145, 0.1465);
141  Variable mean3("kpi_mc_mean3", 0.145404, 0.00001, 0.144, 0.147);
142 
143  Variable sigma1("kpi_mc_sigma1", 0.00010, 0.00001, 0.000001, 0.002);
144  Variable sigma2("kpi_mc_sigma2", 0.00075, 0.00001, 0.000001, 0.005);
145  Variable sigma3("kpi_mc_sigma3", 0.00020, 0.00001, 0.000005, 0.001);
146 
147  Variable pimass("kpi_mc_pimass", 0.13957);
148  Variable aslope("kpi_mc_aslope", -20.0, 1, -100.0, 10.0);
149  Variable apower("kpi_mc_apower", 1.3, 0.1, 0.1, 10.0);
150  Variable gfrac1("kpi_mc_gfrac1", 0.65, 0.01, 0.0, 0.9);
151  Variable gfrac2("kpi_mc_gfrac2", 0.02, 0.001, 0.0, 0.12);
152  Variable afrac("kpi_mc_afrac", 0.005, 0.003, 0.0, 0.1);
153 
154  GaussianPdf gauss1("gauss1", dm, mean1, sigma1);
155  GaussianPdf gauss2("gauss2", dm, mean2, sigma2);
156  GaussianPdf gauss3("gauss3", dm, mean3, sigma3);
157  ArgusPdf argus("argus", dm, pimass, aslope, false, apower);
158 
159  AddPdf resolution{"resolution", {gfrac1, gfrac2, afrac}, {&gauss1, &gauss2, &argus, &gauss3}};
160 
161  resolution.setData(mc_dataset.get());
162 
163  FitManager mcpdf{&resolution};
164 
165  GOOFIT_INFO("Done with collecting MC, starting minimisation");
166  mcpdf.fit();
167 
168  if(plot) {
169  GOOFIT_INFO("Plotting MC");
170  mc_hist->SetLineColor(kBlack);
171  mc_hist->Draw("e");
172 
173  double step = mc_hist->GetXaxis()->GetBinWidth(2);
174  auto tot_hist = resolution.plotToROOT(dm, mc_dataset->getNumEvents() * step);
175  tot_hist->SetLineColor(kGreen);
176 
177  tot_hist->Draw("SAME");
178 
179  foo.SaveAs("MC_plot.png");
180  }
181 
182  // Locking the MC variables
183  mean1.setFixed(true);
184  mean2.setFixed(true);
185  mean3.setFixed(true);
186  sigma1.setFixed(true);
187  sigma2.setFixed(true);
188  sigma3.setFixed(true);
189  pimass.setFixed(true);
190  aslope.setFixed(true);
191  gfrac1.setFixed(true);
192  gfrac2.setFixed(true);
193  afrac.setFixed(true);
194  apower.setFixed(true);
195 
196  Variable dummyzero("kpi_rd_dummyzero", 0);
197  Variable delta("kpi_rd_delta", 0.000002, -0.00005, 0.00005);
198  Variable epsilon("kpi_rd_epsilon", 0.05, -0.1, 0.2);
199 
200  ScaledGaussianPdf resolution1("resolution1", dm, dummyzero, sigma1, delta, epsilon);
201  ScaledGaussianPdf resolution2("resolution2", dm, dummyzero, sigma2, delta, epsilon);
202  ScaledGaussianPdf resolution3("resolution3", dm, dummyzero, sigma3, delta, epsilon);
203 
204  Variable width_bw("kpi_rd_width_bw", 0.0001, 0.00001, 0.0005);
205  KinLimitBWPdf rbw1("rbw1", dm, mean1, width_bw);
206  KinLimitBWPdf rbw2("rbw2", dm, mean2, width_bw);
207  KinLimitBWPdf rbw3("rbw3", dm, mean3, width_bw);
208 
209  ConvolutionPdf signal1{"signal1", dm, &rbw1, &resolution1};
210  ConvolutionPdf signal2{"signal2", dm, &rbw2, &resolution2};
211  ConvolutionPdf signal3{"signal3", dm, &rbw3, &resolution3};
212 
213  signal1.setIntegrationConstants(0.1395, 0.1665, 0.0000027);
214  signal2.setIntegrationConstants(0.1395, 0.1665, 0.0000027);
215  signal3.setIntegrationConstants(0.1395, 0.1665, 0.0000027);
216 
217  AddPdf signal{"signal", {gfrac1, gfrac2, afrac}, {&signal1, &signal2, &argus, &signal3}};
218 
219  Variable slope("kpi_rd_slope", -1.0, 0.1, -35.0, 25.0);
220  ArgusPdf bkg("bkg", dm, pimass, slope, false);
221 
222  Variable bkg_frac("kpi_rd_bkg_frac", 0.03, 0.0, 0.3);
223 
224  TH1D *data_hist = getData(data_dataset.get(), dm, datafile, reduce, "data_hist");
225 
226  AddPdf total("total", {bkg_frac}, {&bkg, &signal});
227 
228  total.setData(data_dataset.get());
229 
230  std::shared_ptr<BinnedChisqFit> chi_control;
231  if(2 == mode) {
232  chi_control.reset(new BinnedChisqFit);
233  total.setFitControl(chi_control);
234  }
235 
236  FitManager datapdf{&total};
237 
238  GOOFIT_INFO("Starting fit");
239 
240  datapdf.fit();
241 
242  if(plot) {
243  GOOFIT_INFO("Plotting results");
244 
245  data_hist->SetLineColor(kBlack);
246  data_hist->Draw("e");
247 
248  double scale = data_hist->GetXaxis()->GetBinWidth(2) * data_dataset->getNumEvents();
249 
250  auto sig_hist = signal.plotToROOT(dm, (1 - bkg_frac.getValue()) * scale);
251  sig_hist->SetLineColor(kBlue);
252  auto back_hist = bkg.plotToROOT(dm, bkg_frac.getValue() * scale);
253  back_hist->SetLineColor(kRed);
254  auto tot_hist = total.plotToROOT(dm, scale);
255  tot_hist->SetLineColor(kGreen);
256 
257  tot_hist->Draw("SAME");
258  sig_hist->Draw("SAME");
259  back_hist->Draw("SAME");
260 
261  foo.SaveAs("ResultFit.png");
262  }
263 
264  return datapdf;
265 }
Thrown when a general error is encountered.
Definition: Error.h:10
size_t getNumEvents() const
Definition: DataSet.h:47
void setFixed(bool fix)
Set the fixedness of a variable.
Definition: Variable.h:205
void setNumBins(size_t num)
Set the number of bins.
Definition: Variable.h:128
#define GOOFIT_INFO(...)
Definition: Log.h:10
__host__ __device__ OutputType reduce(goofit_policy &exec, InputIterator first, InputIterator last, OutputType init, BinaryFunction binary_op)
Special class for observables. Used in DataSets.
Definition: Variable.h:109
UnbinnedDataSet * data
__host__ void setData(DataSet *data)
TH1D * getData(DataSet *data, Observable var, std::string filename, size_t reduce=1, std::string keyw="data_hist")
Definition: zachFit.cpp:32
TCanvas foo
Definition: chisquare.cpp:22
fptype getValue() const
Get the value.
Definition: Variable.h:68
#define GOOFIT_PARSE(app,...)
Definition: Application.h:11
virtual void addEvent()=0
int main(int argc, char **argv)
Definition: zachFit.cpp:64