Belle II Software  release-08-01-10
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 <tracking/calibration/BeamSpotStandAlone.h>
46 #include <tracking/calibration/Splitter.h>
47 #include <tracking/calibration/tools.h>
48 #else
49 #include <BeamSpotStandAlone.h>
50 #include <Splitter.h>
51 #include <tools.h>
52 #endif
53 
54 using Eigen::VectorXd;
55 using Eigen::Vector3d;
56 using Eigen::MatrixXd;
57 using Eigen::Matrix3d;
58 
59 
60 namespace 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::rn
TString rn()
Get random string.
Definition: tools.h:38
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::vec2vec
Eigen::VectorXd vec2vec(std::vector< double > vec)
std vector -> ROOT vector
Definition: tools.h:51
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::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::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::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::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::getStats
std::vector< double > getStats()
Get basic stats.
Definition: BeamSpotStandAlone.cc:258
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