1 added + 5 modified, total 6 files
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLStripClusterData.java 2014-09-17 17:45:13 UTC (rev 1035)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLStripClusterData.java 2014-09-17 18:01:56 UTC (rev 1036)
@@ -10,6 +10,7 @@
*
* @author phansson
*
+ * @version $Id:
*/
public class GBLStripClusterData implements GenericObject {
@@ -61,6 +62,23 @@
public GBLStripClusterData(int id) {
setId(id);
}
+
+ /*
+ * Constructor from GenericObject
+ * TODO add size checks for backwards compatability
+ */
+ public GBLStripClusterData(GenericObject o)
+ {
+ for(int i=0; i<GBLINT.BANK_INT_SIZE; ++i)
+ {
+ bank_int[i] = o.getIntVal(i);
+ }
+ for(int i=0; i<GBLDOUBLE.BANK_DOUBLE_SIZE; ++i)
+ {
+ bank_double[i] = o.getDoubleVal(i);
+ }
+
+ }
/**
* @param set track id to val
@@ -263,11 +281,6 @@
return getDoubleVal(GBLDOUBLE.MSANGLE);
}
-
-
-
-
-
/*
* The functions below are all overide from
* @see org.lcsim.event.GenericObject#getNInt()
@@ -301,7 +314,4 @@
return false;
}
-
-
-
}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java 2014-09-17 17:45:13 UTC (rev 1035)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLTrackData.java 2014-09-17 18:01:56 UTC (rev 1036)
@@ -3,6 +3,11 @@
import org.hps.recon.tracking.gbl.GBLOutput.PerigeeParams;
import org.lcsim.event.GenericObject;
+/**
+ *
+ *
+ * @version $Id:
+ */
public class GBLTrackData implements GenericObject {
/*
@@ -34,6 +39,23 @@
public GBLTrackData(int id) {
setTrackId(id);
}
+
+ /*
+ * Constructor from GenericObject
+ * TODO add size checks for backwards compatability
+ */
+ public GBLTrackData(GenericObject o)
+ {
+ for(int i=0; i<GBLINT.BANK_INT_SIZE; ++i)
+ {
+ bank_int[i] = o.getIntVal(i);
+ }
+ for(int i=0; i<GBLDOUBLE.BANK_DOUBLE_SIZE; ++i)
+ {
+ bank_double[i] = o.getDoubleVal(i);
+ }
+
+ }
/**
* @param set track id to val
@@ -104,5 +126,4 @@
return false;
}
-
-}
+}
\ No newline at end of file
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblData.java 2014-09-17 17:45:13 UTC (rev 1035)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblData.java 2014-09-17 18:01:56 UTC (rev 1036)
@@ -2,9 +2,9 @@
import java.util.ArrayList;
import java.util.List;
+import org.hps.recon.tracking.gbl.matrix.Matrix;
import org.hps.recon.tracking.gbl.matrix.VVector;
-
/**
*
* @author Norman A Graf
@@ -40,92 +40,109 @@
thePrediction = 0.;
}
-///// Add derivatives from measurement.
-///**
-// * Add (non-zero) derivatives to data block. Fill list of labels of used fit parameters.
-// * \param [in] iRow Row index (0-4) in up to 5D measurement
-// * \param [in] labDer Labels for derivatives
-// * \param [in] matDer Derivatives (matrix) 'measurement vs track fit parameters'
-// * \param [in] iOff Offset for row index for additional parameters
-// * \param [in] derLocal Derivatives (matrix) for additional local parameters
-// * \param [in] labGlobal Labels for additional global (MP-II) parameters
-// * \param [in] derGlobal Derivatives (matrix) for additional global (MP-II) parameters
-// * \param [in] extOff Offset for external parameters
-// * \param [in] extDer Derivatives for external Parameters
-// */
-//void addDerivatives(unsigned int iRow,
-// const std::vector<unsigned int> &labDer, const SMatrix55 &matDer,
-// unsigned int iOff, const TMatrixD &derLocal,
-// const std::vector<int> &labGlobal, const TMatrixD &derGlobal,
-// unsigned int extOff, const TMatrixD &extDer) {
-//
-// unsigned int nParMax = 5 + derLocal.GetNcols() + extDer.GetNcols();
+/// Add derivatives from measurement.
+ /**
+ * Add (non-zero) derivatives to data block. Fill list of labels of used fit
+ * parameters. \param [in] iRow Row index (0-4) in up to 5D measurement
+ * \param [in] labDer Labels for derivatives \param [in] matDer Derivatives
+ * (matrix) 'measurement vs track fit parameters' \param [in] iOff Offset
+ * for row index for additional parameters \param [in] derLocal Derivatives
+ * (matrix) for additional local parameters \param [in] labGlobal Labels for
+ * additional global (MP-II) parameters \param [in] derGlobal Derivatives
+ * (matrix) for additional global (MP-II) parameters \param [in] extOff
+ * Offset for external parameters \param [in] extDer Derivatives for
+ * external Parameters
+ */
+ void addDerivatives(int iRow,
+ List<Integer> labDer, Matrix matDer,
+ int iOff, Matrix derLocal,
+ List<Integer> labGlobal, Matrix derGlobal,
+ int extOff, Matrix extDer)
+ {
+ int nLocal = 0;
+ int nExt = 0;
+ if (derLocal != null) {
+ nLocal = derLocal.getColumnDimension();
+ }
+ if (extDer != null) {
+ nExt = extDer.getColumnDimension();
+ }
+ int nParMax = 5 + nLocal + nExt;
// theParameters.reserve(nParMax); // have to be sorted
// theDerivatives.reserve(nParMax);
-//
-// for (int i = 0; i < derLocal.GetNcols(); ++i) // local derivatives
-// {
-// if (derLocal(iRow - iOff, i)) {
-// theParameters.push_back(i + 1);
-// theDerivatives.push_back(derLocal(iRow - iOff, i));
-// }
-// }
-//
-// for (int i = 0; i < extDer.GetNcols(); ++i) // external derivatives
-// {
-// if (extDer(iRow - iOff, i)) {
-// theParameters.push_back(extOff + i + 1);
-// theDerivatives.push_back(extDer(iRow - iOff, i));
-// }
-// }
-//
-// for (unsigned int i = 0; i < 5; ++i) // curvature, offset derivatives
-// {
-// if (labDer[i] && matDer(iRow, i)) {
-// theParameters.push_back(labDer[i]);
-// theDerivatives.push_back(matDer(iRow, i));
-// }
-// }
-//
-// globalLabels = labGlobal;
-// for (int i = 0; i < derGlobal.GetNcols(); ++i) // global derivatives
-// globalDerivatives.push_back(derGlobal(iRow - iOff, i));
-//}
-//
-///// Add derivatives from kink.
-///**
-// * Add (non-zero) derivatives to data block. Fill list of labels of used fit parameters.
-// * \param [in] iRow Row index (0-1) in 2D kink
-// * \param [in] labDer Labels for derivatives
-// * \param [in] matDer Derivatives (matrix) 'kink vs track fit parameters'
-// * \param [in] extOff Offset for external parameters
-// * \param [in] extDer Derivatives for external Parameters
-// */
-//void addDerivatives(unsigned int iRow,
-// const std::vector<unsigned int> &labDer, const SMatrix27 &matDer,
-// unsigned int extOff, const TMatrixD &extDer) {
-//
-// unsigned int nParMax = 7 + extDer.GetNcols();
+
+ if (derLocal != null) {
+ for (int i = 0; i < derLocal.getColumnDimension(); ++i) // local derivatives
+ {
+ if (derLocal.get(iRow - iOff, i) != 0) {
+ theParameters.add(i + 1);
+ theDerivatives.add(derLocal.get(iRow - iOff, i));
+ }
+ }
+ }
+
+ if (extDer != null) {
+ for (int i = 0; i < extDer.getColumnDimension(); ++i) // external derivatives
+ {
+ if (extDer.get(iRow - iOff, i) > 0) {
+ theParameters.add(extOff + i + 1);
+ theDerivatives.add(extDer.get(iRow - iOff, i));
+ }
+ }
+ }
+ for (int i = 0; i < 5; ++i) // curvature, offset derivatives
+ {
+ if (labDer.get(i) != 0 && matDer.get(iRow, i) != 0) {
+ theParameters.add(labDer.get(i));
+ theDerivatives.add(matDer.get(iRow, i));
+ }
+ }
+
+ globalLabels = labGlobal;
+ for (int i = 0; i < derGlobal.getColumnDimension(); ++i) // global derivatives
+ {
+ globalDerivatives.add(derGlobal.get(iRow - iOff, i));
+ }
+ }
+
+/// Add derivatives from kink.
+ /**
+ * Add (non-zero) derivatives to data block. Fill list of labels of used fit
+ * parameters. \param [in] iRow Row index (0-1) in 2D kink \param [in]
+ * labDer Labels for derivatives \param [in] matDer Derivatives (matrix)
+ * 'kink vs track fit parameters' \param [in] extOff Offset for external
+ * parameters \param [in] extDer Derivatives for external Parameters
+ */
+ void addDerivatives(int iRow,
+ List<Integer> labDer, Matrix matDer,
+ int extOff, Matrix extDer)
+ {
+ int nExtDer = 0;
+ if (extDer != null) {
+ nExtDer = extDer.getColumnDimension();
+ }
+ int nParMax = 7 + nExtDer;
// theParameters.reserve(nParMax); // have to be sorted
// theDerivatives.reserve(nParMax);
-//
-// for (int i = 0; i < extDer.GetNcols(); ++i) // external derivatives
-// {
-// if (extDer(iRow, i)) {
-// theParameters.push_back(extOff + i + 1);
-// theDerivatives.push_back(extDer(iRow, i));
-// }
-// }
-//
-// for (unsigned int i = 0; i < 7; ++i) // curvature, offset derivatives
-// {
-// if (labDer[i] && matDer(iRow, i)) {
-// theParameters.push_back(labDer[i]);
-// theDerivatives.push_back(matDer(iRow, i));
-// }
-// }
-//}
-//
+
+ if (extDer != null) {
+ for (int i = 0; i < extDer.getColumnDimension(); ++i) // external derivatives
+ {
+ if (extDer.get(iRow, i) != 0) {
+ theParameters.add(extOff + i + 1);
+ theDerivatives.add(extDer.get(iRow, i));
+ }
+ }
+ }
+ for (int i = 0; i < 7; ++i) // curvature, offset derivatives
+ {
+ if (labDer.get(i) != 0 && matDer.get(iRow, i) != 0) {
+ theParameters.add(labDer.get(i));
+ theDerivatives.add(matDer.get(iRow, i));
+ }
+ }
+ }
+
///// Add derivatives from external seed.
///**
// * Add (non-zero) derivatives to data block. Fill list of labels of used fit parameters.
@@ -191,23 +208,23 @@
return aDiff * aDiff * thePrecision * theDownWeight;
}
-///// Print data block.
-//void printData() const {
-//
-// std::cout << " measurement at label " << theLabel << ": " << theValue
-// << ", " << thePrecision << std::endl;
-// std::cout << " param " << theParameters.size() << ":";
-// for (unsigned int i = 0; i < theParameters.size(); ++i)
-// std::cout << " " << theParameters[i];
-// }
-// std::cout << std::endl;
-// std::cout << " deriv " << theDerivatives.size() << ":";
-// for (unsigned int i = 0; i < theDerivatives.size(); ++i) {
-// std::cout << " " << theDerivatives[i];
-// }
-// std::cout << std::endl;
-//}
-//
+/// Print data block.
+ void printData()
+ {
+ System.out.println(" measurement at label " + theLabel + ": " + theValue
+ + ", " + thePrecision);
+ System.out.print(" param " + theParameters.size() + ":");
+ for (int i = 0; i < theParameters.size(); ++i) {
+ System.out.print(" " + theParameters.get(i));
+ }
+ System.out.println("\n");
+ System.out.print(" deriv " + theDerivatives.size() + ":");
+ for (int i = 0; i < theDerivatives.size(); ++i) {
+ System.out.print(" " + theDerivatives.get(i));
+ }
+ System.out.println("");
+ }
+
public String toString()
{
StringBuffer sb = new StringBuffer(" measurement at label " + theLabel + ": " + theValue + ", " + thePrecision + "\n");
@@ -216,6 +233,10 @@
sb.append(" " + theParameters.get(i));
}
sb.append("\n");
+ for (int i = 0; i < theDerivatives.size(); ++i) {
+ sb.append(" " + theDerivatives.get(i));
+ }
+ sb.append("\n");
return sb.toString();
}
@@ -236,16 +257,17 @@
void getLocalData(double[] retVal, int[] indLocal, double[] derLocal)
{
-
retVal[0] = theValue;
retVal[1] = thePrecision * theDownWeight;
- indLocal = new int[theParameters.size()];
- derLocal = new double[theParameters.size()];
for (int i = 0; i < theParameters.size(); ++i) {
indLocal[i] = theParameters.get(i);
derLocal[i] = theDerivatives.get(i);
}
+ }
+ int getNumParameters()
+ {
+ return theParameters.size();
}
//
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblPoint.java 2014-09-17 17:45:13 UTC (rev 1035)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblPoint.java 2014-09-17 18:01:56 UTC (rev 1036)
@@ -14,9 +14,11 @@
*
* @version $Id:
*/
-public class GblPoint {
+public class GblPoint
+{
- public GblPoint(hep.physics.matrix.BasicMatrix jacPointToPoint) {
+ public GblPoint(hep.physics.matrix.BasicMatrix jacPointToPoint)
+ {
theLabel = 0;
theOffset = 0;
measDim = 0;
@@ -34,70 +36,65 @@
}
}
- public void addMeasurement(hep.physics.matrix.Matrix proL2m, hep.physics.matrix.BasicMatrix meas, hep.physics.matrix.BasicMatrix measPrec) {
+ public void addMeasurement(hep.physics.matrix.Matrix proL2m, hep.physics.matrix.BasicMatrix meas, hep.physics.matrix.BasicMatrix measPrec)
+ {
int ncols = proL2m.getNColumns();
int nrows = proL2m.getNRows();
- System.out.println("proL2m has "+nrows+" rows and "+ ncols+ "columns");
+ System.out.println("proL2m has " + nrows + " rows and " + ncols + "columns");
Matrix a = new Matrix(nrows, ncols);
- for(int i=0; i<nrows; ++i)
- {
- for(int j=0; j<ncols; ++j)
- {
- a.set(i,j,proL2m.e(i, j));
+ for (int i = 0; i < nrows; ++i) {
+ for (int j = 0; j < ncols; ++j) {
+ a.set(i, j, proL2m.e(i, j));
}
}
System.out.println("GblPoint add matrix: ");
a.print(10, 6);
-
-
-
+
int measnrows = meas.getNRows();
int measncols = meas.getNColumns();
- System.out.println("meas has "+measnrows+" rows and "+measncols+" columns");
-
+ System.out.println("meas has " + measnrows + " rows and " + measncols + " columns");
+
Vector measvec = new Vector(measncols);
- for(int i=0; i<measnrows; ++i)
- {
+ for (int i = 0; i < measnrows; ++i) {
measvec.set(i, meas.e(0, i));
}
System.out.println("GblPoint add meas: ");
measvec.print(10, 6);
-
-
+
int measPrecnrows = measPrec.getNRows();
int measPrecncols = measPrec.getNColumns();
-
- System.out.println("measPrec has "+measPrecnrows+" rows and "+measPrecncols+" columns");
+
+ System.out.println("measPrec has " + measPrecnrows + " rows and " + measPrecncols + " columns");
Vector measPrecvec = new Vector(measPrecncols);
- for(int i=0; i<measPrecnrows; ++i)
- {
+ for (int i = 0; i < measPrecnrows; ++i) {
measPrecvec.set(i, measPrec.e(0, i));
}
System.out.println("GblPoint add measPrec: ");
- measPrecvec.print(10, 6);
-
+ measPrecvec.print(10, 6);
+
addMeasurement(a, measvec, measPrecvec, 0.);
}
- public void addScatterer(hep.physics.matrix.BasicMatrix scat, hep.physics.matrix.BasicMatrix scatPrec) {
+ public void addScatterer(hep.physics.matrix.BasicMatrix scat, hep.physics.matrix.BasicMatrix scatPrec)
+ {
// TODO Auto-generated method stub
}
private int theLabel; ///< Label identifying point
private int theOffset; ///< Offset number at point if not negative (else interpolation needed)
- private Matrix p2pJacobian = new SymMatrix(5); ///< Point-to-point jacobian from previous point
- private Matrix prevJacobian = new SymMatrix(5); ///< Jacobian to previous scatterer (or first measurement)
- private Matrix nextJacobian = new SymMatrix(5); ///< Jacobian to next scatterer (or last measurement)
+ private Matrix p2pJacobian = new Matrix(5, 5); ///< Point-to-point jacobian from previous point
+ private Matrix prevJacobian = new Matrix(5, 5); ///< Jacobian to previous scatterer (or first measurement)
+ private Matrix nextJacobian = new Matrix(5, 5); ///< Jacobian to next scatterer (or last measurement)
private int measDim; ///< Dimension of measurement (1-5), 0 indicates absence of measurement
- private SymMatrix measProjection = new SymMatrix(5); ///< Projection from measurement to local system
+ private Matrix measProjection = new Matrix(5, 5); ///< Projection from measurement to local system
private Vector measResiduals = new Vector(5); ///< Measurement residuals
private Vector measPrecision = new Vector(5); ///< Measurement precision (diagonal of inverse covariance matrix)
private boolean transFlag; ///< Transformation exists?
private Matrix measTransformation; ///< Transformation of diagonalization (of meas. precision matrix)
private boolean scatFlag; ///< Scatterer present?
- private SymMatrix scatTransformation = new SymMatrix(2); ///< Transformation of diagonalization (of scat. precision matrix)
+ private Matrix scatTransformation = new Matrix(2, 2); ///< Transformation of diagonalization (of scat. precision matrix)
private Vector scatResiduals = new Vector(2); ///< Scattering residuals (initial kinks if iterating)
private Vector scatPrecision = new Vector(2); ///< Scattering precision (diagonal of inverse covariance matrix)
private Matrix localDerivatives = new Matrix(0, 0); ///< Derivatives of measurement vs additional local (fit) parameters
@@ -110,7 +107,8 @@
* previous point. \param [in] aJacobian Transformation jacobian from
* previous point
*/
- public GblPoint(Matrix aJacobian) {
+ public GblPoint(Matrix aJacobian)
+ {
theLabel = 0;
theOffset = 0;
@@ -129,7 +127,8 @@
}
}
- public GblPoint(SymMatrix aJacobian) {
+ public GblPoint(SymMatrix aJacobian)
+ {
theLabel = 0;
theOffset = 0;
p2pJacobian = new SymMatrix(aJacobian);
@@ -153,8 +152,9 @@
* Minimal precision to accept measurement
*/
public void addMeasurement(Matrix aProjection,
- Vector aResiduals, Vector aPrecision,
- double minPrecision) {
+ Vector aResiduals, Vector aPrecision,
+ double minPrecision)
+ {
measDim = aResiduals.getRowDimension();
int iOff = 5 - measDim;
for (int i = 0; i < measDim; ++i) {
@@ -176,8 +176,9 @@
* (matrix) \param [in] minPrecision Minimal precision to accept measurement
*/
public void addMeasurement(Matrix aProjection,
- Vector aResiduals, SymMatrix aPrecision,
- double minPrecision) {
+ Vector aResiduals, SymMatrix aPrecision,
+ double minPrecision)
+ {
measDim = aResiduals.getRowDimension();
EigenvalueDecomposition measEigen = new EigenvalueDecomposition(aPrecision);
measTransformation.ResizeTo(measDim, measDim);
@@ -206,7 +207,8 @@
* minPrecision Minimal precision to accept measurement
*/
public void addMeasurement(Vector aResiduals,
- Vector aPrecision, double minPrecision) {
+ Vector aPrecision, double minPrecision)
+ {
measDim = aResiduals.getRowDimension();
int iOff = 5 - measDim;
for (int i = 0; i < measDim; ++i) {
@@ -214,7 +216,7 @@
measPrecision.set(iOff + i,
aPrecision.get(i) >= minPrecision ? aPrecision.get(i) : 0.);
}
- measProjection.setToIdentity();
+ measProjection.UnitMatrix();// setToIdentity();
}
/// Add a measurement to a point.
@@ -226,7 +228,8 @@
* (matrix) \param [in] minPrecision Minimal precision to accept measurement
*/
public void addMeasurement(Vector aResiduals,
- SymMatrix aPrecision, double minPrecision) {
+ SymMatrix aPrecision, double minPrecision)
+ {
measDim = aResiduals.getRowDimension();
EigenvalueDecomposition measEigen = new EigenvalueDecomposition(aPrecision);
measTransformation.ResizeTo(measDim, measDim);
@@ -249,7 +252,8 @@
/**
* Get dimension of measurement (0 = none). \return measurement dimension
*/
- int hasMeasurement() {
+ int hasMeasurement()
+ {
return measDim;
}
@@ -259,18 +263,20 @@
* local system \param [out] aResiduals Measurement residuals \param [out]
* aPrecision Measurement precision (diagonal)
*/
- public void getMeasurement(SymMatrix aProjection, Vector aResiduals,
- Vector aPrecision) {
- aProjection = measProjection;
- aResiduals = measResiduals;
- aPrecision = measPrecision;
+ public void getMeasurement(Matrix aProjection, Vector aResiduals,
+ Vector aPrecision)
+ {
+ aProjection.placeAt(measProjection, 0, 0);
+ aResiduals.placeInCol(measResiduals, 0, 0);
+ aPrecision.placeInCol(measPrecision, 0, 0);
}
/// Get measurement transformation (from diagonalization).
/**
* \param [out] aTransformation Transformation matrix
*/
- public void getMeasTransformation(Matrix aTransformation) {
+ public void getMeasTransformation(Matrix aTransformation)
+ {
aTransformation.ResizeTo(measDim, measDim);
if (transFlag) {
aTransformation = measTransformation;
@@ -288,13 +294,14 @@
* Scatterer precision (diagonal of inverse covariance matrix)
*/
public void addScatterer(Vector aResiduals,
- Vector aPrecision) {
+ Vector aPrecision)
+ {
scatFlag = true;
scatResiduals.set(0, aResiduals.get(0));
scatResiduals.set(1, aResiduals.get(1));
scatPrecision.set(0, aPrecision.get(0));
scatPrecision.set(1, aPrecision.get(1));
- scatTransformation.setToIdentity();
+ scatTransformation.UnitMatrix();// setToIdentity();
}
/// Add a (thin) scatterer to a point.
@@ -313,7 +320,8 @@
* Scatterer precision (matrix)
*/
public void addScatterer(Vector aResiduals,
- SymMatrix aPrecision) {
+ SymMatrix aPrecision)
+ {
scatFlag = true;
EigenvalueDecomposition scatEigen = new EigenvalueDecomposition(aPrecision);
Matrix aTransformation = scatEigen.getEigenVectors();
@@ -330,7 +338,8 @@
}
/// Check for scatterer at a point.
- boolean hasScatterer() {
+ boolean hasScatterer()
+ {
return scatFlag;
}
@@ -340,18 +349,20 @@
* diagonalization \param [out] aResiduals Scatterer residuals \param [out]
* aPrecision Scatterer precision (diagonal)
*/
- public void getScatterer(SymMatrix aTransformation, Vector aResiduals,
- Vector aPrecision) {
- aTransformation = scatTransformation;
- aResiduals = scatResiduals;
- aPrecision = scatPrecision;
+ public void getScatterer(Matrix aTransformation, Vector aResiduals,
+ Vector aPrecision)
+ {
+ aTransformation.placeAt(scatTransformation, 0, 0);
+ aResiduals.placeAt(scatResiduals, 0, 0);
+ aPrecision.placeAt(scatPrecision, 0, 0);
}
/// Get scatterer transformation (from diagonalization).
/**
* \param [out] aTransformation Transformation matrix
*/
- public void getScatTransformation(Matrix aTransformation) {
+ public void getScatTransformation(Matrix aTransformation)
+ {
aTransformation.ResizeTo(2, 2);
if (scatFlag) {
for (int i = 0; i < 2; ++i) {
@@ -369,7 +380,8 @@
* Point needs to have a measurement. \param [in] aDerivatives Local
* derivatives (matrix)
*/
- public void addLocals(Matrix aDerivatives) {
+ public void addLocals(Matrix aDerivatives)
+ {
if (measDim != 0) {
localDerivatives.ResizeTo(aDerivatives.getRowDimension(), aDerivatives.getColumnDimension());
if (transFlag) {
@@ -381,12 +393,14 @@
}
/// Retrieve number of local derivatives from a point.
- int getNumLocals() {
+ int getNumLocals()
+ {
return localDerivatives.getColumnDimension();
}
/// Retrieve local derivatives from a point.
- Matrix getLocalDerivatives() {
+ Matrix getLocalDerivatives()
+ {
return localDerivatives;
}
@@ -396,7 +410,8 @@
* labels \param [in] aDerivatives Global derivatives (matrix)
*/
public void addGlobals(List<Integer> aLabels,
- Matrix aDerivatives) {
+ Matrix aDerivatives)
+ {
if (measDim != 0) {
globalLabels = aLabels;
globalDerivatives.ResizeTo(aDerivatives.getRowDimension(), aDerivatives.getColumnDimension());
@@ -410,17 +425,20 @@
}
/// Retrieve number of global derivatives from a point.
- int getNumGlobals() {
+ int getNumGlobals()
+ {
return globalDerivatives.getColumnDimension();
}
/// Retrieve global derivatives labels from a point.
- List<Integer> getGlobalLabels() {
+ List<Integer> getGlobalLabels()
+ {
return globalLabels;
}
/// Retrieve global derivatives from a point.
- Matrix getGlobalDerivatives() {
+ Matrix getGlobalDerivatives()
+ {
return globalDerivatives;
}
@@ -428,12 +446,14 @@
/**
* \param [in] aLabel Label identifying point
*/
- public void setLabel(int aLabel) {
+ public void setLabel(int aLabel)
+ {
theLabel = aLabel;
}
/// Retrieve label of point
- int getLabel() {
+ int getLabel()
+ {
return theLabel;
}
@@ -441,17 +461,20 @@
/**
* \param [in] anOffset Offset number
*/
- public void setOffset(int anOffset) {
+ public void setOffset(int anOffset)
+ {
theOffset = anOffset;
}
/// Retrieve offset for point
- int getOffset() {
+ int getOffset()
+ {
return theOffset;
}
/// Retrieve point-to-(previous)point jacobian
- Matrix getP2pJacobian() {
+ Matrix getP2pJacobian()
+ {
return p2pJacobian;
}
@@ -459,7 +482,8 @@
/**
* \param [in] aJac Jacobian
*/
- public void addPrevJacobian(Matrix aJac) {
+ public void addPrevJacobian(Matrix aJac)
+ {
int ifail = 0;
// to optimize: need only two last rows of inverse
// prevJacobian = aJac.InverseFast(ifail);
@@ -476,7 +500,8 @@
/**
* \param [in] aJac Jacobian
*/
- public void addNextJacobian(Matrix aJac) {
+ public void addNextJacobian(Matrix aJac)
+ {
nextJacobian = aJac;
}
@@ -489,21 +514,23 @@
* std::overflow_error : matrix S is singular.
*/
public void getDerivatives(int aDirection, Matrix matW, Matrix matWJ,
- Vector vecWd) {
+ Vector vecWd)
+ {
Matrix matJ;
+ Matrix matWt;
Vector vecd;
if (aDirection < 1) {
matJ = prevJacobian.sub(2, 3, 3);
- matW = prevJacobian.sub(2, 3, 1).uminus();
+ matWt = prevJacobian.sub(2, 3, 1).uminus();
vecd = prevJacobian.subCol(2, 0, 3);
} else {
matJ = nextJacobian.sub(2, 3, 3);
- matW = nextJacobian.sub(2, 3, 1);
+ matWt = nextJacobian.sub(2, 3, 1);
vecd = nextJacobian.subCol(2, 0, 3);
}
- matW = matW.inverse();
+ matW.placeAt(matWt.inverse(), 0, 0);
// if (!matW.InvertFast()) {
// std::cout << " getDerivatives failed to invert matrix: "
// << matW << "\n";
@@ -512,8 +539,8 @@
// << "\n";
// throw std::overflow_error("Singular matrix inversion exception");
// }
- matWJ = matW.times(matJ);
- vecWd = matW.times(vecd);
+ matWJ.placeAt(matW.times(matJ), 0, 0);
+ vecWd.placeAt(matW.times(vecd), 0, 0);
}
//
@@ -522,7 +549,8 @@
/**
* \param [in] level print level (0: minimum, >0: more)
*/
- public void printPoint(int level) {
+ public void printPoint(int level)
+ {
StringBuffer sb = new StringBuffer("GblPoint");
if (theLabel != 0) {
sb.append(", label " + theLabel);
@@ -586,5 +614,4 @@
}
System.out.println(sb.toString());
}
-
}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectory.java 2014-09-17 17:45:13 UTC (rev 1035)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblTrajectory.java 2014-09-17 18:01:56 UTC (rev 1036)
@@ -1,43 +1,48 @@
package org.hps.recon.tracking.gbl;
/**
- *
- *
+ *
+ *
* @author Per Hansson Adrian <[log in to unmask]>
*
* @author Norman A Graf
*
* @version $Id:
- *
+ *
*/
-
+import static java.lang.Math.abs;
import static java.lang.Math.max;
import java.util.ArrayList;
import java.util.List;
+import org.apache.commons.math3.util.Pair;
import org.hps.recon.tracking.gbl.matrix.BorderedBandMatrix;
import org.hps.recon.tracking.gbl.matrix.Matrix;
import org.hps.recon.tracking.gbl.matrix.SymMatrix;
import org.hps.recon.tracking.gbl.matrix.VVector;
import org.hps.recon.tracking.gbl.matrix.Vector;
+public class GblTrajectory
+{
-public class GblTrajectory {
+ private boolean debug = true;
- public GblTrajectory(List<GblPoint> listOfPoints) {
- this(listOfPoints, 0, new SymMatrix(0), false, false, false);
- }
+ public GblTrajectory(List<GblPoint> listOfPoints)
+ {
+ this(listOfPoints, true, true, true);
+ }
+ public void fit(double m_chi2, int m_ndf, int m_lost_weight)
+ {
+ // TODO Auto-generated method stub
- public void fit(double m_chi2, int m_ndf, int m_lost_weight) {
- // TODO Auto-generated method stub
-
- }
+ }
- public void milleOut(MilleBinary mille) {
- // TODO Auto-generated method stub
-
- }
+ public void milleOut(MilleBinary mille)
+ {
+ // TODO Auto-generated method stub
+ }
+
/**
* \mainpage General information
*
@@ -142,7 +147,7 @@
List<List<GblPoint>> thePoints = new ArrayList<List<GblPoint>>(); ///< (list of) List of points on trajectory
List<GblData> theData = new ArrayList<GblData>(); ///< List of data blocks
List<Integer> measDataIndex = new ArrayList<Integer>(); ///< mapping points to data blocks from measurements
- List<Integer> scatDataIndex; ///< mapping points to data blocks from scatterers
+ List<Integer> scatDataIndex = new ArrayList<Integer>(); ///< mapping points to data blocks from scatterers
List<Integer> externalIndex; ///< List of fit parameters used by external seed
SymMatrix externalSeed; ///< Precision (inverse covariance matrix) of external seed
List<Matrix> innerTransformations; ///< Transformations at innermost points of
@@ -160,38 +165,31 @@
// * [in] flagCurv Use q/p \param [in] flagU1dir Use in u1 direction \param
// * [in] flagU2dir Use in u2 direction
// */
-// GblTrajectory(List<GblPoint> aPointList,
-// boolean flagCurv, boolean flagU1dir, boolean flagU2dir)
-// {
-// numAllPoints = aPointList.size();
-// //numPoints()
-// numOffsets = 0;
-// numInnerTrans = 0;
-// numCurvature = (flagCurv ? 1 : 0);
-// numParameters = 0;
-// numLocals = 0;
-// numMeasurements = 0;
-// externalPoint = 0;
-// //externalSeed()
-// //innerTransformations()
-// //externalDerivatives()
-// //externalMeasurements()
-// // externalPrecisions()
-//
-// if (flagU1dir)
-// {
-// theDimension.add(0);
-// }
-// if (flagU2dir)
-// {
-// theDimension.add(1);
-// }
+ GblTrajectory(List<GblPoint> aPointList,
+ boolean flagCurv, boolean flagU1dir, boolean flagU2dir)
+ {
+ numAllPoints = aPointList.size();
+ numOffsets = 0;
+ numInnerTrans = 0;
+ numCurvature = (flagCurv ? 1 : 0);
+ numParameters = 0;
+ numLocals = 0;
+ numMeasurements = 0;
+ externalPoint = 0;
+
+ if (flagU1dir) {
+ theDimension.add(0);
+ }
+ if (flagU2dir) {
+ theDimension.add(1);
+ }
// // simple (single) trajectory
-// thePoints.add(aPointList);
-// numPoints.add(numAllPoints);
-//// construct(); // construct trajectory
-// }
+ thePoints.add(aPointList);
+ numPoints.add(numAllPoints);
+ construct(); // construct trajectory
+ }
/// Create new (simple) trajectory from list of points with external seed.
+
/**
* Curved trajectory in space (default) or without curvature (q/p) or in one
* plane (u-direction) only. \param [in] aPointList List of points \param
@@ -201,36 +199,21 @@
* flagU1dir Use in u1 direction \param [in] flagU2dir Use in u2 direction
*/
GblTrajectory(List<GblPoint> aPointList,
- int aLabel, SymMatrix aSeed, boolean flagCurv,
- boolean flagU1dir, boolean flagU2dir)
+ int aLabel, SymMatrix aSeed, boolean flagCurv,
+ boolean flagU1dir, boolean flagU2dir)
{
numAllPoints = aPointList.size();
- //numPoints(),
numOffsets = 0;
numInnerTrans = 0;
numCurvature = (flagCurv ? 1 : 0);
numParameters = 0;
numLocals = 0;
numMeasurements = 0;
- //externalPoint(aLabel)
- //theDimension(),
- //thePoints(),
- //theData(),
- //measDataIndex(),
- //scatDataIndex(),
- //externalIndex(),
- //externalSeed(aSeed),
- //innerTransformations(),
- //externalDerivatives(),
- //externalMeasurements(),
- //externalPrecisions() {
- if (flagU1dir)
- {
+ if (flagU1dir) {
theDimension.add(0);
}
- if (flagU2dir)
- {
+ if (flagU2dir) {
theDimension.add(1);
}
// simple (single) trajectory
@@ -239,62 +222,6 @@
construct(); // construct trajectory
}
-///// Create new composed trajectory from list of points and transformations.
-///**
-// * Composed of curved trajectory in space.
-// * \param [in] aPointsAndTransList List containing pairs with list of points and transformation (at inner (first) point)
-// */
-//GblTrajectory(
-// const List<std::pair<List<GblPoint>, TMatrixD> > &aPointsAndTransList) :
-// numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
-// aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
-// 0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalIndex(), externalSeed(), innerTransformations(), externalDerivatives(), externalMeasurements(), externalPrecisions() {
-//
-// for ( int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
-// thePoints.add(aPointsAndTransList[iTraj].first);
-// numPoints.add(thePoints.back().size());
-// numAllPoints += numPoints.back();
-// innerTransformations.add(aPointsAndTransList[iTraj].second);
-// }
-// theDimension.add(0);
-// theDimension.add(1);
-// numCurvature = innerTransformations[0].GetNcols();
-// construct(); // construct (composed) trajectory
-//}
-//
-///// Create new composed trajectory from list of points and transformations with (independent) external measurements.
-///**
-// * Composed of curved trajectory in space.
-// * \param [in] aPointsAndTransList List containing pairs with list of points and transformation (at inner (first) point)
-// * \param [in] extDerivatives Derivatives of external measurements vs external parameters
-// * \param [in] extMeasurements External measurements (residuals)
-// * \param [in] extPrecisions Precision of external measurements
-// */
-//GblTrajectory(
-// const List<std::pair<List<GblPoint>, TMatrixD> > &aPointsAndTransList,
-// const TMatrixD &extDerivatives, const TVectorD &extMeasurements,
-// const TVectorD &extPrecisions) :
-// numAllPoints(), numPoints(), numOffsets(0), numInnerTrans(
-// aPointsAndTransList.size()), numParameters(0), numLocals(0), numMeasurements(
-// 0), externalPoint(0), theDimension(0), thePoints(), theData(), measDataIndex(), scatDataIndex(), externalIndex(), externalSeed(), innerTransformations(), externalDerivatives(
-// extDerivatives), externalMeasurements(extMeasurements), externalPrecisions(
-// extPrecisions) {
-//
-// for ( int iTraj = 0; iTraj < aPointsAndTransList.size(); ++iTraj) {
-// thePoints.add(aPointsAndTransList[iTraj].first);
-// numPoints.add(thePoints.back().size());
-// numAllPoints += numPoints.back();
-// innerTransformations.add(aPointsAndTransList[iTraj].second);
-// }
-// theDimension.add(0);
-// theDimension.add(1);
-// numCurvature = innerTransformations[0].GetNcols();
-// construct(); // construct (composed) trajectory
-//}
-//
-//~GblTrajectory() {
-//}
-//
/// Retrieve validity of trajectory
boolean isValid()
{
@@ -320,20 +247,8 @@
int aLabel = 0;
// loop over trajectories
numTrajectories = thePoints.size();
- //std::cout << " numTrajectories: " << numTrajectories << ", "<< innerTransformations.size() << std::endl;
-// for ( int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
-// List<GblPoint>::iterator itPoint;
-// for (itPoint = thePoints[iTraj].begin();
-// itPoint < thePoints[iTraj].end(); ++itPoint) {
-// numLocals = std::max(numLocals, itPoint->getNumLocals());
-// numMeasurements += itPoint->hasMeasurement();
-// itPoint->setLabel(++aLabel);
-// }
-// }
- for (List<GblPoint> list : thePoints)
- {
- for (GblPoint p : list)
- {
+ for (List<GblPoint> list : thePoints) {
+ for (GblPoint p : list) {
numLocals = max(numLocals, p.getNumLocals());
numMeasurements += p.hasMeasurement();
p.setLabel(++aLabel);
@@ -341,13 +256,7 @@
}
defineOffsets();
calcJacobians();
-// try {
- prepare();
-// } catch (std::overflow_error &e) {
-// std::cout << " GblTrajectory construction failed: " << e.what()
-// << std::endl;
-// return;
-// }
+ prepare();
constructOK = true;
// number of fit parameters
numParameters = (numOffsets - 2 * numInnerTrans) * theDimension.size()
@@ -364,36 +273,20 @@
{
// loop over trajectories
- for (int iTraj = 0; iTraj < numTrajectories; ++iTraj)
- {
+ for (int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
List<GblPoint> list = thePoints.get(iTraj);
int size = list.size();
// first point is offset
-// thePoints[iTraj].front().setOffset(numOffsets++);
list.get(0).setOffset(numOffsets++); // intermediate scatterers are offsets
-// List<GblPoint>::iterator itPoint;
-// for (itPoint = thePoints[iTraj].begin() + 1;
-// itPoint < thePoints[iTraj].end() - 1; ++itPoint) {
-// if (itPoint->hasScatterer()) {
-// itPoint->setOffset(numOffsets++);
-// } else {
-// itPoint->setOffset(-numOffsets);
-// }
-// }
- for (int i = 1; i < size - 1; ++i)
- {
+ for (int i = 1; i < size - 1; ++i) {
GblPoint p = list.get(i);
- if (p.hasScatterer())
- {
+ if (p.hasScatterer()) {
p.setOffset(numOffsets++);
- } else
- {
+ } else {
p.setOffset(-numOffsets);
}
-
}
// last point is offset
-// thePoints[iTraj].back().setOffset(numOffsets++);
list.get(size - 1).setOffset(numOffsets++);
}
}
@@ -404,70 +297,31 @@
Matrix scatJacobian = new Matrix(5, 5);
// loop over trajectories
- for (int iTraj = 0; iTraj < numTrajectories; ++iTraj)
- {
+ for (int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
List<GblPoint> list = thePoints.get(iTraj);
int size = list.size();
// forward propagation (all)
GblPoint previousPoint = list.get(0);
int numStep = 0;
- for (int i = 1; i < size; ++i)
- {
+ for (int i = 1; i < size; ++i) {
GblPoint p = list.get(i);
- if (numStep == 0)
- {
+ if (numStep == 0) {
scatJacobian = p.getP2pJacobian();
- } else
- {
+ } else {
scatJacobian = p.getP2pJacobian().times(scatJacobian);
}
numStep++;
p.addPrevJacobian(scatJacobian);// iPoint -> previous scatterer
- if (p.getOffset() >= 0)
- {
+ if (p.getOffset() >= 0) {
previousPoint.addNextJacobian(scatJacobian); // lastPoint -> next scatterer
numStep = 0;
previousPoint = p;
}
}
-// List<GblPoint>::iterator itPoint;
-// for (itPoint = thePoints[iTraj].begin() + 1;
-// itPoint < thePoints[iTraj].end(); ++itPoint)
-// {
-// if (numStep == 0)
-// {
-// scatJacobian = itPoint -> getP2pJacobian();
-// } else
-// {
-// scatJacobian = itPoint -> getP2pJacobian() * scatJacobian;
-// }
-// numStep++;
-// itPoint -> addPrevJacobian(scatJacobian); // iPoint -> previous scatterer
-// if (itPoint -> getOffset() >= 0)
-// {
-// previousPoint -> addNextJacobian(scatJacobian); // lastPoint -> next scatterer
-// numStep = 0;
-// previousPoint = & ( * itPoint);
-// }
-// }
-// // backward propagation (without scatterers)
-// for (itPoint = thePoints[iTraj].end() - 1;
-// itPoint > thePoints[iTraj].begin(); --itPoint)
-// {
-// if (itPoint -> getOffset() >= 0)
-// {
-// scatJacobian = itPoint -> getP2pJacobian();
-// continue; // skip offsets
-// }
-// itPoint -> addNextJacobian(scatJacobian); // iPoint -> next scatterer
-// scatJacobian = scatJacobian * itPoint -> getP2pJacobian();
-// }
// backward propagation (without scatterers)
- for (int i = size - 1; i > 0; --i)
- {
+ for (int i = size - 1; i > 0; --i) {
GblPoint p = list.get(i);
- if (p.getOffset() >= 0)
- {
+ if (p.getOffset() >= 0) {
scatJacobian = p.getP2pJacobian();
continue; // skip offsets
}
@@ -476,82 +330,85 @@
}
}
}
-//
-///// Get jacobian for transformation from fit to track parameters at point.
-///**
-// * Jacobian broken lines (q/p,..,u_i,u_i+1..) to track (q/p,u',u) parameters
-// * including additional local parameters.
-// * \param [in] aSignedLabel (Signed) label of point for external seed
-// * (<0: in front, >0: after point, slope changes at scatterer!)
-// * \return List of fit parameters with non zero derivatives and
-// * corresponding transformation matrix
-// */
-//std::pair<List< int>, TMatrixD> getJacobian(
-// int aSignedLabel) const {
-//
-// int nDim = theDimension.size();
-// int nCurv = numCurvature;
-// int nLocals = numLocals;
-// int nBorder = nCurv + nLocals;
-// int nParBRL = nBorder + 2 * nDim;
-// int nParLoc = nLocals + 5;
-// List< int> anIndex;
-// anIndex.reserve(nParBRL);
-// TMatrixD aJacobian(nParLoc, nParBRL);
-// aJacobian.Zero();
-//
-// int aLabel = abs(aSignedLabel);
-// int firstLabel = 1;
-// int lastLabel = 0;
-// int aTrajectory = 0;
-// // loop over trajectories
-// for ( int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
-// aTrajectory = iTraj;
-// lastLabel += numPoints[iTraj];
-// if (aLabel <= lastLabel)
-// break;
-// if (iTraj < numTrajectories - 1)
-// firstLabel += numPoints[iTraj];
-// }
-// int nJacobian; // 0: prev, 1: next
-// // check consistency of (index, direction)
-// if (aSignedLabel > 0) {
-// nJacobian = 1;
-// if (aLabel >= lastLabel) {
-// aLabel = lastLabel;
-// nJacobian = 0;
-// }
-// } else {
-// nJacobian = 0;
-// if (aLabel <= firstLabel) {
-// aLabel = firstLabel;
-// nJacobian = 1;
-// }
-// }
-// const GblPoint aPoint = thePoints[aTrajectory][aLabel - firstLabel];
-// List< int> labDer(5);
-// SMatrix55 matDer;
-// getFitToLocalJacobian(labDer, matDer, aPoint, 5, nJacobian);
-//
-// // from local parameters
-// for ( int i = 0; i < nLocals; ++i) {
-// aJacobian(i + 5, i) = 1.0;
-// anIndex.add(i + 1);
-// }
-// // from trajectory parameters
-// int iCol = nLocals;
-// for ( int i = 0; i < 5; ++i) {
-// if (labDer[i] > 0) {
-// anIndex.add(labDer[i]);
-// for ( int j = 0; j < 5; ++j) {
-// aJacobian(j, iCol) = matDer(j, i);
-// }
-// ++iCol;
-// }
-// }
-// return std::make_pair(anIndex, aJacobian);
-//}
+/// Get jacobian for transformation from fit to track parameters at point.
+ /**
+ * Jacobian broken lines (q/p,..,u_i,u_i+1..) to track (q/p,u',u) parameters
+ * including additional local parameters. \param [in] aSignedLabel (Signed)
+ * label of point for external seed (<0: in front, >0: after point, slope
+ * changes at scatterer!) \return List of fit parameters with non zero
+ * derivatives and corresponding transformation matrix
+ */
+ Pair<List< Integer>, Matrix> getJacobian(int aSignedLabel)
+ {
+
+ int nDim = theDimension.size();
+ int nCurv = numCurvature;
+ int nLocals = numLocals;
+ int nBorder = nCurv + nLocals;
+ int nParBRL = nBorder + 2 * nDim;
+ int nParLoc = nLocals + 5;
+ List< Integer> anIndex = new ArrayList<Integer>();
+
+ Matrix aJacobian = new Matrix(nParLoc, nParBRL);
+
+ int aLabel = abs(aSignedLabel);
+ int firstLabel = 1;
+ int lastLabel = 0;
+ int aTrajectory = 0;
+ // loop over trajectories
+ for (int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
+ aTrajectory = iTraj;
+ lastLabel += numPoints.get(iTraj);
+ if (aLabel <= lastLabel) {
+ break;
+ }
+ if (iTraj < numTrajectories - 1) {
+ firstLabel += numPoints.get(iTraj);
+ }
+ }
+ int nJacobian; // 0: prev, 1: next
+ // check consistency of (index, direction)
+ if (aSignedLabel > 0) {
+ nJacobian = 1;
+ if (aLabel >= lastLabel) {
+ aLabel = lastLabel;
+ nJacobian = 0;
+ }
+ } else {
+ nJacobian = 0;
+ if (aLabel <= firstLabel) {
+ aLabel = firstLabel;
+ nJacobian = 1;
+ }
+ }
+ GblPoint aPoint = thePoints.get(aTrajectory).get(aLabel - firstLabel);
+ List< Integer> labDer = new ArrayList<Integer>();
+ for (int i = 0; i < 5; ++i) {
+ labDer.add(0);
+ }
+ Matrix matDer = new Matrix(5, 5);
+ getFitToLocalJacobian(labDer, matDer, aPoint, 5, nJacobian);
+
+ // from local parameters
+ for (int i = 0; i < nLocals; ++i) {
+ aJacobian.set(i + 5, i, 1.0);
+ anIndex.add(i + 1);
+ }
+ // from trajectory parameters
+ int iCol = nLocals;
+ for (int i = 0; i < 5; ++i) {
+ if (labDer.get(i) > 0) {
+ anIndex.add(labDer.get(i));
+ for (int j = 0; j < 5; ++j) {
+ aJacobian.set(j, iCol, matDer.get(j, i));
+ }
+ ++iCol;
+ }
+ }
+ return new Pair(anIndex, aJacobian);
+ }
+
/// Get (part of) jacobian for transformation from (trajectory) fit to track parameters at point.
/**
* Jacobian broken lines (q/p,..,u_i,u_i+1..) to local (q/p,u',u)
@@ -563,8 +420,8 @@
* offset)
*/
void getFitToLocalJacobian(List<Integer> anIndex,
- Matrix aJacobian, GblPoint aPoint, int measDim,
- int nJacobian)
+ Matrix aJacobian, GblPoint aPoint, int measDim,
+ int nJacobian)
{
int nDim = theDimension.size();
@@ -596,35 +453,18 @@
int iOff = nDim * (-nOffset - 1) + nLocals + nCurv + 1; // first offset ('i' in u_i)
// local offset
- if (nCurv > 0)
- {
+ if (nCurv > 0) {
Vector negDiff = prevNd.minus(nextNd).uminus();
aJacobian.placeInCol(negDiff, 3, 0); // from curvature
anIndex.set(0, nLocals + 1);
}
aJacobian.placeAt(prevNW, 3, 1); // from 1st Offset
aJacobian.placeAt(nextNW, 3, 3); // from 2nd Offset
- for (int i = 0; i < nDim; ++i)
- {
+ for (int i = 0; i < nDim; ++i) {
anIndex.set(1 + theDimension.get(i), iOff + i);
anIndex.set(3 + theDimension.get(i), iOff + nDim + i);
}
-//
-//// // local slope and curvature
-//// if (measDim > 2) {
-//// // derivatives for u'_int
-//// const SMatrix22 prevWPN(nextWJ * prevNW); // W+ * J+ * N * W-
-//// const SMatrix22 nextWPN(prevWJ * nextNW); // W- * J- * N * W+
-//// const SVector2 prevWNd(nextWJ * prevNd); // W+ * J+ * N * W- * d-
-//// const SVector2 nextWNd(prevWJ * nextNd); // W- * J- * N * W+ * d+
-//// if (nCurv > 0) {
-//// aJacobian(0, 0) = 1.0;
-//// aJacobian.Place_in_col(prevWNd - nextWNd, 1, 0); // from curvature
-//// }
-//// aJacobian.Place_at(-prevWPN, 1, 1); // from 1st Offset
-//// aJacobian.Place_at(nextWPN, 1, 3); // from 2nd Offset
-//// }
-// } else { // at point
+ } else { // at point
// anIndex must be sorted
// forward : iOff2 = iOff1 + nDim, index1 = 1, index2 = 3
// backward: iOff2 = iOff1 - nDim, index1 = 3, index2 = 1
@@ -635,29 +475,25 @@
// local offset
aJacobian.set(3, index1, 1.0); // from 1st Offset
aJacobian.set(4, index1 + 1, 1.0);
- for (int i = 0; i < nDim; ++i)
- {
+ for (int i = 0; i < nDim; ++i) {
anIndex.set(index1 + theDimension.get(i), iOff1 + i);
}
// local slope and curvature
- if (measDim > 2)
- {
+ if (measDim > 2) {
Matrix matW = new Matrix(2, 2);
Matrix matWJ = new Matrix(2, 2);
Vector vecWd = new Vector(2);
aPoint.getDerivatives(nJacobian, matW, matWJ, vecWd); // W, W * J, W * d
double sign = (nJacobian > 0) ? 1. : -1.;
- if (nCurv > 0)
- {
+ if (nCurv > 0) {
aJacobian.set(0, 0, 1.0);
aJacobian.placeInCol(vecWd.timesScalar(-sign), 1, 0); // from curvature
anIndex.set(0, nLocals + 1);
}
aJacobian.placeAt(matWJ.times(-sign), 1, index1); // from 1st Offset
aJacobian.placeAt(matW.times(sign), 1, index2); // from 2nd Offset
- for (int i = 0; i < nDim; ++i)
- {
+ for (int i = 0; i < nDim; ++i) {
anIndex.set(index2 + theDimension.get(i), iOff2 + i);
}
}
@@ -672,7 +508,7 @@
* \param [in] aPoint Point to use
*/
void getFitToKinkJacobian(List< Integer> anIndex,
- Matrix aJacobian, GblPoint aPoint)
+ Matrix aJacobian, GblPoint aPoint)
{
//nb aJacobian has dimension 2,7
@@ -696,431 +532,187 @@
int iOff = (nOffset - 1) * nDim + nCurv + nLocals + 1; // first offset ('i' in u_i)
// local offset
- if (nCurv > 0)
- {
+ if (nCurv > 0) {
aJacobian.placeInCol(sumWd.uminus(), 0, 0); // from curvature
anIndex.set(0, nLocals + 1);
}
aJacobian.placeAt(prevW, 0, 1); // from 1st Offset
aJacobian.placeAt(sumWJ.uminus(), 0, 3); // from 2nd Offset
aJacobian.placeAt(nextW, 0, 5); // from 1st Offset
- for (int i = 0; i < nDim; ++i)
- {
+ for (int i = 0; i < nDim; ++i) {
anIndex.set(1 + theDimension.get(i), iOff + i);
anIndex.set(3 + theDimension.get(i), iOff + nDim + i);
anIndex.set(5 + theDimension.get(i), iOff + nDim * 2 + i);
}
}
-///// Get fit results at point.
-///**
-// * Get corrections and covariance matrix for local track and additional parameters
-// * in forward or backward direction.
-// * \param [in] aSignedLabel (Signed) label of point on trajectory
-// * (<0: in front, >0: after point, slope changes at scatterer!)
-// * \param [out] localPar Corrections for local parameters
-// * \param [out] localCov Covariance for local parameters
-// * \return error code (non-zero if trajectory not fitted successfully)
-// */
-// int getResults(int aSignedLabel, TVectorD &localPar,
-// SymMatrix &localCov) const {
-// if (! fitOK)
-// return 1;
-// std::pair<List< int>, TMatrixD> indexAndJacobian =
-// getJacobian(aSignedLabel);
-// int nParBrl = indexAndJacobian.first.size();
-// TVectorD aVec(nParBrl); // compressed vector
-// for ( int i = 0; i < nParBrl; ++i) {
-// aVec[i] = theVector(indexAndJacobian.first[i] - 1);
-// }
-// SymMatrix aMat = theMatrix.getBlockMatrix(indexAndJacobian.first); // compressed matrix
-// localPar = indexAndJacobian.second * aVec;
-// localCov = aMat.Similarity(indexAndJacobian.second);
-// return 0;
-//}
-//
-///// Get residuals at point from measurement.
-///**
-// * Get (diagonalized) residual, error of measurement and residual and down-weighting
-// * factor for measurement at point
-// *
-// * \param [in] aLabel Label of point on trajectory
-// * \param [out] numData Number of data blocks from measurement at point
-// * \param [out] aResiduals Measurements-Predictions
-// * \param [out] aMeasErrors Errors of Measurements
-// * \param [out] aResErrors Errors of Residuals (including correlations from track fit)
-// * \param [out] aDownWeights Down-Weighting factors
-// * \return error code (non-zero if trajectory not fitted successfully)
-// */
-// int getMeasResults( int aLabel,
-// int &numData, TVectorD &aResiduals, TVectorD &aMeasErrors,
-// TVectorD &aResErrors, TVectorD &aDownWeights) {
-// numData = 0;
-// if (! fitOK)
-// return 1;
-//
-// int firstData = measDataIndex[aLabel - 1]; // first data block with measurement
-// numData = measDataIndex[aLabel] - firstData; // number of data blocks
-// for ( int i = 0; i < numData; ++i) {
-// getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
-// aResErrors[i], aDownWeights[i]);
-// }
-// return 0;
-//}
-//
-///// Get (kink) residuals at point from scatterer.
-///**
-// * Get (diagonalized) residual, error of measurement and residual and down-weighting
-// * factor for scatterering kinks at point
-// *
-// * \param [in] aLabel Label of point on trajectory
-// * \param [out] numData Number of data blocks from scatterer at point
-// * \param [out] aResiduals (kink)Measurements-(kink)Predictions
-// * \param [out] aMeasErrors Errors of (kink)Measurements
-// * \param [out] aResErrors Errors of Residuals (including correlations from track fit)
-// * \param [out] aDownWeights Down-Weighting factors
-// * \return error code (non-zero if trajectory not fitted successfully)
-// */
-// int getScatResults( int aLabel,
-// int &numData, TVectorD &aResiduals, TVectorD &aMeasErrors,
-// TVectorD &aResErrors, TVectorD &aDownWeights) {
-// numData = 0;
-// if (! fitOK)
-// return 1;
-//
-// int firstData = scatDataIndex[aLabel - 1]; // first data block with scatterer
-// numData = scatDataIndex[aLabel] - firstData; // number of data blocks
-// for ( int i = 0; i < numData; ++i) {
-// getResAndErr(firstData + i, aResiduals[i], aMeasErrors[i],
-// aResErrors[i], aDownWeights[i]);
-// }
-// return 0;
-//}
-//
-///// Get (list of) labels of points on (simple) trajectory
-///**
-// * \param [out] aLabelList List of labels (aLabelList[i] = i+1)
-// */
-//void getLabels(List< int> &aLabelList) {
-// int aLabel = 0;
-// int nPoint = thePoints[0].size();
-// aLabelList.resize(nPoint);
-// for ( i = 0; i < nPoint; ++i) {
-// aLabelList[i] = ++aLabel;
-// }
-//}
-//
-///// Get (list of lists of) labels of points on (composed) trajectory
-///**
-// * \param [out] aLabelList List of of lists of labels
-// */
-//void getLabels(
-// List<List< int> > &aLabelList) {
-// int aLabel = 0;
-// aLabelList.resize(numTrajectories);
-// for ( int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
-// int nPoint = thePoints[iTraj].size();
-// aLabelList[iTraj].resize(nPoint);
-// for ( i = 0; i < nPoint; ++i) {
-// aLabelList[iTraj][i] = ++aLabel;
-// }
-// }
-//}
-//
-///// Get residual and errors from data block.
-///**
-// * Get residual, error of measurement and residual and down-weighting
-// * factor for (single) data block
-// * \param [in] aData Label of data block
-// * \param [out] aResidual Measurement-Prediction
-// * \param [out] aMeasError Error of Measurement
-// * \param [out] aResError Error of Residual (including correlations from track fit)
-// * \param [out] aDownWeight Down-Weighting factor
-// */
-//void getResAndErr( int aData, double &aResidual,
-// double &aMeasError, double &aResError, double &aDownWeight) {
-//
-// double aMeasVar;
-// List< int>* indLocal;
-// List<double>* derLocal;
-// theData[aData].getResidual(aResidual, aMeasVar, aDownWeight, indLocal,
-// derLocal);
-// int nParBrl = (*indLocal).size();
-// TVectorD aVec(nParBrl); // compressed vector of derivatives
-// for ( int j = 0; j < nParBrl; ++j) {
-// aVec[j] = (*derLocal)[j];
-// }
-// SymMatrix aMat = theMatrix.getBlockMatrix(*indLocal); // compressed (covariance) matrix
-// double aFitVar = aMat.Similarity(aVec); // variance from track fit
-// aMeasError = sqrt(aMeasVar); // error of measurement
-// aResError = (aFitVar < aMeasVar ? sqrt(aMeasVar - aFitVar) : 0.); // error of residual
-//}
+/// Get fit results at point.
+ /**
+ * Get corrections and covariance matrix for local track and additional
+ * parameters in forward or backward direction. \param [in] aSignedLabel
+ * (Signed) label of point on trajectory (<0: in front, >0: after point,
+ * slope changes at scatterer!) \param [out] localPar Corrections for local
+ * parameters \param [out] localCov Covariance for local parameters \return
+ * error code (non-zero if trajectory not fitted successfully)
+ */
+ int getResults(int aSignedLabel, Vector localPar,
+ SymMatrix localCov)
+ {
+ if (!fitOK) {
+ return 1;
+ }
+ Pair<List< Integer>, Matrix> indexAndJacobian = getJacobian(aSignedLabel);
+ int nParBrl = indexAndJacobian.getFirst().size();
+ Vector aVec = new Vector(nParBrl); // compressed vector
+ for (int i = 0; i < nParBrl; ++i) {
+ aVec.set(i, theVector.get(indexAndJacobian.getFirst().get(i) - 1));
+ }
+ SymMatrix aMat = theMatrix.getBlockMatrix(indexAndJacobian.getFirst()); // compressed matrix
+ localPar.placeAt(indexAndJacobian.getSecond().times(aVec), 0, 0);
+ localCov.placeAt(aMat.Similarity(indexAndJacobian.getSecond()), 0, 0);
+ return 0;
+ }
/// Build linear equation system from data (blocks).
+
void buildLinearEquationSystem()
{
int nBorder = numCurvature + numLocals;
-// theVector. resize(numParameters);
+ theVector = new VVector(numParameters);
theMatrix.resize(numParameters, nBorder, 5);
double[] retVals = new double[2];
double aValue, aWeight;
- int[] indLocal = null;
- double[] derLocal = null;
- for (GblData d : theData)
- {
+ int nData = 0;
+ for (GblData d : theData) {
+ int size = d.getNumParameters();
+ int[] indLocal = new int[size];
+ double[] derLocal = new double[size];
d.getLocalData(retVals, indLocal, derLocal);
aValue = retVals[0];
aWeight = retVals[1];
- for (int j = 0; j < indLocal.length; ++j)
- {
+ for (int j = 0; j < size; ++j) {
theVector.addTo(indLocal[j] - 1, derLocal[j] * aWeight * aValue);
}
theMatrix.addBlockMatrix(aWeight, indLocal, derLocal);
+ nData++;
}
}
//}
/// Prepare fit for simple or composed trajectory
-/**
- * Generate data (blocks) from measurements, kinks, external seed and measurements.
- */
-void prepare() {
- int nDim = theDimension.size();
- // upper limit
- int maxData = numMeasurements + nDim * (numOffsets - 2);
-//cng + externalSeed.getRowDimension();
-// theData.reserve(maxData);
-// measDataIndex.resize(numAllPoints + 3); // include external seed and measurements
- //cng
- for(int i = 0; i<numAllPoints + 3; ++i) measDataIndex.add(i,0);
- //cng
-// scatDataIndex.resize(numAllPoints + 1);
- int nData = 0;
- List<Matrix> innerTransDer = new ArrayList<Matrix>();
- List<List<Integer> > innerTransLab = new ArrayList<List<Integer>>();
-// // composed trajectory ?
-// if (numInnerTrans > 0) {
-// //std::cout << "composed trajectory" << std::endl;
-// for ( int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
-// // innermost point
-// GblPoint* innerPoint = &thePoints[iTraj].front();
-// // transformation fit to local track parameters
-// List< int> firstLabels(5);
-// SMatrix55 matFitToLocal, matLocalToFit;
-// getFitToLocalJacobian(firstLabels, matFitToLocal, *innerPoint, 5);
-// // transformation local track to fit parameters
-// int ierr;
-// matLocalToFit = matFitToLocal.Inverse(ierr);
-// TMatrixD localToFit(5, 5);
-// for ( int i = 0; i < 5; ++i) {
-// for ( int j = 0; j < 5; ++j) {
-// localToFit(i, j) = matLocalToFit(i, j);
-// }
-// }
-// // transformation external to fit parameters at inner (first) point
-// innerTransDer.add(localToFit * innerTransformations[iTraj]);
-// innerTransLab.add(firstLabels);
-// }
-// }
- // measurements
- Matrix matP = new Matrix(5,5);
-// // loop over trajectories
-// List<GblPoint>::iterator itPoint;
-// for ( int iTraj = 0; iTraj < numTrajectories; ++iTraj) {
-// for (itPoint = thePoints[iTraj].begin();
-// itPoint < thePoints[iTraj].end(); ++itPoint) {
-// SVector5 aMeas, aPrec;
-// int nLabel = itPoint->getLabel();
-// int measDim = itPoint->hasMeasurement();
-// if (measDim) {
-// const TMatrixD localDer = itPoint->getLocalDerivatives();
-// const List<int> globalLab = itPoint->getGlobalLabels();
-// const TMatrixD globalDer = itPoint->getGlobalDerivatives();
-// TMatrixD transDer;
-// itPoint->getMeasurement(matP, aMeas, aPrec);
-// int iOff = 5 - measDim; // first active component
-// List< int> labDer(5);
-// SMatrix55 matDer, matPDer;
-// int nJacobian =
-// (itPoint < thePoints[iTraj].end() - 1) ? 1 : 0; // last point needs backward propagation
-// getFitToLocalJacobian(labDer, matDer, *itPoint, measDim,
-// nJacobian);
-// if (measDim > 2) {
-// matPDer = matP * matDer;
-// } else { // 'shortcut' for position measurements
-// matPDer.Place_at(
-// matP.Sub<SMatrix22>(3, 3)
-// * matDer.Sub<SMatrix25>(3, 0), 3, 0);
-// }
-//
-//// if (numInnerTrans > 0) {
-//// // transform for external parameters
-//// TMatrixD proDer(measDim, 5);
-//// // match parameters
-//// int ifirst = 0;
-//// int ilabel = 0;
-//// while (ilabel < 5) {
-//// if (labDer[ilabel] > 0) {
-//// while (innerTransLab[iTraj][ifirst]
-//// != labDer[ilabel] && ifirst < 5) {
-//// ++ifirst;
-//// }
-//// if (ifirst >= 5) {
-//// labDer[ilabel] -= 2 * nDim * (iTraj + 1); // adjust label
-//// } else {
-//// // match
-//// labDer[ilabel] = 0; // mark as related to external parameters
-//// for ( int k = iOff; k < 5; ++k) {
-//// proDer(k - iOff, ifirst) = matPDer(k,
-//// ilabel);
-//// }
-//// }
-//// }
-//// ++ilabel;
-//// }
-//// transDer.ResizeTo(measDim, numCurvature);
-//// transDer = proDer * innerTransDer[iTraj];
-//// }
-// for ( int i = iOff; i < 5; ++i) {
-// if (aPrec(i) > 0.) {
-// GblData aData(nLabel, aMeas(i), aPrec(i));
-// aData.addDerivatives(i, labDer, matPDer, iOff, localDer,
-// globalLab, globalDer, numLocals, transDer);
-// theData.add(aData);
-// nData++;
-// }
-// }
-//
-// }
-// measDataIndex[nLabel] = nData;
-// }
-// } // end loop over trajectories
+ /**
+ * Generate data (blocks) from measurements, kinks, external seed and
+ * measurements.
+ */
+ void prepare()
+ {
+ int nDim = theDimension.size();
+ // upper limit
[truncated at 1000 lines; 484 more skipped]
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java (rev 0)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java 2014-09-17 18:01:56 UTC (rev 1036)
@@ -0,0 +1,418 @@
+package org.hps.recon.tracking.gbl;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.GenericObject;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.Track;
+import org.lcsim.geometry.Detector;
+import org.lcsim.util.Driver;
+
+import static java.lang.Math.sin;
+import static java.lang.Math.sqrt;
+import static java.lang.Math.abs;
+
+import org.hps.recon.tracking.gbl.matrix.Matrix;
+import org.hps.recon.tracking.gbl.matrix.SymMatrix;
+import org.hps.recon.tracking.gbl.matrix.Vector;
+
+/**
+ * A Driver which refits tracks using GBL
+ * Modeled on the hps-dst code written by Per Hansson and Omar Moreno
+ * Requires the GBL Collections and Relations to be present in the event.
+ *
+ * @author Norman A Graf
+ *
+ * @version $Id:
+ */
+public class HpsGblRefitter extends Driver
+{
+ private boolean _debug = false;
+ private final String trackCollectionName = "MatchedTracks";
+ private final String track2GblTrackRelationName = "TrackToGBLTrack";
+ private final String gblTrack2StripRelationName = "GBLTrackToStripData";
+
+ public void setDebug(boolean debug)
+ {
+ _debug = debug;
+ }
+
+ @Override
+ protected void process(EventHeader event)
+ {
+ Hep3Vector bfield = event.getDetector().getFieldMap().getField(new BasicHep3Vector(0., 0., 1.));
+ double By = bfield.y();
+ double bfac = 0.0002998 * By;
+ // get the tracks
+// List<Track> tracks = null;
+ if (event.hasCollection(Track.class, trackCollectionName)) {
+// tracks = event.get(Track.class, trackCollectionName);
+ if (_debug) {
+// System.out.printf("%s: Event %d has %d tracks\n", this.getClass().getSimpleName(), event.getEventNumber(), tracks.size());
+ }
+ } else {
+ if (_debug) {
+ System.out.printf("%s: No tracks in Event %d \n", this.getClass().getSimpleName(), event.getEventNumber());
+ }
+ return;
+ }
+// System.out.println("Tracks from event " + event.getRunNumber());
+// for (Track t : tracks) {
+// System.out.println(t);
+// System.out.println(t.getTrackStates().get(0));
+// }
+ //get the relations to the GBLtracks
+ if(!event.hasItem(track2GblTrackRelationName))
+ {
+ System.out.println("Need Relations "+track2GblTrackRelationName);
+ return;
+ }
+ // and strips
+ if(!event.hasItem(gblTrack2StripRelationName))
+ {
+ System.out.println("Need Relations "+gblTrack2StripRelationName);
+ return;
+ }
+
+
+ List<LCRelation> track2GblTrackRelations = event.get(LCRelation.class, track2GblTrackRelationName);
+ //need a map of GBLTrackData keyed on the Generic object from which it created
+ Map<GenericObject, GBLTrackData> gblObjMap = new HashMap<GenericObject, GBLTrackData>();
+
+ // loop over the relations
+ for (LCRelation relation : track2GblTrackRelations) {
+ Track t = (Track) relation.getFrom();
+ GenericObject gblTrackObject = (GenericObject) relation.getTo();
+ GBLTrackData gblT = new GBLTrackData(gblTrackObject);
+ gblObjMap.put(gblTrackObject, gblT);
+ } //end of loop over tracks
+
+ //get the strip hit relations
+ List<LCRelation> gblTrack2StripRelations = event.get(LCRelation.class, gblTrack2StripRelationName);
+
+ // need a map of lists of strip data keyed by the gblTrack to which they correspond
+ Map<GBLTrackData, List<GBLStripClusterData>> stripsGblMap = new HashMap<GBLTrackData, List<GBLStripClusterData>>();
+ for (LCRelation relation : gblTrack2StripRelations) {
+ //from GBLTrackData to GBLStripClusterData
+ GenericObject gblTrackObject = (GenericObject) relation.getFrom();
+ //Let's get the GBLTrackData that corresponds to this object...
+ GBLTrackData gblT = gblObjMap.get(gblTrackObject);
+ GBLStripClusterData sd = new GBLStripClusterData((GenericObject) relation.getTo());
+ if (stripsGblMap.containsKey(gblT)) {
+ stripsGblMap.get(gblT).add(sd);
+ } else {
+ stripsGblMap.put(gblT, new ArrayList<GBLStripClusterData>());
+ stripsGblMap.get(gblT).add(sd);
+ }
+ }
+
+ Set<GBLTrackData> keys = stripsGblMap.keySet();
+ int trackNum = 0;
+ for (GBLTrackData t : keys) {
+ int stat = fit(t, stripsGblMap.get(t), bfac);
+ ++trackNum;
+ }
+
+ }
+
+ private int fit(GBLTrackData track, List<GBLStripClusterData> hits, double bfac)
+ {
+ // path length along trajectory
+ double s = 0.;
+
+ // jacobian to transport errors between points along the path
+ Matrix jacPointToPoint = new Matrix(5, 5);
+ jacPointToPoint.UnitMatrix();
+ // Option to use uncorrelated MS errors
+ // This is similar to what is done in lcsim seedtracker
+ // The msCov below holds the MS errors
+ // This is for testing purposes only.
+ boolean useUncorrMS = false;
+ Matrix msCov = new Matrix(5, 5);
+ Matrix measMsCov = new Matrix(2, 2);
+ // Vector of the strip clusters used for the GBL fit
+ List<GblPoint> listOfPoints = new ArrayList<GblPoint>();
+
+ // Store the projection from local to measurement frame for each strip cluster
+ Map< Integer, Matrix> proL2m_list = new HashMap<Integer, Matrix>();
+ // Save the association between strip cluster and label
+
+ //start trajectory at refence point (s=0) - this point has no measurement
+ GblPoint ref_point = new GblPoint(jacPointToPoint);
+ listOfPoints.add(ref_point);
+
+ // Loop over strips
+ int n_strips = hits.size();
+ for (int istrip = 0; istrip != n_strips; ++istrip) {
+ GBLStripClusterData strip = hits.get(istrip);
+ if (_debug) {
+ System.out.println("HpsGblFitter: Processing strip " + istrip + " with id/layer " + strip.getId());
+ }
+ // Path length step for this cluster
+ double step = strip.getPath3D() - s;
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + "Path length step " + step + " from " + s + " to " + strip.getPath());
+ }
+ // Measurement direction (perpendicular and parallel to strip direction)
+ Matrix mDir = new Matrix(2, 3);
+ mDir.set(0, 0, strip.getUx());
+ mDir.set(0, 1, strip.getUy());
+ mDir.set(0, 2, strip.getUz());
+ mDir.set(1, 0, strip.getVx());
+ mDir.set(1, 1, strip.getVy());
+ mDir.set(1, 2, strip.getVz());
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + "mDir");
+ mDir.print(4, 6);
+ }
+ Matrix mDirT = mDir.copy().transpose();
+
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + "mDirT");
+ mDirT.print(4, 6);
+ }
+
+ // Track direction
+ double sinLambda = sin(strip.getTrackLambda());//->GetLambda());
+ double cosLambda = sqrt(1.0 - sinLambda * sinLambda);
+ double sinPhi = sin(strip.getTrackPhi());//->GetPhi());
+ double cosPhi = sqrt(1.0 - sinPhi * sinPhi);
+
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + "Track direction sinLambda=" + sinLambda + " sinPhi=" + sinPhi);
+ }
+
+ // Track direction in curvilinear frame (U,V,T)
+ // U = Z x T / |Z x T|, V = T x U
+ Matrix uvDir = new Matrix(2, 3);
+ uvDir.set(0, 0, -sinPhi);
+ uvDir.set(0, 1, cosPhi);
+ uvDir.set(0, 2, 0.);
+ uvDir.set(1, 0, -sinLambda * cosPhi);
+ uvDir.set(1, 1, -sinLambda * sinPhi);
+ uvDir.set(1, 2, cosLambda);
+
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + "uvDir");
+ uvDir.print(6, 4);
+ }
+
+ // projection from measurement to local (curvilinear uv) directions (duv/dm)
+ Matrix proM2l = uvDir.times(mDirT);
+
+ //projection from local (uv) to measurement directions (dm/duv)
+ Matrix proL2m = proM2l.copy();
+ proL2m = proL2m.inverse();
+ proL2m_list.put(strip.getId(), proL2m.copy()); // is a copy needed or is that just a C++/root thing?
+
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + "proM2l:");
+ proM2l.print(4, 6);
+ System.out.println("HpsGblFitter: " + "proL2m:");
+ proL2m.print(4, 6);
+ System.out.println("HpsGblFitter: " + "proM2l*proL2m (should be unit matrix):");
+ (proM2l.times(proL2m)).print(4, 6);
+ }
+
+ // measurement/residual in the measurement system
+ // only 1D measurement in u-direction, set strip measurement direction to zero
+ Vector meas = new Vector(2);
+// double uRes = strip->GetUmeas() - strip->GetTrackPos().x(); // how can this be correct?
+ double uRes = strip.getMeas() - strip.getTrackPos().x();
+ meas.set(0, uRes);
+ meas.set(1, 0.);
+// //meas[0][0] += deltaU[iLayer] # misalignment
+ Vector measErr = new Vector(2);
+ measErr.set(0, strip.getMeasErr());
+ measErr.set(1, 0.);
+ Vector measPrec = new Vector(2);
+ measPrec.set(0, 1.0 / (measErr.get(0) * measErr.get(0)));
+ measPrec.set(1, 0.);
+
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + "meas: ");
+ meas.print(4, 6);
+ System.out.println("HpsGblFitter: " + "measErr:");
+ measErr.print(4, 6);
+ System.out.println("HpsGblFitter: " + "measPrec:");
+ measPrec.print(4, 6);
+ }
+
+ //Find the Jacobian to be able to propagate the covariance matrix to this strip position
+ jacPointToPoint = gblSimpleJacobianLambdaPhi(step, cosLambda, abs(bfac));
+
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + "jacPointToPoint to extrapolate to this point:");
+ jacPointToPoint.print(4, 6);
+ }
+
+ // Get the transpose of the Jacobian
+ Matrix jacPointToPointTransposed = jacPointToPoint.copy().transpose();
+
+ // Propagate the MS covariance matrix (in the curvilinear frame) to this strip position
+ msCov = msCov.times(jacPointToPointTransposed);
+ msCov = jacPointToPoint.times(msCov);
+
+ // Get the MS covariance for the measurements in the measurement frame
+ Matrix proL2mTransposed = proL2m.copy().transpose();
+
+ measMsCov = proL2m.times((msCov.getMatrix(3, 4, 3, 4)).times(proL2mTransposed));
+
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + " msCov at this point:");
+ msCov.print(4, 6);
+ System.out.println("HpsGblFitter: " + "measMsCov at this point:");
+ measMsCov.print(4, 6);
+ }
+
+ GblPoint point = new GblPoint(jacPointToPoint);
+ //Add measurement to the point
+ point.addMeasurement(proL2m, meas, measPrec, 0.);
+ //Add scatterer in curvilinear frame to the point
+ // no direction in this frame
+ Vector scat = new Vector(2);
+
+ // Scattering angle in the curvilinear frame
+ //Note the cosLambda to correct for the projection in the phi direction
+ Vector scatErr = new Vector(2);
+ scatErr.set(0, strip.getScatterAngle());
+ scatErr.set(1, strip.getScatterAngle() / cosLambda);
+ Vector scatPrec = new Vector(2);
+ scatPrec.set(0, 1.0 / (scatErr.get(0) * scatErr.get(0)));
+ scatPrec.set(1, 1.0 / (scatErr.get(1) * scatErr.get(1)));
+
+ // add scatterer if not using the uncorrelated MS covariances for testing
+ if (!useUncorrMS) {
+ point.addScatterer(scat, scatPrec);
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + "adding scatError to this point:");
+ scatErr.print(4, 6);
+ }
+ }
+
+ // Add this GBL point to list that will be used in fit
+ listOfPoints.add(point);
+ int iLabel = listOfPoints.size();
+
+ // Update MS covariance matrix
+ msCov.set(1, 1, msCov.get(1, 1) + scatErr.get(0) * scatErr.get(0));
+ msCov.set(2, 2, msCov.get(2, 2) + scatErr.get(1) * scatErr.get(1));
+
+ /*
+
+ #####
+ ## Calculate global derivatives for this point
+ # track direction in tracking/global frame
+ tDirGlobal = np.array( [ [cosPhi * cosLambda, sinPhi * cosLambda, sinLambda] ] )
+ # Cross-check that the input is consistent
+ if( np.linalg.norm( tDirGlobal - strip.tDir) > 0.00001):
+ print 'ERROR: tDirs are not consistent!'
+ sys.exit(1)
+ # rotate track direction to measurement frame
+ tDirMeas = np.dot( tDirGlobal, np.array([strip.u, strip.v, strip.w]) )
+ #tDirMeas = utils.rotateGlToMeas(strip,tDirGlobal)
+ normalMeas = np.dot( strip.w , np.array([strip.u, strip.v, strip.w]) )
+ #normalMeas = utils.rotateGlToMeas(strip,strip.w)
+ # non-measured directions
+ vmeas = 0.
+ wmeas = 0.
+ # calculate and add derivatives to point
+ glDers = utils.globalDers(strip.layer,strip.meas,vmeas,wmeas,tDirMeas,normalMeas)
+ ders = glDers.getDers(track.isTop())
+ labGlobal = ders['labels']
+ addDer = ders['ders']
+ if debug:
+ print 'global derivatives:'
+ print labGlobal
+ print addDer
+ point.addGlobals(labGlobal, addDer)
+ #####
+
+ */
+ if (_debug) {
+ System.out.println("HpsGblFitter: " + "uRes " + strip.getId() + " uRes " + uRes + " pred (" + strip.getTrackPos().x() + "," + strip.getTrackPos().y() + "," + strip.getTrackPos().z() + ") s(3D) " + strip.getPath3D());
+ }
+
+ //go to next point
+ s += step;
+
+ } //strips
+
+ //create the trajectory
+ GblTrajectory traj = new GblTrajectory(listOfPoints); //,seedLabel, clSeed);
+
+ if (!traj.isValid()) {
+ System.out.println("HpsGblFitter: " + " Invalid GblTrajectory -> skip");
+ return 1;//INVALIDTRAJ;
+ }
+
+ // print the trajectory
+ if (_debug) {
+ System.out.println("%%%% Gbl Trajectory ");
+ traj.printTrajectory(1);
+ traj.printData();
+ traj.printPoints(4);
+ }
+ // fit trajectory
+ double[] dVals = new double[2];
+ int[] iVals = new int[1];
+ traj.fit(dVals, iVals, "");
+ if (_debug) {
+ System.out.println("HpsGblFitter: Chi2 " + " Fit: " + dVals[0] + ", " + iVals[0] + ", " + dVals[1]);
+ }
+
+ Vector aCorrection = new Vector(5);
+ SymMatrix aCovariance = new SymMatrix(5);
+ traj.getResults(1, aCorrection, aCovariance);
+ System.out.println(" cor ");
+ aCorrection.print(6, 4);
+ System.out.println(" cov ");
+ aCovariance.print(6, 4);
+
+// // write to MP binary file
+// traj.milleOut(mille);
+//
+ return 0;
+ }
+
+ @Override
+ protected void detectorChanged(Detector detector)
+ {
+ }
+
+ private Matrix gblSimpleJacobianLambdaPhi(double ds, double cosl, double bfac)
+ {
+ /**
+ * Simple jacobian: quadratic in arc length difference. using lambda phi
+ * as directions
+ *
+ * @param ds: arc length difference
+ * @type ds: float
+ * @param cosl: cos(lambda)
+ * @type cosl: float
+ * @param bfac: Bz*c
+ * @type bfac: float
+ * @return: jacobian to move by 'ds' on trajectory
+ * @rtype: matrix(float) ajac(1,1)= 1.0D0 ajac(2,2)= 1.0D0
+ * ajac(3,1)=-DBLE(bfac*ds) ajac(3,3)= 1.0D0
+ * ajac(4,1)=-DBLE(0.5*bfac*ds*ds*cosl) ajac(4,3)= DBLE(ds*cosl)
+ * ajac(4,4)= 1.0D0 ajac(5,2)= DBLE(ds) ajac(5,5)= 1.0D0 ''' jac =
+ * np.eye(5) jac[2, 0] = -bfac * ds jac[3, 0] = -0.5 * bfac * ds * ds *
+ * cosl jac[3, 2] = ds * cosl jac[4, 1] = ds return jac
+ */
+ Matrix jac = new Matrix(5, 5);
+ jac.UnitMatrix();
+ jac.set(2, 0, -bfac * ds);
+ jac.set(3, 0, -0.5 * bfac * ds * ds * cosl);
+ jac.set(3, 2, ds * cosl);
+ jac.set(4, 1, ds);
+ return jac;
+ }
+
+}
SVNspam 0.1