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
41#include <framework/particledb/EvtGenDatabasePDG.h>
42
43
44//if compiled within BASF2
45#ifdef _PACKAGE_
46#include <reconstruction/calibration/BeamSpotBoostInvMass/InvariantMassBhadStandAlone.h>
47#include <reconstruction/calibration/BeamSpotBoostInvMass/Splitter.h>
48#include <reconstruction/calibration/BeamSpotBoostInvMass/tools.h>
49#else
50#include <InvariantMassBhadStandAlone.h>
51#include <Splitter.h>
52#include <tools.h>
53#endif
54
55using Eigen::MatrixXd;
56using Eigen::VectorXd;
57
58
59namespace 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