TrkRUtil Namespace Reference

Functions
void	PredictXY (const double *x1, const double *sig1, const double *x2, const double *sig2, double *x3, double *sig3)
	Predicts what XY position of a straight line would be if points x1 and x2 are on a straight line.
double	FitLineFromWires (int n, const double *u, const double *w, const double *z, const double *cosA, const double *sinA, double z1, double z2, double par[4], double parErr[4])
	Fits a straight line to a set of wires.
double	FitCurvedTrack (int n, const double *u, const double *w, const double *z, const double *cosA, const double *sinA, double z0, const double *lambdaX, const double *lambdaY, double par[5], double parErr[5])
	This function does a non-iterative fit of particle momentum by assuming that integral of integral of Bdl over the particle flight path is known a priori.
double	FitCurvedTrack (int n, const double *u, const double *w, const double *z, const double *cosA, const double *sinA, double z0, const double *lambda, double par[5], double parErr[5])
	This function does a non-iterative fit of particle momentum by assuming that integral of integral of Bdl over the particle flight path is known a priori.
bool	CalcLambda (double z0, double lx[][4], double ly[][4], double *par, int minCham, int maxCham)
	This function calculates coefficients lambda for wire plane z-positions.

---

Function Documentation

void TrkRUtil::PredictXY	(	const double *	x1,
		const double *	sig1,
		const double *	x2,
		const double *	sig2,
		double *	x3,
		double *	sig3
	)

Predicts what XY position of a straight line would be if points x1 and x2 are on a straight line.

We write the slope in each coordinate as a=(x1-x2)/(z1-z2) x3 = a*(z-z1) + x1=(x1*(z-z2)-x2*(z-z1))/(z1-z2); sig3^2 = ((z-z2)/(z1-z2))^2*sig1^2 + ((z-z1)/(z1-z2))^2*sig2^2

Definition at line 29 of file TrkRUtil.cxx.

Referenced by TrkSegBuilder::AddPoint(), AlignBC::FillNt(), TrackSeg::GetXY(), main(), ResolveDAF(), and TrkSegBuilder::TryFit().

{
  // First compute x position
  double zDiff = 1/(x1[2] - x2[2]);
  double zDiff1 = (x3[2] - x1[2]) * zDiff;
  double zDiff2 = (x3[2] - x2[2]) * zDiff;

  for (int i = 0; i < 2; ++i) {
    x3[i] = x1[i]*zDiff2 - x2[i]*zDiff1;
    if (sig3) sig3[i] = sqrt(pow(zDiff2*sig1[i], 2) + 
                 pow(zDiff1*sig2[i], 2));
  }
}

double TrkRUtil::FitLineFromWires	(	int	n,
		const double *	u,
		const double *	w,
		const double *	z,
		const double *	cosA,
		const double *	sinA,
		double	z1,
		double	z2,
		double	par[4],
		double	parErr[4]
	)

Fits a straight line to a set of wires.

The line should be given by 4 parameters: point (x1, y2) at z1 and point (x2, y2) at z2. Array par will be returned as (x1, y2, x2, y2) The function will return chi squared of the fit.

Parameters:

	n	- number of chamber wires
	u	- position of every wire in its coordinate (see GChamGeo)
	w	- weight of every wire
	z	- z position of every wire
	cosA	- cosine of the wire angle w.r.t. vertical
	sinA	- sine of the wire angle w.r.t. vertical
	z1	- position of the first point on the track
	z2	- position of the last point on the track
	par	- returned fit parameters
	parErr	- returned fit parameter errors

Definition at line 61 of file TrkRUtil.cxx.

References ok.

Referenced by AlignChamZNoBF::ComputeChi2F(), TrackBeamPart::FillDiffHisto(), TrackSeg::Fit(), TrkSegFitter::FitTrack(), main(), TrkSegBuilder::MakeTrack(), ResolveDAF(), ResolveDumb(), and TrkSegBuilder::TryFit().

