60 #include "../include/GblFitStatus.h"
75 #include <Math/SMatrix.h>
77 #include <TVectorDfwd.h>
85 using namespace genfit;
90 GblFitter::~GblFitter() {
91 if (m_segmentController) {
92 delete m_segmentController;
93 m_segmentController =
nullptr;
99 if (m_segmentController) {
100 delete m_segmentController;
101 m_segmentController =
nullptr;
103 m_segmentController = controler;
108 cleanGblInfo(trk, rep);
115 bool fitQoverP =
true;
118 if (!(Bfield > 1.
e-16))
125 double lostWeight = 0.;
139 constructGblInfo(trk, rep)
146 if (m_externalIterations < 1)
151 unsigned int nFailed = 0;
154 std::vector<std::string> gblIterations;
155 gblIterations.push_back(m_gblInternalIterations);
159 for (
unsigned int iIter = 0; iIter < m_externalIterations; iIter++) {
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++) {
174 fitRes = traj.
fit(Chi2, Ndf, lostWeight, (iIter == m_externalIterations - 1) ? m_gblInternalIterations :
"");
177 updateGblInfo(traj, trk, rep);
181 if (m_recalcJacobians > iIter) {
184 for (
unsigned int ip = 0; ip < trk->
getNumPoints(); ip++) {
188 prevFitterInfo = currFitterInfo;
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();
214 cout <<
" (" << ndummy <<
" dummy) ";
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;
255 double arcLenPos = 0;
258 for (
int ipoint_meas = 0; ipoint_meas < npoints_meas - 1; ipoint_meas++) {
276 std::vector<gbl::GblPoint> thePoints;
280 for (
unsigned int ip = 0; ip < trk->
getNumPoints(); ip++) {
296 for (
unsigned int ip = 0; ip < trk->
getNumPoints(); ip++) {
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;
335 for (
unsigned int i = 0;
i < steps.size();
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.0
e-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, scatTheta,
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;