Author: [log in to unmask]
Date: Tue Mar 22 17:14:50 2016
New Revision: 4324
Log:
clean up vertexing code
Modified:
java/trunk/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java
java/trunk/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java
java/trunk/recon/src/main/java/org/hps/recon/vertexing/BilliorVertexer.java
Modified: java/trunk/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java
=============================================================================
--- java/trunk/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java (original)
+++ java/trunk/recon/src/main/java/org/hps/recon/particle/HpsReconParticleDriver.java Tue Mar 22 17:14:50 2016
@@ -242,6 +242,7 @@
BilliorVertexer vtxFitter = new BilliorVertexer(bField);
// TODO: The beam size should come from the conditions database.
vtxFitter.setBeamSize(beamSize);
+ vtxFitter.setDebug(debug);
// Perform the vertexing based on the specified constraint.
switch (constraint) {
Modified: java/trunk/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java
=============================================================================
--- java/trunk/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java (original)
+++ java/trunk/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java Tue Mar 22 17:14:50 2016
@@ -573,7 +573,7 @@
/**
* Indicates whether debug text should be output or not.
*/
- private boolean debug = false;
+ protected boolean debug = false;
/**
* The simple name of the class used for debug print statements.
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 Tue Mar 22 17:14:50 2016
@@ -16,7 +16,8 @@
/**
* @version $Id: BilliorVertexer.java,v 1.3 2013/03/13 19:24:20 mgraham Exp $
- * @version Vertex tracks using least-squares method laid out by billior etal used in the HPS Java package.
+ * @version Vertex tracks using least-squares method laid out by billior etal
+ * used in the HPS Java package.
*/
public class BilliorVertexer {
// the value of the magnetic field in the vicinity of the vertex
@@ -51,184 +52,193 @@
private Matrix _constrainedFit;
private Matrix _constrainedCov;
private double _chiSq;
- private String _constraintType="Unspecified";
+ private String _constraintType = "Unspecified";
+
// constructor
+
public BilliorVertexer() {
}
public BilliorVertexer(double bField) {
_bField = bField;
- _constraintType="Unconstrained";
- _beamspotConstraint =false;
+ _constraintType = "Unconstrained";
+ _beamspotConstraint = false;
_targetConstraint = false;
}
-
- public BilliorVertexer(double bField,boolean bsConst, boolean constToBS) {
+
+ public BilliorVertexer(double bField, boolean bsConst, boolean constToBS) {
_bField = bField;
- _beamspotConstraint =bsConst;
+ _beamspotConstraint = bsConst;
_targetConstraint = constToBS;
- if(_beamspotConstraint&&_targetConstraint)
+ if (_beamspotConstraint && _targetConstraint) {
System.out.println("BilliorVertexer::Warning!!! Setting both _beamspotConstraint and _targetConstraint to true!");
- if(_beamspotConstraint)
- _constraintType="BeamspotConstrained";
- if(_targetConstraint)
- _constraintType="TargetConstrained";
+ }
+ if (_beamspotConstraint) {
+ _constraintType = "BeamspotConstrained";
+ }
+ if (_targetConstraint) {
+ _constraintType = "TargetConstrained";
+ }
+ }
+
+ public void setDebug(boolean debug) {
+ _debug = debug;
}
public BilliorVertex fitVertex(List<BilliorTrack> tracks) {
_ntracks = tracks.size();
follow1985Paper(tracks);
- if (_beamspotConstraint)
- addV0fromBSConstraint();
- else if (_targetConstraint)
- constrainV0toBS();
- Map<Integer,Hep3Vector> pFitMap=new HashMap<Integer,Hep3Vector>();
- for(int i=0;i<tracks.size();i++){
- Hep3Vector pFit=new BasicHep3Vector(this.getFittedMomentum(i));
- pFitMap.put(i, pFit);
- }
- Hep3Vector vert=new BasicHep3Vector(_vertexPosition.e(0, 0),_vertexPosition.e(1, 0),_vertexPosition.e(2, 0));
- Hep3Vector vertDet=CoordinateTransformations.transformVectorToDetector(vert);
- SymmetricMatrix covVtxDet=CoordinateTransformations.transformCovarianceToDetector(new SymmetricMatrix( _covVtx));
- return new BilliorVertex(vertDet,covVtxDet,_chiSq,getInvMass(),pFitMap,_constraintType);
+ if (_beamspotConstraint) {
+ addV0fromBSConstraint(true);
+ } else if (_targetConstraint) {
+ addV0fromBSConstraint(false);
+ }
+ Map<Integer, Hep3Vector> pFitMap = new HashMap<Integer, Hep3Vector>();
+ for (int i = 0; i < tracks.size(); i++) {
+ Hep3Vector pFit = new BasicHep3Vector(this.getFittedMomentum(i));
+ pFitMap.put(i, pFit);
+ }
+ Hep3Vector vert = new BasicHep3Vector(_vertexPosition.e(0, 0), _vertexPosition.e(1, 0), _vertexPosition.e(2, 0));
+ Hep3Vector vertDet = CoordinateTransformations.transformVectorToDetector(vert);
+ SymmetricMatrix covVtxDet = CoordinateTransformations.transformCovarianceToDetector(new SymmetricMatrix(_covVtx));
+ 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);
- }
- }
-
+ 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);
@@ -274,8 +284,9 @@
firstTrack = false;
}
_covVtx = MatrixOp.inverse(sumwi);
- if (_debug)
+ if (_debug) {
System.out.println("fastVertex::_covVtx matrix " + _covVtx.toString());
+ }
_vertexPosition = (BasicMatrix) MatrixOp.mult(_covVtx, sumwiXi);
_chiSq = 0;
//get the chisq
@@ -312,80 +323,80 @@
}
}
- 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);
- }
-
- }
+// 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
* this method is based on progressive least squares fit
@@ -394,16 +405,17 @@
* all notation is taken from:
* W. Hulsbergen, NIM 552 (2005) 566-575
*/
- private void addV0fromBSConstraint() {
- BasicMatrix Hk = new BasicMatrix(3 * (_ntracks + 1), 3);
+ private void addV0fromBSConstraint(boolean bscon) {
+ String methodName = bscon ? "addV0fromBSConstraint" : "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("addV0fromBSConstraint::Track " + n + " covVtxMom : " + covVtxMom.toString());
+ if (_debug) {
+ System.out.println(methodName + "::Track " + n + " covVtxMom : " + covVtxMom.toString());
+ }
MatrixOp.setSubMatrix(Ckm1, covVtxMom, 0, 3 * n);
MatrixOp.setSubMatrix(Ckm1, MatrixOp.transposed(covVtxMom), 3 * n, 0);
n++;
@@ -411,10 +423,12 @@
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("addV0fromBSConstraint::Track " + i + " p : " + pi.toString());
- for (int j = 0; j < _ntracks; j++)
+ if (_debug) {
+ System.out.println(methodName + "::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.
@@ -443,56 +457,10 @@
pztot += pz;
}
//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);
- rk.setElement(0, 0, 0);
- rk.setElement(1, 0, 0 - (Vy - pytot / pxtot * Vx));
- rk.setElement(2, 0, 0 - (Vz - pztot / pxtot * Vx));
-
-// ok, can set the derivitives wrt to V
- Hk.setElement(0, 0, 0);
- Hk.setElement(0, 1, pytot / pxtot);
- Hk.setElement(0, 2, pztot / pxtot);
- Hk.setElement(1, 0, 0);
- Hk.setElement(1, 1, 1);
- Hk.setElement(1, 2, 0);
- Hk.setElement(2, 0, 0);
- Hk.setElement(2, 1, 0);
- Hk.setElement(2, 2, 1);
-//ok, loop over tracks again to set the derivitives wrt track momenta (theta,phi,rho)
- 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);
- //derivities wrt theta
- Hk.setElement(3 * (i + 1), 0, 0);
- Hk.setElement(3 * (i + 1), 1, 0);
- Hk.setElement(3 * (i + 1), 2, -Pt / Math.pow(sin(theta), 2) * Vx);
- //derivities wrt phi
- Hk.setElement(3 * (i + 1) + 1, 0, 0);
- 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);
- //derivities wrt rho
- Hk.setElement(3 * (i + 1) + 2, 0, 0);
-// Hk.setElement(3 * (i + 1) + 2, 1,
-// (pytot / pxtot - 1) * (Pt / rho) * (1 / pxtot) * Vx);
-// Hk.setElement(3 * (i + 1) + 2, 2,
-// (pztot / pxtot - 1) * (Pt / rho) * (1 / pxtot) * Vx);
- Hk.setElement(3 * (i + 1) + 2, 1,
- (cos(phiv) * pytot / pxtot - sin(phiv)) * (Pt / rho) * (1 / pxtot) * Vx);
- Hk.setElement(3 * (i + 1) + 2, 2,
- (cos(phiv) * pztot / pxtot - sin(phiv)) * (Pt / rho) * (1 / pxtot) * Vx);
- // if(_debug)System.out.println("pxtot = "+pxtot+"; rho = "+rho+"; Pt = "+Pt);
- // if(_debug)System.out.println("cos(phiv)*pytot / pxtot - sin(phiv) = "+(cos(phiv)*pytot / pxtot - sin(phiv)));
- // if(_debug)System.out.println("Pt/(rho*pxtot) = "+(Pt / rho) * (1 / pxtot));
- }
+ BasicMatrix rk = makeRk(Vx, Vy, Vz, pxtot, pytot, pztot, bscon);
+
+ BasicMatrix Hk = makeHk(_ntracks, pxtot, pytot, pztot, bscon);
+
// the beam covariance
BasicMatrix Vk = new BasicMatrix(3, 3);
Vk.setElement(0, 0, _beamSize[0] * _beamSize[0]);
@@ -501,15 +469,18 @@
//now do the matrix operations to get the constrained parameters
BasicMatrix Hkt = (BasicMatrix) MatrixOp.transposed(Hk);
- if (_debug)
- System.out.println("addV0fromBSConstraint::Ckm1Hk = " + MatrixOp.mult(Ckm1, Hk));
+ if (_debug) {
+ System.out.println(methodName + "::Ckm1Hk = " + MatrixOp.mult(Ckm1, Hk));
+ }
BasicMatrix Rk = (BasicMatrix) MatrixOp.mult(Hkt, MatrixOp.mult(Ckm1, Hk));
- if (_debug)
+ if (_debug) {
System.out.println("Pre Vk: Rk = " + Rk.toString());
+ }
Rk = (BasicMatrix) MatrixOp.add(Rk, Vk);
- if (_debug)
+ if (_debug) {
System.out.println("Post Vk: Rk = " + Rk.toString());
+ }
BasicMatrix Rkinv = (BasicMatrix) MatrixOp.inverse(Rk);
BasicMatrix Kk = (BasicMatrix) MatrixOp.mult(Ckm1, MatrixOp.mult(Hk, Rkinv));
@@ -547,70 +518,138 @@
// if(_debug)System.out.println("Unconstrained chi^2 = "+_chiSq);
//ok...add to the chi^2
- if (_debug)
+ 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 void constrainV0toBS() {
- BasicMatrix Hk = new BasicMatrix(3 * (_ntracks + 1), 3);
- 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)
+// 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);
- //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;
- }
- //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);
-
+// 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, 1);
- Hk.setElement(0, 1, 0);
- Hk.setElement(0, 2, 0);
+ Hk.setElement(0, 0, 0);
+ if (bscon) {
+ Hk.setElement(0, 1, pytot / pxtot);
+ Hk.setElement(0, 2, pztot / pxtot);
+ } else {
+ Hk.setElement(0, 1, 0);
+ Hk.setElement(0, 2, 0);
+ }
Hk.setElement(1, 0, 0);
Hk.setElement(1, 1, 1);
Hk.setElement(1, 2, 0);
@@ -618,92 +657,71 @@
Hk.setElement(2, 1, 0);
Hk.setElement(2, 2, 1);
//ok, loop over tracks again to set the derivitives wrt track momenta (theta,phi,rho)
- for (int i = 0; i < _ntracks; i++) {
+ 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);
+// double px = Pt * Math.cos(phiv);
+// double py = Pt * Math.sin(phiv);
+// double pz = Pt * 1 / Math.tan(theta);
//derivities wrt theta
Hk.setElement(3 * (i + 1), 0, 0);
Hk.setElement(3 * (i + 1), 1, 0);
- Hk.setElement(3 * (i + 1), 2, 0);
+ if (bscon) {
+ Hk.setElement(3 * (i + 1), 2, -Pt / Math.pow(sin(theta), 2) * Vx);
+ } else {
+ Hk.setElement(3 * (i + 1), 2, 0);
+ }
//derivities wrt phi
Hk.setElement(3 * (i + 1) + 1, 0, 0);
- Hk.setElement(3 * (i + 1) + 1, 1,
- 0);
- Hk.setElement(3 * (i + 1) + 1, 2, 0);
+ 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);
+ } else {
+ Hk.setElement(3 * (i + 1) + 1, 1, 0);
+ Hk.setElement(3 * (i + 1) + 1, 2, 0);
+ }
//derivities wrt rho
Hk.setElement(3 * (i + 1) + 2, 0, 0);
// Hk.setElement(3 * (i + 1) + 2, 1,
// (pytot / pxtot - 1) * (Pt / rho) * (1 / pxtot) * Vx);
// Hk.setElement(3 * (i + 1) + 2, 2,
// (pztot / pxtot - 1) * (Pt / rho) * (1 / pxtot) * Vx);
- Hk.setElement(3 * (i + 1) + 2, 1,
- 0);
- Hk.setElement(3 * (i + 1) + 2, 2,
- 0);
+ if (bscon) {
+ Hk.setElement(3 * (i + 1) + 2, 1,
+ (cos(phiv) * pytot / pxtot - sin(phiv)) * (Pt / rho) * (1 / pxtot) * Vx);
+ Hk.setElement(3 * (i + 1) + 2, 2,
+ (cos(phiv) * pztot / pxtot - sin(phiv)) * (Pt / rho) * (1 / pxtot) * Vx);
+ } else {
+ Hk.setElement(3 * (i + 1) + 2, 1, 0);
+ Hk.setElement(3 * (i + 1) + 2, 2, 0);
+ }
// if(_debug)System.out.println("pxtot = "+pxtot+"; rho = "+rho+"; Pt = "+Pt);
// if(_debug)System.out.println("cos(phiv)*pytot / pxtot - sin(phiv) = "+(cos(phiv)*pytot / pxtot - sin(phiv)));
// if(_debug)System.out.println("Pt/(rho*pxtot) = "+(Pt / rho) * (1 / pxtot));
}
- // 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);
+ return Hk;
+ }
+
+ public BasicMatrix makeRk(double Vx, double Vy, double Vz, double pxtot, double pytot, double pztot, boolean bscon) {
+ //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);
+ }
+ 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));
+ } else {
+ rk.setElement(0, 0, -Vx);
+ rk.setElement(1, 0, -Vy);
+ rk.setElement(2, 0, -Vz);
+ }
+ return rk;
}
public void setV0(double[] v0) {
@@ -718,13 +736,13 @@
public void doBeamSpotConstraint(boolean bsconst) {
_beamspotConstraint = bsconst;
- _constraintType="BeamspotConstrained";
-
+ _constraintType = "BeamspotConstrained";
+
}
public void doTargetConstraint(boolean bsconst) {
_targetConstraint = bsconst;
- _constraintType="TargetConstrained";
+ _constraintType = "TargetConstrained";
}
public double getChiSq() {
@@ -742,12 +760,12 @@
mom[0] = Pt * Math.cos(phiv);
mom[1] = Pt * Math.sin(phiv);
mom[2] = Pt * 1 / Math.tan(theta);
- if (_debug){
- System.out.println("getFittedMomentum:: "+mom[0] + "; " + mom[1] + "; " + mom[2]);
-
- System.out.println("pT= "+Pt+"; phi = "+phiv+"; B = "+ _bField);
- }
- return mom;
+ if (_debug) {
+ System.out.println("getFittedMomentum:: " + mom[0] + "; " + mom[1] + "; " + mom[2]);
+
+ System.out.println("pT= " + Pt + "; phi = " + phiv + "; B = " + _bField);
+ }
+ return mom;
}
private double getInvMass() {
@@ -772,10 +790,11 @@
double psum = Math.sqrt(pxsum * pxsum + pysum * pysum + pzsum * pzsum);
double evtmass = esum * esum - psum * psum;
- if (evtmass > 0)
+ if (evtmass > 0) {
return Math.sqrt(evtmass);
- else
+ } else {
return -99;
+ }
}
public String toString() {
@@ -815,8 +834,6 @@
double rho = par[4];
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
double vv = v0.e(1, 0) * cos(par[3]) - v0.e(0, 0) * sin(par[3]);//R
@@ -833,10 +850,10 @@
BasicMatrix q0 = new BasicMatrix(3, 1);
/* this looks just wrong...
- q0.setElement(0, 0, Pt * Math.cos(phiv));
- q0.setElement(1, 0, Pt * Math.sin(phiv));
- q0.setElement(2, 0, Pt * 1 / Math.tan(theta));
- q0s.add(q0);
+ q0.setElement(0, 0, Pt * Math.cos(phiv));
+ q0.setElement(1, 0, Pt * Math.sin(phiv));
+ q0.setElement(2, 0, Pt * 1 / Math.tan(theta));
+ q0s.add(q0);
*/
q0.setElement(0, 0, theta);
q0.setElement(1, 0, phiVert);
@@ -875,10 +892,11 @@
BasicMatrix tmpG = (BasicMatrix) MatrixOp.inverse(bt.covariance());
Gs.add(tmpG);
- if (firstTrack)
+ if (firstTrack) {
D0 = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(tmpA), MatrixOp.mult(tmpG, tmpA));
- else
+ } else {
D0 = (BasicMatrix) MatrixOp.add(D0, MatrixOp.mult(MatrixOp.transposed(tmpA), MatrixOp.mult(tmpG, tmpA)));
+ }
BasicMatrix tmpDi = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(tmpA), MatrixOp.mult(tmpG, tmpB));
BasicMatrix tmpEi = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(tmpB), MatrixOp.mult(tmpG, tmpB));
@@ -904,15 +922,17 @@
BasicMatrix beIbTg = (BasicMatrix) MatrixOp.mult(b, MatrixOp.mult(MatrixOp.inverse(e), MatrixOp.mult(MatrixOp.transposed(b), g)));
BasicMatrix MinusaTgbeIbTg = (BasicMatrix) MatrixOp.mult(-1, MatrixOp.mult(aTg, beIbTg));
- if (firstTrack)
+ if (firstTrack) {
bigsum = (BasicMatrix) MatrixOp.mult(MatrixOp.add(aTg, MinusaTgbeIbTg), p);
- else
+ } else {
bigsum = (BasicMatrix) MatrixOp.add(bigsum, MatrixOp.mult(MatrixOp.add(aTg, MinusaTgbeIbTg), p));
+ }
}
BasicMatrix covVtx = (BasicMatrix) MatrixOp.inverse(tmpCovVtx);
BasicMatrix xtilde = (BasicMatrix) MatrixOp.mult(covVtx, bigsum);
- if (_debug)
+ if (_debug) {
System.out.println("follow1985Paper::Vertex at : \nx= " + xtilde.e(0, 0) + " +/- " + Math.sqrt(covVtx.e(0, 0)) + "\ny= " + xtilde.e(1, 0) + " +/- " + Math.sqrt(covVtx.e(1, 1)) + "\nz= " + xtilde.e(2, 0) + " +/- " + Math.sqrt(covVtx.e(2, 2)));
+ }
//ok, now the momentum
List<Matrix> qtildes = new ArrayList<Matrix>();
@@ -938,18 +958,24 @@
BasicMatrix ptilde = (BasicMatrix) MatrixOp.add(MatrixOp.mult(a, xtilde), MatrixOp.mult(b, qtilde));
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)
+ if (_debug) {
System.out.println("\n\nfollow1985Paper::Track #" + j);
- if (_debug)
+ }
+ if (_debug) {
System.out.println("eps(meas) = " + p.e(0, 0) + " eps(fit) =" + ptilde.e(0, 0));
- if (_debug)
+ }
+ if (_debug) {
System.out.println("zp(meas) = " + p.e(1, 0) + " zp(fit) =" + ptilde.e(1, 0));
- if (_debug)
+ }
+ if (_debug) {
System.out.println("theta(meas) = " + p.e(2, 0) + " theta(fit) =" + ptilde.e(2, 0));
- if (_debug)
+ }
+ if (_debug) {
System.out.println("phi(meas) = " + p.e(3, 0) + " phi(fit) =" + ptilde.e(3, 0));
- if (_debug)
+ }
+ if (_debug) {
System.out.println("rho(meas) = " + p.e(4, 0) + " rho(fit) =" + ptilde.e(4, 0));
+ }
BasicMatrix tmpC0j = (BasicMatrix) MatrixOp.mult(-1, MatrixOp.mult(covVtx, MatrixOp.mult(d, MatrixOp.inverse(e))));
C0j.add(tmpC0j);
@@ -970,8 +996,9 @@
pfit.add(tmppfit);
}
- if (_debug)
+ if (_debug) {
System.out.println("follow1985Paper::chi^2 = " + chisq);
+ }
_chiSq = chisq;
_covVtx = covVtx;
|
