Author: [log in to unmask]
Date: Wed Mar 23 18:20:27 2016
New Revision: 4328
Log:
fix bug in residual calculation
Modified:
java/trunk/recon/src/main/java/org/hps/recon/vertexing/BilliorVertexer.java
Modified: java/trunk/recon/src/main/java/org/hps/recon/vertexing/BilliorVertexer.java
=============================================================================
--- java/trunk/recon/src/main/java/org/hps/recon/vertexing/BilliorVertexer.java (original)
+++ java/trunk/recon/src/main/java/org/hps/recon/vertexing/BilliorVertexer.java Wed Mar 23 18:20:27 2016
@@ -24,40 +24,25 @@
// default is a constant field along the z axis
private boolean _debug = false;
- private double _bField;
- private boolean _beamspotConstraint = true;
- private boolean _targetConstraint = false;
- private double[] _beamSize = {0.001, 0.01, 0.01}; //10um in y and z
+ private final double _bField;
+ private boolean _beamspotConstraint;
+ private boolean _targetConstraint;
+ private String _constraintType;
+ private final double[] _beamSize = {0.001, 0.01, 0.01}; //10um in y and z
+ private final double[] _beamPosition = {0.0, 0.0, 0.0}; //origin
private int _ntracks;
- private List<Matrix> paramList = new ArrayList<Matrix>();
- private List<Matrix> WList = new ArrayList<Matrix>();
- private List<Matrix> DList = new ArrayList<Matrix>();
- private List<Matrix> EList = new ArrayList<Matrix>();
- private Matrix A;
- private Matrix T;
- private List<Matrix> BList = new ArrayList<Matrix>();
- private List<Matrix> CinvList = new ArrayList<Matrix>();
- private List<Matrix> CList = new ArrayList<Matrix>();
- private List<Matrix> UList = new ArrayList<Matrix>();
- private List<Matrix> dqList = new ArrayList<Matrix>();
- private double[] _v0 = {0.0, 0.0, 0.0};
+ private double[] _v0 = {0.0, 0.0, 0.0}; //initial guess for unconstrained vertex fit
// private double[] _vertexPosition = {0., 0.0, 0.0};
- private Matrix _vertexPosition = new BasicMatrix(3, 1);
- private Matrix _covVtx = new BasicMatrix(3, 3);
- private List<Matrix> _pFit = new ArrayList<Matrix>();
+ private Matrix _vertexPosition;
+ private Matrix _covVtx;
+ private List<Matrix> _pFit;
;//theta,phi_v,rho
- private List<Matrix> covVtxMomList = new ArrayList<Matrix>();
- private Matrix[][] covMomList = new Matrix[2][2];//max 2 tracks...just make this bigger for more
+ private List<Matrix> covVtxMomList;
+ private Matrix[][] covMomList;//max 2 tracks...just make this bigger for more
private Matrix _constrainedFit;
private Matrix _constrainedCov;
private double _chiSq;
- private String _constraintType = "Unspecified";
-
- // constructor
-
- public BilliorVertexer() {
- }
public BilliorVertexer(double bField) {
_bField = bField;
@@ -89,9 +74,9 @@
_ntracks = tracks.size();
follow1985Paper(tracks);
if (_beamspotConstraint) {
- addV0fromBSConstraint(true);
+ applyBSconstraint(true);
} else if (_targetConstraint) {
- addV0fromBSConstraint(false);
+ applyBSconstraint(false);
}
Map<Integer, Hep3Vector> pFitMap = new HashMap<Integer, Hep3Vector>();
for (int i = 0; i < tracks.size(); i++) {
@@ -104,309 +89,15 @@
return new BilliorVertex(vertDet, covVtxDet, _chiSq, getInvMass(), pFitMap, _constraintType);
}
- public BilliorVertex fitFastVertex(List<BilliorTrack> tracks) {
- _ntracks = tracks.size();
- fastVertex(tracks);
- Hep3Vector vert = new BasicHep3Vector(_vertexPosition.e(0, 0), _vertexPosition.e(1, 0), _vertexPosition.e(2, 0));
- return new BilliorVertex((Hep3Vector) _vertexPosition, _covVtx, _chiSq, getInvMass());
- }
-
-// private void calculateCovariance() {
-// for (int i = 0; i < _ntracks; i++) {
-// BasicMatrix b = (BasicMatrix) BList.get(i);
-// BasicMatrix cinv = (BasicMatrix) CinvList.get(i);
-// BasicMatrix bt = (BasicMatrix) MatrixOp.transposed(b);
-// covVtxMomList.add((MatrixOp.mult(-1, MatrixOp.mult(_covVtx, MatrixOp.mult(b, cinv)))));
-// for (int j = 0; j < _ntracks; j++) {
-// BasicMatrix bj = (BasicMatrix) BList.get(j);
-// BasicMatrix cjinv = (BasicMatrix) CinvList.get(j);
-// BasicMatrix tmp = (BasicMatrix) MatrixOp.mult(cinv, MatrixOp.mult(bt, MatrixOp.mult(_covVtx, MatrixOp.mult(bj, cjinv))));
-// if (i == j)
-// tmp = (BasicMatrix) MatrixOp.add(tmp, cinv);
-// covMomList[i][j] = tmp;
-// }
-// }
-// }
-//
-// private void calculateMomenta() {
-//
-// for (int i = 0; i < _ntracks; i++) {
-// BasicMatrix params = (BasicMatrix) paramList.get(i);
-// BasicMatrix b = (BasicMatrix) BList.get(i);
-// BasicMatrix cinv = (BasicMatrix) CinvList.get(i);
-// BasicMatrix u = (BasicMatrix) UList.get(i);
-// //not sure following line is correct...mg 10/21/10
-// BasicMatrix CinvU = (BasicMatrix) MatrixOp.mult(cinv, u);
-// BasicMatrix CinvBTdV = (BasicMatrix) MatrixOp.mult(-1, MatrixOp.mult(cinv, MatrixOp.mult(MatrixOp.transposed(b), _vertexPosition)));
-//// if(_debug)System.out.println(" B = "+b.toString());
-//// if(_debug)System.out.println(" cinv = "+cinv.toString());
-//// if(_debug)System.out.println(" u = "+u.toString());
-//// if(_debug)System.out.println(" CinvU = "+CinvU.toString());
-//// if(_debug)System.out.println(" CinvBTdV = "+CinvBTdV.toString());
-// BasicMatrix tmpP = (BasicMatrix) MatrixOp.add(CinvBTdV, CinvU);
-// tmpP.setElement(0, 0, tmpP.e(0, 0) + params.e(2, 0));
-// tmpP.setElement(1, 0, tmpP.e(1, 0) + params.e(3, 0));
-// tmpP.setElement(2, 0, tmpP.e(2, 0) + params.e(4, 0));
-// _pFit.add(tmpP);
-//// if(_debug)System.out.println("Track "+i+" orig parameters = "+params);
-//// if(_debug)System.out.println("Track "+i+" deltaP = "+MatrixOp.add(CinvBTdV, CinvU));
-//// if(_debug)System.out.println("Track " + i + " _pFit = " + tmpP);
-// }
-// }
-//
-// private void calculateVertexPosition() {
-// BasicMatrix tmpcov = new BasicMatrix(3, 3);
-// BasicMatrix tmp = new BasicMatrix(3, 1);
-// for (int i = 0; i < _ntracks; i++) {
-// BasicMatrix b = (BasicMatrix) BList.get(i);
-// BasicMatrix cinv = (BasicMatrix) CinvList.get(i);
-// BasicMatrix u = (BasicMatrix) UList.get(i);
-//// if(_debug)System.out.println("Cinv matrix " + cinv.toString());
-//// if(_debug)System.out.println("B matrix " + b.toString());
-//// if(_debug)System.out.println("U matrix " + u.toString());
-// BasicMatrix bt = (BasicMatrix) MatrixOp.transposed(b);
-// // if(_debug)System.out.println("Adding this to tmpcov : " + MatrixOp.mult(-1, MatrixOp.mult(b, MatrixOp.mult(cinv, bt))));
-// if (i == 0) {
-// tmpcov = (BasicMatrix) MatrixOp.mult(-1, MatrixOp.mult(b, MatrixOp.mult(cinv, bt)));
-// tmp = (BasicMatrix) MatrixOp.mult(-1, MatrixOp.mult(b, MatrixOp.mult(cinv, u)));
-// } else {
-// tmpcov = (BasicMatrix) MatrixOp.add(tmpcov, MatrixOp.mult(-1, MatrixOp.mult(b, MatrixOp.mult(cinv, bt))));
-// tmp = (BasicMatrix) MatrixOp.add(tmp, MatrixOp.mult(-1, MatrixOp.mult(b, MatrixOp.mult(cinv, u))));
-// }
-//// if(_debug)System.out.println("tmpCov matrix " + tmpcov.toString());
-//// if(_debug)System.out.println("tmp matrix " + tmp.toString());
-// }
-////
-//// if(_debug)System.out.println("A matrix " + A.toString());
-//// if(_debug)System.out.println("tmpCov matrix " + tmpcov.toString());
-//// if(_debug)System.out.println("sum of A and tmpCov = " + MatrixOp.add(A, tmpcov).toString());
-// _covVtx = MatrixOp.inverse(MatrixOp.add(A, tmpcov));
-//// if(_debug)System.out.println("_covVtx matrix " + _covVtx.toString());
-//// if(_debug)System.out.println("T matrix " + T.toString());
-// _vertexPosition = (BasicMatrix) MatrixOp.mult(_covVtx, MatrixOp.add(T, tmp));
-//
-// }
-//
-// private void makeOtherMatrices() {
-// BasicMatrix tmpA = new BasicMatrix(3, 3);
-// BasicMatrix tmpT = new BasicMatrix(3, 1);
-//
-// for (int i = 0; i < _ntracks; i++) {
-// BasicMatrix tmpD = (BasicMatrix) DList.get(i);
-// BasicMatrix tmpE = (BasicMatrix) EList.get(i);
-// BasicMatrix dq = (BasicMatrix) dqList.get(i);
-// BasicMatrix tmpW = (BasicMatrix) WList.get(i);
-//
-// if (i == 0) {
-// tmpA = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(tmpD), MatrixOp.mult(tmpW, tmpD));
-// tmpT = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(tmpD), MatrixOp.mult(tmpW, dq));
-// } else {
-// tmpT = (BasicMatrix) MatrixOp.add(tmpT, MatrixOp.mult(MatrixOp.transposed(tmpD), MatrixOp.mult(tmpW, dq)));
-// tmpA = (BasicMatrix) MatrixOp.add(tmpA, MatrixOp.mult(MatrixOp.transposed(tmpD), MatrixOp.mult(tmpW, tmpD)));
-// }
-// BList.add(MatrixOp.mult(MatrixOp.transposed(tmpD), MatrixOp.mult(tmpW, tmpE)));
-// BasicMatrix tmpC = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(tmpE), MatrixOp.mult(tmpW, tmpE));
-// CList.add(tmpC);
-// CinvList.add(MatrixOp.inverse(tmpC));
-// UList.add(MatrixOp.mult(MatrixOp.transposed(tmpE), MatrixOp.mult(tmpW, dq)));
-//
-// }
-// A = tmpA;
-// T = tmpT;
-// }
-//
-// private void calculateChisq() {
-// _chiSq = 0;
-// for (int i = 0; i < _ntracks; i++) {
-// BasicMatrix params = (BasicMatrix) paramList.get(i);
-// BasicMatrix d = (BasicMatrix) DList.get(i);
-// BasicMatrix e = (BasicMatrix) EList.get(i);
-// BasicMatrix w = (BasicMatrix) WList.get(i);
-// BasicMatrix pi = (BasicMatrix) _pFit.get(i);
-// BasicMatrix Vtilde = (BasicMatrix) MatrixOp.mult(d, _vertexPosition);
-// BasicMatrix Trtilde = (BasicMatrix) MatrixOp.mult(e, pi);
-// BasicMatrix ptilde = (BasicMatrix) MatrixOp.add(Vtilde, Trtilde);
-// // if(_debug)System.out.println("Vtilde = "+Vtilde);
-// // if(_debug)System.out.println("Trtilde = "+Trtilde);
-// BasicMatrix resid = (BasicMatrix) MatrixOp.add(params, MatrixOp.mult(-1, ptilde));
-// BasicMatrix residT = (BasicMatrix) MatrixOp.transposed(resid);
-//// if(_debug)System.out.println("ptilde = "+ptilde);
-//// if(_debug)System.out.println("params = "+params);
-//// if(_debug)System.out.println("resid = "+resid);
-//// if(_debug)System.out.println("Covariance = "+MatrixOp.inverse(w));
-//// if(_debug)System.out.println("Weight = "+w);
-// _chiSq = _chiSq + (MatrixOp.mult(residT, MatrixOp.mult(w, resid))).e(0, 0);
-//// if(_debug)System.out.println("_chiSq = "+_chiSq);
-// }
-// }
- private void fastVertex(List<BilliorTrack> tracks) {
- boolean firstTrack = true;
- BasicMatrix sumwi = new BasicMatrix(3, 3);
- BasicMatrix sumwiXi = new BasicMatrix(3, 1);
- BasicMatrix dX = new BasicMatrix(3, 1);
-
- for (BilliorTrack bt : tracks) {
- double[] par = bt.parameters();
-// if(_debug)System.out.println("Track parameters = (" + par[0] + ", " + par[1] + ", " + par[2] + ", " + par[3] + ", " + par[4] + ")");
- double cotth = 1. / tan(par[2]);
- double phiv = par[3];
- double cosf = cos(phiv);
- double sinf = sin(phiv);
-
- double xi = par[0] * sin(par[3]);
- double yi = -par[0] * cos(par[3]);
- double zi = par[1];
-
- dX.setElement(0, 0, xi);
- dX.setElement(1, 0, yi);
- dX.setElement(2, 0, zi);
-
- BasicMatrix tmpD = new BasicMatrix(2, 3);
- tmpD.setElement(0, 0, sinf);
- tmpD.setElement(0, 1, -cosf);
- tmpD.setElement(1, 0, -cotth * cosf);
- tmpD.setElement(1, 1, -cotth * sinf);
- tmpD.setElement(1, 2, 1);
- BasicMatrix trkCov = new BasicMatrix(2, 2);
- trkCov.setElement(0, 0, bt.covariance().e(0, 0));
- trkCov.setElement(0, 1, bt.covariance().e(0, 1));
- trkCov.setElement(1, 0, bt.covariance().e(1, 0));
- trkCov.setElement(1, 1, bt.covariance().e(1, 1));
- BasicMatrix tmpW = (BasicMatrix) MatrixOp.inverse(trkCov);
- BasicMatrix wi = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(tmpD), MatrixOp.mult(tmpW, tmpD));
- if (firstTrack) {
- sumwi = wi;
- sumwiXi = (BasicMatrix) MatrixOp.mult(wi, dX);
- } else {
- sumwi = (BasicMatrix) MatrixOp.add(sumwi, wi);
- sumwiXi = (BasicMatrix) MatrixOp.add(sumwiXi, MatrixOp.mult(wi, dX));
- }
- firstTrack = false;
- }
- _covVtx = MatrixOp.inverse(sumwi);
- if (_debug) {
- System.out.println("fastVertex::_covVtx matrix " + _covVtx.toString());
- }
- _vertexPosition = (BasicMatrix) MatrixOp.mult(_covVtx, sumwiXi);
- _chiSq = 0;
- //get the chisq
- for (BilliorTrack bt : tracks) {
- double[] par = bt.parameters();
-// if(_debug)System.out.println("Track parameters = (" + par[0] + ", " + par[1] + ", " + par[2] + ", " + par[3] + ", " + par[4] + ")");
- double cotth = 1. / tan(par[2]);
- double phiv = par[3];
- double cosf = cos(phiv);
- double sinf = sin(phiv);
-
- double xi = par[0] * sin(par[3]);
- double yi = -par[0] * cos(par[3]);
- double zi = par[1];
- //this is xi - fitted vertex now
- dX.setElement(0, 0, xi - _vertexPosition.e(0, 0));
- dX.setElement(1, 0, yi - _vertexPosition.e(1, 0));
- dX.setElement(2, 0, zi - _vertexPosition.e(2, 0));
-
- BasicMatrix tmpD = new BasicMatrix(2, 3);
- tmpD.setElement(0, 0, sinf);
- tmpD.setElement(0, 1, -cosf);
- tmpD.setElement(1, 0, -cotth * cosf);
- tmpD.setElement(1, 1, -cotth * sinf);
- tmpD.setElement(1, 2, 1);
- BasicMatrix trkCov = new BasicMatrix(2, 2);
- trkCov.setElement(0, 0, bt.covariance().e(0, 0));
- trkCov.setElement(0, 1, bt.covariance().e(0, 1));
- trkCov.setElement(1, 0, bt.covariance().e(1, 0));
- trkCov.setElement(1, 1, bt.covariance().e(1, 1));
- BasicMatrix tmpW = (BasicMatrix) MatrixOp.inverse(trkCov);
- BasicMatrix wi = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(tmpD), MatrixOp.mult(tmpW, tmpD));
- _chiSq += MatrixOp.mult(MatrixOp.transposed(dX), MatrixOp.mult(wi, dX)).e(0, 0);
- }
- }
-
-// private void makeDerivativeMatrices(List<BilliorTrack> tracks) {
-//
-// //DList.clear();
-// //EList.clear();
-// //paramList.clear();
-// //dqList.clear();
-// //WList.clear();
-// BasicMatrix dq = new BasicMatrix(5, 1);
-// BasicMatrix tmpW = new BasicMatrix(5, 5);
-// for (BilliorTrack bt : tracks) {
-// double[] par = bt.parameters();
-// BasicMatrix tmpPar = new BasicMatrix(5, 1);
-// tmpPar.setElement(0, 0, par[0]);
-// tmpPar.setElement(1, 0, par[1]);
-// tmpPar.setElement(2, 0, par[2]);
-// tmpPar.setElement(3, 0, par[3]);
-// tmpPar.setElement(4, 0, par[4]);
-// paramList.add(tmpPar);
-// double cotth = 1. / tan(par[2]);
-// double uu = _v0[0] * cos(par[3]) + _v0[1] * sin(par[3]);//Q
-// double vv = _v0[1] * cos(par[3]) - _v0[0] * sin(par[3]);//R
-// double eps = -vv - .5 * uu * uu * par[4];
-// double zp = _v0[2] - uu * (1 - vv * par[4]) * cotth;
-// // * phi at vertex with these parameters
-// double phiv = par[3] + uu * par[4];
-// double cosf = cos(phiv);
-// double sinf = sin(phiv);
-//
-// BasicMatrix tmpD = new BasicMatrix(5, 3);
-// tmpD.setElement(0, 0, sinf);
-// tmpD.setElement(0, 1, -cosf);
-// tmpD.setElement(1, 0, -cotth * cosf);
-// tmpD.setElement(1, 1, -cotth * sinf);
-// tmpD.setElement(1, 2, 1);
-// tmpD.setElement(3, 0, -par[4] * cosf);
-// tmpD.setElement(3, 1, -par[4] * sinf);
-//
-// BasicMatrix tmpE = new BasicMatrix(5, 3);
-// tmpE.setElement(0, 1, uu);
-// tmpE.setElement(0, 2, -uu * uu / 2);
-// tmpE.setElement(1, 0, uu * (1 + cotth * cotth));
-// tmpE.setElement(1, 1, -vv * cotth);
-// tmpE.setElement(1, 2, uu * vv * cotth);
-// tmpE.setElement(3, 1, 1);
-// tmpE.setElement(3, 2, -uu);
-// tmpE.setElement(2, 0, 1); //partial(theta)/dtheta
-// tmpE.setElement(4, 2, 1); //partial (rho)/drho
-// DList.add(tmpD);
-// EList.add(tmpE);
-//
-// double deps = par[0] - eps;
-// double dzp = par[1] - zp;
-// double dphi = par[3] - phiv;
-//
-// dq.setElement(0, 0, deps);
-// dq.setElement(1, 0, dzp);
-// dq.setElement(3, 0, dphi);
-// dqList.add(dq);
-// tmpW = (BasicMatrix) MatrixOp.inverse(bt.covariance());
-// WList.add(tmpW);
-//
-// if (_debug)
-// System.out.println("makeDerivativeMatrices::Params = \n" + tmpPar);
-// if (_debug)
-// System.out.println("D = \n" + tmpD);
-// if (_debug)
-// System.out.println("E = \n" + tmpE);
-// if (_debug)
-// System.out.println("dq = \n" + dq);
-// if (_debug)
-// System.out.println("W = \n" + tmpW);
-// }
-//
-// }
-
- /* Add the constraint that V0 points back to beamspot
+ /* Add the constraint that V0 is at/points back to beamspot
* this method is based on progressive least squares fit
* using the unconstrained fit result as the (k-1) fit
*
* all notation is taken from:
* W. Hulsbergen, NIM 552 (2005) 566-575
*/
- private void addV0fromBSConstraint(boolean bscon) {
- String methodName = bscon ? "addV0fromBSConstraint" : "constrainV0toBS";
+ private void applyBSconstraint(boolean pointback) {
+ String methodName = pointback ? "constrainV0toBS" : "constrainV0toTarget";
BasicMatrix Ckm1 = new BasicMatrix(3 * (_ntracks + 1), 3 * (_ntracks + 1));
BasicMatrix Xkm1 = new BasicMatrix(3 * (_ntracks + 1), 1);
MatrixOp.setSubMatrix(Ckm1, _covVtx, 0, 0);
@@ -457,9 +148,15 @@
pztot += pz;
}
//calculate the position of the A' at X=0
- BasicMatrix rk = makeRk(Vx, Vy, Vz, pxtot, pytot, pztot, bscon);
-
- BasicMatrix Hk = makeHk(_ntracks, pxtot, pytot, pztot, bscon);
+ BasicMatrix rk = makeRk(Vx, Vy, Vz, pxtot, pytot, pztot, pointback);
+ if (_debug) {
+ System.out.println(methodName + "::rk = " + rk);
+ }
+
+ BasicMatrix Hk = makeHk(_ntracks, pxtot, pytot, pztot, pointback);
+ if (_debug) {
+ System.out.println(methodName + "::Hk = " + Hk);
+ }
// the beam covariance
BasicMatrix Vk = new BasicMatrix(3, 3);
@@ -511,6 +208,10 @@
// if(_debug)System.out.println("With Constraint : " + _vertexPosition.toString());
// if(_debug)System.out.println("With Constraint : " + _covVtx.toString());
+ if (_debug) {
+ System.out.println("Constrained vertex: " + _vertexPosition);
+ }
+
for (int i = 0; i < _ntracks; i++) {
BasicMatrix ptmp = (BasicMatrix) MatrixOp.getSubMatrix(_constrainedFit, 3 * (i + 1), 0, 3, 1);
_pFit.set(i, ptmp);
@@ -519,134 +220,24 @@
// if(_debug)System.out.println("Unconstrained chi^2 = "+_chiSq);
//ok...add to the chi^2
if (_debug) {
- System.out.println(MatrixOp.mult(MatrixOp.transposed(rk), MatrixOp.mult(Rkinv, rk)));
+ System.out.println("Chisq contribution: " + MatrixOp.mult(MatrixOp.transposed(rk), MatrixOp.mult(Rkinv, rk)));
}
_chiSq += MatrixOp.mult(MatrixOp.transposed(rk), MatrixOp.mult(Rkinv, rk)).e(0, 0);
// if(_debug)System.out.println("Constrained chi^2 = "+_chiSq);
}
-// private void constrainV0toBS() {
-// BasicMatrix Ckm1 = new BasicMatrix(3 * (_ntracks + 1), 3 * (_ntracks + 1));
-// BasicMatrix Xkm1 = new BasicMatrix(3 * (_ntracks + 1), 1);
-// MatrixOp.setSubMatrix(Ckm1, _covVtx, 0, 0);
-// MatrixOp.setSubMatrix(Xkm1, _vertexPosition, 0, 0);
-//
-// int n = 1;
-// for (Matrix covVtxMom : covVtxMomList) {
-// if (_debug)
-// System.out.println("constrainV0toBS::Track " + n + " covVtxMom : " + covVtxMom.toString());
-// MatrixOp.setSubMatrix(Ckm1, covVtxMom, 0, 3 * n);
-// MatrixOp.setSubMatrix(Ckm1, MatrixOp.transposed(covVtxMom), 3 * n, 0);
-// n++;
-// }
-// for (int i = 0; i < _ntracks; i++) {
-// BasicMatrix pi = (BasicMatrix) _pFit.get(i);
-// MatrixOp.setSubMatrix(Xkm1, pi, 3 * (i + 1), 0);
-// // if(_debug)System.out.println("Track "+i+" p : " + pi.toString());
-// for (int j = 0; j < _ntracks; j++)
-// MatrixOp.setSubMatrix(Ckm1, covMomList[i][j], 3 * (i + 1), 3 * (j + 1));
-// }
-// // now calculate the derivative matrix for the beam constraint.
-// // the beamspot is assumed to be at bvec=(0,0,0)
-// // the V0 production position is Vbvec=(0,0,0)
-// // where ptot=sum_i (pi)
-// // need derivites wrt to the vertex position and momentum (theta,phi_v,rho)
-// double Vx = _vertexPosition.e(0, 0);
-// double Vy = _vertexPosition.e(1, 0);
-// double Vz = _vertexPosition.e(2, 0);
-// //first, get the sum of momenta...
-// double pxtot = 0;
-// double pytot = 0;
-// double pztot = 0;
-// for (int i = 0; i < _ntracks; i++) {
-// BasicMatrix pi = (BasicMatrix) _pFit.get(i);
-// double theta = pi.e(0, 0);
-// double phiv = pi.e(1, 0);
-// double rho = pi.e(2, 0);
-// double Pt = Math.abs((1. / rho) * _bField * Constants.fieldConversion);
-// double px = Pt * Math.cos(phiv);
-// double py = Pt * Math.sin(phiv);
-// double pz = Pt * 1 / Math.tan(theta);
-// pxtot += px;
-// pytot += py;
-// pztot += pz;
-// }
-// BasicMatrix rk = makeRk(Vx,Vy,Vz,pxtot,pytot,pztot,false);
-//// //calculate the position of the A' at X=0
-//// BasicMatrix rk = new BasicMatrix(3, 1);
-//// // if(_debug)System.out.println("Vx = " + Vx + "; Vy = " + Vy + "; Vz = " + Vz + "; pxtot = " + pxtot + "; pytot = " + pytot + "; pztot = " + pztot);
-//// rk.setElement(0, 0, -Vx);
-//// rk.setElement(1, 0, -Vy);
-//// rk.setElement(2, 0, -Vz);
-//
-// BasicMatrix Hk = makeHk(_ntracks,pxtot,pytot,pztot,false);
-//
-// // the beam covariance
-// BasicMatrix Vk = new BasicMatrix(3, 3);
-// Vk.setElement(0, 0, _beamSize[0] * _beamSize[0]);
-// Vk.setElement(1, 1, _beamSize[1] * _beamSize[1]);
-// Vk.setElement(2, 2, _beamSize[2] * _beamSize[2]);
-//
-// //now do the matrix operations to get the constrained parameters
-// BasicMatrix Hkt = (BasicMatrix) MatrixOp.transposed(Hk);
-//// if(_debug)System.out.println("Ckm1Hk = " + MatrixOp.mult(Ckm1, Hk));
-//
-// BasicMatrix Rk = (BasicMatrix) MatrixOp.mult(Hkt, MatrixOp.mult(Ckm1, Hk));
-//// if(_debug)System.out.println("Pre Vk: Rk = " + Rk.toString());
-// Rk = (BasicMatrix) MatrixOp.add(Rk, Vk);
-// BasicMatrix Rkinv = (BasicMatrix) MatrixOp.inverse(Rk);
-// BasicMatrix Kk = (BasicMatrix) MatrixOp.mult(Ckm1, MatrixOp.mult(Hk, Rkinv));
-//
-//// if(_debug)System.out.println("Ckm1 = " + Ckm1.toString());
-//// if(_debug)System.out.println("Hk = " + Hk.toString());
-//// if(_debug)System.out.println("Rk = " + Rk.toString());
-//// if(_debug)System.out.println("Vk = " + Vk.toString());
-//// if(_debug)System.out.println("rk = " + rk.toString());
-//// if(_debug)System.out.println("Kk = " + Kk.toString());
-// _constrainedFit = MatrixOp.mult(Kk, rk);
-// _constrainedFit = MatrixOp.add(_constrainedFit, Xkm1);//Xk
-//
-// //ok, get the new covariance
-// BasicMatrix RkKkt = (BasicMatrix) MatrixOp.mult(Rk, MatrixOp.transposed(Kk));
-// BasicMatrix HkCkm1 = (BasicMatrix) MatrixOp.mult(Hkt, Ckm1);
-// RkKkt = (BasicMatrix) MatrixOp.mult(1, RkKkt);
-// HkCkm1 = (BasicMatrix) MatrixOp.mult(-2, HkCkm1);
-// BasicMatrix sumMatrix = (BasicMatrix) MatrixOp.mult(Kk, MatrixOp.add(HkCkm1, RkKkt));
-// _constrainedCov = (BasicMatrix) MatrixOp.add(Ckm1, sumMatrix);
-//
-// //update the regular parameter names to the constrained result
-//// if(_debug)System.out.println("Without Constraint : " + _vertexPosition.toString());
-//// if(_debug)System.out.println("Without Constraint: x= "+_vertexPosition.e(0,0));
-// // if(_debug)System.out.println(_constrainedFit.toString());
-//// if(_debug)System.out.println("Without Constraint : " + _covVtx.toString());
-// _vertexPosition = (BasicMatrix) MatrixOp.getSubMatrix(_constrainedFit, 0, 0, 3, 1);
-// _covVtx = (BasicMatrix) MatrixOp.getSubMatrix(_constrainedCov, 0, 0, 3, 3);
-//// if(_debug)System.out.println("With Constraint : " + _vertexPosition.toString());
-//// if(_debug)System.out.println("With Constraint : " + _covVtx.toString());
-//
-// for (int i = 0; i < _ntracks; i++) {
-// BasicMatrix ptmp = (BasicMatrix) MatrixOp.getSubMatrix(_constrainedFit, 3 * (i + 1), 0, 3, 1);
-// _pFit.set(i, ptmp);
-// }
-//
-//// if(_debug)System.out.println("Unconstrained chi^2 = "+_chiSq);
-// //ok...add to the chi^2
-// if (_debug)
-// System.out.println(MatrixOp.mult(MatrixOp.transposed(rk), MatrixOp.mult(Rkinv, rk)));
-// _chiSq += MatrixOp.mult(MatrixOp.transposed(rk), MatrixOp.mult(Rkinv, rk)).e(0, 0);
-//// if(_debug)System.out.println("Constrained chi^2 = "+_chiSq);
-// }
private BasicMatrix makeHk(int ntracks, double pxtot, double pytot, double pztot, boolean bscon) {
double Vx = _vertexPosition.e(0, 0);
// double Vy = _vertexPosition.e(1, 0);
// double Vz = _vertexPosition.e(2, 0);
BasicMatrix Hk = new BasicMatrix(3 * (ntracks + 1), 3);
// ok, can set the derivitives wrt to V
- Hk.setElement(0, 0, 0);
if (bscon) {
+ Hk.setElement(0, 0, 0);
Hk.setElement(0, 1, pytot / pxtot);
Hk.setElement(0, 2, pztot / pxtot);
} else {
+ Hk.setElement(0, 0, 1);
Hk.setElement(0, 1, 0);
Hk.setElement(0, 2, 0);
}
@@ -670,7 +261,8 @@
Hk.setElement(3 * (i + 1), 0, 0);
Hk.setElement(3 * (i + 1), 1, 0);
if (bscon) {
- Hk.setElement(3 * (i + 1), 2, -Pt / Math.pow(sin(theta), 2) * Vx);
+ Hk.setElement(3 * (i + 1), 2,
+ -Pt / Math.pow(sin(theta), 2) * Vx);
} else {
Hk.setElement(3 * (i + 1), 2, 0);
}
@@ -679,7 +271,8 @@
if (bscon) {
Hk.setElement(3 * (i + 1) + 1, 1,
(Pt * Pt * cos(phiv) * sin(phiv) / (pxtot * pxtot)) * Vx);
- Hk.setElement(3 * (i + 1) + 1, 2, (Pt * sin(phiv) / (pxtot * pxtot)) * Vx * pztot);
+ Hk.setElement(3 * (i + 1) + 1, 2,
+ (Pt * sin(phiv) / (pxtot * pxtot)) * Vx * pztot);
} else {
Hk.setElement(3 * (i + 1) + 1, 1, 0);
Hk.setElement(3 * (i + 1) + 1, 2, 0);
@@ -710,16 +303,16 @@
//calculate the position of the A' at X=0
BasicMatrix rk = new BasicMatrix(3, 1);
if (_debug) {
- System.out.println("addV0fromBSConstraint::Vx = " + Vx + "; Vy = " + Vy + "; Vz = " + Vz + "; pxtot = " + pxtot + "; pytot = " + pytot + "; pztot = " + pztot);
+ System.out.println("makeRk::Vx = " + Vx + "; Vy = " + Vy + "; Vz = " + Vz + "; pxtot = " + pxtot + "; pytot = " + pytot + "; pztot = " + pztot);
}
if (bscon) {
rk.setElement(0, 0, 0);
- rk.setElement(1, 0, 0 - (Vy - pytot / pxtot * Vx));
- rk.setElement(2, 0, 0 - (Vz - pztot / pxtot * Vx));
+ rk.setElement(1, 0, _beamPosition[1] - (Vy - pytot / pxtot * (Vx - _beamPosition[0])));
+ rk.setElement(2, 0, _beamPosition[2] - (Vz - pztot / pxtot * (Vx - _beamPosition[0])));
} else {
- rk.setElement(0, 0, -Vx);
- rk.setElement(1, 0, -Vy);
- rk.setElement(2, 0, -Vz);
+ rk.setElement(0, 0, _beamPosition[0] - Vx);
+ rk.setElement(1, 0, _beamPosition[1] - Vy);
+ rk.setElement(2, 0, _beamPosition[2] - Vz);
}
return rk;
}
@@ -734,13 +327,21 @@
_beamSize[2] = bs[2];
}
+ public void setBeamPosition(double[] bp) {
+ _beamPosition[0] = bp[0];
+ _beamPosition[1] = bp[1];
+ _beamPosition[2] = bp[2];
+ }
+
public void doBeamSpotConstraint(boolean bsconst) {
_beamspotConstraint = bsconst;
+ _targetConstraint = false;
_constraintType = "BeamspotConstrained";
}
public void doTargetConstraint(boolean bsconst) {
+ _beamspotConstraint = false;
_targetConstraint = bsconst;
_constraintType = "TargetConstrained";
}
@@ -797,9 +398,10 @@
}
}
+ @Override
public String toString() {
- StringBuffer sb = new StringBuffer("Vertex at : \nx= " + _vertexPosition.e(0, 0) + " +/- " + Math.sqrt(_covVtx.e(0, 0)) + "\ny= " + _vertexPosition.e(1, 0) + " +/- " + Math.sqrt(_covVtx.e(1, 1)) + "\nz= " + _vertexPosition.e(2, 0) + " +/- " + Math.sqrt(_covVtx.e(2, 2)));
- return sb.toString();
+ String sb = "Vertex at : \nx= " + _vertexPosition.e(0, 0) + " +/- " + Math.sqrt(_covVtx.e(0, 0)) + "\ny= " + _vertexPosition.e(1, 0) + " +/- " + Math.sqrt(_covVtx.e(1, 1)) + "\nz= " + _vertexPosition.e(2, 0) + " +/- " + Math.sqrt(_covVtx.e(2, 2));
+ return sb;
}
private void follow1985Paper(List<BilliorTrack> tracks) {
@@ -808,8 +410,8 @@
v0.setElement(0, 0, _v0[0]);
v0.setElement(1, 0, _v0[1]);
v0.setElement(2, 0, _v0[2]);
- List<Matrix> params = new ArrayList<Matrix>();
- List<Matrix> q0s = new ArrayList<Matrix>();
+// List<Matrix> params = new ArrayList<Matrix>();
+// List<Matrix> q0s = new ArrayList<Matrix>();
List<Matrix> Gs = new ArrayList<Matrix>();
List<Matrix> Ds = new ArrayList<Matrix>();
List<Matrix> Es = new ArrayList<Matrix>();
@@ -828,11 +430,11 @@
tmpPar.setElement(2, 0, par[2]);
tmpPar.setElement(3, 0, par[3]);
tmpPar.setElement(4, 0, par[4]);
- params.add(tmpPar);
+// params.add(tmpPar);
double theta = par[2];
- double phiv = par[3];
+// double phiv = par[3];
double rho = par[4];
- double Pt = Math.abs((1. / rho) * _bField * Constants.fieldConversion);
+// double Pt = Math.abs((1. / rho) * _bField * Constants.fieldConversion);
double cotth = 1. / tan(par[2]);
double uu = v0.e(0, 0) * cos(par[3]) + v0.e(1, 0) * sin(par[3]);//Q
@@ -858,7 +460,7 @@
q0.setElement(0, 0, theta);
q0.setElement(1, 0, phiVert);
q0.setElement(2, 0, rho);
- q0s.add(q0);
+// q0s.add(q0);
double cosf = cos(phiVert);
double sinf = sin(phiVert);
@@ -935,8 +537,8 @@
}
//ok, now the momentum
- List<Matrix> qtildes = new ArrayList<Matrix>();
- List<Matrix> ptildes = new ArrayList<Matrix>();
+// List<Matrix> qtildes = new ArrayList<Matrix>();
+// List<Matrix> ptildes = new ArrayList<Matrix>();
List<Matrix> C0j = new ArrayList<Matrix>();
List<Matrix> pfit = new ArrayList<Matrix>();
Matrix[][] Cij = new Matrix[2][2];//max 2 tracks...just make this bigger for more
@@ -954,9 +556,9 @@
BasicMatrix second = (BasicMatrix) MatrixOp.mult(MatrixOp.inverse(e), MatrixOp.mult(MatrixOp.transposed(b), g));
second = (BasicMatrix) MatrixOp.mult(second, p);
BasicMatrix qtilde = (BasicMatrix) MatrixOp.add(first, second);
- qtildes.add(qtilde);
+// qtildes.add(qtilde);
BasicMatrix ptilde = (BasicMatrix) MatrixOp.add(MatrixOp.mult(a, xtilde), MatrixOp.mult(b, qtilde));
- ptildes.add(ptilde);
+// ptildes.add(ptilde);
chisq += MatrixOp.mult(MatrixOp.transposed(MatrixOp.add(p, MatrixOp.mult(-1, ptilde))), MatrixOp.mult(g, MatrixOp.add(p, MatrixOp.mult(-1, ptilde)))).e(0, 0);
if (_debug) {
System.out.println("\n\nfollow1985Paper::Track #" + j);
@@ -980,12 +582,12 @@
BasicMatrix tmpC0j = (BasicMatrix) MatrixOp.mult(-1, MatrixOp.mult(covVtx, MatrixOp.mult(d, MatrixOp.inverse(e))));
C0j.add(tmpC0j);
for (int i = 0; i < _ntracks; i++) {
- BasicMatrix ai = (BasicMatrix) As.get(i);
- BasicMatrix bi = (BasicMatrix) Bs.get(i);
- BasicMatrix di = (BasicMatrix) Ds.get(i);
- BasicMatrix ei = (BasicMatrix) Es.get(i);
- BasicMatrix gi = (BasicMatrix) Gs.get(i);
- BasicMatrix pi = (BasicMatrix) pis.get(i);
+// BasicMatrix ai = (BasicMatrix) As.get(i);
+// BasicMatrix bi = (BasicMatrix) Bs.get(i);
+// BasicMatrix di = (BasicMatrix) Ds.get(i);
+// BasicMatrix ei = (BasicMatrix) Es.get(i);
+// BasicMatrix gi = (BasicMatrix) Gs.get(i);
+// BasicMatrix pi = (BasicMatrix) pis.get(i);
BasicMatrix tmpCij = (BasicMatrix) MatrixOp.mult(-1, MatrixOp.mult(MatrixOp.inverse(e), MatrixOp.mult(MatrixOp.transposed(d), tmpC0j)));
Cij[i][j] = tmpCij;
}
