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