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;
TVector3 getDir(const StateOnPlane &state) const
Get the direction vector of a state.
gbl::GblPoint constructGblPoint()
Collect all data and create a GblPoint.
void setScatterer(ThinScatterer *scatterer)
virtual void getForwardJacobianAndNoise(TMatrixD &jacobian, TMatrixDSym &noise, TVectorD &deltaState) const =0
Get the jacobian and noise matrix of the last extrapolation.
bool hasScatterer() const
boost::shared_ptr< genfit::DetPlane > SharedPlanePtr
Shared Pointer to a DetPlane.
void setCurvature(bool useCurvature)
virtual SharedPlanePtr constructPlane(const StateOnPlane &state) const =0
void setSortingParameter(double sortingParameter)
FitStatus for use with GblFitter.
double getMass(const StateOnPlane &state) const
Get tha particle mass in GeV/c^2.
TrackSegmentController for use with GblFitter.
bool hasFitterInfo(const AbsTrackRep *rep) const
void setIsFittedWithReferenceTrack(bool fittedWithReferenceTrack=true)
const MeasuredStateOnPlane & getFittedState(int id=0, const AbsTrackRep *rep=NULL, bool biased=true) const
Shortcut to get FittedStates.
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, the extrapolated state.
void insertPoint(TrackPoint *point, int id=-1)
Insert TrackPoint BEFORE TrackPoint with position id, if id >= 0.
unsigned int getNumPoints() const
unsigned int hasMeasurement() const
Check for measurement at a point.
virtual void setPosMom(StateOnPlane &state, const TVector3 &pos, const TVector3 &mom) const =0
Set position and momentum of state.
TVector3 getFieldVal(const TVector3 &position)
This does NOT use the cache!
Collects information needed and produced by a GblFitter/GBL and is specific to one AbsTrackRep of the...
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...
double getTimeSeed() const
void setNFailedPoints(int nFailedPoints)
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...
AbsMeasurement * getRawMeasurement(int i=0) const
unsigned int getNumPointsWithMeasurement() const
virtual std::vector< genfit::MatStep > getSteps() const =0
Get stepsizes and material properties of crossed materials of the last extrapolation.
static FieldManager * getInstance()
Get singleton instance.
void setIsFitted(bool fitted=true)
Simple struct containing MaterialProperties and stepsize in the material.
virtual void setTime(StateOnPlane &state, double time) const =0
Set time at which the state was defined.
TrackPoint * getPoint(int id) const
bool sort()
Sort TrackPoint and according to their sorting parameters.
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, std::string optionList="")
Perform fit of trajectory.
Namespace for the general broken lines package.
void setTrackLen(double trackLen)
Abstract base class for a track representation.
void updateFitResults(gbl::GblTrajectory &traj)
Update fitter info from GBL fit results.
void setNdf(const double &ndf)
void setNumIterations(unsigned int numIterations)
double getSortingParameter() const
Object containing AbsMeasurement and AbsFitterInfo objects.
double extrapolateToPlane(const SharedPlanePtr &plane, bool stopAtBoundary=false, bool calcJacobianNoise=false)
void setCharge(double charge)
Material properties needed e.g. for material effects calculation.
const SharedPlanePtr & getPlane() const
A state with arbitrary dimension defined in a DetPlane.
TrackPoint * getPointWithMeasurement(int id) const
virtual double getMomMag(const StateOnPlane &state) const =0
get the magnitude of the momentum in GeV.
bool hasRawMeasurements() const
void recalculateJacobian(GblFitterInfo *prevFitterInfo)
Re-extrapolates between prevFitterInfo and this point using forward state to update the Jacobian (if ...
Collection of TrackPoint objects, AbsTrackRep objects and FitStatus objects.
bool isValid() const
Retrieve validity of trajectory.
MaterialProperties materialProperties_
void setIsFitConvergedPartially(bool fitConverged=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 setIsFitConvergedFully(bool fitConverged=true)
AbsFitterInfo * getFitterInfo(const AbsTrackRep *rep=NULL) const
Get fitterInfo for rep. Per default, use cardinal rep.
void setChi2(const double &chi2)
Defines for I/O streams used for error and debug printing.
void setFitStatus(FitStatus *fitStatus, const AbsTrackRep *rep)
void setFitterInfo(genfit::AbsFitterInfo *fitterInfo)
Takes Ownership.
virtual unsigned int getDim() const =0
Get the dimension of the state vector used by the track representation.
void deleteFitterInfo(const AbsTrackRep *rep)
const TVectorD & getStateSeed() const
void setJacobian(TMatrixD jacobian)
Set the Jacobian for further GblPoint construction.