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 + pow(mB, 2)); // 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 s += p.first / pow(p.second, 2);
398 sw += 1. / pow(p.second, 2);
399 }
400 const double mumuMean = s / sw;
401 const double shiftInit = meanInit - mumuMean;
402
403
404 const double cMBp = EvtGenDatabasePDG::Instance()->GetParticle("B+")->Mass();
405 const double cMB0 = EvtGenDatabasePDG::Instance()->GetParticle("B0")->Mass();
406
407 using namespace RooFit;
408
409 RooRealVar eNow("eNow", "E^{*}_{B} [GeV]", cMBp, 5.37);
410
411 std::vector<RooDataSet*> dataE0(limits.size());
412 std::vector<RooDataSet*> dataEp(limits.size());
413
414 for (unsigned i = 0; i < limits.size(); ++i) {
415 dataE0[i] = new RooDataSet(Form("dataE0_%u", i), Form("dataE0_%u", i), RooArgSet(eNow));
416 dataEp[i] = new RooDataSet(Form("dataEp_%u", i), Form("dataEp_%u", i), RooArgSet(eNow));
417 }
418
419
420 TH1D* hDeltaE = new TH1D(rn(), "", 30, -0.05, 0.05);
421 TH1D* hMD = new TH1D(rn(), "", 30, 1.7, 1.9);
422 TH1D* hMB = new TH1D(rn(), "", 30, -0.15, 0.15);
423
424
425 B2ASSERT("Assert the existence of the Y4S particle data", EvtGenDatabasePDG::Instance()->GetParticle("Upsilon(4S)"));
426
427 const double cmsE0 = EvtGenDatabasePDG::Instance()->GetParticle("Upsilon(4S)")->Mass(); //Y4S mass
428
429 for (auto event : evts) {
430 for (auto cand : event.cand) {
431
432 double p = cand.pBcms;
433 double mInv = cand.mB;
434
435 //get mass of B+- or B0
436 double mB = EvtGenDatabasePDG::Instance()->GetParticle(abs(cand.pdg))->Mass();
437
438
439 // Filling the events
440 if (1.830 < cand.mD && cand.mD < 1.894)
441 if (abs(mInv - mB) < 0.05)
442 if (cand.R2 < 0.3)
443 if ((cand.dmDstar < -10) || (0.143 < cand.dmDstar && cand.dmDstar < 0.147)) {
444 double eBC = sqrt(p * p + pow(mB, 2)); // beam constrained energy
445 if (eBC > 5.37) continue;
446
447 for (unsigned i = 0; i < limits.size(); ++i) {
448 if (limits[i].first <= event.t && event.t < limits[i].second) {
449 if (abs(cand.pdg) == 511) {
450 eNow.setVal(eBC);
451 dataE0[i]->add(RooArgSet(eNow));
452 } else {
453 eNow.setVal(eBC);
454 dataEp[i]->add(RooArgSet(eNow));
455 }
456 }
457 }
458
459 hDeltaE->Fill(eBC - cmsE0 / 2);
460 hMD->Fill(cand.mD);
461 hMB->Fill(mInv - mB);
462 }
463 }
464 }
465
466
467
468 RooCategory Bcharge("sample", "sample") ;
469 for (unsigned i = 0; i < limits.size(); ++i) {
470 Bcharge.defineType(Form("B0_%u", i));
471 Bcharge.defineType(Form("Bp_%u", i));
472 }
473
474 RooDataSet* combData = nullptr;
475
476 std::map<std::string, RooDataSet*> dataSetMap;
477 for (unsigned i = 0; i < limits.size(); ++i) {
478 dataSetMap[Form("B0_%u", i)] = dataE0[i];
479 dataSetMap[Form("Bp_%u", i)] = dataEp[i];
480 }
481 combData = new RooDataSet("combData", "combined data", eNow, Index(Bcharge), Import(dataSetMap));
482
483
484 // --- Build Gaussian signal PDF ---
485
486 RooRealVar sigwidth("#sigma", "width of B-meson energy in CMS", sigmaInit, 0.0001, 0.030) ;
487
488 std::vector<RooRealVar*> sigmean(limits.size());
489 std::vector<RooGaussian*> gauss(limits.size());
490
491
492
493 for (unsigned i = 0; i < limits.size(); ++i) {
494 sigmean[i] = new RooRealVar(Form("Mean_%u", i), "mean B-meson energy in CMS", meanInit, 5.27, 5.30) ;
495 gauss[i] = new RooGaussian(Form("gauss_%u", i), "gaussian PDF", eNow, *sigmean[i], sigwidth) ;
496 }
497
498 // --- Build Argus background PDF ---
499
500
501 RooRealVar argpar("Argus_param", "argus shape parameter", argparInit, -300., +50.0) ;
502 RooRealVar endpointBp("EndPointBp", "endPoint parameter", cMBp, 5.27, 5.291) ; //B+ value
503 RooRealVar endpointB0("EndPointB0", "endPoint parameter", cMB0, 5.27, 5.291) ; //B0 value
504 endpointB0.setConstant(kTRUE);
505 endpointBp.setConstant(kTRUE);
506
507
508 RooRealVar zero("zero", "", 0);
509 RooRealVar two("two", "", 2);
510 RooRealVar minus("minus", "", -1);
511
512 //B0 pars
513 RooPolyVar shape2B0("EndPoint2B0", "shape parameter", endpointB0, RooArgSet(zero, two));
514 RooPolyVar eNowDifB0("eNowDifB0", "eNowDifB0", eNow, RooArgSet(shape2B0, minus));
515 RooArgusBG argusB0("argusB0", "Argus PDF", eNowDifB0, endpointB0, argpar) ;
516
517 //Bp pars
518 RooPolyVar shape2Bp("EndPoint2Bp", "shape parameter", endpointBp, RooArgSet(zero, two));
519 RooPolyVar eNowDifBp("eNowDifBp", "eNowDifBp", eNow, RooArgSet(shape2Bp, minus));
520 RooArgusBG argusBp("argusBp", "Argus PDF", eNowDifBp, endpointBp, argpar) ;
521
522
523 std::vector<RooRealVar*> nsigB0(limits.size());
524 std::vector<RooRealVar*> nbkgB0(limits.size());
525
526 std::vector<RooRealVar*> nsigBp(limits.size());
527 std::vector<RooRealVar*> nbkgBp(limits.size());
528
529 std::vector<RooAddPdf*> sumB0(limits.size());
530 std::vector<RooAddPdf*> sumBp(limits.size());
531
532
533 // --- Construct signal+background PDF ---
534 for (unsigned i = 0; i < limits.size(); ++i) {
535 nsigB0[i] = new RooRealVar(Form("nsigB0_%u", i), "#signal events", 100, 0., 100000);
536 nbkgB0[i] = new RooRealVar(Form("nbkgB0_%u", i), "#background events", 100, 0., 500000) ;
537
538 nsigBp[i] = new RooRealVar(Form("nsigBp_%u", i), "#signal events", 100, 0., 100000) ;
539 nbkgBp[i] = new RooRealVar(Form("nbkgBp_%u", i), "#background events", 100, 0., 500000) ;
540
541 sumB0[i] = new RooAddPdf(Form("sumB0_%u", i), "g0+a0", RooArgList(*gauss[i], argusB0), RooArgList(*nsigB0[i], *nbkgB0[i]));
542 sumBp[i] = new RooAddPdf(Form("sumBp_%u", i), "gP+aP", RooArgList(*gauss[i], argusBp), RooArgList(*nsigBp[i], *nbkgBp[i]));
543
544 }
545
546
547 // Construct a simultaneous pdf using category sample as index
548 RooSimultaneous simPdf("simPdf", "simultaneous pdf", Bcharge) ;
549
550 // Associate model with the physics state and model_ctl with the control state
551 for (unsigned i = 0; i < limits.size(); ++i) {
552 simPdf.addPdf(*sumB0[i], Form("B0_%u", i));
553 simPdf.addPdf(*sumBp[i], Form("Bp_%u", i));
554 }
555
556 RooRealVar shift("shift", "shift to mumu", shiftInit, -30e-3, 30e-3);
557
558 std::vector<RooGaussian*> fconstraint(limits.size());
559 std::vector<RooPolyVar*> shiftNow(limits.size());
560
561 for (unsigned i = 0; i < limits.size(); ++i) {
562 shiftNow[i] = new RooPolyVar(Form("shiftNow_%u", i), "shiftShift", shift, RooArgSet(RooConst(mumuVals[i].first), RooConst(1.)));
563
564 fconstraint[i] = new RooGaussian(Form("fconstraint_%u", i), "fconstraint", *sigmean[i], *shiftNow[i],
565 RooConst(mumuVals[i].second / 1.)) ;
566 }
567
568
569 RooArgSet constraintSet;
570 for (auto& c : fconstraint) {
571 constraintSet.add(*c);
572 }
573 simPdf.fitTo(*combData, ExternalConstraints(constraintSet));
574
575
576 std::map<TString, std::pair<double, double>> resMap;
577 for (unsigned i = 0; i < limits.size(); ++i) {
578 resMap[Form("sigmean_%u", i)] = {sigmean[i]->getValV(), sigmean[i]->getError()};
579 resMap[Form("pull_%u", i)] = { (sigmean[i]->getValV() - shift.getValV() - mumuVals[i].first) / mumuVals[i].second, 0};
580 }
581
582 resMap["sigwidth"] = {sigwidth.getValV(), sigwidth.getError()};
583 resMap["argpar"] = {argpar.getValV(), argpar.getError()};
584 resMap["shift"] = {shift.getValV(), shift.getError()};
585
586
587 namespace fs = std::filesystem;
588 fs::create_directories("plotsHadB");
589
590 for (unsigned i = 0; i < limits.size(); ++i) {
591 plotArgusFit(dataE0[i], *sumB0[i], argusB0, *gauss[i], eNow, Form("plotsHadB/B0_%d_%u.pdf", int(round(limits[0].first)), i));
592 plotArgusFit(dataEp[i], *sumBp[i], argusBp, *gauss[i], eNow, Form("plotsHadB/Bp_%d_%u.pdf", int(round(limits[0].first)), i));
593 }
594
595
596 // Delete rooFit objects
597 for (unsigned i = 0; i < limits.size(); ++i) {
598 delete dataE0[i];
599 delete dataEp[i];
600
601 delete sigmean[i];
602 delete gauss[i];
603
604 delete nsigB0[i];
605 delete nbkgB0[i];
606 delete nsigBp[i];
607 delete nbkgBp[i];
608 delete sumB0[i];
609 delete sumBp[i];
610
611 delete shiftNow[i];
612 delete fconstraint[i];
613 }
614
615 delete combData;
616
617 delete hDeltaE;
618 delete hMD;
619 delete hMB;
620
621
622 return resMap;
623 }
624
625
626
627
628
629
630 std::vector<std::vector<double>> doBhadFit(const std::vector<Event>& evts, std::vector<std::pair<double, double>> limits,
631 std::vector<std::pair<double, double>> mumuVals, const std::vector<double>& startPars)
632 {
633
634 auto r = argusFitConstrained(evts, limits, mumuVals, startPars);
635 assert(limits.size() == mumuVals.size());
636
637 std::vector<std::vector<double>> result(limits.size());
638
639
640 for (unsigned i = 0; i < limits.size(); ++i) {
641 std::string n = Form("sigmean_%u", i);
642 std::string np = Form("pull_%u", i);
643
644 //convert to whole eCMS (ecms, ecmsSpread, ecmsShift)
645 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};
646 }
647
648 return result;
649
650 }
651
652
653 std::vector<double> doBhadOnlyFit(const std::vector<Event>& evts, const std::vector<std::pair<double, double>>& limits)
654 {
655
656 auto r = argusFit(evts, limits);
657
658 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};
659
660 }
661
662
663
664}
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