1 added + 10 modified, total 11 files
hps-java/src/main/java/org/lcsim/hps/users/mgraham
diff -u -r1.4 -r1.5
--- DetailedAnalysisDriver.java 19 Sep 2012 22:56:46 -0000 1.4
+++ DetailedAnalysisDriver.java 12 Mar 2013 19:40:11 -0000 1.5
@@ -48,6 +48,7 @@
import org.lcsim.hps.recon.vertexing.StraightLineTrack;
import org.lcsim.hps.recon.tracking.FindableTrack;
import org.lcsim.hps.recon.tracking.TrackAnalysis;
+import org.lcsim.hps.recon.vertexing.BilliorVertexer;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
@@ -101,7 +102,7 @@
public String outputPlots = "myplots.aida";
Map<String, IProfile1D> clsizeMap = new HashMap<String, IProfile1D>();
String[] detNames = {"Tracker"};
- Integer[] nlayers = {10};
+ Integer[] nlayers = {12};
int trk_count = 0;
int nevt = 0;
int _nmcTrk = 0;
@@ -187,15 +188,14 @@
}
}
-
-
@Override
public void detectorChanged(Detector detector) {
// Setup default stereo pairings, which should work for even number of
// modules.
nlayers[0] = detector.getSubdetector("Tracker").getLayering().getLayers().getNumberOfLayers();
- System.out.println("Setting nlayers to "+nlayers[0]);
+ System.out.println("Setting nlayers to " + nlayers[0]);
}
+
@Override
public void process(
EventHeader event) {
@@ -298,7 +298,7 @@
// List<HelicalTrackHit> hthits = event.get(HelicalTrackHit.class, "MatchedHTHits");
-// List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
+ List<HelicalTrackHit> NonRotatedHits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "RotatedHelicalTrackHits");
// List<HelicalTrackHit> axialhits = event.get(HelicalTrackHit.class, "AxialTrackHits");
@@ -331,13 +331,20 @@
for (SiTrackerHitStrip1D stripCluster : stripHits) {
Hep3Vector strCluPos = stripCluster.getPositionAsVector();
+ stripCluster.getMeasuredCoordinate();
+ Hep3Vector measVec = stripCluster.getMeasuredCoordinate();
+ double umeas = VecOp.dot(strCluPos, measVec);
+
double yHit = strCluPos.y();
Set<MCParticle> mcparts = stripCluster.getMCParticles();
aida.cloud1D(occDir + "associated MC Particles").fill(mcparts.size());
List<RawTrackerHit> rthList = stripCluster.getRawHits();
int nhits = rthList.size();
String detlayer = "Foobar";
+ Hep3Vector sthPos = new BasicHep3Vector();
for (RawTrackerHit rth : rthList) {
+ if (rth.getSimTrackerHits().size() > 0)
+ sthPos = rth.getSimTrackerHits().get(0).getPositionVec();
IDetectorElement rhDetE = rth.getDetectorElement();
String rhDetName = rhDetE.getName();
int rhLayer = rth.getLayerNumber();
@@ -345,6 +352,8 @@
if (rhDetName.contains(myname))
detlayer = myname + "_layer" + rhLayer;
}
+ Hep3Vector resid = VecOp.sub(strCluPos, sthPos);
+ double magResid = resid.magnitude();
clsizeMap.get(detlayer).fill(yHit, nhits);
aida.cloud1D(occDir + detlayer + "associated MC Particles").fill(mcparts.size());
aida.cloud1D(occDir + detlayer + " cluster size").fill(nhits);
@@ -748,7 +757,8 @@
// trktomcAxial.add(track, tkanal.getMCParticle());
//
// }
- for (HelicalTrackHit hit : toththits) {
+// for (HelicalTrackHit hit : toththits) {
+ for (HelicalTrackHit hit : NonRotatedHits) {
int nAssHits = hit.getRawHits().size();
aida.cloud1D(debugDir + hit.Detector() + " nAssHits").fill(nAssHits);
@@ -775,6 +785,7 @@
String stripdir = "axial";
double umeas = cluster.umeas();
double charge = cluster.dEdx() * 1000.0;
+ int layer = cluster.layer();
for (RawTrackerHit rhit : rawhits) {
String deName = rhit.getDetectorElement().getName();
@@ -806,19 +817,19 @@
du_stereo = umeas - umc;
if (stripdir.contains("axial"))
du_axial = umeas - umc;
- aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)").fill(umeas - umc);
- aida.cloud1D(debugDir + hit.Detector() + " delta(u)").fill(umeas - umc);
+ aida.histogram1D(debugDir + "layer=" + layer + " delta(u)", 50, -0.03, 0.03).fill(umeas - umc);
+ aida.histogram1D(debugDir + hit.Detector() + " delta(u)", 50, -0.03, 0.03).fill(umeas - umc);
if (nstrips == 1) {
- aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)--1 strip").fill(umeas - umc);
- aida.cloud1D(debugDir + hit.Detector() + " delta(u)--1 strip").fill(umeas - umc);
+ aida.histogram1D(debugDir + "layer=" + layer + " delta(u)--1 strip", 50, -0.03, 0.03).fill(umeas - umc);
+ aida.histogram1D(debugDir + hit.Detector() + " delta(u)--1 strip", 50, -0.03, 0.03).fill(umeas - umc);
}
if (nstrips == 2) {
- aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)--2 strip").fill(umeas - umc);
- aida.cloud1D(debugDir + hit.Detector() + " delta(u)--2 strip").fill(umeas - umc);
+ aida.histogram1D(debugDir + "layer=" + layer + " delta(u)--2 strip", 50, -0.03, 0.03).fill(umeas - umc);
+ aida.histogram1D(debugDir + hit.Detector() + " delta(u)--2 strip", 50, -0.03, 0.03).fill(umeas - umc);
}
if (nstrips == 3) {
- aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)--3 strip").fill(umeas - umc);
- aida.cloud1D(debugDir + hit.Detector() + " delta(u)--3 strip").fill(umeas - umc);
+ aida.histogram1D(debugDir + "layer=" + layer + " delta(u)--3 strip", 50, -0.03, 0.03).fill(umeas - umc);
+ aida.histogram1D(debugDir + hit.Detector() + " delta(u)--3 strip", 50, -0.03, 0.03).fill(umeas - umc);
}
}
@@ -902,64 +913,50 @@
- BilliorVertex bvertexUC = new BilliorVertex(bfield);
- bvertexUC.doBeamSpotConstraint(false);
- bvertexUC.tryNewFormalism(btlist);
-
- BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getVertexPosition();
- BasicMatrix bvtxCovUC = (BasicMatrix) bvertexUC.getVertexCovariance();
- double invMass = bvertexUC.getInvMass();
-
- aida.histogram1D(vertex + "BilliorVertex X -- UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0));
- aida.histogram1D(vertex + "BilliorVertex Y -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0));
- aida.histogram1D(vertex + "BilliorVertex Z -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0));
- aida.histogram1D(vertex + "BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChiSq());
- aida.histogram1D(vertex + "BilliorVertex X Pull -- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) / Math.sqrt(bvtxCovUC.e(0, 0)));
- aida.histogram1D(vertex + "BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1)));
- aida.histogram1D(vertex + "BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2)));
- aida.histogram1D(vertex + "BilliorVertex Mass -- UnConstrained", 250, 0.0, 0.25).fill(bvertexUC.getInvMass());
+ BilliorVertexer bvertexerUC = new BilliorVertexer(bfield);
+ BilliorVertex bvertexUC = bvertexerUC.fitVertex(btlist);
+// bvertexUC.fitVertex(btlist);
+ BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getPosition();
+ SymmetricMatrix bvtxCovUC = bvertexUC.getCovMatrix();
+ double invMassUC = bvertexUC.getParameters().get("invMass");
+// System.out.println("UnConstrained");
+// System.out.println("Vertex Position: " + bvtxPosUC.toString());
+// System.out.println("chisq : " + bvertexUC.getChiSq());
+ aida.histogram1D("BilliorVertex X -- UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0));
+ aida.histogram1D("BilliorVertex Y -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0));
+ aida.histogram1D("BilliorVertex Z -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0));
+ aida.histogram1D("BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChi2());
+ aida.histogram1D("BilliorVertex X Pull -- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) / Math.sqrt(bvtxCovUC.e(0, 0)));
+ aida.histogram1D("BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1)));
+ aida.histogram1D("BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2)));
- aida.cloud1D(apdir + "e+e- Invariant Mass").fill(invMass);
+ aida.cloud1D(apdir + "e+e- Invariant Mass").fill(invMassUC);
if (eleMC != null && posMC != null && ele == eleID && pos == posID)
- aida.cloud1D(apdir + "Matched A' Invariant Mass").fill(invMass);
+ aida.cloud1D(apdir + "Matched A' Invariant Mass").fill(invMassUC);
+
+ BilliorVertexer bconvertexer = new BilliorVertexer(bfield);
+ bconvertexer.setBeamSize(beamsize);
+ bconvertexer.constrainV0toBeamSpot(true);
- BilliorVertex bvertex = new BilliorVertex(bfield);
- bvertex.doBeamSpotConstraint(true);
- bvertex.tryNewFormalism(btlist);
- BasicMatrix bvtxPos = (BasicMatrix) bvertex.getVertexPosition();
- BasicMatrix bvtxCov = (BasicMatrix) bvertex.getVertexCovariance();
-
- aida.histogram1D(vertex + "BilliorVertex X -- Constrained", 100, -10, 20).fill(bvtxPos.e(0, 0));
- aida.histogram1D(vertex + "BilliorVertex Y -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(1, 0));
- aida.histogram1D(vertex + "BilliorVertex Z -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(2, 0));
- aida.histogram1D(vertex + "BilliorVertex ChiSq -- Constrained", 100, -10, 50).fill(bvertex.getChiSq());
- aida.histogram1D(vertex + "BilliorVertex X Pull -- Constrained", 100, -4, 4).fill(bvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
- aida.histogram1D(vertex + "BilliorVertex Y Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
- aida.histogram1D(vertex + "BilliorVertex Z Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
-
- aida.histogram1D(vertex + "BilliorVertex Mass -- Constrained", 250, 0.0, 0.25).fill(bvertex.getInvMass());
+ BilliorVertex bsconfit = bconvertexer.fitVertex(btlist);
+// bvertexUC.fitVertex(btlist);
+ BasicMatrix bsconvtxPos = (BasicMatrix) bsconfit.getPosition();
+ SymmetricMatrix bsconvtxCov = bsconfit.getCovMatrix();
+ double invMassBSCon = bsconfit.getParameters().get("invMass");
+ aida.histogram1D("BilliorVertex X -- BS Constrained", 100, -10, 20).fill(bsconvtxPos.e(0, 0));
+ aida.histogram1D("BilliorVertex Y -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(1, 0));
+ aida.histogram1D("BilliorVertex Z -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(2, 0));
+ aida.histogram1D("BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChi2());
+ aida.histogram1D("BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bsconvtxCov.e(0, 0)));
+ aida.histogram1D("BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bsconvtxCov.e(1, 1)));
+ aida.histogram1D("BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bsconvtxCov.e(2, 2)));
- BilliorVertex bsconfit = new BilliorVertex(bfield);
- bsconfit.setBeamSize(beamsize);
- bsconfit.doBeamSpotConstraint(false);
- bsconfit.constrainV0toBeamSpot(true);
- bsconfit.tryNewFormalism(btlist);
- BasicMatrix bsconvtxPos = (BasicMatrix) bsconfit.getVertexPosition();
- BasicMatrix bsconvtxCov = (BasicMatrix) bsconfit.getVertexCovariance();
-
- aida.histogram1D(vertex + "BilliorVertex X -- BS Constrained", 100, -10, 20).fill(bsconvtxPos.e(0, 0));
- aida.histogram1D(vertex + "BilliorVertex Y -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(1, 0));
- aida.histogram1D(vertex + "BilliorVertex Z -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(2, 0));
- aida.histogram1D(vertex + "BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChiSq());
- aida.histogram1D(vertex + "BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
- aida.histogram1D(vertex + "BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
- aida.histogram1D(vertex + "BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
+ aida.histogram1D("BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(invMassBSCon);
- aida.histogram1D(vertex + "BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(bsconfit.getInvMass());
}
}
}
hps-java/src/main/java/org/lcsim/hps/users/mgraham
diff -u -r1.8 -r1.9
--- FastTrackAnalysisDriver.java 19 Sep 2012 22:56:46 -0000 1.8
+++ FastTrackAnalysisDriver.java 12 Mar 2013 19:40:11 -0000 1.9
@@ -2,6 +2,7 @@
import hep.aida.IAnalysisFactory;
import hep.physics.matrix.BasicMatrix;
+import hep.physics.matrix.SymmetricMatrix;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
@@ -37,6 +38,7 @@
import org.lcsim.hps.recon.vertexing.StraightLineTrack;
import org.lcsim.hps.recon.tracking.FindableTrack;
import org.lcsim.hps.recon.tracking.TrackAnalysis;
+import org.lcsim.hps.recon.vertexing.BilliorVertexer;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
@@ -366,100 +368,57 @@
- BilliorVertex vtxfit = new BilliorVertex(bfield);
- vtxfit.doBeamSpotConstraint(false);
- vtxfit.setBeamSize(beamsize);
-// vtxfit.fitVertex(btlist);
- vtxfit.tryNewFormalism(btlist);
- double[] momE = vtxfit.getFittedMomentum(0);
- double[] momP = vtxfit.getFittedMomentum(1);
- double vtxpxE = momE[0];
- double vtxpyE = momE[1];
- double vtxpzE = momE[2];
- double vtxpxP = momP[0];
- double vtxpyP = momP[1];
- double vtxpzP = momP[2];
- double chisq = vtxfit.getChiSq();
- BasicMatrix vtx = (BasicMatrix) vtxfit.getVertexPosition();
- BasicMatrix vtxcov = (BasicMatrix) vtxfit.getVertexCovariance();
+ BilliorVertexer vtxfitter = new BilliorVertexer(bfield);
+ vtxfitter.doBeamSpotConstraint(false);
+ vtxfitter.setBeamSize(beamsize);
+
+ BilliorVertexer vtxfitterCon = new BilliorVertexer(bfield);
+ vtxfitterCon.doBeamSpotConstraint(true);
+ vtxfitterCon.setBeamSize(beamsize);
+
+ BilliorVertexer vtxfitterBSCon = new BilliorVertexer(bfield);
+ vtxfitterBSCon.constrainV0toBeamSpot(true);
+ vtxfitterBSCon.setBeamSize(beamsize);
+
+
+ BilliorVertex vtxfit = vtxfitter.fitVertex(btlist);
+ Map<String, Double> vtxMap = vtxfit.getParameters();
+
+ double vtxpxE = vtxMap.get("p1X");
+ double vtxpyE = vtxMap.get("p1Y");
+ double vtxpzE = vtxMap.get("p1Z");
+ double vtxpxP = vtxMap.get("p2X");
+ double vtxpyP = vtxMap.get("p2Y");
+ double vtxpzP = vtxMap.get("p2Z");
+ double chisq = vtxfit.getChi2();
+ BasicMatrix vtx = (BasicMatrix) vtxfit.getPosition();
+ SymmetricMatrix vtxcov = vtxfit.getCovMatrix();
double chisqE = ele.getChi2();
double chisqP = pos.getChi2();
- BilliorVertex confit = new BilliorVertex(bfield);
- confit.setBeamSize(beamsize);
- confit.doBeamSpotConstraint(true);
-// confit.fitVertex(btlist);
- confit.tryNewFormalism(btlist);
- double[] conmomE = confit.getFittedMomentum(0);
- double[] conmomP = confit.getFittedMomentum(1);
- double conpxE = conmomE[0];
- double conpyE = conmomE[1];
- double conpzE = conmomE[2];
- double conpxP = conmomP[0];
- double conpyP = conmomP[1];
- double conpzP = conmomP[2];
- double conchisq = confit.getChiSq();
- BasicMatrix conVtx = (BasicMatrix) confit.getVertexPosition();
- BasicMatrix conVtxCov = (BasicMatrix) confit.getVertexCovariance();
-
- // 2nd pass vertex parameters
- BilliorVertex vtxfit2 = new BilliorVertex(bfield);
- vtxfit2.setBeamSize(beamsize);
- double[] vtxfitv0 = {vtx.e(0, 0), vtx.e(1, 0), vtx.e(2, 0)};
- vtxfit2.setV0(vtxfitv0);
- vtxfit2.doBeamSpotConstraint(false);
-// vtxfit.fitVertex(btlist);
- vtxfit2.tryNewFormalism(btlist);
- double[] mom2E = vtxfit2.getFittedMomentum(0);
- double[] mom2P = vtxfit2.getFittedMomentum(1);
- double vtxpx2E = mom2E[0];
- double vtxpy2E = mom2E[1];
- double vtxpz2E = mom2E[2];
- double vtxpx2P = mom2P[0];
- double vtxpy2P = mom2P[1];
- double vtxpz2P = mom2P[2];
- double chisq2 = vtxfit2.getChiSq();
- BasicMatrix vtx2 = (BasicMatrix) vtxfit2.getVertexPosition();
- BasicMatrix vtxcov2 = (BasicMatrix) vtxfit2.getVertexCovariance();
-
- BilliorVertex confit2 = new BilliorVertex(bfield);
- confit2.setBeamSize(beamsize);
-
- double[] confitv0 = {conVtx.e(0, 0), conVtx.e(1, 0), conVtx.e(2, 0)};
- confit2.setV0(confitv0);
- confit2.doBeamSpotConstraint(true);
-// confit.fitVertex(btlist);
- confit2.tryNewFormalism(btlist);
- double[] conmom2E = confit2.getFittedMomentum(0);
- double[] conmom2P = confit2.getFittedMomentum(1);
- double conpx2E = conmom2E[0];
- double conpy2E = conmom2E[1];
- double conpz2E = conmom2E[2];
- double conpx2P = conmom2P[0];
- double conpy2P = conmom2P[1];
- double conpz2P = conmom2P[2];
- double conchisq2 = confit2.getChiSq();
- BasicMatrix conVtx2 = (BasicMatrix) confit2.getVertexPosition();
- BasicMatrix conVtxCov2 = (BasicMatrix) confit2.getVertexCovariance();
-
- BilliorVertex bsconfit = new BilliorVertex(bfield);
- bsconfit.setBeamSize(beamsize);
- bsconfit.doBeamSpotConstraint(false);
- bsconfit.constrainV0toBeamSpot(true);
-// bsconfit.fitVertex(btlist);
- bsconfit.tryNewFormalism(btlist);
- double[] bsconmomE = bsconfit.getFittedMomentum(0);
- double[] bsconmomP = bsconfit.getFittedMomentum(1);
- double bsconpxE = bsconmomE[0];
- double bsconpyE = bsconmomE[1];
- double bsconpzE = bsconmomE[2];
- double bsconpxP = bsconmomP[0];
- double bsconpyP = bsconmomP[1];
- double bsconpzP = bsconmomP[2];
- double bsconchisq = bsconfit.getChiSq();
- BasicMatrix bsconVtx = (BasicMatrix) bsconfit.getVertexPosition();
- BasicMatrix bsconVtxCov = (BasicMatrix) bsconfit.getVertexCovariance();
-
+ BilliorVertex confit = vtxfitterCon.fitVertex(btlist);
+ Map<String, Double> conMap = confit.getParameters();
+ double conpxE = conMap.get("p1X");
+ double conpyE = conMap.get("p1Y");
+ double conpzE = conMap.get("p1Z");
+ double conpxP = conMap.get("p2X");
+ double conpyP = conMap.get("p2Y");
+ double conpzP = conMap.get("p2Z");
+ double conchisq = confit.getChi2();
+ BasicMatrix conVtx = (BasicMatrix) confit.getPosition();
+ SymmetricMatrix conVtxCov = confit.getCovMatrix();
+
+ BilliorVertex bsconfit = vtxfitterBSCon.fitVertex(btlist);
+ Map<String, Double> bsconMap = bsconfit.getParameters();
+ double bsconpxE = bsconMap.get("p1X");
+ double bsconpyE = bsconMap.get("p1Y");
+ double bsconpzE = bsconMap.get("p1Z");
+ double bsconpxP = bsconMap.get("p2X");
+ double bsconpyP = bsconMap.get("p2Y");
+ double bsconpzP = bsconMap.get("p2Z");
+ double bsconchisq = bsconfit.getChi2();
+ BasicMatrix bsconVtx = (BasicMatrix) bsconfit.getPosition();
+ SymmetricMatrix bsconVtxCov = bsconfit.getCovMatrix();
double l1minE = -99;
double l1minP = -99;
@@ -582,23 +541,6 @@
pw.format("%5.5f %5.5f %5.5f %5.5f %5.5f %5.5f ", conVtx.e(0, 0), conErrX, conVtx.e(1, 0), conErrY, conVtx.e(2, 0), conErrZ);
pw.format("%5.5f ", conchisq);
- pw.format("%5.5f %5.5f %5.5f ", vtxpx2E, vtxpy2E, vtxpz2E);
- pw.format("%5.5f %5.5f %5.5f ", vtxpx2P, vtxpy2P, vtxpz2P);
- pw.format("%5.5f %5.5f %5.5f %5.5f %5.5f %5.5f ", vtx2.e(0, 0), Math.sqrt(vtxcov2.e(0, 0)), vtx2.e(1, 0), Math.sqrt(vtxcov2.e(1, 1)), vtx2.e(2, 0), Math.sqrt(vtxcov2.e(2, 2)));
- pw.format("%5.5f ", chisq2);
-
- pw.format("%5.5f %5.5f %5.5f ", conpx2E, conpy2E, conpz2E);
- pw.format("%5.5f %5.5f %5.5f ", conpx2P, conpy2P, conpz2P);
-
- // get the errors on the constrained vertex an make sure they aren't NaN
- // there must be something wrong in the constraint...hopefully just the error calc
- double conErrX2 = getErr(conVtxCov2.e(0, 0));
- double conErrY2 = getErr(conVtxCov2.e(1, 1));
- double conErrZ2 = getErr(conVtxCov2.e(2, 2));
-
- pw.format("%5.5f %5.5f %5.5f %5.5f %5.5f %5.5f ", conVtx2.e(0, 0), conErrX2, conVtx2.e(1, 0), conErrY2, conVtx2.e(2, 0), conErrZ2);
- pw.format("%5.5f ", conchisq2);
-
pw.format("%5.5f %5.5f %5.5f ", bsconpxE, bsconpyE, bsconpzE);
pw.format("%5.5f %5.5f %5.5f ", bsconpxP, bsconpyP, bsconpzP);
hps-java/src/main/java/org/lcsim/hps/users/mgraham
diff -u -r1.4 -r1.5
--- MainJASDriver.java 26 Nov 2012 19:27:39 -0000 1.4
+++ MainJASDriver.java 12 Mar 2013 19:40:11 -0000 1.5
@@ -19,13 +19,16 @@
add(new HPSSVTSensorSetup());
add(new TrackerDigiDriver());
// add(new TrackerReconDriver());
- HelicalTrackHitDriver hth = new HelicalTrackHitDriver();
+// HelicalTrackHitDriver hth = new HelicalTrackHitDriver();
+
+ SingleSensorHelicalTrackHitDriver hth = new SingleSensorHelicalTrackHitDriver();
hth.setMaxSeperation(20.0);
hth.setTolerance(1.0);
add(hth);
TrackerReconDriver trd=new TrackerReconDriver();
- trd.setStrategyResource("/org/lcsim/hps/recon/tracking/strategies/HPS-Full.xml");
+ trd.setStrategyResource("/org/lcsim/hps/recon/tracking/strategies/HPS-SingleSensors.xml");
+// trd.setStrategyResource("/org/lcsim/hps/recon/tracking/strategies/HPS-Full.xml");
// trd.setStrategyResource("/org/lcsim/hps/recon/tracking/strategies/HPS-Test-Lyr50.xml");
add(trd);
add(new DetailedAnalysisDriver(12));
hps-java/src/main/java/org/lcsim/hps/users/mgraham
diff -u -r1.1 -r1.2
--- ExamplePlotter.java 11 Sep 2012 20:30:10 -0000 1.1
+++ ExamplePlotter.java 12 Mar 2013 19:40:11 -0000 1.2
@@ -4,12 +4,22 @@
import hep.aida.IHistogram1D;
import hep.aida.IPlotter;
import hep.aida.IPlotterStyle;
+import hep.physics.vec.Hep3Vector;
+import java.io.IOException;
import java.util.List;
+import java.util.logging.Level;
+import java.util.logging.Logger;
import org.lcsim.event.EventHeader;
import org.lcsim.event.Track;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.geometry.Detector;
import org.lcsim.hps.monitoring.AIDAFrame;
import org.lcsim.hps.monitoring.Resettable;
+import org.lcsim.hps.recon.tracking.HPSTrack;
+import org.lcsim.hps.recon.vertexing.HelixConverter;
+import org.lcsim.hps.recon.vertexing.StraightLineTrack;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.recon.tracking.seedtracker.SeedTrack;
import org.lcsim.util.Driver;
import org.lcsim.util.aida.AIDA;
@@ -24,6 +34,8 @@
IPlotter plotter;
IAnalysisFactory fac = aida.analysisFactory();
private String trackCollectionName = "MatchedTracks";
+ private double zAtConverter = -674.062;//mm
+ private String outputPlots = null;
protected void detectorChanged(Detector detector) {
aida.tree().cd("/");
@@ -35,18 +47,21 @@
IPlotterStyle style = plotter.style();
style.dataStyle().fillStyle().setColor("yellow");
style.dataStyle().errorBarStyle().setVisible(false);
- plotter.createRegions(2, 2);
+ plotter.createRegions(2, 3);
plotterFrame.addPlotter(plotter);
IHistogram1D trkPx = aida.histogram1D("Track Momentum (Px)", 25, -0.25, 0.25);
IHistogram1D trkPy = aida.histogram1D("Track Momentum (Py)", 25, -0.1, 0.1);
IHistogram1D trkPz = aida.histogram1D("Track Momentum (Pz)", 25, 0, 3.5);
IHistogram1D trkChi2 = aida.histogram1D("Track Chi2", 25, 0, 25.0);
-
+ IHistogram1D xAtConvert = aida.histogram1D("X (mm) @ Converter", 50, -50, 50);
+ IHistogram1D yAtConvert = aida.histogram1D("Y (mm) @ Converter", 50, -20, 20);
plotter.region(0).plot(trkPx);
plotter.region(1).plot(trkPy);
plotter.region(2).plot(trkPz);
plotter.region(3).plot(trkChi2);
+ plotter.region(4).plot(xAtConvert);
+ plotter.region(5).plot(yAtConvert);
plotterFrame.pack();
plotterFrame.setVisible(true);
@@ -60,6 +75,17 @@
aida.histogram1D("Track Momentum (Py)").fill(trk.getPZ());
aida.histogram1D("Track Momentum (Pz)").fill(trk.getPX());
aida.histogram1D("Track Chi2").fill(trk.getChi2());
+
+ SeedTrack stEle = (SeedTrack) trk;
+ SeedCandidate seedEle = stEle.getSeedCandidate();
+ HelicalTrackFit ht = seedEle.getHelix();
+ HelixConverter converter = new HelixConverter(0);
+ StraightLineTrack slt = converter.Convert(ht);
+ HPSTrack hpstrack = new HPSTrack(ht);
+ Hep3Vector[] trkatconver = hpstrack.getPositionAtZMap(100, zAtConverter, 1);
+ aida.histogram1D("X (mm) @ Converter").fill(trkatconver[0].x()); // y tracker frame?
+ aida.histogram1D("Y (mm) @ Converter").fill(trkatconver[0].y()); // z tracker frame?
+
}
}
@@ -67,4 +93,19 @@
public void reset() {
throw new UnsupportedOperationException("Not supported yet.");
}
+
+ public void setOutputPlots(String output) {
+ this.outputPlots = output;
+ }
+
+ public void endOfData() {
+ System.out.println("Output");
+ if (outputPlots != null) {
+ try {
+ aida.saveAs(outputPlots);
+ } catch (IOException ex) {
+ Logger.getLogger(ElwinsTrackingRecon.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+ }
}
hps-java/src/main/java/org/lcsim/hps/users/mgraham
diff -u -r1.2 -r1.3
--- JasAnalysisDriver.java 1 Aug 2011 18:25:20 -0000 1.2
+++ JasAnalysisDriver.java 12 Mar 2013 19:40:11 -0000 1.3
@@ -40,6 +40,7 @@
import org.lcsim.hps.recon.vertexing.StraightLineTrack;
import org.lcsim.hps.recon.tracking.FindableTrack;
import org.lcsim.hps.recon.tracking.TrackAnalysis;
+import org.lcsim.hps.recon.vertexing.*;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
@@ -48,8 +49,8 @@
import org.lcsim.util.aida.AIDA;
/**
- *
- * @author partridge
+
+ @author partridge
*/
public class JasAnalysisDriver extends Driver {
@@ -104,7 +105,7 @@
FileWriter fw;
PrintWriter pw;
boolean isBeamConstrain = false;
- double[] beamsize = {0.001, 0.02, 0.02};
+ double[] beamsize = {0.001, 0.02, 0.02};
public JasAnalysisDriver(int layers) {
// Define the efficiency histograms
@@ -405,17 +406,15 @@
for (Track track2 : tracklist) {
if (track1 != track2 && track1.getCharge() > 0 && track2.getCharge() < 0) {
/*
- vlist.clear();
- vlist.add(track1);
- vlist.add(track2);
- Vertex vtx = bFit.fit(vlist, SP, isBeamConstrain);
- System.out.println(vtx.toString());
- double vtxChi2 = vtx._chi2;
- double[] vtxPos = vtx._xyzf;
- aida.histogram1D("Vertex Chi2", 100, 0, 1000).fill(vtxChi2);
- aida.histogram1D("Vertex X", 100, -20, 50).fill(vtxPos[0]);
- aida.histogram1D("Vertex Y", 100, -1, 1).fill(vtxPos[1]);
- aida.histogram1D("Vertex Z", 100, -1, 1).fill(vtxPos[2]);
+ vlist.clear(); vlist.add(track1); vlist.add(track2); Vertex
+ vtx = bFit.fit(vlist, SP, isBeamConstrain);
+ System.out.println(vtx.toString()); double vtxChi2 =
+ vtx._chi2; double[] vtxPos = vtx._xyzf;
+ aida.histogram1D("Vertex Chi2", 100, 0,
+ 1000).fill(vtxChi2); aida.histogram1D("Vertex X", 100, -20,
+ 50).fill(vtxPos[0]); aida.histogram1D("Vertex Y", 100, -1,
+ 1).fill(vtxPos[1]); aida.histogram1D("Vertex Z", 100, -1,
+ 1).fill(vtxPos[2]);
*/
@@ -425,115 +424,120 @@
btlist.add(bt1);
btlist.add(bt2);
/*
- BilliorVertex bvertex = new BilliorVertex(1.0);
- bvertex.fitVertex(btlist);
- BasicMatrix bvtxPos = (BasicMatrix) bvertex.getVertexPosition();
- BasicMatrix bvtxCov = (BasicMatrix) bvertex.getVertexCovariance();
- System.out.println("Constrained");
- System.out.println("Vertex Position: " + bvtxPos.toString());
- System.out.println("chisq : " + bvertex.getChiSq());
-
- aida.histogram1D("BilliorVertex X -- Constrained", 100, -10, 20).fill(bvtxPos.e(0, 0));
- aida.histogram1D("BilliorVertex Y -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(1, 0));
- aida.histogram1D("BilliorVertex Z -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(2, 0));
- aida.histogram1D("BilliorVertex ChiSq -- Constrained", 100, 0, 50).fill(bvertex.getChiSq());
- aida.histogram1D("BilliorVertex X Pull -- Constrained", 100, -4, 4).fill(bvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
- aida.histogram1D("BilliorVertex Y Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
- aida.histogram1D("BilliorVertex Z Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
+ BilliorVertex bvertex = new BilliorVertex(1.0);
+ bvertex.fitVertex(btlist); BasicMatrix bvtxPos =
+ (BasicMatrix) bvertex.getVertexPosition(); BasicMatrix
+ bvtxCov = (BasicMatrix) bvertex.getVertexCovariance();
+ System.out.println("Constrained");
+ System.out.println("Vertex Position: " +
+ bvtxPos.toString()); System.out.println("chisq : " +
+ bvertex.getChiSq());
+
+ aida.histogram1D("BilliorVertex X -- Constrained", 100,
+ -10, 20).fill(bvtxPos.e(0, 0));
+ aida.histogram1D("BilliorVertex Y -- Constrained", 100,
+ -0.4, 0.4).fill(bvtxPos.e(1, 0));
+ aida.histogram1D("BilliorVertex Z -- Constrained", 100,
+ -0.4, 0.4).fill(bvtxPos.e(2, 0));
+ aida.histogram1D("BilliorVertex ChiSq -- Constrained", 100,
+ 0, 50).fill(bvertex.getChiSq());
+ aida.histogram1D("BilliorVertex X Pull -- Constrained",
+ 100, -4, 4).fill(bvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0,
+ 0))); aida.histogram1D("BilliorVertex Y Pull--
+ Constrained", 100, -4, 4).fill(bvtxPos.e(1, 0) /
+ Math.sqrt(bvtxCov.e(1, 1)));
+ aida.histogram1D("BilliorVertex Z Pull-- Constrained", 100,
+ -4, 4).fill(bvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
*/
/*
- BilliorVertex bvertexUC = new BilliorVertex(1.0);
- bvertexUC.doBeamSpotConstraint(false);
- bvertexUC.fitVertex(btlist);
- BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getVertexPosition();
- BasicMatrix bvtxCovUC = (BasicMatrix) bvertexUC.getVertexCovariance();
- System.out.println("UnConstrained");
- System.out.println("Vertex Position: " + bvtxPosUC.toString());
- System.out.println("chisq : " + bvertexUC.getChiSq());
- aida.histogram1D("BilliorVertex X -- UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0));
- aida.histogram1D("BilliorVertex Y -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0));
- aida.histogram1D("BilliorVertex Z -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0));
- aida.histogram1D("BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChiSq());
- aida.histogram1D("BilliorVertex X Pull -- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) / Math.sqrt(bvtxCovUC.e(0, 0)));
- aida.histogram1D("BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1)));
- aida.histogram1D("BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2)));
+ BilliorVertex bvertexUC = new BilliorVertex(1.0);
+ bvertexUC.doBeamSpotConstraint(false);
+ bvertexUC.fitVertex(btlist); BasicMatrix bvtxPosUC =
+ (BasicMatrix) bvertexUC.getVertexPosition(); BasicMatrix
+ bvtxCovUC = (BasicMatrix) bvertexUC.getVertexCovariance();
+ System.out.println("UnConstrained");
+ System.out.println("Vertex Position: " +
+ bvtxPosUC.toString()); System.out.println("chisq : " +
+ bvertexUC.getChiSq()); aida.histogram1D("BilliorVertex X --
+ UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0));
+ aida.histogram1D("BilliorVertex Y -- UnConstrained", 100,
+ -0.4, 0.4).fill(bvtxPosUC.e(1, 0));
+ aida.histogram1D("BilliorVertex Z -- UnConstrained", 100,
+ -0.4, 0.4).fill(bvtxPosUC.e(2, 0));
+ aida.histogram1D("BilliorVertex ChiSq -- UnConstrained",
+ 100, 0, 50).fill(bvertexUC.getChiSq());
+ aida.histogram1D("BilliorVertex X Pull -- UnConstrained",
+ 100, -4, 4).fill(bvtxPosUC.e(0, 0) /
+ Math.sqrt(bvtxCovUC.e(0, 0)));
+ aida.histogram1D("BilliorVertex Y Pull-- UnConstrained",
+ 100, -4, 4).fill(bvtxPosUC.e(1, 0) /
+ Math.sqrt(bvtxCovUC.e(1, 1)));
+ aida.histogram1D("BilliorVertex Z Pull-- UnConstrained",
+ 100, -4, 4).fill(bvtxPosUC.e(2, 0) /
+ Math.sqrt(bvtxCovUC.e(2, 2)));
*/
- BilliorVertex bvertexUC = new BilliorVertex(bfield);
- bvertexUC.doBeamSpotConstraint(false);
- bvertexUC.tryNewFormalism(btlist);
+ BilliorVertexer bvertexerUC = new BilliorVertexer(bfield);
+ BilliorVertex bvertexUC = bvertexerUC.fitVertex(btlist);
// bvertexUC.fitVertex(btlist);
- BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getVertexPosition();
- BasicMatrix bvtxCovUC = (BasicMatrix) bvertexUC.getVertexCovariance();
+ BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getPosition();
+ SymmetricMatrix bvtxCovUC = bvertexUC.getCovMatrix();
+ double invMassUC = bvertexUC.getParameters().get("invMass");
// System.out.println("UnConstrained");
// System.out.println("Vertex Position: " + bvtxPosUC.toString());
// System.out.println("chisq : " + bvertexUC.getChiSq());
aida.histogram1D("BilliorVertex X -- UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0));
aida.histogram1D("BilliorVertex Y -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0));
aida.histogram1D("BilliorVertex Z -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0));
- aida.histogram1D("BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChiSq());
+ aida.histogram1D("BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChi2());
aida.histogram1D("BilliorVertex X Pull -- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) / Math.sqrt(bvtxCovUC.e(0, 0)));
aida.histogram1D("BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1)));
aida.histogram1D("BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2)));
- BilliorVertex bvertex2ndPass = new BilliorVertex(bfield);
- bvertex2ndPass.doBeamSpotConstraint(false);
- double[] v0 = {bvtxPosUC.e(0, 0), bvtxPosUC.e(1, 0), bvtxPosUC.e(2, 0)};
- bvertex2ndPass.setV0(v0);
- bvertex2ndPass.tryNewFormalism(btlist);
-// bvertex2ndPass.fitVertex(btlist);
- BasicMatrix bvtxPos2ndPass = (BasicMatrix) bvertex2ndPass.getVertexPosition();
- BasicMatrix bvtxCov2ndPass = (BasicMatrix) bvertex2ndPass.getVertexCovariance();
-// System.out.println("2ndPass");
-// System.out.println("Vertex Position: " + bvtxPos2ndPass.toString());
-// System.out.println("chisq : " + bvertex2ndPass.getChiSq());
- aida.histogram1D("BilliorVertex X -- 2ndPass", 100, -10, 20).fill(bvtxPos2ndPass.e(0, 0));
- aida.histogram1D("BilliorVertex Y -- 2ndPass", 100, -0.4, 0.4).fill(bvtxPos2ndPass.e(1, 0));
- aida.histogram1D("BilliorVertex Z -- 2ndPass", 100, -0.4, 0.4).fill(bvtxPos2ndPass.e(2, 0));
- aida.histogram1D("BilliorVertex ChiSq -- 2ndPass", 100, 0, 50).fill(bvertex2ndPass.getChiSq());
- aida.histogram1D("BilliorVertex X Pull -- 2ndPass", 100, -4, 4).fill(bvtxPos2ndPass.e(0, 0) / Math.sqrt(bvtxCov2ndPass.e(0, 0)));
- aida.histogram1D("BilliorVertex Y Pull-- 2ndPass", 100, -4, 4).fill(bvtxPos2ndPass.e(1, 0) / Math.sqrt(bvtxCov2ndPass.e(1, 1)));
- aida.histogram1D("BilliorVertex Z Pull-- 2ndPass", 100, -4, 4).fill(bvtxPos2ndPass.e(2, 0) / Math.sqrt(bvtxCov2ndPass.e(2, 2)));
-
-
- BilliorVertex bvertex = new BilliorVertex(bfield);
- bvertex.doBeamSpotConstraint(true);
- bvertex.tryNewFormalism(btlist);
-// bvertex.fitVertex(btlist);
- BasicMatrix bvtxPos = (BasicMatrix) bvertex.getVertexPosition();
- BasicMatrix bvtxCov = (BasicMatrix) bvertex.getVertexCovariance();
+
+ BilliorVertexer bvertexer = new BilliorVertexer(bfield);
+ bvertexer.setBeamSize(beamsize);
+ bvertexer.doBeamSpotConstraint(true);
+
+ BilliorVertex bvertex = bvertexer.fitVertex(btlist);
+ BasicMatrix bvtxPos = (BasicMatrix) bvertex.getPosition();
+ SymmetricMatrix bvtxCov = bvertex.getCovMatrix();
+ double invMass = bvertex.getParameters().get("invMass");
+
// System.out.println("Constrained");
// System.out.println("Vertex Position: " + bvtxPos.toString());
// System.out.println("chisq : " + bvertex.getChiSq());
aida.histogram1D("BilliorVertex X -- Constrained", 100, -10, 20).fill(bvtxPos.e(0, 0));
aida.histogram1D("BilliorVertex Y -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(1, 0));
aida.histogram1D("BilliorVertex Z -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(2, 0));
- aida.histogram1D("BilliorVertex ChiSq -- Constrained", 100, -10, 50).fill(bvertex.getChiSq());
+ aida.histogram1D("BilliorVertex ChiSq -- Constrained", 100, -10, 50).fill(bvertex.getChi2());
aida.histogram1D("BilliorVertex X Pull -- Constrained", 100, -4, 4).fill(bvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
aida.histogram1D("BilliorVertex Y Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
aida.histogram1D("BilliorVertex Z Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
- aida.histogram1D("BilliorVertex Mass -- Constrained", 250, 0.0, 0.25).fill(bvertex.getInvMass());
- aida.histogram1D("BilliorVertex Mass -- UnConstrained", 250, 0.0, 0.25).fill(bvertexUC.getInvMass());
- aida.histogram1D("BilliorVertex Mass -- 2nd Pass", 250, 0.0, 0.25).fill(bvertex2ndPass.getInvMass());
-
-
- BilliorVertex bsconfit = new BilliorVertex(bfield);
- bsconfit.setBeamSize(beamsize);
- bsconfit.doBeamSpotConstraint(false);
- bsconfit.constrainV0toBeamSpot(true);
- bsconfit.tryNewFormalism(btlist);
- BasicMatrix bsconvtxPos = (BasicMatrix) bsconfit.getVertexPosition();
- BasicMatrix bsconvtxCov = (BasicMatrix) bsconfit.getVertexCovariance();
+ aida.histogram1D("BilliorVertex Mass -- Constrained", 250, 0.0, 0.25).fill(invMass);
+ aida.histogram1D("BilliorVertex Mass -- UnConstrained", 250, 0.0, 0.25).fill(invMassUC);
+
+
+ BilliorVertexer bconvertexer = new BilliorVertexer(bfield);
+ bconvertexer.setBeamSize(beamsize);
+ bconvertexer.constrainV0toBeamSpot(true);
+
+ BilliorVertex bsconfit = bconvertexer.fitVertex(btlist);
+// bvertexUC.fitVertex(btlist);
+ BasicMatrix bsconvtxPos = (BasicMatrix) bsconfit.getPosition();
+ SymmetricMatrix bsconvtxCov = bsconfit.getCovMatrix();
+ double invMassBSCon = bsconfit.getParameters().get("invMass");
+
aida.histogram1D("BilliorVertex X -- BS Constrained", 100, -10, 20).fill(bsconvtxPos.e(0, 0));
aida.histogram1D("BilliorVertex Y -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(1, 0));
aida.histogram1D("BilliorVertex Z -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(2, 0));
- aida.histogram1D("BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChiSq());
- aida.histogram1D("BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
- aida.histogram1D("BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
- aida.histogram1D("BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
+ aida.histogram1D("BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChi2());
+ aida.histogram1D("BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bsconvtxCov.e(0, 0)));
+ aida.histogram1D("BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bsconvtxCov.e(1, 1)));
+ aida.histogram1D("BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bsconvtxCov.e(2, 2)));
- aida.histogram1D("BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(bsconfit.getInvMass());
+ aida.histogram1D("BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(invMassBSCon);
}
}
}
@@ -587,7 +591,7 @@
double theta = 180. * Math.acos(cth) / Math.PI;
double eta = -Math.log(Math.tan(Math.atan2(pt, pz) / 2));
double phi = Math.atan2(py, px);
- int pdgid=mcp.getPDGID();
+ int pdgid = mcp.getPDGID();
// Find the number of layers hit by this mc particle
// System.out.println("MC pt=" + pt);
int nhits = findable.LayersHit(mcp);
@@ -700,10 +704,10 @@
ctheffElectrons.fill(cth, wgt);
d0effElectrons.fill(d0, wgt);
z0effElectrons.fill(z0, wgt);
- if (wgt == 0&&pdgid>0) {
+ if (wgt == 0 && pdgid > 0) {
System.out.println("Missed a findable ELECTRON!!!!!");
}
- if (wgt == 0&&pdgid<0) {
+ if (wgt == 0 && pdgid < 0) {
System.out.println("Missed a findable POSITRON!!!!!");
}
hps-java/src/main/java/org/lcsim/hps/examples
diff -u -r1.2 -r1.3
--- JasAnalysisDriver.java 1 Aug 2011 18:25:21 -0000 1.2
+++ JasAnalysisDriver.java 12 Mar 2013 19:40:11 -0000 1.3
@@ -40,6 +40,7 @@
import org.lcsim.hps.recon.vertexing.StraightLineTrack;
import org.lcsim.hps.recon.tracking.FindableTrack;
import org.lcsim.hps.recon.tracking.TrackAnalysis;
+import org.lcsim.hps.recon.vertexing.*;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
@@ -469,82 +470,68 @@
aida.histogram1D("BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1)));
aida.histogram1D("BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2)));
*/
- BilliorVertex bvertexUC = new BilliorVertex(bfield);
- bvertexUC.doBeamSpotConstraint(false);
- bvertexUC.tryNewFormalism(btlist);
+ BilliorVertexer bvertexerUC = new BilliorVertexer(bfield);
+ BilliorVertex bvertexUC=bvertexerUC.fitVertex(btlist);
// bvertexUC.fitVertex(btlist);
- BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getVertexPosition();
- BasicMatrix bvtxCovUC = (BasicMatrix) bvertexUC.getVertexCovariance();
+ BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getPosition();
+ SymmetricMatrix bvtxCovUC = bvertexUC.getCovMatrix();
+ double invMassUC=bvertexUC.getParameters().get("invMass");
// System.out.println("UnConstrained");
// System.out.println("Vertex Position: " + bvtxPosUC.toString());
// System.out.println("chisq : " + bvertexUC.getChiSq());
aida.histogram1D("BilliorVertex X -- UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0));
aida.histogram1D("BilliorVertex Y -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0));
aida.histogram1D("BilliorVertex Z -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0));
- aida.histogram1D("BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChiSq());
+ aida.histogram1D("BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChi2());
aida.histogram1D("BilliorVertex X Pull -- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) / Math.sqrt(bvtxCovUC.e(0, 0)));
aida.histogram1D("BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1)));
aida.histogram1D("BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2)));
- BilliorVertex bvertex2ndPass = new BilliorVertex(bfield);
- bvertex2ndPass.doBeamSpotConstraint(false);
- double[] v0 = {bvtxPosUC.e(0, 0), bvtxPosUC.e(1, 0), bvtxPosUC.e(2, 0)};
- bvertex2ndPass.setV0(v0);
- bvertex2ndPass.tryNewFormalism(btlist);
-// bvertex2ndPass.fitVertex(btlist);
- BasicMatrix bvtxPos2ndPass = (BasicMatrix) bvertex2ndPass.getVertexPosition();
- BasicMatrix bvtxCov2ndPass = (BasicMatrix) bvertex2ndPass.getVertexCovariance();
-// System.out.println("2ndPass");
-// System.out.println("Vertex Position: " + bvtxPos2ndPass.toString());
-// System.out.println("chisq : " + bvertex2ndPass.getChiSq());
- aida.histogram1D("BilliorVertex X -- 2ndPass", 100, -10, 20).fill(bvtxPos2ndPass.e(0, 0));
- aida.histogram1D("BilliorVertex Y -- 2ndPass", 100, -0.4, 0.4).fill(bvtxPos2ndPass.e(1, 0));
- aida.histogram1D("BilliorVertex Z -- 2ndPass", 100, -0.4, 0.4).fill(bvtxPos2ndPass.e(2, 0));
- aida.histogram1D("BilliorVertex ChiSq -- 2ndPass", 100, 0, 50).fill(bvertex2ndPass.getChiSq());
- aida.histogram1D("BilliorVertex X Pull -- 2ndPass", 100, -4, 4).fill(bvtxPos2ndPass.e(0, 0) / Math.sqrt(bvtxCov2ndPass.e(0, 0)));
- aida.histogram1D("BilliorVertex Y Pull-- 2ndPass", 100, -4, 4).fill(bvtxPos2ndPass.e(1, 0) / Math.sqrt(bvtxCov2ndPass.e(1, 1)));
- aida.histogram1D("BilliorVertex Z Pull-- 2ndPass", 100, -4, 4).fill(bvtxPos2ndPass.e(2, 0) / Math.sqrt(bvtxCov2ndPass.e(2, 2)));
-
-
- BilliorVertex bvertex = new BilliorVertex(bfield);
- bvertex.doBeamSpotConstraint(true);
- bvertex.tryNewFormalism(btlist);
-// bvertex.fitVertex(btlist);
- BasicMatrix bvtxPos = (BasicMatrix) bvertex.getVertexPosition();
- BasicMatrix bvtxCov = (BasicMatrix) bvertex.getVertexCovariance();
+
+ BilliorVertexer bvertexer = new BilliorVertexer(bfield);
+ bvertexer.setBeamSize(beamsize);
+ bvertexer.doBeamSpotConstraint(true);
+
+ BilliorVertex bvertex=bvertexer.fitVertex(btlist);
+ BasicMatrix bvtxPos = (BasicMatrix) bvertex.getPosition();
+ SymmetricMatrix bvtxCov = bvertex.getCovMatrix();
+ double invMass=bvertex.getParameters().get("invMass");
+
// System.out.println("Constrained");
// System.out.println("Vertex Position: " + bvtxPos.toString());
// System.out.println("chisq : " + bvertex.getChiSq());
aida.histogram1D("BilliorVertex X -- Constrained", 100, -10, 20).fill(bvtxPos.e(0, 0));
aida.histogram1D("BilliorVertex Y -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(1, 0));
aida.histogram1D("BilliorVertex Z -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(2, 0));
- aida.histogram1D("BilliorVertex ChiSq -- Constrained", 100, -10, 50).fill(bvertex.getChiSq());
+ aida.histogram1D("BilliorVertex ChiSq -- Constrained", 100, -10, 50).fill(bvertex.getChi2());
aida.histogram1D("BilliorVertex X Pull -- Constrained", 100, -4, 4).fill(bvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
aida.histogram1D("BilliorVertex Y Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
aida.histogram1D("BilliorVertex Z Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
- aida.histogram1D("BilliorVertex Mass -- Constrained", 250, 0.0, 0.25).fill(bvertex.getInvMass());
- aida.histogram1D("BilliorVertex Mass -- UnConstrained", 250, 0.0, 0.25).fill(bvertexUC.getInvMass());
- aida.histogram1D("BilliorVertex Mass -- 2nd Pass", 250, 0.0, 0.25).fill(bvertex2ndPass.getInvMass());
+ aida.histogram1D("BilliorVertex Mass -- Constrained", 250, 0.0, 0.25).fill(invMass);
+ aida.histogram1D("BilliorVertex Mass -- UnConstrained", 250, 0.0, 0.25).fill(invMassUC);
- BilliorVertex bsconfit = new BilliorVertex(bfield);
- bsconfit.setBeamSize(beamsize);
- bsconfit.doBeamSpotConstraint(false);
- bsconfit.constrainV0toBeamSpot(true);
- bsconfit.tryNewFormalism(btlist);
- BasicMatrix bsconvtxPos = (BasicMatrix) bsconfit.getVertexPosition();
- BasicMatrix bsconvtxCov = (BasicMatrix) bsconfit.getVertexCovariance();
+ BilliorVertexer bconvertexer = new BilliorVertexer(bfield);
+ bconvertexer.setBeamSize(beamsize);
+ bconvertexer.constrainV0toBeamSpot(true);
+
+ BilliorVertex bsconfit=bconvertexer.fitVertex(btlist);
+// bvertexUC.fitVertex(btlist);
+ BasicMatrix bsconvtxPos = (BasicMatrix) bsconfit.getPosition();
+ SymmetricMatrix bsconvtxCov = bsconfit.getCovMatrix();
+ double invMassBSCon=bsconfit.getParameters().get("invMass");
+
aida.histogram1D("BilliorVertex X -- BS Constrained", 100, -10, 20).fill(bsconvtxPos.e(0, 0));
aida.histogram1D("BilliorVertex Y -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(1, 0));
aida.histogram1D("BilliorVertex Z -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(2, 0));
- aida.histogram1D("BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChiSq());
- aida.histogram1D("BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
- aida.histogram1D("BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
- aida.histogram1D("BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
+ aida.histogram1D("BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChi2());
+ aida.histogram1D("BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bsconvtxCov.e(0, 0)));
+ aida.histogram1D("BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bsconvtxCov.e(1, 1)));
+ aida.histogram1D("BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bsconvtxCov.e(2, 2)));
- aida.histogram1D("BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(bsconfit.getInvMass());
+ aida.histogram1D("BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(invMassBSCon);
}
}
}
hps-java/src/main/java/org/lcsim/hps/examples
diff -u -r1.2 -r1.3
--- DetailedAnalysisDriver.java 1 Aug 2011 18:25:21 -0000 1.2
+++ DetailedAnalysisDriver.java 12 Mar 2013 19:40:11 -0000 1.3
@@ -46,6 +46,7 @@
import org.lcsim.hps.recon.vertexing.StraightLineTrack;
import org.lcsim.hps.recon.tracking.FindableTrack;
import org.lcsim.hps.recon.tracking.TrackAnalysis;
+import org.lcsim.hps.recon.vertexing.BilliorVertexer;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
@@ -53,8 +54,8 @@
import org.lcsim.util.aida.AIDA;
/**
- *
- * @author partridge
+
+ @author partridge
*/
public class DetailedAnalysisDriver extends Driver {
@@ -121,19 +122,21 @@
FileWriter fw;
PrintWriter pw;
double[] beamsize = {0.001, 0.02, 0.02};
- Integer _minLayers=8;
- // flipSign is a kludge...
+ Integer _minLayers = 8;
+ // flipSign is a kludge...
// HelicalTrackFitter doesn't deal with B-fields in -ive Z correctly
// so we set the B-field in +iveZ and flip signs of fitted tracks
// note: this should be -1 for Test configurations and +1 for Full (v3.X and lower) configurations
// this is set by the _config variable (detType in HeavyPhotonDriver)
- int flipSign=1;
- String _config="3pt4";
- public DetailedAnalysisDriver(int trackerLayers,int mintrkLayers, String config) {
+ int flipSign = 1;
+ String _config = "3pt4";
+
+ public DetailedAnalysisDriver(int trackerLayers, int mintrkLayers, String config) {
nlayers[0] = trackerLayers;
- _minLayers=mintrkLayers;
- _config=config;
- if(_config.contains("Test")) flipSign=-1;
+ _minLayers = mintrkLayers;
+ _config = config;
+ if (_config.contains("Test"))
+ flipSign = -1;
// Define the efficiency histograms
IHistogramFactory hf = aida.histogramFactory();
@@ -422,9 +425,9 @@
double slopeErr = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.slopeIndex, HelicalTrackFit.slopeIndex));
double curveErr = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.curvatureIndex, HelicalTrackFit.curvatureIndex));
_nchRec++;
- if (track.getCharge()*flipSign < 0)
+ if (track.getCharge() * flipSign < 0)
_neleRec++;
- if (track.getCharge()*flipSign > 0)
+ if (track.getCharge() * flipSign > 0)
_nposRec++;
SeedTrack stEle = (SeedTrack) track;
@@ -488,9 +491,9 @@
// Make plots for fake, non-fake, and all tracks
if (purity < 0.5) {
- if (track.getCharge()*flipSign < 0)
+ if (track.getCharge() * flipSign < 0)
_neleFake++;
- if (track.getCharge()*flipSign > 0)
+ if (track.getCharge() * flipSign > 0)
_nposFake++;
cthfake.fill(cth, 1.0);
phifake.fill(phi, 1.0);
@@ -499,9 +502,9 @@
fillTrackInfo(trackdir, "fake tracks", track.getChi2(), nhits, p, pperp, px, py, pz, phi, cth, d0, xoca, yoca, z0);
} else {
- if (track.getCharge()*flipSign < 0)
+ if (track.getCharge() * flipSign < 0)
_neleTru++;
- if (track.getCharge()*flipSign > 0)
+ if (track.getCharge() * flipSign > 0)
_nposTru++;
cthfake.fill(cth, 0.0);
phifake.fill(phi, 0.0);
@@ -601,7 +604,7 @@
if (htlayer == 1)
l1DeltaZ.put(track, z - zTr);
- if (purity == 1 && track.getCharge()*flipSign > 0 && nhits == 10) {
+ if (purity == 1 && track.getCharge() * flipSign > 0 && nhits == 10) {
if (clusterlist.get(0).rawhits().size() == 1 && clusterlist.get(1).rawhits().size() == 1) {
aida.cloud1D(hitdir + tkresid + "SingleStrip--Track delta y: Layer " + htlayer).fill(y - yTr);
aida.cloud1D(hitdir + tkresid + "SingleStrip--Track delta z: Layer " + htlayer).fill(z - zTr);
@@ -863,7 +866,7 @@
for (Track track1 : tracklist) {
Track ele = null;
Track pos = null;
- int ch1 = track1.getCharge()*flipSign;
+ int ch1 = track1.getCharge() * flipSign;
int index = tracklist.indexOf(track1);
List<Track> subtracklist = tracklist.subList(index, tracklist.size());
for (Track track2 : subtracklist) {
@@ -887,65 +890,51 @@
-
- BilliorVertex bvertexUC = new BilliorVertex(bfield);
- bvertexUC.doBeamSpotConstraint(false);
- bvertexUC.tryNewFormalism(btlist);
-
- BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getVertexPosition();
- BasicMatrix bvtxCovUC = (BasicMatrix) bvertexUC.getVertexCovariance();
- double invMass = bvertexUC.getInvMass();
-
- aida.histogram1D(vertex + "BilliorVertex X -- UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0));
- aida.histogram1D(vertex + "BilliorVertex Y -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0));
- aida.histogram1D(vertex + "BilliorVertex Z -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0));
- aida.histogram1D(vertex + "BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChiSq());
- aida.histogram1D(vertex + "BilliorVertex X Pull -- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) / Math.sqrt(bvtxCovUC.e(0, 0)));
- aida.histogram1D(vertex + "BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1)));
- aida.histogram1D(vertex + "BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2)));
- aida.histogram1D(vertex + "BilliorVertex Mass -- UnConstrained", 250, 0.0, 0.25).fill(bvertexUC.getInvMass());
+ BilliorVertexer bvertexerUC = new BilliorVertexer(bfield);
+ BilliorVertex bvertexUC = bvertexerUC.fitVertex(btlist);
+// bvertexUC.fitVertex(btlist);
+ BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getPosition();
+ SymmetricMatrix bvtxCovUC = bvertexUC.getCovMatrix();
+ double invMassUC = bvertexUC.getParameters().get("invMass");
+// System.out.println("UnConstrained");
+// System.out.println("Vertex Position: " + bvtxPosUC.toString());
+// System.out.println("chisq : " + bvertexUC.getChiSq());
+ aida.histogram1D("BilliorVertex X -- UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0));
+ aida.histogram1D("BilliorVertex Y -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0));
+ aida.histogram1D("BilliorVertex Z -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0));
+ aida.histogram1D("BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChi2());
+ aida.histogram1D("BilliorVertex X Pull -- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) / Math.sqrt(bvtxCovUC.e(0, 0)));
+ aida.histogram1D("BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1)));
+ aida.histogram1D("BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2)));
- aida.cloud1D(apdir + "e+e- Invariant Mass").fill(invMass);
+ aida.cloud1D(apdir + "e+e- Invariant Mass").fill(invMassUC);
if (eleMC != null && posMC != null && ele == eleID && pos == posID)
- aida.cloud1D(apdir + "Matched A' Invariant Mass").fill(invMass);
+ aida.cloud1D(apdir + "Matched A' Invariant Mass").fill(invMassUC);
+
+
+ BilliorVertexer bconvertexer = new BilliorVertexer(bfield);
+ bconvertexer.setBeamSize(beamsize);
+ bconvertexer.constrainV0toBeamSpot(true);
+ BilliorVertex bsconfit = bconvertexer.fitVertex(btlist);
+// bvertexUC.fitVertex(btlist);
+ BasicMatrix bsconvtxPos = (BasicMatrix) bsconfit.getPosition();
+ SymmetricMatrix bsconvtxCov = bsconfit.getCovMatrix();
+ double invMassBSCon = bsconfit.getParameters().get("invMass");
- BilliorVertex bvertex = new BilliorVertex(bfield);
- bvertex.doBeamSpotConstraint(true);
- bvertex.tryNewFormalism(btlist);
- BasicMatrix bvtxPos = (BasicMatrix) bvertex.getVertexPosition();
- BasicMatrix bvtxCov = (BasicMatrix) bvertex.getVertexCovariance();
-
- aida.histogram1D(vertex + "BilliorVertex X -- Constrained", 100, -10, 20).fill(bvtxPos.e(0, 0));
- aida.histogram1D(vertex + "BilliorVertex Y -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(1, 0));
- aida.histogram1D(vertex + "BilliorVertex Z -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(2, 0));
- aida.histogram1D(vertex + "BilliorVertex ChiSq -- Constrained", 100, -10, 50).fill(bvertex.getChiSq());
- aida.histogram1D(vertex + "BilliorVertex X Pull -- Constrained", 100, -4, 4).fill(bvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
- aida.histogram1D(vertex + "BilliorVertex Y Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
- aida.histogram1D(vertex + "BilliorVertex Z Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
-
- aida.histogram1D(vertex + "BilliorVertex Mass -- Constrained", 250, 0.0, 0.25).fill(bvertex.getInvMass());
+ aida.histogram1D("BilliorVertex X -- BS Constrained", 100, -10, 20).fill(bsconvtxPos.e(0, 0));
+ aida.histogram1D("BilliorVertex Y -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(1, 0));
+ aida.histogram1D("BilliorVertex Z -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(2, 0));
+ aida.histogram1D("BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChi2());
+ aida.histogram1D("BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bsconvtxCov.e(0, 0)));
+ aida.histogram1D("BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bsconvtxCov.e(1, 1)));
+ aida.histogram1D("BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bsconvtxCov.e(2, 2)));
+ aida.histogram1D("BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(invMassBSCon);
- BilliorVertex bsconfit = new BilliorVertex(bfield);
- bsconfit.setBeamSize(beamsize);
- bsconfit.doBeamSpotConstraint(false);
- bsconfit.constrainV0toBeamSpot(true);
- bsconfit.tryNewFormalism(btlist);
- BasicMatrix bsconvtxPos = (BasicMatrix) bsconfit.getVertexPosition();
- BasicMatrix bsconvtxCov = (BasicMatrix) bsconfit.getVertexCovariance();
-
- aida.histogram1D(vertex + "BilliorVertex X -- BS Constrained", 100, -10, 20).fill(bsconvtxPos.e(0, 0));
- aida.histogram1D(vertex + "BilliorVertex Y -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(1, 0));
- aida.histogram1D(vertex + "BilliorVertex Z -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(2, 0));
- aida.histogram1D(vertex + "BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChiSq());
- aida.histogram1D(vertex + "BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
- aida.histogram1D(vertex + "BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
- aida.histogram1D(vertex + "BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
- aida.histogram1D(vertex + "BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(bsconfit.getInvMass());
}
}
}
@@ -991,7 +980,7 @@
// Check cases where we have multiple tracks associated with this MC particle
Set<Track> trklist = trktomc.allTo(mcp);
int ntrk = trklist.size();
- int ntrkAxial=0;
+ int ntrkAxial = 0;
if (event.hasItem("AxialTracks")) {
Set<Track> trklistAxial = trktomcAxial.allTo(mcp);
ntrkAxial = trklistAxial.size();
hps-java/src/main/java/org/lcsim/hps/examples
diff -u -r1.2 -r1.3
--- FastTrackAnalysisDriver.java 1 Aug 2011 18:25:21 -0000 1.2
+++ FastTrackAnalysisDriver.java 12 Mar 2013 19:40:11 -0000 1.3
@@ -1,6 +1,7 @@
package org.lcsim.hps.examples;
import hep.physics.matrix.BasicMatrix;
+import hep.physics.matrix.SymmetricMatrix;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
@@ -36,19 +37,20 @@
import org.lcsim.hps.recon.vertexing.StraightLineTrack;
import org.lcsim.hps.recon.tracking.FindableTrack;
import org.lcsim.hps.recon.tracking.TrackAnalysis;
+import org.lcsim.hps.recon.vertexing.BilliorVertexer;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
import org.lcsim.util.Driver;
/**
- *
- * @author mgraham
+
+ @author mgraham
*/
public class FastTrackAnalysisDriver extends Driver {
String[] detNames = {"Tracker"};
- Integer _minLayers=8;
+ Integer _minLayers = 8;
Integer[] nlayers = {8};
int nevt = 0;
double xref = 50.0; //mm
@@ -56,22 +58,22 @@
FileWriter fw;
PrintWriter pw;
double[] beamsize = {0.001, 0.02, 0.02};
-
- String _config="3pt4";
+ String _config = "3pt4";
// flipSign is a kludge...
// HelicalTrackFitter doesn't deal with B-fields in -ive Z correctly
// so we set the B-field in +iveZ and flip signs of fitted tracks
// note: this should be -1 for Test configurations and +1 for Full (v3.X and lower) configurations
// this is set by the _config variable (detType in HeavyPhotonDriver)
- int flipSign=1;
+ int flipSign = 1;
+
public FastTrackAnalysisDriver(int trackerLayers, int mintrkLayers, String config) {
nlayers[0] = trackerLayers;
- _minLayers=mintrkLayers;
- _config=config;
- if(_config.contains("Test")) flipSign=-1;
+ _minLayers = mintrkLayers;
+ _config = config;
+ if (_config.contains("Test"))
+ flipSign = -1;
}
-
public void process(
EventHeader event) {
if (nevt == 0)
@@ -85,7 +87,7 @@
// Increment the event counter
nevt++;
-
+
Hep3Vector IP = new BasicHep3Vector(0., 0., 0.);
double bfield = event.getDetector().getFieldMap().getField(IP).z();
@@ -99,7 +101,8 @@
int nAxialHitsTotal = axialhits.size();
int nL1Hits = 0;
for (HelicalTrackHit hth : axialhits)
- if (hth.Layer() == 1) nL1Hits++;
+ if (hth.Layer() == 1)
+ nL1Hits++;
// Create a relational table that maps TrackerHits to MCParticles
@@ -192,10 +195,14 @@
double zerr = Math.sqrt(cross.getCorrectedCovMatrix().e(2, 2));
int htlayer = htc.Layer();
- if (htlayer == 1) zlist[0] = z;
- if (htlayer == 3) zlist[1] = z;
- if (htlayer == 5) zlist[2] = z;
- if (htlayer == 1) l1DeltaZ.put(track, z - zTr);
+ if (htlayer == 1)
+ zlist[0] = z;
+ if (htlayer == 3)
+ zlist[1] = z;
+ if (htlayer == 5)
+ zlist[2] = z;
+ if (htlayer == 1)
+ l1DeltaZ.put(track, z - zTr);
for (HelicalTrackStrip cl : clusterlist) {
int layer = cl.layer();
HelicalTrackStrip nearest = getNearestHit(cl, toththits);
@@ -260,11 +267,11 @@
for (Track track1 : tracklist) {
Track ele = null;
Track pos = null;
- int ch1 = track1.getCharge()*flipSign;
+ int ch1 = track1.getCharge() * flipSign;
int index = tracklist.indexOf(track1);
List<Track> subtracklist = tracklist.subList(index, tracklist.size());
for (Track track2 : subtracklist) {
- int ch2 = track2.getCharge()*flipSign;
+ int ch2 = track2.getCharge() * flipSign;
if (track1 != track2 && ch1 == -ch2) {
ele = track1;
pos = track2;
@@ -313,99 +320,57 @@
- BilliorVertex vtxfit = new BilliorVertex(bfield);
- vtxfit.doBeamSpotConstraint(false);
- vtxfit.setBeamSize(beamsize);
-// vtxfit.fitVertex(btlist);
- vtxfit.tryNewFormalism(btlist);
- double[] momE = vtxfit.getFittedMomentum(0);
- double[] momP = vtxfit.getFittedMomentum(1);
- double vtxpxE = momE[0];
- double vtxpyE = momE[1];
- double vtxpzE = momE[2];
- double vtxpxP = momP[0];
- double vtxpyP = momP[1];
- double vtxpzP = momP[2];
- double chisq = vtxfit.getChiSq();
- BasicMatrix vtx = (BasicMatrix) vtxfit.getVertexPosition();
- BasicMatrix vtxcov = (BasicMatrix) vtxfit.getVertexCovariance();
+ BilliorVertexer vtxfitter = new BilliorVertexer(bfield);
+ vtxfitter.doBeamSpotConstraint(false);
+ vtxfitter.setBeamSize(beamsize);
+
+ BilliorVertexer vtxfitterCon = new BilliorVertexer(bfield);
+ vtxfitterCon.doBeamSpotConstraint(true);
+ vtxfitterCon.setBeamSize(beamsize);
+
+ BilliorVertexer vtxfitterBSCon = new BilliorVertexer(bfield);
+ vtxfitterBSCon.constrainV0toBeamSpot(true);
+ vtxfitterBSCon.setBeamSize(beamsize);
+
+
+ BilliorVertex vtxfit = vtxfitter.fitVertex(btlist);
+ Map<String, Double> vtxMap = vtxfit.getParameters();
+
+ double vtxpxE = vtxMap.get("p1X");
+ double vtxpyE = vtxMap.get("p1Y");
+ double vtxpzE = vtxMap.get("p1Z");
+ double vtxpxP = vtxMap.get("p2X");
+ double vtxpyP = vtxMap.get("p2Y");
+ double vtxpzP = vtxMap.get("p2Z");
+ double chisq = vtxfit.getChi2();
+ BasicMatrix vtx = (BasicMatrix) vtxfit.getPosition();
+ SymmetricMatrix vtxcov = vtxfit.getCovMatrix();
double chisqE = ele.getChi2();
double chisqP = pos.getChi2();
- BilliorVertex confit = new BilliorVertex(bfield);
- confit.setBeamSize(beamsize);
- confit.doBeamSpotConstraint(true);
-// confit.fitVertex(btlist);
- confit.tryNewFormalism(btlist);
- double[] conmomE = confit.getFittedMomentum(0);
- double[] conmomP = confit.getFittedMomentum(1);
- double conpxE = conmomE[0];
- double conpyE = conmomE[1];
- double conpzE = conmomE[2];
- double conpxP = conmomP[0];
- double conpyP = conmomP[1];
- double conpzP = conmomP[2];
- double conchisq = confit.getChiSq();
- BasicMatrix conVtx = (BasicMatrix) confit.getVertexPosition();
- BasicMatrix conVtxCov = (BasicMatrix) confit.getVertexCovariance();
-
- // 2nd pass vertex parameters
- BilliorVertex vtxfit2 = new BilliorVertex(bfield);
- vtxfit2.setBeamSize(beamsize);
- double[] vtxfitv0 = {vtx.e(0, 0), vtx.e(1, 0), vtx.e(2, 0)};
- vtxfit2.setV0(vtxfitv0);
- vtxfit2.doBeamSpotConstraint(false);
-// vtxfit.fitVertex(btlist);
- vtxfit2.tryNewFormalism(btlist);
- double[] mom2E = vtxfit2.getFittedMomentum(0);
- double[] mom2P = vtxfit2.getFittedMomentum(1);
- double vtxpx2E = mom2E[0];
- double vtxpy2E = mom2E[1];
- double vtxpz2E = mom2E[2];
- double vtxpx2P = mom2P[0];
- double vtxpy2P = mom2P[1];
- double vtxpz2P = mom2P[2];
- double chisq2 = vtxfit2.getChiSq();
- BasicMatrix vtx2 = (BasicMatrix) vtxfit2.getVertexPosition();
- BasicMatrix vtxcov2 = (BasicMatrix) vtxfit2.getVertexCovariance();
-
- BilliorVertex confit2 = new BilliorVertex(bfield);
- confit2.setBeamSize(beamsize);
-
- double[] confitv0 = {conVtx.e(0, 0), conVtx.e(1, 0), conVtx.e(2, 0)};
- confit2.setV0(confitv0);
- confit2.doBeamSpotConstraint(true);
-// confit.fitVertex(btlist);
- confit2.tryNewFormalism(btlist);
- double[] conmom2E = confit2.getFittedMomentum(0);
- double[] conmom2P = confit2.getFittedMomentum(1);
- double conpx2E = conmom2E[0];
- double conpy2E = conmom2E[1];
- double conpz2E = conmom2E[2];
- double conpx2P = conmom2P[0];
- double conpy2P = conmom2P[1];
- double conpz2P = conmom2P[2];
- double conchisq2 = confit2.getChiSq();
- BasicMatrix conVtx2 = (BasicMatrix) confit2.getVertexPosition();
- BasicMatrix conVtxCov2 = (BasicMatrix) confit2.getVertexCovariance();
-
- BilliorVertex bsconfit = new BilliorVertex(bfield);
- bsconfit.setBeamSize(beamsize);
- bsconfit.doBeamSpotConstraint(false);
- bsconfit.constrainV0toBeamSpot(true);
-// bsconfit.fitVertex(btlist);
- bsconfit.tryNewFormalism(btlist);
- double[] bsconmomE = bsconfit.getFittedMomentum(0);
- double[] bsconmomP = bsconfit.getFittedMomentum(1);
- double bsconpxE = bsconmomE[0];
- double bsconpyE = bsconmomE[1];
- double bsconpzE = bsconmomE[2];
- double bsconpxP = bsconmomP[0];
- double bsconpyP = bsconmomP[1];
- double bsconpzP = bsconmomP[2];
- double bsconchisq = bsconfit.getChiSq();
- BasicMatrix bsconVtx = (BasicMatrix) bsconfit.getVertexPosition();
- BasicMatrix bsconVtxCov = (BasicMatrix) bsconfit.getVertexCovariance();
+ BilliorVertex confit = vtxfitterCon.fitVertex(btlist);
+ Map<String, Double> conMap = confit.getParameters();
+ double conpxE = conMap.get("p1X");
+ double conpyE = conMap.get("p1Y");
+ double conpzE = conMap.get("p1Z");
+ double conpxP = conMap.get("p2X");
+ double conpyP = conMap.get("p2Y");
+ double conpzP = conMap.get("p2Z");
+ double conchisq = confit.getChi2();
+ BasicMatrix conVtx = (BasicMatrix) confit.getPosition();
+ SymmetricMatrix conVtxCov = confit.getCovMatrix();
+
+ BilliorVertex bsconfit = vtxfitterBSCon.fitVertex(btlist);
+ Map<String, Double> bsconMap = bsconfit.getParameters();
+ double bsconpxE = bsconMap.get("p1X");
+ double bsconpyE = bsconMap.get("p1Y");
+ double bsconpzE = bsconMap.get("p1Z");
+ double bsconpxP = bsconMap.get("p2X");
+ double bsconpyP = bsconMap.get("p2Y");
+ double bsconpzP = bsconMap.get("p2Z");
+ double bsconchisq = bsconfit.getChi2();
+ BasicMatrix bsconVtx = (BasicMatrix) bsconfit.getPosition();
+ SymmetricMatrix bsconVtxCov = bsconfit.getCovMatrix();
double l1minE = -99;
@@ -445,8 +410,10 @@
int eleFromAp = 0;
int posFromAp = 0;
- if (eleMC != null && ele == eleID) eleFromAp = 1;
- if (posMC != null && pos == posID) posFromAp = 1;
+ if (eleMC != null && ele == eleID)
+ eleFromAp = 1;
+ if (posMC != null && pos == posID)
+ posFromAp = 1;
MCParticle mcEle = tkanalEle.getMCParticle();
MCParticle mcPos = tkanalPos.getMCParticle();
double[] pmcEle = {-99, -99, -99};
@@ -482,8 +449,10 @@
l123KinkE = l123KinkAngle.get(ele);
if (l123KinkAngle.get(pos) != null)
l123KinkP = l123KinkAngle.get(pos);
- if (!tkanalEle.hasLayerOne()) nElectron = -nElectron;
- if (!tkanalPos.hasLayerOne()) nPositron = -nPositron;
+ if (!tkanalEle.hasLayerOne())
+ nElectron = -nElectron;
+ if (!tkanalPos.hasLayerOne())
+ nPositron = -nPositron;
pw.format("%d %5.5f %5.5f %5.5f ", nevt, pxE, pyE, pzE);
pw.format("%5.5f %5.5f %5.5f %5.5f %5.5f ", d0E, z0E, slopeE, phi0E, RE);
pw.format("%5.5f %5.5f %5.5f %5.5f ", chisqE, l1minE, l1dzE, l123KinkE);
@@ -508,25 +477,7 @@
double conErrZ = getErr(conVtxCov.e(2, 2));
pw.format("%5.5f %5.5f %5.5f %5.5f %5.5f %5.5f ", conVtx.e(0, 0), conErrX, conVtx.e(1, 0), conErrY, conVtx.e(2, 0), conErrZ);
- pw.format("%5.5f ", conchisq);
-
- pw.format("%5.5f %5.5f %5.5f ", vtxpx2E, vtxpy2E, vtxpz2E);
- pw.format("%5.5f %5.5f %5.5f ", vtxpx2P, vtxpy2P, vtxpz2P);
- pw.format("%5.5f %5.5f %5.5f %5.5f %5.5f %5.5f ", vtx2.e(0, 0), Math.sqrt(vtxcov2.e(0, 0)), vtx2.e(1, 0), Math.sqrt(vtxcov2.e(1, 1)), vtx2.e(2, 0), Math.sqrt(vtxcov2.e(2, 2)));
- pw.format("%5.5f ", chisq2);
-
- pw.format("%5.5f %5.5f %5.5f ", conpx2E, conpy2E, conpz2E);
- pw.format("%5.5f %5.5f %5.5f ", conpx2P, conpy2P, conpz2P);
-
- // get the errors on the constrained vertex an make sure they aren't NaN
- // there must be something wrong in the constraint...hopefully just the error calc
- double conErrX2 = getErr(conVtxCov2.e(0, 0));
- double conErrY2 = getErr(conVtxCov2.e(1, 1));
- double conErrZ2 = getErr(conVtxCov2.e(2, 2));
-
- pw.format("%5.5f %5.5f %5.5f %5.5f %5.5f %5.5f ", conVtx2.e(0, 0), conErrX2, conVtx2.e(1, 0), conErrY2, conVtx2.e(2, 0), conErrZ2);
- pw.format("%5.5f ", conchisq2);
-
+ pw.format("%5.5f ", conchisq);
pw.format("%5.5f %5.5f %5.5f ", bsconpxE, bsconpyE, bsconpzE);
pw.format("%5.5f %5.5f %5.5f ", bsconpxP, bsconpyP, bsconpzP);
@@ -564,24 +515,26 @@
if (mcp.getPDGID() == 622)
mcApr = mcp;
if (mcp.getParents().size() == 1 && mcp.getParents().get(0).getPDGID() == 622) {
- if (mcp.getPDGID() == -11) mcPos = mcp;
- if (mcp.getPDGID() == 11) mcEle = mcp;
+ if (mcp.getPDGID() == -11)
+ mcPos = mcp;
+ if (mcp.getPDGID() == 11)
+ mcEle = mcp;
}
}
//should probably check that each MC particles are found, but they should be...
- double pxMCE=-99;
- double pyMCE=-99;
- double pzMCE=-99;
- double pxMCP=-99;
- double pyMCP=-99;
- double pzMCP=-99;
+ double pxMCE = -99;
+ double pyMCE = -99;
+ double pzMCE = -99;
+ double pxMCP = -99;
+ double pyMCP = -99;
+ double pzMCP = -99;
int findableE = 0;
int findableP = 0;
- int foundE=0;
- int foundP=0;
+ int foundE = 0;
+ int foundP = 0;
pxMCE = mcEle.getPX();
- pyMCE = mcEle.getPY();
- pzMCE = mcEle.getPZ();
+ pyMCE = mcEle.getPY();
+ pzMCE = mcEle.getPZ();
if (findable.InnerTrackerIsFindable(mcEle, _minLayers))
findableE = 1;
Set<Track> trklistE = trktomc.allTo(mcEle);
@@ -616,7 +569,6 @@
}
-
public void endOfData() {
pw.close();
@@ -652,7 +604,8 @@
double dzdx1 = slt1.dzdx();
double s1sq = 1 + 1 / (dydx1 * dydx1) + (dzdx1 * dzdx1) / (dydx1 * dydx1);
double truep1y = Math.sqrt(p1mag2 / s1sq);
- if (dydx1 < 0) truep1y = -truep1y;
+ if (dydx1 < 0)
+ truep1y = -truep1y;
double truep1x = truep1y / dydx1;
double truep1z = dzdx1 * truep1x;
@@ -666,7 +619,8 @@
double dzdx2 = slt2.dzdx();
double s2sq = 1 + 1 / (dydx2 * dydx2) + (dzdx2 * dzdx2) / (dydx2 * dydx2);
double truep2y = Math.sqrt(p2mag2 / s2sq);
- if (dydx2 < 0) truep2y = -truep2y;
+ if (dydx2 < 0)
+ truep2y = -truep2y;
double truep2x = truep2y / dydx2;
double truep2z = dzdx2 * truep2x;
@@ -701,7 +655,6 @@
return doca;
}
-
//find the XOCA to the beamline extrpolating linearly from the reference point
private double findXoca(double y, double z, double px, double py, double pz) {
double xoca = 0;
@@ -754,7 +707,8 @@
double dzdx1 = slt1.dzdx();
double s1sq = 1 + 1 / (dydx1 * dydx1) + (dzdx1 * dzdx1) / (dydx1 * dydx1);
truep[1] = Math.sqrt(p1mag2 / s1sq);
- if (dydx1 < 0) truep[1] = -truep[1];
+ if (dydx1 < 0)
+ truep[1] = -truep[1];
truep[0] = truep[1] / dydx1;
truep[2] = dzdx1 * truep[0];
return new BasicHep3Vector(truep[0], truep[1], truep[2]);
@@ -846,9 +800,9 @@
}
private double getErr(double errSquared) {
- if (errSquared > 0) return Math.sqrt(errSquared);
- else return -99;
+ if (errSquared > 0)
+ return Math.sqrt(errSquared);
+ else
+ return -99;
}
}
-
-
hps-java/src/main/java/org/lcsim/hps/users/mgraham/jlabrotation
diff -u -r1.2 -r1.3
--- DetailedAnalysisDriver.java 1 Aug 2011 18:25:21 -0000 1.2
+++ DetailedAnalysisDriver.java 12 Mar 2013 19:40:11 -0000 1.3
@@ -53,8 +53,8 @@
import org.lcsim.util.aida.AIDA;
/**
- *
- * @author partridge
+
+ @author partridge
*/
public class DetailedAnalysisDriver extends Driver {
@@ -205,7 +205,7 @@
// Get the magnetic field
Hep3Vector IP = new BasicHep3Vector(0., 0., 0.1);
// double bfield = event.getDetector().getFieldMap().getField(IP).y();
- double bfield =0.5;
+ double bfield = 0.5;
List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, "RawTrackerHitMaker_RawTrackerHits");
List<SiTrackerHitStrip1D> stripHits = event.get(SiTrackerHitStrip1D.class, "StripClusterer_SiTrackerHitStrip1D");
@@ -872,67 +872,100 @@
btlist.add(btEle);
btlist.add(btPos);
+ /*
-
- BilliorVertex bvertexUC = new BilliorVertex(bfield);
- bvertexUC.doBeamSpotConstraint(false);
- bvertexUC.tryNewFormalism(btlist);
-
- BasicMatrix bvtxPosUC = (BasicMatrix) bvertexUC.getVertexPosition();
- BasicMatrix bvtxCovUC = (BasicMatrix) bvertexUC.getVertexCovariance();
- double invMass = bvertexUC.getInvMass();
-
- aida.histogram1D(vertex + "BilliorVertex X -- UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0));
- aida.histogram1D(vertex + "BilliorVertex Y -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0));
- aida.histogram1D(vertex + "BilliorVertex Z -- UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0));
- aida.histogram1D(vertex + "BilliorVertex ChiSq -- UnConstrained", 100, 0, 50).fill(bvertexUC.getChiSq());
- aida.histogram1D(vertex + "BilliorVertex X Pull -- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) / Math.sqrt(bvtxCovUC.e(0, 0)));
- aida.histogram1D(vertex + "BilliorVertex Y Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1)));
- aida.histogram1D(vertex + "BilliorVertex Z Pull-- UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) / Math.sqrt(bvtxCovUC.e(2, 2)));
- aida.histogram1D(vertex + "BilliorVertex Mass -- UnConstrained", 250, 0.0, 0.25).fill(bvertexUC.getInvMass());
-
-
- aida.cloud1D(apdir + "e+e- Invariant Mass").fill(invMass);
- if (eleMC != null && posMC != null && ele == eleID && pos == posID)
- aida.cloud1D(apdir + "Matched A' Invariant Mass").fill(invMass);
-
-
- BilliorVertex bvertex = new BilliorVertex(bfield);
- bvertex.doBeamSpotConstraint(true);
- bvertex.tryNewFormalism(btlist);
- BasicMatrix bvtxPos = (BasicMatrix) bvertex.getVertexPosition();
- BasicMatrix bvtxCov = (BasicMatrix) bvertex.getVertexCovariance();
-
- aida.histogram1D(vertex + "BilliorVertex X -- Constrained", 100, -10, 20).fill(bvtxPos.e(0, 0));
- aida.histogram1D(vertex + "BilliorVertex Y -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(1, 0));
- aida.histogram1D(vertex + "BilliorVertex Z -- Constrained", 100, -0.4, 0.4).fill(bvtxPos.e(2, 0));
- aida.histogram1D(vertex + "BilliorVertex ChiSq -- Constrained", 100, -10, 50).fill(bvertex.getChiSq());
- aida.histogram1D(vertex + "BilliorVertex X Pull -- Constrained", 100, -4, 4).fill(bvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
- aida.histogram1D(vertex + "BilliorVertex Y Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
- aida.histogram1D(vertex + "BilliorVertex Z Pull-- Constrained", 100, -4, 4).fill(bvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
-
- aida.histogram1D(vertex + "BilliorVertex Mass -- Constrained", 250, 0.0, 0.25).fill(bvertex.getInvMass());
-
-
-
- BilliorVertex bsconfit = new BilliorVertex(bfield);
- bsconfit.setBeamSize(beamsize);
- bsconfit.doBeamSpotConstraint(false);
- bsconfit.constrainV0toBeamSpot(true);
- bsconfit.tryNewFormalism(btlist);
- BasicMatrix bsconvtxPos = (BasicMatrix) bsconfit.getVertexPosition();
- BasicMatrix bsconvtxCov = (BasicMatrix) bsconfit.getVertexCovariance();
-
- aida.histogram1D(vertex + "BilliorVertex X -- BS Constrained", 100, -10, 20).fill(bsconvtxPos.e(0, 0));
- aida.histogram1D(vertex + "BilliorVertex Y -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(1, 0));
- aida.histogram1D(vertex + "BilliorVertex Z -- BS Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(2, 0));
- aida.histogram1D(vertex + "BilliorVertex ChiSq -- BS Constrained", 100, -10, 50).fill(bsconfit.getChiSq());
- aida.histogram1D(vertex + "BilliorVertex X Pull -- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) / Math.sqrt(bvtxCov.e(0, 0)));
- aida.histogram1D(vertex + "BilliorVertex Y Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
- aida.histogram1D(vertex + "BilliorVertex Z Pull-- BS Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) / Math.sqrt(bvtxCov.e(2, 2)));
-
- aida.histogram1D(vertex + "BilliorVertex Mass -- BS Constrained", 100, 0.08, 0.12).fill(bsconfit.getInvMass());
+ BilliorVertex bvertexUC = new BilliorVertex(bfield);
+ bvertexUC.doBeamSpotConstraint(false);
+ bvertexUC.tryNewFormalism(btlist);
+
+ BasicMatrix bvtxPosUC = (BasicMatrix)
+ bvertexUC.getVertexPosition(); BasicMatrix bvtxCovUC =
+ (BasicMatrix) bvertexUC.getVertexCovariance(); double
+ invMass = bvertexUC.getInvMass();
+
+ aida.histogram1D(vertex + "BilliorVertex X --
+ UnConstrained", 100, -10, 20).fill(bvtxPosUC.e(0, 0));
+ aida.histogram1D(vertex + "BilliorVertex Y --
+ UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(1, 0));
+ aida.histogram1D(vertex + "BilliorVertex Z --
+ UnConstrained", 100, -0.4, 0.4).fill(bvtxPosUC.e(2, 0));
+ aida.histogram1D(vertex + "BilliorVertex ChiSq --
+ UnConstrained", 100, 0, 50).fill(bvertexUC.getChiSq());
+ aida.histogram1D(vertex + "BilliorVertex X Pull --
+ UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(0, 0) /
+ Math.sqrt(bvtxCovUC.e(0, 0))); aida.histogram1D(vertex +
+ "BilliorVertex Y Pull-- UnConstrained", 100, -4,
+ 4).fill(bvtxPosUC.e(1, 0) / Math.sqrt(bvtxCovUC.e(1, 1)));
+ aida.histogram1D(vertex + "BilliorVertex Z Pull--
+ UnConstrained", 100, -4, 4).fill(bvtxPosUC.e(2, 0) /
+ Math.sqrt(bvtxCovUC.e(2, 2))); aida.histogram1D(vertex +
+ "BilliorVertex Mass -- UnConstrained", 250, 0.0,
+ 0.25).fill(bvertexUC.getInvMass());
+
+
+ aida.cloud1D(apdir + "e+e- Invariant Mass").fill(invMass);
+ if (eleMC != null && posMC != null && ele == eleID && pos
+ == posID) aida.cloud1D(apdir + "Matched A' Invariant
+ Mass").fill(invMass);
+
+
+ BilliorVertex bvertex = new BilliorVertex(bfield);
+ bvertex.doBeamSpotConstraint(true);
+ bvertex.tryNewFormalism(btlist); BasicMatrix bvtxPos =
+ (BasicMatrix) bvertex.getVertexPosition(); BasicMatrix
+ bvtxCov = (BasicMatrix) bvertex.getVertexCovariance();
+
+ aida.histogram1D(vertex + "BilliorVertex X -- Constrained",
+ 100, -10, 20).fill(bvtxPos.e(0, 0));
+ aida.histogram1D(vertex + "BilliorVertex Y -- Constrained",
+ 100, -0.4, 0.4).fill(bvtxPos.e(1, 0));
+ aida.histogram1D(vertex + "BilliorVertex Z -- Constrained",
+ 100, -0.4, 0.4).fill(bvtxPos.e(2, 0));
+ aida.histogram1D(vertex + "BilliorVertex ChiSq --
+ Constrained", 100, -10, 50).fill(bvertex.getChiSq());
+ aida.histogram1D(vertex + "BilliorVertex X Pull --
+ Constrained", 100, -4, 4).fill(bvtxPos.e(0, 0) /
+ Math.sqrt(bvtxCov.e(0, 0))); aida.histogram1D(vertex +
+ "BilliorVertex Y Pull-- Constrained", 100, -4,
+ 4).fill(bvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
+ aida.histogram1D(vertex + "BilliorVertex Z Pull--
+ Constrained", 100, -4, 4).fill(bvtxPos.e(2, 0) /
+ Math.sqrt(bvtxCov.e(2, 2)));
+
+ aida.histogram1D(vertex + "BilliorVertex Mass --
+ Constrained", 250, 0.0, 0.25).fill(bvertex.getInvMass());
+
+
+
+ BilliorVertex bsconfit = new BilliorVertex(bfield);
+ bsconfit.setBeamSize(beamsize);
+ bsconfit.doBeamSpotConstraint(false);
+ bsconfit.constrainV0toBeamSpot(true);
+ bsconfit.tryNewFormalism(btlist); BasicMatrix bsconvtxPos =
+ (BasicMatrix) bsconfit.getVertexPosition(); BasicMatrix
+ bsconvtxCov = (BasicMatrix) bsconfit.getVertexCovariance();
+
+ aida.histogram1D(vertex + "BilliorVertex X -- BS
+ Constrained", 100, -10, 20).fill(bsconvtxPos.e(0, 0));
+ aida.histogram1D(vertex + "BilliorVertex Y -- BS
+ Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(1, 0));
+ aida.histogram1D(vertex + "BilliorVertex Z -- BS
+ Constrained", 100, -0.4, 0.4).fill(bsconvtxPos.e(2, 0));
+ aida.histogram1D(vertex + "BilliorVertex ChiSq -- BS
+ Constrained", 100, -10, 50).fill(bsconfit.getChiSq());
+ aida.histogram1D(vertex + "BilliorVertex X Pull -- BS
+ Constrained", 100, -4, 4).fill(bsconvtxPos.e(0, 0) /
+ Math.sqrt(bvtxCov.e(0, 0))); aida.histogram1D(vertex +
+ "BilliorVertex Y Pull-- BS Constrained", 100, -4,
+ 4).fill(bsconvtxPos.e(1, 0) / Math.sqrt(bvtxCov.e(1, 1)));
+ aida.histogram1D(vertex + "BilliorVertex Z Pull-- BS
+ Constrained", 100, -4, 4).fill(bsconvtxPos.e(2, 0) /
+ Math.sqrt(bvtxCov.e(2, 2)));
+
+ aida.histogram1D(vertex + "BilliorVertex Mass -- BS
+ Constrained", 100, 0.08, 0.12).fill(bsconfit.getInvMass());
+ */
}
}
}
hps-java/src/main/java/org/lcsim/hps/recon/vertexing
diff -N BilliorVertexer.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ BilliorVertexer.java 12 Mar 2013 19:40:11 -0000 1.1
@@ -0,0 +1,971 @@
+package org.lcsim.hps.recon.vertexing;
+
+import hep.physics.matrix.BasicMatrix;
+import hep.physics.matrix.Matrix;
+import hep.physics.matrix.MatrixOp;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import static java.lang.Math.*;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import org.lcsim.constants.Constants;
+
+/**
+ * @version $Id: BilliorVertexer.java,v 1.1 2013/03/12 19:40:11 mgraham Exp $
+ * @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
+ // default is a constant field along the z axis
+
+ private boolean _debug = false;
+ private double _bField;
+ private boolean _beamspotConstraint = true;
+ private boolean _constrainToBS = false;
+ private double[] _beamSize = {0.001, 0.01, 0.01}; //10um in y and z
+ 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[] _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>();
+
+ ;//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 Matrix _constrainedFit;
+ private Matrix _constrainedCov;
+ private double _chiSq;
+
+ // constructor
+ public BilliorVertexer() {
+ }
+
+ public BilliorVertexer(double bField) {
+ _bField = bField;
+ _beamspotConstraint =false;
+ _constrainToBS = false;
+ }
+
+ public BilliorVertexer(double bField,boolean bsConst, boolean constToBS) {
+ _bField = bField;
+ _beamspotConstraint =bsConst;
+ _constrainToBS = constToBS;
+ if(_beamspotConstraint&&_constrainToBS)
+ System.out.println("BilliorVertexer::Warning!!! Setting both _beamspotConstraint and _constrainToBS to true!");
+ }
+
+ public BilliorVertex fitVertex(List<BilliorTrack> tracks) {
+ _ntracks = tracks.size();
+ follow1985Paper(tracks);
+ if (_beamspotConstraint)
+ addV0fromBSConstraint();
+ else if (_constrainToBS)
+ 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);
+ }
+ return new BilliorVertex((Hep3Vector)_vertexPosition,_covVtx,_chiSq,getInvMass(),pFitMap);
+ }
+
+ public BilliorVertex fitFastVertex(List<BilliorTrack> tracks) {
+ _ntracks = tracks.size();
+ fastVertex(tracks);
+
+ 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
+ * 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() {
+ 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("addV0fromBSConstraint::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("addV0fromBSConstraint::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,-(ptot_y)/(ptot_x)*Vx+Vy, -(ptot_z)/(ptot_x)*Vx+Vz)
+ // 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;
+ }
+ //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));
+ }
+ // 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[1]);
+
+ //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));
+
+ 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);
+ 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));
+
+// 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 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)
+ 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);
+
+// 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(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, 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);
+ //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(_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[1]);
+
+ //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);
+ }
+
+ public void setV0(double[] v0) {
+ _v0 = v0;
+ }
+
+ public void setBeamSize(double[] bs) {
+ _beamSize[0] = bs[0];
+ _beamSize[1] = bs[1];
+ _beamSize[2] = bs[2];
+ }
+
+ public void doBeamSpotConstraint(boolean bsconst) {
+ _beamspotConstraint = bsconst;
+ }
+
+ public void constrainV0toBeamSpot(boolean bsconst) {
+ _constrainToBS = bsconst;
+ }
+
+ public double getChiSq() {
+ return _chiSq;
+ }
+ //in theta/phi/rho for now...should fix this to return Hep3Vector
+
+ private double[] getFittedMomentum(int index) {
+ BasicMatrix pi = (BasicMatrix) _pFit.get(index);
+ double[] mom = {0, 0, 0};
+ 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);
+ 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;
+ }
+
+ private double getInvMass() {
+ double esum = 0.;
+ double pxsum = 0.;
+ double pysum = 0.;
+ double pzsum = 0.;
+ double me = 0.000511;
+ // Loop over tracks
+ for (int i = 0; i < _ntracks; i++) {
+ double[] p = getFittedMomentum(i);
+ double p1x = p[0];
+ double p1y = p[1];
+ double p1z = p[2];
+ double p1mag2 = p1x * p1x + p1y * p1y + p1z * p1z;
+ double e1 = Math.sqrt(p1mag2 + me * me);
+ pxsum += p1x;
+ pysum += p1y;
+ pzsum += p1z;
+ esum += e1;
+ }
+ double psum = Math.sqrt(pxsum * pxsum + pysum * pysum + pzsum * pzsum);
+ double evtmass = esum * esum - psum * psum;
+
+ if (evtmass > 0)
+ return Math.sqrt(evtmass);
+ else
+ return -99;
+ }
+
+ 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();
+ }
+
+ private void follow1985Paper(List<BilliorTrack> tracks) {
+
+ BasicMatrix v0 = new BasicMatrix(3, 1);
+ 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> Gs = new ArrayList<Matrix>();
+ List<Matrix> Ds = new ArrayList<Matrix>();
+ List<Matrix> Es = new ArrayList<Matrix>();
+ List<Matrix> As = new ArrayList<Matrix>();
+ List<Matrix> Bs = new ArrayList<Matrix>();
+ List<Matrix> pis = new ArrayList<Matrix>();
+ List<Matrix> cis = new ArrayList<Matrix>();
+
+ BasicMatrix D0 = new BasicMatrix(3, 3);
+ boolean firstTrack = true;
+ 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]);
+ params.add(tmpPar);
+ double theta = par[2];
+ double phiv = par[3];
+ 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
+ double eps = -vv - .5 * uu * uu * par[4];
+ double zp = v0.e(2, 0) - uu * (1 - vv * par[4]) * cotth;
+ // * phi at vertex with these parameters
+ double phiVert = par[3] + uu * par[4];
+ BasicMatrix p0 = new BasicMatrix(5, 1);
+ p0.setElement(0, 0, eps);
+ p0.setElement(1, 0, zp);
+ p0.setElement(2, 0, theta);
+ p0.setElement(3, 0, phiVert);
+ p0.setElement(4, 0, rho);
+
+ 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, theta);
+ q0.setElement(1, 0, phiVert);
+ q0.setElement(2, 0, rho);
+ q0s.add(q0);
+
+ double cosf = cos(phiVert);
+ double sinf = sin(phiVert);
+ BasicMatrix tmpA = new BasicMatrix(5, 3);
+ tmpA.setElement(0, 0, sinf);
+ tmpA.setElement(0, 1, -cosf);
+ tmpA.setElement(1, 0, -cotth * cosf);
+ tmpA.setElement(1, 1, -cotth * sinf);
+ tmpA.setElement(1, 2, 1);
+ tmpA.setElement(3, 0, -par[4] * cosf);
+ tmpA.setElement(3, 1, -par[4] * sinf);
+
+ BasicMatrix tmpB = new BasicMatrix(5, 3);
+ tmpB.setElement(0, 1, uu);
+ tmpB.setElement(0, 2, -uu * uu / 2);
+ tmpB.setElement(1, 0, uu * (1 + cotth * cotth));
+ tmpB.setElement(1, 1, -vv * cotth);
+ tmpB.setElement(1, 2, uu * vv * cotth);
+ tmpB.setElement(3, 1, 1);
+ tmpB.setElement(3, 2, -uu);
+ tmpB.setElement(2, 0, 1); //partial(theta)/dtheta
+ tmpB.setElement(4, 2, 1); //partial (rho)/drho
+ As.add(tmpA);
+ Bs.add(tmpB);
+
+ BasicMatrix ci = (BasicMatrix) MatrixOp.add(p0, MatrixOp.mult(-1, MatrixOp.mult(tmpA, v0)));
+ ci = (BasicMatrix) MatrixOp.add(ci, MatrixOp.mult(-1, MatrixOp.mult(tmpB, q0)));
+ cis.add(ci);
+ pis.add(MatrixOp.add(tmpPar, MatrixOp.mult(-1, ci)));
+
+ BasicMatrix tmpG = (BasicMatrix) MatrixOp.inverse(bt.covariance());
+ Gs.add(tmpG);
+
+ if (firstTrack)
+ D0 = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(tmpA), MatrixOp.mult(tmpG, tmpA));
+ 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));
+ Ds.add(tmpDi);
+ Es.add(tmpEi);
+ firstTrack = false;
+ }
+
+ //ok, now compute the vertex fit.
+ BasicMatrix tmpCovVtx = D0;
+ BasicMatrix bigsum = new BasicMatrix(3, 1);
+ for (int i = 0; i < _ntracks; i++) {
+ BasicMatrix a = (BasicMatrix) As.get(i);
+ BasicMatrix b = (BasicMatrix) Bs.get(i);
+ BasicMatrix d = (BasicMatrix) Ds.get(i);
+ BasicMatrix e = (BasicMatrix) Es.get(i);
+ BasicMatrix g = (BasicMatrix) Gs.get(i);
+ BasicMatrix p = (BasicMatrix) pis.get(i);
+ BasicMatrix sub = (BasicMatrix) MatrixOp.mult(d, MatrixOp.mult(MatrixOp.inverse(e), MatrixOp.transposed(d)));
+ tmpCovVtx = (BasicMatrix) MatrixOp.add(tmpCovVtx, MatrixOp.mult(-1, sub));
+
+ BasicMatrix aTg = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(a), g);
+ 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)
+ bigsum = (BasicMatrix) MatrixOp.mult(MatrixOp.add(aTg, MinusaTgbeIbTg), p);
+ 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)
+ 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>();
+ 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
+ double chisq = 0;
+ for (int j = 0; j < _ntracks; j++) {
+ BasicMatrix a = (BasicMatrix) As.get(j);
+ BasicMatrix b = (BasicMatrix) Bs.get(j);
+ BasicMatrix d = (BasicMatrix) Ds.get(j);
+ BasicMatrix e = (BasicMatrix) Es.get(j);
+ BasicMatrix g = (BasicMatrix) Gs.get(j);
+ BasicMatrix p = (BasicMatrix) pis.get(j);
+ BasicMatrix c = (BasicMatrix) cis.get(j);
+ BasicMatrix first = (BasicMatrix) MatrixOp.mult(-1, MatrixOp.mult(MatrixOp.inverse(e), MatrixOp.transposed(d)));
+ first = (BasicMatrix) MatrixOp.mult(first, xtilde);
+ 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);
+ 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)
+ System.out.println("\n\nfollow1985Paper::Track #" + j);
+ if (_debug)
+ System.out.println("eps(meas) = " + p.e(0, 0) + " eps(fit) =" + ptilde.e(0, 0));
+ if (_debug)
+ System.out.println("zp(meas) = " + p.e(1, 0) + " zp(fit) =" + ptilde.e(1, 0));
+ if (_debug)
+ System.out.println("theta(meas) = " + p.e(2, 0) + " theta(fit) =" + ptilde.e(2, 0));
+ if (_debug)
+ System.out.println("phi(meas) = " + p.e(3, 0) + " phi(fit) =" + ptilde.e(3, 0));
+ 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);
+ 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 tmpCij = (BasicMatrix) MatrixOp.mult(-1, MatrixOp.mult(MatrixOp.inverse(e), MatrixOp.mult(MatrixOp.transposed(d), tmpC0j)));
+ Cij[i][j] = tmpCij;
+ }
+ BasicMatrix tmppfit = new BasicMatrix(3, 1);
+ tmppfit.setElement(0, 0, qtilde.e(0, 0) + c.e(2, 0));
+ tmppfit.setElement(1, 0, qtilde.e(1, 0) + c.e(3, 0));
+ tmppfit.setElement(2, 0, qtilde.e(2, 0) + c.e(4, 0));
+ pfit.add(tmppfit);
+ }
+
+ if (_debug)
+ System.out.println("follow1985Paper::chi^2 = " + chisq);
+
+ _chiSq = chisq;
+ _covVtx = covVtx;
+ _vertexPosition = xtilde;
+ _pFit = pfit;
+ covMomList = Cij;
+ covVtxMomList = C0j;
+
+ }
+}
hps-java/src/main/java/org/lcsim/hps/recon/vertexing
diff -u -r1.3 -r1.4
--- BilliorVertex.java 29 Aug 2012 15:40:46 -0000 1.3
+++ BilliorVertex.java 12 Mar 2013 19:40:11 -0000 1.4
@@ -4,971 +4,113 @@
import hep.physics.matrix.Matrix;
import java.util.List;
import hep.physics.matrix.MatrixOp;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.Hep3Vector;
import java.util.ArrayList;
import org.lcsim.constants.Constants;
import static java.lang.Math.sin;
import static java.lang.Math.cos;
import static java.lang.Math.tan;
+import java.util.HashMap;
+import java.util.Map;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.Vertex;
/**
- * @version $Id: BilliorVertex.java,v 1.3 2012/08/29 15:40:46 mgraham Exp $
- * @version Vertex tracks using least-squares method laid out by billior etal used in the HPS Java package.
+ @version $Id: BilliorVertex.java,v 1.4 2013/03/12 19:40:11 mgraham Exp $
+ @version Vertex tracks using least-squares method laid out by billior etal used
+ in the HPS Java package.
*/
-public class BilliorVertex {
+public class BilliorVertex implements Vertex {
// the value of the magnetic field in the vicinity of the vertex
// default is a constant field along the z axis
private boolean _debug = false;
- private double _bField;
- private boolean _beamspotConstraint = true;
- private boolean _constrainToBS = false;
- private double[] _beamSize = {0.001, 0.01, 0.01}; //10um in y and z
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[] _vertexPosition = {0., 0.0, 0.0};
- private Matrix _vertexPosition = new BasicMatrix(3, 1);
+ private Hep3Vector _vertexPosition;
private Matrix _covVtx = new BasicMatrix(3, 3);
- private List<Matrix> _pFit = new ArrayList<Matrix>();
-
- ;//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 Matrix _constrainedFit;
- private Matrix _constrainedCov;
+// private List<Matrix> _pFit = new ArrayList<Matrix>();
+// private List<Matrix> covVtxMomList = new ArrayList<Matrix>();
private double _chiSq;
+ private double _invMass;
+ private List<BilliorTrack> _tracks;
+ private Map<Integer, Hep3Vector> _fittedMomentum = new HashMap<Integer, Hep3Vector>();
// constructor
public BilliorVertex() {
}
- public BilliorVertex(double bField) {
- _bField = bField;
- }
-
- public void fitVertex(List<BilliorTrack> tracks) {
- _ntracks = tracks.size();
- makeDerivativeMatrices(tracks);
- makeOtherMatrices();
- calculateVertexPosition();
- calculateMomenta();
- calculateCovariance();
- calculateChisq();
- if (_beamspotConstraint)
- addV0fromBSConstraint();
- else if (_constrainToBS)
- constrainV0toBS();
-
- }
-
- public void tryNewFormalism(List<BilliorTrack> tracks) {
- _ntracks = tracks.size();
- follow1985Paper(tracks);
- if (_beamspotConstraint)
- addV0fromBSConstraint();
- else if (_constrainToBS)
- constrainV0toBS();
- }
-
- public void fitFastVertex(List<BilliorTrack> tracks) {
- _ntracks = tracks.size();
- fastVertex(tracks);
- }
-
- 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
- * 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() {
- 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("addV0fromBSConstraint::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("addV0fromBSConstraint::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,-(ptot_y)/(ptot_x)*Vx+Vy, -(ptot_z)/(ptot_x)*Vx+Vz)
- // 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;
- }
- //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));
- }
- // 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[1]);
-
- //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));
-
- 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);
- 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));
-
-// 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 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)
- 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);
-
-// 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(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, 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);
- //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(_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[1]);
-
- //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);
- }
-
- public Matrix getVertexPosition() {
- return _vertexPosition;
+ BilliorVertex(Hep3Vector vtxPos, Matrix covVtx, double chiSq, double invMass, Map<Integer, Hep3Vector> pFitMap) {
+ _chiSq = chiSq;
+ _covVtx = covVtx;
+ _vertexPosition = vtxPos;
+ _invMass = invMass;
+ _fittedMomentum = pFitMap;
}
- public Matrix getVertexCovariance() {
- return _covVtx;
- }
+ BilliorVertex(Hep3Vector vtxPos, Matrix covVtx, double chiSq, double invMass) {
+ _chiSq = chiSq;
+ _covVtx = covVtx;
+ _vertexPosition = vtxPos;
+ _invMass = invMass;
- public void setV0(double[] v0) {
- _v0 = v0;
}
- public void setBeamSize(double[] bs) {
- _beamSize[0] = bs[0];
- _beamSize[1] = bs[1];
- _beamSize[2] = bs[2];
+ public String toString() {
+ StringBuffer sb = new StringBuffer("Vertex at : \nx= " + _vertexPosition.x() + " +/- " + Math.sqrt(_covVtx.e(0, 0)) + "\ny= " + _vertexPosition.y() + " +/- " + Math.sqrt(_covVtx.e(1, 1)) + "\nz= " + _vertexPosition.z() + " +/- " + Math.sqrt(_covVtx.e(2, 2)));
+ return sb.toString();
}
- public void doBeamSpotConstraint(boolean bsconst) {
- _beamspotConstraint = bsconst;
+ @Override
+ public boolean isPrimary() {
+ throw new UnsupportedOperationException("Not supported yet.");
}
- public void constrainV0toBeamSpot(boolean bsconst) {
- _constrainToBS = bsconst;
+ @Override
+ public String getAlgorithmType() {
+ return ("Billior");
}
- public double getChiSq() {
+ @Override
+ public double getChi2() {
return _chiSq;
}
- //in theta/phi/rho for now...should fix this to return Hep3Vector
- public double[] getFittedMomentum(int index) {
- BasicMatrix pi = (BasicMatrix) _pFit.get(index);
- double[] mom = {0, 0, 0};
- 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);
- 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;
+ @Override
+ public double getProbability() {
+ throw new UnsupportedOperationException("Not supported yet.");
}
- public double getInvMass() {
- double esum = 0.;
- double pxsum = 0.;
- double pysum = 0.;
- double pzsum = 0.;
- double me = 0.000511;
- // Loop over tracks
- for (int i = 0; i < _ntracks; i++) {
- double[] p = getFittedMomentum(i);
- double p1x = p[0];
- double p1y = p[1];
- double p1z = p[2];
- double p1mag2 = p1x * p1x + p1y * p1y + p1z * p1z;
- double e1 = Math.sqrt(p1mag2 + me * me);
- pxsum += p1x;
- pysum += p1y;
- pzsum += p1z;
- esum += e1;
- }
- double psum = Math.sqrt(pxsum * pxsum + pysum * pysum + pzsum * pzsum);
- double evtmass = esum * esum - psum * psum;
-
- if (evtmass > 0)
- return Math.sqrt(evtmass);
- else
- return -99;
+ @Override
+ public Hep3Vector getPosition() {
+ return (Hep3Vector) _vertexPosition;
}
- 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();
+ @Override
+ public SymmetricMatrix getCovMatrix() {
+ return (SymmetricMatrix) _covVtx;
}
- private void follow1985Paper(List<BilliorTrack> tracks) {
-
- BasicMatrix v0 = new BasicMatrix(3, 1);
- 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> Gs = new ArrayList<Matrix>();
- List<Matrix> Ds = new ArrayList<Matrix>();
- List<Matrix> Es = new ArrayList<Matrix>();
- List<Matrix> As = new ArrayList<Matrix>();
- List<Matrix> Bs = new ArrayList<Matrix>();
- List<Matrix> pis = new ArrayList<Matrix>();
- List<Matrix> cis = new ArrayList<Matrix>();
-
- BasicMatrix D0 = new BasicMatrix(3, 3);
- boolean firstTrack = true;
- 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]);
- params.add(tmpPar);
- double theta = par[2];
- double phiv = par[3];
- 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
- double eps = -vv - .5 * uu * uu * par[4];
- double zp = v0.e(2, 0) - uu * (1 - vv * par[4]) * cotth;
- // * phi at vertex with these parameters
- double phiVert = par[3] + uu * par[4];
- BasicMatrix p0 = new BasicMatrix(5, 1);
- p0.setElement(0, 0, eps);
- p0.setElement(1, 0, zp);
- p0.setElement(2, 0, theta);
- p0.setElement(3, 0, phiVert);
- p0.setElement(4, 0, rho);
-
- 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, theta);
- q0.setElement(1, 0, phiVert);
- q0.setElement(2, 0, rho);
- q0s.add(q0);
-
- double cosf = cos(phiVert);
- double sinf = sin(phiVert);
- BasicMatrix tmpA = new BasicMatrix(5, 3);
- tmpA.setElement(0, 0, sinf);
- tmpA.setElement(0, 1, -cosf);
- tmpA.setElement(1, 0, -cotth * cosf);
- tmpA.setElement(1, 1, -cotth * sinf);
- tmpA.setElement(1, 2, 1);
- tmpA.setElement(3, 0, -par[4] * cosf);
- tmpA.setElement(3, 1, -par[4] * sinf);
-
- BasicMatrix tmpB = new BasicMatrix(5, 3);
- tmpB.setElement(0, 1, uu);
- tmpB.setElement(0, 2, -uu * uu / 2);
- tmpB.setElement(1, 0, uu * (1 + cotth * cotth));
- tmpB.setElement(1, 1, -vv * cotth);
- tmpB.setElement(1, 2, uu * vv * cotth);
- tmpB.setElement(3, 1, 1);
- tmpB.setElement(3, 2, -uu);
- tmpB.setElement(2, 0, 1); //partial(theta)/dtheta
- tmpB.setElement(4, 2, 1); //partial (rho)/drho
- As.add(tmpA);
- Bs.add(tmpB);
-
- BasicMatrix ci = (BasicMatrix) MatrixOp.add(p0, MatrixOp.mult(-1, MatrixOp.mult(tmpA, v0)));
- ci = (BasicMatrix) MatrixOp.add(ci, MatrixOp.mult(-1, MatrixOp.mult(tmpB, q0)));
- cis.add(ci);
- pis.add(MatrixOp.add(tmpPar, MatrixOp.mult(-1, ci)));
-
- BasicMatrix tmpG = (BasicMatrix) MatrixOp.inverse(bt.covariance());
- Gs.add(tmpG);
-
- if (firstTrack)
- D0 = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(tmpA), MatrixOp.mult(tmpG, tmpA));
- 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));
- Ds.add(tmpDi);
- Es.add(tmpEi);
- firstTrack = false;
- }
-
- //ok, now compute the vertex fit.
- BasicMatrix tmpCovVtx = D0;
- BasicMatrix bigsum = new BasicMatrix(3, 1);
- for (int i = 0; i < _ntracks; i++) {
- BasicMatrix a = (BasicMatrix) As.get(i);
- BasicMatrix b = (BasicMatrix) Bs.get(i);
- BasicMatrix d = (BasicMatrix) Ds.get(i);
- BasicMatrix e = (BasicMatrix) Es.get(i);
- BasicMatrix g = (BasicMatrix) Gs.get(i);
- BasicMatrix p = (BasicMatrix) pis.get(i);
- BasicMatrix sub = (BasicMatrix) MatrixOp.mult(d, MatrixOp.mult(MatrixOp.inverse(e), MatrixOp.transposed(d)));
- tmpCovVtx = (BasicMatrix) MatrixOp.add(tmpCovVtx, MatrixOp.mult(-1, sub));
-
- BasicMatrix aTg = (BasicMatrix) MatrixOp.mult(MatrixOp.transposed(a), g);
- 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)
- bigsum = (BasicMatrix) MatrixOp.mult(MatrixOp.add(aTg, MinusaTgbeIbTg), p);
- 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)
- 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>();
- 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
- double chisq = 0;
- for (int j = 0; j < _ntracks; j++) {
- BasicMatrix a = (BasicMatrix) As.get(j);
- BasicMatrix b = (BasicMatrix) Bs.get(j);
- BasicMatrix d = (BasicMatrix) Ds.get(j);
- BasicMatrix e = (BasicMatrix) Es.get(j);
- BasicMatrix g = (BasicMatrix) Gs.get(j);
- BasicMatrix p = (BasicMatrix) pis.get(j);
- BasicMatrix c = (BasicMatrix) cis.get(j);
- BasicMatrix first = (BasicMatrix) MatrixOp.mult(-1, MatrixOp.mult(MatrixOp.inverse(e), MatrixOp.transposed(d)));
- first = (BasicMatrix) MatrixOp.mult(first, xtilde);
- 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);
- 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)
- System.out.println("\n\nfollow1985Paper::Track #" + j);
- if (_debug)
- System.out.println("eps(meas) = " + p.e(0, 0) + " eps(fit) =" + ptilde.e(0, 0));
- if (_debug)
- System.out.println("zp(meas) = " + p.e(1, 0) + " zp(fit) =" + ptilde.e(1, 0));
- if (_debug)
- System.out.println("theta(meas) = " + p.e(2, 0) + " theta(fit) =" + ptilde.e(2, 0));
- if (_debug)
- System.out.println("phi(meas) = " + p.e(3, 0) + " phi(fit) =" + ptilde.e(3, 0));
- 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);
- 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);
[truncated at 1000 lines; 40 more skipped]
CVSspam 0.2.12