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(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(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 += pow(spls[i].vals[k] - s, 2);
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 += pow(x - m, 2);
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, (pow(SizeX, 2) - pow(SizeY, 2)));
277 return angle;
278 }
279
280
287 std::pair<double, double> getSizeMinMax(double SizeX, double SizeY, double SizeXY)
288 {
289 double A = pow(SizeX, 2) + pow(SizeY, 2);
290 double B = pow(SizeX, 2) * pow(SizeY, 2) - pow(SizeXY, 4);
291 double D = pow(A, 2) - 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 - pow(SizeXY, 2);
317 double D = pow(A, 2) - 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(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 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 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 (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 = ((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(), [](Event e1, 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 = pow((r(i) - dataH(i)) / (1 - Ahat), 2);
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 = 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 pars[0] = eigHigh * pow(cos(phi), 2) + s2MinLimit * pow(sin(phi), 2);
864 pars[1] = eigHigh * pow(sin(phi), 2) + s2MinLimit * pow(cos(phi), 2);
865 pars[2] = sin(phi) * cos(phi) * (eigHigh - s2MinLimit);
866
867
868 return pars;
869
870 }
871
872
873
874
876 TH1D* getResolution(TH2D* hRes)
877 {
878 TH1D* hSigmaAll = new TH1D(rn(), "", 50, -M_PI, M_PI);
879 for (int i = 1; i <= hRes->GetNbinsX(); ++i) {
880 TH1D* hProj = hRes->ProjectionY(rn(), i, i);
881 double rms = hProj->GetRMS();
882 double rmsErr = hProj->GetRMSError();
883 hSigmaAll->SetBinContent(i, rms * rms); //from cm2 to um2
884 hSigmaAll->SetBinError(i, 2 * abs(rms)*rmsErr);
885 }
886 return hSigmaAll;
887 }
888
890 TH1D* getMean(const TH2D* hRes)
891 {
892 TH1D* hMean = new TH1D(rn(), "", 50, -M_PI, M_PI);
893 for (int i = 1; i <= hRes->GetNbinsX(); ++i) {
894 TH1D* hProj = hRes->ProjectionY(rn(), i, i);
895 double mean = hProj->GetMean();
896 double meanErr = hProj->GetMeanError();
897 hMean->SetBinContent(i, mean);
898 hMean->SetBinError(i, meanErr);
899 }
900 return hMean;
901 }
902
904 double getD12th(Event e, std::vector<double> sizesXY)
905 {
906 double sxx = sizesXY[0];
907 double syy = sizesXY[1];
908 double sxy = sizesXY[2];
909
910 double cc = cos(e.mu0.phi0) * cos(e.mu1.phi0);
911 double ss = sin(e.mu0.phi0) * sin(e.mu1.phi0);
912 double sc = -(sin(e.mu0.phi0) * cos(e.mu1.phi0) + sin(e.mu1.phi0) * cos(e.mu0.phi0));
913
914 return ss * sxx + cc * syy + sc * sxy;
915 }
916
918 double getZ12th(Event e, std::vector<double> sizesXY)
919 {
920 double sxx = sizesXY[0];
921 double syy = sizesXY[1];
922
923 double corr = e.mu0.tanlambda * e.mu1.tanlambda * (sxx * cos(e.mu0.phi0) * cos(e.mu1.phi0) + syy * sin(e.mu0.phi0) * sin(
924 e.mu1.phi0) +
925 + (sin(e.mu0.phi0) * cos(e.mu1.phi0) + cos(e.mu0.phi0) * sin(e.mu1.phi0)));
926 return corr;
927 }
928
929
930
931
932
933
935 void plotSpotPositionFit(const std::vector<Event>& evts, SpotParam par, TString fName)
936 {
937 TGraph* gr = new TGraph();
938 TProfile* dProf = new TProfile(rn(), "dProf", 100, -M_PI, M_PI, "S");
939 TProfile* dProfRes = new TProfile(rn(), "dProfRes", 100, -M_PI, M_PI, "S");
940
941 for (auto e : evts) {
942 if (!e.isSig) continue;
943
944 double d1 = getDtimeConst(e.mu0, e.t, par);
945 double d2 = getDtimeConst(e.mu1, e.t, par);
946
947 gr->SetPoint(gr->GetN(), e.mu0.phi0, d1);
948 gr->SetPoint(gr->GetN(), e.mu1.phi0, d2);
949
950 dProf->Fill(e.mu0.phi0, d1);
951 dProf->Fill(e.mu1.phi0, d2);
952
953
954 double d1c = getCorrD(e.mu0, e.t, par);
955 double d2c = getCorrD(e.mu1, e.t, par);
956
957 dProfRes->Fill(e.mu0.phi0, d1c);
958 dProfRes->Fill(e.mu1.phi0, d2c);
959 }
960 TF1* f = new TF1(rn(), "[0]*sin(x) - [1]*cos(x)", -M_PI, M_PI);
961 f->SetParameters(par.x.val(par.x.center()), par.y.val(par.y.center()));
962
963 TCanvas* c = new TCanvas(rn(), "");
964 gr->Draw("a p");
965 gr->GetXaxis()->SetRangeUser(-M_PI, M_PI);
966 gr->SetMaximum(+1.3 * f->GetMaximum());
967 gr->SetMinimum(-1.3 * f->GetMaximum());
968
969 gr->GetXaxis()->SetTitle("#phi_{0} [rad]");
970 gr->GetYaxis()->SetTitle("d_{0} [#mum]");
971
972 f->Draw("same");
973 c->SaveAs(fName + "_dots.pdf");
974
975
976 TCanvas* d = new TCanvas(rn(), "");
977 gStyle->SetOptStat(0);
978 dProf->Draw();
979 dProf->GetXaxis()->SetRangeUser(-M_PI, M_PI);
980 dProf->SetMaximum(+1.3 * f->GetMaximum());
981 dProf->SetMinimum(-1.3 * f->GetMaximum());
982
983 dProf->GetXaxis()->SetTitle("#phi_{0} [rad]");
984 dProf->GetYaxis()->SetTitle("d_{0} [#mum]");
985
986 f->Draw("same");
987
988 B2INFO("Saving " << fName << " prof ");
989 d->SaveAs(fName + "_prof.pdf");
990
991
992 TCanvas* e = new TCanvas(rn(), "");
993 gStyle->SetOptStat(0);
994 dProfRes->Draw();
995 dProfRes->GetXaxis()->SetRangeUser(-M_PI, M_PI);
996
997 dProfRes->GetXaxis()->SetTitle("#phi_{0} [rad]");
998 dProfRes->GetYaxis()->SetTitle("d_{0} res [#mum]");
999
1000 TH1D* errP = new TH1D(rn(), "dErrP", 100, -M_PI, M_PI);
1001 TH1D* errM = new TH1D(rn(), "dErrM", 100, -M_PI, M_PI);
1002 for (int i = 1; i <= errP->GetNbinsX(); ++i) {
1003 errP->SetBinContent(i, dProfRes->GetBinError(i));
1004 errM->SetBinContent(i, -dProfRes->GetBinError(i));
1005 }
1006
1007 errP->Draw("hist same");
1008 errM->Draw("hist same");
1009
1010 f->SetParameters(0, 0);
1011 f->Draw("same");
1012
1013 B2INFO("Saving " << fName << " profRes ");
1014 e->SaveAs(fName + "_profRes.pdf");
1015
1016 }
1017
1018
1020 void plotSpotZPositionFit(const std::vector<Event>& evts, SpotParam par, TString fName)
1021 {
1022 TProfile* zProf = new TProfile(rn(), "dProf", 100, -M_PI, M_PI, "S");
1023 TGraph* gr = new TGraph();
1024 for (auto e : evts) {
1025 if (!e.isSig) continue;
1026
1027
1028 double z1ip = getZIPest(e.mu0, e.t, par);
1029 double z2ip = getZIPest(e.mu1, e.t, par);
1030
1031 double zipT = par.z.val(e.t);
1032 double zipM = par.z.val(par.z.center());
1033
1034 double val1 = z1ip - (zipT - zipM);
1035 double val2 = z2ip - (zipT - zipM);
1036
1037 gr->SetPoint(gr->GetN(), e.mu0.phi0, val1);
1038 gr->SetPoint(gr->GetN(), e.mu1.phi0, val2);
1039
1040 zProf->Fill(e.mu0.phi0, val1);
1041 zProf->Fill(e.mu1.phi0, val2);
1042 }
1043 TF1* f = new TF1(rn(), "[0]", -M_PI, M_PI);
1044 f->SetParameter(0, par.z.val(par.z.center()));
1045
1046 TCanvas* c = new TCanvas(rn(), "");
1047 c->SetLeftMargin(0.12);
1048 gr->Draw("a p");
1049 gr->GetXaxis()->SetRangeUser(-M_PI, M_PI);
1050 gr->SetMaximum(1000);
1051 gr->SetMinimum(-1000);
1052
1053 gr->GetXaxis()->SetTitle("#phi_{0} [rad]");
1054 gr->GetYaxis()->SetTitle("z_{IP} estimate [#mum]");
1055
1056 f->Draw("same");
1057 c->SaveAs(fName + "_dots.pdf");
1058
1059
1060 TCanvas* d = new TCanvas(rn(), "");
1061 gStyle->SetOptStat(0);
1062 d->SetLeftMargin(0.12);
1063 zProf->Draw();
1064 zProf->GetXaxis()->SetRangeUser(-M_PI, M_PI);
1065 zProf->SetMaximum(1000);
1066 zProf->SetMinimum(-1000);
1067
1068 zProf->GetXaxis()->SetTitle("#phi_{0} [rad]");
1069 zProf->GetYaxis()->SetTitle("z_{IP} estimate [#mum]");
1070
1071 f->Draw("same");
1072 d->SaveAs(fName + +"_prof.pdf");
1073
1074 }
1075
1076
1077
1078
1079
1081 void plotSpotPositionPull(const std::vector<Event>& evts, const SpotParam& par, TString fName, double cut = 70)
1082 {
1083 TH1D* hPull = new TH1D(rn(), "", 200, -200, 200);
1084
1085 for (auto& e : evts) {
1086 if (!e.isSig) continue;
1087
1088 double d0 = getCorrD(e.mu0, e.t, par);
1089 double d1 = getCorrD(e.mu1, e.t, par);
1090
1091 hPull->Fill(d0);
1092 hPull->Fill(d1);
1093 }
1094
1095 TCanvas* c = new TCanvas(rn(), "");
1096 gStyle->SetOptStat(2210);
1097 hPull->Draw();
1098
1099 hPull->GetXaxis()->SetTitle("pull [#mum]");
1100 hPull->GetYaxis()->SetTitle("#tracks");
1101
1102 TLine* l = new TLine;
1103 l->SetLineColor(kRed);
1104 l->DrawLine(-cut, 0, -cut, 500);
1105 l->DrawLine(+cut, 0, +cut, 500);
1106
1107 c->SaveAs(fName + ".pdf");
1108 }
1109
1110
1112 void plotKxKyFit(const std::vector<Event>& evts, SpotParam par, TString fName)
1113 {
1114 TProfile* profRes = new TProfile(rn(), "dProf", 100, -800, 800, "S");
1115 TProfile* profResKx = new TProfile(rn(), "dProfKx", 100, -800, 800, "S");
1116 TProfile* profResKy = new TProfile(rn(), "dProfKy", 100, -800, 800, "S");
1117
1118 SpotParam parNoKx = par;
1119 SpotParam parNoKy = par;
1120 parNoKx.kX.vals = {0};
1121 parNoKy.kY.vals = {0};
1122 parNoKx.kX.nodes = {};
1123 parNoKy.kY.nodes = {};
1124
1125 for (auto& e : evts) {
1126 if (!e.isSig) continue;
1127
1128 double zDiff1 = getZIPest(e.mu0, e.t, par) - par.z.val(e.t);
1129 double zDiff2 = getZIPest(e.mu1, e.t, par) - par.z.val(e.t);
1130
1131
1132 double d1 = getCorrD(e.mu0, e.t, par);
1133 double d2 = getCorrD(e.mu1, e.t, par);
1134
1135 double d1KX = getCorrD(e.mu0, e.t, parNoKx);
1136 double d2KX = getCorrD(e.mu1, e.t, parNoKx);
1137
1138
1139 profResKx->Fill(zDiff1 * sin(e.mu0.phi0), d1KX);
1140 profResKx->Fill(zDiff2 * sin(e.mu1.phi0), d2KX);
1141
1142 double d1KY = getCorrD(e.mu0, e.t, parNoKy);
1143 double d2KY = getCorrD(e.mu1, e.t, parNoKy);
1144 profResKy->Fill(-zDiff1 * cos(e.mu0.phi0), d1KY);
1145 profResKy->Fill(-zDiff2 * cos(e.mu1.phi0), d2KY);
1146
1147 profRes->Fill(zDiff1, d1);
1148 profRes->Fill(zDiff2, d2);
1149
1150
1151 }
1152
1153 TCanvas* cX = new TCanvas(rn(), "");
1154 gStyle->SetOptStat(0);
1155 profResKx->Draw();
1156
1157 profResKx->GetXaxis()->SetTitle("(z_{IP} - z_{IP}^{0}) sin #phi_{0} [#mum]");
1158 profResKx->GetYaxis()->SetTitle("d_{0} res [#mum]");
1159
1160 TF1* f = new TF1(rn(), "[0]*x", -800, 800);
1161 f->SetParameter(0, par.kX.val(par.kX.center()));
1162 f->Draw("same");
1163
1164 cX->SaveAs(fName + "_kX.pdf");
1165
1166 TCanvas* cY = new TCanvas(rn(), "");
1167 gStyle->SetOptStat(0);
1168 profResKy->Draw();
1169
1170 profResKy->GetXaxis()->SetTitle("-(z_{IP} - z_{IP}^{0}) cos #phi_{0} [#mum]");
1171 profResKy->GetYaxis()->SetTitle("d_{0} res [#mum]");
1172
1173 f->SetParameter(0, par.kY.val(par.kY.center()));
1174 f->Draw("same");
1175
1176 cY->SaveAs(fName + "_kY.pdf");
1177
1178
1179 }
1180
1182 void plotXYtimeDep(const std::vector<Event>& evts, SpotParam par, TString fName)
1183 {
1184 TProfile* profRes = new TProfile(rn(), "dProf", 50, -0.5, 0.5);
1185 TProfile* profResTx = new TProfile(rn(), "dProfTx", 50, -0.5, 0.5);
1186 TProfile* profResTy = new TProfile(rn(), "dProfTy", 50, -0.5, 0.5);
1187
1188 SpotParam parNoTx = par;
1189 SpotParam parNoTy = par;
1190 parNoTx.x.nodes = {};
1191 parNoTx.x.vals = {par.x.val(par.x.center())};
1192 parNoTy.y.nodes = {};
1193 parNoTy.y.vals = {par.y.val(par.y.center())};
1194
1195 for (auto& e : evts) {
1196 if (!e.isSig) continue;
1197
1198 double tDiff = (e.t - par.x.val(par.x.center()));
1199 //double tDiff = e.t;
1200
1201
1202 double d1 = getCorrD(e.mu0, e.t, par);
1203 double d2 = getCorrD(e.mu1, e.t, par);
1204
1205 double d1TX = getCorrD(e.mu0, e.t, parNoTx);
1206 double d2TX = getCorrD(e.mu1, e.t, parNoTx);
1207
1208 profResTx->Fill(tDiff * sin(e.mu0.phi0), d1TX);
1209 profResTx->Fill(tDiff * sin(e.mu1.phi0), d2TX);
1210
1211 double d1TY = getCorrD(e.mu0, e.t, parNoTy);
1212 double d2TY = getCorrD(e.mu1, e.t, parNoTy);
1213
1214 profResTy->Fill(-tDiff * cos(e.mu0.phi0), d1TY);
1215 profResTy->Fill(-tDiff * cos(e.mu1.phi0), d2TY);
1216
1217
1218 profRes->Fill(tDiff * sin(e.mu0.phi0), d1);
1219 profRes->Fill(tDiff * sin(e.mu1.phi0), d2);
1220
1221
1222 }
1223
1224 TCanvas* cX = new TCanvas(rn(), "");
1225 gStyle->SetOptStat(0);
1226 profResTx->Draw();
1227
1228
1229 TF1* f = new TF1(rn(), "[0]*x", -1, 1);
1230 f->SetParameter(0, (par.x.val(1) - par.x.val(0)));
1231 f->Draw("same");
1232 B2INFO("Table value " << par.x.val(1) - par.x.val(0));
1233
1234 cX->SaveAs(fName + "_tX.pdf");
1235
1236
1237 }
1238
1239
1240
1241
1242
1244 void plotSpotZpositionPull(const std::vector<Event>& evts, const SpotParam& par, TString fName, double cut = 1000)
1245 {
1246 TH1D* hPull = new TH1D(rn(), "", 200, -2000, 2000);
1247
1248 for (auto& e : evts) {
1249 if (!e.isSig) continue;
1250
1251 double z0 = getCorrZ(e.mu0, e.t, par);
1252 double z1 = getCorrZ(e.mu1, e.t, par);
1253
1254 hPull->Fill(z0);
1255 hPull->Fill(z1);
1256 }
1257
1258 gStyle->SetOptStat(2210);
1259 TCanvas* c = new TCanvas(rn(), "");
1260 hPull->Draw();
1261
1262 hPull->GetXaxis()->SetTitle("pull [#mum]");
1263 hPull->GetYaxis()->SetTitle("#tracks");
1264
1265 TLine* l = new TLine;
1266 l->SetLineColor(kRed);
1267 l->DrawLine(-cut, 0, -cut, 500);
1268 l->DrawLine(+cut, 0, +cut, 500);
1269
1270 c->SaveAs(fName + ".pdf");
1271 }
1272
1273
1275 void removeSpotPositionOutliers(std::vector<Event>& evts, const SpotParam& par, double cut = 70)
1276 {
1277 int nRem = 0;
1278 int nAll = 0;
1279 for (auto& e : evts) {
1280 if (!e.isSig) continue;
1281
1282 double d0 = getCorrD(e.mu0, e.t, par);
1283 double d1 = getCorrD(e.mu1, e.t, par);
1284
1285 e.isSig = abs(d0) < cut && abs(d1) < cut;
1286 nRem += !e.isSig;
1287 ++nAll;
1288 }
1289 B2INFO("Removed fraction Position " << nRem / (nAll + 0.));
1290 }
1291
1292
1294 void removeSpotZpositionOutliers(std::vector<Event>& evts, const SpotParam& par, double cut = 1000)
1295 {
1296 int nRem = 0;
1297 int nAll = 0;
1298 for (auto& e : evts) {
1299 if (!e.isSig) continue;
1300
1301 double z0 = getCorrZ(e.mu0, e.t, par);
1302 double z1 = getCorrZ(e.mu1, e.t, par);
1303
1304 e.isSig = abs(z0) < cut && abs(z1) < cut;
1305 nRem += !e.isSig;
1306 ++nAll;
1307 }
1308 B2INFO("Removed fraction Position " << nRem / (nAll + 0.));
1309 }
1310
1311
1312
1314 std::vector<std::vector<double>> fillSplineBasesLinear(const std::vector<Event>& evts, std::vector<double> spl,
1315 std::function<double(Track, double)> fun)
1316 {
1317 int n = spl.size(); //number of params
1318 if (n == 0 || (n == 2 && spl[0] > spl[1]))
1319 n = 1;
1320
1321 std::vector<std::vector<double>> vecs(n);
1322
1323 if (n == 1) { //no time dependence
1324 for (const auto& e : evts) {
1325 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1326 vecs[0].push_back(1 * fun(e.mu0, e.t));
1327 vecs[0].push_back(1 * fun(e.mu1, e.t));
1328 }
1329 }
1330 } else {
1331 for (int k = 0; k < n; ++k) {
1332 double xCnt = spl[k];
1333 double xLow = (k == 0) ? spl[0] : spl[k - 1];
1334 double xHigh = (k == n - 1) ? spl[n - 1] : spl[k + 1];
1335
1336 for (const auto& e : evts) {
1337 double x = e.t;
1338 double v = 0;
1339 if (xLow <= x && x < xCnt)
1340 v = (x - xLow) / (xCnt - xLow);
1341 else if (xCnt < x && x <= xHigh)
1342 v = (xHigh - x) / (xHigh - xCnt);
1343
1344
1345 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1346 vecs[k].push_back(v * fun(e.mu0, e.t));
1347 vecs[k].push_back(v * fun(e.mu1, e.t));
1348 }
1349 }
1350 }
1351 }
1352
1353 return vecs;
1354 }
1355
1356
1358 std::vector<std::vector<double>> fillSplineBasesZero(const std::vector<Event>& evts, std::vector<double> spl,
1359 std::function<double(Track, double)> fun)
1360 {
1361 int n = spl.size() + 1; //number of params
1362
1363 std::vector<std::vector<double>> vecs(n);
1364
1365 if (n == 1) { //no time dependence
1366 for (const auto& e : evts) {
1367 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1368 vecs[0].push_back(1 * fun(e.mu0, e.t));
1369 vecs[0].push_back(1 * fun(e.mu1, e.t));
1370 }
1371 }
1372 } else { // at least two parameters
1373 for (int k = 0; k < n; ++k) { //loop over spline parameters
1374 double xLow = -1e30;
1375 double xHigh = +1e30;
1376
1377 if (k == 0) {
1378 xHigh = spl[0];
1379 } else if (k == n - 1) {
1380 xLow = spl.back();
1381 } else {
1382 xLow = spl[k - 1];
1383 xHigh = spl[k];
1384 }
1385
1386 for (const auto& e : evts) {
1387 double x = e.t;
1388 double v = 0;
1389 if (xLow <= x && x < xHigh)
1390 v = 1;
1391
1392 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1393 vecs[k].push_back(v * fun(e.mu0, e.t));
1394 vecs[k].push_back(v * fun(e.mu1, e.t));
1395 }
1396 }
1397 }
1398 }
1399
1400 return vecs;
1401 }
1402
1403
1404
1405
1406
1407
1409 double compareSplines(const Spline& spl1, const Spline& spl2)
1410 {
1411 double sum = 0;
1412
1413 double step = 0.001;
1414 for (double x = 0; x <= 1 + step / 2; x += step) {
1415 double v1 = spl1.val(x);
1416 double e1 = spl1.err(x);
1417 double v2 = spl2.val(x);
1418 double e2 = spl2.err(x);
1419
1420 double d = pow(v2 - v1, 2) / pow(std::max(e1, e2), 2);
1421 sum += d * step;
1422 }
1423 return sum;
1424 }
1425
1427 double fitSpotZwidth(const std::vector<Event>& evts, const SpotParam& spotPar, const std::vector<double>& sizesXY)
1428 {
1429
1430 std::vector<double> dataVec;
1431 std::vector<double> zzVec;
1432
1433
1434 for (auto e : evts) {
1435 double z0 = getCorrZ(e.mu0, e.t, spotPar);
1436 double z1 = getCorrZ(e.mu1, e.t, spotPar);
1437
1438 double corr = getZ12th(e, sizesXY);
1439 double z0z1Corr = z0 * z1 - corr;
1440
1441
1442 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1443 dataVec.push_back(z0z1Corr);
1444 zzVec.push_back(1);
1445 }
1446 }
1447
1448 MatrixXd mat = vecs2mat({zzVec});
1449
1450 std::vector<double> pars, err2;
1451 double err2Mean, err2press, err2pressErr;
1452 std::tie(pars, err2) = linearFitErr(mat, vec2vec(dataVec), err2Mean, err2press, err2pressErr);
1453
1454 return pars[0];
1455
1456 }
1457
1458
1459
1460
1462 SpotParam fitSpotPositionSplines(const std::vector<Event>& evts, const std::vector<double>& splX, const std::vector<double>& splY,
1463 const std::vector<double>& splKX, const std::vector<double>& splKY)
1464 {
1465 std::vector<std::vector<double>> basesX = fillSplineBasesZero(evts, splX, [](Track tr, double) {return sin(tr.phi0);});
1466 std::vector<std::vector<double>> basesY = fillSplineBasesZero(evts, splY, [](Track tr, double) {return -cos(tr.phi0);});
1467
1468 std::vector<std::vector<double>> basesKX = fillSplineBasesZero(evts, splKX, [](Track tr, double) {return sin(tr.phi0) * tr.z0;});
1469 std::vector<std::vector<double>> basesKY = fillSplineBasesZero(evts, splKY, [](Track tr, double) {return -cos(tr.phi0) * tr.z0;});
1470
1471
1472 std::vector<double> dataVec;
1473 for (auto e : evts) {
1474 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1475 dataVec.push_back(e.mu0.d0);
1476 dataVec.push_back(e.mu1.d0);
1477 }
1478 }
1479
1480 std::vector<std::vector<double>> allVecs = merge({basesX, basesY, basesKX, basesKY});
1481
1482 MatrixXd A = vecs2mat(allVecs);
1483
1484
1485 VectorXd vData = vec2vec(dataVec);
1486
1487 std::vector<double> pars(A.cols()), err2(A.cols());
1488 double err2Mean, err2press, err2pressErr;
1489 std::tie(pars, err2) = linearFitErr(A, vData, err2Mean, err2press, err2pressErr);
1490
1491 for (auto& e : err2) e = sqrt(e);
1492 return SpotParam(pars, err2, {splX, splY, splKX, splKY});
1493 }
1494
1496 SpotParam fitSpotPositionSplines(const std::vector<Event>& evts, const std::vector<double>& splX, const std::vector<double>& splY,
1497 const std::vector<double>& splKX, const std::vector<double>& splKY, const SpotParam& spotPars)
1498 {
1499 std::vector<std::vector<double>> basesX = fillSplineBasesZero(evts, splX, [](Track tr, double) {return sin(tr.phi0);});
1500 std::vector<std::vector<double>> basesY = fillSplineBasesZero(evts, splY, [](Track tr, double) {return -cos(tr.phi0);});
1501
1502 std::vector<std::vector<double>> basesKX = fillSplineBasesZero(evts, splKX, [ = ](Track tr, double t) {return sin(tr.phi0) * (getZIPest(tr, t, spotPars) - spotPars.z.val(t));});
1503 std::vector<std::vector<double>> basesKY = fillSplineBasesZero(evts, splKY, [ = ](Track tr, double t) {return -cos(tr.phi0) * (getZIPest(tr, t, spotPars) - spotPars.z.val(t));});
1504
1505
1506 std::vector<double> dataVec;
1507 for (auto e : evts) {
1508 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1509 dataVec.push_back(e.mu0.d0);
1510 dataVec.push_back(e.mu1.d0);
1511 }
1512 }
1513
1514 std::vector<std::vector<double>> allVecs = merge({basesX, basesY, basesKX, basesKY});
1515
1516 MatrixXd A = vecs2mat(allVecs);
1517
1518
1519 VectorXd vData = vec2vec(dataVec);
1520
1521 std::vector<double> pars(A.cols()), err2(A.cols());
1522 double err2Mean, err2press, err2pressErr;
1523 std::tie(pars, err2) = linearFitErr(A, vData, err2Mean, err2press, err2pressErr);
1524
1525 for (auto& e : err2) e = sqrt(e);
1526 auto res = SpotParam(pars, err2, {splX, splY, splKX, splKY});
1527 res.z = spotPars.z;
1528 return res;
1529 }
1530
1531
1532
1533
1534
1536 SpotParam fitSpotPositionSplines(const std::vector<Event>& evts, const std::vector<double>& splX, const std::vector<double>& splY)
1537 {
1538 std::vector<std::vector<double>> basesX = fillSplineBasesZero(evts, splX, [](Track tr, double) {return sin(tr.phi0);});
1539 std::vector<std::vector<double>> basesY = fillSplineBasesZero(evts, splY, [](Track tr, double) {return -cos(tr.phi0);});
1540
1541 std::vector<double> dataVec;
1542 for (auto e : evts) {
1543 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1544 dataVec.push_back(e.mu0.d0);
1545 dataVec.push_back(e.mu1.d0);
1546 }
1547 }
1548
1549 std::vector<std::vector<double>> allVecs = merge({basesX, basesY});
1550
1551 MatrixXd A = vecs2mat(allVecs);
1552
1553 VectorXd vData = vec2vec(dataVec);
1554
1555 std::vector<double> pars(A.cols()), err2(A.cols());
1556 double err2Mean, err2press, err2pressErr;
1557 std::tie(pars, err2) = linearFitErr(A, vData, err2Mean, err2press, err2pressErr);
1558
1559 for (auto& e : err2) e = sqrt(e);
1560 return SpotParam(pars, err2, {splX, splY});
1561 }
1562
1563
1564
1565
1566
1567
1569 SpotParam fitZpositionSplines(const std::vector<Event>& evts, const std::vector<double>& splX, const std::vector<double>& splY,
1570 const std::vector<double>& splKX, const std::vector<double>& splKY,
1571 const std::vector<double>& splZ)
1572 {
1573 std::vector<std::vector<double>> basesX = fillSplineBasesZero(evts, splX, [](Track tr, double) {return -tr.tanlambda * cos(tr.phi0);});
1574 std::vector<std::vector<double>> basesY = fillSplineBasesZero(evts, splY, [](Track tr, double) {return -tr.tanlambda * sin(tr.phi0);});
1575
1576 std::vector<std::vector<double>> basesKX = fillSplineBasesZero(evts, splKX, [](Track tr, double) {return -tr.z0 * tr.tanlambda * cos(tr.phi0);});
1577 std::vector<std::vector<double>> basesKY = fillSplineBasesZero(evts, splKY, [](Track tr, double) {return -tr.z0 * tr.tanlambda * sin(tr.phi0);});
1578
1579 std::vector<std::vector<double>> basesZ = fillSplineBasesZero(evts, splZ, [](Track, double) {return 1;});
1580
1581
1582 std::vector<double> dataVec;
1583 for (auto e : evts) {
1584 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1585 dataVec.push_back(e.mu0.z0);
1586 dataVec.push_back(e.mu1.z0);
1587 }
1588 }
1589
1590 std::vector<std::vector<double>> allVecs = merge({basesX, basesY, basesKX, basesKY, basesZ});
1591
1592 MatrixXd A = vecs2mat(allVecs);
1593
1594 VectorXd vData = vec2vec(dataVec);
1595
1596 std::vector<double> pars(A.cols()), err2(A.cols());
1597 double err2Mean, err2press, err2pressErr;
1598 std::tie(pars, err2) = linearFitErr(A, vData, err2Mean, err2press, err2pressErr);
1599
1600 for (auto& e : err2) e = sqrt(e);
1601 return SpotParam(pars, err2, {splX, splY, splKX, splKY, splZ});
1602 }
1603
1604
1605
1607 SpotParam fitZpositionSplinesSimple(const std::vector<Event>& evts, const std::vector<double>& splZ, const SpotParam& spotPars)
1608 {
1609 std::vector<std::vector<double>> basesZ = fillSplineBasesZero(evts, splZ, [](Track, double) {return 1;});
1610
1611 std::vector<double> dataVec;
1612 for (auto e : evts) {
1613 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1614 double z1 = getZIPest(e.mu0, e.t, spotPars);
1615 double z2 = getZIPest(e.mu1, e.t, spotPars);
1616 dataVec.push_back(z1);
1617 dataVec.push_back(z2);
1618 }
1619 }
1620
1621 MatrixXd A = vecs2mat({basesZ});
1622
1623 VectorXd vData = vec2vec(dataVec);
1624
1625 std::vector<double> pars(A.cols()), err2(A.cols());
1626 double err2Mean, err2press, err2pressErr;
1627 std::tie(pars, err2) = linearFitErr(A, vData, err2Mean, err2press, err2pressErr);
1628
1629 for (auto& e : err2) e = sqrt(e);
1630
1631 SpotParam parsUpd = spotPars;
1632 parsUpd.z.vals = pars;
1633 parsUpd.z.errs = err2;
1634 parsUpd.z.nodes = splZ;
1635
1636 return parsUpd;
1637 }
1638
1639
1640
1642 std::vector<double> fitSpotWidthCMS(const std::vector<Event>& evts, const SpotParam& spotPar)
1643 {
1644
1645 std::vector<double> dataVec, ccVec, ssVec, scVec;
1646
1647
1648 for (auto e : evts) {
1649 double d0 = getCorrD(e.mu0, e.t, spotPar);
1650 double d1 = getCorrD(e.mu1, e.t, spotPar);
1651
1652 for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1653 dataVec.push_back(d0 * d1);
1654
1655 ccVec.push_back(cos(e.mu0.phi0)*cos(e.mu1.phi0));
1656 ssVec.push_back(sin(e.mu0.phi0)*sin(e.mu1.phi0));
1657 scVec.push_back(-(sin(e.mu0.phi0)*cos(e.mu1.phi0) + sin(e.mu1.phi0)*cos(e.mu0.phi0)));
1658 }
1659 }
1660
1661
1662 MatrixXd mat = vecs2mat({ssVec, ccVec, scVec});
1663
1664 // Linear fit with constraint on positiveness
1665 VectorXd resPhys = linearFitPos(mat, vec2vec(dataVec));
1666
1667 return {resPhys(0), resPhys(1), resPhys(2)};
1668 }
1669
1670
1672 void plotSpotSizePull(const std::vector<Event>& evts, const SpotParam& spotPar, const std::vector<double>& sizesXY)
1673 {
1674 TH1D* hPull = new TH1D(rn(), "", 100, -2000, 2000);
1675 for (auto& e : evts) {
1676 if (!e.isSig) continue;
1677
1678 double d0 = getCorrD(e.mu0, e.t, spotPar);
1679 double d1 = getCorrD(e.mu1, e.t, spotPar);
1680
1681 double d12Th = getD12th(e, sizesXY);
1682
1683 hPull->Fill(d0 * d1 - d12Th);
1684 }
1685 TCanvas* c = new TCanvas(rn(), "");
1686 hPull->Draw();
1687 c->SaveAs("pullsSize.pdf");
1688 }
1689
1690
1692 void plotSpotSizeZPull(const std::vector<Event>& evts, const SpotParam& spotPar, const std::vector<double>& sizesXY, double sizeZZ)
1693 {
1694 TH1D* hPull = new TH1D(rn(), "", 100, -300e3, 600e3);
1695 for (auto& e : evts) {
1696 if (!e.isSig) continue;
1697
1698 double z0 = getCorrZ(e.mu0, e.t, spotPar);
1699 double z1 = getCorrZ(e.mu1, e.t, spotPar);
1700
1701 double corr = getZ12th(e, sizesXY);
1702 double res = z0 * z1 - corr - sizeZZ;
1703
1704 hPull->Fill(res);
1705 }
1706
1707 gStyle->SetOptStat(2210);
1708 TCanvas* c = new TCanvas(rn(), "");
1709 hPull->Draw();
1710 B2INFO("zSizeFit mean " << hPull->GetMean());
1711 B2INFO("zSizeFit rms " << hPull->GetRMS());
1712
1713 c->SaveAs("pullsZSize.pdf");
1714 }
1715
1716
1717
1718
1720 void plotSpotSizeFit(const std::vector<Event>& evts, const SpotParam& par, const std::vector<double>& sizeXY)
1721 {
1722 double sxx = sizeXY[0];
1723 double syy = sizeXY[1];
1724 double sxy = sizeXY[2];
1725
1726 gStyle->SetOptStat(0);
1727
1728 TProfile* profSxx = new TProfile(rn(), "", 50, -1, 1);
1729 TProfile* profSyy = new TProfile(rn(), "", 50, -1, 1);
1730 TProfile* profSxy = new TProfile(rn(), "", 50, -1, 1);
1731 for (auto e : evts) {
1732 if (!e.isSig) continue;
1733
1734 double cc = cos(e.mu0.phi0) * cos(e.mu1.phi0);
1735 double ss = sin(e.mu0.phi0) * sin(e.mu1.phi0);
1736 double sc = - (sin(e.mu0.phi0) * cos(e.mu1.phi0) + sin(e.mu1.phi0) * cos(e.mu0.phi0));
1737
1738 double d0 = getCorrD(e.mu0, e.t, par);
1739 double d1 = getCorrD(e.mu1, e.t, par);
1740
1741 double data = d0 * d1;
1742
1743 profSxx->Fill(ss, data - syy * cc - sxy * sc);
1744 profSyy->Fill(cc, data - sxx * ss - sxy * sc);
1745 profSxy->Fill(sc, data - syy * cc - sxx * ss);
1746 }
1747
1748 TCanvas* c = new TCanvas(rn(), "", 1200, 500);
1749 c->Divide(3, 1);
1750 c->cd(1);
1751 profSxx->Draw();
1752 profSxx->GetXaxis()->SetTitle("sin #phi_{1} sin #phi_{2}");
1753 profSxx->GetYaxis()->SetTitle("#LTd_{1} d_{2}#GT - corr_{xx} [#mum^{2}]");
1754 TF1* fxx = new TF1(rn(), "[0]*x", -1, 1);
1755 fxx->SetParameter(0, sxx);
1756 fxx->Draw("same");
1757
1758 c->cd(2);
1759 profSyy->Draw();
1760 profSyy->GetXaxis()->SetTitle("cos #phi_{1} cos #phi_{2}");
1761 profSyy->GetYaxis()->SetTitle("#LTd_{1} d_{2}#GT - corr_{yy} [#mum^{2}]");
1762 TF1* fyy = new TF1(rn(), "[0]*x", -1, 1);
1763 fyy->SetParameter(0, syy);
1764 fyy->Draw("same");
1765
1766 c->cd(3);
1767 profSxy->Draw();
1768 profSxy->GetXaxis()->SetTitle("-(sin #phi_{1} cos #phi_{2} + sin #phi_{2} cos #phi_{1})");
1769 profSxy->GetYaxis()->SetTitle("#LTd_{1} d_{2}#GT - corr_{xy} [#mum^{2}]");
1770 TF1* fxy = new TF1(rn(), "[0]*x", -1, 1);
1771 fxy->SetParameter(0, sxy);
1772 fxy->Draw("same");
1773
1774 c->SaveAs("SizeFit.pdf");
1775 }
1776
1777
1779 void plotSpotZSizeFit(const std::vector<Event>& evts, const SpotParam& par, const std::vector<double>& sizesXY, double sizeZZ)
1780 {
1781
1782 gStyle->SetOptStat(0);
1783
1784
1785 TProfile* zzProfPhi = new TProfile(rn(), "", 100, -M_PI, M_PI);
1786 TProfile* zzProfXX = new TProfile(rn(), "", 100, -M_PI / 4, 2 * M_PI);
1787 TProfile* zzProfYY = new TProfile(rn(), "", 100, -M_PI / 4, 2 * M_PI);
1788 TProfile* zzProfXY = new TProfile(rn(), "", 100, -2 * M_PI, 2 * M_PI);
1789 TProfile* zzProfXZ = new TProfile(rn(), "", 100, -2 * M_PI, 2 * M_PI);
1790 TProfile* zzProfYZ = new TProfile(rn(), "", 100, -2 * M_PI, 2 * M_PI);
1791
1792
1793 for (auto e : evts) {
1794 double z0 = getCorrZ(e.mu0, e.t, par);
1795 double z1 = getCorrZ(e.mu1, e.t, par);
1796
1797 double corr = getZ12th(e, sizesXY);
1798 double z0z1Corr = z0 * z1 - corr;
1799
1800 if (e.isSig) {
1801
1802 double xx = e.mu0.tanlambda * e.mu1.tanlambda * cos(e.mu0.phi0) * cos(e.mu1.phi0);
1803 double yy = e.mu0.tanlambda * e.mu1.tanlambda * sin(e.mu0.phi0) * sin(e.mu1.phi0);
1804 double xy = e.mu0.tanlambda * e.mu1.tanlambda * (sin(e.mu0.phi0) * cos(e.mu1.phi0) + cos(e.mu0.phi0) * sin(e.mu1.phi0));
1805 double xz = - (e.mu0.tanlambda * cos(e.mu0.phi0) + e.mu1.tanlambda * cos(e.mu1.phi0));
1806 double yz = - (e.mu0.tanlambda * sin(e.mu0.phi0) + e.mu1.tanlambda * sin(e.mu1.phi0));
1807
1808
1809 zzProfPhi->Fill(e.mu0.phi0, z0z1Corr);
1810 zzProfPhi->Fill(e.mu1.phi0, z0z1Corr);
1811 zzProfXX->Fill(xx, z0z1Corr);
1812 zzProfYY->Fill(yy, z0z1Corr);
1813 zzProfXY->Fill(xy, z0z1Corr);
1814 zzProfXZ->Fill(xz, z0z1Corr);
1815 zzProfYZ->Fill(yz, z0z1Corr);
1816 }
1817 }
1818
1819 TF1* f = new TF1(rn(), "[0]", -2 * M_PI, 2 * M_PI);
1820 f->SetParameter(0, sizeZZ);
1821
1822 TCanvas* c = new TCanvas(rn(), "", 1200, 500);
1823 c->Divide(3, 2);
1824 c->cd(1);
1825 zzProfPhi->Draw();
1826 zzProfPhi->GetXaxis()->SetTitle("#phi_{0} [rad]");
1827 zzProfPhi->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1828 f->Draw("same");
1829
1830 c->cd(2);
1831 zzProfXX->Draw();
1832 zzProfXX->GetXaxis()->SetTitle("xx sensitive");
1833 zzProfXX->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1834 f->Draw("same");
1835
1836 c->cd(3);
1837 zzProfYY->Draw();
1838 zzProfYY->GetXaxis()->SetTitle("yy sensitive");
1839 zzProfYY->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1840 f->Draw("same");
1841
1842 c->cd(4);
1843 zzProfXY->Draw();
1844 zzProfXY->GetXaxis()->SetTitle("xy sensitive");
1845 zzProfXY->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1846 f->Draw("same");
1847
1848 c->cd(5);
1849 zzProfXZ->Draw();
1850 zzProfXZ->GetXaxis()->SetTitle("xz sensitive");
1851 zzProfXZ->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1852 f->Draw("same");
1853
1854 c->cd(6);
1855 zzProfYZ->Draw();
1856 zzProfYZ->GetXaxis()->SetTitle("yz sensitive");
1857 zzProfYZ->GetYaxis()->SetTitle("#LTz_{1} z_{2}#GT - corr [#mum^{2}]");
1858 f->Draw("same");
1859
1860 c->SaveAs("SizeZFit.pdf");
1861 }
1862
1863
1864
1866 void removeSpotSizeOutliers(std::vector<Event>& evts, const SpotParam& spotPar, const std::vector<double>& sizesXY,
1867 double cut = 1500)
1868 {
1869
1870 int nRem = 0;
1871 int nAll = 0;
1872 for (auto& e : evts) {
1873 if (!e.isSig) continue;
1874
1875 double d0 = getCorrD(e.mu0, e.t, spotPar);
1876 double d1 = getCorrD(e.mu1, e.t, spotPar);
1877 double d12Th = getD12th(e, sizesXY);
1878
1879 e.isSig = abs(d0 * d1 - d12Th) < cut;
1880 nRem += !e.isSig;
1881 ++nAll;
1882 }
1883 B2INFO("Removed fraction Size " << nRem / (nAll + 0.));
1884 }
1885
1886
1888 void removeSpotSizeZOutliers(std::vector<Event>& evts, const SpotParam& spotPar, const std::vector<double>& sizesXY, double sizeZZ,
1889 double cut = 150000)
1890 {
1891
1892 int nRem = 0;
1893 int nAll = 0;
1894 for (auto& e : evts) {
1895 if (!e.isSig) continue;
1896
1897 double z0 = getCorrZ(e.mu0, e.t, spotPar);
1898 double z1 = getCorrZ(e.mu1, e.t, spotPar);
1899
1900 double corr = getZ12th(e, sizesXY);
1901 double res = z0 * z1 - corr - sizeZZ;
1902
1903
1904 e.isSig = abs(res) < cut;
1905 nRem += !e.isSig;
1906 ++nAll;
1907 }
1908 B2INFO("Removed fraction Size " << nRem / (nAll + 0.));
1909 }
1910
1911
1913 MatrixXd toMat(TRotation rot)
1914 {
1915 MatrixXd rotM(3, 3);
1916 rotM(0, 0) = rot.XX();
1917 rotM(0, 1) = rot.XY();
1918 rotM(0, 2) = rot.XZ();
1919 rotM(1, 0) = rot.YX();
1920 rotM(1, 1) = rot.YY();
1921 rotM(1, 2) = rot.YZ();
1922 rotM(2, 0) = rot.ZX();
1923 rotM(2, 1) = rot.ZY();
1924 rotM(2, 2) = rot.ZZ();
1925
1926 return rotM;
1927 }
1928
1929
1937 MatrixXd getRotatedSizeMatrix(std::vector<double> xySize, double zzSize, double kX, double kY)
1938 {
1939 TRotation rot; // rotation moving eZ=(0,0,1) to (kX, kY, 1)
1940 rot.RotateX(-kY); //x-rot
1941 rot.RotateY(+kX); //y-rot
1942
1943 MatrixXd rotM = toMat(rot);
1944 MatrixXd rotMT = rotM.transpose();
1945
1946 Matrix3d eigenMat = Matrix3d::Zero(); //z-rot included in eigenMat
1947 eigenMat(0, 0) = xySize[0];
1948 eigenMat(1, 1) = xySize[1];
1949 eigenMat(0, 1) = xySize[2];
1950 eigenMat(1, 0) = xySize[2];
1951 eigenMat(2, 2) = zzSize;
1952
1953 return (rotM * eigenMat * rotMT);
1954 }
1955
1956
1957
1958
1959
1960
1961
1962 // Returns tuple with the beamspot parameters
1963 std::tuple<std::vector<VectorXd>, std::vector<MatrixXd>, MatrixXd> runBeamSpotAnalysis(std::vector<Event> evts,
1964 const std::vector<double>& splitPoints)
1965 {
1966 const double xyPosLimit = 70; //um
1967 const double xySize2Limit = pow(40, 2); //um^2
1968 const double zPosLimit = 1200; //um
1969
1970
1971 std::vector<double> indX = splitPoints;
1972 std::vector<double> indY = splitPoints;
1973 std::vector<double> indZ = splitPoints;
1974
1975 //no time dependence, as for beam size
1976 std::vector<double> indKX = {};
1977 std::vector<double> indKY = {};
1978
1979 UnknownPars allPars;
1980 for (int k = 0; k < 1; ++k) { //loop over BootStrap replicas
1981 for (auto& e : evts) e.isSig = true; //reset cuts
1982 if (k != 0) bootStrap(evts);
1983
1984
1985 //simple XY pos fit
1986 auto resTemp = fitSpotPositionSplines(evts, indX, indY);
1987
1988 const int kPlot = -1;
1989 // cppcheck-suppress knownConditionTrueFalse
1990 if (k == kPlot) {
1991 plotSpotPositionFit(evts, resTemp, "positionFitSimpe");
1992 plotSpotPositionPull(evts, resTemp, "pullsPositionSimple", xyPosLimit);
1993 }
1994 removeSpotPositionOutliers(evts, resTemp, xyPosLimit);
1995
1996 //simple XY pos fit (with outliers removed)
1997 auto resFin = fitSpotPositionSplines(evts, indX, indY);
1998 if (k == kPlot) {
1999 plotSpotPositionFit(evts, resFin, "positionFitSimpleC");
2000 plotSpotPositionPull(evts, resFin, "pullsPositionSimpleC", xyPosLimit);
2001 plotXYtimeDep(evts, resFin, "simplePosTimeDep");
2002 }
2003
2004 //Z position fit
2005 auto resZmy = fitZpositionSplinesSimple(evts, indZ, resFin);
2006 if (k == kPlot) {
2007 plotSpotZPositionFit(evts, resZmy, "positionFitSimpleZ");
2008 plotSpotZpositionPull(evts, resZmy, "zPositionPull", zPosLimit);
2009 }
2010
2011 removeSpotZpositionOutliers(evts, resZmy, zPosLimit);
2012
2013 //Z position fit (with outliers removed)
2014 resZmy = fitZpositionSplinesSimple(evts, indZ, resZmy);
2015
2016
2017 //complete XY pos fit
2018 auto resNew = fitSpotPositionSplines(evts, indX, indY, indKX, indKY, resZmy);
2019 if (k == kPlot) {
2020 plotSpotPositionFit(evts, resNew, "positionFitFull");
2021 plotKxKyFit(evts, resNew, "slopes");
2022 }
2023
2024 //Z position fit (iteration 2)
2025 resZmy = fitZpositionSplinesSimple(evts, indZ, resNew);
2026 if (k == kPlot) plotSpotZPositionFit(evts, resZmy, "positionFitSimpleZLast");
2027
2028
2029 //complete XY pos fit (iteration 2)
2030 resNew = fitSpotPositionSplines(evts, indX, indY, indKX, indKY, resZmy);
2031
2032 //XYZ pos fits (iteration 3)
2033 resZmy = fitZpositionSplinesSimple(evts, indZ, resNew);
2034 resNew = fitSpotPositionSplines(evts, indX, indY, indKX, indKY, resZmy);
2035
2036
2037 // fit of XY sizes (original + with outliers removed)
2038 auto vecXY = fitSpotWidthCMS(evts, resNew);
2039 if (k == kPlot) plotSpotSizePull(evts, resNew, vecXY);
2040 removeSpotSizeOutliers(evts, resNew, vecXY, xySize2Limit);
2041 vecXY = fitSpotWidthCMS(evts, resNew);
2042 if (k == kPlot) plotSpotSizeFit(evts, resNew, vecXY);
2043
2044
2045 // fit of Z size
2046 double sizeZZ = fitSpotZwidth(evts, resNew, vecXY);
2047
2048 if (k == kPlot) {
2049 plotSpotZSizeFit(evts, resNew, vecXY, sizeZZ);
2050 plotSpotSizeZPull(evts, resNew, vecXY, sizeZZ);
2051 }
2052
2053 //removeSpotSizeZOutliers(evts, resNew, vecXY, sizeZZ, 150000);
2054 //sizeZZ = fitSpotZwidth(evts, resNew, vecXY);
2055
2056 allPars.add(resNew, sqrtS(vecXY[0]), sqrtS(vecXY[1]), sqrtS(vecXY[2]), sqrtS(sizeZZ));
2057
2058 }
2059
2060 //allPars.printStat();
2061
2062 std::vector<VectorXd> vtxPos;
2063 std::vector<MatrixXd> vtxErr;
2064 MatrixXd sizeMat;
2065
2066 allPars.getOutput(vtxPos, vtxErr, sizeMat);
2067
2068 return std::make_tuple(vtxPos, vtxErr, sizeMat);
2069 }
2070
2071}
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::vec2vec
Eigen::VectorXd vec2vec(std::vector< double > vec)
std vector -> ROOT vector
Definition: tools.h:51
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(std::vector< std::vector< std::vector< double > > > toMerge)
merge { vector<double> a, vector<double> b} into {a, b}
Definition: tools.h:41
Belle2::vecs2mat
Eigen::MatrixXd vecs2mat(std::vector< std::vector< double > > vecs)
merge columns (from std::vectors) into ROOT matrix
Definition: tools.h:72
Belle2::rn
TString rn()
Get random string.
Definition: tools.h:38
VXDTFFilterTest::v2
const std::vector< double > v2
MATLAB generated random vector.
Definition: decorrelationMatrix.cc:42
VXDTFFilterTest::v1
const std::vector< double > v1
MATLAB generated random vector.
Definition: decorrelationMatrix.cc:28
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::UnknowSpline
Spline with uncertainty obtained from the boot-strap replicas.
Definition: BeamSpotStandAlone.cc:151
Belle2::BeamSpotCalib::UnknowSpline::add
void add(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::add
void add(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::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::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
void setBranchSpline(TTree *T, Spline *spl, TString n)
save bootstrap spline to TTree
Definition: BeamSpotStandAlone.cc:486
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::setBranchVal
void setBranchVal(TTree *T, std::vector< double > *vec, TString n)
save bootstrap variable to TTree
Definition: BeamSpotStandAlone.cc:478
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
Spline structure for zero-order & linear splines.
Definition: tools.h:165
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::print
void print(TString tag="")
print the spline
Definition: tools.h:188
Belle2::Spline::center
double center() const
Get center of the spline domain.
Definition: tools.h:177