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