Belle II Software  release-06-02-00
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 
55 using namespace std;
56 
57 using Eigen::VectorXd;
58 using Eigen::Vector3d;
59 using Eigen::MatrixXd;
60 using Eigen::Matrix3d;
61 
62 
63 namespace Belle2::BeamSpotCalib {
64 
65 
66  inline double sqrS(double x) {return x >= 0 ? x * x : -x * x; }
67  inline double sqrtS(double x) {return x >= 0 ? sqrt(x) : -sqrt(-x); }
68 
69  MatrixXd getRotatedSizeMatrix(vector<double> xySize, double zzSize, double kX, double kY);
70 
71 
73  struct SpotParam {
74  Spline x;
75  Spline y;
76  Spline z;
77  Spline kX;
78  Spline kY;
79 
80  SpotParam() {}
81 
83  void print()
84  {
85  B2INFO("x");
86  x.print();
87  B2INFO("y");
88  y.print();
89  B2INFO("z");
90  z.print();
91  B2INFO("kX");
92  kX.print();
93  B2INFO("kY");
94  kY.print();
95  }
96 
101  SpotParam(const vector<double>& vals, const vector<double>& errs, const vector<vector<double>>& spls, int order = 0)
102  {
103  auto getSize = [order](const vector<double>& sp) {
104  if (sp.size() == 0)
105  return 1;
106  else {
107  if (order == 0) {
108  B2ASSERT("There must be least one node at this place", sp.size() >= 1);
109  return int(sp.size() + 1);
110  } else if (order == 1) {
111  B2ASSERT("Must be at least two nodes in lin. spline", sp.size() >= 2);
112  return int(sp.size());
113  } else {
114  B2FATAL("Unknown order");
115  }
116  }
117  };
118 
119  int nx = getSize(spls[0]);
120  int ny = getSize(spls[1]);
121  x.nodes = spls[0];
122  y.nodes = spls[1];
123  x.vals = slice(vals, 0, nx);
124  y.vals = slice(vals, nx, ny);
125  x.errs = slice(errs, 0, nx);
126  y.errs = slice(errs, nx, ny);
127 
128  if (spls.size() >= 4) {
129  int nkx = getSize(spls[2]);
130  int nky = getSize(spls[3]);
131  kX.nodes = spls[2];
132  kY.nodes = spls[3];
133  kX.vals = slice(vals, nx + ny, nkx);
134  kY.vals = slice(vals, nx + ny + nkx, nky);
135 
136  kX.errs = slice(errs, nx + ny, nkx);
137  kY.errs = slice(errs, nx + ny + nkx, nky);
138 
139  if (spls.size() >= 5) {
140  int nz = getSize(spls[4]);
141  z.nodes = spls[4];
142  z.vals = slice(vals, nx + ny + nkx + nky, nz);
143  z.errs = slice(errs, nx + ny + nkx + nky, nz);
144  }
145  }
146  }
147 
148 
149  };
150 
151 
153  struct UnknowSpline {
154  vector<Spline> spls;
155  void add(Spline spl) { spls.push_back(spl); }
156 
158  Spline getMeanSigma()
159  {
160  Spline sAvg;
161  sAvg = spls[0];
162  int nNd = spls[0].vals.size();
163  for (int k = 0; k < nNd; ++k) {
164  double s = 0, ss = 0;
165  for (unsigned i = 0; i < spls.size(); ++i) {
166  s += spls[i].vals[k];
167  }
168  s /= spls.size();
169 
170  for (unsigned i = 0; i < spls.size(); ++i) {
171  ss += pow(spls[i].vals[k] - s, 2);
172  }
173 
174  if (spls.size() > 1)
175  ss = sqrt(ss / (spls.size() - 1));
176  else
177  ss = 0;
178 
179  sAvg.vals[k] = s;
180  sAvg.errs[k] = ss;
181 
182  }
183  return sAvg;
184  }
185 
188  Spline getLimit(double v)
189  {
190  Spline sLim = spls[0];
191 
192  double indx = (spls.size() - 1) * v;
193  int nNd = spls[0].vals.size();
194  for (int k = 0; k < nNd; ++k) {
195  vector<double> vals;
196  for (unsigned i = 0; i < spls.size(); ++i) {
197  vals.push_back(spls[i].vals[k]);
198  }
199  sort(vals.begin(), vals.end());
200 
201  int I = indx;
202  double r = indx - I;
203  sLim.vals[k] = vals[I] * (1 - r) + vals[I + 1] * r;
204  sLim.errs[k] = 0;
205  }
206  return sLim;
207  }
208 
209  };
210 
211 
212 
213 
215  struct UnknowVar {
216  vector<double> vars;
217  void add(double x) { vars.push_back(x); }
218 
220  double getMean()
221  {
222  B2ASSERT("Must be at least one replica", vars.size() >= 1);
223  return accumulate(vars.begin(), vars.end(), 0.) / vars.size();
224  }
225 
227  double getSigma()
228  {
229  B2ASSERT("Must be at least one replica", vars.size() >= 1);
230  double m = getMean();
231  double s = 0;
232  for (auto x : vars)
233  s += pow(x - m, 2);
234  if (vars.size() > 1)
235  return sqrt(s / (vars.size() - 1));
236  else
237  return 0; //dummy unc. for single replica
238  }
239 
241  double getLimit(double v)
242  {
243  B2ASSERT("Must be at least one replica", vars.size() >= 1);
244  double indx = (vars.size() - 1) * v;
245  sort(vars.begin(), vars.end());
246  int I = indx;
247  double r = indx - I;
248  return vars[I] * (1 - r) + vars[I + 1] * r;
249  }
250 
252  void printStat(TString n)
253  {
254  B2ASSERT("Must be at least one replica", vars.size() >= 1);
255  B2INFO(n << " : " << getMean() << "+-" << getSigma() << " : " << getLimit(0.50) << " (" << getLimit(0.16) << " , " << getLimit(
256  1 - 0.16) << " )");
257  }
258 
260  vector<double> getStats()
261  {
262  return {getLimit(0.50), getLimit(0.16), getLimit(1 - 0.16)};
263  }
264  };
265 
266 
273  double getAngle(double SizeX, double SizeY, double SizeXY)
274  {
275  double C = sqrS(SizeXY);
276  // By convention the range of angle is [-pi/2, pi/2],
277  // the reason is the same as with straight line, e.g. it can point to pi/4 but also to -3/4*pi
278  double angle = 1. / 2 * atan2(2 * C, (pow(SizeX, 2) - pow(SizeY, 2)));
279  return angle;
280  }
281 
282 
289  pair<double, double> getSizeMinMax(double SizeX, double SizeY, double SizeXY)
290  {
291  double A = pow(SizeX, 2) + pow(SizeY, 2);
292  double B = pow(SizeX, 2) * pow(SizeY, 2) - pow(SizeXY, 4);
293  double D = pow(A, 2) - 4 * B;
294 
295  if (D < 0) {
296  B2FATAL("Problem with D value : " << D);
297  }
298 
299  double Size2Min = 2 * B / (A + sqrt(D));
300  if (abs(Size2Min) < 1e-7) Size2Min = 0;
301  if (Size2Min < 0) {
302  B2FATAL("Negative BS eigen size : " << Size2Min);
303  }
304  double Size2Max = (A + sqrt(D)) / 2;
305  return {sqrtS(Size2Min), sqrtS(Size2Max)};
306  }
307 
308 
315  double getSize2MinMat(double SizeXX, double SizeYY, double SizeXY)
316  {
317  double A = SizeXX + SizeYY;
318  double B = SizeXX * SizeYY - pow(SizeXY, 2);
319  double D = pow(A, 2) - 4 * B;
320 
321  if (D < 0) {
322  B2FATAL("Problem with D value : " << D);
323  }
324 
325  double Size2Min = 2 * B / (A + sqrt(D));
326 
327  return Size2Min;
328  }
329 
330 
331 
332 
333 
335  struct UnknownPars {
336  UnknowSpline x;
337  UnknowSpline y;
338  UnknowSpline kX;
339  UnknowSpline kY;
340  UnknowSpline z;
341 
342  UnknowVar sizeX;
343  UnknowVar sizeY;
344  UnknowVar sizeXY;
345  UnknowVar sizeMin;
346  UnknowVar sizeMax;
347  UnknowVar xyAngle;
348 
349 
350  UnknowVar sizeZ;
351  UnknowVar crossAngle;
352 
353  UnknowVar matXX;
354  UnknowVar matYY;
355  UnknowVar matZZ;
356  UnknowVar matXY;
357  UnknowVar matXZ;
358  UnknowVar matYZ;
359 
361  void add(SpotParam sPar, double SizeX, double SizeY, double SizeXY, double SizeZ)
362  {
363  x.add(sPar.x);
364  y.add(sPar.y);
365  kX.add(sPar.kX);
366  kY.add(sPar.kY);
367  z.add(sPar.z);
368 
369 
370  sizeX.add(SizeX);
371  sizeY.add(SizeY);
372  sizeXY.add(SizeXY);
373  sizeZ.add(SizeZ);
374 
375 
376  // Calculate the eigen-values
377  double SizeMin, SizeMax;
378  tie(SizeMin, SizeMax) = getSizeMinMax(SizeX, SizeY, SizeXY);
379 
380  sizeMin.add(SizeMin);
381  sizeMax.add(SizeMax);
382 
383  // and angle in mrad
384  double angle = 1e3 * getAngle(SizeX, SizeY, SizeXY);
385 
386 
387  xyAngle.add(angle);
388 
389  //Get whole cov matrix
390  MatrixXd matSize = getRotatedSizeMatrix({sqrS(SizeX), sqrS(SizeY), sqrS(SizeXY)}, sqrS(SizeZ), sPar.kX.val(sPar.kX.center()),
391  sPar.kY.val(sPar.kY.center()));
392 
393  // Store elements in [um]
394  matXX.add(sqrtS(matSize(0, 0)));
395  matYY.add(sqrtS(matSize(1, 1)));
396  matZZ.add(sqrtS(matSize(2, 2)));
397  matXY.add(sqrtS(matSize(0, 1)));
398 
399  matXZ.add(sqrtS(matSize(0, 2)));
400  matYZ.add(sqrtS(matSize(1, 2)));
401 
402 
403  // crossing-angle in mrad
404  double crossAngleVal = 1e3 * 2 * sqrtS(matSize(0, 0)) / sqrtS(matSize(2, 2));
405  crossAngle.add(crossAngleVal);
406  }
407 
409  void printStat()
410  {
411  x.getMeanSigma().print("x");
412  y.getMeanSigma().print("y");
413  kX.getMeanSigma().print("kX");
414  kY.getMeanSigma().print("kY");
415  z.getMeanSigma().print("z");
416 
417  sizeX.printStat("sizeX");
418  sizeY.printStat("sizeY");
419  sizeXY.printStat("sizeXY");
420 
421 
422  sizeMin.printStat("sizeMin");
423  sizeMax.printStat("sizeMax");
424  xyAngle.printStat("xyAngle");
425  sizeZ.printStat("sizeZ");
426  crossAngle.printStat("crossAngle");
427 
428 
429  matXX.printStat("matXX");
430  matYY.printStat("matYY");
431  matZZ.printStat("matZZ");
432  matXY.printStat("matXY");
433  matXZ.printStat("matXZ");
434  matYZ.printStat("matYZ");
435  }
436 
438  void getOutput(vector<VectorXd>& vtxPos, vector<MatrixXd>& vtxErr, MatrixXd& sizeMat)
439  {
440  //Store the vertex position
441  int nVals = x.spls[0].vals.size();
442 
443  vtxPos.clear();
444  vtxErr.clear();
445 
446  const double toCm = 1e-4;
447 
448  for (int i = 0; i < nVals; ++i) {
449  //vertex position
450  Vector3d vtx(x.spls[0].vals[i]*toCm, y.spls[0].vals[i]*toCm, z.spls[0].vals[i]*toCm);
451 
452  //vertex error matrix (symmetric)
453  Matrix3d mS = Matrix3d::Zero();
454  mS(0, 0) = sqrS(x.spls[0].errs[i] * toCm);
455  mS(1, 1) = sqrS(y.spls[0].errs[i] * toCm);
456  mS(2, 2) = sqrS(z.spls[0].errs[i] * toCm);
457 
458  vtxPos.push_back(vtx);
459  vtxErr.push_back(mS);
460  }
461 
462  //BeamSpot size matrix (from boot-strap iteration 0)
463 
464  sizeMat.resize(3, 3);
465  sizeMat(0, 0) = sqrS(matXX.vars[0] * toCm);
466  sizeMat(1, 1) = sqrS(matYY.vars[0] * toCm);
467  sizeMat(2, 2) = sqrS(matZZ.vars[0] * toCm);
468 
469  sizeMat(0, 1) = sqrS(matXY.vars[0] * toCm);
470  sizeMat(0, 2) = sqrS(matXZ.vars[0] * toCm);
471  sizeMat(1, 2) = sqrS(matYZ.vars[0] * toCm);
472 
473  sizeMat(1, 0) = sizeMat(0, 1);
474  sizeMat(2, 0) = sizeMat(0, 2);
475  sizeMat(2, 1) = sizeMat(1, 2);
476  }
477 
478 
480  void setBranchVal(TTree* T, vector<double>* vec, TString n)
481  {
482  T->Branch(n, &vec->at(0), n + "/D");
483  T->Branch(n + "_low", &vec->at(1), n + "_low/D");
484  T->Branch(n + "_high", &vec->at(2), n + "_high/D");
485  }
486 
488  void setBranchSpline(TTree* T, Spline* spl, TString n)
489  {
490  T->Branch(n + "_nodes", &spl->nodes);
491  T->Branch(n + "_vals", &spl->vals);
492  T->Branch(n + "_errs", &spl->errs);
493  }
494 
495 
496 
498  void save2tree(TString fName)
499  {
500 
501  TTree* T = new TTree("runs", "beam conditions of runs");
502 
503  int run = -99; //currently dummy
504  T->Branch("run", &run, "run/I");
505 
506  Spline xAvg = x.getMeanSigma();
507  setBranchSpline(T, &xAvg, "x");
508  Spline yAvg = y.getMeanSigma();
509  setBranchSpline(T, &yAvg, "y");
510  Spline zAvg = z.getMeanSigma();
511  setBranchSpline(T, &zAvg, "z");
512 
513  Spline kxAvg = kX.getMeanSigma();
514  setBranchSpline(T, &kxAvg, "kX");
515  Spline kyAvg = kY.getMeanSigma();
516  setBranchSpline(T, &kyAvg, "kY");
517 
518  vector<double> sizeXVar = sizeX.getStats();
519  setBranchVal(T, &sizeXVar, "sizeX");
520  vector<double> sizeYVar = sizeY.getStats();
521  setBranchVal(T, &sizeYVar, "sizeY");
522  vector<double> sizeXYVar = sizeXY.getStats();
523  setBranchVal(T, &sizeXYVar, "sizeXY");
524  vector<double> sizeZVar = sizeZ.getStats();
525  setBranchVal(T, &sizeZVar, "sizeZ");
526 
527  vector<double> sizeMinVar = sizeMin.getStats();
528  setBranchVal(T, &sizeMinVar, "sizeMin");
529  vector<double> sizeMaxVar = sizeMax.getStats();
530  setBranchVal(T, &sizeMaxVar, "sizeMax");
531  vector<double> xyAngleVar = xyAngle.getStats();
532  setBranchVal(T, &xyAngleVar, "xyAngle");
533 
534  vector<double> crossAngleVar = crossAngle.getStats();
535  setBranchVal(T, &crossAngleVar, "crossAngle");
536 
537 
538  vector<double> matXXVar = matXX.getStats();
539  vector<double> matYYVar = matYY.getStats();
540  vector<double> matZZVar = matZZ.getStats();
541  vector<double> matXYVar = matXY.getStats();
542  vector<double> matXZVar = matXZ.getStats();
543  vector<double> matYZVar = matYZ.getStats();
544 
545  setBranchVal(T, &matXXVar, "matXX");
546  setBranchVal(T, &matYYVar, "matYY");
547  setBranchVal(T, &matZZVar, "matZZ");
548 
549  setBranchVal(T, &matXYVar, "matXY");
550  setBranchVal(T, &matXZVar, "matXZ");
551  setBranchVal(T, &matYZVar, "matYZ");
552 
553 
554  T->Fill();
555  T->SaveAs(fName);
556  }
557  };
558 
559 
560 
569  double getZIPest(const Track& tr, double t, const SpotParam& spotPar, int nestMax = 5, int nest = 0)
570  {
571  double x0, y0;
572  if (nest < nestMax) {
573  double zIP = getZIPest(tr, t, spotPar, nestMax, nest + 1);
574 
575  x0 = spotPar.x.val(t) + spotPar.kX.val(t) * (zIP - spotPar.z.val(t));
576  y0 = spotPar.y.val(t) + spotPar.kY.val(t) * (zIP - spotPar.z.val(t));
577  } else {
578  x0 = spotPar.x.val(t);
579  y0 = spotPar.y.val(t);
580  }
581 
582  double f0 = tr.tanlambda * (x0 * cos(tr.phi0) + y0 * sin(tr.phi0));
583 
584  return (tr.z0 + f0);
585  }
586 
587 
594  double getCorrD(const Track& tr, double t, const SpotParam& spotPar)
595  {
596  double zIP = getZIPest(tr, t, spotPar);
597 
598  double x0 = spotPar.x.val(t) + spotPar.kX.val(t) * (zIP - spotPar.z.val(t));
599  double y0 = spotPar.y.val(t) + spotPar.kY.val(t) * (zIP - spotPar.z.val(t));
600 
601  double f0 = x0 * sin(tr.phi0) - y0 * cos(tr.phi0);
602 
603  return (tr.d0 - f0);
604  }
605 
607  double getDtimeConst(const Track& tr, double t, const SpotParam& spotPar)
608  {
609  double zIP = getZIPest(tr, t, spotPar);
610  double zIPM = getZIPest(tr, spotPar.z.center(), spotPar);
611 
612  double x0 = spotPar.x.val(t) + spotPar.kX.val(t) * (zIP - spotPar.z.val(t));
613  double y0 = spotPar.y.val(t) + spotPar.kY.val(t) * (zIP - spotPar.z.val(t));
614 
615  double xM = spotPar.x.val(spotPar.x.center()) + spotPar.kX.val(spotPar.kX.center()) * (zIPM - spotPar.z.val(spotPar.z.center()));
616  double yM = spotPar.y.val(spotPar.y.center()) + spotPar.kY.val(spotPar.kY.center()) * (zIPM - spotPar.z.val(spotPar.z.center()));
617 
618 
619  double f0 = (x0 - xM) * sin(tr.phi0) - (y0 - yM) * cos(tr.phi0);
620 
621  return (tr.d0 - f0);
622  }
623 
624 
631  double getCorrZ(const Track& tr, double t, const SpotParam& spotPar)
632  {
633  double zIP = getZIPest(tr, t, spotPar);
634 
635  double x0 = spotPar.x.val(t) + spotPar.kX.val(t) * (zIP - spotPar.z.val(t));
636  double y0 = spotPar.y.val(t) + spotPar.kY.val(t) * (zIP - spotPar.z.val(t));
637  double z0 = spotPar.z.val(t);
638 
639  double f0 = z0 - tr.tanlambda * (x0 * cos(tr.phi0) + y0 * sin(tr.phi0));
640 
641  return (tr.z0 - f0);
642  }
643 
644 
645 
646 
647 
652  pair<double, double> getStartStop(const vector<Event>& evts)
653  {
654  double minT = 1e20, maxT = -1e20;
655  for (auto ev : evts) {
656  minT = min(minT, ev.t);
657  maxT = max(maxT, ev.t);
658  }
659  return {minT, maxT};
660  }
661 
663  vector<TString> extractFileNames(TString str)
664  {
665  vector<TString> files;
666  auto tempVec = str.Tokenize(",");
667  for (int i = 0; i < tempVec->GetEntries(); ++i) {
668  TString s = ((TObjString*)tempVec->At(i))->GetString();
669  files.push_back(s.Strip());
670  }
671  return files;
672  }
673 
675  vector<Event> getEvents(TTree* tr)
676  {
677 
678  vector<Event> events;
679  events.reserve(tr->GetEntries());
680 
681  Event evt;
682 
683  tr->SetBranchAddress("run", &evt.run);
684  tr->SetBranchAddress("exp", &evt.exp);
685  tr->SetBranchAddress("event", &evt.evtNo);
686  tr->SetBranchAddress("mu0_d0", &evt.mu0.d0);
687  tr->SetBranchAddress("mu1_d0", &evt.mu1.d0);
688  tr->SetBranchAddress("mu0_z0", &evt.mu0.z0);
689  tr->SetBranchAddress("mu1_z0", &evt.mu1.z0);
690 
691  tr->SetBranchAddress("mu0_tanlambda", &evt.mu0.tanlambda);
692  tr->SetBranchAddress("mu1_tanlambda", &evt.mu1.tanlambda);
693 
694 
695  tr->SetBranchAddress("mu0_phi0", &evt.mu0.phi0);
696  tr->SetBranchAddress("mu1_phi0", &evt.mu1.phi0);
697 
698  tr->SetBranchAddress("time", &evt.t); //time in hours
699 
700 
701  for (int i = 0; i < tr->GetEntries(); ++i) {
702  tr->GetEntry(i);
703  evt.toMicroM();
704 
705  evt.nBootStrap = 1;
706  evt.isSig = true;
707  events.push_back(evt);
708  }
709 
710  //sort by time
711  sort(events.begin(), events.end(), [](Event e1, Event e2) {return e1.t < e2.t;});
712 
713 
714  return events;
715  }
716 
718  void bootStrap(vector<Event>& evts)
719  {
720  for (auto& e : evts)
721  e.nBootStrap = gRandom->Poisson(1);
722  }
723 
724 
725 
726 
732  VectorXd linearFit(MatrixXd mat, VectorXd r)
733  {
734  MatrixXd matT = mat.transpose();
735  MatrixXd A = matT * mat;
736 
737  VectorXd v = matT * r;
738  MatrixXd Ainv = A.inverse();
739 
740  return (Ainv * v);
741  }
742 
751  pair<vector<double>, vector<double>> linearFitErr(MatrixXd m, VectorXd r, double& err2Mean, double& err2press, double& err2pressErr)
752  {
753  MatrixXd mT = m.transpose();
754  MatrixXd mat = mT * m;
755 
756  mat = mat.inverse();
757  MatrixXd A = mat * mT;
758  VectorXd res = A * r;
759  VectorXd dataH = m * res;
760 
761 
762  //errs
763  double err2 = (dataH - r).squaredNorm() / (r.rows() - res.rows());
764 
765 
766  // Get PRESS statistics of the linear fit
767  {
768  //Ahat = m*A;
769  double press = 0;
770  double press2 = 0;
771  for (int i = 0; i < r.rows(); ++i) {
772  double Ahat = 0;
773  for (int k = 0; k < m.cols(); ++k)
774  Ahat += m(i, k) * A(k, i);
775 
776  double v = pow((r(i) - dataH(i)) / (1 - Ahat), 2);
777  press += v;
778  press2 += v * v;
779  }
780  press /= r.rows();
781  press2 /= (r.rows() - 1);
782 
783  err2press = press;
784  err2pressErr = sqrt((press2 - press * press) / r.rows()) / sqrt(r.rows());
785  }
786 
787 
788 
789  MatrixXd AT = A.transpose(); // A.T();
790  MatrixXd errMat = err2 * A * AT;
791  VectorXd errs2(errMat.rows());
792  for (int i = 0; i < errs2.rows(); ++i)
793  errs2(i) = errMat(i, i);
794 
795 
796  err2Mean = err2;
797 
798  return {vec2vec(res), vec2vec(errs2)};
799  }
800 
801 
802 
803 
809  VectorXd linearFitPos(MatrixXd mat, VectorXd r)
810  {
811  const double s2MinLimit = pow(0.05, 2); //Minimal value of the BeamSpot eigenSize
812  B2ASSERT("Matrix size for size fit should be 3", mat.cols() == 3);
813  MatrixXd matT = mat.transpose();
814  MatrixXd A = matT * mat;
815  VectorXd v = matT * r;
816  MatrixXd Ainv = A.inverse();
817 
818  VectorXd pars = Ainv * v;
819 
820  //If everything is OK
821  double s2Min = getSize2MinMat(pars[0], pars[1], pars[2]);
822  if (pars[0] >= 0 && pars[1] >= 0 && s2Min >= s2MinLimit)
823  return pars;
824 
826  //Get the error matrix
828 
829  int nDf = r.rows() - 3;
830  //Calculate average error
831  double err2 = (mat * pars - r).squaredNorm() / nDf;
832 
833  MatrixXd wMat = Ainv * matT;
834  MatrixXd wMatT = wMat.transpose();
835 
836  MatrixXd covMat = err2 * wMat * wMatT;
837  MatrixXd covMatI = covMat.inverse();
838 
840  //Get maximum likelihood
842 
843  //Maximum likelihood over eigenvector and angle
844  TF2 fEig(rn(), [covMatI, pars, s2MinLimit](double * x, double*) {
845  double eig1 = x[0];
846  double eig2 = s2MinLimit;
847  double phi = x[1];
848  double c = cos(phi);
849  double s = sin(phi);
850 
851  VectorXd xVec(3);
852  xVec[0] = eig1 * c * c + eig2 * s * s;
853  xVec[1] = eig1 * s * s + eig2 * c * c;
854  xVec[2] = s * c * (eig1 - eig2);
855 
856  //double res = covMatI.Similarity(xVec - pars);
857  double res = (xVec - pars).transpose() * covMatI * (xVec - pars);
858  return res;
859  }, s2MinLimit, 400, 0, 2 * M_PI, 0);
860 
861  double eigHigh, phi;
862  fEig.GetMinimumXY(eigHigh, phi);
863 
864  pars[0] = eigHigh * pow(cos(phi), 2) + s2MinLimit * pow(sin(phi), 2);
865  pars[1] = eigHigh * pow(sin(phi), 2) + s2MinLimit * pow(cos(phi), 2);
866  pars[2] = sin(phi) * cos(phi) * (eigHigh - s2MinLimit);
867 
868 
869  return pars;
870 
871  }
872 
873 
874 
875 
877  TH1D* getResolution(TH2D* hRes)
878  {
879  TH1D* hSigmaAll = new TH1D(rn(), "", 50, -M_PI, M_PI);
880  for (int i = 1; i <= hRes->GetNbinsX(); ++i) {
881  TH1D* hProj = hRes->ProjectionY(rn(), i, i);
882  double rms = hProj->GetRMS();
883  double rmsErr = hProj->GetRMSError();
884  hSigmaAll->SetBinContent(i, rms * rms); //from cm2 to um2
885  hSigmaAll->SetBinError(i, 2 * abs(rms)*rmsErr);
886  }
887  return hSigmaAll;
888  }
889 
891  TH1D* getMean(const TH2D* hRes)
892  {
893  TH1D* hMean = new TH1D(rn(), "", 50, -M_PI, M_PI);
894  for (int i = 1; i <= hRes->GetNbinsX(); ++i) {
895  TH1D* hProj = hRes->ProjectionY(rn(), i, i);
896  double mean = hProj->GetMean();
897  double meanErr = hProj->GetMeanError();
898  hMean->SetBinContent(i, mean);
899  hMean->SetBinError(i, meanErr);
900  }
901  return hMean;
902  }
903 
905  double getD12th(Event e, vector<double> sizesXY)
906  {
907  double sxx = sizesXY[0];
908  double syy = sizesXY[1];
909  double sxy = sizesXY[2];
910 
911  double cc = cos(e.mu0.phi0) * cos(e.mu1.phi0);
912  double ss = sin(e.mu0.phi0) * sin(e.mu1.phi0);
913  double sc = -(sin(e.mu0.phi0) * cos(e.mu1.phi0) + sin(e.mu1.phi0) * cos(e.mu0.phi0));
914 
915  return ss * sxx + cc * syy + sc * sxy;
916  }
917 
919  double getZ12th(Event e, vector<double> sizesXY)
920  {
921  double sxx = sizesXY[0];
922  double syy = sizesXY[1];
923 
924  double corr = e.mu0.tanlambda * e.mu1.tanlambda * (sxx * cos(e.mu0.phi0) * cos(e.mu1.phi0) + syy * sin(e.mu0.phi0) * sin(
925  e.mu1.phi0) +
926  + (sin(e.mu0.phi0) * cos(e.mu1.phi0) + cos(e.mu0.phi0) * sin(e.mu1.phi0)));
927  return corr;
928  }
929 
930 
931 
932 
933 
934 
936  void plotSpotPositionFit(const vector<Event>& evts, SpotParam par, TString fName)
937  {
938  TGraph* gr = new TGraph();
939  TProfile* dProf = new TProfile(rn(), "dProf", 100, -M_PI, M_PI, "S");
940  TProfile* dProfRes = new TProfile(rn(), "dProfRes", 100, -M_PI, M_PI, "S");
941 
942  for (auto e : evts) {
943  if (!e.isSig) continue;
944 
945  double d1 = getDtimeConst(e.mu0, e.t, par);
946  double d2 = getDtimeConst(e.mu1, e.t, par);
947 
948  gr->SetPoint(gr->GetN(), e.mu0.phi0, d1);
949  gr->SetPoint(gr->GetN(), e.mu1.phi0, d2);
950 
951  dProf->Fill(e.mu0.phi0, d1);
952  dProf->Fill(e.mu1.phi0, d2);
953 
954 
955  double d1c = getCorrD(e.mu0, e.t, par);
956  double d2c = getCorrD(e.mu1, e.t, par);
957 
958  dProfRes->Fill(e.mu0.phi0, d1c);
959  dProfRes->Fill(e.mu1.phi0, d2c);
960  }
961  TF1* f = new TF1(rn(), "[0]*sin(x) - [1]*cos(x)", -M_PI, M_PI);
962  f->SetParameters(par.x.val(par.x.center()), par.y.val(par.y.center()));
963 
964  TCanvas* c = new TCanvas(rn(), "");
965  gr->Draw("a p");
966  gr->GetXaxis()->SetRangeUser(-M_PI, M_PI);
967  gr->SetMaximum(+1.3 * f->GetMaximum());
968  gr->SetMinimum(-1.3 * f->GetMaximum());
969 
970  gr->GetXaxis()->SetTitle("#phi_{0} [rad]");
971  gr->GetYaxis()->SetTitle("d_{0} [#mum]");
972 
973  f->Draw("same");
974  c->SaveAs(fName + "_dots.pdf");
975 
976 
977  TCanvas* d = new TCanvas(rn(), "");
978  gStyle->SetOptStat(0);
979  dProf->Draw();
980  dProf->GetXaxis()->SetRangeUser(-M_PI, M_PI);
981  dProf->SetMaximum(+1.3 * f->GetMaximum());
982  dProf->SetMinimum(-1.3 * f->GetMaximum());
983 
984  dProf->GetXaxis()->SetTitle("#phi_{0} [rad]");
985  dProf->GetYaxis()->SetTitle("d_{0} [#mum]");
986 
987  f->Draw("same");
988 
989  B2INFO("Saving " << fName << " prof ");
990  d->SaveAs(fName + "_prof.pdf");
991 
992 
993  TCanvas* e = new TCanvas(rn(), "");
994  gStyle->SetOptStat(0);
995  dProfRes->Draw();
996  dProfRes->GetXaxis()->SetRangeUser(-M_PI, M_PI);
997 
998  dProfRes->GetXaxis()->SetTitle("#phi_{0} [rad]");
999  dProfRes->GetYaxis()->SetTitle("d_{0} res [#mum]");
1000 
1001  TH1D* errP = new TH1D(rn(), "dErrP", 100, -M_PI, M_PI);
1002  TH1D* errM = new TH1D(rn(), "dErrM", 100, -M_PI, M_PI);
1003  for (int i = 1; i <= errP->GetNbinsX(); ++i) {
1004  errP->SetBinContent(i, dProfRes->GetBinError(i));
1005  errM->SetBinContent(i, -dProfRes->GetBinError(i));
1006  }
1007 
1008  errP->Draw("hist same");
1009  errM->Draw("hist same");
1010 
1011  f->SetParameters(0, 0);
1012  f->Draw("same");
1013 
1014  B2INFO("Saving " << fName << " profRes ");
1015  e->SaveAs(fName + "_profRes.pdf");
1016 
1017  }
1018 
1019 
1021  void plotSpotZPositionFit(const vector<Event>& evts, SpotParam par, TString fName)
1022  {
1023  TProfile* zProf = new TProfile(rn(), "dProf", 100, -M_PI, M_PI, "S");
1024  TGraph* gr = new TGraph();
1025  for (auto e : evts) {
1026  if (!e.isSig) continue;
1027 
1028 
1029  double z1ip = getZIPest(e.mu0, e.t, par);
1030  double z2ip = getZIPest(e.mu1, e.t, par);
1031 
1032  double zipT = par.z.val(e.t);
1033  double zipM = par.z.val(par.z.center());
1034 
1035  double val1 = z1ip - (zipT - zipM);
1036  double val2 = z2ip - (zipT - zipM);
1037 
1038  gr->SetPoint(gr->GetN(), e.mu0.phi0, val1);
1039  gr->SetPoint(gr->GetN(), e.mu1.phi0, val2);
1040 
1041  zProf->Fill(e.mu0.phi0, val1);
1042  zProf->Fill(e.mu1.phi0, val2);
1043  }
1044  TF1* f = new TF1(rn(), "[0]", -M_PI, M_PI);
1045  f->SetParameter(0, par.z.val(par.z.center()));
1046 
1047  TCanvas* c = new TCanvas(rn(), "");
1048  c->SetLeftMargin(0.12);
1049  gr->Draw("a p");
1050  gr->GetXaxis()->SetRangeUser(-M_PI, M_PI);
1051  gr->SetMaximum(1000);
1052  gr->SetMinimum(-1000);
1053 
1054  gr->GetXaxis()->SetTitle("#phi_{0} [rad]");
1055  gr->GetYaxis()->SetTitle("z_{IP} estimate [#mum]");
1056 
1057  f->Draw("same");
1058  c->SaveAs(fName + "_dots.pdf");
1059 
1060 
1061  TCanvas* d = new TCanvas(rn(), "");
1062  gStyle->SetOptStat(0);
1063  d->SetLeftMargin(0.12);
1064  zProf->Draw();
1065  zProf->GetXaxis()->SetRangeUser(-M_PI, M_PI);
1066  zProf->SetMaximum(1000);
1067  zProf->SetMinimum(-1000);
1068 
1069  zProf->GetXaxis()->SetTitle("#phi_{0} [rad]");
1070  zProf->GetYaxis()->SetTitle("z_{IP} estimate [#mum]");
1071 
1072  f->Draw("same");
1073  d->SaveAs(fName + +"_prof.pdf");
1074 
1075  }
1076 
1077 
1078 
1079 
1080 
1082  void plotSpotPositionPull(const vector<Event>& evts, const SpotParam& par, TString fName, double cut = 70)
1083  {
1084  TH1D* hPull = new TH1D(rn(), "", 200, -200, 200);
1085 
1086  for (auto& e : evts) {
1087  if (!e.isSig) continue;
1088 
1089  double d0 = getCorrD(e.mu0, e.t, par);
1090  double d1 = getCorrD(e.mu1, e.t, par);
1091 
1092  hPull->Fill(d0);
1093  hPull->Fill(d1);
1094  }
1095 
1096  TCanvas* c = new TCanvas(rn(), "");
1097  gStyle->SetOptStat(2210);
1098  hPull->Draw();
1099 
1100  hPull->GetXaxis()->SetTitle("pull [#mum]");
1101  hPull->GetYaxis()->SetTitle("#tracks");
1102 
1103  TLine* l = new TLine;
1104  l->SetLineColor(kRed);
1105  l->DrawLine(-cut, 0, -cut, 500);
1106  l->DrawLine(+cut, 0, +cut, 500);
1107 
1108  c->SaveAs(fName + ".pdf");
1109  }
1110 
1111 
1113  void plotKxKyFit(const vector<Event>& evts, SpotParam par, TString fName)
1114  {
1115  TProfile* profRes = new TProfile(rn(), "dProf", 100, -800, 800, "S");
1116  TProfile* profResKx = new TProfile(rn(), "dProfKx", 100, -800, 800, "S");
1117  TProfile* profResKy = new TProfile(rn(), "dProfKy", 100, -800, 800, "S");
1118 
1119  SpotParam parNoKx = par;
1120  SpotParam parNoKy = par;
1121  parNoKx.kX.vals = {0};
1122  parNoKy.kY.vals = {0};
1123  parNoKx.kX.nodes = {};
1124  parNoKy.kY.nodes = {};
1125 
1126  for (auto& e : evts) {
1127  if (!e.isSig) continue;
1128 
1129  double zDiff1 = getZIPest(e.mu0, e.t, par) - par.z.val(e.t);
1130  double zDiff2 = getZIPest(e.mu1, e.t, par) - par.z.val(e.t);
1131 
1132 
1133  double d1 = getCorrD(e.mu0, e.t, par);
1134  double d2 = getCorrD(e.mu1, e.t, par);
1135 
1136  double d1KX = getCorrD(e.mu0, e.t, parNoKx);
1137  double d2KX = getCorrD(e.mu1, e.t, parNoKx);
1138 
1139 
1140  profResKx->Fill(zDiff1 * sin(e.mu0.phi0), d1KX);
1141  profResKx->Fill(zDiff2 * sin(e.mu1.phi0), d2KX);
1142 
1143  double d1KY = getCorrD(e.mu0, e.t, parNoKy);
1144  double d2KY = getCorrD(e.mu1, e.t, parNoKy);
1145  profResKy->Fill(-zDiff1 * cos(e.mu0.phi0), d1KY);
1146  profResKy->Fill(-zDiff2 * cos(e.mu1.phi0), d2KY);
1147 
1148  profRes->Fill(zDiff1, d1);
1149  profRes->Fill(zDiff2, d2);
1150 
1151 
1152  }
1153 
1154  TCanvas* cX = new TCanvas(rn(), "");
1155  gStyle->SetOptStat(0);
1156  profResKx->Draw();
1157 
1158  profResKx->GetXaxis()->SetTitle("(z_{IP} - z_{IP}^{0}) sin #phi_{0} [#mum]");
1159  profResKx->GetYaxis()->SetTitle("d_{0} res [#mum]");
1160 
1161  TF1* f = new TF1(rn(), "[0]*x", -800, 800);
1162  f->SetParameter(0, par.kX.val(par.kX.center()));
1163  f->Draw("same");
1164 
1165  cX->SaveAs(fName + "_kX.pdf");
1166 
1167  TCanvas* cY = new TCanvas(rn(), "");
1168  gStyle->SetOptStat(0);
1169  profResKy->Draw();
1170 
1171  profResKy->GetXaxis()->SetTitle("-(z_{IP} - z_{IP}^{0}) cos #phi_{0} [#mum]");
1172  profResKy->GetYaxis()->SetTitle("d_{0} res [#mum]");
1173 
1174  f->SetParameter(0, par.kY.val(par.kY.center()));
1175  f->Draw("same");
1176 
1177  cY->SaveAs(fName + "_kY.pdf");
1178 
1179 
1180  }
1181 
1183  void plotXYtimeDep(const vector<Event>& evts, SpotParam par, TString fName)
1184  {
1185  TProfile* profRes = new TProfile(rn(), "dProf", 50, -0.5, 0.5);
1186  TProfile* profResTx = new TProfile(rn(), "dProfTx", 50, -0.5, 0.5);
1187  TProfile* profResTy = new TProfile(rn(), "dProfTy", 50, -0.5, 0.5);
1188 
1189  SpotParam parNoTx = par;
1190  SpotParam parNoTy = par;
1191  parNoTx.x.nodes = {};
1192  parNoTx.x.vals = {par.x.val(par.x.center())};
1193  parNoTy.y.nodes = {};
1194  parNoTy.y.vals = {par.y.val(par.y.center())};
1195 
1196  for (auto& e : evts) {
1197  if (!e.isSig) continue;
1198 
1199  double tDiff = (e.t - par.x.val(par.x.center()));
1200  //double tDiff = e.t;
1201 
1202 
1203  double d1 = getCorrD(e.mu0, e.t, par);
1204  double d2 = getCorrD(e.mu1, e.t, par);
1205 
1206  double d1TX = getCorrD(e.mu0, e.t, parNoTx);
1207  double d2TX = getCorrD(e.mu1, e.t, parNoTx);
1208 
1209  profResTx->Fill(tDiff * sin(e.mu0.phi0), d1TX);
1210  profResTx->Fill(tDiff * sin(e.mu1.phi0), d2TX);
1211 
1212  double d1TY = getCorrD(e.mu0, e.t, parNoTy);
1213  double d2TY = getCorrD(e.mu1, e.t, parNoTy);
1214 
1215  profResTy->Fill(-tDiff * cos(e.mu0.phi0), d1TY);
1216  profResTy->Fill(-tDiff * cos(e.mu1.phi0), d2TY);
1217 
1218 
1219  profRes->Fill(tDiff * sin(e.mu0.phi0), d1);
1220  profRes->Fill(tDiff * sin(e.mu1.phi0), d2);
1221 
1222 
1223  }
1224 
1225  TCanvas* cX = new TCanvas(rn(), "");
1226  gStyle->SetOptStat(0);
1227  profResTx->Draw();
1228 
1229 
1230  TF1* f = new TF1(rn(), "[0]*x", -1, 1);
1231  f->SetParameter(0, (par.x.val(1) - par.x.val(0)));
1232  f->Draw("same");
1233  B2INFO("Table value " << par.x.val(1) - par.x.val(0));
1234 
1235  cX->SaveAs(fName + "_tX.pdf");
1236 
1237 
1238  }
1239 
1240 
1241 
1242 
1243 
1245  void plotSpotZpositionPull(const vector<Event>& evts, const SpotParam& par, TString fName, double cut = 1000)
1246  {
1247  TH1D* hPull = new TH1D(rn(), "", 200, -2000, 2000);
1248 
1249  for (auto& e : evts) {
1250  if (!e.isSig) continue;
1251 
1252  double z0 = getCorrZ(e.mu0, e.t, par);
1253  double z1 = getCorrZ(e.mu1, e.t, par);
1254 
1255  hPull->Fill(z0);
1256  hPull->Fill(z1);
1257  }
1258 
1259  gStyle->SetOptStat(2210);
1260  TCanvas* c = new TCanvas(rn(), "");
1261  hPull->Draw();
1262 
1263  hPull->GetXaxis()->SetTitle("pull [#mum]");
1264  hPull->GetYaxis()->SetTitle("#tracks");
1265 
1266  TLine* l = new TLine;
1267  l->SetLineColor(kRed);
1268  l->DrawLine(-cut, 0, -cut, 500);
1269  l->DrawLine(+cut, 0, +cut, 500);
1270 
1271  c->SaveAs(fName + ".pdf");
1272  }
1273 
1274 
1276  void removeSpotPositionOutliers(vector<Event>& evts, const SpotParam& par, double cut = 70)
1277  {
1278  int nRem = 0;
1279  int nAll = 0;
1280  for (auto& e : evts) {
1281  if (!e.isSig) continue;
1282 
1283  double d0 = getCorrD(e.mu0, e.t, par);
1284  double d1 = getCorrD(e.mu1, e.t, par);
1285 
1286  e.isSig = abs(d0) < cut && abs(d1) < cut;
1287  nRem += !e.isSig;
1288  ++nAll;
1289  }
1290  B2INFO("Removed fraction Position " << nRem / (nAll + 0.));
1291  }
1292 
1293 
1295  void removeSpotZpositionOutliers(vector<Event>& evts, const SpotParam& par, double cut = 1000)
1296  {
1297  int nRem = 0;
1298  int nAll = 0;
1299  for (auto& e : evts) {
1300  if (!e.isSig) continue;
1301 
1302  double z0 = getCorrZ(e.mu0, e.t, par);
1303  double z1 = getCorrZ(e.mu1, e.t, par);
1304 
1305  e.isSig = abs(z0) < cut && abs(z1) < cut;
1306  nRem += !e.isSig;
1307  ++nAll;
1308  }
1309  B2INFO("Removed fraction Position " << nRem / (nAll + 0.));
1310  }
1311 
1312 
1313 
1315  vector<vector<double>> fillSplineBasesLinear(const vector<Event>& evts, vector<double> spl,
1316  std::function<double(Track, double)> fun)
1317  {
1318  int n = spl.size(); //number of params
1319  if (n == 0 || (n == 2 && spl[0] > spl[1]))
1320  n = 1;
1321 
1322  vector<vector<double>> vecs(n);
1323 
1324  if (n == 1) { //no time dependence
1325  for (const auto& e : evts) {
1326  for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1327  vecs[0].push_back(1 * fun(e.mu0, e.t));
1328  vecs[0].push_back(1 * fun(e.mu1, e.t));
1329  }
1330  }
1331  } else {
1332  for (int k = 0; k < n; ++k) {
1333  double xCnt = spl[k];
1334  double xLow = (k == 0) ? spl[0] : spl[k - 1];
1335  double xHigh = (k == n - 1) ? spl[n - 1] : spl[k + 1];
1336 
1337  for (const auto& e : evts) {
1338  double x = e.t;
1339  double v = 0;
1340  if (xLow <= x && x < xCnt)
1341  v = (x - xLow) / (xCnt - xLow);
1342  else if (xCnt < x && x <= xHigh)
1343  v = (xHigh - x) / (xHigh - xCnt);
1344 
1345 
1346  for (int i = 0; i < e.nBootStrap * e.isSig; ++i) {
1347  vecs[k].push_back(v * fun(e.mu0, e.t));
1348  vecs[k].push_back(v * fun(e.mu1, e.t));
1349  }
1350  }
1351  }
1352  }
1353 
1354  return vecs;
1355  }
1356 
1357 
1359  vector<vector<double>> fillSplineBasesZero(const vector<Event>& evts, vector<double> spl, std::function<double(Track, double)> fun)
1360  {
1361  int n = spl.size() + 1; //number of params
1362 
1363  vector<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(max(e1, e2), 2);
1421  sum += d * step;
1422  }
1423  return sum;
1424  }
1425 
1427  double fitSpotZwidth(const vector<Event>& evts, const SpotParam& spotPar, const vector<double>& sizesXY)
1428  {
1429 
1430  vector<double> dataVec;
1431  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  vector<double> pars, err2;
1451  double err2Mean, err2press, err2pressErr;
1452  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 vector<Event>& evts, const vector<double>& splX, const vector<double>& splY,
1463  const vector<double>& splKX, const vector<double>& splKY)
1464  {
1465  vector<vector<double>> basesX = fillSplineBasesZero(evts, splX, [](Track tr, double) {return sin(tr.phi0);});
1466  vector<vector<double>> basesY = fillSplineBasesZero(evts, splY, [](Track tr, double) {return -cos(tr.phi0);});
1467 
1468  vector<vector<double>> basesKX = fillSplineBasesZero(evts, splKX, [](Track tr, double) {return sin(tr.phi0) * tr.z0;});
1469  vector<vector<double>> basesKY = fillSplineBasesZero(evts, splKY, [](Track tr, double) {return -cos(tr.phi0) * tr.z0;});
1470 
1471 
1472  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  vector<vector<double>> allVecs = merge({basesX, basesY, basesKX, basesKY});
1481 
1482  MatrixXd A = vecs2mat(allVecs);
1483 
1484 
1485  VectorXd vData = vec2vec(dataVec);
1486 
1487  vector<double> pars(A.cols()), err2(A.cols());
1488  double err2Mean, err2press, err2pressErr;
1489  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 vector<Event>& evts, const vector<double>& splX, const vector<double>& splY,
1497  const vector<double>& splKX, const vector<double>& splKY, const SpotParam& spotPars)
1498  {
1499  vector<vector<double>> basesX = fillSplineBasesZero(evts, splX, [](Track tr, double) {return sin(tr.phi0);});
1500  vector<vector<double>> basesY = fillSplineBasesZero(evts, splY, [](Track tr, double) {return -cos(tr.phi0);});
1501 
1502  vector<vector<double>> basesKX = fillSplineBasesZero(evts, splKX, [ = ](Track tr, double t) {return sin(tr.phi0) * (getZIPest(tr, t, spotPars) - spotPars.z.val(t));});
1503  vector<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  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  vector<vector<double>> allVecs = merge({basesX, basesY, basesKX, basesKY});
1515 
1516  MatrixXd A = vecs2mat(allVecs);
1517 
1518 
1519  VectorXd vData = vec2vec(dataVec);
1520 
1521  vector<double> pars(A.cols()), err2(A.cols());
1522  double err2Mean, err2press, err2pressErr;
1523  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 vector<Event>& evts, const vector<double>& splX, const vector<double>& splY)
1537  {
1538  vector<vector<double>> basesX = fillSplineBasesZero(evts, splX, [](Track tr, double) {return sin(tr.phi0);});
1539  vector<vector<double>> basesY = fillSplineBasesZero(evts, splY, [](Track tr, double) {return -cos(tr.phi0);});
1540 
1541  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  vector<vector<double>> allVecs = merge({basesX, basesY});
1550 
1551  MatrixXd A = vecs2mat(allVecs);
1552 
1553  VectorXd vData = vec2vec(dataVec);
1554 
1555  vector<double> pars(A.cols()), err2(A.cols());
1556  double err2Mean, err2press, err2pressErr;
1557  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 vector<Event>& evts, const vector<double>& splX, const vector<double>& splY,
1570  const vector<double>& splKX, const vector<double>& splKY,
1571  const vector<double>& splZ)
1572  {
1573  vector<vector<double>> basesX = fillSplineBasesZero(evts, splX, [](Track tr, double) {return -tr.tanlambda * cos(tr.phi0);});
1574  vector<vector<double>> basesY = fillSplineBasesZero(evts, splY, [](Track tr, double) {return -tr.tanlambda * sin(tr.phi0);});
1575 
1576  vector<vector<double>> basesKX = fillSplineBasesZero(evts, splKX, [](Track tr, double) {return -tr.z0 * tr.tanlambda * cos(tr.phi0);});
1577  vector<vector<double>> basesKY = fillSplineBasesZero(evts, splKY, [](Track tr, double) {return -tr.z0 * tr.tanlambda * sin(tr.phi0);});
1578 
1579  vector<vector<double>> basesZ = fillSplineBasesZero(evts, splZ, [](Track , double) {return 1;});
1580 
1581 
1582  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  vector<vector<double>> allVecs = merge({basesX, basesY, basesKX, basesKY, basesZ});
1591 
1592  MatrixXd A = vecs2mat(allVecs);
1593 
1594  VectorXd vData = vec2vec(dataVec);
1595 
1596  vector<double> pars(A.cols()), err2(A.cols());
1597  double err2Mean, err2press, err2pressErr;
1598  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 vector<Event>& evts, const vector<double>& splZ, const SpotParam& spotPars)
1608  {
1609  vector<vector<double>> basesZ = fillSplineBasesZero(evts, splZ, [](Track, double) {return 1;});
1610 
1611  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  vector<double> pars(A.cols()), err2(A.cols());
1626  double err2Mean, err2press, err2pressErr;
1627  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  vector<double> fitSpotWidthCMS(const vector<Event>& evts, const SpotParam& spotPar)
1643  {
1644 
1645  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 vector<Event>& evts, const SpotParam& spotPar, const 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 vector<Event>& evts, const SpotParam& spotPar, const 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 vector<Event>& evts, const SpotParam& par, const 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 vector<Event>& evts, const SpotParam& par, const 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(vector<Event>& evts, const SpotParam& spotPar, const vector<double>& sizesXY, double cut = 1500)
1867  {
1868 
1869  int nRem = 0;
1870  int nAll = 0;
1871  for (auto& e : evts) {
1872  if (!e.isSig) continue;
1873 
1874  double d0 = getCorrD(e.mu0, e.t, spotPar);
1875  double d1 = getCorrD(e.mu1, e.t, spotPar);
1876  double d12Th = getD12th(e, sizesXY);
1877 
1878  e.isSig = abs(d0 * d1 - d12Th) < cut;
1879  nRem += !e.isSig;
1880  ++nAll;
1881  }
1882  B2INFO("Removed fraction Size " << nRem / (nAll + 0.));
1883  }
1884 
1885 
1887  void removeSpotSizeZOutliers(vector<Event>& evts, const SpotParam& spotPar, const vector<double>& sizesXY, double sizeZZ,
1888  double cut = 150000)
1889  {
1890 
1891  int nRem = 0;
1892  int nAll = 0;
1893  for (auto& e : evts) {
1894  if (!e.isSig) continue;
1895 
1896  double z0 = getCorrZ(e.mu0, e.t, spotPar);
1897  double z1 = getCorrZ(e.mu1, e.t, spotPar);
1898 
1899  double corr = getZ12th(e, sizesXY);
1900  double res = z0 * z1 - corr - sizeZZ;
1901 
1902 
1903  e.isSig = abs(res) < cut;
1904  nRem += !e.isSig;
1905  ++nAll;
1906  }
1907  B2INFO("Removed fraction Size " << nRem / (nAll + 0.));
1908  }
1909 
1910 
1912  MatrixXd toMat(TRotation rot)
1913  {
1914  MatrixXd rotM(3, 3);
1915  rotM(0, 0) = rot.XX();
1916  rotM(0, 1) = rot.XY();
1917  rotM(0, 2) = rot.XZ();
1918  rotM(1, 0) = rot.YX();
1919  rotM(1, 1) = rot.YY();
1920  rotM(1, 2) = rot.YZ();
1921  rotM(2, 0) = rot.ZX();
1922  rotM(2, 1) = rot.ZY();
1923  rotM(2, 2) = rot.ZZ();
1924 
1925  return rotM;
1926  }
1927 
1928 
1936  MatrixXd getRotatedSizeMatrix(vector<double> xySize, double zzSize, double kX, double kY)
1937  {
1938  TRotation rot; // rotation moving eZ=(0,0,1) to (kX, kY, 1)
1939  rot.RotateX(-kY); //x-rot
1940  rot.RotateY(+kX); //y-rot
1941 
1942  MatrixXd rotM = toMat(rot);
1943  MatrixXd rotMT = rotM.transpose();
1944 
1945  Matrix3d eigenMat = Matrix3d::Zero(); //z-rot included in eigenMat
1946  eigenMat(0, 0) = xySize[0];
1947  eigenMat(1, 1) = xySize[1];
1948  eigenMat(0, 1) = xySize[2];
1949  eigenMat(1, 0) = xySize[2];
1950  eigenMat(2, 2) = zzSize;
1951 
1952  return (rotM * eigenMat * rotMT);
1953  }
1954 
1955 
1956 
1957 
1958 
1959 
1960 
1961 // Returns tuple with the beamspot parameters
1962  tuple<vector<VectorXd>, vector<MatrixXd>, MatrixXd> runBeamSpotAnalysis(vector<Event> evts,
1963  const vector<double>& splitPoints)
1964  {
1965  const double xyPosLimit = 70; //um
1966  const double xySize2Limit = pow(40, 2); //um^2
1967  const double zPosLimit = 1200; //um
1968 
1969 
1970  vector<double> indX = splitPoints;
1971  vector<double> indY = splitPoints;
1972  vector<double> indZ = splitPoints;
1973 
1974  //no time dependence, as for beam size
1975  vector<double> indKX = {};
1976  vector<double> indKY = {};
1977 
1978  UnknownPars allPars;
1979  const int kPlot = -1; //do plots for index kPlot
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  if (k == kPlot) {
1989  plotSpotPositionFit(evts, resTemp, "positionFitSimpe");
1990  plotSpotPositionPull(evts, resTemp, "pullsPositionSimple", xyPosLimit);
1991  }
1992  removeSpotPositionOutliers(evts, resTemp, xyPosLimit);
1993 
1994  //simple XY pos fit (with outliers removed)
1995  auto resFin = fitSpotPositionSplines(evts, indX, indY);
1996  if (k == kPlot) {
1997  plotSpotPositionFit(evts, resFin, "positionFitSimpleC");
1998  plotSpotPositionPull(evts, resFin, "pullsPositionSimpleC", xyPosLimit);
1999  plotXYtimeDep(evts, resFin, "simplePosTimeDep");
2000  }
2001 
2002  //Z position fit
2003  auto resZmy = fitZpositionSplinesSimple(evts, indZ, resFin);
2004  if (k == kPlot) {
2005  plotSpotZPositionFit(evts, resZmy, "positionFitSimpleZ");
2006  plotSpotZpositionPull(evts, resZmy, "zPositionPull", zPosLimit);
2007  }
2008 
2009  removeSpotZpositionOutliers(evts, resZmy, zPosLimit);
2010 
2011  //Z position fit (with outliers removed)
2012  resZmy = fitZpositionSplinesSimple(evts, indZ, resZmy);
2013 
2014 
2015  //complete XY pos fit
2016  auto resNew = fitSpotPositionSplines(evts, indX, indY, indKX, indKY, resZmy);
2017  if (k == kPlot) {
2018  plotSpotPositionFit(evts, resNew, "positionFitFull");
2019  plotKxKyFit(evts, resNew, "slopes");
2020  }
2021 
2022  //Z position fit (iteration 2)
2023  resZmy = fitZpositionSplinesSimple(evts, indZ, resNew);
2024  if (k == kPlot) plotSpotZPositionFit(evts, resZmy, "positionFitSimpleZLast");
2025 
2026 
2027  //complete XY pos fit (iteration 2)
2028  resNew = fitSpotPositionSplines(evts, indX, indY, indKX, indKY, resZmy);
2029 
2030  //XYZ pos fits (iteration 3)
2031  resZmy = fitZpositionSplinesSimple(evts, indZ, resNew);
2032  resNew = fitSpotPositionSplines(evts, indX, indY, indKX, indKY, resZmy);
2033 
2034 
2035  // fit of XY sizes (original + with outliers removed)
2036  auto vecXY = fitSpotWidthCMS(evts, resNew);
2037  if (k == kPlot) plotSpotSizePull(evts, resNew, vecXY);
2038  removeSpotSizeOutliers(evts, resNew, vecXY, xySize2Limit);
2039  vecXY = fitSpotWidthCMS(evts, resNew);
2040  if (k == kPlot) plotSpotSizeFit(evts, resNew, vecXY);
2041 
2042 
2043  // fit of Z size
2044  double sizeZZ = fitSpotZwidth(evts, resNew, vecXY);
2045 
2046  if (k == kPlot) {
2047  plotSpotZSizeFit(evts, resNew, vecXY, sizeZZ);
2048  plotSpotSizeZPull(evts, resNew, vecXY, sizeZZ);
2049  }
2050 
2051  //removeSpotSizeZOutliers(evts, resNew, vecXY, sizeZZ, 150000);
2052  //sizeZZ = fitSpotZwidth(evts, resNew, vecXY);
2053 
2054  allPars.add(resNew, sqrtS(vecXY[0]), sqrtS(vecXY[1]), sqrtS(vecXY[2]), sqrtS(sizeZZ));
2055 
2056  }
2057 
2058  //allPars.printStat();
2059 
2060  vector<VectorXd> vtxPos;
2061  vector<MatrixXd> vtxErr;
2062  MatrixXd sizeMat;
2063 
2064  allPars.getOutput(vtxPos, vtxErr, sizeMat);
2065 
2066  return make_tuple(vtxPos, vtxErr, sizeMat);
2067  }
2068 
2069 }
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:25
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:73
Belle2::BeamSpotCalib::SpotParam::print
void print()
Print BeamSpot parameters.
Definition: BeamSpotStandAlone.cc:83
Belle2::BeamSpotCalib::SpotParam::z
Spline z
spline for BS center position as a function of time (z coordinate)
Definition: BeamSpotStandAlone.cc:76
Belle2::BeamSpotCalib::SpotParam::x
Spline x
spline for BS center position as a function of time (x coordinate)
Definition: BeamSpotStandAlone.cc:74
Belle2::BeamSpotCalib::SpotParam::kX
Spline kX
spline for BS angle in the xz plane as a function of time
Definition: BeamSpotStandAlone.cc:77
Belle2::BeamSpotCalib::SpotParam::SpotParam
SpotParam(const vector< double > &vals, const vector< double > &errs, const vector< vector< double >> &spls, int order=0)
Constructor based output of the linear regression, assuming zero-order splines vals,...
Definition: BeamSpotStandAlone.cc:101
Belle2::BeamSpotCalib::SpotParam::kY
Spline kY
spline for BS angle in the yz plane as a function of time
Definition: BeamSpotStandAlone.cc:78
Belle2::BeamSpotCalib::SpotParam::y
Spline y
spline for BS center position as a function of time (y coordinate)
Definition: BeamSpotStandAlone.cc:75
Belle2::BeamSpotCalib::UnknowSpline
Spline with uncertainty obtained from the boot-strap replicas.
Definition: BeamSpotStandAlone.cc:153
Belle2::BeamSpotCalib::UnknowSpline::spls
vector< Spline > spls
vector with replicas
Definition: BeamSpotStandAlone.cc:154
Belle2::BeamSpotCalib::UnknowSpline::add
void add(Spline spl)
add boot-strap replica
Definition: BeamSpotStandAlone.cc:155
Belle2::BeamSpotCalib::UnknowSpline::getMeanSigma
Spline getMeanSigma()
Get mean and 1-sigma errors of the spline values.
Definition: BeamSpotStandAlone.cc:158
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:188
Belle2::BeamSpotCalib::UnknowVar
variable with uncertainty from boot-strap replicas
Definition: BeamSpotStandAlone.cc:215
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:241
Belle2::BeamSpotCalib::UnknowVar::add
void add(double x)
add value to the replicas
Definition: BeamSpotStandAlone.cc:217
Belle2::BeamSpotCalib::UnknowVar::getStats
vector< double > getStats()
Get basic stats.
Definition: BeamSpotStandAlone.cc:260
Belle2::BeamSpotCalib::UnknowVar::vars
vector< double > vars
vector of variable values for all replicas
Definition: BeamSpotStandAlone.cc:216
Belle2::BeamSpotCalib::UnknowVar::getSigma
double getSigma()
Get standard deviation.
Definition: BeamSpotStandAlone.cc:227
Belle2::BeamSpotCalib::UnknowVar::printStat
void printStat(TString n)
Print variable of name n with stat-info.
Definition: BeamSpotStandAlone.cc:252
Belle2::BeamSpotCalib::UnknownPars
structure including all variables of interest with uncertainties from boot-strap
Definition: BeamSpotStandAlone.cc:335
Belle2::BeamSpotCalib::UnknownPars::setBranchVal
void setBranchVal(TTree *T, vector< double > *vec, TString n)
save bootstrap variable to TTree
Definition: BeamSpotStandAlone.cc:480
Belle2::BeamSpotCalib::UnknownPars::kY
UnknowSpline kY
BS angle in yz plane.
Definition: BeamSpotStandAlone.cc:339
Belle2::BeamSpotCalib::UnknownPars::crossAngle
UnknowVar crossAngle
derived value of the crossing angle of the HER & LER beams
Definition: BeamSpotStandAlone.cc:351
Belle2::BeamSpotCalib::UnknownPars::matZZ
UnknowVar matZZ
ZZ element of BS size cov matrix.
Definition: BeamSpotStandAlone.cc:355
Belle2::BeamSpotCalib::UnknownPars::sizeMin
UnknowVar sizeMin
smallest eigenvalue of the BS size cov matrix (similar to sizeY)
Definition: BeamSpotStandAlone.cc:345
Belle2::BeamSpotCalib::UnknownPars::printStat
void printStat()
Print interesting statistics from boot-strap.
Definition: BeamSpotStandAlone.cc:409
Belle2::BeamSpotCalib::UnknownPars::matYY
UnknowVar matYY
YY element of BS size cov matrix.
Definition: BeamSpotStandAlone.cc:354
Belle2::BeamSpotCalib::UnknownPars::sizeZ
UnknowVar sizeZ
BS size in z direction.
Definition: BeamSpotStandAlone.cc:350
Belle2::BeamSpotCalib::UnknownPars::z
UnknowSpline z
BS position (z coordinate)
Definition: BeamSpotStandAlone.cc:340
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:361
Belle2::BeamSpotCalib::UnknownPars::getOutput
void getOutput(vector< VectorXd > &vtxPos, vector< MatrixXd > &vtxErr, MatrixXd &sizeMat)
get output in Belle2-like format
Definition: BeamSpotStandAlone.cc:438
Belle2::BeamSpotCalib::UnknownPars::sizeY
UnknowVar sizeY
BS size in y direction.
Definition: BeamSpotStandAlone.cc:343
Belle2::BeamSpotCalib::UnknownPars::matXZ
UnknowVar matXZ
XZ element of BS size cov matrix.
Definition: BeamSpotStandAlone.cc:357
Belle2::BeamSpotCalib::UnknownPars::matXY
UnknowVar matXY
XY element of BS size cov matrix.
Definition: BeamSpotStandAlone.cc:356
Belle2::BeamSpotCalib::UnknownPars::y
UnknowSpline y
BS position (y coordinate)
Definition: BeamSpotStandAlone.cc:337
Belle2::BeamSpotCalib::UnknownPars::sizeX
UnknowVar sizeX
BS size in x direction.
Definition: BeamSpotStandAlone.cc:342
Belle2::BeamSpotCalib::UnknownPars::matXX
UnknowVar matXX
XX element of BS size cov matrix.
Definition: BeamSpotStandAlone.cc:353
Belle2::BeamSpotCalib::UnknownPars::save2tree
void save2tree(TString fName)
save everything to TTree
Definition: BeamSpotStandAlone.cc:498
Belle2::BeamSpotCalib::UnknownPars::sizeMax
UnknowVar sizeMax
middle eigenvalue of the BS size cov matrix (similar to sizeX)
Definition: BeamSpotStandAlone.cc:346
Belle2::BeamSpotCalib::UnknownPars::setBranchSpline
void setBranchSpline(TTree *T, Spline *spl, TString n)
save bootstrap spline to TTree
Definition: BeamSpotStandAlone.cc:488
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:344
Belle2::BeamSpotCalib::UnknownPars::x
UnknowSpline x
BS position (x coordinate)
Definition: BeamSpotStandAlone.cc:336
Belle2::BeamSpotCalib::UnknownPars::xyAngle
UnknowVar xyAngle
angle of the BS in xy plane when z' is aligned with z
Definition: BeamSpotStandAlone.cc:347
Belle2::BeamSpotCalib::UnknownPars::matYZ
UnknowVar matYZ
YZ element of BS size cov matrix.
Definition: BeamSpotStandAlone.cc:358
Belle2::BeamSpotCalib::UnknownPars::kX
UnknowSpline kX
BS angle in xz plane.
Definition: BeamSpotStandAlone.cc:338
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