Belle II Software  release-08-01-10
InvariantMassMuMuStandAlone.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 
10 #include <iostream>
11 #include <iomanip>
12 #include <filesystem>
13 #include <vector>
14 #include <tuple>
15 #include <numeric>
16 #include <fstream>
17 
18 #include <TROOT.h>
19 #include <TTree.h>
20 #include <TGraph.h>
21 #include <TRandom3.h>
22 #include <TH1D.h>
23 #include <TGraph.h>
24 #include <TLegend.h>
25 #include <TLine.h>
26 #include <TCanvas.h>
27 #include <TStyle.h>
28 #include <TPad.h>
29 #include <Math/Functor.h>
30 #include <Math/SpecFuncMathCore.h>
31 #include <Math/DistFunc.h>
32 
33 #include <Eigen/Dense>
34 
35 #include <framework/particledb/EvtGenDatabasePDG.h>
36 
37 //if compiled within BASF2
38 #ifdef _PACKAGE_
39 #include <tracking/calibration/InvariantMassMuMuStandAlone.h>
40 #include <tracking/calibration/InvariantMassMuMuIntegrator.h>
41 #include <tracking/calibration/BoostVectorStandAlone.h>
42 #include <tracking/calibration/Splitter.h>
43 #include <tracking/calibration/tools.h>
44 #include <tracking/calibration/ChebFitter.h>
45 #else
46 #include <InvariantMassMuMuStandAlone.h>
47 #include <InvariantMassMuMuIntegrator.h>
48 #include <Splitter.h>
49 #include <tools.h>
50 #include <ChebFitter.h>
51 #endif
52 
53 using Eigen::MatrixXd;
54 using Eigen::VectorXd;
55 
56 
57 namespace Belle2::InvariantMassMuMuCalib {
58 
59 
60 
61 
63  std::vector<Event> getEvents(TTree* tr, bool is4S)
64  {
65 
66  std::vector<Event> events;
67  events.reserve(tr->GetEntries());
68 
69  Event evt;
70 
71  tr->SetBranchAddress("run", &evt.run);
72  tr->SetBranchAddress("exp", &evt.exp);
73  tr->SetBranchAddress("event", &evt.evtNo);
74 
75  B2Vector3D* p0 = nullptr;
76  B2Vector3D* p1 = nullptr;
77 
78  tr->SetBranchAddress("mu0_p", &p0);
79  tr->SetBranchAddress("mu1_p", &p1);
80 
81  tr->SetBranchAddress("mu0_pid", &evt.mu0.pid);
82  tr->SetBranchAddress("mu1_pid", &evt.mu1.pid);
83 
84 
85  tr->SetBranchAddress("time", &evt.t); //time in hours
86 
87  const double mMu = EvtGenDatabasePDG::Instance()->GetParticle("mu-")->Mass(); //muon mass, i.e. around 105.66e-3 [GeV]
88 
89  for (int i = 0; i < tr->GetEntries(); ++i) {
90  tr->GetEntry(i);
91 
92  evt.mu0.p = *p0;
93  evt.mu1.p = *p1;
94 
95  evt.m = sqrt(pow(hypot(evt.mu0.p.Mag(), mMu) + hypot(evt.mu1.p.Mag(), mMu), 2) - (evt.mu0.p + evt.mu1.p).Mag2());
96 
97  evt.nBootStrap = 1;
98  evt.isSig = true;
99  evt.is4S = is4S;
100  events.push_back(evt);
101  }
102 
103  //sort by time
104  sort(events.begin(), events.end(), [](Event e1, Event e2) {return e1.t < e2.t;});
105 
106 
107  return events;
108  }
109 
110 
111 
112 
113 
114  // Numerical integration using Gauss-Konrod algorithm
115  double integrate(std::function<double(double)> f, double a, double b)
116  {
117  static const std::vector<double> nodes = {
118  -0.991455371120813,
119  -0.949107912342759,
120  -0.864864423359769,
121  -0.741531185599394,
122  -0.586087235467691,
123  -0.405845151377397,
124  -0.207784955007898,
125  0.000000000000000
126  };
127 
128  static const std::vector<double> wgts = {
129  0.022935322010529,
130  0.063092092629979,
131  0.104790010322250,
132  0.140653259715525,
133  0.169004726639267,
134  0.190350578064785,
135  0.204432940075298,
136  0.209482141084728
137  };
138 
139  if (b < a) B2FATAL("Wrongly defined integration interval");
140 
141  double m = (b + a) / 2; //middle of the interval
142  double d = (b - a) / 2; //half-width of the interval
143 
144  double sum = 0;
145  for (unsigned i = 0; i < nodes.size() - 1; ++i) {
146  double x1 = m - d * nodes[i];
147  double x2 = m + d * nodes[i];
148  sum += (f(x1) + f(x2)) * wgts[i];
149  }
150 
151  //add the central point (which is not in pair)
152  sum += f(m) * wgts.back();
153 
154 
155  // scale by size of interval, for example a=-1, b=1 -> no scaling
156  return sum * d;
157  }
158 
159 
160 
161 
162 
163 
164 
165 
166 
167 
169  double gausInt(double a, double b, double c, double d)
170  {
171  double res = sqrt(M_PI) / (2 * sqrt(c)) * (TMath::Erf((b * c - d) / sqrt(c)) - TMath::Erf((a * c - d) / sqrt(c)));
172  return res;
173  }
174 
175 
178  double convGausGaus(double sK, double s, double a, double b, double m, double x)
179  {
180  a -= m;
181  b -= m;
182  x -= m;
183 
184 
185  double c = 1. / 2 * (1. / s / s + 1. / sK / sK);
186  double d = 1. / 2 * x / sK / sK;
187 
188  double Const = 1. / (2 * M_PI) * 1. / (s * sK) * exp(-1. / 2 * x * x / (s * s + sK * sK));
189 
190  double res = Const * gausInt(a, b, c, d);
191  assert(isfinite(res));
192  return res;
193  }
194 
195 
196 
197 
201  double convExpGaus(double sK, double tau, double a, double x)
202  {
203  x -= a;
204 
205  double A = 1. / sqrt(2) * (-x / sK + sK / tau);
206  double B = -x / tau + 1. / 2 * pow(sK / tau, 2);
207  double res = 0;
208  if (B > 700 || A > 20) { // safety term to deal with 0 * inf limit
209  res = 1. / (2 * tau) * 1. / sqrt(M_PI) * exp(-A * A + B) * (1 / A - 1 / 2. / pow(A, 3) + 3. / 4 / pow(A, 5));
210  } else {
211  res = 1. / (2 * tau) * TMath::Erfc(A) * exp(B);
212  }
213  assert(isfinite(res));
214  return res;
215  }
216 
217 
218  // convolution of Gaussian with exp tails with the Gaussian smearing kernel
219  double gausExpConv(double mean, double sigma, double bMean, double bDelta, double tauL, double tauR, double sigmaK, double x)
220  {
221  double bL = bMean - bDelta;
222  double bR = bMean + bDelta;
223 
224  double xL = bL * sigma + mean;
225  double xR = bR * sigma + mean;
226 
227  double iGaus = sqrt(2 * M_PI) * sigma * convGausGaus(sigmaK, sigma, xL, xR, mean, x);
228  double iRight = exp(-1. / 2 * pow(bR, 2)) * tauR * convExpGaus(sigmaK, tauR, xR, x);
229  double iLeft = exp(-1. / 2 * pow(bL, 2)) * tauL * convExpGaus(sigmaK, tauL, -xL, -x);
230 
231  return (iGaus + iLeft + iRight);
232  }
233 
235  double gausExpConvRoot(const double* par)
236  {
237  double x = par[0]; // point where the function is evaluated
238  double mean = par[1]; // mean of Gauss in "Crystal Ball" with exp tails instead of pow
239  double sigma = par[2]; // sigma of Gauss in "Crystal Ball" with exp tails instead of pow
240  double bMean = par[3]; // mean of the transition points between Gaus and exp
241  double bDelta = par[4]; // diff/2 of the transition points between Gaus and exp
242  double tauL = par[5]; // decay par of the left exp
243  double tauR = par[6]; // decay par of the right exp
244  double sigmaK = par[7]; // sigma of the Gaussian smearing Kernel
245  double sigmaA = par[8]; // sigma of the added Gaussian
246  double fA = par[9]; // fraction of the added Gaussian
247 
248  //added Gaussian
249  double G = 1. / (sqrt(2 * M_PI) * sigmaA) * exp(-1. / 2 * pow((x - mean) / sigmaA, 2));
250  return (1 - fA) * gausExpConv(mean, sigma, bMean, bDelta, tauL, tauR, sigmaK, x) + fA * G;
251  }
252 
253 
256  double gausExpPowConvRoot(const double* par)
257  {
258  double x = par[0]; // point where the function is evaluated
259  double mean = par[1]; // mean of Gauss in "Crystal Ball" with exp tails instead of pow
260  double sigma = par[2]; // sigma of Gauss in "Crystal Ball" with exp tails instead of pow
261  double bMean = par[3]; // mean of the transition points between Gaus and exp
262  double bDelta = par[4]; // diff/2 of the transition points between Gaus and exp
263  double tauL = par[5]; // decay par of the left exp
264  double tauR = par[6]; // decay par of the right exp
265  double sigmaK = par[7]; // sigma of the Gaussian smearing Kernel
266  double sigmaA = par[8]; // sigma of the added Gaussian
267  double fA = par[9]; // fraction of the added Gaussian
268 
269 
270 
271  double eCMS = par[10]; // center of mass energy of the collisions
272  double slope = par[11]; // power-slope of gen-level spectra from ISR
273  double K = par[12]; // normalisation of the part without photon ISR
274 
275  double step = 0.10; // step size in the integration
276  int N = 800 * 1. / step;
277 
278  double sum = 0;
279 
280  //calculation of the convolution using trapezium rule
281  for (int i = 0; i < N; ++i) {
282 
283  double t = eCMS - i * step;
284 
285  double y = x - t;
286  double G = 1. / (sqrt(2 * M_PI) * sigmaA) * exp(-1. / 2 * pow((y - mean) / sigmaA, 2));
287  double Core = (1 - fA) * gausExpConv(mean, sigma, bMean, bDelta, tauL, tauR, sigmaK, y) + fA * G;
288 
289  double C = (i == 0 || i == N - 1) ? 0.5 : 1;
290 
291  double Kernel;
292  if (i == 0)
293  Kernel = K * pow(step, -slope);
294  else
295  Kernel = pow(eCMS - t, -slope);
296 
297 
298  sum += Kernel * Core * step * C;
299  }
300 
301  return sum;
302  }
303 
304 
305 
306 
307 
308 
311  void plotTest()
312  {
313  double mean = 0;
314  double bMean = 0;
315  double bDelta = 5;
316  double tauL = 30;
317  double tauR = 30;
318  double sigma = 10;
319  double sigmaK = 40;
320 
321  TGraph* gr = new TGraph;
322  for (double x = -300; x < 300; x += 1) {
323 
324  double v = gausExpConv(mean, sigma, bMean, bDelta, tauL, tauR, sigmaK, x);
325 
326  gr->SetPoint(gr->GetN(), x, v);
327  }
328 
329  gr->Draw();
330  }
331 
332 
334  double gausExp(const double* par)
335  {
336  double x = par[0]; // point where the function is evaluated
337  double mean = par[1]; // mean of Gaus
338  double sigma = par[2]; // sigma of Gaus
339  double bMean = par[3]; // mean of the transition points between Gaus and exp
340  double bDelta = par[4]; // diff/2 of the transition points between Gaus and exp
341  double tauL = par[5]; // decay par of the left exp
342  double tauR = par[6]; // decay par of the right exp
343 
344  double bL = bMean - bDelta;
345  double bR = bMean + bDelta;
346 
347 
348  double r = (x - mean) / sigma;
349  if (bL <= r && r <= bR) {
350  return exp(-1. / 2 * pow(r, 2));
351  } else if (r < bL) {
352  double bp = exp(-1. / 2 * pow(bL, 2));
353  double xb = mean + bL * sigma;
354  return exp((x - xb) / tauL) * bp;
355  } else {
356  double bp = exp(-1. / 2 * pow(bR, 2));
357  double xb = mean + bR * sigma;
358  return exp(-(x - xb) / tauR) * bp;
359  }
360  }
361 
362 
363 
364 
365 
366 
369  void plotInt()
370  {
371  double mean = 4;
372  double sigma = 30;
373  double sigmaK = 30;
374  double bMean = 0;
375  double bDelta = 2.6;
376  double tau = 60;
377 
378  double x = 10300;
379 
380  TGraph* gr = new TGraph;
381  TGraph* grE = new TGraph;
382  TGraph* grRat = new TGraph;
383  TGraph* grR = new TGraph;
384 
385  double m0 = 10500;
386  double slope = 0.95;
387  double eps = 0.01;
388 
389  double frac = 0.2;
390  double sigmaE = 30;
391 
392  for (double t = eps; t < 5000; t += 1.00) {
393  double Core = gausExpConv(mean, sigma, bMean, bDelta, tau, tau, sigmaK, x + t - m0);
394  double K = (+t) >= eps ? pow(+ t, -slope) : 0;
395 
396  double fun = Core * K;
397  gr->SetPoint(gr->GetN(), t, fun);
398 
399  double s = exp(-t / tau);
400 
401  double funE = exp(-abs(x - m0 + t) / tau) * 1 / t;
402  grE->SetPoint(grE->GetN(), t, funE);
403  if (funE > 0) {
404  grRat->SetPoint(grE->GetN(), t, fun / funE);
405  grR->SetPoint(grR->GetN(), s, fun / funE);
406  }
407  }
408 
409 
410  InvariantMassMuMuIntegrator integrator;
411 
412  TGraph* grM = new TGraph; // the PDF used in the fit
413 
414  double C = 16;
415 
416  for (double xNow = 10000; xNow <= 11000; xNow += 0.1) {
417 
418  integrator.init(mean,
419  sigma,
420  sigmaK,
421  bMean,
422  bDelta,
423  tau,
424  sigmaE,
425  frac,
426  m0,
427  eps,
428  C,
429  slope,
430  xNow);
431 
432 
433  grM->SetPoint(grM->GetN(), xNow, integrator.integralKronrod(2000));
434  }
435 
436  grM->Draw();
437 
438 
439  delete gr;
440  delete grE;
441  delete grRat;
442  delete grR;
443  delete grM;
444 
445  }
446 
447 
449  double mainFunction(double xx, Pars par)
450  {
451  InvariantMassMuMuIntegrator fun;
452 
453  fun.init(par.at("mean"), // mean
454  par.at("sigma"), // sigma
455  par.at("sigma"), // sigmaK
456  par.at("bMean"), // bMean
457  par.at("bDelta"),// bDelta
458  par.at("tau"), // tau=tauL=tauR
459  par.at("sigma"), // sigmaE
460  par.at("frac"), // frac
461  par.at("m0"), // m0
462  0.1, // eps
463  par.at("C"), // C
464  par.at("slope"), // slope
465  xx); // x
466 
467  return fun.integralKronrod(2000);
468  }
469 
470 
473  std::vector<double> readEvents(const std::vector<Event>& evts, double pidCut, double a, double b)
474  {
475 
476  std::vector<double> vMass;
477  for (const auto& ev : evts) {
478 
479  //Keep only muons
480  if (ev.mu0.pid < pidCut || ev.mu1.pid < pidCut) {
481  continue;
482  }
483 
484  double m = 1e3 * ev.m; // from GeV to MeV
485  if (a < m && m < b)
486  vMass.push_back(m);
487  }
488 
489  return vMass;
490  }
491 
492 
494  static void plotMuMuFitBase(TH1D* hData, TGraph* gr, TH1D* hPull, Pars pars, Eigen::MatrixXd mat, int time)
495  {
496  bool isBatch = gROOT->IsBatch();
497  gROOT->SetBatch(kTRUE);
498 
499  gStyle->SetOptStat(0);
500 
501  TCanvas* can = new TCanvas(Form("canMuMu_%d", time), "");
502 
503  TPad* pad1 = new TPad(Form("pad1_%d", time), "", 0, 0.3, 1, 1.0);
504  TPad* pad2 = new TPad(Form("pad2_%d", time), "", 0, 0, 1, 0.3);
505 
506  pad1->SetBottomMargin(0.05);
507  pad2->SetTopMargin(0.05);
508  pad2->SetBottomMargin(0.35);
509 
510  pad1->Draw();
511  pad2->Draw();
512 
514  // Main plot
516 
517  pad1->cd();
518 
519  hData->SetMarkerStyle(kFullCircle);
520  hData->Draw();
521  gr->SetLineColor(kRed);
522  gr->SetLineWidth(2);
523  gr->Draw("same");
524  hData->GetXaxis()->SetLabelSize(0.0001);
525  hData->GetYaxis()->SetLabelSize(0.05);
526  hData->GetYaxis()->SetTitle("Number of events");
527  hData->GetYaxis()->SetTitleSize(0.05);
528  hData->GetYaxis()->SetTitleOffset(0.9);
529 
530  double mY4S = 10579.4;
531  double y = gr->Eval(mY4S);
532  TLine* line = new TLine(mY4S, 0, mY4S, y);
533  line->SetLineColor(kGreen);
534  line->SetLineWidth(2);
535  line->Draw();
536 
537 
538 
539  TLegend* leg = new TLegend(.15, .4, .35, .87);
540  int i = 0, nPars = 0;
541  for (auto p : pars) {
542  double err = sqrt(mat(i, i));
543  if (err != 0) {
544  int nDig = log10(p.second / err) + 2;
545 
546  TString s = "%s = %." + TString(Form("%d", nDig)) + "g";
547  TString dig = "%." + TString(Form("%d", nDig)) + "g";
548  TString digE = "%.2g";
549  leg->AddEntry((TObject*)0, Form("%s = " + dig + " #pm " + digE, p.first.c_str(), p.second, err), "h");
550  ++nPars;
551  }
552  ++i;
553  }
554  leg->SetTextSize(0.05);
555  leg->SetBorderSize(0);
556  leg->SetFillStyle(0);
557  leg->Draw();
558 
559 
560  double chi2 = 0;
561  for (int j = 1; j <= hPull->GetNbinsX(); ++j)
562  chi2 += pow(hPull->GetBinContent(j), 2);
563  int ndf = hPull->GetNbinsX() - nPars - 1;
564 
565 
566  TLegend* leg2 = new TLegend(.73, .75, .93, .87);
567  leg2->AddEntry((TObject*)0, Form("chi2/ndf = %.2f", chi2 / ndf), "h");
568  leg2->AddEntry((TObject*)0, Form("p = %.2g", TMath::Prob(chi2, ndf)), "h");
569 
570  leg2->SetTextSize(0.05);
571  leg2->SetBorderSize(0);
572  leg2->SetFillStyle(0);
573  leg2->Draw();
574 
575 
576  double mFit = pars.at("m0");
577  double yF = gr->Eval(mFit);
578  TLine* lineR = new TLine(mFit, 0, mFit, yF);
579  lineR->SetLineColor(kRed);
580  lineR->SetLineWidth(2);
581  lineR->Draw();
582 
583 
584 
586  // Ratio plot
588 
589  pad2->cd();
590  hPull->SetMarkerStyle(kFullCircle);
591  hPull->Draw("p");
592 
593  hPull->GetXaxis()->SetTitle("M (#mu#mu) [MeV]");
594  hPull->GetYaxis()->SetTitle("pull");
595  hPull->GetXaxis()->SetTitleSize(0.13);
596  hPull->GetXaxis()->SetTitleOffset(1.25);
597  hPull->GetXaxis()->SetLabelSize(0.13);
598  hPull->GetXaxis()->SetLabelOffset(0.05);
599  hPull->GetXaxis()->SetTickSize(0.07);
600 
601 
602  hPull->GetYaxis()->SetTitleSize(0.13);
603  hPull->GetYaxis()->SetLabelSize(0.13);
604  hPull->GetYaxis()->SetTitleOffset(0.2);
605  hPull->GetYaxis()->CenterTitle();
606 
607 
608  hPull->GetYaxis()->SetNdivisions(404);
609 
610  hPull->GetYaxis()->SetRangeUser(-5, 5);
611 
612  TGraph* grLine = new TGraph(2);
613  grLine->SetPoint(0, hPull->GetBinLowEdge(1), 0);
614  grLine->SetPoint(1, hPull->GetBinLowEdge(hPull->GetNbinsX()) + hPull->GetBinWidth(hPull->GetNbinsX()), 0);
615  grLine->SetLineWidth(2);
616  grLine->SetLineColor(kRed);
617  grLine->Draw("same");
618 
619  std::filesystem::create_directories("plotsMuMu");
620 
621  can->SaveAs(Form("plotsMuMu/mumu_%d.pdf", time));
622 
623 
624  delete leg;
625  delete leg2;
626  delete line;
627  delete lineR;
628  delete grLine;
629 
630  delete pad1;
631  delete pad2;
632  delete can;
633 
634 
635  gROOT->SetBatch(isBatch);
636  }
637 
639  static void plotMuMuFit(const std::vector<double>& data, const Pars& pars, Eigen::MatrixXd mat, double mMin, double mMax, int time)
640  {
641  const int nBins = 100;
642 
643  // Fill the data histogram
644  TH1D::SetDefaultSumw2();
645  TH1D* hData = new TH1D("hData", "", nBins, mMin, mMax);
646  TH1D* hFit = new TH1D("hFit", "", nBins, mMin, mMax);
647  TH1D* hPull = new TH1D("hPull", "", nBins, mMin, mMax);
648  hData->SetDirectory(nullptr);
649  hFit->SetDirectory(nullptr);
650  hPull->SetDirectory(nullptr);
651 
652  // fill histogram with data
653  for (auto d : data)
654  hData->Fill(d);
655 
656 
657  // construct the fitted function
658  TGraph* gr = new TGraph();
659  const double step = (mMax - mMin) / (nBins);
660 
661  for (int i = 0; i <= 2 * nBins; ++i) {
662  double m = mMin + 0.5 * step * i;
663  double V = mainFunction(m, pars);
664  gr->SetPoint(gr->GetN(), m, V);
665  }
666 
667 
668  // Calculate integrals of the fitted function within each bin
669  for (int i = 0; i < nBins; ++i) {
670  double lV = gr->GetPointY(2 * i + 0);
671  double cV = gr->GetPointY(2 * i + 1);
672  double rV = gr->GetPointY(2 * i + 2);
673 
674  double I = step / 6 * (lV + 4 * cV + rV);
675  hFit->SetBinContent(i + 1, I);
676  }
677 
678  //Normalization factor
679  double F = hData->Integral() / hFit->Integral();
680 
681  hFit->Scale(F);
682 
683  // Normalize the curve
684  for (int i = 0; i < gr->GetN(); ++i)
685  gr->SetPointY(i, gr->GetPointY(i) * F * step);
686 
687 
688  // calculate pulls
689  for (int i = 1; i <= nBins; ++i) {
690  double pull = (hData->GetBinContent(i) - hFit->GetBinContent(i)) / sqrt(hFit->GetBinContent(i));
691  hPull->SetBinContent(i, pull);
692  }
693 
694 
695  plotMuMuFitBase(hData, gr, hPull, pars, mat, time);
696 
697  delete hData;
698  delete hFit;
699  delete hPull;
700  delete gr;
701  }
702 
703 
704 
705 
708  std::pair<Pars, MatrixXd> getInvMassPars(const std::vector<Event>& evts, Pars pars, double mMin, double mMax, int bootStrap = 0)
709  {
710  bool is4S = evts[0].is4S;
711 
712  std::vector<double> dataNow = readEvents(evts, 0.9/*PIDcut*/, mMin, mMax);
713 
714 
715  // do bootStrap
716  std::vector<double> data;
717  TRandom3* rand = nullptr;
718  if (bootStrap) rand = new TRandom3(bootStrap);
719  for (auto d : dataNow) {
720  int nP = bootStrap ? rand->Poisson(1) : 1;
721  for (int i = 0; i < nP; ++i)
722  data.push_back(d);
723  }
724 
725  if (bootStrap)
726  delete rand;
727 
728  ChebFitter fitter;
729  fitter.setDataAndFunction(mainFunction, data);
730  fitter.init(256 + 1, mMin, mMax);
731 
732 
733  Pars pars0_4S = {
734  {"C", 15 },
735  {"bDelta", 1.60307 },
736  {"bMean", 0 },
737  {"frac", 0.998051 },
738  {"m0", 10570.2 },
739  {"mean", 4.13917 },
740  {"sigma", 37.0859 },
741  {"slope", 0.876812 },
742  {"tau", 99.4225}
743  };
744 
745  Pars pars0_Off = {
746  {"C", 15 },
747  {"bDelta", 2.11 },
748  {"bMean", 0 },
749  {"frac", 0.9854 },
750  {"m0", mMax - 230 },
751  {"mean", 4.13917 },
752  {"sigma", 36.4 },
753  {"slope", 0.892 },
754  {"tau", 64.9}
755  };
756 
757  if (pars.empty()) {
758  pars = is4S ? pars0_4S : pars0_Off;
759  }
760 
761 
762 
763  Limits limits = {
764  {"mean", std::make_pair(0, 0)},
765  {"sigma", std::make_pair(10, 120)},
766  {"bMean", std::make_pair(0, 0)},
767  {"bDelta", std::make_pair(0.01, 10.)},
768  {"tau", std::make_pair(20, 250)},
769  {"frac", std::make_pair(0.00, 1.0)},
770 
771  {"m0", std::make_pair(10450, 10950)},
772  {"slope", std::make_pair(0.3, 0.999)},
773  {"C", std::make_pair(0, 0)}
774  };
775 
776 
777  return fitter.fitData(pars, limits, true/*useCheb*/);
778 
779  }
780 
781 
782 
783 
784  // Returns tuple with the invariant mass parameters (cmsEnergy in GeV)
785  std::tuple<std::vector<VectorXd>, std::vector<MatrixXd>, MatrixXd> runMuMuInvariantMassAnalysis(std::vector<Event> evts,
786  const std::vector<double>& splitPoints)
787  {
788  int n = splitPoints.size() + 1;
789 
790  std::vector<VectorXd> invMassVec(n);
791  std::vector<MatrixXd> invMassVecUnc(n);
792  MatrixXd invMassVecSpred;
793 
794  std::ofstream mumuTextOut("mumuEcalib.txt", std::ios::app);
795  static int iPrint = 0;
796  if (iPrint == 0)
797  mumuTextOut << "n id t1 t2 exp1 run1 exp2 run2 Ecms EcmsUnc state" << std::endl;
798  ++iPrint;
799 
800  for (int iDiv = 0; iDiv < n; ++iDiv) {
801 
802 
803  invMassVec[iDiv].resize(1); //1D vector for center of the 1D Gauss
804  invMassVecUnc[iDiv].resize(1, 1); //1x1 matrix covariance mat of the center
805  invMassVecSpred.resize(1, 1); //1x1 matrix for spread of the 1D Gauss
806 
807  std::vector<Event> evtsNow;
808  for (auto ev : evts) {
809  double tMin = (iDiv != 0) ? splitPoints[iDiv - 1] : -1e40;
810  double tMax = (iDiv != n - 1) ? splitPoints[iDiv] : 1e40;
811  if (tMin <= ev.t && ev.t < tMax)
812  evtsNow.push_back(ev);
813 
814  }
815 
816 
817 
818  // default fitting range for the mumu invariant mass
819  double mMin = 10.2e3, mMax = 10.8e3;
820 
821  // in case of offResonance runs adjust limits from median
822  if (!evtsNow[0].is4S) {
823  std::vector<double> dataNow;
824  for (const auto& ev : evtsNow)
825  dataNow.push_back(ev.m);
826  double mMedian = 1e3 * Belle2::BoostVectorCalib::median(dataNow.data(), dataNow.size());
827  double est = mMedian + 30;
828  mMax = est + 220;
829  mMin = est - 380;
830  }
831 
832 
833 
834 
835  // number of required successful bootstrap replicas
836  const int nRep = 25;
837 
838 
839  std::vector<double> vals, errs;
840  for (int rep = 0; rep < 200; ++rep) {
841  double errEst = 50. / sqrt(evtsNow.size());
842 
843  Pars resP, inDummy;
844  MatrixXd resM;
845 
846 
847  // fit using bootstrap replica replicas, rep=0 is no replica
848  tie(resP, resM) = getInvMassPars(evtsNow, inDummy, mMin, mMax, rep);
849 
850  int ind = distance(resP.begin(), resP.find("m0"));
851  double mass = resP.at("m0");
852  double err = sqrt(resM(ind, ind));
853 
854 
855  // if there are problems with fit, try again with diffrent bootstrap replica
856  if (!(errEst < err && err < 4 * errEst))
857  continue;
858 
859  vals.push_back(mass);
860  errs.push_back(err);
861 
862  // if now bootstrapping needed, plot & break
863  if (rep == 0) {
864  plotMuMuFit(readEvents(evtsNow, 0.9/*PIDcut*/, mMin, mMax), resP, resM, mMin, mMax, int(round(evtsNow[0].t)));
865  break;
866  }
867 
868  // try fit with different input parameters, but without bootstrapping
869  Pars resP0;
870  MatrixXd resM0;
871  tie(resP0, resM0) = getInvMassPars(evtsNow, resP, mMin, mMax, 0);
872 
873  int ind0 = distance(resP0.begin(), resP0.find("m0"));
874  double mass0 = resP0.at("m0");
875  double err0 = sqrt(resM0(ind0, ind0));
876 
877  // if successful, plot & break
878  if (errEst < err0 && err0 < 4 * errEst) {
879  vals = {mass0};
880  errs = {err0};
881  plotMuMuFit(readEvents(evtsNow, 0.9/*PIDcut*/, mMin, mMax), resP0, resM0, mMin, mMax, int(round(evtsNow[0].t)));
882  break;
883  }
884 
885 
886  // if the fit was sucesfull several times only on replicas, plot & break
887  if (vals.size() >= nRep) {
888  plotMuMuFit(readEvents(evtsNow, 0.9/*PIDcut*/, mMin, mMax), resP, resM, mMin, mMax, int(round(evtsNow[0].t)));
889  break;
890  }
891 
892  }
893 
894  if (vals.size() != 1 && vals.size() != nRep)
895  B2FATAL("Inconsistency of number of results with number of replicas");
896 
897  double meanMass = accumulate(vals.begin(), vals.end(), 0.) / vals.size();
898  double meanMassUnc = accumulate(errs.begin(), errs.end(), 0.) / errs.size();
899 
900  double sum2 = 0;
901  for (auto v : vals)
902  sum2 += pow(v - meanMass, 2);
903  double errBootStrap = vals.size() > 1 ? sqrt(sum2 / (vals.size() - 1)) : 0;
904 
905  mumuTextOut << n << " " << iDiv << " " << std::setprecision(14) << evtsNow.front().t << " " << evtsNow.back().t << " " <<
906  evtsNow.front().exp << " " << evtsNow.front().run << " " << evtsNow.back().exp << " " << evtsNow.back().run << " " << meanMass
907  <<
908  " " << meanMassUnc << " " << errBootStrap << std::endl;
909 
910  // Convert to GeV
911  invMassVec[iDiv](0) = meanMass / 1e3;
912  invMassVecUnc[iDiv](0, 0) = meanMassUnc / 1e3;
913  invMassVecSpred(0, 0) = 0;
914  }
915 
916  mumuTextOut.close();
917 
918  return std::make_tuple(invMassVec, invMassVecUnc, invMassVecSpred);
919  }
920 
921 }
#define K(x)
macro autogenerated by FFTW
static EvtGenDatabasePDG * Instance()
Instance method that loads the EvtGen table.
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
std::map< std::string, double > Pars
values of parameters in ML fit
Definition: ChebFitter.h:25
std::map< std::string, std::pair< double, double > > Limits
limits of parameters in ML fit
Definition: ChebFitter.h:28