Belle II Software development
BeamSpotStandAlone.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// Program obtain the Beam Spot properties from mumu track variables
10// inspired by https://docs.belle2.org/record/1511/files/BELLE2-NOTE-TE-2019-018.pdf
11// and https://arxiv.org/pdf/1405.6569.pdf
12
13
14
15#include <iomanip>
16#include <iostream>
17#include <tuple>
18
19#include <TH1D.h>
20#include <TH2D.h>
21#include <TF1.h>
22#include <TF2.h>
23#include <TProfile.h>
24#include <TTree.h>
25#include <TGraph.h>
26#include <TStyle.h>
27#include <TCanvas.h>
28#include <TLine.h>
29#include <TRandom.h>
30#include <TRotation.h>
31#include <TString.h>
32#include <TObjString.h>
33
34
35#include <functional>
36#include <vector>
37#include <numeric>
38
39
40#include <Eigen/Dense>
41
42
43//if compiled within BASF2
44#ifdef _PACKAGE_
45#include <reconstruction/calibration/BeamSpotBoostInvMass/BeamSpotStandAlone.h>
46#include <reconstruction/calibration/BeamSpotBoostInvMass/Splitter.h>
47#include <reconstruction/calibration/BeamSpotBoostInvMass/tools.h>
48#else
49#include <BeamSpotStandAlone.h>
50#include <Splitter.h>
51#include <tools.h>
52#endif
53
54using Eigen::VectorXd;
55using Eigen::Vector3d;
56using Eigen::MatrixXd;
57using Eigen::Matrix3d;
58
59
60namespace Belle2::BeamSpotCalib {
61
62
63 inline double sqrS(double x) {return x >= 0 ? x * x : -x * x; }
64 inline double sqrtS(double x) {return x >= 0 ? sqrt(x) : -sqrt(-x); }
65
66 MatrixXd getRotatedSizeMatrix(const std::vector<double>& xySize, double zzSize, double kX, double kY);
67
68
70 struct SpotParam {
71 Spline x;
72 Spline y;
73 Spline z;
74 Spline kX;
75 Spline kY;
76
77 SpotParam() {}
78
80 void print()
81 {
82 B2INFO("x");
83 x.print();
84 B2INFO("y");
85 y.print();
86 B2INFO("z");
87 z.print();
88 B2INFO("kX");
89 kX.print();
90 B2INFO("kY");
91 kY.print();
92 }
93
98 SpotParam(const std::vector<double>& vals, const std::vector<double>& errs, const std::vector<std::vector<double>>& spls,
99 int order = 0)
100 {
101 auto getSize = [order](const std::vector<double>& sp) {
102 if (sp.size() == 0)
103 return 1;
104 else {
105 if (order == 0) {
106 B2ASSERT("There must be least one node at this place", sp.size() >= 1);
107 return int(sp.size() + 1);
108 } else if (order == 1) {
109 B2ASSERT("Must be at least two nodes in lin. spline", sp.size() >= 2);
110 return int(sp.size());
111 } else {
112 B2FATAL("Unknown order");
113 }
114 }
115 };
116
117 int nx = getSize(spls[0]);
118 int ny = getSize(spls[1]);
119 x.nodes = spls[0];
120 y.nodes = spls[1];
121 x.vals = slice(vals, 0, nx);
122 y.vals = slice(vals, nx, ny);
123 x.errs = slice(errs, 0, nx);
124 y.errs = slice(errs, nx, ny);
125
126 if (spls.size() >= 4) {
127 int nkx = getSize(spls[2]);
128 int nky = getSize(spls[3]);
129 kX.nodes = spls[2];
130 kY.nodes = spls[3];
131 kX.vals = slice(vals, nx + ny, nkx);
132 kY.vals = slice(vals, nx + ny + nkx, nky);
133
134 kX.errs = slice(errs, nx + ny, nkx);
135 kY.errs = slice(errs, nx + ny + nkx, nky);
136
137 if (spls.size() >= 5) {
138 int nz = getSize(spls[4]);
139 z.nodes = spls[4];
140 z.vals = slice(vals, nx + ny + nkx + nky, nz);
141 z.errs = slice(errs, nx + ny + nkx + nky, nz);
142 }
143 }
144 }
145
146
147 };
148
149
151 struct UnknowSpline {
152 std::vector<Spline> spls;
153 void add(const Spline& spl) { spls.push_back(spl); }
154
156 Spline getMeanSigma()
157 {
158 Spline sAvg;
159 sAvg = spls[0];
160 int nNd = spls[0].vals.size();
161 for (int k = 0; k < nNd; ++k) {
162 double s = 0, ss = 0;
163 for (unsigned i = 0; i < spls.size(); ++i) {
164 s += spls[i].vals[k];
165 }
166 s /= spls.size();
167
168 for (unsigned i = 0; i < spls.size(); ++i) {
169 ss += (spls[i].vals[k] - s) * (spls[i].vals[k] - s);
170 }
171
172 if (spls.size() > 1)
173 ss = sqrt(ss / (spls.size() - 1));
174 else
175 ss = 0;
176
177 sAvg.vals[k] = s;
178 sAvg.errs[k] = ss;
179
180 }
181 return sAvg;
182 }
183
186 Spline getLimit(double v)
187 {
188 Spline sLim = spls[0];
189
190 double indx = (spls.size() - 1) * v;
191 int nNd = spls[0].vals.size();
192 for (int k = 0; k < nNd; ++k) {
193 std::vector<double> vals;
194 for (unsigned i = 0; i < spls.size(); ++i) {
195 vals.push_back(spls[i].vals[k]);
196 }
197 sort(vals.begin(), vals.end());
198
199 int I = indx;
200 double r = indx - I;
201 sLim.vals[k] = vals[I] * (1 - r) + vals[I + 1] * r;
202 sLim.errs[k] = 0;
203 }
204 return sLim;
205 }
206
207 };
208
209
210
211
213 struct UnknowVar {
214 std::vector<double> vars;
215 void add(double x) { vars.push_back(x); }
216
218 double getMean()
219 {
220 B2ASSERT("Must be at least one replica", vars.size() >= 1);
221 return accumulate(vars.begin(), vars.end(), 0.) / vars.size();
222 }
223
225 double getSigma()
226 {
227 B2ASSERT("Must be at least one replica", vars.size() >= 1);
228 double m = getMean();
229 double s = 0;
230 for (auto x : vars)
231 s += (x - m) * (x - m);
232 if (vars.size() > 1)
233 return sqrt(s / (vars.size() - 1));
234 else
235 return 0; //dummy unc. for single replica
236 }
237
239 double getLimit(double v)
240 {
241 B2ASSERT("Must be at least one replica", vars.size() >= 1);
242 double indx = (vars.size() - 1) * v;
243 sort(vars.begin(), vars.end());
244 int I = indx;
245 double r = indx - I;
246 return vars[I] * (1 - r) + vars[I + 1] * r;
247 }
248
250 void printStat(TString n)
251 {
252 B2ASSERT("Must be at least one replica", vars.size() >= 1);
253 B2INFO(n << " : " << getMean() << "+-" << getSigma() << " : " << getLimit(0.50) << " (" << getLimit(0.16) << " , " << getLimit(
254 1 - 0.16) << " )");
255 }
256
258 std::vector<double> getStats()
259 {
260 return {getLimit(0.50), getLimit(0.16), getLimit(1 - 0.16)};
261 }
262 };
263
264
271 double getAngle(double SizeX, double SizeY, double SizeXY)
272 {
273 double C = sqrS(SizeXY);
274 // By convention the range of angle is [-pi/2, pi/2],
275 // the reason is the same as with straight line, e.g. it can point to pi/4 but also to -3/4*pi
276 double angle = 1. / 2 * atan2(2 * C, ((SizeX * SizeX) - (SizeY * SizeY)));
277 return angle;
278 }
279
280
287 std::pair<double, double> getSizeMinMax(double SizeX, double SizeY, double SizeXY)
288 {
289 double A = (SizeX * SizeX) + (SizeY * SizeY);
290 double B = (SizeX * SizeX) * (SizeY * SizeY) - (SizeXY * SizeXY * SizeXY * SizeXY);
291 double D = A * A - 4 * B;
292
293 if (D < 0) {
294 B2FATAL("Problem with D value : " << D);
295 }
296
297 double Size2Min = 2 * B / (A + sqrt(D));
298 if (abs(Size2Min) < 1e-7) Size2Min = 0;
299 if (Size2Min < 0) {
300 B2FATAL("Negative BS eigen size : " << Size2Min);
301 }
302 double Size2Max = (A + sqrt(D)) / 2;
303 return {sqrtS(Size2Min), sqrtS(Size2Max)};
304 }
305
306
313 double getSize2MinMat(double SizeXX, double SizeYY, double SizeXY)
314 {
315 double A = SizeXX + SizeYY;
316 double B = SizeXX * SizeYY - (SizeXY * SizeXY);
317 double D = A * A - 4 * B;
318
319 if (D < 0) {
320 B2FATAL("Problem with D value : " << D);
321 }
322
323 double Size2Min = 2 * B / (A + sqrt(D));
324
325 return Size2Min;
326 }
327
328
329
330
331
333 struct UnknownPars {
334 UnknowSpline x;
335 UnknowSpline y;
336 UnknowSpline kX;
337 UnknowSpline kY;
338 UnknowSpline z;
339
340 UnknowVar sizeX;
341 UnknowVar sizeY;
342 UnknowVar sizeXY;
343 UnknowVar sizeMin;
344 UnknowVar sizeMax;
345 UnknowVar xyAngle;
346
347
348 UnknowVar sizeZ;
349 UnknowVar crossAngle;
350
351 UnknowVar matXX;
352 UnknowVar matYY;
353 UnknowVar matZZ;
354 UnknowVar matXY;
355 UnknowVar matXZ;
356 UnknowVar matYZ;
357
359 void add(const SpotParam& sPar, double SizeX, double SizeY, double SizeXY, double SizeZ)
360 {
361 x.add(sPar.x);
362 y.add(sPar.y);
363 kX.add(sPar.kX);
364 kY.add(sPar.kY);
365 z.add(sPar.z);
366
367
368 sizeX.add(SizeX);
369 sizeY.add(SizeY);
370 sizeXY.add(SizeXY);
371 sizeZ.add(SizeZ);
372
373
374 // Calculate the eigen-values
375 double SizeMin, SizeMax;
376 std::tie(SizeMin, SizeMax) = getSizeMinMax(SizeX, SizeY, SizeXY);
377
378 sizeMin.add(SizeMin);
379 sizeMax.add(SizeMax);
380
381 // and angle in mrad
382 double angle = 1e3 * getAngle(SizeX, SizeY, SizeXY);
383
384
385 xyAngle.add(angle);
386
387 //Get whole cov matrix
388 MatrixXd matSize = getRotatedSizeMatrix({sqrS(SizeX), sqrS(SizeY), sqrS(SizeXY)}, sqrS(SizeZ), sPar.kX.val(sPar.kX.center()),
389 sPar.kY.val(sPar.kY.center()));
390
391 // Store elements in [um]
392 matXX.add(sqrtS(matSize(0, 0)));
393 matYY.add(sqrtS(matSize(1, 1)));
394 matZZ.add(sqrtS(matSize(2, 2)));
395 matXY.add(sqrtS(matSize(0, 1)));
396
397 matXZ.add(sqrtS(matSize(0, 2)));
398 matYZ.add(sqrtS(matSize(1, 2)));
399
400
401 // crossing-angle in mrad
402 double crossAngleVal = 1e3 * 2 * sqrtS(matSize(0, 0)) / sqrtS(matSize(2, 2));
403 crossAngle.add(crossAngleVal);
404 }
405
407 void printStat()
408 {
409 x.getMeanSigma().print("x");
410 y.getMeanSigma().print("y");
411 kX.getMeanSigma().print("kX");
412 kY.getMeanSigma().print("kY");
413 z.getMeanSigma().print("z");
414
415 sizeX.printStat("sizeX");
416 sizeY.printStat("sizeY");
417 sizeXY.printStat("sizeXY");
418
419
420 sizeMin.printStat("sizeMin");
421 sizeMax.printStat("sizeMax");
422 xyAngle.printStat("xyAngle");
423 sizeZ.printStat("sizeZ");
424 crossAngle.printStat("crossAngle");
425
426
427 matXX.printStat("matXX");
428 matYY.printStat("matYY");
429 matZZ.printStat("matZZ");
430 matXY.printStat("matXY");
431 matXZ.printStat("matXZ");
432 matYZ.printStat("matYZ");
433 }
434
436 void getOutput(std::vector<VectorXd>& vtxPos, std::vector<MatrixXd>& vtxErr, MatrixXd& sizeMat)
437 {
438 //Store the vertex position
439 int nVals = x.spls[0].vals.size();
440
441 vtxPos.clear();
442 vtxErr.clear();
443
444 const double toCm = 1e-4;
445
446 for (int i = 0; i < nVals; ++i) {
447 //vertex position
448 Vector3d vtx(x.spls[0].vals[i]*toCm, y.spls[0].vals[i]*toCm, z.spls[0].vals[i]*toCm);
449
450 //vertex error matrix (symmetric)
451 Matrix3d mS = Matrix3d::Zero();
452 mS(0, 0) = sqrS(x.spls[0].errs[i] * toCm);
453 mS(1, 1) = sqrS(y.spls[0].errs[i] * toCm);
454 mS(2, 2) = sqrS(z.spls[0].errs[i] * toCm);
455
456 vtxPos.push_back(vtx);
457 vtxErr.push_back(mS);
458 }
459
460 //BeamSpot size matrix (from boot-strap iteration 0)
461
462 sizeMat.resize(3, 3);
463 sizeMat(0, 0) = sqrS(matXX.vars[0] * toCm);
464 sizeMat(1, 1) = sqrS(matYY.vars[0] * toCm);
465 sizeMat(2, 2) = sqrS(matZZ.vars[0] * toCm);
466
467 sizeMat(0, 1) = sqrS(matXY.vars[0] * toCm);
468 sizeMat(0, 2) = sqrS(matXZ.vars[0] * toCm);
469 sizeMat(1, 2) = sqrS(matYZ.vars[0] * toCm);
470
471 sizeMat(1, 0) = sizeMat(0, 1);
472 sizeMat(2, 0) = sizeMat(0, 2);
473 sizeMat(2, 1) = sizeMat(1, 2);
474 }
475
476
478 static void setBranchVal(TTree* T, std::vector<double>* vec, TString n)
479 {
480 T->Branch(n, &vec->at(0), n + "/D");
481 T->Branch(n + "_low", &vec->at(1), n + "_low/D");
482 T->Branch(n + "_high", &vec->at(2), n + "_high/D");
483 }
484
486 static void setBranchSpline(TTree* T, Spline* spl, TString n)
487 {
488 T->Branch(n + "_nodes", &spl->nodes);
489 T->Branch(n + "_vals", &spl->vals);
490 T->Branch(n + "_errs", &spl->errs);
491 }
492
493
494
496 void save2tree(TString fName)
497 {
498
499 TTree* T = new TTree("runs", "beam conditions of runs");
500
501 int run = -99; //currently dummy
502 T->Branch("run", &run, "run/I");
503
504 Spline xAvg = x.getMeanSigma();
505 setBranchSpline(T, &xAvg, "x");
506 Spline yAvg = y.getMeanSigma();
507 setBranchSpline(T, &yAvg, "y");
508 Spline zAvg = z.getMeanSigma();
509 setBranchSpline(T, &zAvg, "z");
510
511 Spline kxAvg = kX.getMeanSigma();
512 setBranchSpline(T, &kxAvg, "kX");
513 Spline kyAvg = kY.getMeanSigma();
514 setBranchSpline(T, &kyAvg, "kY");
515
516 std::vector<double> sizeXVar = sizeX.getStats();
517 setBranchVal(T, &sizeXVar, "sizeX");
518 std::vector<double> sizeYVar = sizeY.getStats();
519 setBranchVal(T, &sizeYVar, "sizeY");
520 std::vector<double> sizeXYVar = sizeXY.getStats();
521 setBranchVal(T, &sizeXYVar, "sizeXY");
522 std::vector<double> sizeZVar = sizeZ.getStats();
523 setBranchVal(T, &sizeZVar, "sizeZ");
524
525 std::vector<double> sizeMinVar = sizeMin.getStats();
526 setBranchVal(T, &sizeMinVar, "sizeMin");
527 std::vector<double> sizeMaxVar = sizeMax.getStats();
528 setBranchVal(T, &sizeMaxVar, "sizeMax");
529 std::vector<double> xyAngleVar = xyAngle.getStats();
530 setBranchVal(T, &xyAngleVar, "xyAngle");
531
532 std::vector<double> crossAngleVar = crossAngle.getStats();
533 setBranchVal(T, &crossAngleVar, "crossAngle");
534
535
536 std::vector<double> matXXVar = matXX.getStats();
537 std::vector<double> matYYVar = matYY.getStats();
538 std::vector<double> matZZVar = matZZ.getStats();
539 std::vector<double> matXYVar = matXY.getStats();
540 std::vector<double> matXZVar = matXZ.getStats();
541 std::vector<double> matYZVar = matYZ.getStats();
542
543 setBranchVal(T, &matXXVar, "matXX");
544 setBranchVal(T, &matYYVar, "matYY");
545 setBranchVal(T, &matZZVar, "matZZ");
546
547 setBranchVal(T, &matXYVar, "matXY");
548 setBranchVal(T, &matXZVar, "matXZ");
549 setBranchVal(T, &matYZVar, "matYZ");
550
551
552 T->Fill();
553 T->SaveAs(fName);
554 }
555 };
556
557
558
567 double getZIPest(const Track& tr, double t, const SpotParam& spotPar, int nestMax = 5, int nest = 0)
568 {
569 double x0, y0;
570 if (nest < nestMax) {
571 double zIP = getZIPest(tr, t, spotPar, nestMax, nest + 1);
572
573 x0 = spotPar.x.val(t) + spotPar.kX.val(t) * (zIP - spotPar.z.val(t));
574 y0 = spotPar.y.val(t) + spotPar.kY.val(t) * (zIP - spotPar.z.val(t));
575 } else {
576 x0 = spotPar.x.val(t);
577 y0 = spotPar.y.val(t);
578 }
579
580 double f0 = tr.tanlambda * (x0 * cos(tr.phi0) + y0 * sin(tr.phi0));
581
582 return (tr.z0 + f0);
583 }
584
585
592 double getCorrD(const Track& tr, double t, const SpotParam& spotPar)
593 {
594 double zIP = getZIPest(tr, t, spotPar);
595
596 double x0 = spotPar.x.val(t) + spotPar.kX.val(t) * (zIP - spotPar.z.val(t));
597 double y0 = spotPar.y.val(t) + spotPar.kY.val(t) * (zIP - spotPar.z.val(t));
598
599 double f0 = x0 * sin(tr.phi0) - y0 * cos(tr.phi0);
600
601 return (tr.d0 - f0);
602 }
603
605 double getDtimeConst(const Track& tr, double t, const SpotParam& spotPar)
606 {
607 double zIP = getZIPest(tr, t, spotPar);
608 double zIPM = getZIPest(tr, spotPar.z.center(), spotPar);
609
610 double x0 = spotPar.x.val(t) + spotPar.kX.val(t) * (zIP - spotPar.z.val(t));
611 double y0 = spotPar.y.val(t) + spotPar.kY.val(t) * (zIP - spotPar.z.val(t));
612
613 double xM = spotPar.x.val(spotPar.x.center()) + spotPar.kX.val(spotPar.kX.center()) * (zIPM - spotPar.z.val(spotPar.z.center()));
614 double yM = spotPar.y.val(spotPar.y.center()) + spotPar.kY.val(spotPar.kY.center()) * (zIPM - spotPar.z.val(spotPar.z.center()));
615
616
617 double f0 = (x0 - xM) * sin(tr.phi0) - (y0 - yM) * cos(tr.phi0);
618
619 return (tr.d0 - f0);
620 }
621
622
629 double getCorrZ(const Track& tr, double t, const SpotParam& spotPar)
630 {
631 double zIP = getZIPest(tr, t, spotPar);
632
633 double x0 = spotPar.x.val(t) + spotPar.kX.val(t) * (zIP - spotPar.z.val(t));
634 double y0 = spotPar.y.val(t) + spotPar.kY.val(t) * (zIP - spotPar.z.val(t));
635 double z0 = spotPar.z.val(t);
636
637 double f0 = z0 - tr.tanlambda * (x0 * cos(tr.phi0) + y0 * sin(tr.phi0));
638
639 return (tr.z0 - f0);
640 }
641
642
643
644
645
650 std::pair<double, double> getStartStop(const std::vector<Event>& evts)
651 {
652 double minT = 1e20, maxT = -1e20;
653 for (const auto& ev : evts) {
654 minT = std::min(minT, ev.t);
655 maxT = std::max(maxT, ev.t);
656 }
657 return {minT, maxT};
658 }
659
661 std::vector<TString> extractFileNames(TString str)
662 {
663 std::vector<TString> files;
664 auto tempVec = str.Tokenize(",");
665 for (int i = 0; i < tempVec->GetEntries(); ++i) {
666 TString s = (static_cast<TObjString*>(tempVec->At(i)))->GetString();
667 files.push_back(s.Strip());
668 }
669 return files;
670 }
671
673 std::vector<Event> getEvents(TTree* tr)
674 {
675
676 std::vector<Event> events;
677 events.reserve(tr->GetEntries());
678
679 Event evt;
680
681 tr->SetBranchAddress("run", &evt.run);
682 tr->SetBranchAddress("exp", &evt.exp);
683 tr->SetBranchAddress("event", &evt.evtNo);
684 tr->SetBranchAddress("mu0_d0", &evt.mu0.d0);
685 tr->SetBranchAddress("mu1_d0", &evt.mu1.d0);
686 tr->SetBranchAddress("mu0_z0", &evt.mu0.z0);
687 tr->SetBranchAddress("mu1_z0", &evt.mu1.z0);
688
689 tr->SetBranchAddress("mu0_tanlambda", &evt.mu0.tanlambda);
690 tr->SetBranchAddress("mu1_tanlambda", &evt.mu1.tanlambda);
691
692
693 tr->SetBranchAddress("mu0_phi0", &evt.mu0.phi0);
694 tr->SetBranchAddress("mu1_phi0", &evt.mu1.phi0);
695
696 tr->SetBranchAddress("time", &evt.t); //time in hours
697
698
699 for (int i = 0; i < tr->GetEntries(); ++i) {
700 tr->GetEntry(i);
701 evt.toMicroM();
702
703 evt.nBootStrap = 1;
704 evt.isSig = true;
705 events.push_back(evt);
706 }
707
708 //sort by time
709 sort(events.begin(), events.end(), [](const Event & e1, const Event & e2) {return e1.t < e2.t;});
710
711
712 return events;
713 }
714
716 void bootStrap(std::vector<Event>& evts)
717 {
718 for (auto& e : evts)
719 e.nBootStrap = gRandom->Poisson(1);
720 }
721
722
723
724
730 VectorXd linearFit(MatrixXd mat, VectorXd r)
731 {
732 MatrixXd matT = mat.transpose();
733 MatrixXd A = matT * mat;
734
735 VectorXd v = matT * r;
736 MatrixXd Ainv = A.inverse();
737
738 return (Ainv * v);
739 }
740
749 std::pair<std::vector<double>, std::vector<double>> linearFitErr(MatrixXd m, VectorXd r, double& err2Mean, double& err2press,
750 double& err2pressErr)
751 {
752 MatrixXd mT = m.transpose();
753 MatrixXd mat = mT * m;
754
755 mat = mat.inverse();
756 MatrixXd A = mat * mT;
757 VectorXd res = A * r;
758 VectorXd dataH = m * res;
759
760
761 //errs
762 double err2 = (dataH - r).squaredNorm() / (r.rows() - res.rows());
763
764
765 // Get PRESS statistics of the linear fit
766 {
767 //Ahat = m*A;
768 double press = 0;
769 double press2 = 0;
770 for (int i = 0; i < r.rows(); ++i) {
771 double Ahat = 0;
772 for (int k = 0; k < m.cols(); ++k)
773 Ahat += m(i, k) * A(k, i);
774
775 double v = ((r(i) - dataH(i)) / (1 - Ahat)) * ((r(i) - dataH(i)) / (1 - Ahat));
776 press += v;
777 press2 += v * v;
778 }
779 press /= r.rows();
780 press2 /= (r.rows() - 1);
781
782 err2press = press;
783 err2pressErr = sqrt((press2 - press * press) / r.rows()) / sqrt(r.rows());
784 }
785
786
787
788 MatrixXd AT = A.transpose(); // A.T();
789 MatrixXd errMat = err2 * A * AT;
790 VectorXd errs2(errMat.rows());
791 for (int i = 0; i < errs2.rows(); ++i)
792 errs2(i) = errMat(i, i);
793
794
795 err2Mean = err2;
796
797 return {vec2vec(res), vec2vec(errs2)};
798 }
799
800
801
802
808 VectorXd linearFitPos(MatrixXd mat, VectorXd r)
809 {
810 const double s2MinLimit = 0.0025; // = pow(0.05, 2), Minimal value of the BeamSpot eigenSize
811 B2ASSERT("Matrix size for size fit should be 3", mat.cols() == 3);
812 MatrixXd matT = mat.transpose();
813 MatrixXd A = matT * mat;
814 VectorXd v = matT * r;
815 MatrixXd Ainv = A.inverse();
816
817 VectorXd pars = Ainv * v;
818
819 //If everything is OK
820 double s2Min = getSize2MinMat(pars[0], pars[1], pars[2]);
821 if (pars[0] >= 0 && pars[1] >= 0 && s2Min >= s2MinLimit)
822 return pars;
823
825 //Get the error matrix
827
828 int nDf = r.rows() - 3;
829 //Calculate average error
830 double err2 = (mat * pars - r).squaredNorm() / nDf;
831
832 MatrixXd wMat = Ainv * matT;
833 MatrixXd wMatT = wMat.transpose();
834
835 MatrixXd covMat = err2 * wMat * wMatT;
836 MatrixXd covMatI = covMat.inverse();
837
839 //Get maximum likelihood
841
842 //Maximum likelihood over eigenvector and angle
843 TF2 fEig(rn(), [covMatI, pars, s2MinLimit](const double * x, const double*) {
844 double eig1 = x[0];
845 double eig2 = s2MinLimit;
846 double phi = x[1];
847 double c = cos(phi);
848 double s = sin(phi);
849
850 VectorXd xVec(3);
851 xVec[0] = eig1 * c * c + eig2 * s * s;
852 xVec[1] = eig1 * s * s + eig2 * c * c;
853 xVec[2] = s * c * (eig1 - eig2);
854
855 //double res = covMatI.Similarity(xVec - pars);
856 double res = (xVec - pars).transpose() * covMatI * (xVec - pars);
857 return res;
858 }, s2MinLimit, 400, 0, 2 * M_PI, 0);
859
860 double eigHigh, phi;
861 fEig.GetMinimumXY(eigHigh, phi);
862
863 const double cPhi = cos(phi);
864 const double sPhi = sin(phi);
865
866 pars[0] = eigHigh * (cPhi * cPhi) + s2MinLimit * (sPhi * sPhi);
867 pars[1] = eigHigh * (sPhi * sPhi) + s2MinLimit * (cPhi * cPhi);
868 pars[2] = sPhi * cPhi * (eigHigh - s2MinLimit);
869
870
871 return pars;
872
873 }
874
875
876
877
879 TH1D* getResolution(TH2D* hRes)
880 {
881 TH1D* hSigmaAll = new TH1D(rn(), "", 50, -M_PI, M_PI);
882 for (int i = 1; i <= hRes->GetNbinsX(); ++i) {
883 TH1D* hProj = hRes->ProjectionY(rn(), i, i);
884 double rms = hProj->GetRMS();
885 double rmsErr = hProj->GetRMSError();
886 hSigmaAll->SetBinContent(i, rms * rms); //from cm2 to um2
887 hSigmaAll->SetBinError(i, 2 * abs(rms)*rmsErr);
888 }
889 return hSigmaAll;
890 }
891
893 TH1D* getMean(const TH2D* hRes)
894 {
895 TH1D* hMean = new TH1D(rn(), "", 50, -M_PI, M_PI);
896 for (int i = 1; i <= hRes->GetNbinsX(); ++i) {
897 TH1D* hProj = hRes->ProjectionY(rn(), i, i);
898 double mean = hProj->GetMean();
899 double meanErr = hProj->GetMeanError();
900 hMean->SetBinContent(i, mean);
901 hMean->SetBinError(i, meanErr);
902 }
903 return hMean;
904 }
905
907 double getD12th(const Event& e, const std::vector<double>& sizesXY)
908 {
909 double sxx = sizesXY[0];
910 double syy = sizesXY[1];
911 double sxy = sizesXY[2];
912
913 double cc = cos(e.mu0.phi0) * cos(e.mu1.phi0);
914 double ss = sin(e.mu0.phi0) * sin(e.mu1.phi0);
915 double sc = -(sin(e.mu0.phi0) * cos(e.mu1.phi0) + sin(e.mu1.phi0) * cos(e.mu0.phi0));
916
917 return ss * sxx + cc * syy + sc * sxy;
918 }
919
921 double getZ12th(const Event& e, const std::vector<double>& sizesXY)
922 {
923 double sxx = sizesXY[0];
924 double syy = sizesXY[1];
925
926 double corr = e.mu0.tanlambda * e.mu1.tanlambda * (sxx * cos(e.mu0.phi0) * cos(e.mu1.phi0) + syy * sin(e.mu0.phi0) * sin(
927 e.mu1.phi0) +
928 + (sin(e.mu0.phi0) * cos(e.mu1.phi0) + cos(e.mu0.phi0) * sin(e.mu1.phi0)));
929 return corr;
930 }
931
932
933
934
935
936
938 void plotSpotPositionFit(const std::vector<Event>& evts, const SpotParam& par, TString fName)
939 {
940 TGraph* gr = new TGraph();
941 TProfile* dProf = new TProfile(rn(), "dProf", 100, -M_PI, M_PI, "S");
942 TProfile* dProfRes = new TProfile(rn(), "dProfRes", 100, -M_PI, M_PI, "S");
943
944 for (auto e : evts) {
945 if (!e.isSig) continue;
946
947 double d1 = getDtimeConst(e.mu0, e.t, par);
948 double d2 = getDtimeConst(e.mu1, e.t, par);
949
950 gr->SetPoint(gr->GetN(), e.mu0.phi0, d1);
951 gr->SetPoint(gr->GetN(), e.mu1.phi0, d2);
952
953 dProf->Fill(e.mu0.phi0, d1);
954 dProf->Fill(e.mu1.phi0, d2);
955
956
957 double d1c = getCorrD(e.mu0, e.t, par);
958 double d2c = getCorrD(e.mu1, e.t, par);
959
960 dProfRes->Fill(e.mu0.phi0, d1c);
961 dProfRes->Fill(e.mu1.phi0, d2c);
962 }
963 TF1* f = new TF1(rn(), "[0]*sin(x) - [1]*cos(x)", -M_PI, M_PI);
964 f->SetParameters(par.x.val(par.x.center()), par.y.val(par.y.center()));
965
966 TCanvas* c = new TCanvas(rn(), "");
967 gr->Draw("a p");
968 gr->GetXaxis()->SetRangeUser(-M_PI, M_PI);
969 gr->SetMaximum(+1.3 * f->GetMaximum());
970 gr->SetMinimum(-1.3 * f->GetMaximum());
971
972 gr->GetXaxis()->SetTitle("#phi_{0} [rad]");
973 gr->GetYaxis()->SetTitle("d_{0} [#mum]");
974
975 f->Draw("same");
976 c->SaveAs(fName + "_dots.pdf");
977
978
979 TCanvas* d = new TCanvas(rn(), "");
980 gStyle->SetOptStat(0);
981 dProf->Draw();
982 dProf->GetXaxis()->SetRangeUser(-M_PI, M_PI);
983 dProf->SetMaximum(+1.3 * f->GetMaximum());
984 dProf->SetMinimum(-1.3 * f->GetMaximum());
985
986 dProf->GetXaxis()->SetTitle("#phi_{0} [rad]");
987 dProf->GetYaxis()->SetTitle("d_{0} [#mum]");
988
989 f->Draw("same");
990
991 B2INFO("Saving " << fName << " prof ");
992 d->SaveAs(fName + "_prof.pdf");
993
994
995 TCanvas* e = new TCanvas(rn(), "");
996 gStyle->SetOptStat(0);
997 dProfRes->Draw();
998 dProfRes->GetXaxis()->SetRangeUser(-M_PI, M_PI);
999
1000 dProfRes->GetXaxis()->SetTitle("#phi_{0} [rad]");
1001 dProfRes->GetYaxis()->SetTitle("d_{0} res [#mum]");
1002
1003 TH1D* errP = new TH1D(rn(), "dErrP", 100, -M_PI, M_PI);
1004 TH1D* errM = new TH1D(rn(), "dErrM", 100, -M_PI, M_PI);
1005 for (int i = 1; i <= errP->GetNbinsX(); ++i) {
1006 errP->SetBinContent(i, dProfRes->GetBinError(i));
1007 errM->SetBinContent(i, -dProfRes->GetBinError(i));
1008 }
1009
1010 errP->Draw("hist same");
1011 errM->Draw("hist same");
1012
1013 f->SetParameters(0, 0);
1014 f->Draw("same");
1015
1016 B2INFO("Saving " << fName << " profRes ");
1017 e->SaveAs(fName + "_profRes.pdf");
1018
1019 }
1020
1021
1023 void plotSpotZPositionFit(const std::vector<Event>& evts, const SpotParam& par, TString fName)
1024 {
1025 TProfile* zProf = new TProfile(rn(), "dProf", 100, -M_PI, M_PI, "S");
1026 TGraph* gr = new TGraph();
1027 for (auto e : evts) {
1028 if (!e.isSig) continue;
1029
1030
1031 double z1ip = getZIPest(e.mu0, e.t, par);
1032 double z2ip = getZIPest(e.mu1, e.t, par);
1033
1034 double zipT = par.z.val(e.t);
1035 double zipM = par.z.val(par.z.center());
1036
1037 double val1 = z1ip - (zipT - zipM);
1038 double val2 = z2ip - (zipT - zipM);
1039
1040 gr->SetPoint(gr->GetN(), e.mu0.phi0, val1);
1041 gr->SetPoint(gr->GetN(), e.mu1.phi0, val2);
1042
1043 zProf->Fill(e.mu0.phi0, val1);
1044 zProf->Fill(e.mu1.phi0, val2);
1045 }
1046 TF1* f = new TF1(rn(), "[0]", -M_PI, M_PI);
1047 f->SetParameter(0, par.z.val(par.z.center()));
1048
1049 TCanvas* c = new TCanvas(rn(), "");
1050 c->SetLeftMargin(0.12);
1051 gr->Draw("a p");
1052 gr->GetXaxis()->SetRangeUser(-M_PI, M_PI);
1053 gr->SetMaximum(1000);
1054 gr->SetMinimum(-1000);
1055
1056 gr->GetXaxis()->SetTitle("#phi_{0} [rad]");
1057 gr->GetYaxis()->SetTitle("z_{IP} estimate [#mum]");
1058
1059 f->Draw("same");
1060 c->SaveAs(fName + "_dots.pdf");
1061
1062
1063 TCanvas* d = new TCanvas(rn(), "");
1064 gStyle->SetOptStat(0);
1065 d->SetLeftMargin(0.12);
1066 zProf->Draw();
1067 zProf->GetXaxis()->SetRangeUser(-M_PI, M_PI);
1068 zProf->SetMaximum(1000);
1069 zProf->SetMinimum(-1000);
1070
1071 zProf->GetXaxis()->SetTitle("#phi_{0} [rad]");
1072 zProf->GetYaxis()->SetTitle("z_{IP} estimate [#mum]");
1073
1074 f->Draw("same");
1075 d->SaveAs(fName + +"_prof.pdf");
1076
1077 }
1078
1079
1080
1081
1082
1084 void plotSpotPositionPull(const std::vector<Event>& evts, const SpotParam& par, TString fName, double cut = 70)
1085 {
1086 TH1D* hPull = new TH1D(rn(), "", 200, -200, 200);
1087
1088 for (const auto& e : evts) {
1089 if (!e.isSig) continue;
1090
1091 double d0 = getCorrD(e.mu0, e.t, par);
1092 double d1 = getCorrD(e.mu1, e.t, par);
1093
1094 hPull->Fill(d0);
1095 hPull->Fill(d1);
1096 }
1097
1098 TCanvas* c = new TCanvas(rn(), "");
1099 gStyle->SetOptStat(2210);
1100 hPull->Draw();
1101
1102 hPull->GetXaxis()->SetTitle("pull [#mum]");
1103 hPull->GetYaxis()->SetTitle("#tracks");
1104
1105 TLine* l = new TLine;
1106 l->SetLineColor(kRed);
1107 l->DrawLine(-cut, 0, -cut, 500);
1108 l->DrawLine(+cut, 0, +cut, 500);
1109
1110 c->SaveAs(fName + ".pdf");
1111 }
1112
1113
1115 void plotKxKyFit(const std::vector<Event>& evts, const SpotParam& par, TString fName)
1116 {
1117 TProfile* profRes = new TProfile(rn(), "dProf", 100, -800, 800, "S");
1118 TProfile* profResKx = new TProfile(rn(), "dProfKx", 100, -800, 800, "S");
1119 TProfile* profResKy = new TProfile(rn(), "dProfKy", 100, -800, 800, "S");
1120
1121 SpotParam parNoKx = par;
1122 SpotParam parNoKy = par;
1123 parNoKx.kX.vals = {0};
1124 parNoKy.kY.vals = {0};
1125 parNoKx.kX.nodes = {};
1126 parNoKy.kY.nodes = {};
1127
1128 for (const auto& e : evts) {
1129 if (!e.isSig) continue;
1130
1131 double zDiff1 = getZIPest(e.mu0, e.t, par) - par.z.val(e.t);
1132 double zDiff2 = getZIPest(e.mu1, e.t, par) - par.z.val(e.t);
1133
1134
1135 double d1 = getCorrD(e.mu0, e.t, par);
1136 double d2 = getCorrD(e.mu1, e.t, par);
1137
1138 double d1KX = getCorrD(e.mu0, e.t, parNoKx);
1139 double d2KX = getCorrD(e.mu1, e.t, parNoKx);
1140
1141
1142 profResKx->Fill(zDiff1 * sin(e.mu0.phi0), d1KX);
1143 profResKx->Fill(zDiff2 * sin(e.mu1.phi0), d2KX);
1144
1145 double d1KY = getCorrD(e.mu0, e.t, parNoKy);
1146 double d2KY = getCorrD(e.mu1, e.t, parNoKy);
1147 profResKy->Fill(-zDiff1 * cos(e.mu0.phi0), d1KY);
1148 profResKy->Fill(-zDiff2 * cos(e.mu1.phi0), d2KY);
1149
1150 profRes->Fill(zDiff1, d1);
1151 profRes->Fill(zDiff2, d2);
1152
1153
1154 }
1155
1156 TCanvas* cX = new TCanvas(rn(), "");
1157 gStyle->SetOptStat(0);
1158 profResKx->Draw();
1159
1160 profResKx->GetXaxis()->SetTitle("(z_{IP} - z_{IP}^{0}) sin #phi_{0} [#mum]");
1161 profResKx->GetYaxis()->SetTitle("d_{0} res [#mum]");
1162
1163 TF1* f = new TF1(rn(), "[0]*x", -800, 800);
1164 f->SetParameter(0, par.kX.val(par.kX.center()));
1165 f->Draw("same");
1166
1167 cX->SaveAs(fName + "_kX.pdf");
1168
1169 TCanvas* cY = new TCanvas(rn(), "");
1170 gStyle->SetOptStat(0);
1171 profResKy->Draw();
1172
1173 profResKy->GetXaxis()->SetTitle("-(z_{IP} - z_{IP}^{0}) cos #phi_{0} [#mum]");
1174 profResKy->GetYaxis()->SetTitle("d_{0} res [#mum]");
1175
1176 f->SetParameter(0, par.kY.val(par.kY.center()));
1177 f->Draw("same");
1178
1179 cY->SaveAs(fName + "_kY.pdf");
1180
1181
1182 }
1183
1185 void plotXYtimeDep(const std::vector<Event>& evts, const SpotParam& par, TString fName)
1186 {
1187 TProfile* profRes = new TProfile(rn(), "dProf", 50, -0.5, 0.5);
1188 TProfile* profResTx = new TProfile(rn(), "dProfTx", 50, -0.5, 0.5);
1189 TProfile* profResTy = new TProfile(rn(), "dProfTy", 50, -0.5, 0.5);
1190
1191 SpotParam parNoTx = par;
1192 SpotParam parNoTy = par;
1193 parNoTx.x.nodes = {};
1194 parNoTx.x.vals = {par.x.val(par.x.center())};
1195 parNoTy.y.nodes = {};
1196 parNoTy.y.vals = {par.y.val(par.y.center())};
1197
1198 for (const auto& e : evts) {
1199 if (!e.isSig) continue;
1200
1201 double tDiff = (e.t - par.x.val(par.x.center()));
1202 //double tDiff = e.t;
1203
1204
1205 double d1 = getCorrD(e.mu0, e.t, par);
1206 double d2 = getCorrD(e.mu1, e.t, par);
1207
1208 double d1TX = getCorrD(e.mu0, e.t, parNoTx);
1209 double d2TX = getCorrD(e.mu1, e.t, parNoTx);
1210
1211 profResTx->Fill(tDiff * sin(e.mu0.phi0), d1TX);
1212 profResTx->Fill(tDiff * sin(e.mu1.phi0), d2TX);
1213
1214 double d1TY = getCorrD(e.mu0, e.t, parNoTy);
1215 double d2TY = getCorrD(e.mu1, e.t, parNoTy);
1216
1217 profResTy->Fill(-tDiff * cos(e.mu0.phi0), d1TY);
1218 profResTy->Fill(-tDiff * cos(e.mu1.phi0), d2TY);
1219
1220
1221 profRes->Fill(tDiff * sin(e.mu0.phi0), d1);
1222 profRes->Fill(tDiff * sin(e.mu1.phi0), d2);
1223
1224
1225 }
1226
1227 TCanvas* cX = new TCanvas(rn(), "");
1228 gStyle->SetOptStat(0);
1229 profResTx->Draw();
1230
1231
1232 TF1* f = new TF1(rn(), "[0]*x", -1, 1);
1233 f->SetParameter(0, (par.x.val(1) - par.x.val(0)));
1234 f->Draw("same");
1235 B2INFO("Table value " << par.x.val(1) - par.x.val(0));
1236
1237 cX->SaveAs(fName + "_tX.pdf");
1238
1239
1240 }
1241
1242
1243
1244
1245
1247 void plotSpotZpositionPull(const std::vector<Event>& evts, const SpotParam& par, TString fName, double cut = 1000)
1248 {
1249 TH1D* hPull = new TH1D(rn(), "", 200, -2000, 2000);
1250
1251 for (const auto& e : evts) {
1252 if (!e.isSig) continue;
1253
1254 double z0 = getCorrZ(e.mu0, e.t, par);
1255 double z1 = getCorrZ(e.mu1, e.t, par);
1256
1257 hPull->Fill(z0);
1258 hPull->Fill(z1);
1259 }
1260
1261 gStyle->SetOptStat(2210);
1262 TCanvas* c = new TCanvas(rn(), "");
1263 hPull->Draw();
1264
1265 hPull->GetXaxis()->SetTitle("pull [#mum]");
1266 hPull->GetYaxis()->SetTitle("#tracks");
1267
1268 TLine* l = new TLine;
1269 l->SetLineColor(kRed);
1270 l->DrawLine(-cut, 0, -cut, 500);
1271 l->DrawLine(+cut, 0, +cut, 500);
1272
1273 c->SaveAs(fName + ".pdf");
1274 }
1275
1276
1278 void removeSpotPositionOutliers(std::vector<Event>& evts, const SpotParam& par, double cut = 70)
1279 {
1280 int nRem = 0;
1281 int nAll = 0;
1282 for (auto& e : evts) {
1283 if (!e.isSig) continue;
1284
1285 double d0 = getCorrD(e.mu0, e.t, par);
1286 double d1 = getCorrD(e.mu1, e.t, par);
1287
1288 e.isSig = abs(d0) < cut && abs(d1) < cut;
1289 nRem += !e.isSig;
1290 ++nAll;
1291 }
1292 B2INFO("Removed fraction Position " << nRem / (nAll + 0.));
1293 }
1294
1295
1297 void removeSpotZpositionOutliers(std::vector<Event>& evts, const SpotParam& par, double cut = 1000)
1298 {
1299 int nRem = 0;
1300 int nAll = 0;
1301 for (auto& e : evts) {
1302 if (!e.isSig) continue;
1303
1304 double z0 = getCorrZ(e.mu0, e.t, par);
1305 double z1 = getCorrZ(e.mu1, e.t, par);
1306
1307 e.isSig = abs(z0) < cut && abs(z1) < cut;
1308 nRem += !e.isSig;
1309 ++nAll;
1310 }
1311 B2INFO("Removed fraction Position " << nRem / (nAll + 0.));
1312 }
1313
1314
1315
1317 std::vector<std::vector<double>> fillSplineBasesLinear(const std::vector<Event>& evts, const std::vector<double>& spl,
1318 std::function<double(Track, double)> fun)
1319 {
1320 int n = spl.size(); //number of params
1321 if (n == 0 || (n == 2 && spl[0] > spl[1]))
1322 n = 1;
1323
1324 std::vector<std::vector<double>> vecs(n);
1325
1326 if (n == 1) { //no time dependence
1327 for (const auto& e : evts) {
1328 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1329 vecs[0].push_back(1 * fun(e.mu0, e.t));
1330 vecs[0].push_back(1 * fun(e.mu1, e.t));
1331 }
1332 }
1333 } else {
1334 for (int k = 0; k < n; ++k) {
1335 double xCnt = spl[k];
1336 double xLow = (k == 0) ? spl[0] : spl[k - 1];
1337 double xHigh = (k == n - 1) ? spl[n - 1] : spl[k + 1];
1338
1339 for (const auto& e : evts) {
1340 double x = e.t;
1341 double v = 0;
1342 if (xLow <= x && x < xCnt)
1343 v = (x - xLow) / (xCnt - xLow);
1344 else if (xCnt < x && x <= xHigh)
1345 v = (xHigh - x) / (xHigh - xCnt);
1346
1347
1348 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1349 vecs[k].push_back(v * fun(e.mu0, e.t));
1350 vecs[k].push_back(v * fun(e.mu1, e.t));
1351 }
1352 }
1353 }
1354 }
1355
1356 return vecs;
1357 }
1358
1359
1361 std::vector<std::vector<double>> fillSplineBasesZero(const std::vector<Event>& evts, const std::vector<double>& spl,
1362 std::function<double(Track, double)> fun)
1363 {
1364 int n = spl.size() + 1; //number of params
1365
1366 std::vector<std::vector<double>> vecs(n);
1367
1368 if (n == 1) { //no time dependence
1369 for (const auto& e : evts) {
1370 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1371 vecs[0].push_back(1 * fun(e.mu0, e.t));
1372 vecs[0].push_back(1 * fun(e.mu1, e.t));
1373 }
1374 }
1375 } else { // at least two parameters
1376 for (int k = 0; k < n; ++k) { //loop over spline parameters
1377 double xLow = -1e30;
1378 double xHigh = +1e30;
1379
1380 if (k == 0) {
1381 xHigh = spl[0];
1382 } else if (k == n - 1) {
1383 xLow = spl.back();
1384 } else {
1385 xLow = spl[k - 1];
1386 xHigh = spl[k];
1387 }
1388
1389 for (const auto& e : evts) {
1390 double x = e.t;
1391 double v = 0;
1392 if (xLow <= x && x < xHigh)
1393 v = 1;
1394
1395 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1396 vecs[k].push_back(v * fun(e.mu0, e.t));
1397 vecs[k].push_back(v * fun(e.mu1, e.t));
1398 }
1399 }
1400 }
1401 }
1402
1403 return vecs;
1404 }
1405
1406
1407
1408
1409
1410
1412 double compareSplines(const Spline& spl1, const Spline& spl2)
1413 {
1414 double sum = 0;
1415
1416 double step = 0.001;
1417 for (double x = 0; x <= 1 + step / 2; x += step) {
1418 double v1 = spl1.val(x);
1419 double e1 = spl1.err(x);
1420 double v2 = spl2.val(x);
1421 double e2 = spl2.err(x);
1422
1423 const double maxe = std::max(e1, e2);
1424 double d = (v2 - v1) * (v2 - v1) / (maxe * maxe);
1425 sum += d * step;
1426 }
1427 return sum;
1428 }
1429
1431 double fitSpotZwidth(const std::vector<Event>& evts, const SpotParam& spotPar, const std::vector<double>& sizesXY)
1432 {
1433
1434 std::vector<double> dataVec;
1435 std::vector<double> zzVec;
1436
1437
1438 for (const auto& e : evts) {
1439 double z0 = getCorrZ(e.mu0, e.t, spotPar);
1440 double z1 = getCorrZ(e.mu1, e.t, spotPar);
1441
1442 double corr = getZ12th(e, sizesXY);
1443 double z0z1Corr = z0 * z1 - corr;
1444
1445
1446 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1447 dataVec.push_back(z0z1Corr);
1448 zzVec.push_back(1);
1449 }
1450 }
1451
1452 MatrixXd mat = vecs2mat({zzVec});
1453
1454 std::vector<double> pars, err2;
1455 double err2Mean, err2press, err2pressErr;
1456 std::tie(pars, err2) = linearFitErr(mat, vec2vec(dataVec), err2Mean, err2press, err2pressErr);
1457
1458 return pars[0];
1459
1460 }
1461
1462
1463
1464
1466 SpotParam fitSpotPositionSplines(const std::vector<Event>& evts, const std::vector<double>& splX, const std::vector<double>& splY,
1467 const std::vector<double>& splKX, const std::vector<double>& splKY)
1468 {
1469 std::vector<std::vector<double>> basesX = fillSplineBasesZero(evts, splX, [](const Track & tr, double) {return sin(tr.phi0);});
1470 std::vector<std::vector<double>> basesY = fillSplineBasesZero(evts, splY, [](const Track & tr, double) {return -cos(tr.phi0);});
1471
1472 std::vector<std::vector<double>> basesKX = fillSplineBasesZero(evts, splKX, [](const Track & tr, double) {return sin(tr.phi0) * tr.z0;});
1473 std::vector<std::vector<double>> basesKY = fillSplineBasesZero(evts, splKY, [](const Track & tr, double) {return -cos(tr.phi0) * tr.z0;});
1474
1475
1476 std::vector<double> dataVec;
1477 for (const auto& e : evts) {
1478 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1479 dataVec.push_back(e.mu0.d0);
1480 dataVec.push_back(e.mu1.d0);
1481 }
1482 }
1483
1484 std::vector<std::vector<double>> allVecs = merge({basesX, basesY, basesKX, basesKY});
1485
1486 MatrixXd A = vecs2mat(allVecs);
1487
1488
1489 VectorXd vData = vec2vec(dataVec);
1490
1491 std::vector<double> pars(A.cols()), err2(A.cols());
1492 double err2Mean, err2press, err2pressErr;
1493 std::tie(pars, err2) = linearFitErr(A, vData, err2Mean, err2press, err2pressErr);
1494
1495 for (auto& e : err2) e = sqrt(e);
1496 return SpotParam(pars, err2, {splX, splY, splKX, splKY});
1497 }
1498
1500 SpotParam fitSpotPositionSplines(const std::vector<Event>& evts, const std::vector<double>& splX, const std::vector<double>& splY,
1501 const std::vector<double>& splKX, const std::vector<double>& splKY, const SpotParam& spotPars)
1502 {
1503 std::vector<std::vector<double>> basesX = fillSplineBasesZero(evts, splX, [](const Track & tr, double) {return sin(tr.phi0);});
1504 std::vector<std::vector<double>> basesY = fillSplineBasesZero(evts, splY, [](const Track & tr, double) {return -cos(tr.phi0);});
1505
1506 std::vector<std::vector<double>> basesKX = fillSplineBasesZero(evts, splKX, [ = ](const Track & tr, double t) {return sin(tr.phi0) * (getZIPest(tr, t, spotPars) - spotPars.z.val(t));});
1507 std::vector<std::vector<double>> basesKY = fillSplineBasesZero(evts, splKY, [ = ](const Track & tr, double t) {return -cos(tr.phi0) * (getZIPest(tr, t, spotPars) - spotPars.z.val(t));});
1508
1509
1510 std::vector<double> dataVec;
1511 for (const auto& e : evts) {
1512 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1513 dataVec.push_back(e.mu0.d0);
1514 dataVec.push_back(e.mu1.d0);
1515 }
1516 }
1517
1518 std::vector<std::vector<double>> allVecs = merge({basesX, basesY, basesKX, basesKY});
1519
1520 MatrixXd A = vecs2mat(allVecs);
1521
1522
1523 VectorXd vData = vec2vec(dataVec);
1524
1525 std::vector<double> pars(A.cols()), err2(A.cols());
1526 double err2Mean, err2press, err2pressErr;
1527 std::tie(pars, err2) = linearFitErr(A, vData, err2Mean, err2press, err2pressErr);
1528
1529 for (auto& e : err2) e = sqrt(e);
1530 auto res = SpotParam(pars, err2, {splX, splY, splKX, splKY});
1531 res.z = spotPars.z;
1532 return res;
1533 }
1534
1535
1536
1537
1538
1540 SpotParam fitSpotPositionSplines(const std::vector<Event>& evts, const std::vector<double>& splX, const std::vector<double>& splY)
1541 {
1542 std::vector<std::vector<double>> basesX = fillSplineBasesZero(evts, splX, [](const Track & tr, double) {return sin(tr.phi0);});
1543 std::vector<std::vector<double>> basesY = fillSplineBasesZero(evts, splY, [](const Track & tr, double) {return -cos(tr.phi0);});
1544
1545 std::vector<double> dataVec;
1546 for (const auto& e : evts) {
1547 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1548 dataVec.push_back(e.mu0.d0);
1549 dataVec.push_back(e.mu1.d0);
1550 }
1551 }
1552
1553 std::vector<std::vector<double>> allVecs = merge({basesX, basesY});
1554
1555 MatrixXd A = vecs2mat(allVecs);
1556
1557 VectorXd vData = vec2vec(dataVec);
1558
1559 std::vector<double> pars(A.cols()), err2(A.cols());
1560 double err2Mean, err2press, err2pressErr;
1561 std::tie(pars, err2) = linearFitErr(A, vData, err2Mean, err2press, err2pressErr);
1562
1563 for (auto& e : err2) e = sqrt(e);
1564 return SpotParam(pars, err2, {splX, splY});
1565 }
1566
1567
1568
1569
1570
1571
1573 SpotParam fitZpositionSplines(const std::vector<Event>& evts, const std::vector<double>& splX, const std::vector<double>& splY,
1574 const std::vector<double>& splKX, const std::vector<double>& splKY,
1575 const std::vector<double>& splZ)
1576 {
1577 std::vector<std::vector<double>> basesX = fillSplineBasesZero(evts, splX, [](const Track & tr, double) {return -tr.tanlambda * cos(tr.phi0);});
1578 std::vector<std::vector<double>> basesY = fillSplineBasesZero(evts, splY, [](const Track & tr, double) {return -tr.tanlambda * sin(tr.phi0);});
1579
1580 std::vector<std::vector<double>> basesKX = fillSplineBasesZero(evts, splKX, [](const Track & tr, double) {return -tr.z0 * tr.tanlambda * cos(tr.phi0);});
1581 std::vector<std::vector<double>> basesKY = fillSplineBasesZero(evts, splKY, [](const Track & tr, double) {return -tr.z0 * tr.tanlambda * sin(tr.phi0);});
1582
1583 std::vector<std::vector<double>> basesZ = fillSplineBasesZero(evts, splZ, [](Track, double) {return 1;});
1584
1585
1586 std::vector<double> dataVec;
1587 for (const auto& e : evts) {
1588 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1589 dataVec.push_back(e.mu0.z0);
1590 dataVec.push_back(e.mu1.z0);
1591 }
1592 }
1593
1594 std::vector<std::vector<double>> allVecs = merge({basesX, basesY, basesKX, basesKY, basesZ});
1595
1596 MatrixXd A = vecs2mat(allVecs);
1597
1598 VectorXd vData = vec2vec(dataVec);
1599
1600 std::vector<double> pars(A.cols()), err2(A.cols());
1601 double err2Mean, err2press, err2pressErr;
1602 std::tie(pars, err2) = linearFitErr(A, vData, err2Mean, err2press, err2pressErr);
1603
1604 for (auto& e : err2) e = sqrt(e);
1605 return SpotParam(pars, err2, {splX, splY, splKX, splKY, splZ});
1606 }
1607
1608
1609
1611 SpotParam fitZpositionSplinesSimple(const std::vector<Event>& evts, const std::vector<double>& splZ, const SpotParam& spotPars)
1612 {
1613 std::vector<std::vector<double>> basesZ = fillSplineBasesZero(evts, splZ, [](Track, double) {return 1;});
1614
1615 std::vector<double> dataVec;
1616 for (const auto& e : evts) {
1617 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1618 double z1 = getZIPest(e.mu0, e.t, spotPars);
1619 double z2 = getZIPest(e.mu1, e.t, spotPars);
1620 dataVec.push_back(z1);
1621 dataVec.push_back(z2);
1622 }
1623 }
1624
1625 MatrixXd A = vecs2mat({basesZ});
1626
1627 VectorXd vData = vec2vec(dataVec);
1628
1629 std::vector<double> pars(A.cols()), err2(A.cols());
1630 double err2Mean, err2press, err2pressErr;
1631 std::tie(pars, err2) = linearFitErr(A, vData, err2Mean, err2press, err2pressErr);
1632
1633 for (auto& e : err2) e = sqrt(e);
1634
1635 SpotParam parsUpd = spotPars;
1636 parsUpd.z.vals = pars;
1637 parsUpd.z.errs = err2;
1638 parsUpd.z.nodes = splZ;
1639
1640 return parsUpd;
1641 }
1642
1643
1644
1646 std::vector<double> fitSpotWidthCMS(const std::vector<Event>& evts, const SpotParam& spotPar)
1647 {
1648
1649 std::vector<double> dataVec, ccVec, ssVec, scVec;
1650
1651
1652 for (const auto& e : evts) {
1653 double d0 = getCorrD(e.mu0, e.t, spotPar);
1654 double d1 = getCorrD(e.mu1, e.t, spotPar);
1655
1656 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1657 dataVec.push_back(d0 * d1);
1658
1659 ccVec.push_back(cos(e.mu0.phi0)*cos(e.mu1.phi0));
1660 ssVec.push_back(sin(e.mu0.phi0)*sin(e.mu1.phi0));
1661 scVec.push_back(-(sin(e.mu0.phi0)*cos(e.mu1.phi0) + sin(e.mu1.phi0)*cos(e.mu0.phi0)));
1662 }
1663 }
1664
1665
1666 MatrixXd mat = vecs2mat({ssVec, ccVec, scVec});
1667
1668 // Linear fit with constraint on positiveness
1669 VectorXd resPhys = linearFitPos(mat, vec2vec(dataVec));
1670
1671 return {resPhys(0), resPhys(1), resPhys(2)};
1672 }
1673
1674
1676 void plotSpotSizePull(const std::vector<Event>& evts, const SpotParam& spotPar, const std::vector<double>& sizesXY)
1677 {
1678 TH1D* hPull = new TH1D(rn(), "", 100, -2000, 2000);
1679 for (const auto& e : evts) {
1680 if (!e.isSig) continue;
1681
1682 double d0 = getCorrD(e.mu0, e.t, spotPar);
1683 double d1 = getCorrD(e.mu1, e.t, spotPar);
1684
1685 double d12Th = getD12th(e, sizesXY);
1686
1687 hPull->Fill(d0 * d1 - d12Th);
1688 }
1689 TCanvas* c = new TCanvas(rn(), "");
1690 hPull->Draw();
1691 c->SaveAs("pullsSize.pdf");
1692 }
1693
1694
1696 void plotSpotSizeZPull(const std::vector<Event>& evts, const SpotParam& spotPar, const std::vector<double>& sizesXY, double sizeZZ)
1697 {
1698 TH1D* hPull = new TH1D(rn(), "", 100, -300e3, 600e3);
1699 for (const auto& e : evts) {
1700 if (!e.isSig) continue;
1701
1702 double z0 = getCorrZ(e.mu0, e.t, spotPar);
1703 double z1 = getCorrZ(e.mu1, e.t, spotPar);
1704
1705 double corr = getZ12th(e, sizesXY);
1706 double res = z0 * z1 - corr - sizeZZ;
1707
1708 hPull->Fill(res);
1709 }
1710
1711 gStyle->SetOptStat(2210);
1712 TCanvas* c = new TCanvas(rn(), "");
1713 hPull->Draw();
1714 B2INFO("zSizeFit mean " << hPull->GetMean());
1715 B2INFO("zSizeFit rms " << hPull->GetRMS());
1716
1717 c->SaveAs("pullsZSize.pdf");
1718 }
1719
1720
1721
1722
1724 void plotSpotSizeFit(const std::vector<Event>& evts, const SpotParam& par, const std::vector<double>& sizeXY)
1725 {
1726 double sxx = sizeXY[0];
1727 double syy = sizeXY[1];
1728 double sxy = sizeXY[2];
1729
1730 gStyle->SetOptStat(0);
1731
1732 TProfile* profSxx = new TProfile(rn(), "", 50, -1, 1);
1733 TProfile* profSyy = new TProfile(rn(), "", 50, -1, 1);
1734 TProfile* profSxy = new TProfile(rn(), "", 50, -1, 1);
1735 for (const auto& e : evts) {
1736 if (!e.isSig) continue;
1737
1738 double cc = cos(e.mu0.phi0) * cos(e.mu1.phi0);
1739 double ss = sin(e.mu0.phi0) * sin(e.mu1.phi0);
1740 double sc = - (sin(e.mu0.phi0) * cos(e.mu1.phi0) + sin(e.mu1.phi0) * cos(e.mu0.phi0));
1741
1742 double d0 = getCorrD(e.mu0, e.t, par);
1743 double d1 = getCorrD(e.mu1, e.t, par);
1744
1745 double data = d0 * d1;
1746
1747 profSxx->Fill(ss, data - syy * cc - sxy * sc);
1748 profSyy->Fill(cc, data - sxx * ss - sxy * sc);
1749 profSxy->Fill(sc, data - syy * cc - sxx * ss);
1750 }
1751
1752 TCanvas* c = new TCanvas(rn(), "", 1200, 500);
1753 c->Divide(3, 1);
1754 c->cd(1);
1755 profSxx->Draw();
1756 profSxx->GetXaxis()->SetTitle("sin #phi_{1} sin #phi_{2}");
1757 profSxx->GetYaxis()->SetTitle("#LTd_{1} d_{2}#GT - corr_{xx} [#mum^{2}]");
1758 TF1* fxx = new TF1(rn(), "[0]*x", -1, 1);
1759 fxx->SetParameter(0, sxx);
1760 fxx->Draw("same");
1761
1762 c->cd(2);
1763 profSyy->Draw();
1764 profSyy->GetXaxis()->SetTitle("cos #phi_{1} cos #phi_{2}");
1765 profSyy->GetYaxis()->SetTitle("#LTd_{1} d_{2}#GT - corr_{yy} [#mum^{2}]");
1766 TF1* fyy = new TF1(rn(), "[0]*x", -1, 1);
1767 fyy->SetParameter(0, syy);
1768 fyy->Draw("same");
1769
1770 c->cd(3);
1771 profSxy->Draw();
1772 profSxy->GetXaxis()->SetTitle("-(sin #phi_{1} cos #phi_{2} + sin #phi_{2} cos #phi_{1})");
1773 profSxy->GetYaxis()->SetTitle("#LTd_{1} d_{2}#GT - corr_{xy} [#mum^{2}]");
1774 TF1* fxy = new TF1(rn(), "[0]*x", -1, 1);
1775 fxy->SetParameter(0, sxy);
1776 fxy->Draw("same");
1777
1778 c->SaveAs("SizeFit.pdf");
1779 }
1780
1781
1783 void plotSpotZSizeFit(const std::vector<Event>& evts, const SpotParam& par, const std::vector<double>& sizesXY, double sizeZZ)
1784 {
1785
1786 gStyle->SetOptStat(0);
1787
1788
1789 TProfile* zzProfPhi = new TProfile(rn(), "", 100, -M_PI, M_PI);
1790 TProfile* zzProfXX = new TProfile(rn(), "", 100, -M_PI / 4, 2 * M_PI);
1791 TProfile* zzProfYY = new TProfile(rn(), "", 100, -M_PI / 4, 2 * M_PI);
1792 TProfile* zzProfXY = new TProfile(rn(), "", 100, -2 * M_PI, 2 * M_PI);
1793 TProfile* zzProfXZ = new TProfile(rn(), "", 100, -2 * M_PI, 2 * M_PI);
1794 TProfile* zzProfYZ = new TProfile(rn(), "", 100, -2 * M_PI, 2 * M_PI);
1795
1796
1797 for (const auto& e : evts) {
1798 if (e.isSig) {
1799 double z0 = getCorrZ(e.mu0, e.t, par);
1800 double z1 = getCorrZ(e.mu1, e.t, par);
1801
1802 double corr = getZ12th(e, sizesXY);
1803 double z0z1Corr = z0 * z1 - corr;
1804
1805 double xx = e.mu0.tanlambda * e.mu1.tanlambda * cos(e.mu0.phi0) * cos(e.mu1.phi0);
1806 double yy = e.mu0.tanlambda * e.mu1.tanlambda * sin(e.mu0.phi0) * sin(e.mu1.phi0);
1807 double xy = e.mu0.tanlambda * e.mu1.tanlambda * (sin(e.mu0.phi0) * cos(e.mu1.phi0) + cos(e.mu0.phi0) * sin(e.mu1.phi0));
1808 double xz = - (e.mu0.tanlambda * cos(e.mu0.phi0) + e.mu1.tanlambda * cos(e.mu1.phi0));
1809 double yz = - (e.mu0.tanlambda * sin(e.mu0.phi0) + e.mu1.tanlambda * sin(e.mu1.phi0));
1810
1811
1812 zzProfPhi->Fill(e.mu0.phi0, z0z1Corr);
1813 zzProfPhi->Fill(e.mu1.phi0, z0z1Corr);
1814 zzProfXX->Fill(xx, z0z1Corr);
1815 zzProfYY->Fill(yy, z0z1Corr);
1816 zzProfXY->Fill(xy, z0z1Corr);
1817 zzProfXZ->Fill(xz, z0z1Corr);
1818 zzProfYZ->Fill(yz, z0z1Corr);
1819 }
1820 }
1821
1822 TF1* f = new TF1(rn(), "[0]", -2 * M_PI, 2 * M_PI);
1823 f->SetParameter(0, sizeZZ);
1824
1825 TCanvas* c = new TCanvas(rn(), "", 1200, 500);
1826 c->Divide(3, 2);
1827 c->cd(1);
1828 zzProfPhi->Draw();
1829 zzProfPhi->GetXaxis()->SetTitle("#phi_{0} [rad]");
1830 zzProfPhi->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1831 f->Draw("same");
1832
1833 c->cd(2);
1834 zzProfXX->Draw();
1835 zzProfXX->GetXaxis()->SetTitle("xx sensitive");
1836 zzProfXX->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1837 f->Draw("same");
1838
1839 c->cd(3);
1840 zzProfYY->Draw();
1841 zzProfYY->GetXaxis()->SetTitle("yy sensitive");
1842 zzProfYY->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1843 f->Draw("same");
1844
1845 c->cd(4);
1846 zzProfXY->Draw();
1847 zzProfXY->GetXaxis()->SetTitle("xy sensitive");
1848 zzProfXY->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1849 f->Draw("same");
1850
1851 c->cd(5);
1852 zzProfXZ->Draw();
1853 zzProfXZ->GetXaxis()->SetTitle("xz sensitive");
1854 zzProfXZ->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1855 f->Draw("same");
1856
1857 c->cd(6);
1858 zzProfYZ->Draw();
1859 zzProfYZ->GetXaxis()->SetTitle("yz sensitive");
1860 zzProfYZ->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1861 f->Draw("same");
1862
1863 c->SaveAs("SizeZFit.pdf");
1864 }
1865
1866
1867
1869 void removeSpotSizeOutliers(std::vector<Event>& evts, const SpotParam& spotPar, const std::vector<double>& sizesXY,
1870 double cut = 1500)
1871 {
1872
1873 int nRem = 0;
1874 int nAll = 0;
1875 for (auto& e : evts) {
1876 if (!e.isSig) continue;
1877
1878 double d0 = getCorrD(e.mu0, e.t, spotPar);
1879 double d1 = getCorrD(e.mu1, e.t, spotPar);
1880 double d12Th = getD12th(e, sizesXY);
1881
1882 e.isSig = abs(d0 * d1 - d12Th) < cut;
1883 nRem += !e.isSig;
1884 ++nAll;
1885 }
1886 B2INFO("Removed fraction Size " << nRem / (nAll + 0.));
1887 }
1888
1889
1891 void removeSpotSizeZOutliers(std::vector<Event>& evts, const SpotParam& spotPar, const std::vector<double>& sizesXY, double sizeZZ,
1892 double cut = 150000)
1893 {
1894
1895 int nRem = 0;
1896 int nAll = 0;
1897 for (auto& e : evts) {
1898 if (!e.isSig) continue;
1899
1900 double z0 = getCorrZ(e.mu0, e.t, spotPar);
1901 double z1 = getCorrZ(e.mu1, e.t, spotPar);
1902
1903 double corr = getZ12th(e, sizesXY);
1904 double res = z0 * z1 - corr - sizeZZ;
1905
1906
1907 e.isSig = abs(res) < cut;
1908 nRem += !e.isSig;
1909 ++nAll;
1910 }
1911 B2INFO("Removed fraction Size " << nRem / (nAll + 0.));
1912 }
1913
1914
1916 MatrixXd toMat(TRotation rot)
1917 {
1918 MatrixXd rotM(3, 3);
1919 rotM(0, 0) = rot.XX();
1920 rotM(0, 1) = rot.XY();
1921 rotM(0, 2) = rot.XZ();
1922 rotM(1, 0) = rot.YX();
1923 rotM(1, 1) = rot.YY();
1924 rotM(1, 2) = rot.YZ();
1925 rotM(2, 0) = rot.ZX();
1926 rotM(2, 1) = rot.ZY();
1927 rotM(2, 2) = rot.ZZ();
1928
1929 return rotM;
1930 }
1931
1932
1940 MatrixXd getRotatedSizeMatrix(const std::vector<double>& xySize, double zzSize, double kX, double kY)
1941 {
1942 TRotation rot; // rotation moving eZ=(0,0,1) to (kX, kY, 1)
1943 rot.RotateX(-kY); //x-rot
1944 rot.RotateY(+kX); //y-rot
1945
1946 MatrixXd rotM = toMat(rot);
1947 MatrixXd rotMT = rotM.transpose();
1948
1949 Matrix3d eigenMat = Matrix3d::Zero(); //z-rot included in eigenMat
1950 eigenMat(0, 0) = xySize[0];
1951 eigenMat(1, 1) = xySize[1];
1952 eigenMat(0, 1) = xySize[2];
1953 eigenMat(1, 0) = xySize[2];
1954 eigenMat(2, 2) = zzSize;
1955
1956 return (rotM * eigenMat * rotMT);
1957 }
1958
1959
1960
1961
1962
1963
1964
1965 // Returns tuple with the beamspot parameters
1966 std::tuple<std::vector<VectorXd>, std::vector<MatrixXd>, MatrixXd> runBeamSpotAnalysis(std::vector<Event> evts,
1967 const std::vector<double>& splitPoints)
1968 {
1969 const double xyPosLimit = 70; //um
1970 const double xySize2Limit = 1600; //40 x 40 um^2
1971 const double zPosLimit = 1200; //um
1972
1973
1974 std::vector<double> indX = splitPoints;
1975 std::vector<double> indY = splitPoints;
1976 std::vector<double> indZ = splitPoints;
1977
1978 //no time dependence, as for beam size
1979 std::vector<double> indKX = {};
1980 std::vector<double> indKY = {};
1981
1982 UnknownPars allPars;
1983 // Only a single BootStrap replica is used and the debugging plots are disabled (kPlot = -1),
1984 // hence the conditions on k below are known to be always false.
1985 // cppcheck-suppress-begin knownConditionTrueFalse
1986 for (int k = 0; k < 1; ++k) { //loop over BootStrap replicas
1987 for (auto& e : evts) e.isSig = true; //reset cuts
1988 if (k != 0) bootStrap(evts);
1989
1990
1991 //simple XY pos fit
1992 auto resTemp = fitSpotPositionSplines(evts, indX, indY);
1993
1994 const int kPlot = -1;
1995 if (k == kPlot) {
1996 plotSpotPositionFit(evts, resTemp, "positionFitSimpe");
1997 plotSpotPositionPull(evts, resTemp, "pullsPositionSimple", xyPosLimit);
1998 }
1999 removeSpotPositionOutliers(evts, resTemp, xyPosLimit);
2000
2001 //simple XY pos fit (with outliers removed)
2002 auto resFin = fitSpotPositionSplines(evts, indX, indY);
2003 if (k == kPlot) {
2004 plotSpotPositionFit(evts, resFin, "positionFitSimpleC");
2005 plotSpotPositionPull(evts, resFin, "pullsPositionSimpleC", xyPosLimit);
2006 plotXYtimeDep(evts, resFin, "simplePosTimeDep");
2007 }
2008
2009 //Z position fit
2010 auto resZmy = fitZpositionSplinesSimple(evts, indZ, resFin);
2011 if (k == kPlot) {
2012 plotSpotZPositionFit(evts, resZmy, "positionFitSimpleZ");
2013 plotSpotZpositionPull(evts, resZmy, "zPositionPull", zPosLimit);
2014 }
2015
2016 removeSpotZpositionOutliers(evts, resZmy, zPosLimit);
2017
2018 //Z position fit (with outliers removed)
2019 resZmy = fitZpositionSplinesSimple(evts, indZ, resZmy);
2020
2021
2022 //complete XY pos fit
2023 auto resNew = fitSpotPositionSplines(evts, indX, indY, indKX, indKY, resZmy);
2024 if (k == kPlot) {
2025 plotSpotPositionFit(evts, resNew, "positionFitFull");
2026 plotKxKyFit(evts, resNew, "slopes");
2027 }
2028
2029 //Z position fit (iteration 2)
2030 resZmy = fitZpositionSplinesSimple(evts, indZ, resNew);
2031 if (k == kPlot) plotSpotZPositionFit(evts, resZmy, "positionFitSimpleZLast");
2032
2033
2034 //complete XY pos fit (iteration 2)
2035 resNew = fitSpotPositionSplines(evts, indX, indY, indKX, indKY, resZmy);
2036
2037 //XYZ pos fits (iteration 3)
2038 resZmy = fitZpositionSplinesSimple(evts, indZ, resNew);
2039 resNew = fitSpotPositionSplines(evts, indX, indY, indKX, indKY, resZmy);
2040
2041
2042 // fit of XY sizes (original + with outliers removed)
2043 auto vecXY = fitSpotWidthCMS(evts, resNew);
2044 if (k == kPlot) plotSpotSizePull(evts, resNew, vecXY);
2045 removeSpotSizeOutliers(evts, resNew, vecXY, xySize2Limit);
2046 vecXY = fitSpotWidthCMS(evts, resNew);
2047 if (k == kPlot) plotSpotSizeFit(evts, resNew, vecXY);
2048
2049
2050 // fit of Z size
2051 double sizeZZ = fitSpotZwidth(evts, resNew, vecXY);
2052
2053 if (k == kPlot) {
2054 plotSpotZSizeFit(evts, resNew, vecXY, sizeZZ);
2055 plotSpotSizeZPull(evts, resNew, vecXY, sizeZZ);
2056 }
2057
2058 //removeSpotSizeZOutliers(evts, resNew, vecXY, sizeZZ, 150000);
2059 //sizeZZ = fitSpotZwidth(evts, resNew, vecXY);
2060
2061 allPars.add(resNew, sqrtS(vecXY[0]), sqrtS(vecXY[1]), sqrtS(vecXY[2]), sqrtS(sizeZZ));
2062
2063 }
2064 // cppcheck-suppress-end knownConditionTrueFalse
2065
2066 //allPars.printStat();
2067
2068 std::vector<VectorXd> vtxPos;
2069 std::vector<MatrixXd> vtxErr;
2070 MatrixXd sizeMat;
2071
2072 allPars.getOutput(vtxPos, vtxErr, sizeMat);
2073
2074 return std::make_tuple(vtxPos, vtxErr, sizeMat);
2075 }
2076
2077}
Belle2::Track
Class that bundles various TrackFitResults.
Definition Track.h:25
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Belle2::vecs2mat
Eigen::MatrixXd vecs2mat(const std::vector< std::vector< double > > &vecs)
merge columns (from std::vectors) into ROOT matrix
Definition tools.h:72
Belle2::slice
std::vector< Atom > slice(std::vector< Atom > vec, int s, int e)
Slice the vector to contain only elements with indexes s .. e (included)
Definition Splitter.h:85
Belle2::merge
std::vector< std::vector< double > > merge(const std::vector< std::vector< std::vector< double > > > &toMerge)
merge { vector<double> a, vector<double> b} into {a, b}
Definition tools.h:41
Belle2::vec2vec
Eigen::VectorXd vec2vec(const std::vector< double > &vec)
std vector -> ROOT vector
Definition tools.h:51
Belle2::rn
TString rn()
Get random string.
Definition tools.h:38
Belle2::BeamSpotCalib::Event
Event containing two tracks.
Definition BeamSpotStandAlone.h:35
Belle2::BeamSpotCalib::SpotParam
structure containing most of the beam spot parameters
Definition BeamSpotStandAlone.cc:70
Belle2::BeamSpotCalib::SpotParam::print
void print()
Print BeamSpot parameters.
Definition BeamSpotStandAlone.cc:80
Belle2::BeamSpotCalib::SpotParam::z
Spline z
spline for BS center position as a function of time (z coordinate)
Definition BeamSpotStandAlone.cc:73
Belle2::BeamSpotCalib::SpotParam::x
Spline x
spline for BS center position as a function of time (x coordinate)
Definition BeamSpotStandAlone.cc:71
Belle2::BeamSpotCalib::SpotParam::SpotParam
SpotParam(const std::vector< double > &vals, const std::vector< double > &errs, const std::vector< std::vector< double > > &spls, int order=0)
Constructor based output of the linear regression, assuming zero-order splines vals,...
Definition BeamSpotStandAlone.cc:98
Belle2::BeamSpotCalib::SpotParam::kX
Spline kX
spline for BS angle in the xz plane as a function of time
Definition BeamSpotStandAlone.cc:74
Belle2::BeamSpotCalib::SpotParam::kY
Spline kY
spline for BS angle in the yz plane as a function of time
Definition BeamSpotStandAlone.cc:75
Belle2::BeamSpotCalib::SpotParam::y
Spline y
spline for BS center position as a function of time (y coordinate)
Definition BeamSpotStandAlone.cc:72
Belle2::BeamSpotCalib::Track
track parameters (neglecting curvature)
Definition BeamSpotStandAlone.h:27
Belle2::BeamSpotCalib::UnknowSpline
Spline with uncertainty obtained from the boot-strap replicas.
Definition BeamSpotStandAlone.cc:151
Belle2::BeamSpotCalib::UnknowSpline::add
void add(const Spline &spl)
add boot-strap replica
Definition BeamSpotStandAlone.cc:153
Belle2::BeamSpotCalib::UnknowSpline::spls
std::vector< Spline > spls
vector with replicas
Definition BeamSpotStandAlone.cc:152
Belle2::BeamSpotCalib::UnknowSpline::getMeanSigma
Spline getMeanSigma()
Get mean and 1-sigma errors of the spline values.
Definition BeamSpotStandAlone.cc:156
Belle2::BeamSpotCalib::UnknowSpline::getLimit
Spline getLimit(double v)
quantile of all points in spline, v=0.5 : median, v=0.16: lower 68% bound, v=0.84 : upper 68% bound
Definition BeamSpotStandAlone.cc:186
Belle2::BeamSpotCalib::UnknowVar
variable with uncertainty from boot-strap replicas
Definition BeamSpotStandAlone.cc:213
Belle2::BeamSpotCalib::UnknowVar::getStats
std::vector< double > getStats()
Get basic stats.
Definition BeamSpotStandAlone.cc:258
Belle2::BeamSpotCalib::UnknowVar::vars
std::vector< double > vars
vector of variable values for all replicas
Definition BeamSpotStandAlone.cc:214
Belle2::BeamSpotCalib::UnknowVar::getLimit
double getLimit(double v)
Get quantile (v=0.5 -> median, v=0.16,v=0.84 68% confidence interval)
Definition BeamSpotStandAlone.cc:239
Belle2::BeamSpotCalib::UnknowVar::add
void add(double x)
add value to the replicas
Definition BeamSpotStandAlone.cc:215
Belle2::BeamSpotCalib::UnknowVar::getSigma
double getSigma()
Get standard deviation.
Definition BeamSpotStandAlone.cc:225
Belle2::BeamSpotCalib::UnknowVar::printStat
void printStat(TString n)
Print variable of name n with stat-info.
Definition BeamSpotStandAlone.cc:250
Belle2::BeamSpotCalib::UnknownPars
structure including all variables of interest with uncertainties from boot-strap
Definition BeamSpotStandAlone.cc:333
Belle2::BeamSpotCalib::UnknownPars::kY
UnknowSpline kY
BS angle in yz plane.
Definition BeamSpotStandAlone.cc:337
Belle2::BeamSpotCalib::UnknownPars::crossAngle
UnknowVar crossAngle
derived value of the crossing angle of the HER & LER beams
Definition BeamSpotStandAlone.cc:349
Belle2::BeamSpotCalib::UnknownPars::matZZ
UnknowVar matZZ
ZZ element of BS size cov matrix.
Definition BeamSpotStandAlone.cc:353
Belle2::BeamSpotCalib::UnknownPars::sizeMin
UnknowVar sizeMin
smallest eigenvalue of the BS size cov matrix (similar to sizeY)
Definition BeamSpotStandAlone.cc:343
Belle2::BeamSpotCalib::UnknownPars::printStat
void printStat()
Print interesting statistics from boot-strap.
Definition BeamSpotStandAlone.cc:407
Belle2::BeamSpotCalib::UnknownPars::matYY
UnknowVar matYY
YY element of BS size cov matrix.
Definition BeamSpotStandAlone.cc:352
Belle2::BeamSpotCalib::UnknownPars::sizeZ
UnknowVar sizeZ
BS size in z direction.
Definition BeamSpotStandAlone.cc:348
Belle2::BeamSpotCalib::UnknownPars::z
UnknowSpline z
BS position (z coordinate)
Definition BeamSpotStandAlone.cc:338
Belle2::BeamSpotCalib::UnknownPars::sizeY
UnknowVar sizeY
BS size in y direction.
Definition BeamSpotStandAlone.cc:341
Belle2::BeamSpotCalib::UnknownPars::matXZ
UnknowVar matXZ
XZ element of BS size cov matrix.
Definition BeamSpotStandAlone.cc:355
Belle2::BeamSpotCalib::UnknownPars::getOutput
void getOutput(std::vector< VectorXd > &vtxPos, std::vector< MatrixXd > &vtxErr, MatrixXd &sizeMat)
get output in Belle2-like format
Definition BeamSpotStandAlone.cc:436
Belle2::BeamSpotCalib::UnknownPars::matXY
UnknowVar matXY
XY element of BS size cov matrix.
Definition BeamSpotStandAlone.cc:354
Belle2::BeamSpotCalib::UnknownPars::y
UnknowSpline y
BS position (y coordinate)
Definition BeamSpotStandAlone.cc:335
Belle2::BeamSpotCalib::UnknownPars::sizeX
UnknowVar sizeX
BS size in x direction.
Definition BeamSpotStandAlone.cc:340
Belle2::BeamSpotCalib::UnknownPars::setBranchVal
static void setBranchVal(TTree *T, std::vector< double > *vec, TString n)
save bootstrap variable to TTree
Definition BeamSpotStandAlone.cc:478
Belle2::BeamSpotCalib::UnknownPars::matXX
UnknowVar matXX
XX element of BS size cov matrix.
Definition BeamSpotStandAlone.cc:351
Belle2::BeamSpotCalib::UnknownPars::save2tree
void save2tree(TString fName)
save everything to TTree
Definition BeamSpotStandAlone.cc:496
Belle2::BeamSpotCalib::UnknownPars::sizeMax
UnknowVar sizeMax
middle eigenvalue of the BS size cov matrix (similar to sizeX)
Definition BeamSpotStandAlone.cc:344
Belle2::BeamSpotCalib::UnknownPars::setBranchSpline
static void setBranchSpline(TTree *T, Spline *spl, TString n)
save bootstrap spline to TTree
Definition BeamSpotStandAlone.cc:486
Belle2::BeamSpotCalib::UnknownPars::add
void add(const SpotParam &sPar, double SizeX, double SizeY, double SizeXY, double SizeZ)
add next boot-strap replica of the BS parameters
Definition BeamSpotStandAlone.cc:359
Belle2::BeamSpotCalib::UnknownPars::sizeXY
UnknowVar sizeXY
off-diagonal component of BS size cov matrix in frame where z' is aligned with z
Definition BeamSpotStandAlone.cc:342
Belle2::BeamSpotCalib::UnknownPars::x
UnknowSpline x
BS position (x coordinate)
Definition BeamSpotStandAlone.cc:334
Belle2::BeamSpotCalib::UnknownPars::xyAngle
UnknowVar xyAngle
angle of the BS in xy plane when z' is aligned with z
Definition BeamSpotStandAlone.cc:345
Belle2::BeamSpotCalib::UnknownPars::matYZ
UnknowVar matYZ
YZ element of BS size cov matrix.
Definition BeamSpotStandAlone.cc:356
Belle2::BeamSpotCalib::UnknownPars::kX
UnknowSpline kX
BS angle in xz plane.
Definition BeamSpotStandAlone.cc:336
Belle2::Spline::val
double val(double x) const
get value of spline at point x
Definition tools.h:171
Belle2::Spline::nodes
std::vector< double > nodes
vector of spline nodes
Definition tools.h:166
Belle2::Spline::errs
std::vector< double > errs
vector of spline errors
Definition tools.h:168
Belle2::Spline::vals
std::vector< double > vals
vector of spline values
Definition tools.h:167
Belle2::Spline::center
double center() const
Get center of the spline domain.
Definition tools.h:177