Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GblFitter.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GblFitter.cc
1 //-*-mode: C++; c-basic-offset: 2; -*-
2 /* Copyright 2013-2014
3  * Authors: Sergey Yashchenko and Tadeas Bilka
4  *
5  * This is an interface to General Broken Lines
6  *
7  * Version: 6 --------------------------------------------------------------
8  * - complete rewrite to GblFitter using GblFitterInfo and GblFitStatus
9  * - mathematics should be the same except for additional iterations
10  * (not possible before)
11  * - track is populated with scatterers + additional points with
12  * scatterers and no measurement (optional)
13  * - final track contains GblFitStatus and fitted states from GBL prediction
14  * - 1D/2D hits supported (pixel, single strip, combined strips(2D), wire)
15  * - At point: Only the very first raw measurement is used and from
16  * that, constructed MeasurementOnPlane with heighest weight
17  * Version: 5 --------------------------------------------------------------
18  * - several bug-fixes:
19  * - Scatterers at bad points
20  * - Jacobians at a point before they should be (code reorganized)
21  * - Change of sign of residuals
22  * Version: 4 --------------------------------------------------------------
23  * Fixed calculation of equvivalent scatterers (solution by C. Kleinwort)
24  * Now a scatterer is inserted at each measurement (except last) and
25  * between each two measurements.
26  * TrueHits/Clusters. Ghost (1D) hits ignored. With or
27  * without magnetic field.
28  * Version: 3 --------------------------------------------------------------
29  * This version now supports both TrueHits and Clusters for VXD.
30  * It can be used for arbitrary material distribution between
31  * measurements. Moments of scattering distribution are computed
32  * and translated into two equivalent thin GBL scatterers placed
33  * at computed positions between measurement points.
34  * Version: 2 ... never published -----------------------------------------
35  * Scatterer at each boundary (tooo many scatterers). TrueHits/Clusters.
36  * Without global der.&MP2 output.
37  * Version: 1 --------------------------------------------------------------
38  * Scatterers at measurement planes. TrueHits
39  * Version: 0 --------------------------------------------------------------
40  * Without scatterers. Genfit 1.
41  * -------------------------------------------------------------------------
42  *
43  * This file is part of GENFIT.
44  *
45  * GENFIT is free software: you can redistribute it and/or modify
46  * it under the terms of the GNU Lesser General Public License as published
47  * by the Free Software Foundation, either version 3 of the License, or
48  * (at your option) any later version.
49  *
50  * GENFIT is distributed in the hope that it will be useful,
51  * but WITHOUT ANY WARRANTY; without even the implied warranty of
52  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
53  * GNU Lesser General Public License for more details.
54  *
55  * You should have received a copy of the GNU Lesser General Public License
56  * along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
57  */
58 
59 #include "GblFitter.h"
60 #include "../include/GblFitStatus.h"
61 #include "GblFitterInfo.h"
62 #include "GblTrajectory.h"
63 #include "GblPoint.h"
65 
66 #include "Track.h"
67 #include <TFile.h>
68 #include <TH1F.h>
69 #include <TTree.h>
70 #include <string>
71 #include <list>
72 #include <FieldManager.h>
73 #include <HMatrixU.h>
74 #include <HMatrixV.h>
75 #include <Math/SMatrix.h>
76 #include <TMatrixD.h>
77 #include <TVectorDfwd.h>
78 #include <TMatrixT.h>
79 #include <TVector3.h>
80 
81 //#define DEBUG
82 
83 using namespace gbl;
84 using namespace std;
85 using namespace genfit;
86 
90 GblFitter::~GblFitter() {
91  if (m_segmentController) {
92  delete m_segmentController;
93  m_segmentController = nullptr;
94  }
95 }
96 
97 void GblFitter::setTrackSegmentController(GblTrackSegmentController* controler)
98 {
99  if (m_segmentController) {
100  delete m_segmentController;
101  m_segmentController = nullptr;
102  }
103  m_segmentController = controler;
104 }
105 
106 void GblFitter::processTrackWithRep(Track* trk, const AbsTrackRep* rep, bool resortHits)
107 {
108  cleanGblInfo(trk, rep);
109 
110  if (resortHits)
111  sortHits(trk, rep);
112 
113  // This flag enables/disables fitting of q/p parameter in GBL
114  // It is switched off automatically if no B-field at (0,0,0) is detected.
115  bool fitQoverP = true;
116  //TODO: Use clever way to determine zero B-field
117  double Bfield = genfit::FieldManager::getInstance()->getFieldVal(TVector3(0., 0., 0.)).Mag();
118  if (!(Bfield > 1.e-16))
119  fitQoverP = false;
120  // degrees of freedom after fit
121  int Ndf = 0;
122  // Chi2 after fit
123  double Chi2 = 0.;
124  //FIXME: d-w's not used so far...
125  double lostWeight = 0.;
126 
127  // Preparation of points (+add reference states) for GBL fit
128  // -----------------------------------------------------------------
130  trk->setFitStatus(gblfs, rep);
131  gblfs->setCurvature(fitQoverP);
132  //
133  // Propagate reference seed, create scattering points, calc Jacobians
134  // and store everything in fitter infos. (ready to collect points and fit)
135  //
136  //
137  gblfs->setTrackLen(
138  //
139  constructGblInfo(trk, rep)
140  //
141  );
142  //
143  //
144  gblfs->setIsFittedWithReferenceTrack(true);
145  gblfs->setNumIterations(0); //default value, still valid, No GBL iteration
146  if (m_externalIterations < 1)
147  return;
148  // -----------------------------------------------------------------
149 
150  // cppcheck-suppress unreadVariable
151  unsigned int nFailed = 0;
152  // cppcheck-suppress unreadVariable
153  int fitRes = 0;
154  std::vector<std::string> gblIterations;
155  gblIterations.push_back(m_gblInternalIterations);
156 
157  // Iterations and updates of fitter infos and fit status
158  // -------------------------------------------------------------------
159  for (unsigned int iIter = 0; iIter < m_externalIterations; iIter++) {
160  // GBL refit (1st of reference, then refit of GBL trajectory itself)
161  int nscat = 0, nmeas = 0, ndummy = 0;
162  std::vector<gbl::GblPoint> points = collectGblPoints(trk, rep);
163  for(unsigned int ip = 0;ip<points.size(); ip++) {
164  GblPoint & p = points.at(ip);
165  if (p.hasScatterer())
166  nscat++;
167  if (p.hasMeasurement())
168  nmeas++;
169  if(!p.hasMeasurement()&&!p.hasScatterer())
170  ndummy++;
171  }
172  gbl::GblTrajectory traj(points, gblfs->hasCurvature());
173 
174  fitRes = traj.fit(Chi2, Ndf, lostWeight, (iIter == m_externalIterations - 1) ? m_gblInternalIterations : "");
175 
176  // Update fit results in fitterinfos
177  updateGblInfo(traj, trk, rep);
178 
179  // This repropagates to get new Jacobians,
180  // if planes changed, predictions are extrapolated to new planes
181  if (m_recalcJacobians > iIter) {
182  GblFitterInfo* prevFitterInfo = 0;
183  GblFitterInfo* currFitterInfo = 0;
184  for (unsigned int ip = 0; ip < trk->getNumPoints(); ip++) {
185  if (trk->getPoint(ip)->hasFitterInfo(rep) && (currFitterInfo = dynamic_cast<GblFitterInfo*>(trk->getPoint(ip)->getFitterInfo(rep)))) {
186 
187  currFitterInfo->recalculateJacobian(prevFitterInfo);
188  prevFitterInfo = currFitterInfo;
189  }
190  }
191  }
192 
193  gblfs->setIsFitted(true);
194  gblfs->setIsFitConvergedPartially(fitRes == 0);
195  nFailed = trk->getNumPointsWithMeasurement() - nmeas;
196  gblfs->setNFailedPoints(nFailed);
197  gblfs->setIsFitConvergedFully(fitRes == 0 && nFailed == 0);
198  gblfs->setNumIterations(iIter + 1);
199  gblfs->setChi2(Chi2);
200  gblfs->setNdf(Ndf);
201  gblfs->setCharge(trk->getFittedState().getCharge());
202 
203  #ifdef DEBUG
204  int npoints_meas = trk->getNumPointsWithMeasurement();
205  int npoints_all = trk->getNumPoints();
206 
207  cout << "-------------------------------------------------------" << endl;
208  cout << " GBL processed genfit::Track " << endl;
209  cout << "-------------------------------------------------------" << endl;
210  cout << " # Track Points : " << npoints_all << endl;
211  cout << " # Meas. Points : " << npoints_meas << endl;
212  cout << " # GBL points all : " << traj.getNumPoints();
213  if (ndummy)
214  cout << " (" << ndummy << " dummy) ";
215  cout << endl;
216  cout << " # GBL points meas : " << nmeas << endl;
217  cout << " # GBL points scat : " << nscat << endl;
218  cout << "-------------- GBL Fit Results ----------- Iteration " << iIter+1 << " " << ((iIter < gblIterations.size()) ? gblIterations[iIter] : "") << endl;
219  cout << " Fit q/p parameter : " << (gblfs->hasCurvature() ? ("True") : ("False")) << endl;
220  cout << " Valid trajectory : " << ((traj.isValid()) ? ("True") : ("False")) << endl;
221  cout << " Fit result : " << fitRes << " (0 for success)" << endl;
222  cout << " GBL track NDF : " << Ndf << " (-1 for failure)" << endl;
223  cout << " GBL track Chi2 : " << Chi2 << endl;
224  cout << " GBL track P-value : " << TMath::Prob(Chi2, Ndf) << endl;
225  cout << "-------------------------------------------------------" << endl;
226  #endif
227 
228  }
229  // -------------------------------------------------------------------
230 
231 }
232 
233 void GblFitter::cleanGblInfo(Track* trk, const AbsTrackRep* rep) const {
234 
235  for (int ip = trk->getNumPoints() - 1; ip >=0; ip--) {
236  trk->getPoint(ip)->setScatterer(nullptr);
237  trk->getPoint(ip)->deleteFitterInfo(rep);
238  //TODO
239  if (!trk->getPoint(ip)->hasRawMeasurements())
240  trk->deletePoint(ip);
241  }
242 }
243 
244 void GblFitter::sortHits(Track* trk, const AbsTrackRep* rep) const {
245  // All measurement points in ref. track
246  int npoints_meas = trk->getNumPointsWithMeasurement();
247  // Prepare state for extrapolation of track seed
248  StateOnPlane reference(rep);
249  rep->setTime(reference, trk->getTimeSeed());
250  rep->setPosMom(reference, trk->getStateSeed());
251  // Take the state to first plane
252  SharedPlanePtr firstPlane(trk->getPointWithMeasurement(0)->getRawMeasurement(0)->constructPlane(reference));
253  reference.extrapolateToPlane(firstPlane);
254  //1st point is at arc-len=0
255  double arcLenPos = 0;
256 
257  // Loop only between meas. points
258  for (int ipoint_meas = 0; ipoint_meas < npoints_meas - 1; ipoint_meas++) {
259  // current measurement point
260  TrackPoint* point_meas = trk->getPointWithMeasurement(ipoint_meas);
261  // Current detector plane
262  SharedPlanePtr plane = point_meas->getRawMeasurement(0)->constructPlane(reference);
263  // Get the next plane
264  SharedPlanePtr nextPlane(trk->getPointWithMeasurement(ipoint_meas + 1)->getRawMeasurement(0)->constructPlane(reference));
265 
266  point_meas->setSortingParameter(arcLenPos);
267  arcLenPos += reference.extrapolateToPlane(nextPlane);
268 
269  } // end of loop over track points with measurement
270  trk->getPointWithMeasurement(npoints_meas - 1)->setSortingParameter(arcLenPos);
271  trk->sort();
272 }
273 
274 std::vector<gbl::GblPoint> GblFitter::collectGblPoints(genfit::Track* trk, const genfit::AbsTrackRep* rep) {
275  //TODO store collected points in in fit status? need streamer for GblPoint (or something like that)
276  std::vector<gbl::GblPoint> thePoints;
277  thePoints.clear();
278 
279  // Collect points from track and fitterInfo(rep)
280  for (unsigned int ip = 0; ip < trk->getNumPoints(); ip++) {
281  GblFitterInfo * gblfi = dynamic_cast<GblFitterInfo*>(trk->getPoint(ip)->getFitterInfo(rep));
282  if (!gblfi)
283  continue;
284  thePoints.push_back(gblfi->constructGblPoint());
285  }
286  return thePoints;
287 }
288 
289 void GblFitter::updateGblInfo(gbl::GblTrajectory& traj, genfit::Track* trk, const genfit::AbsTrackRep* rep) {
290  //FIXME
291  if (!traj.isValid())
292  return;
293 
294  // Update points in track and fitterInfo(rep)
295  int igblfi = -1;
296  for (unsigned int ip = 0; ip < trk->getNumPoints(); ip++) {
297  GblFitterInfo * gblfi = dynamic_cast<GblFitterInfo*>(trk->getPoint(ip)->getFitterInfo(rep));
298  if (!gblfi)
299  continue;
300  igblfi++;
301 
302  // The point will calculate its position on the track
303  // (counting fitter infos) which hopefully
304  gblfi->updateFitResults(traj);
305 
306  // This is agains logic. User can do this if he wants and it is recommended usually
307  // so that following fit could reuse the updated seed
308  //if (igblfi == 0) {
309  // trk->setStateSeed( gblfi->getFittedState(true).getPos(), gblfi->getFittedState(true).getMom() );
310  // trk->setCovSeed( gblfi->getFittedState(true).get6DCov() );
311  //}
312  }
313 }
314 
316  double& theta, double& s, double& ds,
317  const double p, const double mass, const double charge,
318  const std::vector<genfit::MatStep>& steps) const {
319  theta = 0.; s = 0.; ds = 0.; length = 0;
320  if (steps.empty()) return;
321 
322  // sum of step lengths
323  double len = 0.;
324  // normalization
325  double sumxx = 0.;
326  // first moment (non-normalized)
327  double sumx2x2 = 0.;
328  // (part of) second moment / variance (non-normalized)
329  double sumx3x3 = 0.;
330 
331  // cppcheck-suppress unreadVariable
332  double xmin = 0.;
333  double xmax = 0.;
334 
335  for (unsigned int i = 0; i < steps.size(); i++) {
336  const MatStep step = steps.at(i);
337  // inverse of material radiation length ... (in 1/cm) ... "density of scattering"
338  double rho = 1. / step.material_.radiationLength;
339  len += fabs(step.stepSize_);
340  xmin = xmax;
341  xmax = xmin + fabs(step.stepSize_);
342  // Compute integrals
343 
344  // integral of rho(x)
345  sumxx += rho * (xmax - xmin);
346  // integral of x*rho(x)
347  sumx2x2 += rho * (xmax * xmax - xmin * xmin) / 2.;
348  // integral of x*x*rho(x)
349  sumx3x3 += rho * (xmax * xmax * xmax - xmin * xmin * xmin) / 3.;
350  }
351  // This ensures PDG formula still gives positive results (but sumxx should be > 1e-4 for it to hold)
352  if (sumxx < 1.0e-10) return;
353  // Calculate theta from total sum of radiation length
354  // instead of summimg squares of individual deflection angle variances
355  // PDG formula:
356  double beta = p / sqrt(p * p + mass * mass);
357  theta = (0.0136 / p / beta) * fabs(charge) * sqrt(sumxx) * (1. + 0.038 * log(sumxx));
358  //theta = (0.015 / p / beta) * fabs(charge) * sqrt(sumxx);
359 
360  // track length
361  length = len;
362  // Normalization factor
363  double N = 1. / sumxx;
364  // First moment
365  s = N * sumx2x2;
366  // Square of second moment (variance)
367  // integral of (x - s)*(x - s)*rho(x)
368  double ds_2 = N * (sumx3x3 - 2. * sumx2x2 * s + sumxx * s * s);
369  ds = sqrt(ds_2);
370 
371  #ifdef DEBUG
372 
373  #endif
374 }
375 
376 double GblFitter::constructGblInfo(Track* trk, const AbsTrackRep* rep)
377 {
378  // All measurement points in ref. track
379  int npoints_meas = trk->getNumPointsWithMeasurement();
380  // Dimesion of representation/state
381  int dim = rep->getDim();
382  // Jacobian for point with measurement = how to propagate from previous point (scat/meas)
383  TMatrixD jacPointToPoint(dim, dim);
384  jacPointToPoint.UnitMatrix();
385 
386  // Prepare state for extrapolation of track seed
387  // Take the state to first plane
388  StateOnPlane reference(rep);
389  rep->setTime(reference, trk->getTimeSeed());
390  rep->setPosMom(reference, trk->getStateSeed());
391 
392  SharedPlanePtr firstPlane(trk->getPointWithMeasurement(0)->getRawMeasurement(0)->constructPlane(reference));
393  reference.extrapolateToPlane(firstPlane);
394 
395  double sumTrackLen = 0;
396  // NOT used but useful
397  TMatrixDSym noise; TVectorD deltaState;
398 
399  // Loop only between meas. points
400  for (int ipoint_meas = 0; ipoint_meas < npoints_meas; ipoint_meas++) {
401  // current measurement point
402  TrackPoint* point_meas = trk->getPointWithMeasurement(ipoint_meas);
403  // Current detector plane
404  SharedPlanePtr plane = point_meas->getRawMeasurement(0)->constructPlane(reference);
405  // track direction at plane (in global coords)
406  TVector3 trackDir = rep->getDir(reference);
407  // track momentum direction vector at plane (in global coords)
408  double trackMomMag = rep->getMomMag(reference);
409  // charge of particle
410  double particleCharge = rep->getCharge(reference);
411  // mass of particle
412  double particleMass = rep->getMass(reference);
413  // Parameters of a thick scatterer between measurements
414  double trackLen = 0., scatTheta = 0., scatSMean = 0., scatDeltaS = 0.;
415  // Parameters of two equivalent thin scatterers
416  double theta1 = 0., theta2 = 0., s1 = 0., s2 = 0.;
417  // jacobian from s1=0 to s2
418  TMatrixD jacMeas2Scat(dim, dim);
419  jacMeas2Scat.UnitMatrix();
420 
421  // Stop here if we are at last point (do not add scatterers to last point),
422  // just the fitter info
423  if (ipoint_meas >= npoints_meas - 1) {
424 
425  // Construction last measurement (no scatterer)
426  // --------------------------------------------
427  // Just add the fitter info of last plane
428  GblFitterInfo* gblfimeas(new GblFitterInfo(point_meas, rep, reference));
429  gblfimeas->setJacobian(jacPointToPoint);
430  point_meas->setFitterInfo(gblfimeas);
431  // --------------------------------------------
432 
433  break;
434  }
435  // Extrapolate to next measurement to get material distribution
436  // Use a temp copy of the StateOnPlane to propage between measurements
437  StateOnPlane refCopy(reference);
438  // Get the next plane
439  SharedPlanePtr nextPlane(trk->getPointWithMeasurement(ipoint_meas + 1)->getRawMeasurement(0)->constructPlane(reference));
440 
441  // Extrapolation for multiple scattering calculation
442  // Extrap to point + 1, do NOT stop at boundary
443  TVector3 segmentEntry = refCopy.getPos();
444  rep->extrapolateToPlane(refCopy, nextPlane, false, false);
445  TVector3 segmentExit = refCopy.getPos();
446 
447  getScattererFromMatList(trackLen,
448  scatTheta,
449  scatSMean,
450  scatDeltaS,
451  trackMomMag,
452  particleMass,
453  particleCharge,
454  rep->getSteps());
455  // Now calculate positions and scattering variance for equivalent scatterers
456  // (Solution from Claus Kleinwort (DESY))
457  s1 = 0.; s2 = scatSMean + scatDeltaS * scatDeltaS / (scatSMean - s1);
458  theta1 = sqrt(scatTheta * scatTheta * scatDeltaS * scatDeltaS / (scatDeltaS * scatDeltaS + (scatSMean - s1) * (scatSMean - s1)));
459  theta2 = sqrt(scatTheta * scatTheta * (scatSMean - s1) * (scatSMean - s1) / (scatDeltaS * scatDeltaS + (scatSMean - s1) * (scatSMean - s1)));
460 
461  // Call segment controller to set MS options:
462  if (m_segmentController)
463  m_segmentController->controlTrackSegment(segmentEntry, segmentExit, scatTheta, this);
464 
465  // Scattering options: OFF / THIN / THICK
466  if (m_enableScatterers && !m_enableIntermediateScatterer) {
467  theta1 = scatTheta;
468  theta2 = 0;
469  } else if (!m_enableScatterers) {
470  theta1 = 0.;
471  theta2 = 0.;
472  }
473 
474  // Construction of measurement (with scatterer)
475  // --------------------------------------------
476 
477  if (theta1 > scatEpsilon)
478  point_meas->setScatterer(new ThinScatterer(plane, Material(theta1, 0., 0., 0., 0.)));
479 
480  GblFitterInfo* gblfimeas(new GblFitterInfo(point_meas, rep, reference));
481  gblfimeas->setJacobian(jacPointToPoint);
482  point_meas->setFitterInfo(gblfimeas);
483  // --------------------------------------------
484 
485 
486  // If not last measurement, extrapolate and get jacobians for scattering points between this and next measurement
487  if (theta2 > scatEpsilon) {
488  // First scatterer has been placed at current measurement point (see above)
489  // theta2 > 0 ... we want second scatterer:
490  // Extrapolate to s2 (we have s1 = 0)
491  rep->extrapolateBy(reference, s2, false, true);
492  rep->getForwardJacobianAndNoise(jacMeas2Scat, noise, deltaState);
493 
494  // Construction of intermediate scatterer
495  // --------------------------------------
496  TrackPoint* scattp = new TrackPoint(trk);
497  scattp->setSortingParameter(point_meas->getSortingParameter() + s2);
498  scattp->setScatterer(new ThinScatterer(reference.getPlane(), Material(theta2, 0., 0., 0., 0.)));
499  // Add point to track before next point
500  int pointIndex = 0;
501  //TODO Deduce this rather than looping over all points
502  for (unsigned int itp = 0; itp < trk->getNumPoints(); itp++) {
503  if (trk->getPoint(itp) == point_meas) {
504  pointIndex = itp;
505  break;
506  }
507  }
508  trk->insertPoint(scattp, pointIndex + 1);
509  // Create and store fitter info
510  GblFitterInfo * gblfiscat(new GblFitterInfo(scattp, rep, reference));
511  gblfiscat->setJacobian(jacMeas2Scat);
512  scattp->setFitterInfo(gblfiscat);
513  // ---------------------------------------
514 
515  // Finish extrapolation to next measurement
516  double nextStep = rep->extrapolateToPlane(reference, nextPlane, false, true);
517  rep->getForwardJacobianAndNoise(jacPointToPoint, noise, deltaState);
518 
519  if (0. > nextStep) {
520  cout << " ERROR: The extrapolation to measurement point " << (ipoint_meas + 2) << " stepped back by " << nextStep << "cm !!! Track will be cut before this point." << endl;
521  // stop trajectory construction here
522  break;
523  }
524 
525  } else {
526  // No scattering: extrapolate whole distance between measurements
527  double nextStep = rep->extrapolateToPlane(reference, nextPlane, false, true);
528  rep->getForwardJacobianAndNoise(jacPointToPoint, noise, deltaState);
529 
530  if (0. > nextStep) {
531  cout << " ERROR: The extrapolation to measurement point " << (ipoint_meas + 2) << " stepped back by " << nextStep << "cm !!! Track will be cut before this point." << endl;
532  // stop trajectory construction here
533  break;
534  }
535  }
536  // Track length up to next point
537  sumTrackLen += trackLen;
538 
539  } // end of loop over track points with measurement
540  return sumTrackLen;
541 }