GENFIT
Rev:NoNumberAvailable
|
Go to the documentation of this file.
60 #include "../include/GblFitStatus.h"
75 #include <Math/SMatrix.h>
77 #include <TVectorDfwd.h>
91 GblFitter::~GblFitter() {
92 if (m_segmentController) {
93 delete m_segmentController;
94 m_segmentController = NULL;
100 if (m_segmentController) {
101 delete m_segmentController;
102 m_segmentController = NULL;
104 m_segmentController = controler;
109 cleanGblInfo(trk, rep);
116 bool fitQoverP =
true;
119 if (!(Bfield > 1.e-16))
126 double lostWeight = 0.;
140 constructGblInfo(trk, rep)
147 if (m_externalIterations < 1)
152 unsigned int nFailed = 0;
154 std::vector<std::string> gblIterations;
155 gblIterations.push_back(m_gblInternalIterations);
160 for (
unsigned int iIter = 0; iIter < m_externalIterations; iIter++) {
162 int nscat = 0, nmeas = 0, ndummy = 0;
163 std::vector<gbl::GblPoint> points = collectGblPoints(trk, rep);
164 for(
unsigned int ip = 0;ip<points.size(); ip++) {
175 fitRes = traj.
fit(Chi2, Ndf, lostWeight, (iIter == m_externalIterations - 1) ? m_gblInternalIterations :
"");
178 updateGblInfo(traj, trk, rep);
182 if (m_recalcJacobians > iIter) {
185 for (
unsigned int ip = 0; ip < trk->
getNumPoints(); ip++) {
189 prevFitterInfo = currFitterInfo;
208 cout <<
"-------------------------------------------------------" << endl;
209 cout <<
" GBL processed genfit::Track " << endl;
210 cout <<
"-------------------------------------------------------" << endl;
211 cout <<
" # Track Points : " << npoints_all << endl;
212 cout <<
" # Meas. Points : " << npoints_meas << endl;
213 cout <<
" # GBL points all : " << traj.getNumPoints();
215 cout <<
" (" << ndummy <<
" dummy) ";
217 cout <<
" # GBL points meas : " << nmeas << endl;
218 cout <<
" # GBL points scat : " << nscat << endl;
219 cout <<
"-------------- GBL Fit Results ----------- Iteration " << iIter+1 <<
" " << ((iIter < gblIterations.size()) ? gblIterations[iIter] :
"") << endl;
220 cout <<
" Fit q/p parameter : " << (gblfs->
hasCurvature() ? (
"True") : (
"False")) << endl;
221 cout <<
" Valid trajectory : " << ((traj.isValid()) ? (
"True") : (
"False")) << endl;
222 cout <<
" Fit result : " << fitRes <<
" (0 for success)" << endl;
223 cout <<
" GBL track NDF : " << Ndf <<
" (-1 for failure)" << endl;
224 cout <<
" GBL track Chi2 : " << Chi2 << endl;
225 cout <<
" GBL track P-value : " << TMath::Prob(Chi2, Ndf) << endl;
226 cout <<
"-------------------------------------------------------" << endl;
256 double arcLenPos = 0;
259 for (
int ipoint_meas = 0; ipoint_meas < npoints_meas - 1; ipoint_meas++) {
277 std::vector<gbl::GblPoint> thePoints;
281 for (
unsigned int ip = 0; ip < trk->
getNumPoints(); ip++) {
297 for (
unsigned int ip = 0; ip < trk->
getNumPoints(); ip++) {
317 double& theta,
double& s,
double& ds,
318 const double p,
const double mass,
const double charge,
319 const std::vector<genfit::MatStep>& steps)
const {
320 theta = 0.; s = 0.; ds = 0.; length = 0;
321 if (steps.empty())
return;
335 for (
unsigned int i = 0; i < steps.size(); i++) {
336 const MatStep step = steps.at(i);
345 sumxx += rho * (xmax - xmin);
347 sumx2x2 += rho * (xmax * xmax - xmin * xmin) / 2.;
349 sumx3x3 += rho * (xmax * xmax * xmax - xmin * xmin * xmin) / 3.;
352 if (sumxx < 1.0e-10)
return;
356 double beta = p / sqrt(p * p + mass * mass);
357 theta = (0.0136 / p / beta) * fabs(charge) * sqrt(sumxx) * (1. + 0.038 * log(sumxx));
363 double N = 1. / sumxx;
368 double ds_2 = N * (sumx3x3 - 2. * sumx2x2 * s + sumxx * s * s);
383 TMatrixD jacPointToPoint(dim, dim);
384 jacPointToPoint.UnitMatrix();
395 double sumTrackLen = 0;
397 TMatrixDSym noise; TVectorD deltaState;
400 for (
int ipoint_meas = 0; ipoint_meas < npoints_meas; ipoint_meas++) {
406 TVector3 trackDir = rep->
getDir(reference);
408 double trackMomMag = rep->
getMomMag(reference);
410 double particleCharge = rep->
getCharge(reference);
412 double particleMass = rep->
getMass(reference);
414 double trackLen = 0., scatTheta = 0., scatSMean = 0., scatDeltaS = 0.;
416 double theta1 = 0., theta2 = 0., s1 = 0., s2 = 0.;
418 TMatrixD jacMeas2Scat(dim, dim);
419 jacMeas2Scat.UnitMatrix();
423 if (ipoint_meas >= npoints_meas - 1) {
443 TVector3 segmentEntry = refCopy.
getPos();
445 TVector3 segmentExit = refCopy.
getPos();
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)));
462 if (m_segmentController)
463 m_segmentController->controlTrackSegment(segmentEntry, segmentExit,
this);
466 if (m_enableScatterers && !m_enableIntermediateScatterer) {
469 }
else if (!m_enableScatterers) {
502 for (
unsigned int itp = 0; itp < trk->
getNumPoints(); itp++) {
503 if (trk->
getPoint(itp) == point_meas) {
520 cout <<
" ERROR: The extrapolation to measurement point " << (ipoint_meas + 2) <<
" stepped back by " << nextStep <<
"cm !!! Track will be cut before this point." << endl;
531 cout <<
" ERROR: The extrapolation to measurement point " << (ipoint_meas + 2) <<
" stepped back by " << nextStep <<
"cm !!! Track will be cut before this point." << endl;
537 sumTrackLen += trackLen;
virtual double getMomMag(const StateOnPlane &state) const =0
get the magnitude of the momentum in GeV.
Collects information needed and produced by a GblFitter/GBL and is specific to one AbsTrackRep of the...
virtual std::vector< genfit::MatStep > getSteps() const =0
Get stepsizes and material properties of crossed materials of the last extrapolation.
boost::shared_ptr< genfit::DetPlane > SharedPlanePtr
Shared Pointer to a DetPlane.
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,...
Collection of TrackPoint objects, AbsTrackRep objects and FitStatus objects.
void recalculateJacobian(GblFitterInfo *prevFitterInfo)
Re-extrapolates between prevFitterInfo and this point using forward state to update the Jacobian (if ...
void getScattererFromMatList(double &length, double &theta, double &s, double &ds, const double p, const double mass, const double charge, const std::vector< MatStep > &steps)
Evaluates moments of radiation length distribution from list of material steps and computes parameter...
TrackSegmentController for use with GblFitter.
void deleteFitterInfo(const AbsTrackRep *rep)
double getTimeSeed() const
void updateFitResults(gbl::GblTrajectory &traj)
Update fitter info from GBL fit results.
unsigned int hasMeasurement() const
Check for measurement at a point.
void setFitStatus(FitStatus *fitStatus, const AbsTrackRep *rep)
Defines for I/O streams used for error and debug printing.
void setTrackLen(double trackLen)
TVector3 getFieldVal(const TVector3 &position)
This does NOT use the cache!
TrackPoint * getPoint(int id) const
void setNFailedPoints(int nFailedPoints)
double getMass(const StateOnPlane &state) const
Get tha particle mass in GeV/c^2.
void setNdf(const double &ndf)
virtual void setTime(StateOnPlane &state, double time) const =0
Set time at which the state was defined.
void setIsFitted(bool fitted=true)
Abstract base class for a track representation.
void setSortingParameter(double sortingParameter)
void setJacobian(TMatrixD jacobian)
Set the Jacobian for further GblPoint construction.
FitStatus for use with GblFitter.
bool isValid() const
Retrieve validity of trajectory.
Namespace for the general broken lines package.
void setNumIterations(unsigned int numIterations)
Simple struct containing MaterialProperties and stepsize in the material.
double extrapolateToPlane(const SharedPlanePtr &plane, bool stopAtBoundary=false, bool calcJacobianNoise=false)
void setFitterInfo(genfit::AbsFitterInfo *fitterInfo)
Takes Ownership.
bool hasFitterInfo(const AbsTrackRep *rep) const
void insertPoint(TrackPoint *point, int id=-1)
Insert TrackPoint BEFORE TrackPoint with position id, if id >= 0.
TrackPoint * getPointWithMeasurement(int id) const
Object containing AbsMeasurement and AbsFitterInfo objects.
A state with arbitrary dimension defined in a DetPlane.
void setCurvature(bool useCurvature)
bool hasRawMeasurements() const
virtual SharedPlanePtr constructPlane(const StateOnPlane &state) const =0
TVector3 getDir(const StateOnPlane &state) const
Get the direction vector of a state.
AbsMeasurement * getRawMeasurement(int i=0) const
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,...
void setCharge(double charge)
void setIsFittedWithReferenceTrack(bool fittedWithReferenceTrack=true)
virtual double getCharge(const StateOnPlane &state) const =0
Get the (fitted) charge of a state. This is not always equal the pdg charge (e.g. if the charge sign ...
void setChi2(const double &chi2)
unsigned int getNumPointsWithMeasurement() const
gbl::GblPoint constructGblPoint()
Collect all data and create a GblPoint.
const SharedPlanePtr & getPlane() const
MaterialProperties materialProperties_
virtual void getForwardJacobianAndNoise(TMatrixD &jacobian, TMatrixDSym &noise, TVectorD &deltaState) const =0
Get the jacobian and noise matrix of the last extrapolation.
virtual void setPosMom(StateOnPlane &state, const TVector3 &pos, const TVector3 &mom) const =0
Set position and momentum of state.
static FieldManager * getInstance()
Get singleton instance.
void setIsFitConvergedPartially(bool fitConverged=true)
const TVectorD & getStateSeed() const
void setIsFitConvergedFully(bool fitConverged=true)
virtual unsigned int getDim() const =0
Get the dimension of the state vector used by the track representation.
const MeasuredStateOnPlane & getFittedState(int id=0, const AbsTrackRep *rep=NULL, bool biased=true) const
Shortcut to get FittedStates.
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, std::string optionList="")
Perform fit of trajectory.
void setScatterer(ThinScatterer *scatterer)
bool sort()
Sort TrackPoint and according to their sorting parameters.
AbsFitterInfo * getFitterInfo(const AbsTrackRep *rep=NULL) const
Get fitterInfo for rep. Per default, use cardinal rep.
unsigned int getNumPoints() const
double getSortingParameter() const
Material properties needed e.g. for material effects calculation.
bool hasScatterer() const
Check for scatterer at a point.