GooFit  v2.1.3
pipipi0DPFit.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 <TH2F.h>
6 #include <TLegend.h>
7 #include <TLine.h>
8 #include <TRandom.h>
9 #include <TRandom3.h>
10 #include <TStyle.h>
11 #include <TText.h>
12 
13 // System stuff
14 #include <cassert>
15 #include <climits>
16 #include <fstream>
17 
18 // GooFit stuff
19 #include <goofit/Application.h>
20 #include <goofit/PDFs/GooPdf.h>
43 #include <goofit/Variable.h>
45 
48 
49 #include <goofit/FunctorWriter.h>
51 #include <goofit/UnbinnedDataSet.h>
52 
53 using namespace std;
54 using namespace GooFit;
55 
56 TCanvas *foo;
57 TCanvas *foodal;
58 
59 UnbinnedDataSet *data = nullptr;
62 TH2F *weightHistogram = nullptr;
63 TH2F *underlyingBins = nullptr;
64 
65 // I hate having to use globals, but this is the best way for now
67 bool minuit1;
68 
69 Observable *m12 = nullptr;
70 Observable *m13 = nullptr;
72 Observable *massd0 = nullptr;
73 Observable *deltam = nullptr;
74 Observable *dtime = nullptr;
75 Observable *sigma = nullptr;
76 Observable *wSig0 = nullptr;
77 Observable *wBkg1 = nullptr;
78 Observable *wBkg2 = nullptr;
79 Observable *wBkg3 = nullptr;
80 Observable *wBkg4 = nullptr;
81 
82 bool fitMasses = false;
83 Variable fixedRhoMass("rho_mass", 0.7758, 0.01, 0.7, 0.8);
84 Variable fixedRhoWidth("rho_width", 0.1503, 0.01, 0.1, 0.2);
85 
86 // Systematic variables
87 double luckyFrac = 0.5;
88 double mesonRad = 1.5;
89 int normBinning = 240;
90 int blindSeed = 4;
91 int mdslices = 1;
92 double md0offset = 0;
95 double deltam_lower_window = -2;
97 double lowerTime = -2;
98 double upperTime = 3;
99 double maxSigma = 0.8;
100 bool polyEff = false;
101 bool saveEffPlot = true;
102 bool useHistogramSigma = false;
105 std::string bkg2Model_str = "sideband";
108 bool makePlots = false;
109 int m23Slices = 6;
110 bool drop_rho_1450 = false;
111 bool drop_rho_1700 = false;
112 bool drop_f0_980 = false;
113 bool drop_f0_1370 = false;
114 bool drop_f0_1500 = false;
115 bool drop_f0_1710 = false;
116 bool drop_f2_1270 = false;
117 bool drop_f0_600 = false;
118 bool gaussBkgTime = false;
119 bool mikhailSetup = false;
120 int bkgHistBins = 80;
121 string paramUp = "";
122 string paramDn = "";
124 
125 const fptype _mD0 = 1.86484;
126 const fptype _mD02 = _mD0 * _mD0;
127 const fptype _mD02inv = 1. / _mD02;
128 const fptype piPlusMass = 0.13957018;
129 const fptype piZeroMass = 0.1349766;
130 enum Resolutions { DplotRes = 1, TimeRes = 2, Efficiency = 4 };
131 
132 size_t maxEvents = 100000;
133 
134 Variable motherM("motherM", 1.86484);
135 Variable chargeM("chargeM", 0.13957018);
136 Variable neutrlM("neutrlM", 0.1349766);
138 
139 // Constants used in more than one PDF component.
140 Variable constantOne("constantOne", 1);
141 Variable constantTwo("constantTwo", 2);
142 Variable constantZero("constantZero", 0);
143 Variable constantMinusOne("constantMinusOne", -1);
144 Variable minDalitzX("minDalitzX", pow(piPlusMass + piZeroMass, 2));
145 Variable maxDalitzX("maxDalitzX", pow(_mD0 - piPlusMass, 2));
146 Variable minDalitzY("minDalitzY", pow(piPlusMass + piZeroMass, 2));
147 Variable maxDalitzY("maxDalitzY", pow(_mD0 - piPlusMass, 2));
148 Variable minDalitzZ("minDalitzZ", pow(piPlusMass + piPlusMass, 2));
149 Variable maxDalitzZ("maxDalitzZ", pow(_mD0 - piZeroMass, 2));
150 
151 std::vector<Variable> weights;
152 std::vector<Observable> obsweights;
153 
154 std::vector<PdfBase *> comps;
155 TH1F *dataTimePlot = nullptr;
156 TH1F *loM23Sigma = nullptr;
157 TH1F *hiM23Sigma = nullptr;
165 GooPdf *sig0_jsugg = nullptr;
166 GooPdf *bkg2_jsugg = nullptr;
167 GooPdf *bkg3_jsugg = nullptr;
168 GooPdf *bkg4_jsugg = nullptr;
169 Variable massSum("massSum", _mD0 *_mD0 + 2 * piPlusMass * piPlusMass + piZeroMass * piZeroMass); // = 3.53481
170 GooPdf *kzero_veto = nullptr;
171 
172 bool doToyStudy = false;
173 float md0_toy_width = 0.0075;
174 float md0_toy_mean = 1.8654;
175 const float toySigFraction = 0.6;
176 const float toyBkgTimeMean = 0.0;
177 const float toyBkgTimeRMS = 0.7;
178 std::string toyFileName;
179 char strbuffer[1000];
180 
181 SmoothHistogramPdf *makeBackgroundHistogram(int bkgnum, std::string overridename = "");
182 void makeToyDalitzPlots(GooPdf *overallSignal, std::string plotdir = "./plots_from_toy_mixfit/");
183 void getBackgroundFile(int bkgType);
184 
185 double intGaus = -1;
186 double calcGauInteg(double x1, double x2) {
187  if(x1 > x2) {
188  double swp = x2;
189  x2 = x1;
190  x1 = swp;
191  }
192 
193  double sum1 = x1, sum2 = x2;
194  double value1 = x1, value2 = x2;
195 
196  for(int i = 1; i < 100; i++) {
197  value1 = value1 * x1 * x1 / (2 * i + 1);
198  sum1 += value1;
199  value2 = value2 * x2 * x2 / (2 * i + 1);
200  sum2 += value2;
201  }
202 
203  return sum2 * exp(-(x2 * x2) / 2) - sum1 * exp(-(x1 * x1) / 2);
204 }
205 
206 double calcToyWeight(double sigratio, double m) {
207  if(intGaus < 0)
208  intGaus = calcGauInteg(0.0075 * md0_lower_window / md0_toy_width, 0.0075 * md0_upper_window / md0_toy_width);
209 
210  double t = (m - md0_toy_mean) / md0_toy_width;
211  double fsig = sigratio * exp(-t * t / 2.) / intGaus;
212  double fbkg = (1 - sigratio) / ((md0_upper_window - md0_lower_window) * 0.0075 / md0_toy_width);
213  return fsig / (fsig + fbkg);
214 }
215 
216 void printMemoryStatus(std::string file, int line);
217 void loadDataFile(std::string fname, UnbinnedDataSet **setToFill = 0, int effSkip = 3);
218 int runBackgroundDalitzFit(int bkgType, bool plots = false);
219 
220 void normalize(TH1F *dat) {
221  double integral = 0;
222 
223  for(int i = 1; i <= dat->GetNbinsX(); ++i) {
224  integral += dat->GetBinContent(i);
225  }
226 
227  integral = 1.0 / integral;
228 
229  for(int i = 1; i <= dat->GetNbinsX(); ++i) {
230  dat->SetBinContent(i, integral * dat->GetBinContent(i));
231  dat->SetBinError(i, integral * dat->GetBinError(i));
232  }
233 }
234 
235 fptype cpuGetM23(fptype massPZ, fptype massPM) {
236  return (_mD02 + piZeroMass * piZeroMass + piPlusMass * piPlusMass + piPlusMass * piPlusMass - massPZ - massPM);
237 }
238 
239 bool cpuDalitz(fptype m12, fptype m13, fptype bigM, fptype dm1, fptype dm2, fptype dm3) {
240  if(m12 < pow(dm1 + dm2, 2))
241  return false; // This m12 cannot exist, it's less than the square of the (1,2) particle mass.
242 
243  if(m12 > pow(bigM - dm3, 2))
244  return false; // This doesn't work either, there's no room for an at-rest 3 daughter.
245 
246  // Calculate energies of 1 and 3 particles in m12 rest frame.
247  fptype e1star = 0.5 * (m12 - dm2 * dm2 + dm1 * dm1) / sqrt(m12);
248  fptype e3star = 0.5 * (bigM * bigM - m12 - dm3 * dm3) / sqrt(m12);
249 
250  // Bounds for m13 at this value of m12.
251  fptype minimum
252  = pow(e1star + e3star, 2) - pow(sqrt(e1star * e1star - dm1 * dm1) + sqrt(e3star * e3star - dm3 * dm3), 2);
253 
254  if(m13 < minimum)
255  return false;
256 
257  fptype maximum
258  = pow(e1star + e3star, 2) - pow(sqrt(e1star * e1star - dm1 * dm1) - sqrt(e3star * e3star - dm3 * dm3), 2);
259 
260  if(m13 > maximum)
261  return false;
262 
263  return true;
264 }
265 
267  if(!loM23Sigma)
268  return;
269 
270  normalize(loM23Sigma);
271  normalize(hiM23Sigma);
272 
273  loM23Sigma->SetLineWidth(3);
274  loM23Sigma->SetLineColor(kBlue);
275  hiM23Sigma->SetLineWidth(3);
276  hiM23Sigma->SetLineColor(kRed);
277 
278  hiM23Sigma->GetXaxis()->SetTitle("#sigma_{t} [ps]");
279  hiM23Sigma->GetYaxis()->SetTitle("Fraction / 8 fs");
280  hiM23Sigma->Draw("l");
281  loM23Sigma->Draw("lsame");
282 
283  TLegend blah(0.7, 0.8, 0.95, 0.95);
284  blah.AddEntry(loM23Sigma, "m_{0}^{2} < 1.5", "l");
285  blah.AddEntry(hiM23Sigma, "m_{0}^{2} > 1.5", "l");
286  blah.Draw();
287 
288  foo->SaveAs("./plots_from_mixfit/sigmacomp.png");
289 }
290 
291 void plotFit(Observable var, UnbinnedDataSet *dat, GooPdf *fit) {
292  int numEvents = dat->getNumEvents();
293  TH1F *dat_hist = new TH1F(
294  (var.getName() + "_dathist").c_str(), "", var.getNumBins(), var.getLowerLimit(), var.getUpperLimit());
295 
296  for(int i = 0; i < numEvents; ++i) {
297  dat_hist->Fill(dat->getValue(var, i));
298  }
299 
300  TH1F *pdf_hist = new TH1F(
301  (var.getName() + "_pdfhist").c_str(), "", var.getNumBins(), var.getLowerLimit(), var.getUpperLimit());
302  std::vector<fptype> values = fit->evaluateAtPoints(var);
303 
304  double totalPdf = 0;
305 
306  for(int i = 0; i < var.getNumBins(); ++i) {
307  totalPdf += values[i];
308  }
309 
310  for(int i = 0; i < var.getNumBins(); ++i) {
311  pdf_hist->SetBinContent(i + 1, values[i] * numEvents / totalPdf);
312  }
313 
314  pdf_hist->SetStats(false);
315  pdf_hist->SetLineColor(kBlue);
316  pdf_hist->SetLineWidth(3);
317 
318  dat_hist->SetStats(false);
319  dat_hist->SetMarkerStyle(8);
320  dat_hist->SetMarkerSize(0.7);
321  dat_hist->Draw("p");
322 
323  pdf_hist->Draw("lsame");
324 
325  foo->SaveAs((var.getName() + "_fit.png").c_str());
326  foo->SetLogy(true);
327  foo->SaveAs((var.getName() + "_logfit.png").c_str());
328  foo->SetLogy(false);
329 }
330 
331 bool readWrapper(std::ifstream &reader, std::string fname = strbuffer) {
332  std::cout << "Now open file " << fname << " for reading" << std::endl;
333  reader.open(fname.c_str());
334  assert(reader.good());
335  return true;
336 }
337 
338 void getToyData(float sigweight = 0.9) {
339  if(!data) {
340  std::vector<Observable> vars;
341  vars.push_back(*m12);
342  vars.push_back(*m13);
343  vars.push_back(*dtime);
344  vars.push_back(*sigma);
345  vars.push_back(*eventNumber);
346  vars.push_back(*wSig0);
347  // vars.push_back(wBkg1);
348  // vars.push_back(wBkg2);
349  data = new UnbinnedDataSet(vars);
350  }
351 
352  std::ifstream reader;
353  readWrapper(reader, toyFileName);
354  std::string buffer;
355 
356  while(!reader.eof()) {
357  reader >> buffer;
358 
359  if(buffer == "====")
360  break;
361 
362  std::cout << buffer;
363  }
364 
365  TRandom3 donram(42);
366 
367  int nsig = 0;
368  double sigprob = 0;
369  double dummy = 0;
370  double md0 = md0_toy_mean;
371 
372  while(reader >> dummy
373 
374  // Ignoring m23, m(pi+ pi-), called m12 in processToyRoot convention.
375  // Already swapped m12 m13 according to D* charge. m12 = m(pi+pi0)
376  >> dummy >> *m12 >> *m13
377 
378  // Errors on Dalitz variables
379  >> dummy >> dummy >> dummy
380 
381  // Remaining two interesting values
382  >> *dtime >> *sigma
383 
384  // Md0, deltaM, ProbSig, Dst charge, Run, Event, Signal and four bkg fractions
385  >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy) {
386  double resolution = donram.Gaus(0, 1);
387  dtime->setValue(dtime->getValue() + resolution * sigma->getValue());
388 
389  eventNumber->setValue(data->getNumEvents());
390 
391  do {
392  md0 = donram.Gaus(md0_toy_mean, md0_toy_width);
393  } while(md0 <= 1.8654 + 0.0075 * md0_lower_window || md0 >= 1.8654 + 0.0075 * md0_upper_window);
394 
395  // wSig0->getValue() = sigweight;
396  wSig0->setValue(calcToyWeight(sigweight, md0));
397  sigprob += wSig0->getValue();
398  data->addEvent();
399  nsig++;
400 
401  if(data->getNumEvents() < 10) {
402  std::cout << data->getNumEvents() << " : " << m12->getValue() << " " << m13->getValue() << " "
403  << dtime->getValue() << " " << sigma->getValue() << " " << std::endl;
404  }
405 
406  if(data->getNumEvents() >= maxEvents)
407  break;
408  }
409 
410  GOOFIT_INFO("Read in {} events", data->getNumEvents())
411 
412  for(int ib = 0; ib < nsig * (1 - sigweight) / sigweight; ib++) {
413  do {
414  m12->setValue(donram.Uniform() * (m12->getUpperLimit() - m12->getLowerLimit()) + m12->getLowerLimit());
415  m13->setValue(donram.Uniform() * (m13->getUpperLimit() - m13->getLowerLimit()) + m13->getLowerLimit());
416  } while(!cpuDalitz(m12->getValue(), m13->getValue(), _mD0, piZeroMass, piPlusMass, piPlusMass));
417 
418  do {
419  dtime->setValue(donram.Gaus(toyBkgTimeMean, toyBkgTimeRMS));
420  } while(!(dtime->getValue() > dtime->getLowerLimit() && dtime->getValue() < dtime->getUpperLimit()));
421 
422  eventNumber->setValue(data->getNumEvents());
423  md0 = donram.Uniform(1.8654 + 0.0075 * md0_lower_window, 1.8654 + 0.0075 * md0_upper_window);
424  // wSig0->getValue() = sigweight;
425  wSig0->setValue(calcToyWeight(sigweight, md0));
426  sigprob += wSig0->getValue();
427  data->addEvent();
428  }
429 
430  reader.close();
431 }
432 
434  Variable effSmoothing("effSmoothing", 1.0, 0.1, 0, 1.25);
435  // Variable* effSmoothing = new Variable("effSmoothing", 0);
436  SmoothHistogramPdf *ret = new SmoothHistogramPdf("efficiency", binEffData, effSmoothing);
437  return ret;
438 }
439 
441  if(kzero_veto)
442  return kzero_veto;
443 
444  VetoInfo kVetoInfo(Variable("veto_min", 0.475 * 0.475), Variable("veto_max", 0.505 * 0.505), PAIR_23);
445 
446  vector<VetoInfo> vetos;
447  vetos.push_back(kVetoInfo);
448  kzero_veto = new DalitzVetoPdf("kzero_veto", *m12, *m13, motherM, neutrlM, chargeM, chargeM, vetos);
449  return kzero_veto;
450 }
451 
453  if(!kzero_veto)
454  makeKzeroVeto();
455 
456  vector<Variable> offsets;
457  vector<Observable> observables;
458  vector<Variable> coefficients;
459  offsets.push_back(constantOne);
460  offsets.push_back(constantOne);
461 
462  observables.push_back(*m12);
463  observables.push_back(*m13);
464 
465  coefficients.push_back(Variable("x0y0", 1.0));
466  coefficients.push_back(Variable("x1y0", 0.07999, 0.01, -0.5, 0.5));
467  coefficients.push_back(Variable("x2y0", -0.23732, 0.01, -0.5, 0.5));
468  coefficients.push_back(Variable("x3y0", 0.10369, 0.01, -1.5, 0.5));
469  coefficients.push_back(Variable("x0y1", 0.10248, 0.01, -1.5, 0.5));
470  coefficients.push_back(Variable("x1y1", -0.28543, 0.01, -0.5, 0.5));
471  coefficients.push_back(Variable("x2y1", 0.15058, 0.01, -1.5, 0.5));
472  coefficients.push_back(Variable("x0y2", -0.20648, 0.01, -0.5, 0.5));
473  coefficients.push_back(Variable("x1y2", 0.14567, 0.01, -1.5, 0.5));
474  coefficients.push_back(Variable("x0y3", 0.06231, 0.01, -0.5, 0.5));
475 
476  PolynomialPdf *poly = new PolynomialPdf("efficiency", observables, coefficients, offsets, 3);
477  poly->setParameterConstantness(true);
478 
479  Variable decXmin("decXmin", 6.22596);
480  Variable decYmin("decYmin", 6.30722);
481  Variable decZmin("decZmin", 10.82390);
482  Variable conXmin("conXmin", 0.65621);
483  Variable conYmin("conYmin", 0.69527);
484  Variable conZmin("conZmin", 0.31764);
485  Variable decXmax("decXmax", 0.79181);
486  Variable decYmax("decYmax", 5.91211);
487  Variable decZmax("decZmax", 1.52031);
488  Variable conXmax("conXmax", 0.80918);
489  Variable conYmax("conYmax", 0.22158);
490  Variable conZmax("conZmax", 0.41866);
491 
492  TrigThresholdPdf *loX = new TrigThresholdPdf("loX", *m12, minDalitzX, decXmin, conXmin, false);
493  TrigThresholdPdf *hiX = new TrigThresholdPdf("hiX", *m12, maxDalitzX, decXmax, conXmax, true);
494 
495  TrigThresholdPdf *loY = new TrigThresholdPdf("loY", *m13, minDalitzY, decYmin, conYmin, false);
496  TrigThresholdPdf *hiY = new TrigThresholdPdf("hiY", *m13, maxDalitzY, decYmax, conYmax, true);
497 
498  TrigThresholdPdf *loZ = new TrigThresholdPdf("loZ", *m12, *m13, minDalitzZ, decZmin, conZmin, massSum, false);
499  TrigThresholdPdf *hiZ = new TrigThresholdPdf("hiZ", *m12, *m13, maxDalitzZ, decZmax, conZmax, massSum, true);
500 
501  comps.clear();
502  comps.push_back(poly);
503  comps.push_back(loX);
504  comps.push_back(hiX);
505  comps.push_back(loY);
506  comps.push_back(hiY);
507  comps.push_back(loZ);
508  comps.push_back(hiZ);
509  comps.push_back(kzero_veto);
510  ProdPdf *ret = new ProdPdf("efficiency_total", comps);
511 
512  // return poly;
513  return ret;
514 }
515 
516 DecayInfo3t dtop0pp{Variable("tau", 0.4101, 0.001, 0.300, 0.500),
517  Variable("xmixing", 0.0016, 0.001, 0, 0),
518  Variable("ymixing", 0.0055, 0.001, 0, 0)};
519 
520 TddpPdf *makeSignalPdf(MixingTimeResolution *resolution = 0, GooPdf *eff = 0) {
526 
527  // dtop0pp._tau = new Variable("tau", 0.4116, 0.001, 0.300, 0.500);
528  // Setting limits causes trouble with MNHESSE - why is this?
529  // dtop0pp._xmixing = new Variable("xmixing", 0.01, 0.001, -0.05, 0.05);
530  // dtop0pp._ymixing = new Variable("ymixing", 0.01, 0.001, -0.05, 0.05);
531 
532  // dtop0pp._tau.fixed = true;
533  // dtop0pp._xmixing.fixed = true;
534  // dtop0pp._ymixing.fixed = true;
535 
536  ResonancePdf *rhop = new Resonances::RBW(
537  "rhop", Variable("rhop_amp_real", 1), Variable("rhop_amp_imag", 0), fixedRhoMass, fixedRhoWidth, 1, PAIR_12);
538 
539  bool fixAmps = false;
540 
541  ResonancePdf *rhom = new Resonances::RBW(
542  "rhom",
543  fixAmps ? Variable("rhom_amp_real", 0.714) : Variable("rhom_amp_real", 0.714, 0.001, 0, 0),
544  fixAmps ? Variable("rhom_amp_imag", -0.025) : Variable("rhom_amp_imag", -0.025, 0.1, 0, 0),
545  fixedRhoMass,
547  1,
548  PAIR_13);
549 
550  ResonancePdf *rho0
551  = new Resonances::GS("rho0",
552  fixAmps ? Variable("rho0_amp_real", 0.565) : Variable("rho0_amp_real", 0.565, 0.001, 0, 0),
553  fixAmps ? Variable("rho0_amp_imag", 0.164) : Variable("rho0_amp_imag", 0.164, 0.1, 0, 0),
554  fixedRhoMass,
556  1,
557  PAIR_23);
558 
559  Variable rho1450Mass("rhop_1450_mass", 1.465, 0.01, 1.0, 2.0);
560  Variable rho1450Width("rhop_1450_width", 0.400, 0.01, 0.01, 5.0);
561 
562  ResonancePdf *rhop_1450 = new Resonances::RBW(
563  "rhop_1450",
564  fixAmps ? Variable("rhop_1450_amp_real", -0.174) : Variable("rhop_1450_amp_real", -0.174, 0.001, 0, 0),
565  fixAmps ? Variable("rhop_1450_amp_imag", -0.117) : Variable("rhop_1450_amp_imag", -0.117, 0.1, 0, 0),
566  rho1450Mass,
567  rho1450Width,
568  1,
569  PAIR_12);
570 
571  ResonancePdf *rho0_1450 = new Resonances::RBW(
572  "rho0_1450",
573  fixAmps ? Variable("rho0_1450_amp_real", 0.325) : Variable("rho0_1450_amp_real", 0.325, 0.001, 0, 0),
574  fixAmps ? Variable("rho0_1450_amp_imag", 0.057) : Variable("rho0_1450_amp_imag", 0.057, 0.1, 0, 0),
575  rho1450Mass,
576  rho1450Width,
577  1,
578  PAIR_23);
579 
580  ResonancePdf *rhom_1450 = new Resonances::RBW(
581  "rhom_1450",
582  fixAmps ? Variable("rhom_1450_amp_real", 0.788) : Variable("rhom_1450_amp_real", 0.788, 0.001, 0, 0),
583  fixAmps ? Variable("rhom_1450_amp_imag", 0.226) : Variable("rhom_1450_amp_imag", 0.226, 0.1, 0, 0),
584  rho1450Mass,
585  rho1450Width,
586  1,
587  PAIR_13);
588 
589  Variable rho1700Mass("rhop_1700_mass", 1.720, 0.01, 1.6, 1.9);
590  Variable rho1700Width("rhop_1700_width", 0.250, 0.01, 0.1, 1.0);
591 
592  ResonancePdf *rhop_1700 = new Resonances::RBW(
593  "rhop_1700",
594  fixAmps ? Variable("rhop_1700_amp_real", 2.151) : Variable("rhop_1700_amp_real", 2.151, 0.001, 0, 0),
595  fixAmps ? Variable("rhop_1700_amp_imag", -0.658) : Variable("rhop_1700_amp_imag", -0.658, 0.1, 0, 0),
596  rho1700Mass,
597  rho1700Width,
598  1,
599  PAIR_12);
600 
601  ResonancePdf *rho0_1700 = new Resonances::RBW(
602  "rho0_1700",
603  fixAmps ? Variable("rho0_1700_amp_real", 2.400) : Variable("rho0_1700_amp_real", 2.400, 0.001, 0, 0),
604  fixAmps ? Variable("rho0_1700_amp_imag", -0.734) : Variable("rho0_1700_amp_imag", -0.734, 0.1, 0, 0),
605  rho1700Mass,
606  rho1700Width,
607  1,
608  PAIR_23);
609 
610  ResonancePdf *rhom_1700 = new Resonances::RBW(
611  "rhom_1700",
612  fixAmps ? Variable("rhom_1700_amp_real", 1.286) : Variable("rhom_1700_amp_real", 1.286, 0.001, 0, 0),
613  fixAmps ? Variable("rhom_1700_amp_imag", -1.532) : Variable("rhom_1700_amp_imag", -1.532, 0.1, 0, 0),
614  rho1700Mass,
615  rho1700Width,
616  1,
617  PAIR_13);
618 
619  ResonancePdf *f0_980 = new Resonances::RBW("f0_980",
620  fixAmps ? Variable("f0_980_amp_real", 0.008 * (-_mD02))
621  : Variable("f0_980_amp_real", 0.008 * (-_mD02), 0.001, 0, 0),
622  fixAmps ? Variable("f0_980_amp_imag", -0.013 * (-_mD02))
623  : Variable("f0_980_amp_imag", -0.013 * (-_mD02), 0.1, 0, 0),
624  Variable("f0_980_mass", 0.980, 0.01, 0.8, 1.2),
625  Variable("f0_980_width", 0.044, 0.001, 0.001, 0.08),
626  0,
627  PAIR_23);
628 
629  ResonancePdf *f0_1370 = new Resonances::RBW("f0_1370",
630  fixAmps ? Variable("f0_1370_amp_real", -0.058 * (-_mD02))
631  : Variable("f0_1370_amp_real", -0.058 * (-_mD02), 0.001, 0, 0),
632  fixAmps ? Variable("f0_1370_amp_imag", 0.026 * (-_mD02))
633  : Variable("f0_1370_amp_imag", 0.026 * (-_mD02), 0.1, 0, 0),
634  Variable("f0_1370_mass", 1.434, 0.01, 1.2, 1.6),
635  Variable("f0_1370_width", 0.173, 0.01, 0.01, 0.4),
636  0,
637  PAIR_23);
638 
639  ResonancePdf *f0_1500 = new Resonances::RBW("f0_1500",
640  fixAmps ? Variable("f0_1500_amp_real", 0.057 * (-_mD02))
641  : Variable("f0_1500_amp_real", 0.057 * (-_mD02), 0.001, 0, 0),
642  fixAmps ? Variable("f0_1500_amp_imag", 0.012 * (-_mD02))
643  : Variable("f0_1500_amp_imag", 0.012 * (-_mD02), 0.1, 0, 0),
644  Variable("f0_1500_mass", 1.507, 0.01, 1.3, 1.7),
645  Variable("f0_1500_width", 0.109, 0.01, 0.01, 0.3),
646  0,
647  PAIR_23);
648 
649  ResonancePdf *f0_1710 = new Resonances::RBW("f0_1710",
650  fixAmps ? Variable("f0_1710_amp_real", 0.070 * (-_mD02))
651  : Variable("f0_1710_amp_real", 0.070 * (-_mD02), 0.001, 0, 0),
652  fixAmps ? Variable("f0_1710_amp_imag", 0.087 * (-_mD02))
653  : Variable("f0_1710_amp_imag", 0.087 * (-_mD02), 0.1, 0, 0),
654  Variable("f0_1710_mass", 1.714, 0.01, 1.5, 2.9),
655  Variable("f0_1710_width", 0.140, 0.01, 0.01, 0.5),
656  0,
657  PAIR_23);
658 
659  ResonancePdf *f2_1270
660  = new Resonances::RBW("f2_1270",
661  fixAmps ? Variable("f2_1270_amp_real", -1.027 * (-_mD02inv))
662  : Variable("f2_1270_amp_real", -1.027 * (-_mD02inv), 0.001, 0, 0),
663  fixAmps ? Variable("f2_1270_amp_imag", -0.162 * (-_mD02inv))
664  : Variable("f2_1270_amp_imag", -0.162 * (-_mD02inv), 0.1, 0, 0),
665  Variable("f2_1270_mass", 1.2754, 0.01, 1.0, 1.5),
666  Variable("f2_1270_width", 0.1851, 0.01, 0.01, 0.4),
667  2,
668  PAIR_23);
669 
670  ResonancePdf *f0_600 = new Resonances::RBW("f0_600",
671  fixAmps ? Variable("f0_600_amp_real", 0.068 * (-_mD02))
672  : Variable("f0_600_amp_real", 0.068 * (-_mD02), 0.001, 0, 0),
673  fixAmps ? Variable("f0_600_amp_imag", 0.010 * (-_mD02))
674  : Variable("f0_600_amp_imag", 0.010 * (-_mD02), 0.1, 0, 0),
675  Variable("f0_600_mass", 0.500, 0.01, 0.3, 0.7),
676  Variable("f0_600_width", 0.400, 0.01, 0.2, 0.6),
677  0,
678  PAIR_23);
679 
680  ResonancePdf *nonr = new Resonances::NonRes(
681  "nonr",
682  fixAmps ? Variable("nonr_amp_real", 0.5595 * (-1)) : Variable("nonr_amp_real", 0.5595 * (-1), 0.001, 0, 0),
683  fixAmps ? Variable("nonr_amp_imag", -0.108761 * (-1)) : Variable("nonr_amp_imag", -0.108761 * (-1), 0.1, 0, 0));
684 
685  dtop0pp.resonances.push_back(nonr);
686  dtop0pp.resonances.push_back(rhop);
687  dtop0pp.resonances.push_back(rho0);
688  dtop0pp.resonances.push_back(rhom);
689 
690  if(!drop_rho_1450) {
691  dtop0pp.resonances.push_back(rhop_1450);
692  dtop0pp.resonances.push_back(rho0_1450);
693  dtop0pp.resonances.push_back(rhom_1450);
694  }
695 
696  if(!drop_rho_1700) {
697  dtop0pp.resonances.push_back(rhop_1700);
698  dtop0pp.resonances.push_back(rho0_1700);
699  dtop0pp.resonances.push_back(rhom_1700);
700  }
701 
702  if(!drop_f0_980)
703  dtop0pp.resonances.push_back(f0_980);
704 
705  if(!drop_f0_1370)
706  dtop0pp.resonances.push_back(f0_1370);
707 
708  if(!drop_f0_1500)
709  dtop0pp.resonances.push_back(f0_1500);
710 
711  if(!drop_f0_1710)
712  dtop0pp.resonances.push_back(f0_1710);
713 
714  if(!drop_f2_1270)
715  dtop0pp.resonances.push_back(f2_1270);
716 
717  if(!drop_f0_600)
718  dtop0pp.resonances.push_back(f0_600);
719 
720  if(!fitMasses) {
721  for(vector<ResonancePdf *>::iterator res = dtop0pp.resonances.begin(); res != dtop0pp.resonances.end(); ++res) {
722  (*res)->setParameterConstantness(true);
723  }
724  }
725 
726  if(!eff) {
727  vector<Variable> offsets;
728  vector<Observable> observables;
729  vector<Variable> coefficients;
730 
731  observables.push_back(*m12);
732  observables.push_back(*m13);
733  offsets.push_back(constantZero);
734  offsets.push_back(constantZero);
735  coefficients.push_back(constantOne);
736  eff = new PolynomialPdf("constantEff", observables, coefficients, offsets, 0);
737  }
738 
739  std::vector<MixingTimeResolution *> resList;
740 
741  if(!resolution) {
742  if(massd0) {
743  for(int i = 0; i < mdslices; ++i) {
744  sprintf(strbuffer, "coreFrac_%i", i);
745  Variable coreFrac(strbuffer, 0.90, 0.001, 0.55, 0.999);
746  sprintf(strbuffer, "coreBias_%i", i);
747  // Variable* coreBias = new Variable(strbuffer, 0.1, 0.001, -0.20, 0.30);
748  Variable coreBias(strbuffer, -0.1, 0.001, -0.20, 0.30);
749  sprintf(strbuffer, "coreScaleFactor_%i", i);
750  Variable coreScaleFactor(strbuffer, 0.96, 0.001, 0.20, 1.50);
751  sprintf(strbuffer, "tailScaleFactor_%i", i);
752  Variable tailScaleFactor(strbuffer, 1.63, 0.001, 0.90, 6.00);
753  resolution = new ThreeGaussResolution(coreFrac,
754  constantOne,
755  coreBias,
756  coreScaleFactor,
757  constantZero,
758  tailScaleFactor,
759  constantZero,
760  constantOne);
761  resList.push_back(resolution);
762  }
763  } else {
764  Variable coreFrac("coreFrac", 0.90, 0.001, 0.35, 0.999);
765  // Variable* tailFrac = new Variable("tailFrac", 0.90, 0.001, 0.80, 1.00);
766  Variable coreBias("coreBias", 0.0, 0.001, -0.20, 0.30);
767  Variable coreScaleFactor("coreScaleFactor", 0.96, 0.001, 0.20, 1.50);
768  // Variable* tailBias = new Variable("tailBias", 0);
769  Variable tailScaleFactor("tailScaleFactor", 1.63, 0.001, 0.90, 6.00);
770  // Variable* outlBias = new Variable("outlBias", 0);
771  // Variable* outlScaleFactor = new Variable("outlScaleFactor", 5.0, 0.01, 0.1, 10.0);
772  // resolution = new ThreeGaussResolution(coreFrac, tailFrac, coreBias, coreScaleFactor, tailBias,
773  // tailScaleFactor, outlBias, outlScaleFactor);
774 
775  // resolution = new ThreeGaussResolution(coreFrac, constantOne, coreBias, coreScaleFactor,
776  // constantZero, tailScaleFactor, constantZero, constantOne);
777  // resolution = new ThreeGaussResolution(coreFrac, constantOne, constantZero, coreScaleFactor,
778  // constantZero, tailScaleFactor, constantZero, constantOne);
779  if(!doToyStudy)
780  resolution = new ThreeGaussResolution(coreFrac,
781  constantOne,
782  coreBias,
783  coreScaleFactor,
784  coreBias,
785  tailScaleFactor,
786  constantZero,
787  constantOne);
788  else {
789  coreBias.setValue(0);
790  coreScaleFactor.setValue(1);
791  coreScaleFactor.setFixed(false);
792  resolution = new ThreeGaussResolution(constantOne,
793  constantOne,
794  coreBias,
795  coreScaleFactor,
796  constantZero,
797  constantOne,
798  constantZero,
799  constantOne);
800  }
801  }
802  }
803 
804  TddpPdf *mixPdf = 0;
805 
806  if(massd0)
807  mixPdf = new TddpPdf("mixPdf", *dtime, *sigma, *m12, *m13, *eventNumber, dtop0pp, resList, eff, *massd0, wBkg1);
808  else
809  mixPdf = new TddpPdf("mixPdf", *dtime, *sigma, *m12, *m13, *eventNumber, dtop0pp, resolution, eff, wBkg1);
810 
811  return mixPdf;
812 }
813 
814 GooPdf *makeFlatBkgDalitzPdf(bool fixem = true) {
815  vector<Variable> offsets;
816  vector<Observable> observables;
817  vector<Variable> coefficients;
818  offsets.push_back(constantZero);
819  offsets.push_back(constantZero);
820  observables.push_back(*m12);
821  observables.push_back(*m13);
822  coefficients.push_back(constantOne);
823 
824  PolynomialPdf *poly = new PolynomialPdf("flatbkgPdf", observables, coefficients, offsets, 0);
825  Variable g_mean("g_mean", toyBkgTimeMean, 0.01, -0.2, 0.5);
826  Variable g_sigma("g_sigma", toyBkgTimeRMS, 0.01, 0.15, 1.5);
827  GooPdf *gt = new GaussianPdf("flatbkg_timepdf", *dtime, g_mean, g_sigma);
828  comps.clear();
829  comps.push_back(poly);
830  comps.push_back(gt);
831  GooPdf *ret = new ProdPdf("flatbkg_total", comps);
832 
833  if(fixem)
834  ret->setParameterConstantness(true);
835 
836  return ret;
837 }
838 
839 int runToyFit(int ifile, int nfile, bool noPlots = true) {
840  if(!nfile || ifile < 0 || ifile >= 100)
841  return 7; // No File or file error
842 
843  doToyStudy = true;
844  // dtime = new Variable("dtime", lowerTime, upperTime);
845  dtime = new Observable("dtime", -3, 5);
846  dtime->setNumBins(floor((upperTime - lowerTime) / 0.05 + 0.5));
847  // dtime->getNumBins() = 200;
848  // sigma = new Variable("sigma", 0, 0.8);
849  sigma = new Observable("sigma", 0.099, 0.101);
850  sigma->setNumBins(1);
851  // Cheating way to avoid Punzi effect for toy MC. The normalisation integral is now a delta function!
852  m12 = new Observable("m12", 0, 3);
853  m12->setNumBins(240);
854  m13 = new Observable("m13", 0, 3);
855  m13->setNumBins(240);
856  eventNumber = new EventNumber("eventNumber", 0, INT_MAX);
857  wSig0 = new Observable("wSig0", 0, 1);
858 
859  for(int i = 0; i < nfile; i++) {
860  // sprintf(strbuffer, "dataFiles/toyPipipi0/dalitz_toyMC_%03d.txt", (i+ifile)%100);
861  sprintf(strbuffer, "dataFiles/toyPipipi0/dalitz_toyMC_%03d.txt", ifile);
862  toyFileName = app_ptr->get_filename(strbuffer, "examples/pipipi0DPFit");
863  getToyData(toySigFraction);
864  }
865 
866  // TruthResolution* dat = new TruthResolution();
867  // TddpPdf* mixPdf = makeSignalPdf(dat);
868  signalDalitz = makeSignalPdf();
869  signalDalitz->setDataSize(data->getNumEvents(), 6); // Default 5 is fine for toys
870  sig0_jsugg = new ExpPdf("sig0_jsugg", *sigma, constantZero);
871  // sig0_jsugg = makeBkg_sigma_strips(0);
872  sig0_jsugg->addSpecialMask(PdfBase::ForceSeparateNorm);
873  sig0_jsugg->setParameterConstantness(true);
874  comps.clear();
875  comps.push_back(signalDalitz);
876  // comps.push_back(sig0_jsugg);
877  std::cout << "Creating overall PDF\n";
878  // ProdPdf* overallSignal = new ProdPdf("overallSignal", comps);
879  GooPdf *bkgPdf = makeFlatBkgDalitzPdf();
880  bkgPdf->setParameterConstantness(true);
881 
882  std::vector<Observable> evtWeights;
883  evtWeights.push_back(*wSig0);
884  // evtWeights.push_back(wBkg2);
885  std::vector<PdfBase *> components;
886  components.push_back(signalDalitz);
887  components.push_back(bkgPdf);
888  EventWeightedAddPdf *mixPdf = new EventWeightedAddPdf("total", evtWeights, components);
889  // GooPdf* mixPdf = overallSignal;
890 
891  mixPdf->setData(data);
892 
893  int retval;
894  if(minuit1) {
895  GooFit::FitManagerMinuit1 datapdf(mixPdf);
896  datapdf.setMaxCalls(64000);
897  datapdf.fit();
898  retval = datapdf;
899  } else {
900  GooFit::FitManagerMinuit2 datapdf(mixPdf);
901  datapdf.setMaxCalls(64000);
902  datapdf.fit();
903  retval = datapdf;
904  }
905 
906  fmt::print("Fit results Toy fit:\n"
907  "tau : ({:.3}) fs\n"
908  "xmixing: ({:.3})\%\n"
909  "ymixing: ({:.3})\%\n",
910  1000 * Uncertain(dtop0pp._tau),
911  100 * Uncertain(dtop0pp._xmixing),
912  100 * Uncertain(dtop0pp._ymixing));
913 
914  if(!noPlots)
915  makeToyDalitzPlots(mixPdf);
916 
917  // makeToyDalitzPlots(signalDalitz);
918  return retval;
919 }
920 
921 void loadDataFile(std::string fname, UnbinnedDataSet **setToFill, int effSkip) {
922  if(!setToFill)
923  setToFill = &data;
924 
925  std::vector<Observable> vars;
926  vars.push_back(*m12);
927  vars.push_back(*m13);
928  vars.push_back(*dtime);
929  vars.push_back(*sigma);
930  vars.push_back(*eventNumber);
931  vars.push_back(*wSig0);
932  vars.push_back(*wBkg1);
933  vars.push_back(*wBkg2);
934  vars.push_back(*wBkg3);
935  vars.push_back(*wBkg4);
936 
937  if(massd0)
938  vars.push_back(*massd0);
939 
940  (*setToFill) = new UnbinnedDataSet(vars);
941  std::ifstream reader;
942  readWrapper(reader, fname.c_str());
943  std::string buffer;
944 
945  while(!reader.eof()) {
946  reader >> buffer;
947 
948  if(buffer == "====")
949  break;
950 
951  std::cout << buffer;
952  }
953 
954  double integralWeights[5] = {0, 0, 0, 0, 0};
955 
956  double dummy = 0;
957  double mass = 0;
958  double delm = 0;
959  int events = 0;
960 
961  while(!reader.eof()) {
962  reader >> dummy;
963 
964  if(reader.eof())
965  break;
966 
967  reader >> dummy; // m23, m(pi+ pi-), called m12 in processToyRoot convention.
968  reader >> *m12; // Already swapped according to D* charge
969  reader >> *m13;
970 
971  // Errors on Dalitz variables
972  reader >> dummy;
973  reader >> dummy;
974  reader >> dummy;
975 
976  reader >> *dtime;
977  reader >> *sigma;
978 
979  if(massd0) {
980  reader >> *massd0; // Md0
981  } else
982  reader >> mass;
983 
984  if(deltam) {
985  reader >> *deltam;
986  } else
987  reader >> delm;
988 
989  reader >> dummy; // ProbSig
990  reader >> dummy; // Dst charge
991  reader >> dummy; // Run
992  reader >> dummy; // Event
993  reader >> *wSig0;
994  reader >> *wBkg1;
995  reader >> *wBkg2;
996  reader >> *wBkg3;
997  reader >> *wBkg4;
998 
999  if(massd0) {
1000  if(massd0->getValue() <= massd0->getLowerLimit())
1001  continue;
1002 
1003  if(massd0->getValue() >= massd0->getUpperLimit())
1004  continue;
1005  } else {
1006  // Enforce signal box on all data sets!
1007  if(mass <= 1.8654 + 0.0075 * md0_lower_window + md0offset)
1008  continue;
1009 
1010  if(mass >= 1.8654 + 0.0075 * md0_upper_window + md0offset)
1011  continue;
1012  }
1013 
1014  if(deltam) {
1015  if(deltam->getValue() >= deltam->getUpperLimit())
1016  continue;
1017 
1018  if(deltam->getValue() <= deltam->getLowerLimit())
1019  continue;
1020  } else {
1021  if(delm >= 0.1454 + 0.0003 * deltam_upper_window)
1022  continue;
1023 
1024  if(delm <= 0.1454 + 0.0003 * deltam_lower_window)
1025  continue;
1026  }
1027 
1028  if(dtime->getValue() < dtime->getLowerLimit())
1029  continue;
1030 
1031  if(dtime->getValue() > dtime->getUpperLimit())
1032  continue;
1033 
1034  if(sigma->getValue() > sigma->getUpperLimit())
1035  continue; // Lower limit is 0, and it can't be lower than that, so whatever.
1036 
1037  integralWeights[0] += wSig0->getValue();
1038  integralWeights[1] += wBkg1->getValue();
1039  integralWeights[2] += wBkg2->getValue();
1040  integralWeights[3] += wBkg3->getValue();
1041  integralWeights[4] += wBkg4->getValue();
1042  eventNumber->setValue((*setToFill)->getNumEvents());
1043 
1044  // See comments in TddpPdf.hh for explanation of this.
1045  double mistag = wSig0->getValue() + wBkg1->getValue() * luckyFrac;
1046  wSig0->setValue(wSig0->getValue() + wBkg1->getValue());
1047  wBkg1->setValue(mistag / wSig0->getValue());
1048 
1049  if((binEffData) && (0 == events % effSkip)) {
1050  double weight = weightHistogram->GetBinContent(weightHistogram->FindBin(m12->getValue(), m13->getValue()));
1051  // weight = 1;
1052  binEffData->addWeightedEvent(weight);
1053 
1054  // binEffData->addEvent();
1055  if(underlyingBins)
1056  underlyingBins->Fill(m12->getValue(), m13->getValue(), weight);
1057  } else
1058  (*setToFill)->addEvent();
1059 
1060  events++;
1061 
1062  if(loM23Sigma) {
1063  double currM23 = cpuGetM23(m12->getValue(), m13->getValue());
1064 
1065  if(currM23 < 1.5)
1066  loM23Sigma->Fill(sigma->getValue());
1067  else
1068  hiM23Sigma->Fill(sigma->getValue());
1069 
1070  // std::cout << "Filled sigma with " << sigma->getValue() << " " << currM23 << std::endl;
1071  }
1072 
1073  if(wSig0->getValue() > 1.0)
1074  std::cout << "Problem with event " << (*setToFill)->getNumEvents() << ", too-large signal weight "
1075  << wSig0->getValue() << std::endl;
1076 
1077  /*
1078  if ((*setToFill)->getNumEvents() < 10) {
1079  std::cout << (*setToFill)->getNumEvents() << " : "
1080  << m12->getValue() << " "
1081  << m13->getValue() << " "
1082  << dtime->getValue() << " "
1083  << sigma->getValue() << " "
1084  << wSig0->getValue() << " "
1085  << wBkg1->getValue() << " "
1086  << wBkg2->getValue() << " "
1087  << wBkg3->getValue() << " "
1088  << wBkg4->getValue() << " "
1089  << std::endl;
1090  }
1091  */
1092  }
1093 
1094  std::cout << "Loaded " << (*setToFill)->getNumEvents() << " events.\n";
1095  std::cout << "Integrals: " << integralWeights[0] << " " << integralWeights[1] << " " << integralWeights[2] << " "
1096  << integralWeights[3] << " " << integralWeights[4] << "\n";
1097 }
1098 
1100  static bool exists = false;
1101 
1102  if(exists)
1103  return;
1104 
1105  exists = true;
1106 
1107  dtime = new Observable("dtime", lowerTime, upperTime);
1108  dtime->setNumBins(floor((upperTime - lowerTime) / 0.05 + 0.5));
1109  sigma = new Observable("sigma", 0, maxSigma);
1110  sigma->setNumBins(floor((maxSigma / 0.01) + 0.5));
1111  m12 = new Observable("m12", 0, 3);
1112  m13 = new Observable("m13", 0, 3);
1113  m12->setNumBins(normBinning);
1114  m13->setNumBins(normBinning);
1115  eventNumber = new EventNumber("eventNumber", 0, INT_MAX);
1116  wSig0 = new Observable("wSig0", 0, 1);
1117  wBkg1 = new Observable("wBkg1", 0, 1);
1118  wBkg2 = new Observable("wBkg2", 0, 1);
1119  wBkg3 = new Observable("wBkg3", 0, 1);
1120  wBkg4 = new Observable("wBkg4", 0, 1);
1121 }
1122 
1123 GooPdf *makeSignalJSU_gg(int idx, bool fixem = true) {
1124  // Values from TSigma, 'Mikhail default'.
1125  static int jsugg_num = -1;
1126  jsugg_num++;
1127 
1128  sprintf(strbuffer, "js_meana_%i", jsugg_num);
1129  Variable js_meana(strbuffer, 0.0593, 0.01, 0, 0.2);
1130  js_meana.setFixed(fixem);
1131  sprintf(strbuffer, "js_sigma_%i", jsugg_num);
1132  Variable js_sigma(strbuffer, 0.000474, 0.0001, 0, 0.001);
1133  js_sigma.setFixed(fixem);
1134  sprintf(strbuffer, "js_gamma_%i", jsugg_num);
1135  Variable js_gamma(strbuffer, -10.1942, 1, -30, 0);
1136  js_gamma.setFixed(fixem);
1137  sprintf(strbuffer, "js_delta_%i", jsugg_num);
1138  Variable js_delta(strbuffer, 1.4907, 0.1, 0.5, 5);
1139  js_delta.setFixed(fixem);
1140  sprintf(strbuffer, "frac_jsu_%i", jsugg_num);
1141  Variable frac_jsu(strbuffer, 0.9516, 0.01, 0.5, 1.0);
1142  frac_jsu.setFixed(fixem);
1143  sprintf(strbuffer, "frac_ga1_%i", jsugg_num);
1144  Variable frac_ga1(strbuffer, 0.001845, 0.0005, 0.0001, 0.3);
1145  frac_ga1.setFixed(fixem);
1146  sprintf(strbuffer, "g1_meana_%i", jsugg_num);
1147  Variable g1_meana(strbuffer, 0.2578, 0.003, 0.1, 0.5);
1148  g1_meana.setFixed(fixem);
1149  sprintf(strbuffer, "g1_sigma_%i", jsugg_num);
1150  Variable g1_sigma(strbuffer, 0.03086, 0.01, 0.005, 0.25);
1151  g1_sigma.setFixed(fixem);
1152  sprintf(strbuffer, "g2_meana_%i", jsugg_num);
1153  Variable g2_meana(strbuffer, 0.32, 0.01, 0.1, 0.5);
1154  g2_meana.setFixed(fixem);
1155  sprintf(strbuffer, "g2_sigma_%i", jsugg_num);
1156  Variable g2_sigma(strbuffer, 0.05825, 0.01, 0.01, 0.1);
1157  g2_sigma.setFixed(fixem);
1158  // Variable* g2_sigma = new Variable("g2_sigma", 0.5825);
1159 
1160  sprintf(strbuffer, "js_%i", jsugg_num);
1161  JohnsonSUPdf *js = new JohnsonSUPdf(strbuffer, *sigma, js_meana, js_sigma, js_gamma, js_delta);
1162  sprintf(strbuffer, "g1_%i", jsugg_num);
1163  GaussianPdf *g1 = new GaussianPdf(strbuffer, *sigma, g1_meana, g1_sigma);
1164  sprintf(strbuffer, "g2_%i", jsugg_num);
1165  // GaussianPdf* g2 = new GaussianPdf(strbuffer, sigma, g2_meana, g2_sigma);
1166 
1167  weights.clear();
1168  weights.push_back(frac_jsu);
1169  // weights.push_back(frac_ga1);
1170  comps.clear();
1171  comps.push_back(js);
1172  comps.push_back(g1);
1173  // comps.push_back(g2);
1174 
1175  // Deal with special indices to get different starting points
1176  switch(idx) {
1177  case 1:
1178  g2_sigma.setUpperLimit(0.15);
1179  break;
1180 
1181  case 2:
1182  frac_jsu.setValue(0.80);
1183  break;
1184 
1185  // case 5:
1186  // g1_sigma->getLowerLimit() = 0.005;
1187  // break;
1188  case 7:
1189  frac_jsu.setValue(0.80);
1190  // frac_ga1->getValue() = 0.05;
1191  break;
1192 
1193  case 11:
1194  frac_ga1.setUpperLimit(0.4);
1195  frac_jsu.setValue(0.80);
1196  break;
1197 
1198  default:
1199  break;
1200  }
1201 
1202  sprintf(strbuffer, "signal_sigma_%i", jsugg_num);
1203  AddPdf *signal_sigma = new AddPdf(strbuffer, weights, comps);
1204  return signal_sigma;
1205 }
1206 
1207 GooPdf *makeMikhailJSU_gg(bool fixem = true) {
1208  // Values from TSigma, 'Mikhail default'.
1209 
1210  Variable *js_meana = new Variable("js_meana", 0.0593279, 0.01, 0, 0.2);
1211  js_meana->setFixed(fixem);
1212  Variable *js_sigma = new Variable("js_sigma", 0.000474171, 0.0001, 0, 0.001);
1213  js_sigma->setFixed(fixem);
1214  Variable *js_gamma = new Variable("js_gamma", -10.1942, 1, -30, 0);
1215  js_gamma->setFixed(fixem);
1216  Variable *js_delta = new Variable("js_delta", 1.4907, 0.1, 0.5, 5);
1217  js_delta->setFixed(fixem);
1218  Variable *frac_jsu = new Variable("frac_jsu", 0.951638, 0.01, 0.5, 1.0);
1219  frac_jsu->setFixed(fixem);
1220  Variable *frac_ga1 = new Variable("frac_ga1", 0.0184522, 0.00001, 0.001, 0.3);
1221  frac_ga1->setFixed(fixem);
1222  Variable *g1_meana = new Variable("g1_meana", 0.257796, 0.003, 0.1, 0.5);
1223  g1_meana->setFixed(fixem);
1224  Variable *g1_sigma = new Variable("g1_sigma", 0.0308619, 0.01, 0.005, 0.25);
1225  g1_sigma->setFixed(fixem);
1226  Variable *g2_meana = new Variable("g2_meana", 0.319952, 0.01, 0.1, 0.5);
1227  g2_meana->setFixed(fixem);
1228  Variable *g2_sigma = new Variable("g2_sigma", 0.0582524, 0.01, 0.01, 0.1);
1229  g2_sigma->setFixed(fixem);
1230 
1231  // Variable* g2_sigma = new Variable("g2_sigma", 0.5825);
1232 
1233  frac_ga1->setValue(frac_ga1->getValue() * (1 - frac_jsu->getValue()));
1234 
1235  JohnsonSUPdf *js = new JohnsonSUPdf("js", *sigma, *js_meana, *js_sigma, *js_gamma, *js_delta);
1236  GaussianPdf *g1 = new GaussianPdf("g1", *sigma, *g1_meana, *g1_sigma);
1237  GaussianPdf *g2 = new GaussianPdf("g2", *sigma, *g2_meana, *g2_sigma);
1238 
1239  weights.clear();
1240  weights.push_back(*frac_jsu);
1241  weights.push_back(*frac_ga1);
1242  comps.clear();
1243  comps.push_back(js);
1244  comps.push_back(g1);
1245  comps.push_back(g2);
1246 
1247  GooPdf *ret = new AddPdf("signal_sigma", weights, comps);
1248  return ret;
1249 }
1250 
1251 const int numSigmaBins = 36;
1252 TH1F **sigma_dat_hists = 0;
1253 TH1F **sigma_pdf_hists = 0;
1255 vector<GooPdf *> jsuList;
1256 
1258  sigma_dat_hists = new TH1F *[numSigmaBins];
1259  sigma_pdf_hists = new TH1F *[numSigmaBins];
1260  sigma_data = new UnbinnedDataSet *[numSigmaBins];
1261 
1262  std::vector<Observable> vars;
1263  vars.push_back(*sigma);
1264 
1265  for(int i = 0; i < numSigmaBins; ++i) {
1266  sprintf(strbuffer, "sigma_data_hist_%i", i);
1267  sigma_dat_hists[i]
1268  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1269  sigma_dat_hists[i]->SetStats(false);
1270  sigma_dat_hists[i]->SetMarkerStyle(8);
1271  sigma_dat_hists[i]->SetMarkerSize(0.7);
1272  sprintf(strbuffer, "sigma_pdf_hist_%i", i);
1273  sigma_pdf_hists[i]
1274  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1275  sigma_pdf_hists[i]->SetStats(false);
1276  sigma_pdf_hists[i]->SetLineWidth(3);
1277  sigma_pdf_hists[i]->SetLineColor(kBlue);
1278 
1279  sigma_data[i] = new UnbinnedDataSet(vars);
1280  }
1281 
1282  int numEvents = data->getNumEvents();
1283 
1284  for(int i = 0; i < numEvents; ++i) {
1285  m12->setValue(data->getValue(*m12, i));
1286  m13->setValue(data->getValue(*m13, i));
1287  sigma->setValue(data->getValue(*sigma, i));
1288 
1289  int xbin = (int)floor(m12->getValue() / 0.5);
1290  int ybin = (int)floor(m13->getValue() / 0.5);
1291  int overallbin = ybin * 6 + xbin;
1292  sigma_dat_hists[overallbin]->Fill(sigma->getValue());
1293  sigma_data[overallbin]->addEvent();
1294  }
1295 
1296  // vector<GooPdf*> jsuList;
1297  for(int i = 0; i < numSigmaBins; ++i) {
1298  GooPdf *js = makeSignalJSU_gg(i, false);
1299  jsuList.push_back(js);
1300 
1301  // int xbin = i % 6;
1302  // int ybin = i / 6;
1303 
1304  /*int inDalitzPlot = 0;
1305  if (cpuDalitz(0.5*(xbin+0), 0.5*(ybin+0), _mD0, piZeroMass, piPlusMass, piPlusMass)) inDalitzPlot++;
1306  if (cpuDalitz(0.5*(xbin+1), 0.5*(ybin+0), _mD0, piZeroMass, piPlusMass, piPlusMass)) inDalitzPlot++;
1307  if (cpuDalitz(0.5*(xbin+0), 0.5*(ybin+1), _mD0, piZeroMass, piPlusMass, piPlusMass)) inDalitzPlot++;
1308  if (cpuDalitz(0.5*(xbin+1), 0.5*(ybin+1), _mD0, piZeroMass, piPlusMass, piPlusMass)) inDalitzPlot++;
1309  */
1310  if(0 == sigma_data[i]->getNumEvents())
1311  js->setParameterConstantness(true);
1312  else {
1313  std::cout << "\n\nAbout to start fit of sigma box " << i << std::endl;
1314  js->setData(sigma_data[i]);
1315  if(minuit1) {
1316  GooFit::FitManagerMinuit1 currpdf(js);
1317  currpdf.fit();
1318  } else {
1319  GooFit::FitManagerMinuit2 currpdf(js);
1320  currpdf.fit();
1321  }
1322  js->setParameterConstantness(true);
1323  // js->clearCurrentFit();
1324  std::cout << "Done with sigma box " << i << "\n";
1325  }
1326  }
1327 
1328  vector<Observable> obses;
1329  obses.push_back(*m12);
1330  obses.push_back(*m13);
1331  vector<double> limits;
1332  limits.push_back(0);
1333  limits.push_back(0);
1334  vector<double> binSizes;
1335  binSizes.push_back(0.5);
1336  binSizes.push_back(0.5);
1337  vector<int> numBins;
1338  numBins.push_back(6);
1339  numBins.push_back(6);
1340  BinTransformPdf *mapFunction = new BinTransformPdf("mapFunction", obses, limits, binSizes, numBins);
1341 
1342  return new MappedPdf("sigmaMap", mapFunction, jsuList);
1343 }
1344 
1346  sigma_dat_hists = new TH1F *[1];
1347  sigma_pdf_hists = new TH1F *[1];
1348  sigma_data = new UnbinnedDataSet *[1];
1349 
1350  std::vector<Observable> vars;
1351  vars.push_back(*sigma);
1352 
1353  for(int i = 0; i < 1; ++i) {
1354  sprintf(strbuffer, "sigma_data_hist_%i", i);
1355  sigma_dat_hists[i]
1356  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1357  sigma_dat_hists[i]->SetStats(false);
1358  sigma_dat_hists[i]->SetMarkerStyle(8);
1359  sigma_dat_hists[i]->SetMarkerSize(0.7);
1360  sprintf(strbuffer, "sigma_pdf_hist_%i", i);
1361  sigma_pdf_hists[i]
1362  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1363  sigma_pdf_hists[i]->SetStats(false);
1364  sigma_pdf_hists[i]->SetLineWidth(3);
1365  sigma_pdf_hists[i]->SetLineColor(kBlue);
1366 
1367  sigma_data[i] = new UnbinnedDataSet(vars);
1368  }
1369 
1370  int numEvents = data->getNumEvents();
1371 
1372  for(int i = 0; i < numEvents; ++i) {
1373  m12->setValue(data->getValue(*m12, i));
1374  m13->setValue(data->getValue(*m13, i));
1375  sigma->setValue(data->getValue(*sigma, i));
1376 
1377  int overallbin = 0;
1378  sigma_dat_hists[overallbin]->Fill(sigma->getValue());
1379  sigma_data[overallbin]->addEvent();
1380  }
1381 
1382  // vector<GooPdf*> jsuList;
1383  for(int i = 0; i < 1; ++i) {
1384  GooPdf *js = makeSignalJSU_gg(i, false);
1385  jsuList.push_back(js);
1386 
1387  if(0 == sigma_data[i]->getNumEvents())
1388  js->setParameterConstantness(true);
1389  else {
1390  std::cout << "\n\nAbout to start fit of sigma box " << i << std::endl;
1391  js->setData(sigma_data[i]);
1392  if(minuit1) {
1393  GooFit::FitManagerMinuit1 currpdf(js);
1394  currpdf.fit();
1395  } else {
1396  GooFit::FitManagerMinuit2 currpdf(js);
1397  currpdf.fit();
1398  }
1399  js->setParameterConstantness(true);
1400  // js->clearCurrentFit();
1401  std::cout << "Done with sigma box " << i << "\n";
1402  }
1403  }
1404 
1405  vector<Observable> obses;
1406  obses.push_back(*m12);
1407  obses.push_back(*m13);
1408  vector<double> limits;
1409  limits.push_back(0);
1410  limits.push_back(0);
1411  vector<double> binSizes;
1412  binSizes.push_back(3);
1413  binSizes.push_back(3);
1414  vector<int> numBins;
1415  numBins.push_back(1);
1416  numBins.push_back(1);
1417  BinTransformPdf *mapFunction = new BinTransformPdf("mapFunction", obses, limits, binSizes, numBins);
1418 
1419  return new MappedPdf("sigmaMap", mapFunction, jsuList);
1420 }
1421 
1423  sigma_dat_hists = new TH1F *[4];
1424  sigma_pdf_hists = new TH1F *[4];
1425  sigma_data = new UnbinnedDataSet *[4];
1426 
1427  std::vector<Observable> vars;
1428  vars.push_back(*sigma);
1429 
1430  for(int i = 0; i < 4; ++i) {
1431  sprintf(strbuffer, "sigma_data_hist_%i", i);
1432  sigma_dat_hists[i]
1433  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1434  sigma_dat_hists[i]->SetStats(false);
1435  sigma_dat_hists[i]->SetMarkerStyle(8);
1436  sigma_dat_hists[i]->SetMarkerSize(0.7);
1437  sprintf(strbuffer, "sigma_pdf_hist_%i", i);
1438  sigma_pdf_hists[i]
1439  = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1440  sigma_pdf_hists[i]->SetStats(false);
1441  sigma_pdf_hists[i]->SetLineWidth(3);
1442  sigma_pdf_hists[i]->SetLineColor(kBlue);
1443 
1444  sigma_data[i] = new UnbinnedDataSet(vars);
1445  }
1446 
1447  int numEvents = data->getNumEvents();
1448 
1449  for(int i = 0; i < numEvents; ++i) {
1450  m12->setValue(data->getValue(*m12, i));
1451  m13->setValue(data->getValue(*m13, i));
1452  sigma->setValue(data->getValue(*sigma, i));
1453 
1454  int xbin = (int)floor(m12->getValue() / 1.5);
1455  int ybin = (int)floor(m13->getValue() / 1.5);
1456  int overallbin = ybin * 2 + xbin;
1457  sigma_dat_hists[overallbin]->Fill(sigma->getValue());
1458  sigma_data[overallbin]->addEvent();
1459  }
1460 
1461  // vector<GooPdf*> jsuList;
1462  for(int i = 0; i < 4; ++i) {
1463  GooPdf *js = makeSignalJSU_gg(i, false);
1464  jsuList.push_back(js);
1465 
1466  if(0 == sigma_data[i]->getNumEvents())
1467  js->setParameterConstantness(true);
1468  else {
1469  std::cout << "\n\nAbout to start fit of sigma box " << i << std::endl;
1470  js->setData(sigma_data[i]);
1471  if(minuit1) {
1472  GooFit::FitManagerMinuit1 currpdf(js);
1473  currpdf.fit();
1474  } else {
1475  GooFit::FitManagerMinuit2 currpdf(js);
1476  currpdf.fit();
1477  }
1478  js->setParameterConstantness(true);
1479  // js->clearCurrentFit();
1480  std::cout << "Done with sigma box " << i << "\n";
1481  }
1482  }
1483 
1484  vector<Observable> obses;
1485  obses.push_back(*m12);
1486  obses.push_back(*m13);
1487  vector<double> limits;
1488  limits.push_back(0);
1489  limits.push_back(0);
1490  vector<double> binSizes;
1491  binSizes.push_back(1.5);
1492  binSizes.push_back(1.5);
1493  vector<int> numBins;
1494  numBins.push_back(2);
1495  numBins.push_back(2);
1496  BinTransformPdf *mapFunction = new BinTransformPdf("mapFunction", obses, limits, binSizes, numBins);
1497 
1498  return new MappedPdf("sigmaMap", mapFunction, jsuList);
1499 }
1500 
1501 void coarseBin(TH2F &dalitzHist, int grain) {
1502  // Move from too-fine-to-see binning down to reasonable binning in Dalitz plots.
1503  for(int i = 1; i < m12->getNumBins(); i += grain) {
1504  for(int j = 1; j < m13->getNumBins(); j += grain) {
1505  double total = 0;
1506 
1507  for(int k = 0; k < grain; ++k) {
1508  for(int l = 0; l < grain; ++l) {
1509  total += dalitzHist.GetBinContent(i + k, j + l);
1510  }
1511  }
1512 
1513  for(int k = 0; k < grain; ++k) {
1514  for(int l = 0; l < grain; ++l) {
1515  dalitzHist.SetBinContent(i + k, j + l, total);
1516  }
1517  }
1518  }
1519  }
1520 }
1521 
1522 struct BigBin {
1523  int xbin;
1524  int ybin;
1525  int width;
1526  int height;
1527  double getContent(TH2F *plot);
1528 };
1529 
1530 double BigBin::getContent(TH2F *plot) {
1531  double ret = 0;
1532 
1533  // std::cout << "getContent with " << width << " " << height << " " << xbin << " " << ybin <<std::endl;
1534  for(unsigned int i = 0; i < width; ++i) {
1535  for(unsigned int j = 0; j < height; ++j) {
1536  // std::cout << i << ", " << j << std::endl;
1537  if(xbin + i > plot->GetNbinsX())
1538  continue;
1539 
1540  if(ybin + j > plot->GetNbinsY())
1541  continue;
1542 
1543  ret += plot->GetBinContent(xbin + i, ybin + j);
1544  }
1545  }
1546 
1547  // std::cout << "Total " << ret << std::endl;
1548  return ret;
1549 }
1550 
1551 struct ChisqInfo {
1552  ChisqInfo();
1553  double chisq;
1554  int dof;
1556 };
1557 
1559  : chisq(0)
1560  , dof(0)
1561  , contribPlot(0) {}
1562 
1563 ChisqInfo *getAdaptiveChisquare(TH2F *datPlot, TH2F *pdfPlot) {
1564  bool acceptable = false;
1565  int binSize = 1;
1566  vector<BigBin> binlist;
1567  double limit = 26;
1568 
1569  while(!acceptable) {
1570  binlist.clear();
1571  std::cout << "Attempting bin generation with size " << binSize << std::endl;
1572 
1573  for(int xbin = 1; xbin <= datPlot->GetNbinsX(); xbin += binSize) {
1574  for(int ybin = 1; ybin <= datPlot->GetNbinsY(); ybin += binSize) {
1575  double lox = datPlot->GetXaxis()->GetBinLowEdge(xbin + 0);
1576  double hix = datPlot->GetXaxis()->GetBinLowEdge(xbin + 1 + binSize);
1577  double loy = datPlot->GetYaxis()->GetBinLowEdge(ybin + 0);
1578  double hiy = datPlot->GetYaxis()->GetBinLowEdge(ybin + 1 + binSize);
1579  bool corner = false;
1580 
1581  if(cpuDalitz(lox, loy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1582  corner = true;
1583  else if(cpuDalitz(lox, hiy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1584  corner = true;
1585  else if(cpuDalitz(hix, loy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1586  corner = true;
1587  else if(cpuDalitz(hix, hiy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1588  corner = true;
1589 
1590  if(!corner)
1591  continue;
1592 
1593  BigBin curr;
1594  curr.xbin = xbin;
1595  curr.ybin = ybin;
1596  curr.width = curr.height = binSize;
1597  binlist.push_back(curr);
1598  }
1599  }
1600 
1601  acceptable = true;
1602 
1603  for(vector<BigBin>::iterator bin = binlist.begin(); bin != binlist.end(); ++bin) {
1604  if((*bin).getContent(datPlot) >= limit)
1605  continue;
1606 
1607  acceptable = false;
1608  binSize *= 2;
1609  std::cout << "Couldn't get good bins, retry.\n";
1610  break;
1611  }
1612  }
1613 
1614  std::cout << "Good bins at size " << binSize << ", beginning splits.\n";
1615 
1616  // Now attempt to split bins.
1617  int numSplits = 1;
1618 
1619  while(0 < numSplits) {
1620  numSplits = 0;
1621  vector<BigBin> newbins;
1622 
1623  for(vector<BigBin>::iterator bin = binlist.begin(); bin != binlist.end(); ++bin) {
1624  if(1 == (*bin).width * (*bin).height) {
1625  newbins.push_back(*bin);
1626  continue;
1627  }
1628 
1629  BigBin lolef;
1630  BigBin lorig;
1631  BigBin hilef;
1632  BigBin hirig;
1633  lolef.xbin = (*bin).xbin;
1634  hilef.xbin = (*bin).xbin;
1635  lorig.xbin = (*bin).xbin + (*bin).width / 2;
1636  hirig.xbin = (*bin).xbin + (*bin).width / 2;
1637  lolef.ybin = (*bin).ybin;
1638  hilef.ybin = (*bin).ybin + (*bin).height / 2;
1639  lorig.ybin = (*bin).ybin;
1640  hirig.ybin = (*bin).ybin + (*bin).height / 2;
1641 
1642  lolef.width = (*bin).width / 2;
1643  lorig.width = (*bin).width / 2;
1644  hilef.width = (*bin).width / 2;
1645  hirig.width = (*bin).width / 2;
1646  lolef.height = (*bin).height / 2;
1647  lorig.height = (*bin).height / 2;
1648  hilef.height = (*bin).height / 2;
1649  hirig.height = (*bin).height / 2;
1650 
1651  int mask = 0;
1652 
1653  if(limit < lolef.getContent(datPlot))
1654  mask += 1;
1655 
1656  if(limit < lorig.getContent(datPlot))
1657  mask += 2;
1658 
1659  if(limit < hilef.getContent(datPlot))
1660  mask += 4;
1661 
1662  if(limit < hirig.getContent(datPlot))
1663  mask += 8;
1664 
1665  if(mask != 15) {
1666  newbins.push_back(*bin);
1667  } else {
1668  newbins.push_back(lolef);
1669  newbins.push_back(lorig);
1670  newbins.push_back(hilef);
1671  newbins.push_back(hirig);
1672  numSplits++;
1673  }
1674  }
1675 
1676  binlist.clear();
1677 
1678  for(vector<BigBin>::iterator i = newbins.begin(); i != newbins.end(); ++i)
1679  binlist.push_back(*i);
1680 
1681  std::cout << "Split " << numSplits << " bins.\n";
1682  }
1683 
1684  ChisqInfo *ret = new ChisqInfo();
1685  ret->dof = binlist.size();
1686  ret->contribPlot = new TH2F("contribPlot",
1687  "",
1688  datPlot->GetNbinsX(),
1689  datPlot->GetXaxis()->GetBinLowEdge(1),
1690  datPlot->GetXaxis()->GetBinLowEdge(datPlot->GetNbinsX() + 1),
1691  datPlot->GetNbinsY(),
1692  datPlot->GetYaxis()->GetBinLowEdge(1),
1693  datPlot->GetYaxis()->GetBinLowEdge(datPlot->GetNbinsY() + 1));
1694 
1695  double totalDat = 0;
1696  double totalPdf = 0;
1697 
1698  for(vector<BigBin>::iterator bin = binlist.begin(); bin != binlist.end(); ++bin) {
1699  double dat = (*bin).getContent(datPlot);
1700  double pdf = (*bin).getContent(pdfPlot);
1701  double term = (dat - pdf) / sqrt(dat);
1702  ret->chisq += term * term;
1703 
1704  /*
1705  std::cout << "Bin (" << (*bin).xbin << ", " << (*bin).ybin << ") "
1706  << (*bin).width << " " << (*bin).height << " : "
1707  << dat << " " << pdf << " "
1708  << term << std::endl;
1709  */
1710  for(int i = 0; i < (*bin).width; ++i) {
1711  for(int j = 0; j < (*bin).height; ++j) {
1712  bool corner = false;
1713  double lox = datPlot->GetXaxis()->GetBinLowEdge((*bin).xbin + i);
1714  double hix = datPlot->GetXaxis()->GetBinLowEdge((*bin).xbin + i + 1);
1715  double loy = datPlot->GetYaxis()->GetBinLowEdge((*bin).ybin + j);
1716  double hiy = datPlot->GetYaxis()->GetBinLowEdge((*bin).ybin + j + 1);
1717 
1718  if(cpuDalitz(lox, loy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1719  corner = true;
1720  else if(cpuDalitz(lox, hiy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1721  corner = true;
1722  else if(cpuDalitz(hix, loy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1723  corner = true;
1724  else if(cpuDalitz(hix, hiy, _mD0, piZeroMass, piPlusMass, piPlusMass))
1725  corner = true;
1726 
1727  if(!corner)
1728  continue;
1729 
1730  ret->contribPlot->SetBinContent((*bin).xbin + i, (*bin).ybin + j, term);
1731  }
1732  }
1733 
1734  totalPdf += pdf;
1735  totalDat += dat;
1736  }
1737 
1738  return ret;
1739 }
1740 
1741 void makeToyDalitzPlots(GooPdf *overallSignal, std::string plotdir) {
1742  std::string call = "mkdir -p " + plotdir;
1743  system(call.c_str());
1744 
1745  foo->cd();
1746 
1747  TH1F dtime_dat_hist("dtime_dat_hist", "", dtime->getNumBins(), dtime->getLowerLimit(), dtime->getUpperLimit());
1748  dtime_dat_hist.SetStats(false);
1749  dtime_dat_hist.SetMarkerStyle(8);
1750  dtime_dat_hist.SetMarkerSize(1.2);
1751  dtime_dat_hist.GetXaxis()->SetTitle("Decay time [ps]");
1752  dtime_dat_hist.GetYaxis()->SetTitle("Events / 50 fs");
1753  TH1F dtime_pdf_hist("dtime_pdf_hist", "", dtime->getNumBins(), dtime->getLowerLimit(), dtime->getUpperLimit());
1754  dtime_pdf_hist.SetStats(false);
1755  dtime_pdf_hist.SetLineColor(kBlue);
1756  dtime_pdf_hist.SetLineWidth(3);
1757  TH1F dtime_sig_hist("dtime_sig_hist", "", dtime->getNumBins(), dtime->getLowerLimit(), dtime->getUpperLimit());
1758  dtime_sig_hist.SetStats(false);
1759  dtime_sig_hist.SetLineColor(kRed);
1760  dtime_sig_hist.SetLineWidth(3);
1761  TH1F dtime_bg_hist("dtime_bg_hist", "", dtime->getNumBins(), dtime->getLowerLimit(), dtime->getUpperLimit());
1762  dtime_bg_hist.SetStats(false);
1763  dtime_bg_hist.SetLineColor(kMagenta);
1764  dtime_bg_hist.SetLineWidth(3);
1765 
1766  double totalPdf = 0;
1767  double totalPdf_sig = 0;
1768  double totalPdf_bg = 0;
1769  double totalDat = 0;
1770  double totalSigProb = 0;
1771  double totalBGProb = 0;
1772 
1773  for(unsigned int evt = 0; evt < data->getNumEvents(); ++evt) {
1774  double currTime = data->getValue(*dtime, evt);
1775  dtime_dat_hist.Fill(currTime);
1776  totalSigProb += data->getValue(*wSig0, evt);
1777  totalBGProb += 1 - data->getValue(*wSig0, evt);
1778  totalDat++;
1779  }
1780 
1781  std::cout << "totalData = " << totalDat << ", totalSigProb = " << totalSigProb << std::endl;
1782  std::vector<Observable> vars;
1783  vars.push_back(*m12);
1784  vars.push_back(*m13);
1785  vars.push_back(*dtime);
1786  vars.push_back(*sigma);
1787  vars.push_back(*eventNumber);
1788  vars.push_back(*wSig0);
1789  UnbinnedDataSet currData(vars);
1790  sigma->setValue(0.1);
1791  wSig0->setValue(totalSigProb / totalDat);
1792  int evtCounter = 0;
1793 
1794  for(int i = 0; i < m12->getNumBins(); ++i) {
1795  m12->setValue(m12->getLowerLimit()
1796  + (m12->getUpperLimit() - m12->getLowerLimit()) * (i + 0.5) / m12->getNumBins());
1797 
1798  for(int j = 0; j < m13->getNumBins(); ++j) {
1799  m13->setValue(m13->getLowerLimit()
1800  + (m13->getUpperLimit() - m13->getLowerLimit()) * (j + 0.5) / m13->getNumBins());
1801 
1803  continue;
1804 
1805  for(int l = 0; l < dtime->getNumBins(); ++l) {
1806  dtime->setValue(dtime->getLowerLimit()
1807  + (dtime->getUpperLimit() - dtime->getLowerLimit()) * (l + 0.5) / dtime->getNumBins());
1808  eventNumber->setValue(evtCounter);
1809  evtCounter++;
1810  currData.addEvent();
1811  }
1812  }
1813  }
1814 
1815  GOOFIT_INFO("Adding {} signal events from toy", currData.getNumEvents());
1816 
1817  overallSignal->setData(&currData);
1818  signalDalitz->setDataSize(currData.getNumEvents(), 6);
1819  std::vector<std::vector<double>> pdfValues = overallSignal->getCompProbsAtDataPoints();
1820 
1821  for(unsigned int j = 0; j < pdfValues[0].size(); ++j) {
1822  double currTime = currData.getValue(*dtime, j);
1823  dtime_sig_hist.Fill(currTime, pdfValues[1][j]);
1824  dtime_bg_hist.Fill(currTime, pdfValues[2][j]);
1825  totalPdf += pdfValues[0][j];
1826  totalPdf_sig += pdfValues[1][j];
1827  totalPdf_bg += pdfValues[2][j];
1828  }
1829 
1830  for(int i = 1; i <= dtime->getNumBins(); ++i) {
1831  dtime_sig_hist.SetBinContent(i, dtime_sig_hist.GetBinContent(i) * totalSigProb / totalPdf_sig);
1832  dtime_bg_hist.SetBinContent(i, dtime_bg_hist.GetBinContent(i) * totalBGProb / totalPdf_bg);
1833  dtime_pdf_hist.SetBinContent(i, dtime_sig_hist.GetBinContent(i) + dtime_bg_hist.GetBinContent(i));
1834  }
1835 
1836  foo->cd();
1837  dtime_dat_hist.Draw("p");
1838  dtime_pdf_hist.Draw("lsame");
1839  dtime_bg_hist.SetLineStyle(2);
1840  dtime_bg_hist.Draw("lsame");
1841  dtime_sig_hist.SetLineStyle(3);
1842  dtime_sig_hist.Draw("lsame");
1843 
1844  foo->SaveAs((plotdir + "/dtime_fit.png").c_str());
1845  foo->SetLogy(true);
1846  foo->SaveAs((plotdir + "/dtime_fit_log.png").c_str());
1847  foo->SetLogy(false);
1848 }
1849 
1850 void makeDalitzPlots(GooPdf *overallSignal, std::string plotdir = "./plots_from_mixfit/") {
1851  std::string mkplotdir{"mkdir " + plotdir};
1852  system(mkplotdir.c_str());
1853  foo->cd();
1854 
1855  TH1F dtime_dat_hist("dtime_dat_hist", "", dtime->getNumBins(), dtime->getLowerLimit(), dtime->getUpperLimit());
1856  dtime_dat_hist.SetStats(false);
1857  dtime_dat_hist.SetMarkerStyle(8);
1858  dtime_dat_hist.SetMarkerSize(1.2);
1859  dtime_dat_hist.GetXaxis()->SetTitle("Decay time [ps]");
1860  dtime_dat_hist.GetYaxis()->SetTitle("Events / 50 fs");
1861  TH1F dtime_pdf_hist("dtime_pdf_hist", "", dtime->getNumBins(), dtime->getLowerLimit(), dtime->getUpperLimit());
1862  dtime_pdf_hist.SetStats(false);
1863  dtime_pdf_hist.SetLineColor(kBlue);
1864  dtime_pdf_hist.SetLineWidth(3);
1865 
1866  TH1F sigma_dat_hist("sigma_dat_hist", "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1867  sigma_dat_hist.SetStats(false);
1868  sigma_dat_hist.SetMarkerStyle(8);
1869  sigma_dat_hist.SetMarkerSize(1.2);
1870  sigma_dat_hist.GetXaxis()->SetTitle("Decay time error [ps]");
1871  sigma_dat_hist.GetYaxis()->SetTitle("Events / 8 fs");
1872  TH1F sigma_pdf_hist("sigma_pdf_hist", "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1873  sigma_pdf_hist.SetStats(false);
1874  sigma_pdf_hist.SetLineColor(kBlue);
1875  sigma_pdf_hist.SetLineWidth(3);
1876 
1877  TH1F m12_dat_hist("m12_dat_hist", "", m12->getNumBins(), m12->getLowerLimit(), m12->getUpperLimit());
1878  m12_dat_hist.SetStats(false);
1879  m12_dat_hist.SetMarkerStyle(8);
1880  m12_dat_hist.SetMarkerSize(1.2);
1881  m12_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{+} #pi^{0}) [GeV]");
1882  m12_dat_hist.GetYaxis()->SetTitle("Events / 12.5 MeV");
1883  TH1F m12_pdf_hist("m12_pdf_hist", "", m12->getNumBins(), m12->getLowerLimit(), m12->getUpperLimit());
1884  m12_pdf_hist.SetStats(false);
1885  m12_pdf_hist.SetLineColor(kBlue);
1886  m12_pdf_hist.SetLineWidth(3);
1887 
1888  TH1F m13_dat_hist("m13_dat_hist", "", m13->getNumBins(), m13->getLowerLimit(), m13->getUpperLimit());
1889  m13_dat_hist.SetStats(false);
1890  m13_dat_hist.SetMarkerStyle(8);
1891  m13_dat_hist.SetMarkerSize(1.2);
1892  m13_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{-} #pi^{0}) [GeV]");
1893  m13_dat_hist.GetYaxis()->SetTitle("Events / 12.5 MeV");
1894  TH1F m13_pdf_hist("m13_pdf_hist", "", m13->getNumBins(), m13->getLowerLimit(), m13->getUpperLimit());
1895  m13_pdf_hist.SetStats(false);
1896  m13_pdf_hist.SetLineColor(kBlue);
1897  m13_pdf_hist.SetLineWidth(3);
1898 
1899  TH1F m23_dat_hist("m23_dat_hist", "", m13->getNumBins(), m13->getLowerLimit(), m13->getUpperLimit());
1900  m23_dat_hist.SetStats(false);
1901  m23_dat_hist.SetMarkerStyle(8);
1902  m23_dat_hist.SetMarkerSize(1.2);
1903  m23_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{-} #pi^{+}) [GeV]");
1904  m23_dat_hist.GetYaxis()->SetTitle("Events / 12.5 MeV");
1905  TH1F m23_pdf_hist("m23_pdf_hist", "", m13->getNumBins(), m13->getLowerLimit(), m13->getUpperLimit());
1906  m23_pdf_hist.SetStats(false);
1907  m23_pdf_hist.SetLineColor(kBlue);
1908  m23_pdf_hist.SetLineWidth(3);
1909 
1910  TH2F dalitzpm_dat_hist("dalitzpm_dat_hist",
1911  "",
1912  m12->getNumBins(),
1913  m12->getLowerLimit(),
1914  m12->getUpperLimit(),
1915  m13->getNumBins(),
1916  m13->getLowerLimit(),
1917  m13->getUpperLimit());
1918  dalitzpm_dat_hist.SetStats(false);
1919  dalitzpm_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{+} #pi^{0}) [GeV]");
1920  dalitzpm_dat_hist.GetYaxis()->SetTitle("m^{2}(#pi^{-} #pi^{0}) [GeV]");
1921  TH2F dalitzpm_pdf_hist("dalitzpm_pdf_hist",
1922  "",
1923  m12->getNumBins(),
1924  m12->getLowerLimit(),
1925  m12->getUpperLimit(),
1926  m13->getNumBins(),
1927  m13->getLowerLimit(),
1928  m13->getUpperLimit());
1929  dalitzpm_pdf_hist.SetStats(false);
1930 
1931  TH2F dalitzp0_dat_hist("dalitzp0_dat_hist",
1932  "",
1933  m12->getNumBins(),
1934  m12->getLowerLimit(),
1935  m12->getUpperLimit(),
1936  m13->getNumBins(),
1937  m13->getLowerLimit(),
1938  m13->getUpperLimit());
1939  dalitzp0_dat_hist.SetStats(false);
1940  dalitzp0_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{+} #pi^{0}) [GeV]");
1941  dalitzp0_dat_hist.GetYaxis()->SetTitle("m^{2}(#pi^{-} #pi^{+}) [GeV]");
1942  TH2F dalitzp0_pdf_hist("dalitzp0_pdf_hist",
1943  "",
1944  m12->getNumBins(),
1945  m12->getLowerLimit(),
1946  m12->getUpperLimit(),
1947  m13->getNumBins(),
1948  m13->getLowerLimit(),
1949  m13->getUpperLimit());
1950  dalitzp0_pdf_hist.SetStats(false);
1951 
1952  TH2F dalitzm0_dat_hist("dalitzm0_dat_hist",
1953  "",
1954  m12->getNumBins(),
1955  m12->getLowerLimit(),
1956  m12->getUpperLimit(),
1957  m13->getNumBins(),
1958  m13->getLowerLimit(),
1959  m13->getUpperLimit());
1960  dalitzm0_dat_hist.SetStats(false);
1961  dalitzm0_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{-} #pi^{0}) [GeV]");
1962  dalitzm0_dat_hist.GetYaxis()->SetTitle("m^{2}(#pi^{+} #pi^{-}) [GeV]");
1963  TH2F dalitzm0_pdf_hist("dalitzm0_pdf_hist",
1964  "",
1965  m12->getNumBins(),
1966  m12->getLowerLimit(),
1967  m12->getUpperLimit(),
1968  m13->getNumBins(),
1969  m13->getLowerLimit(),
1970  m13->getUpperLimit());
1971  dalitzm0_pdf_hist.SetStats(false);
1972 
1973  TH1F *bkg3_pdfs[6];
1974  TH1F *bkg3_data[6];
1975  double num_sigma_dat[6];
1976  double num_sigma_pdf[6];
1977 
1978  for(int i = 0; i < 6; ++i) {
1979  sprintf(strbuffer, "bkg_sigma_pdf_%i", i);
1980  bkg3_pdfs[i] = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1981  sprintf(strbuffer, "bkg_sigma_dat_%i", i);
1982  bkg3_data[i] = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
1983 
1984  num_sigma_dat[i] = 0;
1985  num_sigma_pdf[i] = 0;
1986 
1987  bkg3_data[i]->SetStats(false);
1988  bkg3_data[i]->SetMarkerStyle(8);
1989  bkg3_data[i]->SetMarkerSize(1.2);
1990  bkg3_data[i]->GetXaxis()->SetTitle("Decay time error [ps]");
1991  bkg3_data[i]->GetYaxis()->SetTitle("Events / 8 fs");
1992 
1993  bkg3_pdfs[i]->SetStats(false);
1994  bkg3_pdfs[i]->SetLineColor(kBlue);
1995  bkg3_pdfs[i]->SetLineWidth(3);
1996  }
1997 
1998  double totalPdf = 0;
1999  double totalDat = 0;
2000 
2001  for(unsigned int evt = 0; evt < data->getNumEvents(); ++evt) {
2002  double currTime = data->getValue(*dtime, evt);
2003  dtime_dat_hist.Fill(currTime);
2004 
2005  double currSigma = data->getValue(*sigma, evt);
2006  sigma_dat_hist.Fill(currSigma);
2007 
2008  double currm12 = data->getValue(*m12, evt);
2009  m12_dat_hist.Fill(currm12);
2010 
2011  double currm13 = data->getValue(*m13, evt);
2012  m13_dat_hist.Fill(currm13);
2013 
2014  dalitzpm_dat_hist.Fill(currm12, currm13);
2015 
2016  double currm23 = cpuGetM23(currm12, currm13);
2017  m23_dat_hist.Fill(currm23);
2018 
2019  dalitzp0_dat_hist.Fill(currm12, currm23);
2020  dalitzm0_dat_hist.Fill(currm13, currm23);
2021 
2022  totalDat++;
2023 
2024  int m23bin = (int)floor(currm23 / 0.5);
2025  bkg3_data[m23bin]->Fill(currSigma);
2026  num_sigma_dat[m23bin]++;
2027  }
2028 
2029  double maxBinContent = 0;
2030  int bestI = 0;
2031  int bestJ = 0;
2032 
2033  for(int i = 1; i <= m12->getNumBins(); ++i) {
2034  for(int j = 1; j <= m13->getNumBins(); ++j) {
2035  double curr = dalitzpm_dat_hist.GetBinContent(i, j);
2036 
2037  if(curr < maxBinContent)
2038  continue;
2039 
2040  maxBinContent = curr;
2041  bestI = i;
2042  bestJ = j;
2043  }
2044  }
2045 
2046  std::cout << "Max bin content: " << maxBinContent << " (" << bestI << ", " << bestJ << ")\n";
2047 
2048  bool dependsOnSigma = true;
2049  std::vector<Observable> obses = overallSignal->getObservables();
2050 
2051  if(std::find(obses.begin(), obses.end(), *sigma) == obses.end())
2052  dependsOnSigma = false;
2053 
2054  // overallSignal->setDebugMask(1);
2055 
2056  wBkg1->setValue(0);
2057  const int division = 2;
2058 
2059  for(int half = 0; half < division; ++half) {
2060  std::vector<Observable> vars;
2061  vars.push_back(*m12);
2062  vars.push_back(*m13);
2063  vars.push_back(*dtime);
2064  vars.push_back(*sigma);
2065  vars.push_back(*eventNumber);
2066  vars.push_back(*wBkg1);
2067 
2068  UnbinnedDataSet currData(vars);
2069 
2070  sigma->setValue(0.5);
2071  int evtCounter = 0;
2072 
2073  for(int i = 0; i < m12->getNumBins(); ++i) {
2074  m12->setValue(m12->getLowerLimit()
2075  + (m12->getUpperLimit() - m12->getLowerLimit()) * (i + 0.5) / m12->getNumBins());
2076 
2077  for(int j = 0; j < m13->getNumBins(); ++j) {
2078  m13->setValue(m13->getLowerLimit()
2079  + (m13->getUpperLimit() - m13->getLowerLimit()) * (j + 0.5) / m13->getNumBins());
2080 
2082  continue;
2083 
2084  for(int l = half; l < dtime->getNumBins(); l += division) {
2085  dtime->setValue(dtime->getLowerLimit()
2086  + (dtime->getUpperLimit() - dtime->getLowerLimit()) * (l + 0.5)
2087  / dtime->getNumBins());
2088  eventNumber->setValue(evtCounter);
2089  evtCounter++;
2090  currData.addEvent();
2091  }
2092  }
2093  }
2094 
2095  for(int k = 0; k < sigma->getNumBins(); ++k) {
2096  sigma->setValue(sigma->getLowerLimit()
2097  + (sigma->getUpperLimit() - sigma->getLowerLimit()) * (k + 0.5) / sigma->getNumBins());
2098 
2099  if(0 == k % 10)
2100  std::cout << "sigma iteration " << half << " " << k << std::endl;
2101 
2102  currData.setValueForAllEvents(*sigma);
2103  overallSignal->setData(&currData);
2104 
2105  if(0 == k) {
2106  if(signalDalitz)
2107  signalDalitz->setDataSize(currData.getNumEvents(), 6);
2108 
2109  if(incsum1)
2110  incsum1->setDataSize(currData.getNumEvents(), 6);
2111 
2112  if(incsum2)
2113  incsum2->setDataSize(currData.getNumEvents(), 6);
2114 
2115  if(incsum3)
2116  incsum3->setDataSize(currData.getNumEvents(), 6);
2117 
2118  if(incsum4)
2119  incsum4->setDataSize(currData.getNumEvents(), 6);
2120 
2121  if(incsum5)
2122  incsum5->setDataSize(currData.getNumEvents(), 6);
2123 
2124  if(incsum6)
2125  incsum6->setDataSize(currData.getNumEvents(), 6);
2126  }
2127 
2128  std::vector<std::vector<double>> pdfValues = overallSignal->getCompProbsAtDataPoints();
2129 
2130  for(unsigned int j = 0; j < pdfValues[0].size(); ++j) {
2131  double currTime = currData.getValue(*dtime, j);
2132  dtime_pdf_hist.Fill(currTime, pdfValues[0][j]);
2133 
2134  double currSigma = currData.getValue(*sigma, j);
2135  sigma_pdf_hist.Fill(currSigma, pdfValues[0][j]);
2136 
2137  // Um... these two are switched? Weirdness...
2138  double currm12 = currData.getValue(*m13, j);
2139  m12_pdf_hist.Fill(currm12, pdfValues[0][j]);
2140 
2141  double currm13 = currData.getValue(*m12, j);
2142  m13_pdf_hist.Fill(currm13, pdfValues[0][j]);
2143  dalitzpm_pdf_hist.Fill(currm12, currm13, pdfValues[0][j]);
2144 
2145  double currm23 = cpuGetM23(currm12, currm13);
2146  m23_pdf_hist.Fill(currm23, pdfValues[0][j]);
2147  dalitzp0_pdf_hist.Fill(currm12, currm23, pdfValues[0][j]);
2148  dalitzm0_pdf_hist.Fill(currm13, currm23, pdfValues[0][j]);
2149 
2150  int currM23Bin = (int)(currm23 / 0.5);
2151  bkg3_pdfs[currM23Bin]->Fill(currSigma, pdfValues[0][j]);
2152  num_sigma_pdf[currM23Bin] += pdfValues[0][j];
2153 
2154  totalPdf += pdfValues[0][j];
2155 
2156  if(std::isnan(pdfValues[0][j])) {
2157  std::cout << "Major problem: " << k << " " << j << std::endl;
2158  assert(false);
2159  }
2160 
2161  if(std::isinf(pdfValues[0][j])) {
2162  std::cout << "Infinity " << k << " " << j << std::endl;
2163  assert(false);
2164  }
2165  }
2166 
2167  // If PDF doesn't depend on sigma, don't project from that dimension.
2168  if(!dependsOnSigma)
2169  break;
2170  }
2171  }
2172 
2173  // std::cout << "Final values: " << totalDat << " " << totalPdf << std::endl;
2174 
2175  for(int i = 1; i <= dtime->getNumBins(); ++i) {
2176  dtime_pdf_hist.SetBinContent(i, dtime_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
2177  }
2178 
2179  for(int i = 1; i <= sigma->getNumBins(); ++i) {
2180  sigma_pdf_hist.SetBinContent(i, sigma_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
2181 
2182  for(int j = 0; j < 6; ++j) {
2183  bkg3_pdfs[j]->SetBinContent(i, bkg3_pdfs[j]->GetBinContent(i) * num_sigma_dat[j] / num_sigma_pdf[j]);
2184  }
2185  }
2186 
2187  for(int i = 1; i <= m12->getNumBins(); ++i) {
2188  m12_pdf_hist.SetBinContent(i, m12_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
2189  }
2190 
2191  for(int i = 1; i <= m13->getNumBins(); ++i) {
2192  m13_pdf_hist.SetBinContent(i, m13_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
2193  m23_pdf_hist.SetBinContent(i, m23_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
2194  }
2195 
2196  for(int i = 1; i <= m12->getNumBins(); ++i) {
2197  for(int j = 1; j <= m13->getNumBins(); ++j) {
2198  dalitzpm_pdf_hist.SetBinContent(i, j, dalitzpm_pdf_hist.GetBinContent(i, j) * totalDat / totalPdf);
2199  dalitzp0_pdf_hist.SetBinContent(i, j, dalitzp0_pdf_hist.GetBinContent(i, j) * totalDat / totalPdf);
2200  dalitzm0_pdf_hist.SetBinContent(i, j, dalitzm0_pdf_hist.GetBinContent(i, j) * totalDat / totalPdf);
2201  }
2202  }
2203 
2204  ChisqInfo *chisq = getAdaptiveChisquare(&dalitzpm_dat_hist, &dalitzpm_pdf_hist);
2205 
2206  std::cout << "Chisquare: " << chisq->chisq << " / " << chisq->dof << std::endl;
2207  foodal->cd();
2208  foodal->SetLogz(false);
2209  chisq->contribPlot->SetStats(false);
2210  chisq->contribPlot->Draw("colz");
2211  foodal->SaveAs((plotdir + "/chisq.png").c_str());
2212  foo->cd();
2213 
2214  coarseBin(dalitzpm_pdf_hist, 2);
2215  coarseBin(dalitzpm_dat_hist, 2);
2216  coarseBin(dalitzp0_pdf_hist, 2);
2217  coarseBin(dalitzp0_dat_hist, 2);
2218  coarseBin(dalitzm0_pdf_hist, 2);
2219  coarseBin(dalitzm0_dat_hist, 2);
2220 
2221  dtime_dat_hist.Draw("p");
2222  dtime_pdf_hist.Draw("lsame");
2223  foo->SaveAs((plotdir + "/dtime_fit.png").c_str());
2224  foo->SetLogy(true);
2225  foo->SaveAs((plotdir + "/dtime_fit_log.png").c_str());
2226  foo->SetLogy(false);
2227 
2228  for(int i = 0; i < 6; ++i) {
2229  if(!dependsOnSigma) {
2230  bkg3_data[i]->Draw("p");
2231  } else if(bkg3_data[i]->GetMaximum() > bkg3_pdfs[i]->GetMaximum()) {
2232  bkg3_data[i]->Draw("p");
2233  bkg3_pdfs[i]->Draw("lsame");
2234  } else {
2235  bkg3_pdfs[i]->Draw("l");
2236  bkg3_data[i]->Draw("psame");
2237  }
2238 
2239  sprintf(strbuffer, "%i", i);
2240  TText slicenum;
2241  slicenum.DrawTextNDC(0.2, 0.8, strbuffer);
2242 
2243  foo->SaveAs((plotdir + bkg3_pdfs[i]->GetName() + ".png").c_str());
2244  foo->SetLogy(true);
2245  foo->SaveAs((plotdir + bkg3_pdfs[i]->GetName() + "_log.png").c_str());
2246  foo->SetLogy(false);
2247  }
2248 
2249  m13_dat_hist.Draw("p");
2250  m13_pdf_hist.Draw("lsame");
2251  foo->SaveAs((plotdir + "/m13_fit.png").c_str());
2252  foo->SetLogy(true);
2253  foo->SaveAs((plotdir + "/m13_fit_log.png").c_str());
2254  foo->SetLogy(false);
2255 
2256  sigma_dat_hist.Draw("p");
2257  sigma_pdf_hist.Draw("lsame");
2258  foo->SaveAs((plotdir + "/sigma_fit.png").c_str());
2259  foo->SetLogy(true);
2260  foo->SaveAs((plotdir + "/sigma_fit_log.png").c_str());
2261  foo->SetLogy(false);
2262 
2263  m12_dat_hist.Draw("p");
2264  m12_pdf_hist.Draw("lsame");
2265  foo->SaveAs((plotdir + "/m12_fit.png").c_str());
2266  foo->SetLogy(true);
2267  foo->SaveAs((plotdir + "/m12_fit_log.png").c_str());
2268  foo->SetLogy(false);
2269 
2270  m23_dat_hist.Draw("p");
2271  m23_pdf_hist.Draw("lsame");
2272  foo->SaveAs((plotdir + "/m23_fit.png").c_str());
2273  foo->SetLogy(true);
2274  foo->SaveAs((plotdir + "/m23_fit_log.png").c_str());
2275  foo->SetLogy(false);
2276 
2277  foodal->cd();
2278  dalitzpm_dat_hist.Draw("colz");
2279 
2280  for(int slice = 0; slice < 6; ++slice) {
2281  double line_m12 = cpuGetM23(0, 0.5 * (slice + 1));
2282  TLine sliceLine;
2283  sliceLine.SetLineWidth(2);
2284  sliceLine.DrawLine(0, line_m12, line_m12, 0);
2285  sprintf(strbuffer, "%i", slice);
2286  TText sliceNum;
2287  sliceNum.DrawText(0.25, line_m12 - 0.25, strbuffer);
2288  }
2289 
2290  foodal->SaveAs((plotdir + "/dalitzpm_dat.png").c_str());
2291  foodal->SetLogz(true);
2292  foodal->SaveAs((plotdir + "/dalitzpm_dat_log.png").c_str());
2293  foodal->SetLogz(false);
2294  dalitzp0_dat_hist.Draw("colz");
2295  foodal->SaveAs((plotdir + "/dalitzp0_dat.png").c_str());
2296  foodal->SetLogz(true);
2297  foodal->SaveAs((plotdir + "/dalitzp0_dat_log.png").c_str());
2298  foodal->SetLogz(false);
2299  dalitzm0_dat_hist.Draw("colz");
2300  foodal->SaveAs((plotdir + "/dalitzm0_dat.png").c_str());
2301  foodal->SetLogz(true);
2302  foodal->SaveAs((plotdir + "/dalitzm0_dat_log.png").c_str());
2303  foodal->SetLogz(false);
2304 
2305  dalitzpm_pdf_hist.Draw("colz");
2306  foodal->SaveAs((plotdir + "/dalitzpm_pdf.png").c_str());
2307  foodal->SetLogz(true);
2308  foodal->SaveAs((plotdir + "/dalitzpm_pdf_log.png").c_str());
2309  foodal->SetLogz(false);
2310  dalitzpm_pdf_hist.GetZaxis()->SetRangeUser(0, dalitzpm_dat_hist.GetMaximum());
2311  dalitzpm_pdf_hist.Draw("colz");
2312  foodal->SaveAs((plotdir + "/dalitzpm_pdf_matched.png").c_str());
2313  foodal->SetLogz(true);
2314  foodal->SaveAs((plotdir + "/dalitzpm_pdf_matched_log.png").c_str());
2315  foodal->SetLogz(false);
2316  dalitzp0_pdf_hist.GetZaxis()->SetRangeUser(0, dalitzp0_dat_hist.GetMaximum());
2317  dalitzp0_pdf_hist.Draw("colz");
2318  foodal->SaveAs((plotdir + "/dalitzp0_pdf.png").c_str());
2319  foodal->SetLogz(true);
2320  foodal->SaveAs((plotdir + "/dalitzp0_pdf_log.png").c_str());
2321  foodal->SetLogz(false);
2322  dalitzm0_pdf_hist.GetZaxis()->SetRangeUser(0, dalitzm0_dat_hist.GetMaximum());
2323  dalitzm0_pdf_hist.Draw("colz");
2324  foodal->SaveAs((plotdir + "/dalitzm0_pdf.png").c_str());
2325  foodal->SetLogz(true);
2326  foodal->SaveAs((plotdir + "/dalitzm0_pdf_log.png").c_str());
2327  foodal->SetLogz(false);
2328 
2329  TH1F pull_pm_hist("pull_pm_hist", "", 100, -5, 5);
2330  pull_pm_hist.GetXaxis()->SetTitle("(Data - PDF) / sqrt(Data)");
2331  pull_pm_hist.GetYaxis()->SetTitle("Bins / 0.1");
2332 
2333  for(int i = 1; i <= m12->getNumBins(); ++i) {
2334  double m12loedge
2335  = m12->getLowerLimit() + ((m12->getUpperLimit() - m12->getLowerLimit()) / m12->getNumBins()) * (i - 1);
2336  double m12hiedge
2337  = m12->getLowerLimit() + ((m12->getUpperLimit() - m12->getLowerLimit()) / m12->getNumBins()) * (i);
2338 
2339  for(int j = 1; j <= m13->getNumBins(); ++j) {
2340  double m13loedge
2341  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j - 1);
2342 
2343  if(!cpuDalitz(m12loedge, m13loedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
2344  dalitzpm_dat_hist.SetBinContent(i, j, 0);
2345  continue;
2346  }
2347 
2348  if(!cpuDalitz(m12hiedge, m13loedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
2349  dalitzpm_dat_hist.SetBinContent(i, j, 0);
2350  continue;
2351  }
2352 
2353  double m13hiedge
2354  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j);
2355 
2356  if(!cpuDalitz(m12loedge, m13hiedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
2357  dalitzpm_dat_hist.SetBinContent(i, j, 0);
2358  continue;
2359  }
2360 
2361  if(!cpuDalitz(m12hiedge, m13hiedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
2362  dalitzpm_dat_hist.SetBinContent(i, j, 0);
2363  continue;
2364  }
2365 
2366  double dat = dalitzpm_dat_hist.GetBinContent(i, j);
2367  double pdf = dalitzpm_pdf_hist.GetBinContent(i, j);
2368 
2369  double pullval = (dat - pdf) / sqrt(max(1.0, dat));
2370  dalitzpm_dat_hist.SetBinContent(i, j, pullval);
2371  pull_pm_hist.Fill(pullval);
2372  }
2373 
2374  for(int j = 1; j <= m13->getNumBins(); ++j) {
2375  double m23loedge
2376  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j - 1);
2377 
2378  // To get 12, 23 instead of 12, 13, just exchange 1 and 2.
2379  if(!cpuDalitz(m12loedge, m23loedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2380  dalitzp0_dat_hist.SetBinContent(i, j, 0);
2381  continue;
2382  }
2383 
2384  if(!cpuDalitz(m12hiedge, m23loedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2385  dalitzp0_dat_hist.SetBinContent(i, j, 0);
2386  continue;
2387  }
2388 
2389  double m23hiedge
2390  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j);
2391 
2392  if(!cpuDalitz(m12loedge, m23hiedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2393  dalitzp0_dat_hist.SetBinContent(i, j, 0);
2394  continue;
2395  }
2396 
2397  if(!cpuDalitz(m12hiedge, m23hiedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2398  dalitzp0_dat_hist.SetBinContent(i, j, 0);
2399  continue;
2400  }
2401 
2402  double dat = dalitzp0_dat_hist.GetBinContent(i, j);
2403  double pdf = dalitzp0_pdf_hist.GetBinContent(i, j);
2404 
2405  dalitzp0_dat_hist.SetBinContent(i, j, (dat - pdf) / sqrt(max(1.0, dat)));
2406  }
2407 
2408  // NB, this exploits symmetry 12 and 13 by treating the outer loop as 13.
2409  for(int j = 1; j <= m13->getNumBins(); ++j) {
2410  double m23loedge
2411  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j - 1);
2412 
2413  if(!cpuDalitz(m12loedge, m23loedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2414  dalitzm0_dat_hist.SetBinContent(i, j, 0);
2415  continue;
2416  }
2417 
2418  if(!cpuDalitz(m12hiedge, m23loedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2419  dalitzm0_dat_hist.SetBinContent(i, j, 0);
2420  continue;
2421  }
2422 
2423  double m23hiedge
2424  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j);
2425 
2426  if(!cpuDalitz(m12loedge, m23hiedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2427  dalitzm0_dat_hist.SetBinContent(i, j, 0);
2428  continue;
2429  }
2430 
2431  if(!cpuDalitz(m12hiedge, m23hiedge, _mD0, piPlusMass, piZeroMass, piPlusMass)) {
2432  dalitzm0_dat_hist.SetBinContent(i, j, 0);
2433  continue;
2434  }
2435 
2436  double dat = dalitzm0_dat_hist.GetBinContent(i, j);
2437  double pdf = dalitzm0_pdf_hist.GetBinContent(i, j);
2438 
2439  dalitzm0_dat_hist.SetBinContent(i, j, (dat - pdf) / sqrt(max(1.0, dat)));
2440  }
2441  }
2442 
2443  dalitzpm_dat_hist.GetZaxis()->SetRangeUser(-5, 5);
2444  dalitzpm_dat_hist.Draw("colz");
2445  foodal->SaveAs((plotdir + "/dalitzpm_pull.png").c_str());
2446  dalitzp0_dat_hist.GetZaxis()->SetRangeUser(-5, 5);
2447  dalitzp0_dat_hist.Draw("colz");
2448  foodal->SaveAs((plotdir + "/dalitzp0_pull.png").c_str());
2449  dalitzm0_dat_hist.GetZaxis()->SetRangeUser(-5, 5);
2450  dalitzm0_dat_hist.Draw("colz");
2451  foodal->SaveAs((plotdir + "/dalitzm0_pull.png").c_str());
2452 
2453  foo->cd();
2454 
2455  pull_pm_hist.Draw("");
2456  foo->SaveAs((plotdir + "/pull_pm_hist.png").c_str());
2457 }
2458 
2460  // This is complicated. I want to make a function such that the parameters
2461  // depend on position in the Dalitz plot. I want to use stripes of m23, because
2462  // that seems to be the real shape of the dependence, and I haven't got so much
2463  // data for bkg3; so stripes work better than boxes of m12, m13. So I need a
2464  // transform from m12, m13 to m23. Then I need a transform from m23 to a bin number.
2465  // Finally I need a transform from bin number to function. Keep the tongue straight
2466  // in the mouth, now!
2467 
2468  vector<Observable> obses;
2469  vector<Variable> offsets;
2470  vector<Variable> coefficients;
2471  vector<PdfBase *> components;
2472  vector<double> limits;
2473  vector<double> binSizes;
2474  vector<int> numBins;
2475 
2476  // First the transform to m23.
2477  // m23 = _mD02 + piZeroMass*piZeroMass + piPlusMass*piPlusMass + piPlusMass*piPlusMass - massPZ - massPM
2478 
2479  obses.push_back(*m12);
2480  obses.push_back(*m13);
2481  coefficients.push_back(constantBigM);
2482  coefficients.push_back(constantMinusOne);
2483  coefficients.push_back(constantMinusOne);
2484  offsets.push_back(constantZero);
2485  offsets.push_back(constantZero);
2486  PolynomialPdf *m23_transform = new PolynomialPdf("m23_transform", obses, coefficients, offsets, 1);
2487 
2488  // Now create the BinTransform which will make bins out of m23 values.
2489  obses.clear();
2490  coefficients.clear();
2491  obses.push_back(*m12); // Fake dependence; CompositePdf will create a fake event using this index.
2492  limits.push_back(0); // Bins of m23 start at 0.
2493  binSizes.push_back(3.0 / m23Slices);
2494  numBins.push_back(m23Slices);
2495  BinTransformPdf *m23_binMap = new BinTransformPdf("m23_binMap", obses, limits, binSizes, numBins);
2496 
2497  // Now make a composite, so that the BinTransform takes the Polynomial result as input.
2498  CompositePdf *m23_composite = new CompositePdf("m23_composite", m23_transform, m23_binMap);
2499  return m23_composite;
2500 }
2501 
2503  GooPdf *m23_composite = make_m23_transform();
2504  std::vector<std::unique_ptr<BinnedDataSet>> sigmaHists;
2505 
2506  for(int i = 0; i < m23Slices; ++i)
2507  sigmaHists.emplace_back(new BinnedDataSet(*sigma));
2508 
2509  std::ifstream reader;
2510  std::string fname = app_ptr->get_filename("./dataFiles/signalMC_truth_mm_0.txt", "examples/pipipi0DPFit");
2511  readWrapper(reader, fname);
2512  std::string buffer;
2513 
2514  while(!reader.eof()) {
2515  reader >> buffer;
2516 
2517  if(buffer == "====")
2518  break;
2519 
2520  std::cout << buffer;
2521  }
2522 
2523  double dummy = 0;
2524  double m23 = 0;
2525 
2526  while(!reader.eof()) {
2527  reader >> dummy;
2528  reader >> m23;
2529  reader >> dummy;
2530  reader >> dummy;
2531  reader >> dummy;
2532  reader >> dummy;
2533  reader >> dummy;
2534  reader >> *dtime;
2535  reader >> *sigma;
2536 
2537  reader >> dummy;
2538  reader >> dummy;
2539  reader >> dummy;
2540  reader >> dummy;
2541  reader >> dummy;
2542  reader >> dummy;
2543  reader >> dummy;
2544  reader >> dummy;
2545  reader >> dummy;
2546  reader >> dummy;
2547  reader >> dummy;
2548 
2549  if(dtime->getValue() < dtime->getLowerLimit())
2550  continue;
2551 
2552  if(dtime->getValue() > dtime->getUpperLimit())
2553  continue;
2554 
2555  int bin = (int)floor((m23 / 3.0) * m23Slices);
2556  sigmaHists[bin]->addEvent();
2557  }
2558 
2559  vector<GooPdf *> jsuList;
2560 
2561  for(int i = 0; i < m23Slices; ++i) {
2562  sprintf(strbuffer, "signal_sigma_hist_%i", i);
2563  SmoothHistogramPdf *hist = new SmoothHistogramPdf(strbuffer, sigmaHists[i].get(), constantZero);
2564  jsuList.push_back(hist);
2565  }
2566 
2567  return new MappedPdf("signalSigmaHist", m23_composite, jsuList);
2568 }
2569 
2571  GooPdf *m23_composite = make_m23_transform();
2572 
2573  vector<GooPdf *> jsuList;
2574  vector<ConvolutionPdf *> convList;
2575  bool useShare = false;
2576 
2577  for(int i = 0; i < m23Slices; ++i) {
2578  sprintf(strbuffer, "bkg%i_sigma_slice%i_expalpha", bkgnum, i);
2579  Variable exp_alpha(strbuffer, 7.50, 0.10, 0, 10.00);
2580  sprintf(strbuffer, "bkg%i_sigma_slice%i_gauss_meana", bkgnum, i);
2581  Variable g_meana(strbuffer, 0.20, 0.01, 0.00, 0.80);
2582  sprintf(strbuffer, "bkg%i_sigma_slice%i_gauss_sigma", bkgnum, i);
2583  Variable g_sigma(strbuffer, 0.03, 0.01, 0.01, 0.20);
2584 
2585  sprintf(strbuffer, "bkg%i_sigma_slice%i_conv", bkgnum, i);
2586  ExpGausPdf *expfunc = new ExpGausPdf(strbuffer, *sigma, g_meana, g_sigma, exp_alpha);
2587  jsuList.push_back(expfunc);
2588  }
2589 
2590  if(useShare) {
2591  for(vector<ConvolutionPdf *>::iterator conv = convList.begin(); conv != convList.end(); ++conv) {
2592  (*conv)->registerOthers(convList);
2593  }
2594  }
2595 
2596  sprintf(strbuffer, "bkg%i_sigma_map", bkgnum);
2597  return new MappedPdf(strbuffer, m23_composite, jsuList);
2598 }
2599 
2601  weightHistogram = new TH2F("weightHistogram",
2602  "",
2603  m12->getNumBins(),
2604  m12->getLowerLimit(),
2605  m12->getUpperLimit(),
2606  m13->getNumBins(),
2607  m13->getLowerLimit(),
2608  m13->getUpperLimit());
2609  weightHistogram->SetStats(false);
2610  double step12 = (m12->getUpperLimit() - m12->getLowerLimit()) / m12->getNumBins();
2611  double step13 = (m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins();
2612 
2613  for(int i = 1; i <= m12->getNumBins(); ++i) {
2614  for(int j = 1; j < m13->getNumBins(); ++j) {
2615  double maxCount = 0;
2616  double count = 0;
2617 
2618  for(double currM12 = m12->getLowerLimit() + step12 * (i - 1) + 0.05 * step12;
2619  currM12 < m12->getLowerLimit() + step12 * i;
2620  currM12 += 0.1 * step12) {
2621  for(double currM13 = m13->getLowerLimit() + step13 * (j - 1) + 0.05 * step13;
2622  currM13 < m13->getLowerLimit() + step13 * j;
2623  currM13 += 0.1 * step13) {
2624  maxCount++;
2625 
2626  if(!cpuDalitz(currM12, currM13, _mD0, piZeroMass, piPlusMass, piPlusMass))
2627  continue;
2628 
2629  count++;
2630  }
2631  }
2632 
2633  if(0.1 > maxCount)
2634  continue;
2635 
2636  count /= maxCount;
2637  weightHistogram->SetBinContent(i, j, count);
2638  }
2639  }
2640 
2641  underlyingBins = new TH2F("underlyingBins",
2642  "",
2643  m12->getNumBins(),
2644  m12->getLowerLimit(),
2645  m12->getUpperLimit(),
2646  m13->getNumBins(),
2647  m13->getLowerLimit(),
2648  m13->getUpperLimit());
2649  underlyingBins->SetStats(false);
2650 }
2651 
2653  makeKzeroVeto();
2654 
2655  int oldBins1 = m12->getNumBins();
2656  int oldBins2 = m13->getNumBins();
2657  // Too fine a binning here leads to bad results due to fluctuations.
2658  m12->setNumBins(120);
2659  m13->setNumBins(120);
2660  vector<Observable> lvars;
2661  lvars.push_back(*m12);
2662  lvars.push_back(*m13);
2663  binEffData = new BinnedDataSet(lvars);
2665  std::cout << "Loading efficiency data" << std::endl;
2666  std::string fname = app_ptr->get_filename("./dataFiles/efficiency_flat.txt", "examples/pipipi0DPFit");
2667  loadDataFile(fname, &effdata);
2668 
2669  if(saveEffPlot) {
2670  system("mkdir plots_from_mixfit");
2671  foodal->cd();
2672  underlyingBins->Draw("colz");
2673  foodal->SaveAs("plots_from_mixfit/efficiency_bins.png");
2674  foodal->SetLogz(true);
2675  foodal->SaveAs("plots_from_mixfit/efficiency_bins_log.png");
2676  foo->cd();
2677  }
2678 
2679  GooPdf *eff = 0;
2680 
2681  // Polynomial version
2682  if(polyEff)
2683  eff = makeEfficiencyPoly();
2684  // SmoothHistogram version
2685  else
2686  eff = makeEfficiencyPdf();
2687 
2688  eff->setData(effdata);
2689 
2690  if(minuit1) {
2691  GooFit::FitManagerMinuit1 effpdf(eff);
2692  effpdf.fit();
2693  } else {
2694  GooFit::FitManagerMinuit2 effpdf(eff);
2695  effpdf.fit();
2696  }
2697 
2698  eff->setParameterConstantness(true);
2699  delete binEffData;
2700  binEffData = nullptr;
2701  delete effdata;
2702  effdata = nullptr;
2703 
2704  m12->setNumBins(oldBins1);
2705  m13->setNumBins(oldBins2);
2706 
2707  comps.clear();
2708  comps.push_back(eff);
2709  comps.push_back(kzero_veto);
2710  ProdPdf *effWithVeto = new ProdPdf("effWithVeto", comps);
2711 
2712  std::cout << "Creating signal PDF\n";
2713  signalDalitz = makeSignalPdf(0, effWithVeto);
2714 
2715  std::cout << "Creating sigma PDF\n";
2716 
2717  // sig0_jsugg = makeSigmaMap();
2718  // sig0_jsugg = make1BinSigmaMap();
2719  // sig0_jsugg = make4BinSigmaMap();
2720  // sig0_jsugg = makeMikhailJSU_gg();
2721  if(useHistogramSigma)
2722  sig0_jsugg = makeSigmaHists();
2723  else
2724  sig0_jsugg = makeBkg_sigma_strips(0);
2725 
2726  sig0_jsugg->addSpecialMask(PdfBase::ForceSeparateNorm);
2727  // sig0_jsugg = makeSignalJSU_gg(-1, false);
2728 
2729  /*
2730  sig0_jsugg->setData(data);
2731  FitManager jsupdf(sig0_jsugg);
2732  gettimeofday(&startTime, nullptr);
2733  jsupdf.fit();
2734  gettimeofday(&stopTime, nullptr);
2735  timersub(&stopTime, &startTime, &totalTime);
2736  std::cout << "Time for sigma fit : " << totalTime.tv_sec + totalTime.tv_usec/1000000.0 << " seconds." << std::endl;
2737  */
2738  sprintf(strbuffer, "signal_sigma_%islices_pdf.txt", m23Slices);
2739  fname = app_ptr->get_filename(strbuffer, "examples/pipipi0DPFit");
2740  readFromFile(sig0_jsugg, strbuffer);
2741  sig0_jsugg->setParameterConstantness(true);
2742 
2743  comps.clear();
2744  comps.push_back(signalDalitz);
2745  comps.push_back(sig0_jsugg);
2746  std::cout << "Creating overall PDF\n";
2747  ProdPdf *overallSignal = new ProdPdf("overallSignal", comps);
2748  return overallSignal;
2749 }
2750 
2751 int runTruthMCFit(std::string fname, bool noPlots = true) {
2753 
2754  std::cout << "Loading MC data from " << fname << std::endl;
2755  loadDataFile(fname);
2756  GooPdf *overallSignal = makeOverallSignal();
2757 
2758  signalDalitz->setDataSize(data->getNumEvents()); // Default 5 is ok here, no event weighting
2759  overallSignal->setData(data);
2760  // overallSignal->setDebugMask(1);
2761 
2762  int retval;
2763  if(minuit1) {
2764  GooFit::FitManagerMinuit1 datapdf(overallSignal);
2765  datapdf.fit();
2766  retval = datapdf;
2767  } else {
2768  GooFit::FitManagerMinuit2 datapdf(overallSignal);
2769  datapdf.fit();
2770  retval = datapdf;
2771  }
2772 
2773  // overallSignal->setDebugMask(0);
2774 
2775  fmt::print("Fit results Toy fit TruthMC fit:\n"
2776  "tau : {:.3}\n"
2777  "xmixing: ({:.3})\%\n"
2778  "ymixing: ({:.3})\%\n",
2780  100 * Uncertain(dtop0pp._xmixing),
2781  100 * Uncertain(dtop0pp._ymixing));
2782 
2783  if(noPlots)
2784  return retval;
2785 
2786  makeDalitzPlots(overallSignal, "./plots_from_mixfit/fullMCfit/");
2787  return retval;
2788 }
2789 
2790 int runGeneratedMCFit(std::string fname, int genResolutions, double dplotres) {
2792  std::cout << "Loading (generated) MC data from " << fname << std::endl;
2793  dtime->setUpperLimit(6);
2794  dtime->setLowerLimit(0);
2795  sigma->setUpperLimit(1.0);
2796  loadDataFile(fname);
2797 
2798  TRandom donram(42);
2799  std::vector<Observable> vars = data->getObservables();
2800  UnbinnedDataSet *smearedData = new UnbinnedDataSet(vars);
2801 
2802  if(0 != genResolutions) {
2803  int numEvents = data->getNumEvents();
2804 
2805  for(int i = 0; i < numEvents; ++i) {
2806  data->loadEvent(i);
2807 
2808  double smear1 = 0;
2809  double smear2 = 0;
2810 
2811  if(DplotRes & genResolutions) {
2812  if(100 > smearedData->getNumEvents())
2813  std::cout << "Before smear: " << m12->getValue() << " " << m13->getValue();
2814 
2815  smear1 = donram.Gaus(0, dplotres);
2816  smear2 = donram.Gaus(0, dplotres);
2817  }
2818 
2819  // if (cpuDalitz(m12->getValue() + smear1, m13->getValue() + smear2, _mD0, piZeroMass, piPlusMass,
2820  // piPlusMass)) {
2821  m12->setValue(m12->getValue() + smear1);
2822  m13->setValue(m13->getValue() + smear2);
2823 
2824  //}
2825  if(100 > smearedData->getNumEvents())
2826  std::cout << " After smear: " << m12->getValue() << " " << m13->getValue() << "\n";
2827 
2828  smearedData->addEvent();
2829  }
2830 
2831  // delete data;
2832  // data = smearedData;
2833  } else
2834  smearedData = data;
2835 
2836  /*
2837  vector<Variable*> lvars;
2838  lvars.push_back(m12);
2839  lvars.push_back(m13);
2840  binEffData = new BinnedDataSet(lvars);
2841 
2842  TH2F genEff("genEff", "", m12->getNumBins(), m12->getLowerLimit(), m12->getUpperLimit(), m13->getNumBins(),
2843  m13->getLowerLimit(), m13->getUpperLimit());
2844  TH2F resEff("resEff", "", m12->getNumBins(), m12->getLowerLimit(), m12->getUpperLimit(), m13->getNumBins(),
2845  m13->getLowerLimit(), m13->getUpperLimit());
2846 
2847  double xstep = (m12->getUpperLimit() - m12->getLowerLimit()) / m12->getNumBins();
2848  double ystep = (m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins();
2849  for (int i = 0; i < m12->getNumBins(); ++i) {
2850  double lox = m12->getLowerLimit() + i*xstep;
2851  double hix = m12->getLowerLimit() + (i+1)*xstep;
2852  if (0 == i%10) std::cout << "Generating efficiency for " << i << std::endl;
2853  for (int j = 0; j < m13->getNumBins(); ++j) {
2854  double loy = m13->getLowerLimit() + j*ystep;
2855  double hiy = m13->getLowerLimit() + (j+1)*ystep;
2856 
2857  bool corner = false;
2858  if (cpuDalitz(lox, loy, _mD0, piZeroMass, piPlusMass, piPlusMass)) corner = true;
2859  else if (cpuDalitz(lox, hiy, _mD0, piZeroMass, piPlusMass, piPlusMass)) corner = true;
2860  else if (cpuDalitz(hix, loy, _mD0, piZeroMass, piPlusMass, piPlusMass)) corner = true;
2861  else if (cpuDalitz(hix, hiy, _mD0, piZeroMass, piPlusMass, piPlusMass)) corner = true;
2862  if (!corner) continue;
2863 
2864  for (int k = 0; k < 10000; ++k) {
2865  double genx = donram.Uniform(lox, hix);
2866  double geny = donram.Uniform(loy, hiy);
2867  if (!cpuDalitz(genx, geny, _mD0, piZeroMass, piPlusMass, piPlusMass)) continue;
2868  genEff.Fill(genx, geny);
2869  if (DplotRes & genResolutions) {
2870  genx += donram.Gaus(0, dplotres);
2871  geny += donram.Gaus(0, dplotres);
2872  }
2873  if (!cpuDalitz(genx, geny, _mD0, piZeroMass, piPlusMass, piPlusMass)) continue;
2874  resEff.Fill(genx, geny);
2875  }
2876  }
2877  }
2878 
2879  for (int i = 0; i < m12->getNumBins(); ++i) {
2880  for (int j = 0; j < m13->getNumBins(); ++j) {
2881  double gen = genEff.GetBinContent(i+1, j+1);
2882  if (0.1 > gen) continue;
2883  double res = resEff.GetBinContent(i+1, j+1);
2884  m12->getValue() = m12->getLowerLimit() + (i+0.5)*xstep;
2885  m13->getValue() = m13->getLowerLimit() + (j+0.5)*ystep;
2886  binEffData->setBinContent(binEffData->getBinNumber(), (res/gen));
2887  resEff.SetBinContent(i+1, j+1, (res/gen));
2888  }
2889  }
2890 
2891  Variable* effSmoothing = new Variable("effSmoothing", 0.0);
2892  SmoothHistogramPdf* eff = new SmoothHistogramPdf("efficiency", binEffData, effSmoothing);
2893  foodal->cd();
2894  resEff.SetStats(false);
2895  foodal->SetLogz(true);
2896  resEff.GetZaxis()->SetRangeUser(0.99, 1);
2897  resEff.Draw("colz");
2898  foodal->SaveAs("gen_efficiency.png");
2899  foodal->SetLogz(false);
2900  foo->cd();
2901  */
2902 
2903  int oldBins1 = m12->getNumBins();
2904  int oldBins2 = m13->getNumBins();
2905  m12->setNumBins(120);
2906  m13->setNumBins(120);
2907 
2909  weightHistogram->Draw("colz");
2910  foodal->SaveAs("./plots_from_mixfit/efficiency_weights.png");
2911 
2912  vector<Observable> lvars;
2913  lvars.push_back(*m12);
2914  lvars.push_back(*m13);
2915  binEffData = new BinnedDataSet(lvars);
2916  fname = app_ptr->get_filename("dataFiles/efficiency_gen.txt", "examples/pipipi0DPFit");
2917  loadDataFile(fname, &effdata, 1);
2918  GooPdf *eff = makeEfficiencyPdf();
2919  m12->setNumBins(oldBins1);
2920  m13->setNumBins(oldBins2);
2921 
2922  // eff->setData(effdata);
2923  // FitManager effpdf(eff);
2924  // effpdf.fit();
2925  // eff->setParameterConstantness(true);
2926  // binEffData = 0;
2927  // delete effdata; effdata = 0;
2928 
2929  TruthResolution *res = new TruthResolution();
2930  signalDalitz = makeSignalPdf(res, eff);
2931  signalDalitz->setDataSize(smearedData->getNumEvents()); // Default 5 is ok here, no event weighting
2932  signalDalitz->setData(smearedData);
2933 
2934  /*
2935  std::vector<std::vector<double> > pdfValues1;
2936  signalDalitz->getCompProbsAtDataPoints(pdfValues1);
2937 
2938  signalDalitz->setDataSize(smearedData->getNumEvents()); // Default 5 is ok here, no event weighting
2939  signalDalitz->setData(smearedData);
2940  std::vector<std::vector<double> > pdfValues2;
2941  signalDalitz->getCompProbsAtDataPoints(pdfValues2);
2942 
2943  double nll1 = 0;
2944  double nll2 = 0;
2945  for (int i = 0; i < data->getNumEvents(); ++i) {
2946  nll1 += log(pdfValues1[0][i]);
2947  nll2 += log(pdfValues2[0][i]);
2948  data->loadEvent(i);
2949  if ((100 > i) || (fabs(pdfValues1[0][i] - pdfValues2[0][i]) > 0.5)) {
2950  double eff1 = binEffData->getBinContent(binEffData->getBinNumber());
2951  std::cout << i << ": " << m12->getValue() << " " << m13->getValue() << " -> " << pdfValues1[0][i] << " " << eff1
2952  << " " << binEffData->getBinNumber();
2953  smearedData->loadEvent(i);
2954  eff1 = binEffData->getBinContent(binEffData->getBinNumber());
2955  std::cout << " | " << m12->getValue() << " " << m13->getValue() << " -> " << pdfValues2[0][i] << " " << eff1
2956  << " " << binEffData->getBinNumber();
2957  std::cout << std::endl;
2958  }
2959  }
2960  std::cout << "Final NLLs: " << nll1 << " " << nll2 << std::endl;
2961  */
2962 
2963  int retval;
2964  if(minuit1) {
2965  GooFit::FitManagerMinuit1 datapdf(signalDalitz);
2966  datapdf.setMaxCalls(64000);
2967  datapdf.fit();
2968  retval = datapdf;
2969  } else {
2970  GooFit::FitManagerMinuit2 datapdf(signalDalitz);
2971  datapdf.setMaxCalls(64000);
2972  datapdf.fit();
2973  retval = datapdf;
2974  }
2975 
2976  fmt::print("Fit results Canonical fit:\n"
2977  "tau : {:.3}\n"
2978  "xmixing: ({:.3})\%\n"
2979  "ymixing: ({:.3})\%\n",
2981  100 * Uncertain(dtop0pp._xmixing),
2982  100 * Uncertain(dtop0pp._ymixing));
2983 
2984  // All this relies on exact formatting of the input data files; it's fragile.
2985  double inputx = 1;
2986  double inputy = 1;
2987  std::string::size_type pos = fname.find("mm");
2988 
2989  if(pos != std::string::npos)
2990  inputx = inputy = -1;
2991  else {
2992  pos = fname.find("mp");
2993 
2994  if(pos != std::string::npos)
2995  inputx = -1;
2996  else {
2997  pos = fname.find("pm");
2998 
2999  if(pos != std::string::npos)
3000  inputy = -1;
3001  else {
3002  pos = fname.find("pp");
3003  assert(pos != std::string::npos);
3004  }
3005  }
3006  }
3007 
3008  std::string ident = fname.substr(pos, 4);
3009  sprintf(strbuffer, "result_%s_%f", ident.c_str(), dplotres);
3010  ofstream writer;
3011  writer.open(strbuffer);
3012  writer << inputx << " " << 100 * dtop0pp._xmixing.getValue() << " " << 100 * dtop0pp._xmixing.getError() << " "
3013  << inputy << " " << 100 * dtop0pp._ymixing.getValue() << " " << 100 * dtop0pp._ymixing.getError()
3014  << std::endl;
3015  writer.close();
3016 
3017  return retval;
3018  // makeDalitzPlots(signalDalitz, "plots_from_mixfit/generated/");
3019 }
3020 
3022  // Variable* bkg2_sigma_js_meana = new Variable("bkg2_sigma_js_meana", 0.01, 0.01, -0.30, 0.30);
3023  // Variable* bkg2_sigma_js_sigma = new Variable("bkg2_sigma_js_sigma", 0.09, 0.01, 0, 0.4);
3024  // Variable* bkg2_sigma_js_gamma = new Variable("bkg2_sigma_js_gamma",-5.00, 0.10, -30, 0);
3025  // Variable* bkg2_sigma_js_delta = new Variable("bkg2_sigma_js_delta", 1.49, 0.01, 0.50, 5.00);
3026  // Variable* bkg2_sigma_frac_jsu = new Variable("bkg2_sigma_frac_jsu", 0.85, 0.01, 0.01, 1.00);
3027  // Variable* bkg2_sigma_frac_ga1 = new Variable("bkg2_sigma_frac_ga1", 0.04, 0.01, 0.01, 0.20);
3028 
3029  Variable bkg2_sigma_num_jsu("bkg2_sigma_num_jsu", 9100, 200, 1000, 15000);
3030  Variable bkg2_sigma_num_ga1("bkg2_sigma_num_ga1", 2400, 200, 500, 7000);
3031  Variable bkg2_sigma_num_ga2("bkg2_sigma_num_ga2", 2900, 200, 500, 7000);
3032 
3033  Variable bkg2_sigma_g1_meana("bkg2_sigma_g1_meana", 0.35, 0.01, 0.10, 0.50);
3034  Variable bkg2_sigma_g1_sigma("bkg2_sigma_g1_sigma", 0.30, 0.01, 0.05, 0.55);
3035  Variable bkg2_sigma_g2_meana("bkg2_sigma_g2_meana", 0.80, 0.01, 0.01, 1.50);
3036  Variable bkg2_sigma_g2_sigma("bkg2_sigma_g2_sigma", 0.90, 0.01, 0.01, 2.75);
3037  // JohnsonSUPdf* bkg2_sigma_js = new JohnsonSUPdf("bkg2_sigma_js", sigma, bkg2_sigma_js_meana, bkg2_sigma_js_sigma,
3038  // bkg2_sigma_js_gamma, bkg2_sigma_js_delta);
3039 
3040  Variable bkg2_sigma_js_meana("bkg2_sigma_js_meana", 0.35, 0.01, 0.00, 0.60);
3041  Variable bkg2_sigma_js_sigma("bkg2_sigma_js_sigma", 0.09, 0.01, 0, 0.4);
3042  Variable bkg2_sigma_js_gamma("bkg2_sigma_js_gamma", 2.00, 0.10, 0, 10);
3043  Variable bkg2_sigma_js_delta("bkg2_sigma_js_delta", 2);
3044  CrystalBallPdf *bkg2_sigma_js = new CrystalBallPdf(
3045  "bkg2_sigma_js", *sigma, bkg2_sigma_js_meana, bkg2_sigma_js_sigma, bkg2_sigma_js_gamma, bkg2_sigma_js_delta);
3046 
3047  GaussianPdf *bkg2_sigma_g1 = new GaussianPdf("bkg2_sigma_g1", *sigma, bkg2_sigma_g1_meana, bkg2_sigma_g1_sigma);
3048  GaussianPdf *bkg2_sigma_g2 = new GaussianPdf("bkg2_sigma_g2", *sigma, bkg2_sigma_g2_meana, bkg2_sigma_g2_sigma);
3049 
3050  weights.clear();
3051  weights.push_back(bkg2_sigma_num_jsu);
3052  weights.push_back(bkg2_sigma_num_ga1);
3053  weights.push_back(bkg2_sigma_num_ga2);
3054  comps.clear();
3055  comps.push_back(bkg2_sigma_js);
3056  comps.push_back(bkg2_sigma_g1);
3057  comps.push_back(bkg2_sigma_g2);
3058 
3059  GooPdf *ret = new AddPdf("bkg2_jsugg", weights, comps);
3060  return ret;
3061  // return bkg2_sigma_js;
3062 }
3063 
3065  // Variable* bkg4_sigma_js_meana = new Variable("bkg4_sigma_js_meana", 0.01, 0.01, -0.30, 0.30);
3066  // Variable* bkg4_sigma_js_sigma = new Variable("bkg4_sigma_js_sigma", 0.09, 0.01, 0, 0.4);
3067  // Variable* bkg4_sigma_js_gamma = new Variable("bkg4_sigma_js_gamma",-5.00, 0.10, -30, 0);
3068  // Variable* bkg4_sigma_js_delta = new Variable("bkg4_sigma_js_delta", 1.49, 0.01, 0.50, 5.00);
3069  // Variable* bkg4_sigma_frac_jsu = new Variable("bkg4_sigma_frac_jsu", 0.85, 0.01, 0.01, 1.00);
3070  // Variable* bkg4_sigma_frac_ga1 = new Variable("bkg4_sigma_frac_ga1", 0.04, 0.01, 0.01, 0.20);
3071 
3072  Variable bkg4_sigma_num_jsu("bkg4_sigma_num_jsu", 9100, 200, 1000, 15000);
3073  Variable bkg4_sigma_num_ga1("bkg4_sigma_num_ga1", 2400, 200, 500, 7000);
3074  Variable bkg4_sigma_num_ga2("bkg4_sigma_num_ga2", 2900, 200, 500, 7000);
3075 
3076  Variable bkg4_sigma_g1_meana("bkg4_sigma_g1_meana", 0.35, 0.01, 0.10, 0.50);
3077  Variable bkg4_sigma_g1_sigma("bkg4_sigma_g1_sigma", 0.30, 0.01, 0.05, 0.55);
3078  Variable bkg4_sigma_g2_meana("bkg4_sigma_g2_meana", 0.80, 0.01, 0.01, 1.50);
3079  Variable bkg4_sigma_g2_sigma("bkg4_sigma_g2_sigma", 0.90, 0.01, 0.01, 2.75);
3080  // JohnsonSUPdf* bkg4_sigma_js = new JohnsonSUPdf("bkg4_sigma_js", sigma, bkg4_sigma_js_meana, bkg4_sigma_js_sigma,
3081  // bkg4_sigma_js_gamma, bkg4_sigma_js_delta);
3082 
3083  Variable bkg4_sigma_js_meana("bkg4_sigma_js_meana", 0.35, 0.01, 0.00, 0.60);
3084  Variable bkg4_sigma_js_sigma("bkg4_sigma_js_sigma", 0.09, 0.01, 0, 0.4);
3085  Variable bkg4_sigma_js_gamma("bkg4_sigma_js_gamma", 2.00, 0.10, 0, 10);
3086  Variable bkg4_sigma_js_delta("bkg4_sigma_js_delta", 2);
3087  CrystalBallPdf *bkg4_sigma_js = new CrystalBallPdf(
3088  "bkg4_sigma_js", *sigma, bkg4_sigma_js_meana, bkg4_sigma_js_sigma, bkg4_sigma_js_gamma, bkg4_sigma_js_delta);
3089 
3090  GaussianPdf *bkg4_sigma_g1 = new GaussianPdf("bkg4_sigma_g1", *sigma, bkg4_sigma_g1_meana, bkg4_sigma_g1_sigma);
3091  GaussianPdf *bkg4_sigma_g2 = new GaussianPdf("bkg4_sigma_g2", *sigma, bkg4_sigma_g2_meana, bkg4_sigma_g2_sigma);
3092 
3093  weights.clear();
3094  weights.push_back(bkg4_sigma_num_jsu);
3095  weights.push_back(bkg4_sigma_num_ga1);
3096  weights.push_back(bkg4_sigma_num_ga2);
3097  comps.clear();
3098  comps.push_back(bkg4_sigma_js);
3099  comps.push_back(bkg4_sigma_g1);
3100  comps.push_back(bkg4_sigma_g2);
3101 
3102  GooPdf *ret = new AddPdf("bkg4_jsugg", weights, comps);
3103  return ret;
3104 }
3105 
3107  // Variable* bkg3_sigma_js_meana = new Variable("bkg3_sigma_js_meana", 0.05, 0.01, -0.30, 0.30);
3108  // Variable* bkg3_sigma_js_sigma = new Variable("bkg3_sigma_js_sigma", 0.013, 0.01, 0, 0.2);
3109  // Variable* bkg3_sigma_js_gamma = new Variable("bkg3_sigma_js_gamma",-6.00, 1.00, -30, 0);
3110  // Variable* bkg3_sigma_js_delta = new Variable("bkg3_sigma_js_delta", 1.99, 0.10, 0.50, 5.00);
3111  Variable bkg3_sigma_frac_jsu("bkg3_sigma_frac_jsu", 0.50, 0.01, 0.01, 1.00);
3112  Variable bkg3_sigma_frac_ga1("bkg3_sigma_frac_ga1", 0.04, 0.01, 0.01, 0.20);
3113 
3114  // Variable* bkg3_sigma_num_jsu = new Variable("bkg3_sigma_num_jsu", 11000, 200, 1000, 35000);
3115  // Variable* bkg3_sigma_num_ga1 = new Variable("bkg3_sigma_num_ga1", 9400, 200, 500, 10000);
3116  // Variable* bkg3_sigma_num_ga2 = new Variable("bkg3_sigma_num_ga2", 3900, 200, 500, 17000);
3117  // Variable* bkg3_sigma_num_ga3 = new Variable("bkg3_sigma_num_ga3", 3900, 200, 500, 7000);
3118 
3119  Variable bkg3_sigma_g1_meana("bkg3_sigma_g1_meana", 0.35, 0.01, 0.10, 0.50);
3120  Variable bkg3_sigma_g1_sigma("bkg3_sigma_g1_sigma", 0.10, 0.01, 0.01, 0.55);
3121  Variable bkg3_sigma_g2_meana("bkg3_sigma_g2_meana", 0.20, 0.01, 0.01, 1.50);
3122  Variable bkg3_sigma_g2_sigma("bkg3_sigma_g2_sigma", 0.10, 0.01, 0.001, 0.15);
3123  Variable bkg3_sigma_g2_gamma("bkg3_sigma_g2_gamma", -2.00, 1.00, -10, 10);
3124  Variable bkg3_sigma_g2_delta("bkg3_sigma_g2_delta", 2, 0.10, 0.50, 5.00);
3125  // Variable* bkg3_sigma_g3_meana = new Variable("bkg3_sigma_g3_meana", 0.20, 0.01, 0.01, 1.50);
3126  // Variable* bkg3_sigma_g3_sigma = new Variable("bkg3_sigma_g3_sigma", 0.20, 0.01, 0.01, 0.75);
3127  // Variable* bkg3_sigma_g2_width = new Variable("bkg3_sigma_g2_width", 0.10, 0.01, 0.01, 0.75);
3128  // JohnsonSUPdf* bkg3_sigma_js = new JohnsonSUPdf("bkg3_sigma_js", sigma, bkg3_sigma_js_meana, bkg3_sigma_js_sigma,
3129  // bkg3_sigma_js_gamma, bkg3_sigma_js_delta);
3130 
3131  Variable bkg3_sigma_js_meana("bkg3_sigma_js_meana", 0.35, 0.01, 0.00, 0.60);
3132  Variable bkg3_sigma_js_sigma("bkg3_sigma_js_sigma", 0.09, 0.01, 0, 0.40);
3133  Variable bkg3_sigma_js_gamma("bkg3_sigma_js_gamma", 2.00, 0.10, 0, 10);
3134  Variable bkg3_sigma_js_delta("bkg3_sigma_js_delta", 2);
3135  CrystalBallPdf *bkg3_sigma_js = new CrystalBallPdf(
3136  "bkg3_sigma_js", *sigma, bkg3_sigma_js_meana, bkg3_sigma_js_sigma, bkg3_sigma_js_gamma, bkg3_sigma_js_delta);
3137  // JohnsonSUPdf* bkg3_sigma_js = new JohnsonSUPdf("bkg3_sigma_js", sigma, bkg3_sigma_js_meana, bkg3_sigma_js_sigma,
3138  // bkg3_sigma_js_gamma, bkg3_sigma_js_delta);
3139 
3140  GaussianPdf *bkg3_sigma_g1 = new GaussianPdf("bkg3_sigma_g1", *sigma, bkg3_sigma_g1_meana, bkg3_sigma_g1_sigma);
3141  // GaussianPdf* bkg3_sigma_g2 = new GaussianPdf("bkg3_sigma_g2", sigma, bkg3_sigma_g2_meana, bkg3_sigma_g2_sigma);
3142  // CrystalBallPdf* bkg3_sigma_g2 = new CrystalBallPdf("bkg3_sigma_g2", sigma, bkg3_sigma_g2_meana,
3143  // bkg3_sigma_g2_sigma, bkg3_sigma_g2_gamma, bkg3_sigma_g2_delta);
3144  JohnsonSUPdf *bkg3_sigma_g2 = new JohnsonSUPdf(
3145  "bkg3_sigma_g2", *sigma, bkg3_sigma_g2_meana, bkg3_sigma_g2_sigma, bkg3_sigma_g2_gamma, bkg3_sigma_g2_delta);
3146  // GaussianPdf* bkg3_sigma_g3 = new GaussianPdf("bkg3_sigma_g3", sigma, bkg3_sigma_g3_meana, bkg3_sigma_g3_sigma);
3147  // VoigtianPdf* bkg3_sigma_g2 = new VoigtianPdf("bkg3_sigma_g2", sigma, bkg3_sigma_g2_meana, bkg3_sigma_g2_sigma,
3148  // bkg3_sigma_g2_width);
3149 
3150  weights.clear();
3151  // weights.push_back(bkg3_sigma_num_jsu);
3152  // weights.push_back(bkg3_sigma_num_ga1);
3153  // weights.push_back(bkg3_sigma_num_ga2);
3154  // weights.push_back(bkg3_sigma_num_ga3);
3155  weights.push_back(bkg3_sigma_frac_jsu);
3156  weights.push_back(bkg3_sigma_frac_ga1);
3157  comps.clear();
3158  comps.push_back(bkg3_sigma_js);
3159  comps.push_back(bkg3_sigma_g1);
3160  comps.push_back(bkg3_sigma_g2);
3161  // comps.push_back(bkg3_sigma_g3);
3162 
3163  GooPdf *ret = new AddPdf("bkg3_jsugg", weights, comps);
3164  return ret;
3165 }
3166 
3168  // Gaussians for decay time.
3169  Variable *frac_ga2;
3170  Variable *frac_ga3;
3171  Variable *g1_meana;
3172  Variable *g1_sigma;
3173  Variable *g2_meana;
3174  Variable *g2_sigma;
3175  Variable *g3_meana;
3176  Variable *g3_sigma;
3177 
3178  std::string bkgname = "";
3179 
3180  switch(bkg) {
3181  case 4:
3182  frac_ga2 = new Variable("bkg4_frac_ga2", 0.60780, 0.01, 0.40, 0.80);
3183  frac_ga3 = new Variable("bkg4_frac_ga3", 0.04776, 0.01, 0.00, 0.20);
3184  g1_meana = new Variable("bkg4_g1_meana", 0.10164, 0.01, 0.00, 0.80);
3185  g1_sigma = new Variable("bkg4_g1_sigma", 0.27504, 0.01, 0.10, 0.80);
3186  g2_meana = new Variable("bkg4_g2_meana", 0.19974, 0.01, 0.10, 0.80);
3187  g2_sigma = new Variable("bkg4_g2_sigma", 0.63765, 0.01, 0.40, 0.80);
3188  g3_meana = new Variable("bkg4_g3_meana", 0.45817, 0.01, 0.20, 0.80);
3189  g3_sigma = new Variable("bkg4_g3_sigma", 1.52905, 0.01, 1.40, 1.80);
3190  bkgname = "bkg4";
3191  break;
3192 
3193  case 3:
3194  frac_ga2 = new Variable("bkg3_frac_ga2", 0.51448, 0.01, 0.25, 0.75);
3195  frac_ga3 = new Variable("bkg3_frac_ga3", 0.04169, 0.01, 0.00, 0.40);
3196  g1_meana = new Variable("bkg3_g1_meana", 0.25101, 0.01, 0.01, 0.50);
3197  g1_sigma = new Variable("bkg3_g1_sigma", 0.31953, 0.01, 0.02, 0.50);
3198  g2_meana = new Variable("bkg3_g2_meana", 0.49139, 0.01, 0.25, 0.75);
3199  g2_sigma = new Variable("bkg3_g2_sigma", 0.65443, 0.01, 0.10, 1.00);
3200  g3_meana = new Variable("bkg3_g3_meana", 0.83600, 0.01, 0.50, 1.00);
3201  g3_sigma = new Variable("bkg3_g3_sigma", 1.51839, 0.01, 0.10, 2.00);
3202  bkgname = "bkg3";
3203  break;
3204 
3205  case 2:
3206  default:
3207  frac_ga2 = new Variable("frac_ga2", 0.48994, 0.01, 0.1, 0.6);
3208  frac_ga3 = new Variable("frac_ga3", 0.04721, 0.01, 0.0, 0.1);
3209  g1_meana = new Variable("g1_meana", 0.01216, 0.01, -0.1, 0.1);
3210  g1_sigma = new Variable("g1_sigma", 0.25813, 0.01, 0.15, 0.35);
3211  g2_meana = new Variable("g2_meana", 0.05335, 0.01, -0.1, 0.1);
3212  g2_sigma = new Variable("g2_sigma", 0.58651, 0.01, 0.5, 1.2);
3213  g3_meana = new Variable("g3_meana", 0.17451, 0.01, 0.1, 1.85);
3214  g3_sigma = new Variable("g3_sigma", 1.15125, 0.01, 0.5, 1.3);
3215  bkgname = "bkg2";
3216  break;
3217  }
3218 
3219  GaussianPdf *g1 = new GaussianPdf((bkgname + "_g1").c_str(), *dtime, *g1_meana, *g1_sigma);
3220  GaussianPdf *g2 = new GaussianPdf((bkgname + "_g2").c_str(), *dtime, *g2_meana, *g2_sigma);
3221  GaussianPdf *g3 = new GaussianPdf((bkgname + "_g3").c_str(), *dtime, *g3_meana, *g3_sigma);
3222 
3223  weights.clear();
3224  weights.push_back(*frac_ga2);
3225 
3226  if(3 == bkg)
3227  weights.push_back(*frac_ga3);
3228 
3229  comps.clear();
3230  comps.push_back(g2);
3231 
3232  if(3 == bkg)
3233  comps.push_back(g3);
3234 
3235  comps.push_back(g1);
3236  AddPdf *bkg_dtime = new AddPdf((bkgname + "_dtime").c_str(), weights, comps);
3237  return bkg_dtime;
3238 }
3239 
3241  std::string bkgname = "";
3242 
3243  switch(bkg) {
3244  case 4:
3245  bkgname = "bkg4";
3246  break;
3247 
3248  case 3:
3249  bkgname = "bkg3";
3250  break;
3251 
3252  case 2:
3253  default:
3254  bkgname = "bkg2";
3255  break;
3256  }
3257 
3258  Variable g1_mean((bkgname + "_dtime_gmean1"), 0, 0.01, -0.5, 0.5);
3259  Variable g1_sigm((bkgname + "_dtime_gsigm1"), 0.2, 0.01, 0.01, 0.8);
3260  Variable e1_alph((bkgname + "_dtime_alpha1"), 2.5, 0.01, 0.01, 7.5);
3261 
3262  Variable g2_mean((bkgname + "_dtime_gmean2"), -0.3, 0.01, -0.85, 0.85);
3263  Variable g2_sigm((bkgname + "_dtime_gsigm2"), 0.2, 0.01, 0.01, 0.8);
3264  Variable e2_alph((bkgname + "_dtime_alpha2"), 0.5, 0.01, 0.01, 10.0);
3265 
3266  ExpGausPdf *exp1 = new ExpGausPdf((bkgname + "_exp1"), *dtime, g1_mean, g1_sigm, e1_alph);
3267  ExpGausPdf *exp2 = new ExpGausPdf((bkgname + "_exp2"), *dtime, g2_mean, g2_sigm, e2_alph);
3268 
3269  Variable frac1((bkgname + "_dtime_frac1"), 0.1, 0.01, 0, 0.8);
3270 
3271  GooPdf *ret = new AddPdf((bkgname + "_dtime"), frac1, exp1, exp2);
3272 
3273  return ret;
3274 }
3275 
3276 GooPdf *makeBkg2DalitzPdf(bool fixem = true) {
3277  if(!kzero_veto)
3278  makeKzeroVeto();
3279 
3280  GooPdf *bkg2_dalitz = nullptr;
3281 
3282  if(Parameter == bkg2Model) {
3283  comps.clear();
3284 
3285  vector<Variable> offsets;
3286  vector<Observable> observables;
3287  vector<Variable> coefficients;
3288  offsets.push_back(constantOne);
3289  offsets.push_back(constantOne);
3290  observables.push_back(*m12);
3291  observables.push_back(*m13);
3292  double weightOffset = 3;
3293  // Recurring factor 3 offsets division by total weight in AddPdf.
3294  coefficients.push_back(Variable("bkg2_x0y0", 1.0 * weightOffset));
3295  coefficients.push_back(
3296  Variable("bkg2_x1y0", 0.13184 * weightOffset, 0.01, 0.01 * weightOffset, 0.18 * weightOffset));
3297  coefficients.push_back(
3298  Variable("bkg2_x2y0", 0.02062 * weightOffset, 0.01, 0.00 * weightOffset, 0.17 * weightOffset));
3299  coefficients.push_back(
3300  Variable("bkg2_x3y0", 0.04688 * weightOffset, 0.01, 0.00 * weightOffset, 0.08 * weightOffset));
3301  coefficients.push_back(
3302  Variable("bkg2_x0y1", -0.02568 * weightOffset, 0.01, -0.15 * weightOffset, 0.04 * weightOffset));
3303  coefficients.push_back(
3304  Variable("bkg2_x1y1", 0.06805 * weightOffset, 0.01, 0.02 * weightOffset, 0.10 * weightOffset));
3305  coefficients.push_back(
3306  Variable("bkg2_x2y1", 0.38557 * weightOffset, 0.01, 0.30 * weightOffset, 0.50 * weightOffset));
3307  coefficients.push_back(
3308  Variable("bkg2_x0y2", 0.11252 * weightOffset, 0.01, 0.05 * weightOffset, 0.20 * weightOffset));
3309  coefficients.push_back(
3310  Variable("bkg2_x1y2", 0.24896 * weightOffset, 0.01, 0.20 * weightOffset, 0.50 * weightOffset));
3311  coefficients.push_back(
3312  Variable("bkg2_x0y3", 0.05605 * weightOffset, 0.01, -0.05 * weightOffset, 0.15 * weightOffset));
3313 
3314  PolynomialPdf *poly = new PolynomialPdf("bkg2Pdf", observables, coefficients, offsets, 3);
3315 
3316  Variable bkg2_decZmin("bkg2_decZmin", 3.30411);
3317  Variable bkg2_conZmin("bkg2_conZmin", 0.29909);
3318  TrigThresholdPdf *bkg2_loZ
3319  = new TrigThresholdPdf("bkg2_loZ", *m12, *m13, minDalitzZ, bkg2_decZmin, bkg2_conZmin, massSum, false);
3320  // bkg2_loZ->setDebugMask(1);
3321 
3322  comps.clear();
3323  comps.push_back(poly);
3324  comps.push_back(bkg2_loZ);
3325  comps.push_back(kzero_veto);
3326  // Separate PDF to avoid triggering numerical normalisation over all four observables.
3327  ProdPdf *poly_x_veto = new ProdPdf("poly_x_veto", comps);
3328  // poly_x_veto->setDebugMask(1);
3329 
3330  // One omega->pipipi0 reflection.
3331  // Factor 3 in amplitudes is to offset division by total weight in AddPdf.
3332  DecayInfo3 special_rho_decay;
3333  special_rho_decay.motherMass = _mD0;
3334  special_rho_decay.daug1Mass = piZeroMass;
3335  special_rho_decay.daug2Mass = piPlusMass;
3336  special_rho_decay.daug3Mass = piPlusMass;
3337  special_rho_decay.meson_radius = 1.5;
3338 
3339  ResonancePdf *bkg2_rho_ref = new Resonances::Gauss(
3340  "bkg2_rho_ref",
3341  Variable("bkg2_rho_ref_amp", 0.00896 * weightOffset, 0.001, 0, 0.015 * weightOffset),
3342  constantZero,
3343  Variable("bkg2_rho_ref_mass", 0.53172),
3344  Variable("bkg2_rho_ref_width", 0.06426),
3345  PAIR_13);
3346  special_rho_decay.resonances.push_back(bkg2_rho_ref);
3347 
3348  Variable bkg2_rho_poly_offset("bkg2_rho_poly_offset", 1.64254);
3349  Variable bkg2_rho_poly_linear("bkg2_rho_poly_linear", 0);
3350  Variable bkg2_rho_poly_second("bkg2_rho_poly_second", -0.48166);
3351  weights.clear();
3352  weights.push_back(constantOne);
3353  weights.push_back(bkg2_rho_poly_linear);
3354  weights.push_back(bkg2_rho_poly_second);
3355  PolynomialPdf *bkg2_rho_poly = new PolynomialPdf("bkg2_rho_poly", *m12, weights, bkg2_rho_poly_offset);
3356  comps.clear();
3357  comps.push_back(kzero_veto);
3358  comps.push_back(bkg2_rho_poly);
3359  comps.push_back(bkg2_loZ);
3360  ProdPdf *bkg2_rho_mods = new ProdPdf("bkg2_rho_mods", comps);
3361  incsum1 = new IncoherentSumPdf("incsum1", *m12, *m13, *eventNumber, special_rho_decay, bkg2_rho_mods);
3362 
3363  // Three spin-0 rho resonances to be added incoherently.
3364  DecayInfo3 incoherent_rho0s;
3365  incoherent_rho0s.motherMass = _mD0;
3366  incoherent_rho0s.daug1Mass = piZeroMass;
3367  incoherent_rho0s.daug2Mass = piPlusMass;
3368  incoherent_rho0s.daug3Mass = piPlusMass;
3369  incoherent_rho0s.meson_radius = 0; // Mikhail uses zero radius for incoherent resonances.
3370 
3371  ResonancePdf *bkg2_incRho0 = new Resonances::RBW(
3372  "bkg2_incRho0",
3373  Variable("bkg2_incRho0_amp", 0.00304 * weightOffset, 0.001, 0.0, 0.006 * weightOffset),
3374  constantZero,
3375  fixedRhoMass,
3376  fixedRhoWidth,
3377  0, // Incoherent rho has effective spin 0.
3378  PAIR_23);
3379  incoherent_rho0s.resonances.push_back(bkg2_incRho0);
3380 
3381  ResonancePdf *bkg2_incRhoP = new Resonances::RBW(
3382  "bkg2_incRhoP",
3383  Variable("bkg2_incRhoP_amp", 0.00586 * weightOffset, 0.001, 0.0, 0.012 * weightOffset),
3384  constantZero,
3385  fixedRhoMass,
3386  fixedRhoWidth,
3387  0,
3388  PAIR_12);
3389  incoherent_rho0s.resonances.push_back(bkg2_incRhoP);
3390 
3391  ResonancePdf *bkg2_incRhoM = new Resonances::RBW(
3392  "bkg2_incRhoM",
3393  Variable("bkg2_incRhoM_amp", 0.00635 * weightOffset, 0.001, 0.0, 0.015 * weightOffset),
3394  constantZero,
3395  fixedRhoMass,
3396  fixedRhoWidth,
3397  0,
3398  PAIR_13);
3399  incoherent_rho0s.resonances.push_back(bkg2_incRhoM);
3400 
3401  comps.clear();
3402  comps.push_back(kzero_veto);
3403  comps.push_back(bkg2_loZ);
3404  ProdPdf *bkg2_rho_mods2 = new ProdPdf("bkg2_rho_mods2", comps);
3405 
3406  incsum2 = new IncoherentSumPdf("incsum2", *m12, *m13, *eventNumber, incoherent_rho0s, bkg2_rho_mods2);
3407 
3408  weights.clear();
3409  weights.push_back(constantOne);
3410  weights.push_back(constantOne);
3411  weights.push_back(constantOne);
3412  comps.clear();
3413  comps.push_back(poly_x_veto);
3414  comps.push_back(incsum1);
3415  comps.push_back(incsum2);
3416 
3417  bkg2_dalitz = new AddPdf("bkg2_dalitz", weights, comps);
3418  bkg2_dalitz->addSpecialMask(PdfBase::ForceCommonNorm);
3419  // bkg2_dalitz->setDebugMask(1);
3420 
3421  } else if(Histogram == bkg2Model) {
3422  bkg2_dalitz = makeBackgroundHistogram(2);
3423  } else if(Sideband == bkg2Model) {
3424  comps.clear();
3425  // Originally had massD0 sidebands separated into deltaM high and low,
3426  // but these distributions were extremely similar - collapsed them
3427  // into just massD0 sidebands.
3428  std::string fname = app_ptr->get_filename("dataFiles/sideband1.txt", "examples/pipipi0DPFit");
3429  comps.push_back(makeBackgroundHistogram(101, fname));
3430  fname = app_ptr->get_filename("dataFiles/sideband2.txt", "examples/pipipi0DPFit");
3431  comps.push_back(makeBackgroundHistogram(102, fname));
3432  // comps.push_back(makeBackgroundHistogram(103, "./dataFiles/sideband3.txt"));
3433  // comps.push_back(makeBackgroundHistogram(104, "./dataFiles/sideband4.txt"));
3434  weights.clear();
3435  weights.push_back(Variable("sband1Weight", 300000, 1000, 100, 750000));
3436  weights.push_back(Variable("sband2Weight", 100000, 1000, 100, 500000));
3437  // weights.push_back(new Variable("sband3Weight", 150000, 1000, 100, 500000));
3438  // weights.push_back(new Variable("sband4Weight", 50000, 1000, 100, 500000));
3439  bkg2_dalitz = new AddPdf("bkg2_dalitz", weights, comps);
3440  } else {
3441  // This cannot possibly happen, and if it does something is wrong. Panic!
3442  assert(Sideband == bkg2Model);
3443  }
3444 
3445  GooPdf *bkg2_dtime = nullptr;
3446 
3447  if(gaussBkgTime)
3448  bkg2_dtime = makeGaussianTimePdf(2);
3449  else
3450  bkg2_dtime = makeExpGausTimePdf(2);
3451 
3452  // bkg2_jsugg = sig0_jsugg; // Mikhail uses same sigma distribution for everything.
3453  // Separate sigma_t distribution
3454  // bkg2_jsugg = makeBkg2_sigma();
3455  bkg2_jsugg = makeBkg_sigma_strips(2);
3456  bkg2_jsugg->addSpecialMask(PdfBase::ForceSeparateNorm);
3457 
3458  // Finally create overall product.
3459  comps.clear();
3460  // poly->setDebugMask(1);
3461  // bkg2_dalitz->setDebugMask(1);
3462  // incsum2->setDebugMask(1);
3463  comps.push_back(bkg2_dalitz);
3464  comps.push_back(bkg2_dtime);
3465  comps.push_back(bkg2_jsugg);
3466 
3467  GooPdf *ret = new ProdPdf("bkg2_total", comps);
3468 
3469  if(fixem)
3470  ret->setParameterConstantness(true);
3471 
3472  return ret;
3473 }
3474 
3476  // Smoothed histogram from flat-file data.
3477  // Only 4500 events, so use large bins.
3478 
3479  obsweights.clear();
3480  obsweights.push_back(*m12);
3481  obsweights.push_back(*m13);
3482 
3483  int m12bins = m12->getNumBins();
3484  int m13bins = m13->getNumBins();
3485 
3486  m12->setNumBins(30);
3487  m13->setNumBins(30);
3488  BinnedDataSet *bkg3_eff_data = new BinnedDataSet(obsweights);
3489  std::ifstream reader;
3490  std::string fname = app_ptr->get_filename("dataFiles/efficiency_bkg3_flat.txt", "examples/pipipi0DPFit");
3491  readWrapper(reader, fname);
3492  std::string buffer;
3493 
3494  while(!reader.eof()) {
3495  reader >> buffer;
3496 
3497  if(buffer == "====")
3498  break;
3499 
3500  std::cout << buffer;
3501  }
3502 
3503  double dummy = 0;
3504 
3505  while(!reader.eof()) {
3506  reader >> dummy;
3507 
3508  if(reader.eof())
3509  break;
3510 
3511  reader >> dummy; // m23, m(pi+ pi-), called m12 in processToyRoot convention.
3512  reader >> *m12; // Already swapped according to D* charge
3513  reader >> *m13;
3514 
3515  // Everything else is irrelevant for this purpose!
3516  for(int i = 0; i < 16; ++i)
3517  reader >> dummy;
3518 
3519  bkg3_eff_data->addEvent();
3520  }
3521 
3522  Variable bkg3_eff_smoothing("bkg3_eff_smoothing", 1.0, 0.1, 0, 1.25);
3523  // Variable* bkg3_eff_smoothing = new Variable("bkg3_eff_smoothing", 1.0);
3524  SmoothHistogramPdf *ret = new SmoothHistogramPdf("bkg3_efficiency", bkg3_eff_data, bkg3_eff_smoothing);
3525 
3526  m12->setNumBins(m12bins);
3527  m13->setNumBins(m13bins);
3528 
3529  return ret;
3530 }
3531 
3532 SmoothHistogramPdf *makeBackgroundHistogram(int bkgnum, std::string overridename) {
3533  std::ifstream reader;
3534  sprintf(strbuffer, "./dataFiles/bkgDalitz_%i.txt", bkgnum);
3535 
3536  if(overridename != "")
3537  sprintf(strbuffer, "%s", overridename.c_str());
3538 
3539  std::string fname = app_ptr->get_filename(strbuffer, "examples/pipipi0DPFit");
3540  readWrapper(reader, fname);
3541  std::string buffer;
3542 
3543  while(!reader.eof()) {
3544  reader >> buffer;
3545 
3546  if(buffer == "====")
3547  break;
3548 
3549  std::cout << buffer;
3550  }
3551 
3552  obsweights.clear();
3553  obsweights.push_back(*m12);
3554  obsweights.push_back(*m13);
3555  BinnedDataSet *bkg_binned_data = new BinnedDataSet(obsweights);
3556 
3557  double dummy = 0;
3558 
3559  while(!reader.eof()) {
3560  reader >> dummy;
3561  reader >> dummy; // m23, m(pi+ pi-), called m12 in processToyRoot convention.
3562  reader >> *m12; // Already swapped according to D* charge. m12 = m(pi+pi0)
3563  reader >> *m13;
3564 
3565  // Don't need the rest.
3566  for(int i = 0; i < 16; ++i)
3567  reader >> dummy;
3568 
3569  bkg_binned_data->addEvent();
3570 
3571  // std::cout << m12->getValue() << " " << m13->getValue() << std::endl;
3572  }
3573 
3574  std::cout << "Read " << bkg_binned_data->getNumEvents() << " events for background " << bkgnum << std::endl;
3575  sprintf(strbuffer, "bkg%i_dalitz_smoothing", bkgnum);
3576  Variable smoothing(strbuffer, 1);
3577 
3578  if((-1 != bkgHistRandSeed) && ((3 == bkgnum) || (4 == bkgnum))) {
3579  std::cout << "Shuffling background " << bkgnum << " histogram with random seed " << bkgHistRandSeed
3580  << std::endl;
3581  TRandom donram(bkgHistRandSeed);
3582 
3583  for(unsigned int bin = 0; bin < bkg_binned_data->getNumBins(); ++bin) {
3584  double events = bkg_binned_data->getBinContent(bin);
3585 
3586  if(1 > events)
3587  continue;
3588 
3589  double newEvents = -1;
3590 
3591  while(0 > newEvents)
3592  newEvents = donram.Gaus(events, sqrt(events));
3593 
3594  bkg_binned_data->setBinContent(bin, newEvents);
3595  }
3596  }
3597 
3598  sprintf(strbuffer, "bkg%i_dalitz", bkgnum);
3599  SmoothHistogramPdf *bkg_dalitz = new SmoothHistogramPdf(strbuffer, bkg_binned_data, smoothing);
3600  // bkg_dalitz->setDebugMask(1);
3601  return bkg_dalitz;
3602 }
3603 
3605  // I can't make this thing describe the background 3 data.
3606 
3607  // GooPdf* bkg3_eff = makeBkg3Eff();
3608  weights.clear();
3609 
3610  vector<Variable> offsets;
3611  vector<Observable> observables;
3612  vector<Variable> coefficients;
3613  offsets.push_back(constantOne);
3614  offsets.push_back(constantOne);
3615 
3616  observables.push_back(*m12);
3617  observables.push_back(*m13);
3618  // Recurring factor 3 offsets division by total weight in AddPdf.
3619  double weightOffset = 1;
3620 
3621  coefficients.push_back(Variable("bkg3_x0y0", 1.00 * weightOffset));
3622  // coefficients.push_back(new Variable("bkg3_x0y0", 1.10 * weightOffset, 0.01, 0.01 * weightOffset, 1.50 *
3623  // weightOffset));
3624  coefficients.push_back(
3625  Variable("bkg3_x1y0", -0.36937 * weightOffset, 0.01, -1.50 * weightOffset, 0.00 * weightOffset));
3626  coefficients.push_back(
3627  Variable("bkg3_x2y0", 1.36184 * weightOffset, 0.01, -0.10 * weightOffset, 1.60 * weightOffset));
3628  // coefficients.push_back(new Variable("bkg3_x3y0", -0.43177 * weightOffset, 0.01,-1.60*weightOffset,
3629  // 0.60*weightOffset));
3630  coefficients.push_back(
3631  Variable("bkg3_x0y1", -0.27691 * weightOffset, 0.01, -1.50 * weightOffset, 0.00 * weightOffset));
3632  coefficients.push_back(
3633  Variable("bkg3_x1y1", 2.16029 * weightOffset, 0.01, 0.30 * weightOffset, 4.50 * weightOffset));
3634  // coefficients.push_back(new Variable("bkg3_x2y1", -2.04133 * weightOffset, 0.01,-2.50*weightOffset,
3635  // 1.50*weightOffset));
3636  coefficients.push_back(
3637  Variable("bkg3_x0y2", 1.33100 * weightOffset, 0.01, 1.00 * weightOffset, 2.00 * weightOffset));
3638  // coefficients.push_back(new Variable("bkg3_x1y2", -1.88226 * weightOffset, 0.01,-2.20*weightOffset,
3639  // 1.00*weightOffset));
3640  // coefficients.push_back(new Variable("bkg3_x0y3", -0.58920 * weightOffset, 0.01,-1.00*weightOffset,
3641  // 2.00*weightOffset));
3642  // PolynomialPdf* poly = new PolynomialPdf("bkg3Pdf", observables, coefficients, offsets, 3);
3643  PolynomialPdf *poly = new PolynomialPdf("bkg3Pdf", observables, coefficients, offsets, 2);
3644 
3645  // coefficients.push_back(new Variable("bkg3_x0y0", 0.10 * weightOffset, 0.01, 0.01 * weightOffset, 0.20 *
3646  // weightOffset));
3647  // PolynomialPdf* poly = new PolynomialPdf("bkg3Pdf", observables, coefficients, offsets, 0);
3648 
3649  // Background 3 does not have a trig threshold in Mikhail's fit
3650  // - that is, it has one, but the dec variable is set above the
3651  // threshold, so it doesn't do anything.
3652  // That's Mikhail's fit; I'm putting one in to try to deal with
3653  // the asymmetry in the rho+.
3654  // Didn't work.
3655  // Variable* bkg3_decZmin = new Variable("bkg3_decZmin", 3.30411, 0.1, 1, 5);
3656  // Variable* bkg3_conZmin = new Variable("bkg3_conZmin", 0.29909, 0.01, 0.1, 0.9);
3657  // TrigThresholdPdf* bkg3_loZ = new TrigThresholdPdf("bkg3_loZ", m12, m13, minDalitzZ, bkg3_decZmin, bkg3_conZmin,
3658  // massSum, false);
3659 
3660  comps.clear();
3661  comps.push_back(poly);
3662  comps.push_back(kzero_veto);
3663  // comps.push_back(bkg3_eff);
3664  // comps.push_back(bkg3_loZ);
3665 
3666  ProdPdf *poly_x_veto = new ProdPdf("poly_x_veto", comps);
3667 
3668  // One misIDpi0.
3669  // Factor 3 in amplitudes is to offset division by total weight in AddPdf.
3670  DecayInfo3 special_pi0_decay;
3671  special_pi0_decay.motherMass = _mD0;
3672  special_pi0_decay.daug1Mass = piZeroMass;
3673  special_pi0_decay.daug2Mass = piPlusMass;
3674  special_pi0_decay.daug3Mass = piPlusMass;
3675  special_pi0_decay.meson_radius = 1.5;
3676 
3677  ResonancePdf *bkg3_pi0_ref = new Resonances::Gauss(
3678  "bkg3_pi0_ref",
3679  Variable("bkg3_pi0_ref_amp", 0.01189 * weightOffset, 0.01, 0.00 * weightOffset, 0.25 * weightOffset),
3680  constantZero,
3681  Variable("bkg3_pi0_ref_mass", 1.65766, 0.01, 1.4, 1.8),
3682  Variable("bkg3_pi0_ref_width", 0.05018, 0.01, 0.02, 0.20),
3683  PAIR_23);
3684  special_pi0_decay.resonances.push_back(bkg3_pi0_ref);
3685 
3686  // Mikhail defines 'transverse Z' as y - x - (parameter).
3687  // Variable* bkg3_pi0_transZ_offset = new Variable("bkg3_pi0_transZ_offset", -0.04381, 0.001, -0.5, 0.5);
3688  Variable bkg3_pi0_transZ_offset("bkg3_pi0_transZ_offset", -0.04381);
3689  offsets.clear();
3690  observables.clear();
3691  coefficients.clear();
3692  observables.push_back(*m12);
3693  observables.push_back(*m13);
3694  coefficients.push_back(bkg3_pi0_transZ_offset);
3695  coefficients.push_back(constantMinusOne);
3696  coefficients.push_back(constantOne);
3697  offsets.push_back(constantZero);
3698  offsets.push_back(constantZero);
3699  PolynomialPdf *bkg3_pi0_transZ = new PolynomialPdf("bkg3_pi0_transZ", observables, coefficients, offsets, 1);
3700 
3701  // Now we're going to take (1 - tz^2 * (parameter)) and multiply that into the misID pi0.
3702  // Variable* bkg3_pi0_transZ_quad = new Variable("bkg3_pi0_transZ_quad", 2.12277, 0.01, -1.5, 6.0);
3703  Variable bkg3_pi0_transZ_quad("bkg3_pi0_transZ_quad", 2.12277);
3704  coefficients.clear();
3705  coefficients.push_back(constantOne);
3706  coefficients.push_back(constantZero);
3707  coefficients.push_back(bkg3_pi0_transZ_quad);
3708  // Notice the fake dependence of the polynomial on m12; in fact CompositePdf
3709  // will send it a fake event, we just have to supply a reasonable index.
3710  PolynomialPdf *bkg3_pi0_shell = new PolynomialPdf("bkg3_pi0_shell", *m12, coefficients);
3711  CompositePdf *bkg3_pi0_transZ_total = new CompositePdf("bkg3_pi0_transZ_total", bkg3_pi0_transZ, bkg3_pi0_shell);
3712 
3713  comps.clear();
3714  comps.push_back(kzero_veto);
3715  comps.push_back(bkg3_pi0_transZ_total);
3716  // comps.push_back(bkg3_eff);
3717  // comps.push_back(bkg3_loZ);
3718  // ProdPdf* bkg3_pi0_mods = new ProdPdf("bkg3_pi0_mods", comps);
3719  // incsum3 = new IncoherentSumPdf("incsum3", m12, m13, eventNumber, special_pi0_decay, bkg3_pi0_mods);
3720 
3721  // Three spin-1 rho resonances to be added incoherently.
3722  DecayInfo3 incoherent_rhos;
3723  incoherent_rhos.motherMass = _mD0;
3724  incoherent_rhos.daug1Mass = piZeroMass;
3725  incoherent_rhos.daug2Mass = piPlusMass;
3726  incoherent_rhos.daug3Mass = piPlusMass;
3727  incoherent_rhos.meson_radius = 0; // Mikhail uses zero radius for incoherent resonances.
3728 
3729  ResonancePdf *bkg3_incRho0 = new Resonances::RBW(
3730  "bkg3_incRho0",
3731  Variable("bkg3_incRho0_amp", 0.00807 * weightOffset, 0.01, 0.00 * weightOffset, 0.25 * weightOffset),
3732  constantZero,
3733  Variable("bkg3_incRho0_mass", 0.800, 0.01, 0.6, 1.0),
3734  Variable("bkg3_incRho0_width", 0.15, 0.01, 0.10, 0.40),
3735  1, // These rhos are spin 1, being bad signal.
3736  PAIR_23);
3737  incoherent_rhos.resonances.push_back(bkg3_incRho0);
3738 
3739  ResonancePdf *bkg3_incRhoP = new Resonances::RBW(
3740  "bkg3_incRhoP",
3741  Variable("bkg3_incRhoP_amp", 0.01683 * weightOffset, 0.01, 0.00 * weightOffset, 0.25 * weightOffset),
3742  constantZero,
3743  Variable("bkg3_incRhoP_mass", 0.800, 0.01, 0.6, 1.0),
3744  Variable("bkg3_incRhoP_width", 0.15, 0.01, 0.10, 0.40),
3745  1,
3746  PAIR_12);
3747  incoherent_rhos.resonances.push_back(bkg3_incRhoP);
3748 
3749  ResonancePdf *bkg3_incRhoM = new Resonances::RBW(
3750  "bkg3_incRhoM",
3751  Variable("bkg3_incRhoM_amp", 0.01645 * weightOffset, 0.01, 0.00 * weightOffset, 0.25 * weightOffset),
3752  constantZero,
3753  Variable("bkg3_incRhoM_mass", 0.900, 0.01, 0.6, 1.0),
3754  Variable("bkg3_incRhoM_width", 0.35, 0.01, 0.10, 0.60),
3755  1,
3756  PAIR_13);
3757  incoherent_rhos.resonances.push_back(bkg3_incRhoM);
3758 
3759  comps.clear();
3760  comps.push_back(kzero_veto);
3761  // comps.push_back(bkg3_loZ);
3762  // comps.push_back(bkg3_eff);
3763  // ProdPdf* bkg3_rho_mods = new ProdPdf("bkg3_rho_mods", comps);
3764 
3765  // incsum4 = new IncoherentSumPdf("incsum4", m12, m13, eventNumber, incoherent_rhos, bkg3_rho_mods);
3766  // incsum4 = new IncoherentSumPdf("incsum4", m12, m13, eventNumber, incoherent_rhos, kzero_veto);
3767 
3768  weights.clear();
3769  weights.push_back(constantOne);
3770  // weights.push_back(constantOne);
3771  // weights.push_back(constantOne);
3772  comps.clear();
3773  comps.push_back(poly_x_veto);
3774  // comps.push_back(incsum3);
3775  // comps.push_back(incsum4);
3776 
3777  AddPdf *bkg3_dalitz = new AddPdf("bkg3_dalitz", weights, comps);
3778  bkg3_dalitz->addSpecialMask(PdfBase::ForceCommonNorm);
3779  return bkg3_dalitz;
3780 }
3781 
3783  vector<Variable> offsets;
3784  vector<Observable> observables;
3785  vector<Variable> coefficients;
3786  offsets.push_back(constantOne);
3787  offsets.push_back(constantOne);
3788 
3789  observables.push_back(*m12);
3790  observables.push_back(*m13);
3791  // Recurring factor 3 offsets division by total weight in AddPdf.
3792 
3793  double weightOffset = 3;
3794  // coefficients.push_back(new Variable("bkg4_x0y0", 1.0 * weightOffset, 0.01, 0.50*weightOffset,
3795  // 1.50*weightOffset));
3796  coefficients.push_back(Variable("bkg4_x0y0", 1.0 * weightOffset));
3797  coefficients.push_back(
3798  Variable("bkg4_x1y0", -0.18594 * weightOffset, 0.01, -0.50 * weightOffset, 0.50 * weightOffset));
3799  coefficients.push_back(
3800  Variable("bkg4_x2y0", 0.45459 * weightOffset, 0.01, 0.25 * weightOffset, 0.75 * weightOffset));
3801  coefficients.push_back(
3802  Variable("bkg4_x3y0", -0.20869 * weightOffset, 0.01, -0.50 * weightOffset, 0.50 * weightOffset));
3803  coefficients.push_back(
3804  Variable("bkg4_x0y1", -0.65061 * weightOffset, 0.01, -1.50 * weightOffset, 0.50 * weightOffset));
3805  coefficients.push_back(
3806  Variable("bkg4_x1y1", 0.11000 * weightOffset, 0.01, 0.00 * weightOffset, 0.50 * weightOffset));
3807  coefficients.push_back(
3808  Variable("bkg4_x2y1", 0.42009 * weightOffset, 0.01, 0.25 * weightOffset, 1.00 * weightOffset));
3809  coefficients.push_back(
3810  Variable("bkg4_x0y2", -0.06151 * weightOffset, 0.01, -0.50 * weightOffset, 0.50 * weightOffset));
3811  coefficients.push_back(
3812  Variable("bkg4_x1y2", 0.58508 * weightOffset, 0.01, 0.20 * weightOffset, 1.50 * weightOffset));
3813  coefficients.push_back(
3814  Variable("bkg4_x0y3", 0.54740 * weightOffset, 0.01, 0.20 * weightOffset, 1.50 * weightOffset));
3815 
3816  PolynomialPdf *poly = new PolynomialPdf("bkg4Pdf", observables, coefficients, offsets, 3);
3817 
3818  Variable bkg4_decZmin("bkg4_decZmin", 2.77576);
3819  Variable bkg4_conZmin("bkg4_conZmin", 0.23328);
3820  TrigThresholdPdf *bkg4_loZ
3821  = new TrigThresholdPdf("bkg4_loZ", *m12, *m13, minDalitzZ, bkg4_decZmin, bkg4_conZmin, massSum, false);
3822 
3823  comps.clear();
3824  comps.push_back(poly);
3825  comps.push_back(bkg4_loZ);
3826  comps.push_back(kzero_veto);
3827  // Separate PDF to avoid triggering numerical normalisation over all four observables.
3828  ProdPdf *poly_x_veto = new ProdPdf("poly_x_veto", comps);
3829 
3830  // One pipi bump.
3831  // Factor 3 in amplitudes is to offset division by total weight in AddPdf.
3832  DecayInfo3 special_pipi_decay;
3833  special_pipi_decay.motherMass = _mD0;
3834  special_pipi_decay.daug1Mass = piZeroMass;
3835  special_pipi_decay.daug2Mass = piPlusMass;
3836  special_pipi_decay.daug3Mass = piPlusMass;
3837  special_pipi_decay.meson_radius = 1.5;
3838 
3839  ResonancePdf *bkg4_pipi_ref = new Resonances::Gauss("bkg4_pipi_ref",
3840  Variable("bkg4_pipi_ref_amp", 0.00147 * weightOffset),
3841  constantZero,
3842  Variable("bkg4_pipi_ref_mass", 1.32447),
3843  Variable("bkg4_pipi_ref_width", 0.04675),
3844  PAIR_23);
3845  special_pipi_decay.resonances.push_back(bkg4_pipi_ref);
3846 
3847  // Mikhail defines 'transverse Z' as y - x - (parameter).
3848  Variable bkg4_pipi_transZ_offset("bkg4_pipi_transZ_offset", -0.39877);
3849 
3850  offsets.clear();
3851  observables.clear();
3852  coefficients.clear();
3853  observables.push_back(*m12);
3854  observables.push_back(*m13);
3855  coefficients.push_back(bkg4_pipi_transZ_offset);
3856  coefficients.push_back(constantMinusOne);
3857  coefficients.push_back(constantOne);
3858  offsets.push_back(constantZero);
3859  offsets.push_back(constantZero);
3860  PolynomialPdf *bkg4_pipi_transZ = new PolynomialPdf("bkg4_pipi_transZ", observables, coefficients, offsets, 1);
3861 
3862  // Now we're going to take (1 - tz^2 * (parameter)) and multiply that into the pipi bump.
3863  Variable bkg4_pipi_transZ_quad("bkg4_pipi_transZ_quad", -0.25640);
3864  coefficients.clear();
3865  coefficients.push_back(constantOne);
3866  coefficients.push_back(constantZero);
3867  coefficients.push_back(bkg4_pipi_transZ_quad);
3868  // Notice the fake dependence of the polynomial on m12; in fact CompositePdf
3869  // will send it a fake event, we just have to supply a reasonable index.
3870  PolynomialPdf *bkg4_pipi_shell = new PolynomialPdf("bkg4_pipi_shell", *m12, coefficients);
3871  CompositePdf *bkg4_pipi_transZ_total
3872  = new CompositePdf("bkg4_pipi_transZ_total", bkg4_pipi_transZ, bkg4_pipi_shell);
3873 
3874  comps.clear();
3875  comps.push_back(kzero_veto);
3876  comps.push_back(bkg4_loZ);
3877  comps.push_back(bkg4_pipi_transZ_total);
3878 
3879  ProdPdf *bkg4_pipi_mods = new ProdPdf("bkg4_pipi_mods", comps);
3880  incsum5 = new IncoherentSumPdf("incsum5", *m12, *m13, *eventNumber, special_pipi_decay, bkg4_pipi_mods);
3881 
3882  // Three spin-0 rho resonances to be added incoherently.
3883  DecayInfo3 incoherent_rho0s;
3884  incoherent_rho0s.motherMass = _mD0;
3885  incoherent_rho0s.daug1Mass = piZeroMass;
3886  incoherent_rho0s.daug2Mass = piPlusMass;
3887  incoherent_rho0s.daug3Mass = piPlusMass;
3888  incoherent_rho0s.meson_radius = 0; // Mikhail uses zero radius for incoherent resonances.
3889 
3890  ResonancePdf *bkg4_incRho0 = new Resonances::RBW("bkg4_incRho0",
3891  Variable("bkg4_incRho0_amp", 0.00429 * weightOffset),
3892  constantZero,
3893  fixedRhoMass,
3894  fixedRhoWidth,
3895  0, // These rhos are spin 0.
3896  PAIR_23);
3897  incoherent_rho0s.resonances.push_back(bkg4_incRho0);
3898 
3899  ResonancePdf *bkg4_incRhoP = new Resonances::RBW("bkg4_incRhoP",
3900  Variable("bkg4_incRhoP_amp", 0.00705 * weightOffset),
3901  constantZero,
3902  fixedRhoMass,
3903  fixedRhoWidth,
3904  0,
3905  PAIR_12);
3906  incoherent_rho0s.resonances.push_back(bkg4_incRhoP);
3907 
3908  ResonancePdf *bkg4_incRhoM = new Resonances::RBW("bkg4_incRhoM",
3909  Variable("bkg4_incRhoM_amp", -0.00043 * weightOffset),
3910  constantZero,
3911  fixedRhoMass,
3912  fixedRhoWidth,
3913  0,
3914  PAIR_13);
3915  incoherent_rho0s.resonances.push_back(bkg4_incRhoM);
3916 
3917  comps.clear();
3918  comps.push_back(kzero_veto);
3919  comps.push_back(bkg4_loZ);
3920  ProdPdf *bkg4_incrho_mods = new ProdPdf("bkg4_incrho_mods", comps);
3921  incsum6 = new IncoherentSumPdf("incsum6", *m12, *m13, *eventNumber, incoherent_rho0s, bkg4_incrho_mods);
3922  // incsum6 = new IncoherentSumPdf("incsum6", m12, m13, eventNumber, incoherent_rho0s, kzero_veto);
3923 
3924  weights.clear();
3925  weights.push_back(constantOne);
3926  weights.push_back(constantOne);
3927  weights.push_back(constantOne);
3928  comps.clear();
3929  comps.push_back(poly_x_veto);
3930  comps.push_back(incsum5);
3931  comps.push_back(incsum6);
3932 
3933  AddPdf *bkg4_dalitz = new AddPdf("bkg4_dalitz", weights, comps);
3934  bkg4_dalitz->addSpecialMask(PdfBase::ForceCommonNorm);
3935  return bkg4_dalitz;
3936 }
3937 
3938 GooPdf *makeBkg3DalitzPdf(bool fixem = true) {
3939  if(!kzero_veto)
3940  makeKzeroVeto();
3941 
3942  comps.clear();
3943  weights.clear();
3944 
3945  // Using a histogram for Dalitz description. Notice: The purpose of this
3946  // fit is to get a description for use in fitting actual data. When fitting
3947  // data I can just use the background 3 Monte Carlo histogram. When fitting
3948  // the MC, it doesn't matter what I do, because I'm going to be using the
3949  // histogram. So I load up all the MC data and use it in the histogram, either
3950  // way.
3951 
3952  GooPdf *bkg3_dalitz = nullptr;
3953 
3955  bkg3_dalitz = makeBackgroundHistogram(3);
3956  else
3957  bkg3_dalitz = makeBackground3DalitzParam();
3958 
3959  // bkg3_dalitz->setDebugMask(1);
3960 
3961  GooPdf *bkg3_dtime = nullptr;
3962 
3963  if(gaussBkgTime)
3964  bkg3_dtime = makeGaussianTimePdf(3);
3965  else
3966  bkg3_dtime = makeExpGausTimePdf(3);
3967 
3968  // bkg3_jsugg = makeBkg3_sigma();
3969  // bkg3_jsugg = sig0_jsugg; // Mikhail uses same sigma distribution for everything.
3970  bkg3_jsugg = makeBkg_sigma_strips(3);
3971  bkg3_jsugg->addSpecialMask(PdfBase::ForceSeparateNorm);
3972  // Otherwise ProdPdf tries to use the default overall integration,
3973  // because bkg3_jsugg depends on m12, m13 due to the striping, and that has
3974  // disastrous results for bkg3_dalitz. Note that this doesn't, actually,
3975  // contradict the ForceCommonNorm above, even though it looks like it should,
3976  // because CommonNorm applies to the AddPdf while SeparateNorm
3977  // applies to the ProdPdf.
3978 
3979  comps.clear();
3980  comps.push_back(bkg3_dalitz);
3981  // bkg3_dalitz->setDebugMask(1);
3982  // incsum3->setDebugMask(1);
3983  comps.push_back(bkg3_dtime);
3984  // bkg3_dtime->setDebugMask(1);
3985  // comps.push_back(bkg3_jsugg);
3986  // sig0_jsugg->setDebugMask(1);
3987 
3988  GooPdf *ret = new ProdPdf("bkg3_total", comps);
3989 
3990  if(fixem)
3991  ret->setParameterConstantness(true);
3992 
3993  return ret;
3994 }
3995 
3996 GooPdf *makeBkg4DalitzPdf(bool fixem = true) {
3997  if(!kzero_veto)
3998  makeKzeroVeto();
3999 
4000  comps.clear();
4001  weights.clear();
4002 
4003  GooPdf *bkg4_dalitz = nullptr;
4004 
4006  bkg4_dalitz = makeBackgroundHistogram(4);
4007  else
4008  bkg4_dalitz = makeBackground4DalitzParam();
4009 
4010  GooPdf *bkg4_dtime = nullptr;
4011 
4012  if(gaussBkgTime)
4013  bkg4_dtime = makeGaussianTimePdf(4);
4014  else
4015  bkg4_dtime = makeExpGausTimePdf(4);
4016 
4017  // Separate sigma_t distribution
4018  // bkg4_jsugg = makeBkg4_sigma();
4019  // bkg4_jsugg = sig0_jsugg; // Mikhail uses same sigma distribution for everything.
4020  bkg4_jsugg = makeBkg_sigma_strips(4);
4021  bkg4_jsugg->addSpecialMask(PdfBase::ForceSeparateNorm); // See comments to bkg3_jsugg.
4022 
4023  comps.clear();
4024  comps.push_back(bkg4_dalitz);
4025  comps.push_back(bkg4_dtime);
4026  // comps.push_back(bkg4_jsugg);
4027 
4028  ProdPdf *ret = new ProdPdf("bkg4_total", comps);
4029 
4030  if(fixem)
4031  ret->setParameterConstantness(true);
4032 
4033  return ret;
4034 }
4035 
4036 int runCanonicalFit(std::string fname, bool noPlots = true) {
4038 
4039  if(mdslices > 1)
4040  massd0 = new Observable(
4041  "massd0", 1.8654 + 0.0075 * md0_lower_window + md0offset, 1.8654 + 0.0075 * md0_upper_window + md0offset);
4042 
4043  std::cout << "Loading events from " << fname << std::endl;
4044  loadDataFile(fname);
4045 
4046  std::cout << "Creating overall signal PDF\n";
4047  GooPdf *overallSignal = makeOverallSignal();
4048 
4049  TRandom donram(blindSeed); // The rain and the sun!
4050 
4051  if(0 != blindSeed) {
4052  dtop0pp._xmixing.setBlind(donram.Gaus(0, 0.005));
4053  dtop0pp._ymixing.setBlind(donram.Gaus(0, 0.005));
4054  }
4055 
4056  // overallSignal->setDebugMask(1);
4057 
4058  int oldBins1 = m12->getNumBins();
4059  int oldBins2 = m13->getNumBins();
4060  // Too fine a binning here leads to bad results due to fluctuations.
4061  m12->setNumBins(bkgHistBins);
4062  m13->setNumBins(bkgHistBins);
4063  std::cout << "Creating background PDFs\n";
4064  GooPdf *bkg2Pdf = makeBkg2DalitzPdf();
4065  GooPdf *bkg3Pdf = makeBkg3DalitzPdf();
4066  GooPdf *bkg4Pdf = makeBkg4DalitzPdf();
4067  m12->setNumBins(oldBins1);
4068  m13->setNumBins(oldBins2);
4069 
4070  getBackgroundFile(2);
4071  std::cout << "Reading bkg2 parameters from " << strbuffer << std::endl;
4072  readFromFile(bkg2Pdf, strbuffer);
4073  // NB, background 3 and 4 params do not actually work. Only hists.
4074  getBackgroundFile(3);
4075  std::cout << "Reading bkg3 parameters from " << strbuffer << std::endl;
4076  readFromFile(bkg3Pdf, strbuffer);
4077  getBackgroundFile(4);
4078  std::cout << "Reading bkg4 parameters from " << strbuffer << std::endl;
4079  readFromFile(bkg4Pdf, strbuffer);
4080 
4081  bkg2Pdf->setParameterConstantness(true);
4082  bkg3Pdf->setParameterConstantness(true);
4083  bkg4Pdf->setParameterConstantness(true);
4084 
4085  // bkg3Pdf->setDebugMask(1);
4086 
4087  int eventSize = massd0 ? 11 : 10;
4088  std::cout << "Setting data size " << eventSize << std::endl;
4089  signalDalitz->setDataSize(data->getNumEvents(), eventSize); // Must take into account event weights!
4090 
4091  // bkg2Pdf->setDebugMask(1);
4092  if(incsum1)
4093  incsum1->setDataSize(data->getNumEvents(), eventSize);
4094 
4095  if(incsum2)
4096  incsum2->setDataSize(data->getNumEvents(), eventSize);
4097 
4098  // bkg3Pdf->setDebugMask(1);
4099  if(incsum3)
4100  incsum3->setDataSize(data->getNumEvents(), eventSize);
4101 
4102  if(incsum4)
4103  incsum4->setDataSize(data->getNumEvents(), eventSize);
4104 
4105  if(incsum5)
4106  incsum5->setDataSize(data->getNumEvents(), eventSize);
4107 
4108  if(incsum6)
4109  incsum6->setDataSize(data->getNumEvents(), eventSize);
4110 
4111  std::cout << "Creating overall PDF\n";
4112  std::vector<Observable> evtWeights;
4113  evtWeights.push_back(*wSig0);
4114  evtWeights.push_back(*wBkg2);
4115  evtWeights.push_back(*wBkg3);
4116  evtWeights.push_back(*wBkg4);
4117  std::vector<PdfBase *> components;
4118  components.push_back(overallSignal);
4119  components.push_back(bkg2Pdf);
4120  components.push_back(bkg3Pdf);
4121  components.push_back(bkg4Pdf);
4122  EventWeightedAddPdf *overallPdf = new EventWeightedAddPdf("total", evtWeights, components);
4123  // overallPdf->setDebugMask(1);
4124  std::cout << "Copying data to GPU\n";
4125  overallPdf->setData(data);
4126 
4127  if(paramUp != "") {
4128  Variable *target = overallPdf->getParameterByName(paramUp);
4129  assert(target);
4130  target->setValue(target->getValue() + target->getError());
4131  }
4132 
4133  if(paramDn != "") {
4134  Variable *target = overallPdf->getParameterByName(paramDn);
4135  assert(target);
4136  target->setValue(target->getValue() - target->getError());
4137  }
4138 
4139  int retval;
4140  if(minuit1) {
4141  GooFit::FitManagerMinuit1 datapdf(overallPdf);
4142  datapdf.setMaxCalls(64000);
4143  datapdf.fit();
4144  retval = datapdf;
4145  } else {
4146  GooFit::FitManagerMinuit2 datapdf(overallPdf);
4147  datapdf.setMaxCalls(64000);
4148  datapdf.fit();
4149  retval = datapdf;
4150  }
4151 
4152 #ifdef PROFILING
4153  overallPdf->printProfileInfo();
4154 #endif
4155 
4156  fmt::print("Fit results Canonical fit:\n"
4157  "tau : ({:.3}) fs\n"
4158  "xmixing: ({:.3})\%\n"
4159  "ymixing: ({:.3})\%\n",
4160  1000 * Uncertain(dtop0pp._tau),
4161  100 * Uncertain(dtop0pp._xmixing),
4162  100 * Uncertain(dtop0pp._ymixing));
4163 
4164  /*
4165  std::cout << "Fit results: \n"
4166  << "tau : (" << 1000*dtop0pp._tau.getValue() << " $\\pm$ " << 1000*dtop0pp._tau.getError() << ") fs\n"
4167  << "xmixing: (" << 100*dtop0pp._xmixing.getValue() << " $\\pm$ " << 100*dtop0pp._xmixing.getError() << ")%\n"
4168  << "ymixing: (" << 100*dtop0pp._ymixing.getValue() << " $\\pm$ " << 100*dtop0pp._ymixing.getError() << ")%\n";
4169  */
4170 
4171  /*
4172  double fitx = dtop0pp._xmixing.getValue();
4173  TH1F xscan("xscan", "", 200, fitx - 0.0001 - 0.0000005, fitx + 0.0001 - 0.0000005);
4174  xscan.SetStats(false);
4175  //double fity = dtop0pp._ymixing.getValue();
4176  for (int i = 0; i < 200; ++i) {
4177  dtop0pp._xmixing.getValue() = fitx - 0.0001 + 0.000001*i;
4178  overallPdf->copyParams();
4179  double nll = overallPdf->calculateNLL();
4180  printf("%i: %.10f\n", i, nll);
4181  xscan.SetBinContent(i+1, nll - floor(nll));
4182  }
4183 
4184  foo->cd();
4185  xscan.Draw();
4186  foo->SaveAs("xscan.png");
4187  */
4188  if(!noPlots)
4189  makeDalitzPlots(overallSignal);
4190 
4191  return retval;
4192 }
4193 
4194 int runSigmaFit(const char *fname) {
4196 
4197  loM23Sigma = new TH1F("loM23Sigma", "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
4198  loM23Sigma->SetStats(false);
4199  hiM23Sigma = new TH1F("hiM23Sigma", "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
4200  hiM23Sigma->SetStats(false);
4201 
4202  loadDataFile(fname);
4203  // GooPdf* jsu_gg = makeSignalJSU_gg(-1, false);
4204  // GooPdf* jsu_gg = makeSigmaMap();
4205  GooPdf *jsu_gg = makeBkg_sigma_strips(0);
4206  jsu_gg->setData(data);
4207  // jsu_gg->copyParams();
4208 
4209  int retval;
4210  if(minuit1) {
4211  GooFit::FitManagerMinuit1 datapdf(jsu_gg);
4212  datapdf.fit();
4213  retval = datapdf;
4214  } else {
4215  GooFit::FitManagerMinuit2 datapdf(jsu_gg);
4216  datapdf.fit();
4217  retval = datapdf;
4218  }
4219 
4220  sprintf(strbuffer, "signal_sigma_%islices_pdf.txt", m23Slices);
4221  writeToFile(jsu_gg, strbuffer);
4222 
4223  foo->cd();
4224  plotLoHiSigma();
4225 
4226  std::vector<Observable> gridvars;
4227  gridvars.push_back(*m12);
4228  gridvars.push_back(*m13);
4229  gridvars.push_back(*sigma);
4230  UnbinnedDataSet grid(gridvars);
4231 
4232  TH1F *sigma_pdfs[6];
4233  TH1F *sigma_data[6];
4234  double num_sigma_dat[6];
4235  double num_sigma_pdf[6];
4236 
4237  for(int i = 0; i < 6; ++i) {
4238  sprintf(strbuffer, "sigma_pdf_%i", i);
4239  sigma_pdfs[i] = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
4240  sprintf(strbuffer, "sigma_dat_%i", i);
4241  sigma_data[i] = new TH1F(strbuffer, "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
4242 
4243  num_sigma_dat[i] = 0;
4244  num_sigma_pdf[i] = 0;
4245 
4246  sigma_data[i]->SetStats(false);
4247  sigma_data[i]->SetMarkerStyle(8);
4248  sigma_data[i]->SetMarkerSize(1.2);
4249  sigma_data[i]->GetXaxis()->SetTitle("Decay time error [ps]");
4250  sigma_data[i]->GetYaxis()->SetTitle("Events / 8 fs");
4251 
4252  sigma_pdfs[i]->SetStats(false);
4253  sigma_pdfs[i]->SetLineColor(kBlue);
4254  sigma_pdfs[i]->SetLineWidth(3);
4255  }
4256 
4257  double totalPdf = 0;
4258  double totalDat = 0;
4259 
4260  for(unsigned int evt = 0; evt < data->getNumEvents(); ++evt) {
4261  double currSigma = data->getValue(*sigma, evt);
4262  double currm12 = data->getValue(*m12, evt);
4263  double currm13 = data->getValue(*m13, evt);
4264  double currm23 = cpuGetM23(currm12, currm13);
4265  int m23bin = (int)floor(currm23 / 0.5);
4266  sigma_data[m23bin]->Fill(currSigma);
4267  num_sigma_dat[m23bin]++;
4268  totalDat++;
4269  }
4270 
4271  for(int i = 0; i < m12->getNumBins(); ++i) {
4272  m12->setValue(m12->getLowerLimit() + (i + 0.5) * m12->getBinSize());
4273 
4274  for(int j = 0; j < m13->getNumBins(); ++j) {
4275  m13->setValue(m13->getLowerLimit() + (j + 0.5) * m13->getBinSize());
4276 
4278  continue;
4279 
4280  for(int k = 0; k < sigma->getNumBins(); ++k) {
4281  sigma->setValue(sigma->getLowerLimit() + (k + 0.5) * sigma->getBinSize());
4282  grid.addEvent();
4283  }
4284  }
4285  }
4286 
4287  jsu_gg->setData(&grid);
4288 
4289  std::vector<std::vector<double>> pdfValues = jsu_gg->getCompProbsAtDataPoints();
4290 
4291  for(unsigned int j = 0; j < pdfValues[0].size(); ++j) {
4292  double currM12 = grid.getValue(*m12, j);
4293  double currM13 = grid.getValue(*m13, j);
4294  double currSigma = grid.getValue(*sigma, j);
4295  double currm23 = cpuGetM23(currM12, currM13);
4296  int m23bin = (int)floor(currm23 / 0.5);
4297  sigma_pdfs[m23bin]->Fill(currSigma, pdfValues[0][j]);
4298  num_sigma_pdf[m23bin] += pdfValues[0][j];
4299  totalPdf += pdfValues[0][j];
4300  }
4301 
4302  for(int i = 1; i <= sigma->getNumBins(); ++i) {
4303  for(int j = 0; j < 6; ++j) {
4304  sigma_pdfs[j]->SetBinContent(i, sigma_pdfs[j]->GetBinContent(i) * num_sigma_dat[j] / num_sigma_pdf[j]);
4305  }
4306  }
4307 
4308  std::string plotdir = "./plots_from_mixfit/";
4309 
4310  for(int i = 0; i < 6; ++i) {
4311  if(sigma_data[i]->GetMaximum() > sigma_pdfs[i]->GetMaximum()) {
4312  sigma_data[i]->Draw("p");
4313  sigma_pdfs[i]->Draw("lsame");
4314  } else {
4315  sigma_pdfs[i]->Draw("l");
4316  sigma_data[i]->Draw("psame");
4317  }
4318 
4319  sprintf(strbuffer, "%i", i);
4320  TText slicenum;
4321  slicenum.DrawTextNDC(0.2, 0.8, strbuffer);
4322 
4323  foo->SaveAs((plotdir + sigma_pdfs[i]->GetName() + ".png").c_str());
4324  foo->SetLogy(true);
4325  foo->SaveAs((plotdir + sigma_pdfs[i]->GetName() + "_log.png").c_str());
4326  foo->SetLogy(false);
4327  }
4328 
4329  /*
4330  // This code assumes you're using the PDF from makeSigmaMap.
4331  TCanvas foodal("", "", 600, 600);
4332  foodal.Divide(6, 6, 0, 0);
4333  for (int i = 0; i < numSigmaBins; ++i) {
4334  int xbin = i % 6;
4335  int ybin = i / 6;
4336 
4337  if (0 == sigma_data[i]->numEvents()) continue;
4338 
4339  //m12->getValue() = 0.5*(xbin + 0.5);
4340  //m13->getValue() = 0.5*(ybin + 0.5);
4341 
4342  std::vector<fptype> values;
4343  jsuList[i]->evaluateAtPoints(sigma, values);
4344 
4345  double totalPdf = 0;
4346  for (int bin = 0; bin < sigma->getNumBins(); ++bin) {
4347  totalPdf += values[bin];
4348  }
4349  for (int bin = 0; bin < sigma->getNumBins(); ++bin) {
4350  sigma_pdf_hists[i]->SetBinContent(bin+1, values[bin] * sigma_data[i]->numEvents() / totalPdf);
4351  }
4352 
4353  int padNumber = 31 - 6*ybin;
4354  padNumber += xbin;
4355  foodal.cd(padNumber);
4356  sigma_dat_hists[i]->Draw("p");
4357  sigma_pdf_hists[i]->Draw("lsame");
4358  }
4359  foodal.SaveAs("./plots_from_mixfit/sigma_dalitz.png");
4360  */
4361 
4362  return retval;
4363 }
4364 
4365 int runEfficiencyFit(int which) {
4367 
4368  if(3 == which) {
4369  m12->setNumBins(m12->getNumBins() / 8);
4370  m13->setNumBins(m13->getNumBins() / 8);
4371  }
4372 
4373  vector<Observable *> lvars;
4374  lvars.push_back(m12);
4375  lvars.push_back(m13);
4376  // binEffData = new BinnedDataSet(lvars);
4377  // GooPdf* eff = makeEfficiencyPdf();
4378 
4379  makeKzeroVeto();
4380  // GooPdf* eff = makeEfficiencyPoly();
4381  GooPdf *eff = makeEfficiencyPdf();
4382 
4383  std::string fname_3flat = app_ptr->get_filename("dataFiles/efficiency_bkg3_flat.txt", "examples/pipipi0DPFit");
4384  std::string fname_flat = app_ptr->get_filename("dataFiles/efficiency_flat.txt", "examples/pipipi0DPFit");
4385 
4386  if(3 == which)
4387  loadDataFile(fname_3flat);
4388  else
4389  loadDataFile(fname_flat);
4390 
4391  if(underlyingBins) {
4392  underlyingBins->GetZaxis()->SetRangeUser(10, 40);
4393  underlyingBins->Draw("colz");
4394  foo->SaveAs("./plots_from_mixfit/efficiency_bins.png");
4395  }
4396 
4397  // eff->setDebugMask(1);
4398  eff->setData(data);
4399 
4400  int retval;
4401  if(minuit1) {
4402  GooFit::FitManagerMinuit1 datapdf(eff);
4403  datapdf.fit();
4404  retval = datapdf;
4405  } else {
4406  GooFit::FitManagerMinuit2 datapdf(eff);
4407  datapdf.fit();
4408  retval = datapdf;
4409  }
4410 
4411  // plotFit(sigma, data, jsu_gg);
4412 
4413  TH2F dalitz_dat_hist("dalitz_dat_hist",
4414  "",
4415  m12->getNumBins(),
4416  m12->getLowerLimit(),
4417  m12->getUpperLimit(),
4418  m13->getNumBins(),
4419  m13->getLowerLimit(),
4420  m13->getUpperLimit());
4421  dalitz_dat_hist.SetStats(false);
4422  dalitz_dat_hist.GetXaxis()->SetTitle("m^{2}(#pi^{+} #pi^{0}) [GeV]");
4423  dalitz_dat_hist.GetYaxis()->SetTitle("m^{2}(#pi^{-} #pi^{0}) [GeV]");
4424  TH2F dalitz_pdf_hist("dalitz_pdf_hist",
4425  "",
4426  m12->getNumBins(),
4427  m12->getLowerLimit(),
4428  m12->getUpperLimit(),
4429  m13->getNumBins(),
4430  m13->getLowerLimit(),
4431  m13->getUpperLimit());
4432  dalitz_pdf_hist.SetStats(false);
4433 
4434  double totalPdf = 0;
4435  double totalDat = 0;
4436 
4437  for(unsigned int evt = 0; evt < data->getNumEvents(); ++evt) {
4438  double currval = data->getValue(*m12, evt);
4439  // m12_dat_hist.Fill(currval);
4440  double currval2 = data->getValue(*m13, evt);
4441  // m13_dat_hist.Fill(currval2);
4442  dalitz_dat_hist.Fill(currval, currval2);
4443  totalDat++;
4444  }
4445 
4446  std::vector<Observable> nvars;
4447  nvars.push_back(*m12);
4448  nvars.push_back(*m13);
4449  UnbinnedDataSet currData(nvars);
4450 
4451  for(int i = 0; i < m12->getNumBins(); ++i) {
4452  m12->setValue(m12->getLowerLimit() + (i + 0.5) * m12->getBinSize());
4453 
4454  for(int j = 0; j < m13->getNumBins(); ++j) {
4455  m13->setValue(m13->getLowerLimit() + (j + 0.5) * m13->getBinSize());
4456 
4458  continue;
4459 
4460  currData.addEvent();
4461  }
4462  }
4463 
4464  eff->setData(&currData);
4465 
4466  // eff->setDebugMask(1);
4467  std::vector<std::vector<double>> pdfValues = eff->getCompProbsAtDataPoints();
4468 
4469  for(unsigned int j = 0; j < pdfValues[0].size(); ++j) {
4470  double currVal = currData.getValue(*m12, j);
4471  // m12_pdf_hist.Fill(currVal, pdfValues[0][j]);
4472  double currVal2 = currData.getValue(*m13, j);
4473  // m13_pdf_hist.Fill(currVal, pdfValues[0][j]);
4474  dalitz_pdf_hist.Fill(currVal, currVal2, pdfValues[0][j]);
4475 
4476  totalPdf += pdfValues[0][j];
4477 
4478  if(std::isnan(pdfValues[0][j])) {
4479  std::cout << "Major problem: " << currVal << " " << currVal2 << " " << j << std::endl;
4480  assert(false);
4481  }
4482 
4483  if(std::isinf(pdfValues[0][j])) {
4484  std::cout << "Infinity " << j << std::endl;
4485  assert(false);
4486  }
4487  }
4488 
4489  /*for (int i = 1; i <= m12->getNumBins(); ++i) {
4490  m12_pdf_hist.SetBinContent(i, m12_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
4491  }
4492 
4493  for (int i = 1; i <= m13->getNumBins(); ++i) {
4494  m13_pdf_hist.SetBinContent(i, m13_pdf_hist.GetBinContent(i) * totalDat / totalPdf);
4495  }
4496  */
4497  for(int i = 1; i <= m12->getNumBins(); ++i) {
4498  for(int j = 1; j <= m13->getNumBins(); ++j) {
4499  double currNormEff = dalitz_pdf_hist.GetBinContent(i, j) * totalDat / totalPdf;
4500  dalitz_pdf_hist.SetBinContent(i, j, currNormEff);
4501 
4502  if(currNormEff >= 50)
4503  std::cout << "High efficiency: " << currNormEff << " " << i << " " << j << std::endl;
4504  }
4505  }
4506 
4507  foodal->cd();
4508 
4509  dalitz_dat_hist.GetZaxis()->SetRangeUser(0, 40);
4510  dalitz_dat_hist.Draw("colz");
4511  foodal->SaveAs("./plots_from_mixfit/efficiency_dat.png");
4512 
4513  dalitz_pdf_hist.GetZaxis()->SetRangeUser(0, 40);
4514  dalitz_pdf_hist.Draw("colz");
4515  foodal->SaveAs("./plots_from_mixfit/efficiency_pdf.png");
4516  foo->cd();
4517 
4518  TH1F pullplot("pullplot", "", 100, -5, 5);
4519  TH1F hiM23pullplot("hiM23pullplot", "", 100, -5, 5);
4520  TH1F loM23pullplot("loM23pullplot", "", 100, -5, 5);
4521 
4522  for(int i = 1; i <= m12->getNumBins(); ++i) {
4523  double m12loedge
4524  = m12->getLowerLimit() + ((m12->getUpperLimit() - m12->getLowerLimit()) / m12->getNumBins()) * (i - 1);
4525  double m12hiedge
4526  = m12->getLowerLimit() + ((m12->getUpperLimit() - m12->getLowerLimit()) / m12->getNumBins()) * (i);
4527 
4528  for(int j = 1; j <= m13->getNumBins(); ++j) {
4529  double m13loedge
4530  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j - 1);
4531 
4532  if(!cpuDalitz(m12loedge, m13loedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
4533  dalitz_dat_hist.SetBinContent(i, j, 0);
4534  continue;
4535  }
4536 
4537  if(!cpuDalitz(m12hiedge, m13loedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
4538  dalitz_dat_hist.SetBinContent(i, j, 0);
4539  continue;
4540  }
4541 
4542  double m13hiedge
4543  = m13->getLowerLimit() + ((m13->getUpperLimit() - m13->getLowerLimit()) / m13->getNumBins()) * (j);
4544 
4545  if(!cpuDalitz(m12loedge, m13hiedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
4546  dalitz_dat_hist.SetBinContent(i, j, 0);
4547  continue;
4548  }
4549 
4550  if(!cpuDalitz(m12hiedge, m13hiedge, _mD0, piZeroMass, piPlusMass, piPlusMass)) {
4551  dalitz_dat_hist.SetBinContent(i, j, 0);
4552  continue;
4553  }
4554 
4555  double dat = dalitz_dat_hist.GetBinContent(i, j);
4556  double pdf = dalitz_pdf_hist.GetBinContent(i, j);
4557 
4558  double pull = (dat - pdf) / sqrt(max(1.0, dat));
4559  // if (fabs(pull) > 5) continue;
4560  dalitz_dat_hist.SetBinContent(i, j, pull);
4561  pullplot.Fill(pull);
4562 
4563  double currm12 = dalitz_dat_hist.GetXaxis()->GetBinCenter(i);
4564  double currm13 = dalitz_dat_hist.GetYaxis()->GetBinCenter(j);
4565 
4566  // double currm23 = cpuGetM23(currm12, currm13);
4567  // if (currm23 > 1.5) hiM23pullplot.Fill(pull);
4568  // else loM23pullplot.Fill(pull);
4569  if((currm13 > 2) || (currm12 > 2))
4570  hiM23pullplot.Fill(pull);
4571  else
4572  loM23pullplot.Fill(pull);
4573  }
4574  }
4575 
4576  foodal->cd();
4577  dalitz_dat_hist.GetZaxis()->SetRangeUser(-5, 5);
4578  dalitz_dat_hist.Draw("colz");
4579  foodal->SaveAs("./plots_from_mixfit/efficiency_pull.png");
4580 
4581  foo->cd();
4582  pullplot.Draw();
4583  foo->SaveAs("./plots_from_mixfit/effpull.png");
4584 
4585  hiM23pullplot.Draw();
4586  foo->SaveAs("./plots_from_mixfit/hieffpull.png");
4587 
4588  loM23pullplot.Draw();
4589  foo->SaveAs("./plots_from_mixfit/loeffpull.png");
4590 
4591  return retval;
4592 }
4593 
4594 int runBackgroundDalitzFit(int bkgType, bool plots) {
4596  makeKzeroVeto();
4597 
4598  GooPdf *bkgPdf = 0;
4599 
4600  switch(bkgType) {
4601  default:
4602  case 2:
4603  bkgPdf = makeBkg2DalitzPdf(false);
4604  break;
4605 
4606  case 3:
4607  bkgPdf = makeBkg3DalitzPdf(false);
4608  break;
4609 
4610  case 4:
4611  bkgPdf = makeBkg4DalitzPdf(false);
4612  break;
4613  }
4614 
4615  sprintf(strbuffer, "./dataFiles/bkgDalitz_%i.txt", bkgType);
4616  std::string fname = app_ptr->get_filename(strbuffer, "examples/pipipi0DPFit");
4617  loadDataFile(fname);
4618 
4619  bkgPdf->setData(data);
4620  // bkgPdf->setDebugMask(1);
4621 
4622  // Incoherent-sum components need to know data size; check which ones exist.
4623  int eventSize = 5; // m12 m13 dtime sigma eventNumber
4624 
4625  if(incsum1)
4626  incsum1->setDataSize(data->getNumEvents(), eventSize);
4627 
4628  if(incsum2)
4629  incsum2->setDataSize(data->getNumEvents(), eventSize);
4630 
4631  if(incsum3)
4632  incsum3->setDataSize(data->getNumEvents(), eventSize);
4633 
4634  if(incsum4)
4635  incsum4->setDataSize(data->getNumEvents(), eventSize);
4636 
4637  if(incsum5)
4638  incsum5->setDataSize(data->getNumEvents(), eventSize);
4639 
4640  if(incsum6)
4641  incsum6->setDataSize(data->getNumEvents(), eventSize);
4642 
4643  int retval;
4644  if(minuit1) {
4645  GooFit::FitManagerMinuit1 fitter(bkgPdf);
4646  fitter.setMaxCalls(32000);
4647  fitter.fit();
4648  retval = fitter;
4649  } else {
4650  GooFit::FitManagerMinuit2 fitter(bkgPdf);
4651  fitter.setMaxCalls(32000);
4652  fitter.fit();
4653  retval = fitter;
4654  }
4655 
4656  if(plots) {
4657  // sprintf(strbuffer, "./plots_from_mixfit/bkgdalitz_%i/", bkgType);
4658  // makeDalitzPlots(bkgPdf, strbuffer);
4659  getBackgroundFile(bkgType);
4660  std::string fname = app_ptr->get_filename(strbuffer, "examples/pipipi0DPFit");
4661  writeToFile(bkgPdf, fname.c_str());
4662  }
4663 
4664  return retval;
4665 }
4666 
4667 void getBackgroundFile(int bkgType) {
4668  if(mikhailSetup)
4669  sprintf(strbuffer, "./bkg_%i_mikhail.txt", bkgType);
4670  else {
4671  if(2 == bkgType) {
4672  if(Sideband == bkg2Model)
4673  sprintf(strbuffer, "./bkg_2_pdf_sideband_%islices.txt", m23Slices);
4674  else
4675  sprintf(strbuffer, "./bkg_2_pdf_%islices.txt", m23Slices);
4676  } else {
4677  std::string pdftype;
4678 
4679  if(((3 == bkgType) && (notUseBackground3Hist)) || ((4 == bkgType) && (notUseBackground4Hist)))
4680  pdftype = "_param";
4681 
4682  sprintf(strbuffer, "./bkg_%i_pdf%s.txt", bkgType, pdftype.c_str());
4683  }
4684  }
4685 }
4686 
4687 void makeTimePlots(std::string fname) {
4689  massd0 = new Observable("massd0", 1.8654 + 0.0075 * md0_lower_window, 1.8654 + 0.0075 * md0_upper_window);
4690  massd0->setNumBins(180);
4691  std::cout << "Loading MC data from " << fname << std::endl;
4692  loadDataFile(fname);
4693 
4694  TH1F timeMean("timeMean", "", 6, massd0->getLowerLimit(), massd0->getUpperLimit());
4695  timeMean.SetStats(false);
4696  timeMean.SetLineWidth(3);
4697  timeMean.SetXTitle("#pi#pi#pi^{0} mass [GeV]");
4698  timeMean.SetYTitle("Mean of decay time [ps]");
4699  TH2F timeVsMass("timeVsMass",
4700  "",
4701  massd0->getNumBins(),
4702  massd0->getLowerLimit(),
4703  massd0->getUpperLimit(),
4704  dtime->getNumBins(),
4705  dtime->getLowerLimit(),
4706  dtime->getUpperLimit());
4707  timeVsMass.SetStats(false);
4708  timeVsMass.GetXaxis()->SetTitle("#pi#pi#pi^{0} mass [GeV]");
4709  timeVsMass.GetYaxis()->SetTitle("Decay time [ps]");
4710 
4711  int colors[6] = {kViolet + 1, kBlue, kCyan, kGreen, kYellow, kRed};
4712  TH1F *timePlots[6];
4713  TH1F *massPlots[5];
4714 
4715  for(int i = 0; i < 6; ++i) {
4716  sprintf(strbuffer, "timePlot_%i.png", i);
4717  timePlots[i] = new TH1F(strbuffer, "", dtime->getNumBins(), dtime->getLowerLimit(), dtime->getUpperLimit());
4718  timePlots[i]->SetStats(false);
4719  timePlots[i]->SetXTitle("Decay time [ps]");
4720  timePlots[i]->SetYTitle("Ratio");
4721  timePlots[i]->SetLineWidth(3);
4722  timePlots[i]->SetLineColor(colors[i]);
4723 
4724  if(i == 5)
4725  continue;
4726 
4727  sprintf(strbuffer, "massPlot_%i.png", i);
4728  massPlots[i] = new TH1F(strbuffer, "", massd0->getNumBins(), massd0->getLowerLimit(), massd0->getUpperLimit());
4729  massPlots[i]->SetStats(false);
4730  massPlots[i]->SetLineWidth(3);
4731  massPlots[i]->SetLineColor(colors[i]);
4732  }
4733 
4734  for(int i = 0; i < data->getNumEvents(); ++i) {
4735  data->loadEvent(i);
4736  timeVsMass.Fill(massd0->getValue(), dtime->getValue());
4737 
4738  if(massd0->getValue() >= massd0->getUpperLimit())
4739  continue;
4740 
4741  if(massd0->getValue() < massd0->getLowerLimit())
4742  continue;
4743 
4744  int slice = (int)floor(6 * (massd0->getValue() - massd0->getLowerLimit())
4745  / (massd0->getUpperLimit() - massd0->getLowerLimit()));
4746  timePlots[slice]->Fill(dtime->getValue());
4747 
4748  slice = (int)floor(5 * (dtime->getValue() - dtime->getLowerLimit())
4749  / (dtime->getUpperLimit() - dtime->getLowerLimit()));
4750  massPlots[slice]->Fill(massd0->getValue());
4751  }
4752 
4753  foo->cd();
4754  normalize(timePlots[3]);
4755  timePlots[3]->SetMinimum(0);
4756  timePlots[3]->Draw("hist");
4757 
4758  for(int i = 0; i < 6; ++i) {
4759  normalize(timePlots[i]);
4760  timePlots[i]->Draw("histsame");
4761  timeMean.SetBinContent(i + 1, timePlots[i]->GetMean());
4762  timeMean.SetBinError(i + 1, timePlots[i]->GetMeanError());
4763  }
4764 
4765  foo->SaveAs("timePlots.png");
4766  timeMean.Draw("e");
4767  foo->SaveAs("timeMeanPlot.png");
4768 
4769  // normalize(massPlots[2]);
4770  massPlots[2]->GetYaxis()->SetRangeUser(0, massPlots[2]->GetMaximum() * 1.1);
4771  massPlots[2]->Draw("");
4772 
4773  for(int i = 0; i < 5; ++i) {
4774  // normalize(massPlots[i]);
4775  massPlots[i]->Draw("same");
4776  }
4777 
4778  foo->SaveAs("massPlots.png");
4779 
4780  timeVsMass.Draw("colz");
4781 
4782  for(int i = 0; i < 6; ++i) {
4783  TLine *currLine
4784  = new TLine(massd0->getLowerLimit() + (i + 0) * (massd0->getUpperLimit() - massd0->getLowerLimit()) / 6,
4785  dtime->getLowerLimit() + 0.09,
4786  massd0->getLowerLimit() + (i + 1) * (massd0->getUpperLimit() - massd0->getLowerLimit()) / 6,
4787  dtime->getLowerLimit() + 0.09);
4788  currLine->SetLineWidth(12);
4789  currLine->SetLineColor(colors[i]);
4790  currLine->Draw();
4791 
4792  if(5 == i)
4793  continue;
4794 
4795  currLine = new TLine(massd0->getLowerLimit() + 0.00025,
4796  dtime->getLowerLimit() + (i + 0) * (dtime->getUpperLimit() - dtime->getLowerLimit()) / 5,
4797  massd0->getLowerLimit() + 0.00025,
4798  dtime->getLowerLimit() + (i + 1) * (dtime->getUpperLimit() - dtime->getLowerLimit()) / 5);
4799  currLine->SetLineWidth(12);
4800  currLine->SetLineColor(colors[i]);
4801  // currLine->Draw();
4802  }
4803 
4804  foo->SaveAs("timeVsMass.png");
4805 }
4806 
4807 int runBackgroundSigmaFit(int bkgType) {
4809 
4810  GooPdf *bkgPdf = 0;
4811 
4812  switch(bkgType) {
4813  default:
4814  case 2:
4815  bkgPdf = makeBkg2_sigma();
4816  break;
4817 
4818  case 3:
4819  loM23Sigma = new TH1F("loM23Sigma", "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
4820  loM23Sigma->SetStats(false);
4821  hiM23Sigma = new TH1F("hiM23Sigma", "", sigma->getNumBins(), sigma->getLowerLimit(), sigma->getUpperLimit());
4822  hiM23Sigma->SetStats(false);
4823 
4824  bkgPdf = makeBkg_sigma_strips(3);
4825  break;
4826  }
4827 
4828  sprintf(strbuffer, "./dataFiles/bkgDalitz_%i.txt", bkgType);
4829  std::string fname = app_ptr->get_filename(strbuffer, "examples/pipipi0DPFit");
4830  loadDataFile(fname);
4831  bkgPdf->setData(data);
4832 
4833  int retval;
4834 
4835  if(minuit1) {
4836  GooFit::FitManagerMinuit1 fitter(bkgPdf);
4837  fitter.setMaxCalls(8000);
4838  fitter.fit();
4839  retval = fitter;
4840  } else {
4841  GooFit::FitManagerMinuit2 fitter(bkgPdf);
4842  fitter.setMaxCalls(8000);
4843  fitter.fit();
4844  retval = fitter;
4845  }
4846 
4847  // bkgPdf->setDebugMask(1);
4848  plotFit(*sigma, data, bkgPdf);
4849  plotLoHiSigma();
4850 
4851  return retval;
4852 
4853  // sprintf(strbuffer, "./plots_from_mixfit/bkgdalitz_%i/", bkgType);
4854  // makeDalitzPlots(bkgPdf, strbuffer);
4855 }
4856 
4860  thrust::host_vector<fptype> host_hist;
4861  bkg3->extractHistogram(host_hist);
4862  sprintf(strbuffer, "Bkg%i_dalitzhist.txt", bkg);
4863  std::string fname = app_ptr->get_filename(strbuffer, "examples/pipipi0DPFit");
4864  writeListOfNumbers(host_hist, fname.c_str());
4865 }
4866 
4868  if(bkg2Model_str == "histogram")
4869  bkg2Model = Histogram;
4870  else if(bkg2Model_str == "parameter")
4871  bkg2Model = Parameter;
4872  else if(bkg2Model_str == "sideband")
4873  bkg2Model = Sideband;
4874 
4875  if(mikhailSetup) {
4876  m23Slices = 1;
4877  gaussBkgTime = true;
4878  }
4879 }
4880 
4881 void parseArg(GooFit::App *app) {
4882  app->add_option("--luckyFrac", luckyFrac, "", true);
4883  app->add_option("--mesonRad", mesonRad, "", true);
4884  app->add_option("--normBins", normBinning, "", true);
4885  app->add_option("--blindSeed", blindSeed, "", true);
4886  app->add_option("--mdslices", mdslices, "", true);
4887  app->add_option("--offset", md0offset, "Offest in GeV", true);
4888  // Previously in MeV
4889  app->add_option("--upper_window", md0_upper_window, "", true);
4890  app->add_option("--lower_window", md0_lower_window, "", true);
4891  app->add_option("--upper_delta_window", deltam_upper_window, "", true);
4892  app->add_option("--lower_delta_window", deltam_lower_window, "", true);
4893  app->add_option("--upperTime", upperTime, "", true);
4894  app->add_option("--lowerTime", lowerTime, "", true);
4895  app->add_option("--maxSigma", maxSigma, "", true);
4896  app->add_option("--polyEff", polyEff, "", true);
4897  app->add_option("--m23Slices", m23Slices, "", true);
4898  app->add_option("--bkgRandSeed", bkgHistRandSeed, "", true);
4899 
4900  app->add_flag("--drop-rho_1450", drop_rho_1450);
4901  app->add_flag("--drop-rho_1700", drop_rho_1700);
4902  app->add_flag("--drop-f0_980", drop_f0_980);
4903  app->add_flag("--drop-f0_1370", drop_f0_1370);
4904  app->add_flag("--drop-f0_1500", drop_f0_1500);
4905  app->add_flag("--drop-f0_1710", drop_f0_1710);
4906  app->add_flag("--drop-f2_1270", drop_f2_1270);
4907  app->add_flag("--drop-f0_600", drop_f0_600);
4908 
4909  app->add_flag("--histSigma", useHistogramSigma);
4910  app->add_flag("--makePlots", makePlots);
4911  app->add_set("--mkg2Model", bkg2Model_str, {"histogram", "parameter", "sideband"}, "", true);
4912  app->add_flag("--bkg3Hist", notUseBackground3Hist);
4913  app->add_flag("--bkg4Hist", notUseBackground4Hist);
4914  app->add_option("--bkgHistBins", bkgHistBins, "", true);
4915  app->add_option("--varyParameterUp", paramUp, "", true);
4916  app->add_option("--varyParameterDn", paramDn, "", true);
4917  app->add_flag("--mikhail", mikhailSetup);
4918 }
4919 
4920 int main(int argc, char **argv) {
4921  gStyle->SetCanvasBorderMode(0);
4922  gStyle->SetCanvasColor(10);
4923  gStyle->SetFrameFillColor(10);
4924  gStyle->SetFrameBorderMode(0);
4925  gStyle->SetPadColor(0);
4926  gStyle->SetTitleColor(1);
4927  gStyle->SetStatColor(0);
4928  gStyle->SetFillColor(0);
4929  gStyle->SetFuncWidth(1);
4930  gStyle->SetLineWidth(1);
4931  gStyle->SetLineColor(1);
4932  gStyle->SetPalette(1, 0);
4933  foo = new TCanvas();
4934  foodal = new TCanvas();
4935  foodal->Size(10, 10);
4936 
4937  int retval = 0;
4938 
4939  GooFit::Application app("pipipi0 Dalitz fit example", argc, argv);
4940  app_ptr = &app;
4941  app.require_subcommand();
4942 
4943  app.add_flag("--minuit1", minuit1, "Use Minuit 1 instead of Minuit 2");
4944 
4945  std::string data;
4946  int sample = 0;
4947  int load = 1;
4948  bool plots;
4949  int genResolutions = 0;
4950  double dplotres = 0;
4951 
4952  auto toy = app.add_subcommand("toy", "Toy MC Performance evaluation");
4953  toy->add_option("-s,--sample,sample", sample, "Sample number to use", true);
4954  toy->add_option("-l,--load,load", load, "Number of times to load", true);
4955  toy->add_option("-m,--max", maxEvents, "Maximum number of events to read", true);
4956  toy->add_flag("-p,--plot", plots, "Also make plots");
4957  toy->set_callback([&]() { retval = runToyFit(sample, load, plots); });
4958 
4959  auto truth_fit = app.add_subcommand("truth", "Truth Monte Carlo fit");
4960  truth_fit->add_option("-d,--data,data", data, "Data to use")->required()->check(GooFit::ExistingFile);
4961  truth_fit->set_callback([&]() { retval = runTruthMCFit(data, false); });
4962 
4963  auto sigma_fit = app.add_subcommand("sigma", "Run sigma fit");
4964  sigma_fit->add_option("-d,--data,data", data, "Data to use")->required()->check(GooFit::ExistingFile);
4965  sigma_fit->add_option("-s,--slices,slices", m23Slices, "m23 slices")->required();
4966  sigma_fit->set_callback([&]() { retval = runSigmaFit(data.c_str()); });
4967 
4968  auto efficiency_fit = app.add_subcommand("efficiency", "Run efficiency fit");
4969  efficiency_fit->add_option("-s,--sample,sample", sample, "Sample number to use", true);
4970  efficiency_fit->set_callback([&]() { retval = runEfficiencyFit(sample); });
4971 
4972  auto canonical_fit = app.add_subcommand("canonical", "Run the canonical fit");
4973  canonical_fit->add_option("-d,--data,data", data, "Data to use")->required()->check(GooFit::ExistingFile);
4974  parseArg(canonical_fit);
4975  canonical_fit->set_callback([&]() {
4977  retval = runCanonicalFit(data, !makePlots);
4978  });
4979 
4980  auto background_dalitz_fit = app.add_subcommand("background_dalitz", "Run the background Dalitz fit");
4981  background_dalitz_fit->add_option("-s,--sample,sample", sample, "Sample number to use", true);
4982  parseArg(background_dalitz_fit);
4983  background_dalitz_fit->set_callback([&]() {
4985  retval = runBackgroundDalitzFit(sample, true);
4986  });
4987 
4988  auto background_sigma_fit = app.add_subcommand("background_sigma", "Run background sigma fit");
4989  background_sigma_fit->add_option("-s,--sample,sample", sample, "Sample number to use", true);
4990  background_sigma_fit->set_callback([&]() { retval = runBackgroundSigmaFit(sample); });
4991 
4992  auto write_background_histograms = app.add_subcommand("background_histograms", "Write background histograms");
4993  write_background_histograms->add_option("-s,--sample,sample", sample, "Sample number to use", true);
4994  write_background_histograms->set_callback([&]() { writeBackgroundHistograms(sample); });
4995 
4996  auto run_gen_mc_fit = app.add_subcommand("run_gen_mc", "Run generated Monte Carlo fit");
4997  run_gen_mc_fit->add_option("-d,--data,data", data, "Data to use")->required()->check(GooFit::ExistingFile);
4998  run_gen_mc_fit->add_option("-g,--genres,gen-resolutions", genResolutions)->required();
4999  run_gen_mc_fit->add_option("-p,--dplotres,dplotres", dplotres);
5000  run_gen_mc_fit->set_callback([&]() {
5001  if(!(DplotRes & genResolutions))
5002  dplotres = 0;
5003  retval = runGeneratedMCFit(data, genResolutions, dplotres);
5004  });
5005 
5006  auto make_time_plots = app.add_subcommand("make_time_plots", "Make time plots");
5007  make_time_plots->add_option("-d,--data,data", data, "Data to use")->required()->check(GooFit::ExistingFile);
5008  make_time_plots->set_callback([&]() { makeTimePlots(data); });
5009 
5010  GOOFIT_PARSE(app);
5011 
5012  return retval;
5013 }
IncoherentSumPdf * incsum2
bool drop_rho_1450
size_t getNumEvents() const
Definition: DataSet.h:47
fptype getBinSize() const
Get the bin size, (upper-lower) / bins.
Definition: Variable.h:133
Observable * wBkg3
__host__ void printProfileInfo(bool topLevel=true)
TH2F * underlyingBins
void normalize(TH1F *dat)
void setFixed(bool fix)
Set the fixedness of a variable.
Definition: Variable.h:205
double fptype
std::vector< ResonancePdf * > resonances
void plotLoHiSigma()
vector< GooPdf * > jsuList
void makeFullFitVariables()
double maxSigma
fptype getBinContent(size_t bin) const
Get the content of a bin.
Definition: BinnedDataSet.h:26
Observable * massd0
__host__ void setParameterConstantness(bool constant=true)
void setNumBins(size_t num)
Set the number of bins.
Definition: Variable.h:128
int runGeneratedMCFit(std::string fname, int genResolutions, double dplotres)
Variable maxDalitzY("maxDalitzY", pow(_mD0 - piPlusMass, 2))
Variable constantZero("constantZero", 0)
double lowerTime
int m23Slices
double deltam_lower_window
GooPdf * makeBkg4_sigma()
GooPdf * makeFlatBkgDalitzPdf(bool fixem=true)
bool drop_f0_600
Variable neutrlM("neutrlM", 0.1349766)
void addEvent() override
bool doToyStudy
#define GOOFIT_INFO(...)
Definition: Log.h:10
bool notUseBackground3Hist
Resolutions
Nonresonant constructor.
Definition: ResonancePdf.h:117
Variable massSum("massSum", _mD0 *_mD0+2 *piPlusMass *piPlusMass+piZeroMass *piZeroMass)
__host__ Variable * getParameterByName(std::string n)
Definition: PdfBase.cpp:95
int runSigmaFit(const char *fname)
TddpPdf * makeSignalPdf(MixingTimeResolution *resolution=0, GooPdf *eff=0)
TH1F * loM23Sigma
Variable chargeM("chargeM", 0.13957018)
void makeToyDalitzPlots(GooPdf *overallSignal, std::string plotdir="./plots_from_toy_mixfit/")
Variable minDalitzX("minDalitzX", pow(piPlusMass+piZeroMass, 2))
fptype cpuGetM23(fptype massPZ, fptype massPM)
char strbuffer[1000]
std::vector< Variable > weights
void setMaxCalls(double mxc)
void coarseBin(TH2F &dalitzHist, int grain)
__host__ std::vector< std::vector< fptype > > getCompProbsAtDataPoints()
Produce a list of probabilies at points.
Bkg2Model bkg2Model
GooPdf * makeBkg_sigma_strips(int bkgnum)
Bkg2Model
bool drop_f0_1710
void parseArg(GooFit::App *app)
const fptype piPlusMass
void makeDalitzPlots(GooPdf *overallSignal, std::string plotdir="./plots_from_mixfit/")
void loadDataFile(std::string fname, UnbinnedDataSet **setToFill=0, int effSkip=3)
void writeToFile(PdfBase *pdf, const char *fname)
void setValue(fptype val)
Set the value.
Definition: Variable.h:70
bool cpuDalitz(fptype m12, fptype m13, fptype bigM, fptype dm1, fptype dm2, fptype dm3)
bool mikhailSetup
TH1F ** sigma_pdf_hists
double getContent(TH2F *plot)
const int numSigmaBins
Gaussian constructor.
Definition: ResonancePdf.h:110
GooPdf * makeEfficiencyPdf()
double intGaus
Observable * m12
int runToyFit(int ifile, int nfile, bool noPlots=true)
GooPdf * bkg3_jsugg
Special class for observables. Used in DataSets.
Definition: Variable.h:109
GooPdf * makeBkg3DalitzPdf(bool fixem=true)
const float toySigFraction
ChisqInfo * getAdaptiveChisquare(TH2F *datPlot, TH2F *pdfPlot)
Variable constantTwo("constantTwo", 2)
fptype getValue(const Observable &var, size_t idx) const
Get the value at a specific variable and event number.
Variable fixedRhoWidth("rho_width", 0.1503, 0.01, 0.1, 0.2)
Variable minDalitzY("minDalitzY", pow(piPlusMass+piZeroMass, 2))
Variable constantBigM("constantBigM", _mD02+2 *piPlusMass *piPlusMass+piZeroMass *piZeroMass)
Variable minDalitzZ("minDalitzZ", pow(piPlusMass+piPlusMass, 2))
void setBinContent(unsigned int bin, fptype value)
Definition: BinnedDataSet.h:45
bool notUseBackground4Hist
GooPdf * makeBackground3DalitzParam()
Observable * deltam
std::string toyFileName
DecayInfo3t dtop0pp
__host__ void extractHistogram(thrust::host_vector< fptype > &host_hist)
size_t getNumBins() const
Get the number of bins.
Definition: Variable.h:130
GooPdf * makeSigmaMap()
TH1F * dataTimePlot
UnbinnedDataSet * effdata
int mdslices
TH1F * hiM23Sigma
UnbinnedDataSet * data
size_t maxEvents
TddpPdf * signalDalitz
Variable constantOne("constantOne", 1)
IncoherentSumPdf * incsum3
Observable * wSig0
BinnedDataSet * binEffData
TCanvas * foodal
double mesonRad
bool drop_f2_1270
string paramUp
Gounaris-Sakurai.
Definition: ResonancePdf.h:89
const fptype _mD0
__host__ void setDataSize(unsigned int dataSize, unsigned int evtSize=5)
bool polyEff
GooPdf * make4BinSigmaMap()
__host__ void setData(DataSet *data)
Observable * m13
GooPdf * bkg4_jsugg
int runCanonicalFit(std::string fname, bool noPlots=true)
Observable * sigma
bool drop_f0_980
void setUpperLimit(fptype val)
Set the upper limit.
Definition: Variable.h:75
const float toyBkgTimeRMS
std::vector< PdfBase * > comps
int runBackgroundDalitzFit(int bkgType, bool plots=false)
void createWeightHistogram()
GooPdf * makeBkg2DalitzPdf(bool fixem=true)
double deltam_upper_window
int runBackgroundSigmaFit(int bkgType)
void getBackgroundFile(int bkgType)
const fptype _mD02inv
void writeListOfNumbers(thrust::host_vector< fptype > &target, const char *fname)
bool minuit1
GooPdf * makeMikhailJSU_gg(bool fixem=true)
const fptype piZeroMass
bool gaussBkgTime
int blindSeed
std::string bkg2Model_str
double upperTime
const float toyBkgTimeMean
GooPdf * makeOverallSignal()
EventNumber * eventNumber
virtual __host__ std::vector< Observable > getObservables() const
Definition: PdfBase.cpp:111
__host__ std::vector< fptype > evaluateAtPoints(Observable var)
GooPdf * makeBkg3_sigma()
GooPdf * makeBkg3Eff()
int bkgHistBins
bool useHistogramSigma
void plotFit(Observable var, UnbinnedDataSet *dat, GooPdf *fit)
fptype getLowerLimit() const
Get the lower limit.
Definition: Variable.h:78
GooPdf * make1BinSigmaMap()
Observable * wBkg4
const std::vector< Observable > & getObservables() const
Definition: DataSet.cpp:49
int md0_lower_window
int normBinning
Variable fixedRhoMass("rho_mass", 0.7758, 0.01, 0.7, 0.8)
double luckyFrac
fptype getValue() const
Get the value.
Definition: Variable.h:68
Variable constantMinusOne("constantMinusOne", -1)
fptype getUpperLimit() const
Get the upper limit.
Definition: Variable.h:73
bool drop_rho_1700
int main(int argc, char **argv)
void addWeightedEvent(double weight) override
int md0_upper_window
bool drop_f0_1500
bool makePlots
void setBlind(fptype val)
Hides the number; the real value is the result minus this value. Cannot be retreived once set...
Definition: Variable.h:211
double calcGauInteg(double x1, double x2)
float md0_toy_mean
Variable maxDalitzX("maxDalitzX", pow(_mD0 - piPlusMass, 2))
double md0offset
Variable motherM("motherM", 1.86484)
Variable maxDalitzZ("maxDalitzZ", pow(_mD0 - piZeroMass, 2))
GooPdf * make_m23_transform()
void setMaxCalls(unsigned int max_calls=0)
Set the maximum number of calls. 0 for Minuit2 default.
void getToyData(float sigweight=0.9)
void set_bkg_model_from_string()
const fptype _mD02
GooPdf * makeKzeroVeto()
Observable * wBkg2
void setValueForAllEvents(const Observable &var)
Set all entries to a constant value (note: this is kind of ugly)
GooPdf * makeGaussianTimePdf(int bkg)
void readFromFile(PdfBase *pdf, const char *fname)
GooPdf * makeExpGausTimePdf(int bkg)
Observable * wBkg1
ROOT::Minuit2::FunctionMinimum fit()
This runs the fit.
#define GOOFIT_PARSE(app,...)
Definition: Application.h:11
string paramDn
size_t getNumBins() const
Relativistic Breit-Wigner.
Definition: ResonancePdf.h:68
IncoherentSumPdf * incsum5
IncoherentSumPdf * incsum1
bool saveEffPlot
void makeTimePlots(std::string fname)