{
  const int nParam = 4;
  if (n < nParam) {
    cout << "TrkRUtil::FitLineFromWires() will not work with less "
     << " than " << nParam << " points" << endl;
    return 1e9;
  }
  double mat[nParam * nParam];
  memset(mat, 0, sizeof(mat));
  double vec[nParam] = {0, 0, 0, 0};

  // Calculate matrix elements, matrix is symmetric, so off-diagonal
  // elements need to be calculated only once
  double dz = 1/ (z1 - z2);
  double dz1, dz2;
  double c[nParam];
  double chi2 = 0;
  for (int pt = 0; pt < n; ++pt) {
    dz1 = (z[pt] - z1) * dz;
    dz2 = (z[pt] - z2) * dz;
    chi2 += w[pt] * u[pt] * u[pt];
    
    // Initialize geometric constants
    c[0] = dz2 * cosA[pt];
    c[1] = - dz2 * sinA[pt];
    c[2] = - dz1 * cosA[pt];
    c[3] = dz1 * sinA[pt];

    for (int i = 0; i < nParam; ++i) {
      vec[i] += w[pt] * u[pt] * c[i];

      for (int j = i; j < nParam; ++j) {
    mat[i * nParam + j] += w[pt] * c[i] * c[j];
    if (i != j) mat[j * nParam + i] = mat[i * nParam + j];
      }
    }
  }
  
  Bool_t ok = true;

  for (int j = 0; j < nParam; ++j) {
    if (mat[j * nParam + j] < 0) {
      ok = false;
      break;
    }
  }

  if (!ok) {
    cout << "TrkRUtil::FitLineFromWires() : got a goofy matrix" << endl;
    for (int j = 0; j < nParam; ++j) {
      for (int k = 0; k < nParam; ++k) cout << mat[j * nParam + k] << " ";
      cout << endl;
    }
    cout << " **** From " << z1 << " to " << z2 << " **** " << endl;
    for (int i = 0; i < n; ++i) {
      cout << i << ": w=" << w[i] << ", u=" << u[i] << ", z=" << z[i]
       << ", cos=" << cosA[i] << ", sin=" << sinA[i] << endl;
    }
    cout << "***********************************************" << endl;
  }

  // Define matrices so as to make ROOT invert the matrix for us
  TMatrixDSym matrix(nParam, mat);
  TVectorD v(nParam, vec);
  TDecompLU lu(matrix);
  if (!lu.Decompose() || !lu.Solve(v)) {
    cout << __FILE__ << ":" << __LINE__ << " error Solving matrix"
     << endl;
    for (int i = 0; i < 4; ++i) {
      par[i] = 1e9;
      if (parErr) parErr[i] = 1e9;
    }
    return 1e9;
  }
  if (parErr) {
    TMatrixD matInv(nParam, nParam);
    lu.Invert(matInv);

    for (int i = 0; i < nParam; ++i) {
      parErr[i] = matInv[i][i];
      if (parErr[i] < 0) {
    cout << "TrkRUtil::FitLineFromWires(): parErr[" << i << "]: "
         << parErr[i] << " this shouldn't happen. Matrix" << endl;
    for (int j = 0; j < nParam; ++j) {
      for (int k = 0; k < nParam; ++k) cout << mat[j * nParam + k] << " ";
      cout << endl;
    }
    cout << " **** From " << z1 << " to " << z2 << " **** " << endl;
    for (int j = 0; j < n; ++j) {
      cout << j << ": w=" << w[j] << ", u=" << u[j] << ", z=" << z[j]
           << ", cos=" << cosA[j] << ", sin=" << sinA[j] << endl;
    }
    cout << "***********************************************" << endl;
    parErr[i] = 0;
      }
      else parErr[i] = sqrt(parErr[i]);
    }
  }

  // Calculate chi squared
  for (int i = 0; i < nParam; ++i) {
    par[i] = v[i];
    chi2 -= vec[i] * par[i];
  }

  // Check for floating point errors
  if (chi2 < 0) return 0;
  
  return chi2;
}

