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