Author: [log in to unmask]
Date: Thu Sep 3 19:18:59 2015
New Revision: 3520
Log:
dedupe tracks with identical hit content before making GBL info
Modified:
java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblFitter.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java
java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java Thu Sep 3 19:18:59 2015
@@ -940,10 +940,15 @@
return meanTime;
}
+ public static List<TrackerHit> getStripHits(Track track, RelationalTable hitToStrips, RelationalTable hitToRotated) {
+ List<TrackerHit> hits = new ArrayList<TrackerHit>();
+ for (TrackerHit hit : track.getTrackerHits())
+ hits.addAll(hitToStrips.allFrom(hitToRotated.from(hit)));
+ return hits;
+ }
+
public static boolean hasSharedStrips(Track track1, Track track2, RelationalTable hitToStrips, RelationalTable hitToRotated) {
- Set<TrackerHit> track1hits = new HashSet<TrackerHit>();
- for (TrackerHit hit : track1.getTrackerHits())
- track1hits.addAll(hitToStrips.allFrom(hitToRotated.from(hit)));
+ Set<TrackerHit> track1hits = new HashSet<TrackerHit>(getStripHits(track1, hitToStrips, hitToRotated));
for (TrackerHit hit : track2.getTrackerHits())
for (TrackerHit hts : (Set<TrackerHit>) hitToStrips.allFrom(hitToRotated.from(hit)))
if (track1hits.contains(hts))
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutput.java Thu Sep 3 19:18:59 2015
@@ -15,14 +15,11 @@
import java.util.Map;
import org.hps.recon.tracking.CoordinateTransformations;
import org.hps.recon.tracking.MaterialSupervisor;
-import org.hps.recon.tracking.MaterialSupervisor.DetectorPlane;
-import org.hps.recon.tracking.MaterialSupervisor.ScatteringDetectorVolume;
import org.hps.recon.tracking.MultipleScattering;
import org.hps.recon.tracking.MultipleScattering.ScatterPoint;
import org.hps.recon.tracking.MultipleScattering.ScatterPoints;
import org.hps.recon.tracking.TrackUtils;
import org.hps.recon.tracking.TrackerHitUtils;
-import org.lcsim.constants.Constants;
import org.lcsim.detector.IDetectorElement;
import org.lcsim.detector.tracker.silicon.HpsSiSensor;
import org.lcsim.event.MCParticle;
@@ -38,7 +35,6 @@
import org.lcsim.geometry.Detector;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
import org.lcsim.recon.tracking.digitization.sisim.TrackerHitType;
-import org.lcsim.recon.tracking.seedtracker.ScatterAngle;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
@@ -197,13 +193,6 @@
}
}
- //GBLDATA
- for (int row = 0; row < perToClPrj.getNRows(); ++row) {
- for (int col = 0; col < perToClPrj.getNColumns(); ++col) {
- gtd.setPrjPerToCl(row, col, perToClPrj.e(row, col));
- }
- }
-
// print chi2 of fit
if (textFile != null) {
textFile.printChi2(htf.chisq(), htf.ndf());
@@ -226,8 +215,7 @@
}
// dummy cov matrix for CL parameters
- BasicMatrix clCov = new BasicMatrix(5, 5);
- initUnit(clCov);
+ BasicMatrix clCov = GblUtils.unitMatrix(5, 5);
clCov = (BasicMatrix) MatrixOp.mult(0.1 * 0.1, clCov);
if (textFile != null) {
@@ -657,26 +645,6 @@
// return j;
//
// }
- private void initUnit(BasicMatrix mat) {
- for (int row = 0; row != mat.getNRows(); row++) {
- for (int col = 0; col != mat.getNColumns(); col++) {
- if (row != col) {
- mat.setElement(row, col, 0);
- } else {
- mat.setElement(row, col, 1);
- }
- }
- }
- }
-
- private void initZero(BasicMatrix mat) {
- for (int row = 0; row != mat.getNRows(); row++) {
- for (int col = 0; col != mat.getNColumns(); col++) {
- mat.setElement(row, col, 0);
- }
- }
- }
-
/**
* Transform MCParticle into a Helix object. Note that it produces the helix
* parameters at nominal x=0 and assumes that there is no field at x<0
@@ -723,13 +691,13 @@
p.setElement(0, 0, perPar.getD0());
p.setElement(0, 1, perPar.getPhi());
p.setElement(0, 2, perPar.getKappa());
- p.setElement(0, 0, perPar.getZ0());
+ p.setElement(0, 3, perPar.getZ0());
p.setElement(0, 4, Math.tan(Math.PI / 2.0 - perPar.getTheta()));
BasicMatrix pt = new BasicMatrix(1, 5);
pt.setElement(0, 0, perParTruth.getD0());
pt.setElement(0, 1, perParTruth.getPhi());
pt.setElement(0, 2, perParTruth.getKappa());
- pt.setElement(0, 0, perParTruth.getZ0());
+ pt.setElement(0, 3, perParTruth.getZ0());
pt.setElement(0, 4, Math.tan(Math.PI / 2.0 - perParTruth.getTheta()));
Matrix error_matrix = MatrixOp.inverse(covariance);
BasicMatrix res = (BasicMatrix) MatrixOp.sub(p, pt);
@@ -887,6 +855,28 @@
trans.setElement(2, 1, VecOp.dot(J, T));
trans.setElement(2, 2, VecOp.dot(K, T));
return trans;
+
+ /*
+ Hep3Vector B = new BasicHep3Vector(0, 0, 1); // TODO sign convention?
+ Hep3Vector H = VecOp.mult(1 / bfield, B);
+ Hep3Vector T = HelixUtils.Direction(helix, 0.);
+ Hep3Vector HcrossT = VecOp.cross(H, T);
+ double alpha = HcrossT.magnitude(); // this should be Bvec cross TrackDir/|B|
+ double Q = Math.abs(bfield) * q / p;
+ Hep3Vector Z = new BasicHep3Vector(0, 0, 1);
+ Hep3Vector J = VecOp.mult(1. / VecOp.cross(T, Z).magnitude(), VecOp.cross(T, Z));
+ Hep3Vector K = Z;
+ Hep3Vector U = VecOp.mult(-1, J);
+ Hep3Vector V = VecOp.cross(T, U);
+ Hep3Vector I = VecOp.cross(J, K);
+ Hep3Vector N = VecOp.mult(1 / alpha, VecOp.cross(H, T)); //-cross(T,H)/alpha = -cross(T,Z) = -J
+ double UdotI = VecOp.dot(U, I); // 0,0
+ double NdotV = VecOp.dot(N, V); // 1,1?
+ double NdotU = VecOp.dot(N, U); // 0,1?
+ double TdotI = VecOp.dot(T, I); // 2,0
+ double VdotI = VecOp.dot(V, I); // 1,0
+ double VdotK = VecOp.dot(V, K); // 1,2
+*/
}
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GBLOutputDriver.java Thu Sep 3 19:18:59 2015
@@ -2,20 +2,24 @@
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
-
import java.io.IOException;
import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.HashSet;
import java.util.List;
+import java.util.Map;
+import java.util.Set;
import java.util.logging.Level;
import java.util.logging.Logger;
-
import org.hps.recon.tracking.EventQuality;
import org.hps.recon.tracking.TrackUtils;
import org.lcsim.event.EventHeader;
import org.lcsim.event.LCRelation;
import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
import org.lcsim.event.SimTrackerHit;
import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
import org.lcsim.event.base.MyLCRelation;
import org.lcsim.geometry.Detector;
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
@@ -23,8 +27,7 @@
import org.lcsim.util.aida.AIDA;
/**
- * This driver class is used to
- * 1) write lcio collection of GBL info objects OR
+ * This driver class is used to 1) write lcio collection of GBL info objects OR
* 2) write GBL info into a unstructures text-based output
*
* It uses a helper class that does the actual work. We will port GBL to java
@@ -37,13 +40,13 @@
*/
public class GBLOutputDriver extends Driver {
- private AIDA aida = AIDA.defaultInstance();
+ private final AIDA aida = AIDA.defaultInstance();
int nevt = 0;
GBLOutput gbl = null;
TruthResiduals truthRes = null;
private String gblFileName = "";
private String outputPlotFileName = "";
- private String MCParticleCollectionName = "MCParticle";
+ private final String MCParticleCollectionName = "MCParticle";
private int _debug = 0;
private boolean isMC = false;
private int totalTracks = 0;
@@ -72,46 +75,63 @@
}
@Override
- public void process(EventHeader event) {
+ public void process(EventHeader event) {
// Check if the event contains a collection of the type Track. If it
// doesn't skip the event.
- if (!event.hasCollection(Track.class))
+ if (!event.hasCollection(Track.class)) {
return;
+ }
// Get all collections of the type Track from the event. This is
// required since the event contains a track collection for each of the
// different tracking strategies.
List<List<Track>> trackCollections = event.get(Track.class);
- if (_debug > 0)
+ if (_debug > 0) {
System.out.printf("%s: Event %d has %d tracks\n", this.getClass().getSimpleName(), event.getEventNumber(), trackCollections.size());
+ }
List<SiTrackerHitStrip1D> stripHits = event.get(SiTrackerHitStrip1D.class, "StripClusterer_SiTrackerHitStrip1D");
- if (_debug > 0)
+ if (_debug > 0) {
System.out.printf("%s: Got %d SiTrackerHitStrip1D in this event\n", stripHits.size());
+ }
List<MCParticle> mcParticles = new ArrayList<MCParticle>();
- if (event.hasCollection(MCParticle.class, this.MCParticleCollectionName))
+ if (event.hasCollection(MCParticle.class, this.MCParticleCollectionName)) {
mcParticles = event.get(MCParticle.class, this.MCParticleCollectionName);
+ }
List<SimTrackerHit> simTrackerHits = new ArrayList<SimTrackerHit>();
- if (event.hasCollection(SimTrackerHit.class, "TrackerHits"))
+ if (event.hasCollection(SimTrackerHit.class, "TrackerHits")) {
simTrackerHits = event.get(SimTrackerHit.class, "TrackerHits");
-
- if (isMC)
- if (truthRes != null)
+ }
+
+ if (isMC) {
+ if (truthRes != null) {
truthRes.processSim(mcParticles, simTrackerHits);
+ }
+ }
+
+ RelationalTable hitToStrips = TrackUtils.getHitToStripsTable(event);
+ RelationalTable hitToRotated = TrackUtils.getHitToRotatedTable(event);
// Loop over each of the track collections retrieved from the event
- List<Track> selected_tracks = new ArrayList<Track>();
- for (List<Track> tracklist : trackCollections) {
+ Map<Set<TrackerHit>, List<Track>> hitsToTracksMap = new HashMap<Set<TrackerHit>, List<Track>>();
+ for (List<Track> tracklist : trackCollections) {
for (Track trk : tracklist) {
totalTracks++;
if (TrackUtils.isGoodTrack(trk, tracklist, EventQuality.Quality.NONE)) {
- if (_debug > 0)
- System.out.printf("%s: Track failed selection\n", this.getClass().getSimpleName());
- selected_tracks.add(trk);
+ Set<TrackerHit> trackHits = new HashSet<TrackerHit>(TrackUtils.getStripHits(trk, hitToStrips, hitToRotated));
+ List<Track> matchingTracks = hitsToTracksMap.get(trackHits);
+ if (matchingTracks == null) {
+ matchingTracks = new ArrayList<Track>();
+ hitsToTracksMap.put(trackHits, matchingTracks);
+ }
+ matchingTracks.add(trk);
+ } else if (_debug > 0) {
+ System.out.printf("%s: Track failed selection\n", this.getClass().getSimpleName());
}
}
}
+
// GBLData
// containers and data
List<GBLEventData> gblEventData = new ArrayList<GBLEventData>();
@@ -125,9 +145,13 @@
gbl.printNewEvent(iEvent, gbl.get_B().z());
iTrack = 0;
- for (Track trk : selected_tracks) {
- if (_debug > 0)
+
+ for (List<Track> matchingTracks : hitsToTracksMap.values()) {
+ Track trk = matchingTracks.get(0);//arbitrarily pick one track to use to generate GBL data
+
+ if (_debug > 0) {
System.out.printf("%s: Print GBL output for this track\n", this.getClass().getSimpleName());
+ }
//GBLDATA
GBLTrackData gblTrackData = new GBLTrackData(iTrack);
@@ -139,7 +163,9 @@
//GBLDATA
//add relation to normal track object
- trackToGBLTrackRelationListAll.add(new MyLCRelation(trk, gblTrackData));
+ for (Track matchingTrack : matchingTracks) {
+ trackToGBLTrackRelationListAll.add(new MyLCRelation(matchingTrack, gblTrackData));
+ }
// add strip clusters to lists
for (GBLStripClusterData gblStripClusterData : gblStripDataList) {
// add all strip clusters from this track to output list
@@ -153,12 +179,12 @@
totalTracksProcessed++;
++iTrack;
}
-
- event.put("GBLEventData" , gblEventData, GBLEventData.class, 0);
- event.put("GBLTrackData" , gblTrackDataList, GBLTrackData.class, 0);
- event.put("GBLStripClusterData" , gblStripDataListAll, GBLStripClusterData.class, 0);
- event.put("GBLTrackToStripData" , gblTrackToStripClusterRelationListAll, LCRelation.class, 0);
- event.put("TrackToGBLTrack" , trackToGBLTrackRelationListAll, LCRelation.class, 0);
+
+ event.put("GBLEventData", gblEventData, GBLEventData.class, 0);
+ event.put("GBLTrackData", gblTrackDataList, GBLTrackData.class, 0);
+ event.put("GBLStripClusterData", gblStripDataListAll, GBLStripClusterData.class, 0);
+ event.put("GBLTrackToStripData", gblTrackToStripClusterRelationListAll, LCRelation.class, 0);
+ event.put("TrackToGBLTrack", trackToGBLTrackRelationListAll, LCRelation.class, 0);
++iEvent;
@@ -166,18 +192,20 @@
@Override
public void endOfData() {
- if (gbl != null)
+ if (gbl != null) {
gbl.close();
-
- if (!"".equals(outputPlotFileName))
+ }
+
+ if (!"".equals(outputPlotFileName)) {
try {
aida.saveAs(outputPlotFileName);
} catch (IOException ex) {
Logger.getLogger(GBLOutputDriver.class.getName()).log(Level.SEVERE, "Couldn't save aida plots to file " + outputPlotFileName, ex);
}
+ }
System.out.println(this.getClass().getSimpleName() + ": Total Number of Events = " + iEvent);
System.out.println(this.getClass().getSimpleName() + ": Total Number of Tracks = " + totalTracks);
- System.out.println(this.getClass().getSimpleName() + ": Total Number of Tracks Processed = " + totalTracksProcessed);
+ System.out.println(this.getClass().getSimpleName() + ": Total Number of Unique Tracks Processed = " + totalTracksProcessed);
}
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/GblUtils.java Thu Sep 3 19:18:59 2015
@@ -17,20 +17,10 @@
*/
public class GblUtils {
- private static GblUtils INSTANCE = null;
-
private GblUtils() {
}
- public static GblUtils getInstance() {
- if (INSTANCE == null) {
- return new GblUtils();
- } else {
- return INSTANCE;
- }
- }
-
- public BasicMatrix gblSimpleJacobianLambdaPhi(double ds, double cosl, double bfac) {
+ public static BasicMatrix gblSimpleJacobianLambdaPhi(double ds, double cosl, double bfac) {
/*
Simple jacobian: quadratic in arc length difference.
using lambda phi as directions
@@ -68,7 +58,7 @@
return mat;
}
- public BasicMatrix unitMatrix(int rows, int cols) {
+ public static BasicMatrix unitMatrix(int rows, int cols) {
BasicMatrix mat = new BasicMatrix(rows, cols);
for (int row = 0; row != mat.getNRows(); row++) {
for (int col = 0; col != mat.getNColumns(); col++) {
@@ -82,7 +72,7 @@
return mat;
}
- public BasicMatrix zeroMatrix(int rows, int cols) {
+ public static BasicMatrix zeroMatrix(int rows, int cols) {
BasicMatrix mat = new BasicMatrix(rows, cols);
for (int row = 0; row != mat.getNRows(); row++) {
for (int col = 0; col != mat.getNColumns(); col++) {
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblFitter.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblFitter.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblFitter.java Thu Sep 3 19:18:59 2015
@@ -10,20 +10,15 @@
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
-import java.util.HashMap;
import java.util.List;
-import java.util.Map;
import org.hps.recon.tracking.MaterialSupervisor;
-import org.hps.recon.tracking.MaterialSupervisor.DetectorPlane;
-import org.hps.recon.tracking.MaterialSupervisor.ScatteringDetectorVolume;
import org.hps.recon.tracking.MultipleScattering;
import org.hps.recon.tracking.MultipleScattering.ScatterPoint;
import org.hps.recon.tracking.MultipleScattering.ScatterPoints;
import org.hps.recon.tracking.TrackUtils;
import org.hps.recon.tracking.TrackerHitUtils;
import org.lcsim.constants.Constants;
-import org.lcsim.detector.IDetectorElement;
import org.lcsim.event.RawTrackerHit;
import org.lcsim.event.Track;
import org.lcsim.fit.helicaltrack.HelicalTrackCross;
@@ -31,7 +26,6 @@
import org.lcsim.fit.helicaltrack.HelicalTrackHit;
import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
import org.lcsim.fit.helicaltrack.HelixUtils;
-import org.lcsim.recon.tracking.seedtracker.ScatterAngle;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
@@ -45,7 +39,7 @@
*/
public class HpsGblFitter {
- private boolean _debug = true;
+ private final boolean _debug = true;
private double _B = 0.;
private double _bfac = 0.;
private boolean isMC = false;
@@ -104,7 +98,7 @@
// path length along trajectory
double s = 0.;
// jacobian to transport errors between points along the path
- BasicMatrix jacPointToPoint = GblUtils.getInstance().unitMatrix(5, 5);
+ BasicMatrix jacPointToPoint = GblUtils.unitMatrix(5, 5);
// Option to use uncorrelated MS errors
// This is similar to what is done in lcsim seedtracker
// The msCov below holds the MS errors
@@ -139,6 +133,7 @@
// TODO use actual path length and not layer id!
//Collections.sort(stripClusters, new HelicalTrackStripComparer());
Collections.sort(stripClusters, new Comparator<HelicalTrackStrip>() {
+ @Override
public int compare(HelicalTrackStrip o1, HelicalTrackStrip o2) {
return o1.layer() < o2.layer() ? -1 : o1.layer() > o2.layer() ? 1 : 0;
}
@@ -268,7 +263,7 @@
}
//Find the Jacobian to be able to propagate the covariance matrix to this strip position
- jacPointToPoint = GblUtils.getInstance().gblSimpleJacobianLambdaPhi(step, cosLambda, Math.abs(_bfac));
+ jacPointToPoint = GblUtils.gblSimpleJacobianLambdaPhi(step, cosLambda, Math.abs(_bfac));
if (_debug) {
System.out.printf("%s: jacPointToPoint \n%s\n", this.getClass().getSimpleName(), jacPointToPoint.toString());
@@ -296,7 +291,7 @@
//Add scatterer in curvilinear frame to the point
// no direction in this frame as it moves along the track
- BasicMatrix scat = GblUtils.getInstance().zeroMatrix(0, 2);
+ BasicMatrix scat = GblUtils.zeroMatrix(0, 2);
// find scattering angle
ScatterPoint scatter = scatters.getScatterPoint(((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement());
@@ -421,6 +416,7 @@
public static class HelicalTrackStripComparer implements Comparator<HelicalTrackStrip> {
+ @Override
public int compare(HelicalTrackStrip o1, HelicalTrackStrip o2) {
// TODO Change this to path length!?
return compare(o1.layer(), o2.layer());
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/HpsGblRefitter.java Thu Sep 3 19:18:59 2015
@@ -1,10 +1,8 @@
package org.hps.recon.tracking.gbl;
-import hep.physics.matrix.SymmetricMatrix;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
-import static java.lang.Math.abs;
import static java.lang.Math.abs;
import static java.lang.Math.sin;
import static java.lang.Math.sqrt;
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/gbl/MakeGblTracks.java Thu Sep 3 19:18:59 2015
@@ -216,7 +216,7 @@
// Strandlie, Wittek, NIMA 566, 2006
Matrix covariance_gbl = new Matrix(5, 5);
//helpers
- double Bz = Constants.fieldConversion * Math.abs(bfield); // TODO sign convention and should it be it scaled from Telsa?
+ double Bz = -Constants.fieldConversion * Math.abs(bfield); // TODO sign convention and should it be it scaled from Telsa?
double p = Math.abs(1 / qOverP_gbl);
double q = Math.signum(qOverP_gbl);
double tanLambda = Math.tan(lambda_gbl);
@@ -240,35 +240,46 @@
double TdotI = VecOp.dot(T, I);
double VdotI = VecOp.dot(V, I);
double VdotK = VecOp.dot(V, K);
+ covariance_gbl.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), VdotK / TdotI);
+ covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1);
+ covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -alpha * Q * UdotI * NdotU / (cosLambda * TdotI));
+ covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -alpha * Q * VdotI * NdotU / (cosLambda * TdotI));
covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), -1 * Bz / cosLambda);
+// covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 0);
covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1 * q * Bz * tanLambda / (p * cosLambda));
- covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 0);
covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), q * Bz * alpha * Q * tanLambda * UdotI * NdotV / (p * cosLambda * TdotI));
covariance_gbl.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), q * Bz * alpha * Q * tanLambda * VdotI * NdotV / (p * cosLambda * TdotI));
+ covariance_gbl.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -1 / TdotI);
covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), -1);
covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), alpha * Q * UdotI * NdotV / TdotI);
covariance_gbl.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), alpha * Q * VdotI * NdotV / TdotI);
- covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1);
- covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -alpha * Q * UdotI * NdotU / (cosLambda * TdotI));
- covariance_gbl.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -alpha * Q * VdotI * NdotU / (cosLambda * TdotI));
- covariance_gbl.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), VdotK / TdotI);
- covariance_gbl.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -1 / TdotI);
-
+
+// System.out.println(clToPerPrj);
+
+// covariance_gbl.print(15, 13);
+//// System.out.println(-alpha * Q * UdotI * NdotU / (cosLambda * TdotI));
+// System.out.println(-alpha * Q * VdotI * NdotU / (cosLambda * TdotI));
+// System.out.format("%f %f %f %f %f %f\n", -alpha, Q, VdotI, NdotU, cosLambda, TdotI);
+// System.out.format("%f %f %f %f %f %f\n", -Math.cos(lambda)/Math.abs(bfield), Q, perToClPrj.e(1, 0), perToClPrj.e(0, 1), Math.cos(lambda_gbl), perToClPrj.e(2, 0));
+//
+//// System.out.println(q * Bz * alpha * Q * tanLambda * UdotI * NdotV / (p * cosLambda * TdotI));
+// System.out.println(q * Bz * alpha * Q * tanLambda * VdotI * NdotV / (p * cosLambda * TdotI));
+//// System.out.println(alpha * Q * UdotI * NdotV / TdotI);
+// System.out.println(alpha * Q * VdotI * NdotV / TdotI);
// Sho's magic
Matrix jacobian = new Matrix(5, 5);
jacobian.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.XT.getValue(), -clToPerPrj.e(1, 0));
jacobian.set(HelicalTrackFit.dcaIndex, FittedGblTrajectory.GBLPARIDX.YT.getValue(), -clToPerPrj.e(1, 1));
+ jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1.0);
+ jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), Constants.fieldConversion * Math.abs(bfield) / Math.cos(lambda_gbl));
+ jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl * Math.tan(lambda_gbl) / Math.cos(lambda_gbl));
jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.XT.getValue(), clToPerPrj.e(2, 0));
jacobian.set(HelicalTrackFit.z0Index, FittedGblTrajectory.GBLPARIDX.YT.getValue(), clToPerPrj.e(2, 1));
- jacobian.set(HelicalTrackFit.phi0Index, FittedGblTrajectory.GBLPARIDX.XTPRIME.getValue(), 1.0);
jacobian.set(HelicalTrackFit.slopeIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Math.pow(Math.cos(lambda_gbl), -2.0));
- jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.QOVERP.getValue(), Constants.fieldConversion * Math.abs(bfield) / Math.cos(lambda_gbl));
- jacobian.set(HelicalTrackFit.curvatureIndex, FittedGblTrajectory.GBLPARIDX.YTPRIME.getValue(), Constants.fieldConversion * Math.abs(bfield) * qOverP_gbl * Math.tan(lambda_gbl) / Math.cos(lambda_gbl));
-
-// covariance_gbl.print(10, 8);
-//
-// jacobian.print(10, 8);
-
+
+// jacobian.print(15, 13);
+// System.out.println(-clToPerPrj.e(1, 1));
+// System.out.println(clToPerPrj.e(2, 0));
Matrix helixCovariance = jacobian.times(locCov.times(jacobian.transpose()));
SymmetricMatrix cov = new SymmetricMatrix(5);
for (int i = 0; i < 5; i++) {