double TrkRUtil::FitCurvedTrack	(	int	n,
		const double *	u,
		const double *	w,
		const double *	z,
		const double *	cosA,
		const double *	sinA,
		double	z0,
		const double *	lambdaX,
		const double *	lambdaY,
		double	par[5],
		double	parErr[5]
	)

This function does a non-iterative fit of particle momentum by assuming that integral of integral of Bdl over the particle flight path is known a priori.

Parameters:

	n	- number of measurement points
	u	- measurement value in respective coordinate
	w	- weight of the measurement
	z	- z position of the measurement
	cosA	- cosine of the wire angle w.r.t. vertical
	sinA	- sine of the wire angle w.r.t. vertical
	z0	- z coordinate where momentum direction is fitted
	lambdaBx	- x-displacement from straight line projection
	lambdaBy	- y-displacement from straight line projection
	par	- returned fit parameters: x, y, px, py, q/p at z=z0
	parErr	- returned fit parameter errors

Definition at line 248 of file TrkRUtil.cxx.

References FitLine().

Referenced by TrkCandBuilder::AddBC(), ChamFCN(), AlignChamZ::ComputeChi2F(), TrackBeamPart::ComputeF(), AlignRotB::ComputeF(), AlignTPCDrift2::ComputeTrkChi2(), TrackBeamPart::FillDiffHisto(), Fit(), TrkCand::Fit(), SPTrk::Fit(), AlignRotB::Fit(), AlignChamZ::Fit(), AlignChamW0::Fit(), FitEqualWeights(), FitMSC(), main(), TrkCandBuilder::MakeTrack(), and TrkCandBuilder::TryBCFit().

{
  const int nParam = 5;
  if (n < nParam) {
    cout << "TrkRUtil::FitCurvedTrack() will not work with "
     << n << " points; need at least " << nParam << endl;
    return 1e9;
  }
  double mat[nParam * nParam];
  memset(mat, 0, sizeof(mat));
  double vec[nParam] = {0, 0, 0, 0, 0};

  double chi2 = 0;
  double c[nParam]; // Geometric constants for each parameter

  double minZ = z[0];
  double maxZ = z[0];
  double maxLam = 0;

  const double kScale[nParam] = {1, 1, 1e-2, 1e-2, 1e-2};

  for (int pt = 0; pt < n; ++pt) {
    if (z[pt] < minZ) minZ = z[pt];
    else if (z[pt] > maxZ) maxZ = z[pt];

    // Initialize geometric constants
    c[0] = cosA[pt] * kScale[0];
    c[1] = -sinA[pt] * kScale[1];
    c[2] = cosA[pt] * (z[pt] - z0) * kScale[2];
    c[3] = -sinA[pt] * (z[pt] - z0) * kScale[3];
    c[4] = (lambdaX[pt] * cosA[pt] - lambdaY[pt] * sinA[pt]) * kScale[4];
    if (fabs(c[4]) > maxLam) maxLam = fabs(c[4]);

    chi2 += w[pt] * u[pt] * u[pt];

    for (int i = 0; i < nParam; ++i) {
      vec[i] += w[pt] * u[pt] * c[i];

      for (int j = i; j < nParam; ++j) {
    mat[i * nParam + j] += w[pt] * c[i] * c[j];
    if (i != j) mat[j * nParam + i] = mat[i * nParam + j];
      }
    }
  }

  double maxEnt = 0, minEnt = 1e15;
  for (int i = 0; i < 25; ++i) {
    if (fabs(maxEnt) < fabs(mat[i])) maxEnt = mat[i];
    if (fabs(minEnt) > fabs(mat[i])) minEnt = mat[i];
  }

  if (maxLam < 1e-6) return FitLine(chi2, mat, vec, par, parErr);

  // Use ROOT linear algebra to solve matrix inversion
  TMatrixDSym matrix(nParam, mat);
  TVectorD v(nParam, vec);
  TDecompLU lu(matrix);
  if (!lu.Decompose()) return FitLine(chi2, mat, vec, par, parErr);
  double d1, d2;
  lu.Det(d1, d2);
  if (d1 <= 0) {
//     cout << "TrkRUtil.cxx:" << __LINE__
//   << " -- Determinant is less than 0" << endl;
    return FitLine(chi2, mat, vec, par, parErr);
  }
  if (!lu.Solve(v)) {
    for (int i = 0; i < nParam; ++i) par[i] = -1e9;
    return 1e9;
  }

  bool goodres = true;

  if (parErr) {
    TMatrixD matInv(nParam, nParam);
    lu.Invert(matInv);

    for (int i = 0; i < nParam; ++i) {
      parErr[i] = matInv[i][i] * kScale[i];
      if (parErr[i] < 0) {
    cout << __FILE__ << ":" << __LINE__ << " parErr[" << i << "]: "
         << parErr[i] << " this shouldn't happen!" << endl;
    parErr[i] = 0;
    goodres = false;
      }
      else parErr[i] = sqrt(parErr[i]);

    }
    if (!goodres) {
      cout << "Matrix:\n";
      for (int i = 0; i < 25; ++i) {
    printf(" %+.4e, ", mat[i]);
    if (i % 5 == 4) cout << "\n";
      }
      cout << "Vector: ";
      for (int i = 0; i < 5; ++i) printf("%+.4e, ", vec[i]);
      cout << endl;

      cout << "Solution: ";
      for (int i = 0; i < 5; ++i) printf("%+.4e, ", v[i]);
      cout << endl;

      cout << "Data [" << n << "]\n";
      for (int i = 0; i < n; ++i) {
    printf("%+.3e: %+.3e %+.2f %+.2f %+.3e %+.3e\n",
           z[i], u[i], cosA[i], sinA[i], lambdaX[i], lambdaY[i]);
      }
    }
  }

  for (int i = 0; i < nParam; ++i) {
    par[i] = v[i];
    chi2 -= vec[i] * par[i];
    par[i] *= kScale[i];
  }
  if (chi2 < 0) return 0;
  return chi2;

}

