Belle II Software  release-08-01-10
InvariantMassBhadStandAlone.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 
11 #include <iostream>
12 #include <iomanip>
13 #include <vector>
14 #include <filesystem>
15 
16 #include <Eigen/Dense>
17 
18 #include <TROOT.h>
19 #include <TTree.h>
20 #include <TCanvas.h>
21 #include <TH1D.h>
22 #include <TLine.h>
23 #include <Math/Functor.h>
24 #include <Math/SpecFuncMathCore.h>
25 #include <Math/DistFunc.h>
26 
27 
28 #include <RooRealVar.h>
29 #include <RooDataSet.h>
30 #include <RooGaussian.h>
31 #include <RooPlot.h>
32 #include <RooHist.h>
33 #include <RooAddPdf.h>
34 #include <RooPolyVar.h>
35 #include <RooArgusBG.h>
36 #include <RooSimultaneous.h>
37 #include <RooCategory.h>
38 #include <RooArgSet.h>
39 #include <RooConstVar.h>
40 
41 #include <framework/particledb/EvtGenDatabasePDG.h>
42 
43 
44 //if compiled within BASF2
45 #ifdef _PACKAGE_
46 #include <tracking/calibration/InvariantMassBhadStandAlone.h>
47 #include <tracking/calibration/Splitter.h>
48 #include <tracking/calibration/tools.h>
49 #else
50 #include <InvariantMassBhadStandAlone.h>
51 #include <Splitter.h>
52 #include <tools.h>
53 #endif
54 
55 using Eigen::MatrixXd;
56 using Eigen::VectorXd;
57 
58 
59 namespace Belle2::InvariantMassBhadCalib {
60 
61 
62 
64  std::vector<Event> getEvents(TTree* tr)
65  {
66 
67  std::vector<Event> events;
68  events.reserve(tr->GetEntries());
69 
70  Event evt;
71 
72  tr->SetBranchAddress("run", &evt.run);
73  tr->SetBranchAddress("exp", &evt.exp);
74  tr->SetBranchAddress("event", &evt.evtNo);
75  tr->SetBranchAddress("time", &evt.t); //time in hours
76 
77 
78  std::vector<double>* pBcms = new std::vector<double>;
79  std::vector<double>* mB = new std::vector<double>;
80  std::vector<int>* pdg = new std::vector<int>;
81  std::vector<int>* mode = new std::vector<int>;
82  std::vector<double>* Kpid = new std::vector<double>;
83  std::vector<double>* R2 = new std::vector<double>;
84  std::vector<double>* mD = new std::vector<double>;
85  std::vector<double>* dmDstar = new std::vector<double>;
86 
87 
88 
89  tr->SetBranchAddress("pBcms", &pBcms);
90  tr->SetBranchAddress("mB", &mB);
91  tr->SetBranchAddress("pdg", &pdg);
92  tr->SetBranchAddress("mode", &mode);
93  tr->SetBranchAddress("Kpid", &Kpid);
94  tr->SetBranchAddress("R2", &R2);
95  tr->SetBranchAddress("mD", &mD);
96  tr->SetBranchAddress("dmDstar", &dmDstar);
97 
98 
99 
100 
101  for (int i = 0; i < tr->GetEntries(); ++i) {
102  tr->GetEntry(i);
103 
104  int nCand = mode->size();
105  evt.cand.resize(nCand);
106 
107  for (int j = 0; j < nCand; ++j) {
108  evt.cand[j].pBcms = pBcms->at(j);
109  evt.cand[j].mB = mB->at(j);
110  evt.cand[j].pdg = pdg->at(j);
111  evt.cand[j].mode = mode->at(j);
112  evt.cand[j].Kpid = Kpid->at(j);
113  evt.cand[j].R2 = R2->at(j);
114  evt.cand[j].mD = mD->at(j);
115  evt.cand[j].dmDstar = dmDstar->at(j);
116 
117  evt.cand[j].isSig = true;
118  }
119 
120  evt.nBootStrap = 1;
121  //evt.isSig = true;
122  events.push_back(evt);
123  }
124 
125  //sort by time
126  sort(events.begin(), events.end(), [](Event e1, Event e2) {return e1.t < e2.t;});
127 
128 
129  return events;
130  }
131 
132 
133 
134  void plotArgusFit(RooDataSet* dataE0, RooAddPdf& sumB0, RooArgusBG& argus, RooGaussian& gauss, RooRealVar& eNow,
135  TString fName = "")
136  {
137  // switch to the batch mode and store the current setup
138  bool isBatch = gROOT->IsBatch();
139  gROOT->SetBatch(kTRUE);
140 
141  using namespace RooFit ;
142 
143  // --- Plot toy data and composite PDF overlaid ---
144  RooPlot* mbcframe = eNow.frame(40) ;
145  dataE0->plotOn(mbcframe) ;
146 
147  sumB0.paramOn(mbcframe, dataE0);
148 
149 
150  sumB0.plotOn(mbcframe, Components(argus), LineStyle(kDashed)) ;
151  sumB0.plotOn(mbcframe, Components(gauss), LineStyle(kDashed), LineColor(kRed)) ;
152 
153  sumB0.plotOn(mbcframe);
154 
155  mbcframe->GetXaxis()->SetTitleSize(0.0001);
156  mbcframe->GetXaxis()->SetLabelSize(0.0001);
157  mbcframe->SetTitleSize(0.0001);
158  mbcframe->GetYaxis()->SetTitleSize(0.06);
159  mbcframe->GetYaxis()->SetTitleOffset(0.7);
160  mbcframe->GetYaxis()->SetLabelSize(0.06) ;
161 
162  mbcframe->SetTitle("");
163 
164 
165  RooHist* hpull = mbcframe->pullHist() ;
166  hpull->Print();
167  RooPlot* frame3 = eNow.frame(Title(".")) ;
168  frame3->GetYaxis()->SetTitle("Pull") ;
169  frame3->GetYaxis()->SetTitleSize(0.13) ;
170 
171  frame3->GetYaxis()->SetNdivisions(504) ;
172  frame3->GetYaxis()->SetLabelSize(0.15) ;
173  frame3->GetXaxis()->SetTitleSize(0.15) ;
174  frame3->GetXaxis()->SetLabelSize(0.15) ;
175 
176  frame3->GetYaxis()->SetTitleOffset(0.3) ;
177  frame3->addPlotable(hpull, "x0 P E1") ;
178  frame3->SetMaximum(5.);
179  frame3->SetMinimum(-5.);
180 
181 
182  TCanvas* c1 = new TCanvas(rn()) ;
183  TPad* pad1 = new TPad("pad1", "pad1", 0, 0.3, 1, 1.0);
184  pad1->Draw(); // Draw the upper pad: pad1
185  pad1->cd();
186  mbcframe->Draw() ;
187 
188 
189  TLine* ll = new TLine;
190  double mRev = 10579.4e-3 / 2;
191  ll->SetLineColor(kGreen);
192  ll->DrawLine(mRev, 0, mRev, mbcframe->GetMaximum());
193 
194 
195  c1->cd(); // Go back to the main canvas before defining pad2
196  TPad* pad2 = new TPad("pad2", "pad2", 0, 0.00, 1, 0.3);
197  pad2->SetTopMargin(0.05);
198  pad2->SetBottomMargin(0.35);
199  pad2->Draw();
200  pad2->cd();
201  frame3->Draw() ;
202  c1->Update();
203 
204  TLine* l = new TLine(5279.34e-3, 0.0, 5.37, 0.0);
205  l->SetLineColor(kBlue);
206  l->SetLineWidth(3);
207  l->Draw();
208 
209 
210  if (fName != "") c1->SaveAs(fName);
211 
212 
213  delete hpull;
214  delete ll;
215  delete l;
216  delete pad1;
217  delete pad2;
218 
219  delete c1;
220 
221 
222  gROOT->SetBatch(isBatch);
223  }
224 
225 
226 
227  std::map<TString, std::pair<double, double>> argusFit(const std::vector<Event>& evts,
228  std::vector<std::pair<double, double>> limits)
229  {
230 
231  const double cMBp = EvtGenDatabasePDG::Instance()->GetParticle("B+")->Mass();
232  const double cMB0 = EvtGenDatabasePDG::Instance()->GetParticle("B0")->Mass();
233 
234  using namespace RooFit;
235 
236  RooRealVar eNow("eNow", "E^{*}_{B} [GeV]", cMBp, 5.37);
237 
238 
239  RooDataSet* dataE0 = new RooDataSet("dataE0", "dataE0", RooArgSet(eNow));
240  RooDataSet* dataEp = new RooDataSet("dataEp", "dataEp", RooArgSet(eNow));
241 
242 
243  TH1D* hDeltaE = new TH1D(rn(), "", 30, -0.05, 0.05);
244  TH1D* hMD = new TH1D(rn(), "", 30, 1.7, 1.9);
245  TH1D* hMB = new TH1D(rn(), "", 30, -0.15, 0.15);
246 
247 
248 
249  B2ASSERT("Assert the existence of the Y4S particle data", EvtGenDatabasePDG::Instance()->GetParticle("Upsilon(4S)"));
250 
251  const double cmsE0 = EvtGenDatabasePDG::Instance()->GetParticle("Upsilon(4S)")->Mass(); //Y4S mass
252 
253  for (auto event : evts) {
254  for (auto cand : event.cand) {
255 
256  double p = cand.pBcms;
257  double mInv = cand.mB;
258 
259  //get mass of B+- or B0
260  double mB = EvtGenDatabasePDG::Instance()->GetParticle(abs(cand.pdg))->Mass();
261 
262 
263  // Filling the events
264  if (1.830 < cand.mD && cand.mD < 1.894)
265  if (abs(mInv - mB) < 0.05)
266  if (cand.R2 < 0.3)
267  if ((cand.dmDstar < -10) || (0.143 < cand.dmDstar && cand.dmDstar < 0.147)) {
268  double eBC = sqrt(p * p + pow(mB, 2)); // beam constrained energy
269  if (eBC > 5.37) continue;
270 
271  for (unsigned i = 0; i < limits.size(); ++i) {
272  if (limits[i].first <= event.t && event.t < limits[i].second) {
273  if (abs(cand.pdg) == 511) {
274  eNow.setVal(eBC);
275  dataE0->add(RooArgSet(eNow));
276  } else {
277  eNow.setVal(eBC);
278  dataEp->add(RooArgSet(eNow));
279  }
280  }
281  }
282 
283  hDeltaE->Fill(eBC - cmsE0 / 2);
284  hMD->Fill(cand.mD);
285  hMB->Fill(mInv - mB);
286  }
287  }
288  }
289 
290 
291 
292  RooCategory Bcharge("sample", "sample") ;
293  Bcharge.defineType("B0") ;
294  Bcharge.defineType("Bp") ;
295 
296  RooDataSet combData("combData", "combined data", eNow, Index(Bcharge), Import("B0", *dataE0), Import("Bp", *dataEp)) ;
297 
298 
299  // --- Build Gaussian signal PDF ---
300  RooRealVar sigmean("Mean", "B^{#pm} mass", 5.29, 5.27, 5.30) ;
301  RooRealVar sigwidth("#sigma", "B^{#pm} width", 0.00237, 0.0001, 0.030) ;
302  //sigwidth.setConstant(kTRUE);
303 
304 
305  RooPolyVar sigmeanNow("sigmeanNow", "shape parameter", sigmean, RooArgSet(RooConst(0.000), RooConst(1.)));
306  RooGaussian gauss("gauss", "gaussian PDF", eNow, sigmeanNow, sigwidth) ;
307 
308  // --- Build Argus background PDF ---
309  RooRealVar argpar("Argus_param", "argus shape parameter", -150.7, -300., +50.0) ;
310  RooRealVar endpointBp("EndPointBp", "endPoint parameter", cMBp, 5.27, 5.291) ; //B+ value
311  RooRealVar endpointB0("EndPointB0", "endPoint parameter", cMB0, 5.27, 5.291) ; //B0 value
312  endpointB0.setConstant(kTRUE);
313  endpointBp.setConstant(kTRUE);
314 
315 
316  //B0 pars
317  RooPolyVar shape2B0("EndPoint2B0", "shape parameter", endpointB0, RooArgSet(RooConst(0.), RooConst(2.)));
318  RooPolyVar eNowDifB0("eNowDifB0", "eNowDifB0", eNow, RooArgSet(shape2B0, RooConst(-1.)));
319  RooArgusBG argusB0("argusB0", "Argus PDF", eNowDifB0, endpointB0, argpar) ;
320 
321  //Bp pars
322  RooPolyVar shape2Bp("EndPoint2Bp", "shape parameter", endpointBp, RooArgSet(RooConst(0.), RooConst(2.)));
323  RooPolyVar eNowDifBp("eNowDifBp", "eNowDifBp", eNow, RooArgSet(shape2Bp, RooConst(-1.)));
324  RooArgusBG argusBp("argusBp", "Argus PDF", eNowDifBp, endpointBp, argpar) ;
325 
326 
327  // --- Construct signal+background PDF ---
328  RooRealVar nsigB0("nsigB0", "#signal events", 100, 0., 100000) ;
329  RooRealVar nbkgB0("nbkgB0", "#background events", 100, 0., 500000) ;
330 
331  RooRealVar nsigBp("nsigBp", "#signal events", 100, 0., 100000) ;
332  RooRealVar nbkgBp("nbkgBp", "#background events", 100, 0., 500000) ;
333 
334 
335 
336  RooAddPdf sumB0("sumB0", "g0+a0", RooArgList(gauss, argusB0), RooArgList(nsigB0, nbkgB0)) ;
337  RooAddPdf sumBp("sumBp", "gP+aP", RooArgList(gauss, argusBp), RooArgList(nsigBp, nbkgBp)) ;
338 
339 
340  // Construct a simultaneous pdf using category sample as index
341  RooSimultaneous simPdf("simPdf", "simultaneous pdf", Bcharge) ;
342 
343  // Associate model with the physics state and model_ctl with the control state
344  simPdf.addPdf(sumB0, "B0") ;
345  simPdf.addPdf(sumBp, "Bp") ;
346 
347  simPdf.fitTo(combData);
348 
349 
350  std::map<TString, std::pair<double, double>> resMap;
351  resMap["sigmean"] = {sigmean.getValV(), sigmean.getError()};
352  resMap["sigwidth"] = {sigwidth.getValV(), sigwidth.getError()};
353  resMap["argpar"] = {argpar.getValV(), argpar.getError()};
354 
355  std::filesystem::create_directories("plotsHadBonly");
356 
357  plotArgusFit(dataE0, sumB0, argusB0, gauss, eNow, Form("plotsHadBonly/B0Single_%d.pdf", int(round(limits[0].first))));
358  plotArgusFit(dataEp, sumBp, argusBp, gauss, eNow, Form("plotsHadBonly/BpSingle_%d.pdf", int(round(limits[0].first))));
359 
360 
361  // Delete rooFit objects
362  delete dataE0;
363  delete dataEp;
364 
365 
366  delete hDeltaE;
367  delete hMD;
368  delete hMB;
369 
370  return resMap;
371  }
372 
373 
374 
375 
376 
377 
378 // Analysis itself
379  std::map<TString, std::pair<double, double>> argusFitConstrained(const std::vector<Event>& evts,
380  std::vector<std::pair<double, double>> limits,
381  std::vector<std::pair<double, double>> mumuVals, std::vector<double> startPars)
382  {
383  // Calculate eCMS/2
384  for (auto& el : mumuVals) {
385  el.first /= 2;
386  el.second /= 2;
387  }
388 
389  //starting values are from the Bonly fit
390  const double meanInit = startPars[0] / 2;
391  const double sigmaInit = startPars[2] / 2;
392  const double argparInit = startPars[4];
393 
394  double s = 0, sw = 0;
395  for (auto p : mumuVals) {
396  s += p.first / pow(p.second, 2);
397  sw += 1. / pow(p.second, 2);
398  }
399  const double mumuMean = s / sw;
400  const double shiftInit = meanInit - mumuMean;
401 
402 
403  const double cMBp = EvtGenDatabasePDG::Instance()->GetParticle("B+")->Mass();
404  const double cMB0 = EvtGenDatabasePDG::Instance()->GetParticle("B0")->Mass();
405 
406  using namespace RooFit;
407 
408  RooRealVar eNow("eNow", "E^{*}_{B} [GeV]", cMBp, 5.37);
409 
410  std::vector<RooDataSet*> dataE0(limits.size());
411  std::vector<RooDataSet*> dataEp(limits.size());
412 
413  for (unsigned i = 0; i < limits.size(); ++i) {
414  dataE0[i] = new RooDataSet(Form("dataE0_%u", i), Form("dataE0_%u", i), RooArgSet(eNow));
415  dataEp[i] = new RooDataSet(Form("dataEp_%u", i), Form("dataEp_%u", i), RooArgSet(eNow));
416  }
417 
418 
419  TH1D* hDeltaE = new TH1D(rn(), "", 30, -0.05, 0.05);
420  TH1D* hMD = new TH1D(rn(), "", 30, 1.7, 1.9);
421  TH1D* hMB = new TH1D(rn(), "", 30, -0.15, 0.15);
422 
423 
424  B2ASSERT("Assert the existence of the Y4S particle data", EvtGenDatabasePDG::Instance()->GetParticle("Upsilon(4S)"));
425 
426  const double cmsE0 = EvtGenDatabasePDG::Instance()->GetParticle("Upsilon(4S)")->Mass(); //Y4S mass
427 
428  for (auto event : evts) {
429  for (auto cand : event.cand) {
430 
431  double p = cand.pBcms;
432  double mInv = cand.mB;
433 
434  //get mass of B+- or B0
435  double mB = EvtGenDatabasePDG::Instance()->GetParticle(abs(cand.pdg))->Mass();
436 
437 
438  // Filling the events
439  if (1.830 < cand.mD && cand.mD < 1.894)
440  if (abs(mInv - mB) < 0.05)
441  if (cand.R2 < 0.3)
442  if ((cand.dmDstar < -10) || (0.143 < cand.dmDstar && cand.dmDstar < 0.147)) {
443  double eBC = sqrt(p * p + pow(mB, 2)); // beam constrained energy
444  if (eBC > 5.37) continue;
445 
446  for (unsigned i = 0; i < limits.size(); ++i) {
447  if (limits[i].first <= event.t && event.t < limits[i].second) {
448  if (abs(cand.pdg) == 511) {
449  eNow.setVal(eBC);
450  dataE0[i]->add(RooArgSet(eNow));
451  } else {
452  eNow.setVal(eBC);
453  dataEp[i]->add(RooArgSet(eNow));
454  }
455  }
456  }
457 
458  hDeltaE->Fill(eBC - cmsE0 / 2);
459  hMD->Fill(cand.mD);
460  hMB->Fill(mInv - mB);
461  }
462  }
463  }
464 
465 
466 
467  RooCategory Bcharge("sample", "sample") ;
468  for (unsigned i = 0; i < limits.size(); ++i) {
469  Bcharge.defineType(Form("B0_%u", i));
470  Bcharge.defineType(Form("Bp_%u", i));
471  }
472 
473  RooDataSet* combData = nullptr;
474 
475  std::map<std::string, RooDataSet*> dataSetMap;
476  for (unsigned i = 0; i < limits.size(); ++i) {
477  dataSetMap[Form("B0_%u", i)] = dataE0[i];
478  dataSetMap[Form("Bp_%u", i)] = dataEp[i];
479  }
480  combData = new RooDataSet("combData", "combined data", eNow, Index(Bcharge), Import(dataSetMap));
481 
482 
483  // --- Build Gaussian signal PDF ---
484 
485  RooRealVar sigwidth("#sigma", "width of B-meson energy in CMS", sigmaInit, 0.0001, 0.030) ;
486 
487  std::vector<RooRealVar*> sigmean(limits.size());
488  std::vector<RooGaussian*> gauss(limits.size());
489 
490 
491 
492  for (unsigned i = 0; i < limits.size(); ++i) {
493  sigmean[i] = new RooRealVar(Form("Mean_%u", i), "mean B-meson energy in CMS", meanInit, 5.27, 5.30) ;
494  gauss[i] = new RooGaussian(Form("gauss_%u", i), "gaussian PDF", eNow, *sigmean[i], sigwidth) ;
495  }
496 
497  // --- Build Argus background PDF ---
498 
499 
500  RooRealVar argpar("Argus_param", "argus shape parameter", argparInit, -300., +50.0) ;
501  RooRealVar endpointBp("EndPointBp", "endPoint parameter", cMBp, 5.27, 5.291) ; //B+ value
502  RooRealVar endpointB0("EndPointB0", "endPoint parameter", cMB0, 5.27, 5.291) ; //B0 value
503  endpointB0.setConstant(kTRUE);
504  endpointBp.setConstant(kTRUE);
505 
506 
507  RooRealVar zero("zero", "", 0);
508  RooRealVar two("two", "", 2);
509  RooRealVar minus("minus", "", -1);
510 
511  //B0 pars
512  RooPolyVar shape2B0("EndPoint2B0", "shape parameter", endpointB0, RooArgSet(zero, two));
513  RooPolyVar eNowDifB0("eNowDifB0", "eNowDifB0", eNow, RooArgSet(shape2B0, minus));
514  RooArgusBG argusB0("argusB0", "Argus PDF", eNowDifB0, endpointB0, argpar) ;
515 
516  //Bp pars
517  RooPolyVar shape2Bp("EndPoint2Bp", "shape parameter", endpointBp, RooArgSet(zero, two));
518  RooPolyVar eNowDifBp("eNowDifBp", "eNowDifBp", eNow, RooArgSet(shape2Bp, minus));
519  RooArgusBG argusBp("argusBp", "Argus PDF", eNowDifBp, endpointBp, argpar) ;
520 
521 
522  std::vector<RooRealVar*> nsigB0(limits.size());
523  std::vector<RooRealVar*> nbkgB0(limits.size());
524 
525  std::vector<RooRealVar*> nsigBp(limits.size());
526  std::vector<RooRealVar*> nbkgBp(limits.size());
527 
528  std::vector<RooAddPdf*> sumB0(limits.size());
529  std::vector<RooAddPdf*> sumBp(limits.size());
530 
531 
532  // --- Construct signal+background PDF ---
533  for (unsigned i = 0; i < limits.size(); ++i) {
534  nsigB0[i] = new RooRealVar(Form("nsigB0_%u", i), "#signal events", 100, 0., 100000);
535  nbkgB0[i] = new RooRealVar(Form("nbkgB0_%u", i), "#background events", 100, 0., 500000) ;
536 
537  nsigBp[i] = new RooRealVar(Form("nsigBp_%u", i), "#signal events", 100, 0., 100000) ;
538  nbkgBp[i] = new RooRealVar(Form("nbkgBp_%u", i), "#background events", 100, 0., 500000) ;
539 
540  sumB0[i] = new RooAddPdf(Form("sumB0_%u", i), "g0+a0", RooArgList(*gauss[i], argusB0), RooArgList(*nsigB0[i], *nbkgB0[i]));
541  sumBp[i] = new RooAddPdf(Form("sumBp_%u", i), "gP+aP", RooArgList(*gauss[i], argusBp), RooArgList(*nsigBp[i], *nbkgBp[i]));
542 
543  }
544 
545 
546  // Construct a simultaneous pdf using category sample as index
547  RooSimultaneous simPdf("simPdf", "simultaneous pdf", Bcharge) ;
548 
549  // Associate model with the physics state and model_ctl with the control state
550  for (unsigned i = 0; i < limits.size(); ++i) {
551  simPdf.addPdf(*sumB0[i], Form("B0_%u", i));
552  simPdf.addPdf(*sumBp[i], Form("Bp_%u", i));
553  }
554 
555  RooRealVar shift("shift", "shift to mumu", shiftInit, -30e-3, 30e-3);
556 
557  std::vector<RooGaussian*> fconstraint(limits.size());
558  std::vector<RooPolyVar*> shiftNow(limits.size());
559 
560  for (unsigned i = 0; i < limits.size(); ++i) {
561  shiftNow[i] = new RooPolyVar(Form("shiftNow_%u", i), "shiftShift", shift, RooArgSet(RooConst(mumuVals[i].first), RooConst(1.)));
562 
563  fconstraint[i] = new RooGaussian(Form("fconstraint_%u", i), "fconstraint", *sigmean[i], *shiftNow[i],
564  RooConst(mumuVals[i].second / 1.)) ;
565  }
566 
567 
568  RooArgSet constraintSet;
569  for (auto& c : fconstraint) {
570  constraintSet.add(*c);
571  }
572  simPdf.fitTo(*combData, ExternalConstraints(constraintSet));
573 
574 
575  std::map<TString, std::pair<double, double>> resMap;
576  for (unsigned i = 0; i < limits.size(); ++i) {
577  resMap[Form("sigmean_%u", i)] = {sigmean[i]->getValV(), sigmean[i]->getError()};
578  resMap[Form("pull_%u", i)] = { (sigmean[i]->getValV() - shift.getValV() - mumuVals[i].first) / mumuVals[i].second, 0};
579  }
580 
581  resMap["sigwidth"] = {sigwidth.getValV(), sigwidth.getError()};
582  resMap["argpar"] = {argpar.getValV(), argpar.getError()};
583  resMap["shift"] = {shift.getValV(), shift.getError()};
584 
585 
586  namespace fs = std::filesystem;
587  fs::create_directories("plotsHadB");
588 
589  for (unsigned i = 0; i < limits.size(); ++i) {
590  plotArgusFit(dataE0[i], *sumB0[i], argusB0, *gauss[i], eNow, Form("plotsHadB/B0_%d_%u.pdf", int(round(limits[0].first)), i));
591  plotArgusFit(dataEp[i], *sumBp[i], argusBp, *gauss[i], eNow, Form("plotsHadB/Bp_%d_%u.pdf", int(round(limits[0].first)), i));
592  }
593 
594 
595  // Delete rooFit objects
596  for (unsigned i = 0; i < limits.size(); ++i) {
597  delete dataE0[i];
598  delete dataEp[i];
599 
600  delete sigmean[i];
601  delete gauss[i];
602 
603  delete nsigB0[i];
604  delete nbkgB0[i];
605  delete nsigBp[i];
606  delete nbkgBp[i];
607  delete sumB0[i];
608  delete sumBp[i];
609 
610  delete shiftNow[i];
611  delete fconstraint[i];
612  }
613 
614  delete combData;
615 
616  delete hDeltaE;
617  delete hMD;
618  delete hMB;
619 
620 
621  return resMap;
622  }
623 
624 
625 
626 
627 
628 
629  std::vector<std::vector<double>> doBhadFit(const std::vector<Event>& evts, std::vector<std::pair<double, double>> limits,
630  std::vector<std::pair<double, double>> mumuVals, const std::vector<double>& startPars)
631  {
632 
633  auto r = argusFitConstrained(evts, limits, mumuVals, startPars);
634  assert(limits.size() == mumuVals.size());
635 
636  std::vector<std::vector<double>> result(limits.size());
637 
638 
639  for (unsigned i = 0; i < limits.size(); ++i) {
640  std::string n = Form("sigmean_%u", i);
641  std::string np = Form("pull_%u", i);
642 
643  //convert to whole eCMS (ecms, ecmsSpread, ecmsShift)
644  result[i] = {2 * r.at(n).first, 2 * r.at(n).second, 2 * r.at("sigwidth").first, 2 * r.at("sigwidth").second, 2 * r.at("shift").first, 2 * r.at("shift").second, r.at(np).first};
645  }
646 
647  return result;
648 
649  }
650 
651 
652  std::vector<double> doBhadOnlyFit(const std::vector<Event>& evts, const std::vector<std::pair<double, double>>& limits)
653  {
654 
655  auto r = argusFit(evts, limits);
656 
657  return {2 * r.at("sigmean").first, 2 * r.at("sigmean").second, 2 * r.at("sigwidth").first, 2 * r.at("sigwidth").second, r.at("argpar").first};
658 
659  }
660 
661 
662 
663 }
Belle2::EvtGenDatabasePDG::Instance
static EvtGenDatabasePDG * Instance()
Instance method that loads the EvtGen table.
Definition: EvtGenDatabasePDG.cc:25
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Belle2::rn
TString rn()
Get random string.
Definition: tools.h:38