Belle II Software  release-05-01-25
GFGbl.cc
1 //-*-mode: C++; c-basic-offset: 2; -*-
2 /* Copyright 2013
3  * Authors: Sergey Yashchenko and Tadeas Bilka
4  *
5  * This is an interface to General Broken Lines
6  *
7  * Version: 5 (Tadeas)
8  * - several bug-fixes:
9  * - Scatterers at bad points
10  * - Jacobians at a point before they should be (code reorganized)
11  * - Change of sign of residuals
12  * Version: 4 (Tadeas)
13  * Fixed calculation of equvivalent scatterers (solution by C. Kleinwort)
14  * Now a scatterer is inserted at each measurement (except last) and between each two measurements.
15  * TrueHits/Clusters. Ghost (1D) hits ignored. With or without magnetic field.
16  * Version: 3 (Tadeas)
17  * This version now supports both TrueHits and Clusters for VXD.
18  * It can be used for arbitrary material distribution between
19  * measurements. Moments of scattering distribution are computed
20  * and translated into two equivalent thin GBL scatterers placed
21  * at computed positions between measurement points.
22  * Version: 2 ... never published (Tadeas)
23  * Scatterer at each boundary (tooo many scatterers). TrueHits/Clusters. Without global der.&MP2 output.
24  * Version: 1 (Sergey & Tadeas)
25  * Scatterers at measurement planes. TrueHits
26  * Version 0: (Sergey)
27  * Without scatterers. Genfit 1.
28  *
29  * This file is part of GENFIT.
30  *
31  * GENFIT is free software: you can redistribute it and/or modify
32  * it under the terms of the GNU Lesser General Public License as published
33  * by the Free Software Foundation, either version 3 of the License, or
34  * (at your option) any later version.
35  *
36  * GENFIT is distributed in the hope that it will be useful,
37  * but WITHOUT ANY WARRANTY; without even the implied warranty of
38  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
39  * GNU Lesser General Public License for more details.
40  *
41  * You should have received a copy of the GNU Lesser General Public License
42  * along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
43  */
44 
45 #include "GFGbl.h"
46 #include "GblTrajectory.h"
47 #include "GblPoint.h"
48 #include "MyDebugTools.h"
49 
50 #include "AbsMeasurement.h"
51 #include "PlanarMeasurement.h"
52 #include "KalmanFitterInfo.h"
53 
54 #include "Track.h"
55 #include <TFile.h>
56 #include <TH1F.h>
57 #include <TTree.h>
58 #include <string>
59 #include <list>
60 #include <FieldManager.h>
61 #include <HMatrixU.h>
62 #include <HMatrixV.h>
63 #include <Math/SMatrix.h>
64 #include <TMatrixD.h>
65 #include <TVectorDfwd.h>
66 #include <TMatrixT.h>
67 
68 #include <TVector3.h>
69 
70 //#define DEBUG
71 //#define OUTPUT
72 
73 
74 #ifdef DEBUG
75 //ofstream debug("gbl.debug");
76 #endif
77 
78 #ifdef OUTPUT
79 
80 std::string rootFileName = "gbl.root";
81 
82 
83 TFile* diag;
84 TH1F* resHistosU[12];
85 TH1F* resHistosV[12];
86 TH1F* mhistosU[12];
87 TH1F* mhistosV[12];
88 TH1F* ghistosU[12];
89 TH1F* ghistosV[12];
90 TH1F* downWeightsHistosU[12];
91 TH1F* downWeightsHistosV[12];
92 TH1F* localPar1[12];
93 TH1F* localPar2[12];
94 TH1F* localPar3[12];
95 TH1F* localPar4[12];
96 TH1F* localPar5[12];
97 TH1F* localCov1[12];
98 TH1F* localCov2[12];
99 TH1F* localCov3[12];
100 TH1F* localCov4[12];
101 TH1F* localCov5[12];
102 TH1F* localCov12[12];
103 TH1F* localCov13[12];
104 TH1F* localCov14[12];
105 TH1F* localCov15[12];
106 TH1F* localCov23[12];
107 TH1F* localCov24[12];
108 TH1F* localCov25[12];
109 TH1F* localCov34[12];
110 TH1F* localCov35[12];
111 TH1F* localCov45[12];
112 TH1I* fitResHisto;
113 TH1I* ndfHisto;
114 TH1F* chi2OndfHisto;
115 TH1F* pValueHisto;
116 TH1I* stats;
117 
118 
119 
120 bool writeHistoDataForLabel(double label, TVectorD res, TVectorD measErr, TVectorD resErr, TVectorD downWeights, TVectorD localPar, TMatrixDSym localCov)
121 {
122  if (label < 1.) return false;
123 
124  unsigned int id = floor(label);
125  // skip segment (5 bits)
126  id = id >> 5;
127  unsigned int sensor = id & 7;
128  id = id >> 3;
129  unsigned int ladder = id & 31;
130  id = id >> 5;
131  unsigned int layer = id & 7;
132  if (layer == 7 && ladder == 2) {
133  label = sensor;
134  } else if (layer == 7 && ladder == 3) {
135  label = sensor + 9 - 3;
136  } else {
137  label = layer + 3;
138  }
139 
140  if (label > 12.) return false;
141 
142  int i = int(label);
143 
144  #ifdef OUTPUT
145  resHistosU[i - 1]->Fill(res[0]);
146  resHistosV[i - 1]->Fill(res[1]);
147  mhistosU[i - 1]->Fill(res[0] / measErr[0]);
148  mhistosV[i - 1]->Fill(res[1] / measErr[1]);
149  ghistosU[i - 1]->Fill(res[0] / resErr[0]);
150  ghistosV[i - 1]->Fill(res[1] / resErr[1]);
151  downWeightsHistosU[i - 1]->Fill(downWeights[0]);
152  downWeightsHistosV[i - 1]->Fill(downWeights[1]);
153  localPar1[i - 1]->Fill(localPar(0));
154  localPar2[i - 1]->Fill(localPar(1));
155  localPar3[i - 1]->Fill(localPar(2));
156  localPar4[i - 1]->Fill(localPar(3));
157  localPar5[i - 1]->Fill(localPar(4));
158  #endif
159 
160 
161  return true;
162 }
163 #endif
164 
165 // Millepede Binary File for output of GBL trajectories for alignment
166 gbl::MilleBinary* milleFile;
167 // Minimum scattering sigma (will be squared and inverted...)
168 const double scatEpsilon = 1.e-8;
169 
170 
171 using namespace gbl;
172 using namespace std;
173 using namespace genfit;
174 
175 GFGbl::GFGbl() :
176 AbsFitter(), m_milleFileName("millefile.dat"), m_gblInternalIterations("THC"), m_pValueCut(0.), m_minNdf(1),
177 m_chi2Cut(0.),
178 m_enableScatterers(true),
179 m_enableIntermediateScatterer(true)
180 {
181 
182 }
183 
184 void GFGbl::beginRun()
185 {
186  milleFile = new MilleBinary(m_milleFileName);
187 
188  #ifdef OUTPUT
189  diag = new TFile(rootFileName.c_str(), "RECREATE");
190  char name[20];
191 
192  for (int i = 0; i < 12; i++) {
193  sprintf(name, "res_u_%i", i + 1);
194  resHistosU[i] = new TH1F(name, "Residual (U)", 1000, -0.1, 0.1);
195  sprintf(name, "res_v_%i", i + 1);
196  resHistosV[i] = new TH1F(name, "Residual (V)", 1000, -0.1, 0.1);
197  sprintf(name, "meas_pull_u_%i", i + 1);
198  mhistosU[i] = new TH1F(name, "Res/Meas.Err. (U)", 1000, -20., 20.);
199  sprintf(name, "meas_pull_v_%i", i + 1);
200  mhistosV[i] = new TH1F(name, "Res/Meas.Err. (V)", 1000, -20., 20.);
201  sprintf(name, "pull_u_%i", i + 1);
202  ghistosU[i] = new TH1F(name, "Res/Res.Err. (U)", 1000, -20., 20.);
203  sprintf(name, "pull_v_%i", i + 1);
204  ghistosV[i] = new TH1F(name, "Res/Res.Err. (V)", 1000, -20., 20.);
205  sprintf(name, "downWeights_u_%i", i + 1);
206  downWeightsHistosU[i] = new TH1F(name, "Down-weights (U)", 1000, 0., 1.);
207  sprintf(name, "downWeights_v_%i", i + 1);
208  downWeightsHistosV[i] = new TH1F(name, "Down-weights (V)", 1000, 0., 1.);
209  sprintf(name, "localPar1_%i", i + 1);
210 
211  localPar1[i] = new TH1F(name, "Residual (U)", 1000, -0.1, 0.1);
212  sprintf(name, "localPar2_%i", i + 1);
213  localPar2[i] = new TH1F(name, "Residual (U)", 1000, -0.1, 0.1);
214  sprintf(name, "localPar3_%i", i + 1);
215  localPar3[i] = new TH1F(name, "Residual (U)", 1000, -0.1, 0.1);
216  sprintf(name, "localPar4_%i", i + 1);
217  localPar4[i] = new TH1F(name, "Residual (U)", 1000, -0.1, 0.1);
218  sprintf(name, "localPar5_%i", i + 1);
219  localPar5[i] = new TH1F(name, "Residual (U)", 1000, -0.1, 0.1);
220  }
221  fitResHisto = new TH1I("fit_result", "GBL Fit Result", 21, -1, 20);
222  ndfHisto = new TH1I("ndf", "GBL Track NDF", 41, -1, 40);
223  chi2OndfHisto = new TH1F("chi2_ndf", "Track Chi2/NDF", 100, 0., 10.);
224  pValueHisto = new TH1F("p_value", "Track P-value", 100, 0., 1.);
225 
226  stats = new TH1I("stats", "0: NDF>0 | 1: fTel&VXD | 2: all | 3: VXD | 4: SVD | 5: all - PXD | 6: fTel&SVD | 7: bTel", 10, 0, 10);
227 
228  #endif
229 }
230 
231 void GFGbl::endRun()
232 {
233  #ifdef OUTPUT
234  diag->cd();
235  diag->Write();
236  diag->Close();
237  #endif
238  // This is needed to close the file before alignment starts
239  if (milleFile)
240  delete milleFile;
241 }
242 
263 void getScattererFromMatList(double& length, double& theta, double& s, double& ds, const double p, const double mass, const double charge, const std::vector<MatStep>& steps)
264 {
265  theta = 0.; s = 0.; ds = 0.;
266  if (steps.empty()) return;
267 
268  // sum of step lengths
269  double len = 0.;
270  // normalization
271  double sumxx = 0.;
272  // first moment (non-normalized)
273  double sumx2x2 = 0.;
274  // (part of) second moment / variance (non-normalized)
275  double sumx3x3 = 0.;
276 
277  double xmin = 0.;
278  double xmax = 0.;
279 
280  for (unsigned int i = 0; i < steps.size(); i++) {
281  const MatStep step = steps.at(i);
282  // inverse of material radiation length ... (in 1/cm) ... "density of scattering"
283  double rho = 1. / step.material_.radiationLength;
284  len += fabs(step.stepSize_);
285  xmin = xmax;
286  xmax = xmin + fabs(step.stepSize_);
287  // Compute integrals
288 
289  // integral of rho(x)
290  sumxx += rho * (xmax - xmin);
291  // integral of x*rho(x)
292  sumx2x2 += rho * (xmax * xmax - xmin * xmin) / 2.;
293  // integral of x*x*rho(x)
294  sumx3x3 += rho * (xmax * xmax * xmax - xmin * xmin * xmin) / 3.;
295  }
296  // This ensures PDG formula still gives positive results (but sumxx should be > 1e-4 for it to hold)
297  if (sumxx < 1.0e-10) return;
298  // Calculate theta from total sum of radiation length
299  // instead of summimg squares of individual deflection angle variances
300  // PDG formula:
301  double beta = p / sqrt(p * p + mass * mass);
302  theta = (0.0136 / p / beta) * fabs(charge) * sqrt(sumxx) * (1. + 0.038 * log(sumxx));
303  //theta = (0.015 / p / beta) * fabs(charge) * sqrt(sumxx);
304 
305  // track length
306  length = len;
307  // Normalization factor
308  double N = 1. / sumxx;
309  // First moment
310  s = N * sumx2x2;
311  // Square of second moment (variance)
312  // integral of (x - s)*(x - s)*rho(x)
313  double ds_2 = N * (sumx3x3 - 2. * sumx2x2 * s + sumxx * s * s);
314  ds = sqrt(ds_2);
315 
316  #ifdef DEBUG
317  //cout << "Thick scatterer parameters:" << endl;
318  //cout << "Variance of theta: " << theta << endl;
319  //cout << "Mean s : " << s << endl;
320  //cout << "Variance of s : " << ds << endl;
321 
322  #endif
323 }
324 
325 
326 void GFGbl::processTrackWithRep(Track* trk, const AbsTrackRep* rep, bool resortHits)
327 {
328  // Flag used to mark error in raw measurement combination
329  // measurement won't be considered, but scattering yes
330  bool skipMeasurement = false;
331  // Chi2 of Reference Track
332  double trkChi2 = 0.;
333  // This flag enables/disables fitting of q/p parameter in GBL
334  // It is switched off automatically if no B-field at (0,0,0) is detected.
335  bool fitQoverP = true;
336  //TODO: Use clever way to determine zero B-field
337  double Bfield = genfit::FieldManager::getInstance()->getFieldVal(TVector3(0., 0., 0.)).Mag();
338  if (!(Bfield > 0.))
339  fitQoverP = false;
340 
341  // Dimesion of repository/state
342  int dim = rep->getDim();
343  // current measurement point
344  TrackPoint* point_meas;
345  // current raw measurement
346  AbsMeasurement* raw_meas;
347 
348  // We assume no initial kinks, this will be reused several times
349  TVectorD scatResidual(2);
350  scatResidual.Zero();
351 
352  // All measurement points in ref. track
353  int npoints_meas = trk->getNumPointsWithMeasurement();
354 
355  #ifdef DEBUG
356  int npoints_all = trk->getNumPoints();
357 
358  if (resortHits)
359  cout << "WARNING: Hits resorting in GBL interface not supported." << endl;
360 
361  cout << "-------------------------------------------------------" << endl;
362  cout << " GBL processing genfit::Track " << endl;
363  cout << "-------------------------------------------------------" << endl;
364  cout << " # Ref. Track Points : " << npoints_all << endl;
365  cout << " # Meas. Points : " << npoints_meas << endl;
366 
367  #endif
368  // List of prepared GBL points for GBL trajectory construction
369  std::vector<GblPoint> listOfPoints;
370 
371  std::vector<double> listOfLayers;
372  //TODO: Add internal/external seed (from CDC) option in the future
373  // index of point with seed information (0 for none)
374  unsigned int seedLabel = 0;
375  // Seed covariance
376  // TMatrixDSym clCov(dim);
377  // Seed state
378  TMatrixDSym clSeed(dim);
379 
380  // propagation Jacobian to next point from current measurement point
381  TMatrixD jacPointToPoint(dim, dim);
382  jacPointToPoint.UnitMatrix();
383 
384  int n_gbl_points = 0;
385  int n_gbl_meas_points = 0;
386  int Ndf = 0;
387  double Chi2 = 0.;
388  double lostWeight = 0.;
389 
390  // Momentum of track at current plane
391  double trackMomMag = 0.;
392  // Charge of particle at current plane :-)
393  double particleCharge = 1.;
394 
395  for (int ipoint_meas = 0; ipoint_meas < npoints_meas; ipoint_meas++) {
396  point_meas = trk->getPointWithMeasurement(ipoint_meas);
397 
398  if (!point_meas->hasFitterInfo(rep)) {
399  cout << " ERROR: Measurement point does not have a fitter info. Track will be skipped." << endl;
400  return;
401  }
402  // Fitter info which contains Reference state and plane
403  KalmanFitterInfo* fi = dynamic_cast<KalmanFitterInfo*>(point_meas->getFitterInfo(rep));
404  if (!fi) {
405  cout << " ERROR: KalmanFitterInfo (with reference state) for measurement does not exist. Track will be skipped." << endl;
406  return;
407  }
408  // Current detector plane
409  SharedPlanePtr plane = fi->getPlane();
410  if (!fi->hasReferenceState()) {
411  cout << " ERROR: Fitter info does not contain reference state. Track will be skipped." << endl;
412  return;
413  }
414  // Reference StateOnPlane for extrapolation
415  ReferenceStateOnPlane* reference = new ReferenceStateOnPlane(*fi->getReferenceState());//(dynamic_cast<const genfit::ReferenceStateOnPlane&>(*fi->getReferenceState()));
416  // Representation state at plane
417  TVectorD state = reference->getState();
418  // track direction at plane (in global coords)
419  TVector3 trackDir = rep->getDir(*reference);
420  // track momentum vector at plane (in global coords)
421  trackMomMag = rep->getMomMag(*reference);
422  // charge of particle
423  particleCharge = rep->getCharge(*reference);
424  // mass of particle
425  double particleMass = rep->getMass(*reference);
426 
427  // Parameters of a thick scatterer between measurements
428  double trackLen = 0.;
429  double scatTheta = 0.;
430  double scatSMean = 0.;
431  double scatDeltaS = 0.;
432  // Parameters of two equivalent thin scatterers
433  double theta1 = 0.;
434  double theta2 = 0.;
435  double s1 = 0.;
436  double s2 = 0.;
437 
438  TMatrixDSym noise;
439  TVectorD deltaState;
440  // jacobian from s2 to M2
441  TMatrixD jacMeas2Scat(dim, dim);
442  jacMeas2Scat.UnitMatrix();
443 
444 
445  // Now get measurement. First have a look if we need to combine SVD clusters...
446  // Load Jacobian from previous extrapolation
447  // rep->getForwardJacobianAndNoise(jacPointToPoint, noise, deltaState);
448  // Try to get VxdId of current plane
449  int sensorId = 0;
450  PlanarMeasurement* measPlanar = dynamic_cast<PlanarMeasurement*>(point_meas->getRawMeasurement(0));
451  if (measPlanar) sensorId = measPlanar->getPlaneId();
452 
453  //WARNING: Now we only support 2D measurements. If 2 raw measurements are stored at the track
454  // point, these are checked if they correspond to "u" and "v" measurement (for SVD clusters) and these
455  // measurements are combined. SLANTED SENSORS NOT YET SUPPORTED!!!
456  if (point_meas->getRawMeasurement(0)->getDim() != 2
457  && trk->getPointWithMeasurement(ipoint_meas)->getNumRawMeasurements() == 2
458  && point_meas->getRawMeasurement(0)->getDim() == 1
459  && point_meas->getRawMeasurement(1)->getDim() == 1) {
460  AbsMeasurement* raw_measU = 0;
461  AbsMeasurement* raw_measV = 0;
462 
463  // cout << " Two 1D Measurements encountered. " << endl;
464 
465  int sensorId2 = -1;
466  PlanarMeasurement* measPlanar2 = dynamic_cast<PlanarMeasurement*>(point_meas->getRawMeasurement(0));
467  if (measPlanar2) sensorId2 = measPlanar->getPlaneId();
468 
469  // We only try to combine if at same sensor id (should be always, but who knows)
470  // otherwise ignore this point
471  if (sensorId != sensorId2) {
472  skipMeasurement = true;
473  cout << " ERROR: Incompatible sensorIDs at measurement point " << ipoint_meas << ". Measurement will be skipped." << endl;
474  }
475 
476  // We have to combine two SVD 1D Clusters at the same plane into one 2D recohit
477  AbsMeasurement* raw_meas1 = point_meas->getRawMeasurement(0);
478  AbsMeasurement* raw_meas2 = point_meas->getRawMeasurement(1);
479  // Decide which cluster is u and which v based on H-matrix
480  if (raw_meas1->constructHMatrix(rep)->isEqual(genfit::HMatrixU())
481  && raw_meas2->constructHMatrix(rep)->isEqual(genfit::HMatrixV())) {
482  // right order U, V
483  raw_measU = raw_meas1;
484  raw_measV = raw_meas2;
485  } else if (raw_meas2->constructHMatrix(rep)->isEqual(genfit::HMatrixU())
486  && raw_meas1->constructHMatrix(rep)->isEqual(genfit::HMatrixV())) {
487  // inversed order V, U
488  raw_measU = raw_meas2;
489  raw_measV = raw_meas1;
490  } else {
491  // Incompatible measurements ... skip track ... I saw this once and just before a segfault ...
492  cout << " ERROR: Incompatible 1D measurements at meas. point " << ipoint_meas << ". Track will be skipped." << endl;
493  return;
494  }
495  // Combine raw measurements
496  TVectorD _raw_coor(2);
497  _raw_coor(0) = raw_measU->getRawHitCoords()(0);
498  _raw_coor(1) = raw_measV->getRawHitCoords()(0);
499  // Combine covariance matrix
500  TMatrixDSym _raw_cov(2);
501  _raw_cov.Zero();
502  _raw_cov(0, 0) = raw_measU->getRawHitCov()(0, 0);
503  _raw_cov(1, 1) = raw_measV->getRawHitCov()(0, 0);
504  // Create new combined measurement
505  raw_meas = new PlanarMeasurement(_raw_coor, _raw_cov, raw_measU->getDetId(), raw_measU->getHitId(), point_meas);
506  } else {
507  // Default behavior
508  raw_meas = point_meas->getRawMeasurement(0);
509  }
510  //TODO: We only support 2D measurements in GBL (ot two 1D combined above)
511  if (raw_meas->getRawHitCoords().GetNoElements() != 2) {
512  skipMeasurement = true;
513  #ifdef DEBUG
514  cout << " WARNING: Measurement " << (ipoint_meas + 1) << " is not 2D. Measurement Will be skipped. " << endl;
515  #endif
516  }
517 
518  // Now we have all necessary information, so lets insert current measurement point
519  // if we don't want to skip it (e.g. ghost SVD hit ... just 1D information)
520  if (!skipMeasurement) {
521  // 2D hit coordinates
522  TVectorD raw_coor = raw_meas->getRawHitCoords();
523  // Covariance matrix of measurement
524  TMatrixDSym raw_cov = raw_meas->getRawHitCov();
525  // Projection matrix from repository state to measurement coords
526  std::unique_ptr<const AbsHMatrix> HitHMatrix(raw_meas->constructHMatrix(rep));
527  // Residual between measured position and reference track position
528  TVectorD residual = -1. * (raw_coor - HitHMatrix->Hv(state));
529 
530  trkChi2 += residual(0) * residual(0) / raw_cov(0, 0) + residual(1) * residual(1) / raw_cov(1, 1);
531 
532  // Measurement point
533  GblPoint measPoint(jacPointToPoint);
534  // Projection from local (state) coordinates to measurement coordinates (inverted)
535  // 2x2 matrix ... last block of H matrix (2 rows x 5 columns)
536  //TMatrixD proL2m = HitHMatrix->getMatrix().GetSub(0, 1, 3, 4);
537  TMatrixD proL2m(2, 2);
538  proL2m.Zero();
539  proL2m(0, 0) = HitHMatrix->getMatrix()(0, 3);
540  proL2m(0, 1) = HitHMatrix->getMatrix()(0, 4);
541  proL2m(1, 1) = HitHMatrix->getMatrix()(1, 4);
542  proL2m(1, 0) = HitHMatrix->getMatrix()(1, 3);
543  proL2m.Invert();
544  //raw_cov *= 100.;
545  //proL2m.Print();
546  measPoint.addMeasurement(proL2m, residual, raw_cov.Invert());
547 
548  //Add global derivatives to the point
549 
550  // sensor label = sensorID * 10, then last digit is label for global derivative for the sensor
551  int label = sensorId * 10;
552  // values for global derivatives
553  //TMatrixD derGlobal(2, 6);
554  TMatrixD derGlobal(2, 12);
555 
556  // labels for global derivatives
557  std::vector<int> labGlobal;
558 
559  // track direction in global coords
560  TVector3 tDir = trackDir;
561  // sensor u direction in global coords
562  TVector3 uDir = plane->getU();
563  // sensor v direction in global coords
564  TVector3 vDir = plane->getV();
565  // sensor normal direction in global coords
566  TVector3 nDir = plane->getNormal();
567  //file << sensorId << endl;
568  //outputVector(uDir, "U");
569  //outputVector(vDir, "V");
570  //outputVector(nDir, "Normal");
571  // track direction in local sensor system
572  TVector3 tLoc = TVector3(uDir.Dot(tDir), vDir.Dot(tDir), nDir.Dot(tDir));
573 
574  // track u-slope in local sensor system
575  double uSlope = tLoc[0] / tLoc[2];
576  // track v-slope in local sensor system
577  double vSlope = tLoc[1] / tLoc[2];
578 
579  // Measured track u-position in local sensor system
580  double uPos = raw_coor[0];
581  // Measured track v-position in local sensor system
582  double vPos = raw_coor[1];
583 
584  //Global derivatives for alignment in sensor local coordinates
585  derGlobal(0, 0) = 1.0;
586  derGlobal(0, 1) = 0.0;
587  derGlobal(0, 2) = - uSlope;
588  derGlobal(0, 3) = vPos * uSlope;
589  derGlobal(0, 4) = -uPos * uSlope;
590  derGlobal(0, 5) = vPos;
591 
592  derGlobal(1, 0) = 0.0;
593  derGlobal(1, 1) = 1.0;
594  derGlobal(1, 2) = - vSlope;
595  derGlobal(1, 3) = vPos * vSlope;
596  derGlobal(1, 4) = -uPos * vSlope;
597  derGlobal(1, 5) = -uPos;
598 
599  labGlobal.push_back(label + 1); // u
600  labGlobal.push_back(label + 2); // v
601  labGlobal.push_back(label + 3); // w
602  labGlobal.push_back(label + 4); // alpha
603  labGlobal.push_back(label + 5); // beta
604  labGlobal.push_back(label + 6); // gamma
605 
606  // Global derivatives for movement of whole detector system in global coordinates
607  //TODO: Usage of this requires Hierarchy Constraints to be provided to MP2
608 
609  // sensor centre position in global system
610  TVector3 detPos = plane->getO();
611  //cout << "detPos" << endl;
612  //detPos.Print();
613 
614  // global prediction from raw measurement
615  TVector3 pred = detPos + uPos * uDir + vPos * vDir;
616  //cout << "pred" << endl;
617  //pred.Print();
618 
619  double xPred = pred[0];
620  double yPred = pred[1];
621  double zPred = pred[2];
622 
623  // scalar product of sensor normal and track direction
624  double tn = tDir.Dot(nDir);
625  //cout << "tn" << endl;
626  //cout << tn << endl;
627 
628  // derivatives of local residuals versus measurements
629  TMatrixD drdm(3, 3);
630  drdm.UnitMatrix();
631  for (int row = 0; row < 3; row++)
632  for (int col = 0; col < 3; col++)
633  drdm(row, col) -= tDir[row] * nDir[col] / tn;
634 
635  //cout << "drdm" << endl;
636  //drdm.Print();
637 
638  // derivatives of measurements versus global alignment parameters
639  TMatrixD dmdg(3, 6);
640  dmdg.Zero();
641  dmdg(0, 0) = 1.; dmdg(0, 4) = -zPred; dmdg(0, 5) = yPred;
642  dmdg(1, 1) = 1.; dmdg(1, 3) = zPred; dmdg(1, 5) = -xPred;
643  dmdg(2, 2) = 1.; dmdg(2, 3) = -yPred; dmdg(2, 4) = xPred;
644 
645  //cout << "dmdg" << endl;
646  //dmdg.Print();
647 
648  // derivatives of local residuals versus global alignment parameters
649  TMatrixD drldrg(3, 3);
650  drldrg.Zero();
651  drldrg(0, 0) = uDir[0]; drldrg(0, 1) = uDir[1]; drldrg(0, 2) = uDir[2];
652  drldrg(1, 0) = vDir[0]; drldrg(1, 1) = vDir[1]; drldrg(1, 2) = vDir[2];
653 
654  //cout << "drldrg" << endl;
655  //drldrg.Print();
656 
657  //cout << "drdm * dmdg" << endl;
658  //(drdm * dmdg).Print();
659 
660  // derivatives of local residuals versus rigid body parameters
661  TMatrixD drldg(3, 6);
662  drldg = drldrg * (drdm * dmdg);
663 
664  //cout << "drldg" << endl;
665  //drldg.Print();
666 
667  // offset to determine labels for sensor sets or individual layers
668  // 0: PXD, TODO 1: SVD, or individual layers
669  // offset 0 is for alignment of whole setup
670  int offset = 0;
671  //if (sensorId > 16704) offset = 20; // svd, but needs to introduce new 6 constraints: sum rot&shifts of pxd&svd = 0
672 
673  derGlobal(0, 6) = drldg(0, 0); labGlobal.push_back(offset + 1);
674  derGlobal(0, 7) = drldg(0, 1); labGlobal.push_back(offset + 2);
675  derGlobal(0, 8) = drldg(0, 2); labGlobal.push_back(offset + 3);
676  derGlobal(0, 9) = drldg(0, 3); labGlobal.push_back(offset + 4);
677  derGlobal(0, 10) = drldg(0, 4); labGlobal.push_back(offset + 5);
678  derGlobal(0, 11) = drldg(0, 5); labGlobal.push_back(offset + 6);
679 
680  derGlobal(1, 6) = drldg(1, 0);
681  derGlobal(1, 7) = drldg(1, 1);
682  derGlobal(1, 8) = drldg(1, 2);
683  derGlobal(1, 9) = drldg(1, 3);
684  derGlobal(1, 10) = drldg(1, 4);
685  derGlobal(1, 11) = drldg(1, 5);
686 
687 
688 
689  measPoint.addGlobals(labGlobal, derGlobal);
690  listOfPoints.push_back(measPoint);
691  listOfLayers.push_back((unsigned int) sensorId);
692  n_gbl_points++;
693  n_gbl_meas_points++;
694  } else {
695  // Incompatible measurement, store point without measurement
696  GblPoint dummyPoint(jacPointToPoint);
697  listOfPoints.push_back(dummyPoint);
698  listOfLayers.push_back((unsigned int) sensorId);
699  n_gbl_points++;
700  skipMeasurement = false;
701  #ifdef DEBUG
702  cout << " Dummy point inserted. " << endl;
703  #endif
704  }
705 
706 
707  //cout << " Starting extrapolation..." << endl;
708  try {
709 
710  // Extrapolate to next measurement to get material distribution
711  if (ipoint_meas < npoints_meas - 1) {
712  // Check if fitter info is in place
713  if (!trk->getPoint(ipoint_meas + 1)->hasFitterInfo(rep)) {
714  cout << " ERROR: Measurement point does not have a fitter info. Track will be skipped." << endl;
715  return;
716  }
717  // Fitter of next point info which is only used now to get the plane
718  KalmanFitterInfo* fi_i_plus_1 = dynamic_cast<KalmanFitterInfo*>(trk->getPoint(ipoint_meas + 1)->getFitterInfo(rep));
719  if (!fi_i_plus_1) {
720  cout << " ERROR: KalmanFitterInfo (with reference state) for measurement does not exist. Track will be skipped." << endl;
721  return;
722  }
723  StateOnPlane refCopy(*reference);
724  // Extrap to point + 1, do NOT stop at boundary
725  rep->extrapolateToPlane(refCopy, trk->getPoint(ipoint_meas + 1)->getFitterInfo(rep)->getPlane(), false, false);
726  getScattererFromMatList(trackLen,
727  scatTheta,
728  scatSMean,
729  scatDeltaS,
730  trackMomMag,
731  particleMass,
732  particleCharge,
733  rep->getSteps());
734  // Now calculate positions and scattering variance for equivalent scatterers
735  // (Solution from Claus Kleinwort (DESY))
736  s1 = 0.;
737  s2 = scatSMean + scatDeltaS * scatDeltaS / (scatSMean - s1);
738  theta1 = sqrt(scatTheta * scatTheta * scatDeltaS * scatDeltaS / (scatDeltaS * scatDeltaS + (scatSMean - s1) * (scatSMean - s1)));
739  theta2 = sqrt(scatTheta * scatTheta * (scatSMean - s1) * (scatSMean - s1) / (scatDeltaS * scatDeltaS + (scatSMean - s1) * (scatSMean - s1)));
740 
741  if (m_enableScatterers && !m_enableIntermediateScatterer) {
742  theta1 = scatTheta;
743  theta2 = 0;
744  } else if (!m_enableScatterers) {
745  theta1 = 0.;
746  theta2 = 0.;
747  }
748 
749  if (s2 < 1.e-4 || s2 >= trackLen - 1.e-4 || s2 <= 1.e-4) {
750  cout << " WARNING: GBL points will be too close. GBLTrajectory construction might fail. Let's try it..." << endl;
751  }
752 
753  }
754  // Return back to state on current plane
755  delete reference;
756  reference = new ReferenceStateOnPlane(*fi->getReferenceState());
757 
758  // If not last measurement, extrapolate and get jacobians for scattering points between this and next measurement
759  if (ipoint_meas < npoints_meas - 1) {
760  if (theta2 > scatEpsilon) {
761  // First scatterer will be placed at current measurement point (see bellow)
762 
763  // theta2 > 0 ... we want second scatterer:
764  // Extrapolate to s2 (remember we have s1 = 0)
765  rep->extrapolateBy(*reference, s2, false, true);
766  rep->getForwardJacobianAndNoise(jacMeas2Scat, noise, deltaState);
767  // Finish extrapolation to next measurement
768  double nextStep = rep->extrapolateToPlane(*reference, trk->getPoint(ipoint_meas + 1)->getFitterInfo(rep)->getPlane(), false, true);
769  if (0. > nextStep) {
770  cout << " ERROR: The extrapolation to measurement point " << (ipoint_meas + 2) << " stepped back by " << nextStep << "cm !!! Track will be skipped." << endl;
771  return;
772  }
773  rep->getForwardJacobianAndNoise(jacPointToPoint, noise, deltaState);
774 
775  } else {
776  // No scattering: extrapolate whole distance between measurements
777  rep->extrapolateToPlane(*reference, trk->getPoint(ipoint_meas + 1)->getFitterInfo(rep)->getPlane(), false, true);
778  //NOTE: we will load the jacobian from this extrapolation in next loop into measurement point
779  //jacPointToPoint.Print();
780  rep->getForwardJacobianAndNoise(jacPointToPoint, noise, deltaState);
781  //jacPointToPoint.Print();
782  }
783  }
784  } catch (...) {
785  cout << " ERROR: Extrapolation failed. Track will be skipped." << endl;
786  return;
787  }
788  //cout << " Extrapolation finished." << endl;
789 
790 
791  // Now store scatterers if not last measurement and if we decided
792  // there should be scatteres, otherwise the jacobian in measurement
793  // stored above is already correct
794  if (ipoint_meas < npoints_meas - 1) {
795 
796  if (theta1 > scatEpsilon) {
797  // We have to insert first scatterer at measurement point
798  // Therefore (because state is perpendicular to plane, NOT track)
799  // we have non-diaonal matrix of multiple scattering covariance
800  // We have to project scattering into plane coordinates
801  double c1 = trackDir.Dot(plane->getU());
802  double c2 = trackDir.Dot(plane->getV());
803  TMatrixDSym scatCov(2);
804  scatCov(0, 0) = 1. - c2 * c2;
805  scatCov(1, 1) = 1. - c1 * c1;
806  scatCov(0, 1) = c1 * c2;
807  scatCov(1, 0) = c1 * c2;
808  scatCov *= theta1 * theta1 / (1. - c1 * c1 - c2 * c2) / (1. - c1 * c1 - c2 * c2) ;
809 
810  // last point is the just inserted measurement (or dummy point)
811  GblPoint& lastPoint = listOfPoints.at(n_gbl_points - 1);
812  lastPoint.addScatterer(scatResidual, scatCov.Invert());
813 
814  }
815 
816  if (theta2 > scatEpsilon) {
817  // second scatterer is somewhere between measurements
818  // TrackRep state is perpendicular to track direction if using extrapolateBy (I asked Johannes Rauch),
819  // therefore scattering covariance is diagonal (and both elements are equal)
820  TMatrixDSym scatCov(2);
821  scatCov.Zero();
822  scatCov(0, 0) = theta2 * theta2;
823  scatCov(1, 1) = theta2 * theta2;
824 
825  GblPoint scatPoint(jacMeas2Scat);
826  scatPoint.addScatterer(scatResidual, scatCov.Invert());
827  listOfPoints.push_back(scatPoint);
828  listOfLayers.push_back(((unsigned int) sensorId) + 0.5);
829  n_gbl_points++;
830  }
831 
832 
833  }
834  // Free memory on the heap
835  delete reference;
836  }
837  // We should have at least two measurement points to fit anything
838  if (n_gbl_meas_points > 1) {
839  int fitRes = -1;
840  double pvalue = 0.;
841  GblTrajectory* traj = 0;
842  try {
843  // Construct the GBL trajectory, seed not used
844  traj = new GblTrajectory(listOfPoints, seedLabel, clSeed, fitQoverP);
845  // Fit the trajectory
846  fitRes = traj->fit(Chi2, Ndf, lostWeight, m_gblInternalIterations);
847  if (fitRes != 0) {
848  //#ifdef DEBUG
849  //traj->printTrajectory(100);
850  //traj->printData();
851  //traj->printPoints(100);
852  //#endif
853  }
854  } catch (...) {
855  // Gbl failed critically (usually GblPoint::getDerivatives ... singular matrix inversion)
856  return;
857  }
858 
859  pvalue = TMath::Prob(Chi2, Ndf);
860 
861  //traj->printTrajectory(100);
862  //traj->printData();
863  //traj->printPoints(100);
864 
865  #ifdef OUTPUT
866  // Fill histogram with fit result
867  fitResHisto->Fill(fitRes);
868  ndfHisto->Fill(Ndf);
869  #endif
870 
871  #ifdef DEBUG
872  cout << " Ref. Track Chi2 : " << trkChi2 << endl;
873  cout << " Ref. end momentum : " << trackMomMag << " GeV/c ";
874  if (abs(trk->getCardinalRep()->getPDG()) == 11) {
875  if (particleCharge < 0.)
876  cout << "(electron)";
877  else
878  cout << "(positron)";
879  }
880  cout << endl;
881  cout << "------------------ GBL Fit Results --------------------" << endl;
882  cout << " Fit q/p parameter : " << ((fitQoverP) ? ("True") : ("False")) << endl;
883  cout << " Valid trajectory : " << ((traj->isValid()) ? ("True") : ("False")) << endl;
884  cout << " Fit result : " << fitRes << " (0 for success)" << endl;
885  cout << " # GBL meas. points : " << n_gbl_meas_points << endl;
886  cout << " # GBL all points : " << n_gbl_points << endl;
887  cout << " GBL track NDF : " << Ndf << " (-1 for failure)" << endl;
888  cout << " GBL track Chi2 : " << Chi2 << endl;
889  cout << " GBL track P-value : " << pvalue;
890  if (pvalue < m_pValueCut)
891  cout << " < p-value cut " << m_pValueCut;
892  cout << endl;
893  cout << "-------------------------------------------------------" << endl;
894  #endif
895 
896  #ifdef OUTPUT
897  bool hittedLayers[12];
898  for (int hl = 0; hl < 12; hl++) {
899  hittedLayers[hl] = false;
900  }
901  #endif
902 
903  // GBL fit succeded if Ndf >= 0, but Ndf = 0 is useless
904  //TODO: Here should be some track quality check
905  // if (Ndf > 0 && fitRes == 0) {
906  if (traj->isValid() && pvalue >= m_pValueCut && Ndf >= m_minNdf) {
907 
908  // In case someone forgot to use beginRun and dind't reset mille file name to ""
909  if (!milleFile && m_milleFileName != "")
910  milleFile = new MilleBinary(m_milleFileName);
911 
912  // Loop over all GBL points
913  for (unsigned int p = 0; p < listOfPoints.size(); p++) {
914  unsigned int label = p + 1;
915  unsigned int numRes;
916  TVectorD residuals(2);
917  TVectorD measErrors(2);
918  TVectorD resErrors(2);
919  TVectorD downWeights(2);
920  //TODO: now we only provide info about measurements, not kinks
921  if (!listOfPoints.at(p).hasMeasurement())
922  continue;
923 
924  #ifdef OUTPUT
925  // Decode VxdId and get layer in TB setup
926  unsigned int l = 0;
927  unsigned int id = listOfLayers.at(p);
928  // skip segment (5 bits)
929  id = id >> 5;
930  unsigned int sensor = id & 7;
931  id = id >> 3;
932  unsigned int ladder = id & 31;
933  id = id >> 5;
934  unsigned int layer = id & 7;
935 
936  if (layer == 7 && ladder == 2) {
937  l = sensor;
938  } else if (layer == 7 && ladder == 3) {
939  l = sensor + 9 - 3;
940  } else {
941  l = layer + 3;
942  }
943 
944  hittedLayers[l - 1] = true;
945  #endif
946  TVectorD localPar(5);
947  TMatrixDSym localCov(5);
948  traj->getResults(label, localPar, localCov);
949  // Get GBL fit results
950  traj->getMeasResults(label, numRes, residuals, measErrors, resErrors, downWeights);
951  if (m_chi2Cut && (fabs(residuals[0] / resErrors[0]) > m_chi2Cut || fabs(residuals[1] / resErrors[1]) > m_chi2Cut))
952  return;
953  // Write layer-wise data
954  #ifdef OUTPUT
955  if (!writeHistoDataForLabel(listOfLayers.at(p), residuals, measErrors, resErrors, downWeights, localPar, localCov))
956  return;
957  #endif
958 
959  } // end for points
960 
961  // Write binary data to mille binary
962  if (milleFile && m_milleFileName != "" && pvalue >= m_pValueCut && Ndf >= m_minNdf) {
963  traj->milleOut(*milleFile);
964  #ifdef DEBUG
965  cout << " GBL Track written to Millepede II binary file." << endl;
966  cout << "-------------------------------------------------------" << endl;
967  #endif
968  }
969 
970  #ifdef OUTPUT
971  // Fill histograms
972  chi2OndfHisto->Fill(Chi2 / Ndf);
973  pValueHisto->Fill(TMath::Prob(Chi2, Ndf));
974  // track counting
975  stats->Fill(0);
976  // hitted sensors statistics
977  if (
978  hittedLayers[0] &&
979  hittedLayers[1] &&
980  hittedLayers[2] &&
981  hittedLayers[3] &&
982  hittedLayers[4] &&
983  hittedLayers[5] &&
984  hittedLayers[6] &&
985  hittedLayers[7] &&
986  hittedLayers[8]
987  ) {
988  // front tel + pxd + svd
989  stats->Fill(1);
990  }
991 
992  if (
993  hittedLayers[0] &&
994  hittedLayers[1] &&
995  hittedLayers[2] &&
996  hittedLayers[3] &&
997  hittedLayers[4] &&
998  hittedLayers[5] &&
999  hittedLayers[6] &&
1000  hittedLayers[7] &&
1001  hittedLayers[8] &&
1002  hittedLayers[9] &&
1003  hittedLayers[10] &&
1004  hittedLayers[11]
1005  ) {
1006  // all layers
1007  stats->Fill(2);
1008  }
1009  if (
1010  hittedLayers[3] &&
1011  hittedLayers[4] &&
1012  hittedLayers[5] &&
1013  hittedLayers[6] &&
1014  hittedLayers[7] &&
1015  hittedLayers[8]
1016  ) {
1017  // vxd
1018  stats->Fill(3);
1019  }
1020  if (
1021  hittedLayers[5] &&
1022  hittedLayers[6] &&
1023  hittedLayers[7] &&
1024  hittedLayers[8]
1025  ) {
1026  // svd
1027  stats->Fill(4);
1028  }
1029  if (
1030  hittedLayers[0] &&
1031  hittedLayers[1] &&
1032  hittedLayers[2] &&
1033 
1034  hittedLayers[5] &&
1035  hittedLayers[6] &&
1036  hittedLayers[7] &&
1037  hittedLayers[8] &&
1038  hittedLayers[9] &&
1039  hittedLayers[10] &&
1040  hittedLayers[11]
1041  ) {
1042  // all except pxd
1043  stats->Fill(5);
1044  }
1045  if (
1046  hittedLayers[0] &&
1047  hittedLayers[1] &&
1048  hittedLayers[2] &&
1049 
1050  hittedLayers[5] &&
1051  hittedLayers[6] &&
1052  hittedLayers[7] &&
1053  hittedLayers[8]
1054  ) {
1055  // front tel + svd
1056  stats->Fill(6);
1057  }
1058  if (
1059  hittedLayers[9] &&
1060  hittedLayers[10] &&
1061  hittedLayers[11]
1062  ) {
1063  // backward tel
1064  stats->Fill(7);
1065  }
1066  #endif
1067  }
1068 
1069  // Free memory
1070  delete traj;
1071  }
1072 }
1073 
gbl::GblTrajectory::getMeasResults
unsigned int getMeasResults(unsigned int aLabel, unsigned int &numRes, TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors, TVectorD &aDownWeights)
Get residuals from fit at point for measurement.
Definition: GblTrajectory.cc:661
gbl::GblTrajectory::getResults
unsigned int getResults(int aSignedLabel, TVectorD &localPar, TMatrixDSym &localCov) const
Get fit results at point.
Definition: GblTrajectory.cc:631
gbl::GblTrajectory::milleOut
void milleOut(MilleBinary &aMille)
Write valid trajectory to Millepede-II binary file.
Definition: GblTrajectory.cc:1091
genfit::TrackPoint
Object containing AbsMeasurement and AbsFitterInfo objects.
Definition: TrackPoint.h:46
genfit::SharedPlanePtr
std::shared_ptr< genfit::DetPlane > SharedPlanePtr
Shared Pointer to a DetPlane.
Definition: SharedPlanePtr.h:40
genfit::AbsTrackRep::extrapolateToPlane
virtual double extrapolateToPlane(StateOnPlane &state, const genfit::SharedPlanePtr &plane, bool stopAtBoundary=false, bool calcJacobianNoise=false) const =0
Extrapolates the state to plane, and returns the extrapolation length and, via reference,...
genfit::AbsTrackRep::getDir
TVector3 getDir(const StateOnPlane &state) const
Get the direction vector of a state.
Definition: AbsTrackRep.h:246
genfit::AbsTrackRep::getPDG
int getPDG() const
Get the pdg code.
Definition: AbsTrackRep.h:272
genfit::ReferenceStateOnPlane
#StateOnPlane with linearized transport to that #ReferenceStateOnPlane from previous and next #Refere...
Definition: ReferenceStateOnPlane.h:43
genfit::AbsTrackRep::getMomMag
virtual double getMomMag(const StateOnPlane &state) const =0
get the magnitude of the momentum in GeV.
genfit::AbsHMatrix::getMatrix
virtual const TMatrixD & getMatrix() const =0
Get the actual matrix representation.
genfit::StateOnPlane
A state with arbitrary dimension defined in a DetPlane.
Definition: StateOnPlane.h:47
genfit
Defines for I/O streams used for error and debug printing.
Definition: AlignablePXDRecoHit.h:19
genfit::GFGbl::processTrackWithRep
void processTrackWithRep(Track *trk, const AbsTrackRep *rep, bool resortHits=false) override
Performs fit on a Track.
Definition: GFGbl.cc:326
genfit::Track
Collection of TrackPoint objects, AbsTrackRep objects and FitStatus objects.
Definition: Track.h:71
gbl::GblPoint::addScatterer
void addScatterer(const TVectorD &aResiduals, const TVectorD &aPrecision)
Add a (thin) scatterer to a point.
Definition: GblPoint.cc:210
genfit::FieldManager::getFieldVal
TVector3 getFieldVal(const TVector3 &position)
This does NOT use the cache!
Definition: FieldManager.h:63
gbl::GblPoint::addGlobals
void addGlobals(const std::vector< int > &aLabels, const TMatrixD &aDerivatives)
Add global derivatives to a point.
Definition: GblPoint.cc:320
genfit::AbsTrackRep
Abstract base class for a track representation.
Definition: AbsTrackRep.h:66
genfit::AbsTrackRep::extrapolateBy
virtual double extrapolateBy(StateOnPlane &state, double step, bool stopAtBoundary=false, bool calcJacobianNoise=false) const =0
Extrapolates the state by step (cm) and returns the extrapolation length and, via reference,...
genfit::AbsTrackRep::getCharge
virtual double getCharge(const StateOnPlane &state) const =0
Get the (fitted) charge of a state.
gbl::GblTrajectory
GBL trajectory.
Definition: GblTrajectory.h:48
genfit::PlanarMeasurement
Measurement class implementing a planar hit geometry (1 or 2D).
Definition: PlanarMeasurement.h:44
genfit::HMatrixV
AbsHMatrix implementation for one-dimensional MeasurementOnPlane and RKTrackRep parameterization.
Definition: HMatrixV.h:37
gbl
Namespace for the general broken lines package.
Definition: BorderedBandMatrix.h:44
genfit::TrackPoint::getFitterInfo
AbsFitterInfo * getFitterInfo(const AbsTrackRep *rep=nullptr) const
Get fitterInfo for rep. Per default, use cardinal rep.
Definition: TrackPoint.cc:170
genfit::AbsMeasurement
Contains the measurement and covariance in raw detector coordinates.
Definition: AbsMeasurement.h:42
genfit::KalmanFitterInfo
Collects information needed and produced by a AbsKalmanFitter implementations and is specific to one ...
Definition: KalmanFitterInfo.h:44
GblPoint.h
genfit::HMatrixU
AbsHMatrix implementation for one-dimensional MeasurementOnPlane and RKTrackRep parameterization.
Definition: HMatrixU.h:37
genfit::AbsTrackRep::getMass
double getMass(const StateOnPlane &state) const
Get tha particle mass in GeV/c^2.
Definition: AbsTrackRep.cc:100
genfit::AbsTrackRep::getDim
virtual unsigned int getDim() const =0
Get the dimension of the state vector used by the track representation.
genfit::AbsMeasurement::constructHMatrix
virtual const AbsHMatrix * constructHMatrix(const AbsTrackRep *) const =0
Returns a new AbsHMatrix object.
genfit::GFGbl::beginRun
void beginRun()
Creates the mille binary file for output of data for Millepede II alignment, can be set by setMP2Opti...
Definition: GFGbl.cc:184
genfit::AbsTrackRep::getSteps
virtual std::vector< genfit::MatStep > getSteps() const =0
Get stepsizes and material properties of crossed materials of the last extrapolation.
genfit::MatStep
Simple struct containing MaterialProperties and stepsize in the material.
Definition: AbsTrackRep.h:42
genfit::GFGbl::endRun
void endRun()
Required to write and close ROOT file with debug output.
Definition: GFGbl.cc:231
genfit::AbsHMatrix::Hv
virtual TVectorD Hv(const TVectorD &v) const
H*v.
Definition: AbsHMatrix.h:49
gbl::GblTrajectory::isValid
bool isValid() const
Retrieve validity of trajectory.
Definition: GblTrajectory.cc:289
GblTrajectory.h
genfit::FieldManager::getInstance
static FieldManager * getInstance()
Get singleton instance.
Definition: FieldManager.h:119
genfit::Track::getCardinalRep
AbsTrackRep * getCardinalRep() const
Get cardinal track representation.
Definition: Track.h:145
gbl::MilleBinary
Millepede-II (binary) record.
Definition: MilleBinary.h:68
gbl::GblPoint::addMeasurement
void addMeasurement(const TMatrixD &aProjection, const TVectorD &aResiduals, const TVectorD &aPrecision, double minPrecision=0.)
Add a measurement to a point.
Definition: GblPoint.cc:69
gbl::GblPoint
Point on trajectory.
Definition: GblPoint.h:68
genfit::AbsFitter
Abstract base class for fitters.
Definition: AbsFitter.h:35
genfit::AbsTrackRep::getForwardJacobianAndNoise
virtual void getForwardJacobianAndNoise(TMatrixD &jacobian, TMatrixDSym &noise, TVectorD &deltaState) const =0
Get the jacobian and noise matrix of the last extrapolation.
gbl::GblTrajectory::fit
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, std::string optionList="")
Perform fit of (valid) trajectory.
Definition: GblTrajectory.cc:1043