double TrkRUtil::FitCurvedTrack	(	int	n,
		const double *	u,
		const double *	w,
		const double *	z,
		const double *	cosA,
		const double *	sinA,
		double	z0,
		const double *	lambda,
		double	par[5],
		double	parErr[5]
	)

This function does a non-iterative fit of particle momentum by assuming that integral of integral of Bdl over the particle flight path is known a priori.

Parameters:

	n	- number of measurement points
	u	- measurement value in respective coordinate
	w	- weight of the measurement
	z	- z position of the measurement
	cosA	- cosine of the wire angle w.r.t. vertical
	sinA	- sine of the wire angle w.r.t. vertical
	z0	- z coordinate where momentum direction is fitted
	lambda	- displacement from straight-line projection in that wire direction
	par	- returned fit parameters: x, y, px, py, q/p at z=z0
	parErr	- returned fit parameter errors

Definition at line 390 of file TrkRUtil.cxx.

References FitLine().

{
  const int nParam = 5;
  if (n < nParam) {
    cout << "TrkRUtil::FitCurvedTrack() will not work with "
     << n << " points; need at least " << nParam << endl;
    return 1e9;
  }
  double mat[nParam * nParam];
  memset(mat, 0, sizeof(mat));
  double vec[nParam] = {0, 0, 0, 0, 0};

  double chi2 = 0;
  double c[nParam]; // Geometric constants for each parameter

  double minZ = z[0];
  double maxZ = z[0];
  double maxLam = 0;

  for (int pt = 0; pt < n; ++pt) {
    if (z[pt] < minZ) minZ = z[pt];
    else if (z[pt] > maxZ) maxZ = z[pt];

    // Initialize geometric constants
    c[0] = cosA[pt];
    c[1] = -sinA[pt];
    c[2] = cosA[pt] * (z[pt] - z0);
    c[3] = -sinA[pt] * (z[pt] - z0);
    c[4] = lambda[pt];
    if (fabs(c[4]) > maxLam) maxLam = fabs(c[4]);

    chi2 += w[pt] * u[pt] * u[pt];

    for (int i = 0; i < nParam; ++i) {
      vec[i] += w[pt] * u[pt] * c[i];

      for (int j = i; j < nParam; ++j) {
    mat[i * nParam + j] += w[pt] * c[i] * c[j];
    if (i != j) mat[j * nParam + i] = mat[i * nParam + j];
      }
    }
  }

  if (maxLam < 1e-6) return FitLine(chi2, mat, vec, par, parErr);

  // Use ROOT linear algebra to solve matrix inversion
  TMatrixDSym matrix(nParam, mat);
  TVectorD v(nParam, vec);
  TDecompLU lu(matrix);
  if (!lu.Decompose() || !lu.Solve(v)) {
    for (int i = 0; i < nParam; ++i) par[i] = -1e9;
    return 1e9;
  }

  if (parErr) {
    TMatrixD matInv(nParam, nParam);
    lu.Invert(matInv);

    for (int i = 0; i < nParam; ++i) {
      parErr[i] = matInv[i][i];
      if (parErr[i] < 0) {
    cout << __FILE__ << ":" << __LINE__ << " parErr[" << i << "]: "
         << parErr[i] << " this shouldn't happen!" << endl;
    parErr[i] = 0;
      }
      else parErr[i] = sqrt(parErr[i]);

    }
  }

  for (int i = 0; i < nParam; ++i) {
    par[i] = v[i];
    chi2 -= vec[i] * par[i];
  }
  if (chi2 < 0) return 0;
  return chi2;

}

bool TrkRUtil::CalcLambda	(	double	z0,
		double	lx[][4],
		double	ly[][4],
		double *	par = 0,
		int	minCham = 0,
		int	maxCham = 8
	)

This function calculates coefficients lambda for wire plane z-positions.

These coefficients are needed for track template and can be thought of as integral of integral of Bdl or equivalent displacement from straight-line trajectory for 1 GeV/c particle. The values are stored for each plane of each wire chamber.

Definition at line 482 of file TrkRUtil.cxx.

References geo, SwimMIPP::GetXY(), GMIPPGeo::Instance(), kNPlane, mom, p, and SwimMIPP::Swim().

Referenced by ChamFCN(), AlignChamZ::ComputeChi2F(), TrackBeamPart::ComputeF(), AlignRotB::ComputeF(), AlignChamZ::ComputeF(), AlignChamZ::EndRun(), TrackBeamPart::FillDiffHisto(), Fit(), AlignRotB::Fit(), AlignChamZ::Fit(), AlignChamW0::Fit(), FitEqualWeights(), FitMSC(), TrkCand::Init(), main(), TrkCandBuilder::MakeTrack(), TrackBeamPart::NewRun(), and TrkCandBuilder::TryBCFit().

{
  static SwimMIPP swim;

  double mom = 120;
  if (par) {
    if (fabs(par[4]) < 1e-3) par[4] = 1e-3;
    mom = 1 / par[4];
  }

  if (!par) swim.Swim(1, 120, 0, 0, z0, 0, 0);
  else swim.Swim(z0, par);    

  for (int c = minCham; c <= maxCham; ++c) {
    GChamGeo& geo = GMIPPGeo::Instance().Cham(c);
    for (int p = 0; p < kNPlane; ++p) {
      // Get (x,y) coordinate prediction at chamber plane
      if (!swim.GetXY(geo.Z(p + 1), &lx[c][p], &ly[c][p])) return false;
    
      if (par) {
    // Subtract initial parameters if we have them
    lx[c][p] -= par[0] + par[2] * (geo.Z(p + 1) - z0);
    ly[c][p] -= par[1] + par[3] * (geo.Z(p + 1) - z0);
      }
      // Multiply by momentum to obtain coefficients necessary for
      // track template fitting
      lx[c][p] *= mom; 
      ly[c][p] *= mom;
    }
  }
  return true;
}

---