1 added + 5 modified, total 6 files
lcsim/src/org/lcsim/contrib/JanStrube/standalone
diff -N NickTrack_FastMCTrackComparison.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ NickTrack_FastMCTrackComparison.java 7 Dec 2006 06:11:24 -0000 1.1
@@ -0,0 +1,351 @@
+import static java.lang.Math.abs;
+import static java.lang.Math.sqrt;
+
+import org.lcsim.event.*;
+import org.lcsim.event.util.*;
+import org.lcsim.event.base.*;
+import org.lcsim.geometry.Detector;
+import org.lcsim.contrib.JanStrube.tracking.HelixSwimmer;
+import org.lcsim.contrib.JanStrube.tracking.NewFastMCTrackFactory;
+import org.lcsim.contrib.JanStrube.tracking.NewTrack;
+import org.lcsim.contrib.JanStrube.tracking.Track;
+import org.lcsim.contrib.NickSinev.tracking.wmfitter.*;
+import org.lcsim.contrib.NickSinev.tracking.util.*;
+import org.lcsim.spacegeom.SpacePoint;
+import org.lcsim.spacegeom.SpaceVector;
+import org.lcsim.util.Driver;
+import org.lcsim.contrib.NickSinev.ztracking.*;
+import org.lcsim.contrib.NickSinev.event.util.*;
+import org.lcsim.constants.Constants;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.matrix.*;
+import org.lcsim.fit.circle.*;
+import java.text.*;
+import java.util.*;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.swim.Helix;
+import org.lcsim.util.swim.HelixSwimmer.SpatialParameters;
+
+import Jama.*;
+
+public class NickTrack_FastMCTrackComparison extends Driver {
+ int firstev = 0;
+ int lastev = 9999;
+ boolean real = true;
+ double minPt = 0.1;
+ double maxPt = 500.;
+ double minTl = 0.;
+ double maxTl = 5.2;
+ double[] minTllms = { 0.01, 0.4, 0.9, 1.4, 2.2, 3.2 };
+ double[] maxTllms = { 0.4, 0.9, 1.4, 2.2, 3.2, 5. };
+ private double vtxbrphires = 0.0035;
+ private double vtxbzres = 0.035;
+ private double vtxecrpres = 0.0035;
+ private double vtxecrres = 0.0035;
+ private double trkbrpres = 0.007;
+ private double trkbzres = 30.;
+ private double trkecrpres = 0.01;
+ private double trkecrres = 0.01;
+ private SLDFitTrack fitter = new SLDFitTrack();
+ private DecimalFormat df = new DecimalFormat();
+ private AIDA aida = AIDA.defaultInstance();
+ private SmearTrackerHits smTh = new SmearTrackerHits();
+ private TrackingCheaterSmHits cheat = new TrackingCheaterSmHits();
+ private int evnum = 0;
+ double field = 5. * Constants.fieldConversion;
+ int ttn = 0;
+ int ngdft = 0;
+ double trtrlpa = 0.;
+
+ public NickTrack_FastMCTrackComparison() {
+ if (real) {
+ smTh.setVtxBarrRPResolution(vtxbrphires);
+ smTh.setVtxBarrZResolution(vtxbzres);
+ smTh.setVtxECRPResolution(vtxecrpres);
+ smTh.setVtxECRResolution(vtxecrres);
+ smTh.setTrkBarrRPResolution(trkbrpres);
+ smTh.setTrkBarrZResolution(trkbzres);
+ smTh.setTrkECRPResolution(trkecrpres);
+ smTh.setTrkECRResolution(trkecrres);
+ }
+ cheat.useChargedOnly = true;
+ smTh.setDoSmear(true); // ++++++++++++++++++++++++++++++++++turn on
+ // smearing
+ add(smTh);
+ add(cheat);
+ df.setMaximumFractionDigits(2);
+ fitter.setDoHist(false);
+ }
+
+ public void process(EventHeader event) {
+ NewFastMCTrackFactory factory = new NewFastMCTrackFactory(event, false);
+ if (!((evnum > firstev - 1) && (evnum < lastev + 1)))
+ return;
+ if (evnum % 100 == 0) {
+ System.out.println("Processing event " + evnum
+ + " accepted chit tracks " + ttn + " fitted " + ngdft);
+ }
+ super.process(event);
+ if (evnum == firstev) {
+ Detector det = event.getDetector();
+ fitter.setDetector(det);
+ }
+ List<CheatedBaseTrack> tracks = null;
+ try {
+ tracks = event.get(CheatedBaseTrack.class, "CheatedBaseTracks");
+ } catch (Exception trnotf) {
+ System.out.println("No CheatedBaseTracks in the event ! \n");
+ } finally {
+ }
+ if (tracks == null)
+ return;
+
+ int nprts = tracks.size();
+// System.out.println("Tracking cheater has found " + nprts
+// + " BaseTracks in event " + evnum);
+ int tn = 0;
+ double[] ortpar = new double[5];
+ List<Track> vtxCandidate = new ArrayList<Track>();
+ for (CheatedBaseTrack track : tracks) {// 4
+ double[] opar = track.getTrackParameters();
+ for (int i = 0; i < 5; i++) {// 5
+ ortpar[i] = opar[i];
+ }// 4
+ MCParticle mcp = track.getMCParticle();
+ int pdgid = -1;
+ double mfmo = 0.;
+ double mcpmas = 0.;
+ if (mcp == null) {
+ System.out.println("No Particle");
+ continue;
+ }
+ Track fastMCTrack = factory.getTrack(new SpaceVector(mcp.getMomentum()), new SpacePoint(mcp.getOrigin()), new SpacePoint(), (int)mcp.getCharge(), new Random(), true);
+ Track fastMCTRackUnsmeared = factory.getTrack(new SpaceVector(mcp.getMomentum()), new SpacePoint(mcp.getOrigin()), new SpacePoint(), (int)mcp.getCharge(), new Random(), false);
+ Hep3Vector po = mcp.getOrigin();
+ Hep3Vector pm = mcp.getMomentum();
+ HelixSwimmer new_swimmer = new HelixSwimmer(5);
+ new_swimmer.setTrack(fastMCTrack);
+ double tracklength = new_swimmer.getTrackLengthToPoint(new SpacePoint(po));
+ SpaceVector newp = new_swimmer.getMomentumAtLength(tracklength);
+ aida.cloud1D("NewFastMCTrack/px").fill(newp.x());
+ aida.cloud1D("NewFastMCTrack/py").fill(newp.y());
+ aida.cloud1D("NewFastMCTrack/pz").fill(newp.z());
+
+ aida.cloud1D("NewFastMCTrack/px_partpx").fill(newp.x()-pm.x());
+ aida.cloud1D("NewFastMCTrack/py_partpy").fill(newp.y()-pm.y());
+ aida.cloud1D("NewFastMCTrack/pz_partpz").fill(newp.z()-pm.z());
+ aida.cloud1D("NewFastMCTrack/d0 unsmeared").fill(fastMCTRackUnsmeared.getParameters().getValues()[0]);
+ aida.cloud1D("NewFastMCTrack/phi0 unsmeared").fill(fastMCTRackUnsmeared.getParameters().getValues()[1]);
+ aida.cloud1D("NewFastMCTrack/omega unsmeared").fill(fastMCTRackUnsmeared.getParameters().getValues()[2]);
+ aida.cloud1D("NewFastMCTrack/z0 unsmeared").fill(fastMCTRackUnsmeared.getParameters().getValues()[3]);
+ aida.cloud1D("NewFastMCTrack/tanLambda unsmeared").fill(fastMCTRackUnsmeared.getParameters().getValues()[4]);
+ aida.cloud1D("NewFastMCTrack/d0 smeared").fill(fastMCTrack.getParameters().getValues()[0]);
+ aida.cloud1D("NewFastMCTrack/phi0 smeared").fill(fastMCTrack.getParameters().getValues()[1]);
+ aida.cloud1D("NewFastMCTrack/omega smeared").fill(fastMCTrack.getParameters().getValues()[2]);
+ aida.cloud1D("NewFastMCTrack/z0 smeared").fill(fastMCTrack.getParameters().getValues()[3]);
+ aida.cloud1D("NewFastMCTrack/tanLambda smeared").fill(fastMCTrack.getParameters().getValues()[4]);
+
+ double mphi0 = Math.atan2(pm.y(), pm.x());
+ double mPt = Math.sqrt(pm.x() * pm.x() + pm.y() * pm.y());
+ mfmo = Math.sqrt(mPt * mPt + pm.z() * pm.z());
+ double mtlam = pm.z() / mPt;
+ pdgid = mcp.getType().getPDGID();
+ mcpmas = mcp.getMass();
+ double mcpip = Math.sqrt(po.x() * po.x() + po.y() * po.y() + po.z()
+ * po.z());
+ int iq = (int) mcp.getCharge();
+ double[] tm = track.getMomentum();
+ double Pt = Math.sqrt(tm[0] * tm[0] + tm[1] * tm[1]);
+ double d0 = Math.abs(track.getTrackParameter(0));
+ double phi0 = track.getTrackParameter(1);
+ double z0 = track.getTrackParameter(3);
+ double tl = track.getTrackParameter(4);
+ double momega = track.getTrackParameter(2);
+// System.out.printf("PT: %.3f mcpip: %.3f tl: %.3f ip: %d", Pt,
+// mcpip, tl, iq);
+ // if ((Pt > minPt) && (mcpip < 0.001) && (Math.abs(tl) < maxTl)
+ // && (iq != 0) && (Pt < maxPt) && (Math.abs(tl) >= minTl)) {
+ if (Pt > 1.0) {
+ ttn++;
+ fitter.setTrack(track);
+ fitter.fit(); // **************************************************
+// System.out.println(track.fitSuccess());
+ if (!track.fitSuccess()) {
+// System.out.println("No fit success");
+ continue;
+ }
+
+ // hist
+
+ double[] ftpar = track.getTrackParameters();
+ double ftl = ftpar[4];
+ aida.cloud2D("am. of mat vs TanLam").fill(Math.abs(ftl),
+ trtrlpa * 100.);
+ double chi2 = track.getChi2();
+ int ndf = track.getNDF();
+ SymmetricMatrix trcovm = track.getErrorMatrix();
+ double dd0exp = Math.sqrt(trcovm.diagonal(0));
+ double dphi0exp = Math.sqrt(trcovm.diagonal(1));
+ double domexp = Math.sqrt(trcovm.diagonal(2));
+ boolean badFit = false;
+ if (domexp < 0.00000001) {
+// System.out.println("domexp is too small");
+ badFit = true;
+ }
+ if (Double.isNaN(domexp))
+ badFit = true;
+ if (badFit) {
+// System.out.println("bad fit");
+ continue;
+ }
+ double dz0exp = Math.sqrt(trcovm.diagonal(3));
+ double dtlexp = Math.sqrt(trcovm.diagonal(4));
+ double chi2aft = chi2 / ndf;
+// System.out.println("Event:" + evnum
+// + "covariance matrix omega error is: " + domexp
+// + " divided by omega: " + Math.abs(domexp / ortpar[2]));
+ ngdft++;
+ aida.cloud1D("Tan Lambda of fitted tracks").fill(ftpar[4]);
+ aida.cloud1D("Phi of fitted tracks").fill(ftpar[1]);
+ double[] fastMCTrackParams = fastMCTrack.getParameters().getValues();
+ aida.cloud1D("Comparison/d0").fill(fastMCTrackParams[0] - ftpar[0]);
+ aida.cloud1D("Comparison/phi").fill(fastMCTrackParams[1] - ftpar[1]);
+ aida.cloud1D("Comparison/omega").fill(fastMCTrackParams[2] - ftpar[2]);
+ aida.cloud1D("Comparison/z0").fill(fastMCTrackParams[3] - ftpar[3]);
+ aida.cloud1D("Comparison/tanLambda").fill(fastMCTrackParams[4] - ftpar[4]);
+
+ SymmetricMatrix fastMCErrorMatrix = fastMCTrack.getErrorMatrix();
+ aida.cloud1D("Comparison/d0 Error").fill(Math.sqrt(fastMCErrorMatrix.diagonal(0)) - dd0exp);
+ aida.cloud1D("Comparison/phi Error").fill(Math.sqrt(fastMCErrorMatrix.diagonal(1)) - dphi0exp);
+ aida.cloud1D("Comparison/omega Error").fill(Math.sqrt(fastMCErrorMatrix.diagonal(2)) - domexp);
+ aida.cloud1D("Comparison/z0 Error").fill(Math.sqrt(fastMCErrorMatrix.diagonal(3)) - dz0exp);
+ aida.cloud1D("Comparison/tanLambda Error").fill(Math.sqrt(fastMCErrorMatrix.diagonal(4)) - dtlexp);
+ org.lcsim.util.swim.HelixSwimmer old_swimmer = new org.lcsim.util.swim.HelixSwimmer(5);
+ old_swimmer.setTrack(track);
+ double tracklen = old_swimmer.getDistanceToPoint(po);
+ Hep3Vector oldp = ((Helix)old_swimmer.getTrajectory()).getTangentAtDistance(tracklen);
+ double m = pm.magnitude();
+ aida.cloud1D("FittedTrack diff px").fill(m*oldp.x() - pm.x());
+ aida.cloud1D("FittedTrack diff py").fill(m*oldp.y() - pm.y());
+ aida.cloud1D("FittedTrack diff pz").fill(m*oldp.z() - pm.z());
+ aida.cloud1D("FittedTrack unswum diff px").fill(track.getPX() - pm.x());
+ aida.cloud1D("FittedTrack unswum diff py").fill(track.getPY() - pm.y());
+ aida.cloud1D("FittedTrack unswum diff pz").fill(track.getPZ() - pm.z());
+
+ double[] fitres = new double[5];
+ for (int i = 0; i < 5; i++) {
+ fitres[i] = ftpar[i] - ortpar[i];
+ }
+ if (fitres[1] > Math.PI)
+ fitres[1] -= 2. * Math.PI;
+ if (fitres[1] < -Math.PI)
+ fitres[1] += 2. * Math.PI;
+ double[] tmf = track.getMomentum();
+ double Ptf = field / Math.abs(ftpar[2]);
+ double bfom = ftpar[2];
+ double ftP = Ptf * Math.sqrt(1. + ftpar[4] * ftpar[4]);
+ double relPer = (ftP - mfmo) / mfmo;
+ aida.cloud1D(
+ "delta full mom over full mom for fitted trks, Tl "
+ + df.format(minTl) + " - " + df.format(maxTl))
+ .fill(relPer);
+ double d0res = fitres[0];
+ aida.cloud1D("number of degrees of freedom for fitted tracks")
+ .fill(ndf);
+ aida.cloud1D(
+ "d0 residual for fitted tracks, Tl " + df.format(minTl)
+ + " - " + df.format(maxTl)).fill(d0res);
+ aida.cloud1D(
+ "phi0 residual for fitted tracks, Tl "
+ + df.format(minTl) + " - " + df.format(maxTl))
+ .fill(fitres[1]);
+ double domovom = (fitres[2]) / Math.abs(ortpar[2]);
+ aida.cloud1D(
+ "domega over omega for fitted tracks, Tl "
+ + df.format(minTl) + " - " + df.format(maxTl))
+ .fill(domovom);
+ domovom /= Pt;
+ aida.cloud1D(
+ "dPt over Pt**2 for fitted traks, Tl "
+ + df.format(minTl) + " - " + df.format(maxTl))
+ .fill(domovom);
+ double noromr = fitres[0] / dd0exp;
+ aida.cloud1D(
+ "normalised to cov.mat d0 residual, Tl "
+ + df.format(minTl) + " - " + df.format(maxTl))
+ .fill(noromr);
+ aida.cloud2D("Track cov.mat. exp d0 error vs tanlam").fill(
+ Math.abs(ftl), dd0exp);
+ noromr = fitres[1] / dphi0exp;
+ if (Math.abs(noromr) < 5.) {
+ aida.cloud1D(
+ "normalised to cov.mat Phi0 residual, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(noromr);
+ }
+ aida.cloud2D("Track cov.mat. exp phi0 error vs tanlam").fill(
+ Math.abs(ftl), dphi0exp);
+ noromr = fitres[2] / domexp;
+ if (Math.abs(noromr) < 5.) {
+ aida.cloud1D(
+ "normalised to cov.mat omega residual, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(noromr);
+ aida
+ .cloud1D(
+ "normalised to cov.mat omega residual for all EC tracks")
+ .fill(noromr);
+ }
+ aida.cloud2D("Track cov.mat. rel. omega error vs tanlam").fill(
+ Math.abs(ftl), domexp / Math.abs(ortpar[2]));
+ noromr = fitres[3] / dz0exp;
+ if (Math.abs(noromr) < 5.) {
+ aida.cloud1D(
+ "normalised to cov.mat z0 residual, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(noromr);
+ aida.cloud1D(
+ "z0 residual for fitted tracks, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(fitres[3]);
+ }
+ aida.cloud2D("Track cov.mat. exp z0 error vs tanlam").fill(
+ Math.abs(ftl), dz0exp);
+ noromr = fitres[4] / dtlexp;
+ if (Math.abs(noromr) < 5.) {
+ aida.cloud1D(
+ "normalised to cov.mat Tanlam residual, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(noromr);
+ aida.cloud1D(
+ "tanl residual for fitted tracks, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(fitres[4]);
+ }
+ aida.cloud2D("Track cov.mat. exp tanlam error vs tanlam").fill(
+ Math.abs(ftl), dtlexp);
+ noromr = chi2 / ndf;
+ if (noromr < 10.) {
+ aida.cloud1D(
+ "Chi2 per dof for fitted tracks, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(chi2 / ndf);
+ }
+ aida.cloud1D(
+ "Pt of fitted tracks, Tl " + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(mPt);
+ aida.cloud1D(
+ "tan lambda of fitted tracks, Tl " + df.format(minTl)
+ + " - " + df.format(maxTl)).fill(tl);
+ if ((chi2 / ndf > 50.) || (chi2 < 0.))
+// System.out.println(" Event " + evnum + " track " + tn
+// + " has chi2: " + chi2);
+ tn++;
+ vtxCandidate.add(new NewTrack(track));
+ }
+ }
+ evnum++;
+ }
+}
lcsim/src/org/lcsim/contrib/JanStrube/standalone
diff -u -r1.6 -r1.7
--- MainLoop.java 14 Oct 2006 08:45:21 -0000 1.6
+++ MainLoop.java 7 Dec 2006 06:11:24 -0000 1.7
@@ -23,19 +23,21 @@
LCSimLoop loop = new LCSimLoop();
// URL location = new URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/ILC500/Zgamma/stdhep/pythia/pythiaZgamma.stdhep");
// URL location = new URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/stdhep/K0S_pipi_Theta45-135_5-25Gev.stdhep");
- URL location = new URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/stdhep/psi_mumu_Theta4-176_5-100GeV.stdhep");
+ //URL location = new URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/stdhep/psi_mumu_Theta4-176_5-100GeV.stdhep");
// URL location = new URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/ILC500/mumu/stdhep/pandorapythia/panpymumu.stdhep");
+ URL location = new URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/sidaug05/slcio/slic/psi_mumu_Theta4-176_5-100GeV_SLIC_v1r9p3_sidaug05.slcio");
FileCache cache = new FileCache();
File trackFile = cache.getCachedFile(location);
- loop.setStdhepRecordSource(trackFile, "sidaug05");
-// loop.setLCIORecordSource(input);
+// loop.setStdhepRecordSource(trackFile, "sidaug05");
+ loop.setLCIORecordSource(trackFile);
// loop.add(new MCFast(false, false));
- loop.add(new FitterTestDriver());
+// loop.add(new FitterTestDriver());
// loop.add(new FitterTestDriver_ReconTrack());
// File output = new File("exampleAnalysisJava.slcio");
// loop.add(new LCIODriver(output));
// System.out.println("starting:");
// loop.add(new TrackComparison());
+ loop.add(new NickTrack_FastMCTrackComparison());
loop.loop(-1);
loop.dispose();
AIDA.defaultInstance().saveAs("exampleAnalysisJava.aida");
lcsim/src/org/lcsim/contrib/JanStrube/standalone
diff -u -r1.5 -r1.6
--- MainLoop.py 29 Nov 2006 02:32:14 -0000 1.5
+++ MainLoop.py 7 Dec 2006 06:11:24 -0000 1.6
@@ -23,14 +23,14 @@
loop = LCSimLoop()
# location = URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/ILC500/Zgamma/stdhep/pythia/pythiaZgamma.stdhep")
# location = URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/stdhep/psi_mumu_Theta4-176_5-100GeV.stdhep")
- location = URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/stdhep/K0S_pipi_Theta45-135_5-25Gev.stdhep")
+# location = URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/stdhep/K0S_pipi_Theta45-135_5-25Gev.stdhep")
# location = URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/sidaug05/slcio/slic/tau_5pi_Theta20-90_10-200GeV_SLIC_v1r9p3_sidaug05.slcio")
# location = URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/sidaug05/slcio/slic/tau_3pi_Theta20-90_10-200GeV_SLIC_v1r9p3_sidaug05.slcio")
-# location = URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/sidaug05/slcio/slic/psi_mumu_Theta4-176_5-100GeV_SLIC_v1r9p3_sidaug05.slcio")
+ location = URL("ftp://ftp-lcd.slac.stanford.edu/lcd/ILC/singleParticle/sidaug05/slcio/slic/psi_mumu_Theta4-176_5-100GeV_SLIC_v1r9p3_sidaug05.slcio")
cache = FileCache()
trackFile = cache.getCachedFile(location)
- loop.setStdhepRecordSource(trackFile, "sidaug05")
-# loop.setLCIORecordSource(trackFile)
+# loop.setStdhepRecordSource(trackFile, "sidaug05")
+ loop.setLCIORecordSource(trackFile)
# no beamspotConstraint, no simple smearing
loop.add(MCFast(False, False))
# loop.add(FitterTestDriver())
lcsim/src/org/lcsim/contrib/JanStrube/standalone
diff -u -r1.1 -r1.2
--- NickTrackVtxFitter.java 28 Oct 2006 01:11:33 -0000 1.1
+++ NickTrackVtxFitter.java 7 Dec 2006 06:11:24 -0000 1.2
@@ -7,8 +7,6 @@
import org.lcsim.geometry.Detector;
import org.lcsim.contrib.JanStrube.tracking.NewTrack;
import org.lcsim.contrib.JanStrube.tracking.Track;
-import org.lcsim.contrib.JanStrube.vtxFitter.Fitter;
-import org.lcsim.contrib.JanStrube.vtxFitter.Vertex;
import org.lcsim.contrib.NickSinev.tracking.wmfitter.*;
import org.lcsim.contrib.NickSinev.tracking.util.*;
import org.lcsim.spacegeom.SpacePoint;
@@ -27,321 +25,300 @@
import Jama.*;
public class NickTrackVtxFitter extends Driver {
- int firstev = 0;
-
- int lastev = 9999;
-
- boolean real = true;
-
- double minPt = 0.1;
-
- double maxPt = 500.;
-
- double minTl = 0.;
-
- double maxTl = 5.2;
-
- double[] minTllms = { 0.01, 0.4, 0.9, 1.4, 2.2, 3.2 };
-
- double[] maxTllms = { 0.4, 0.9, 1.4, 2.2, 3.2, 5. };
-
- private double vtxbrphires = 0.0035;
-
- private double vtxbzres = 0.035;
-
- private double vtxecrpres = 0.0035;
-
- private double vtxecrres = 0.0035;
-
- private double trkbrpres = 0.007;
-
- private double trkbzres = 30.;
-
- private double trkecrpres = 0.01;
-
- private double trkecrres = 0.01;
-
- private SLDFitTrack fitter = new SLDFitTrack();
-
- private DecimalFormat df = new DecimalFormat();
-
- private AIDA aida = AIDA.defaultInstance();
-
- private SmearTrackerHits smTh = new SmearTrackerHits();
-
- private TrackingCheaterSmHits cheat = new TrackingCheaterSmHits();
-
- private int evnum = 0;
-
- double field = 5. * Constants.fieldConversion;
-
- int ttn = 0;
-
- int ngdft = 0;
-
- double trtrlpa = 0.;
-
- public NickTrackVtxFitter() {
- if (real) {
- smTh.setVtxBarrRPResolution(vtxbrphires);
- smTh.setVtxBarrZResolution(vtxbzres);
- smTh.setVtxECRPResolution(vtxecrpres);
- smTh.setVtxECRResolution(vtxecrres);
- smTh.setTrkBarrRPResolution(trkbrpres);
- smTh.setTrkBarrZResolution(trkbzres);
- smTh.setTrkECRPResolution(trkecrpres);
- smTh.setTrkECRResolution(trkecrres);
- }
- cheat.useChargedOnly = true;
- smTh.setDoSmear(true); // ++++++++++++++++++++++++++++++++++turn on
- // smearing
- add(smTh);
- add(cheat);
- df.setMaximumFractionDigits(2);
- fitter.setDoHist(false);
- }
-
- public void process(EventHeader event) {
- if (!((evnum > firstev - 1) && (evnum < lastev + 1)))
- return;
- if (evnum % 100 == 0) {
- System.out.println("Processing event " + evnum
- + " accepted chit tracks " + ttn + " fitted " + ngdft);
- }
- super.process(event);
- if (evnum == firstev) {
- Detector det = event.getDetector();
- fitter.setDetector(det);
- }
- List<CheatedBaseTrack> tracks = null;
- try {
- tracks = event.get(CheatedBaseTrack.class, "CheatedBaseTracks");
- } catch (Exception trnotf) {
- System.out.println("No CheatedBaseTracks in the event ! \n");
- } finally {
- }
- if (tracks == null)
- return;
-
- int nprts = tracks.size();
- System.out.println("Tracking cheater has found " + nprts
- + " BaseTracks in event " + evnum);
- int tn = 0;
- double[] ortpar = new double[5];
- List<Track> vtxCabdidate = new ArrayList<Track>();
- for (CheatedBaseTrack track : tracks) {// 4
- double[] opar = track.getTrackParameters();
- for (int i = 0; i < 5; i++) {// 5
- ortpar[i] = opar[i];
- }// 4
- MCParticle mcp = track.getMCParticle();
- Hep3Vector po = new BasicHep3Vector(1., 1., 1.);
- int pdgid = -1;
- double mfmo = 0.;
- double mcpmas = 0.;
- if (mcp == null) {
- System.out.println("No Particle");
- continue;
- }
- po = mcp.getOrigin();
- Hep3Vector pm = mcp.getMomentum();
- double mphi0 = Math.atan2(pm.y(), pm.x());
- double mPt = Math.sqrt(pm.x() * pm.x() + pm.y() * pm.y());
- mfmo = Math.sqrt(mPt * mPt + pm.z() * pm.z());
- double mtlam = pm.z() / mPt;
- pdgid = mcp.getType().getPDGID();
- mcpmas = mcp.getMass();
- double mcpip = Math.sqrt(po.x() * po.x() + po.y() * po.y() + po.z()
- * po.z());
- int iq = (int) mcp.getCharge();
- double[] tm = track.getMomentum();
- double Pt = Math.sqrt(tm[0] * tm[0] + tm[1] * tm[1]);
- double d0 = Math.abs(track.getTrackParameter(0));
- double phi0 = track.getTrackParameter(1);
- double z0 = track.getTrackParameter(3);
- double tl = track.getTrackParameter(4);
- double momega = track.getTrackParameter(2);
- System.out.printf("PT: %.3f mcpip: %.3f tl: %.3f ip: %d", Pt,
- mcpip, tl, iq);
- // if ((Pt > minPt) && (mcpip < 0.001) && (Math.abs(tl) < maxTl)
- // && (iq != 0) && (Pt < maxPt) && (Math.abs(tl) >= minTl)) {
- if (Pt > 1.0) {
- ttn++;
- fitter.setTrack(track);
- fitter.fit(); // **************************************************
- System.out.println(track.fitSuccess());
- if (!track.fitSuccess()) {
- System.out.println("No fit success");
- continue;
- }
-
- // hist
+ int firstev = 0;
+ int lastev = 9999;
+ boolean real = true;
+ double minPt = 0.1;
+ double maxPt = 500.;
+ double minTl = 0.;
+ double maxTl = 5.2;
+ double[] minTllms = { 0.01, 0.4, 0.9, 1.4, 2.2, 3.2 };
+ double[] maxTllms = { 0.4, 0.9, 1.4, 2.2, 3.2, 5. };
+ private double vtxbrphires = 0.0035;
+ private double vtxbzres = 0.035;
+ private double vtxecrpres = 0.0035;
+ private double vtxecrres = 0.0035;
+ private double trkbrpres = 0.007;
+ private double trkbzres = 30.;
+ private double trkecrpres = 0.01;
+ private double trkecrres = 0.01;
+ private SLDFitTrack fitter = new SLDFitTrack();
+ private DecimalFormat df = new DecimalFormat();
+ private AIDA aida = AIDA.defaultInstance();
+ private SmearTrackerHits smTh = new SmearTrackerHits();
+ private TrackingCheaterSmHits cheat = new TrackingCheaterSmHits();
+ private int evnum = 0;
+ double field = 5. * Constants.fieldConversion;
+ int ttn = 0;
+ int ngdft = 0;
+ double trtrlpa = 0.;
+
+ public NickTrackVtxFitter() {
+ if (real) {
+ smTh.setVtxBarrRPResolution(vtxbrphires);
+ smTh.setVtxBarrZResolution(vtxbzres);
+ smTh.setVtxECRPResolution(vtxecrpres);
+ smTh.setVtxECRResolution(vtxecrres);
+ smTh.setTrkBarrRPResolution(trkbrpres);
+ smTh.setTrkBarrZResolution(trkbzres);
+ smTh.setTrkECRPResolution(trkecrpres);
+ smTh.setTrkECRResolution(trkecrres);
+ }
+ cheat.useChargedOnly = true;
+ smTh.setDoSmear(true); // ++++++++++++++++++++++++++++++++++turn on
+ // smearing
+ add(smTh);
+ add(cheat);
+ df.setMaximumFractionDigits(2);
+ fitter.setDoHist(false);
+ }
+
+ public void process(EventHeader event) {
+ if (!((evnum > firstev - 1) && (evnum < lastev + 1)))
+ return;
+ if (evnum % 100 == 0) {
+ System.out.println("Processing event " + evnum
+ + " accepted chit tracks " + ttn + " fitted " + ngdft);
+ }
+ super.process(event);
+ if (evnum == firstev) {
+ Detector det = event.getDetector();
+ fitter.setDetector(det);
+ }
+ List<CheatedBaseTrack> tracks = null;
+ try {
+ tracks = event.get(CheatedBaseTrack.class, "CheatedBaseTracks");
+ } catch (Exception trnotf) {
+ System.out.println("No CheatedBaseTracks in the event ! \n");
+ } finally {
+ }
+ if (tracks == null)
+ return;
- double[] ftpar = track.getTrackParameters();
- double ftl = ftpar[4];
- aida.cloud2D("am. of mat vs TanLam").fill(Math.abs(ftl),
- trtrlpa * 100.);
- double chi2 = track.getChi2();
- int ndf = track.getNDF();
- SymmetricMatrix trcovm = track.getErrorMatrix();
- double dd0exp = Math.sqrt(trcovm.diagonal(0));
- double dphi0exp = Math.sqrt(trcovm.diagonal(1));
- double domexp = Math.sqrt(trcovm.diagonal(2));
- boolean badFit = false;
- if (domexp < 0.00000001) {
- System.out.println("domexp is too small");
- badFit = true;
- }
- if (Double.toString(domexp).equals("NaN"))
- badFit = true;
- double dz0exp = Math.sqrt(trcovm.diagonal(3));
- double dtlexp = Math.sqrt(trcovm.diagonal(4));
- double chi2aft = chi2 / ndf;
- if (badFit) {
- System.out.println("bad fit");
- continue;
- }
-// System.out.println("Event:" + evnum
-// + "covariance matrix omega error is: " + domexp
-// + " divaded by omega: " + Math.abs(domexp / ortpar[2]));
- ngdft++;
-// aida.cloud1D("Tan Lambda of fitted tracks").fill(ftpar[4]);
-// aida.cloud1D("Phi of fitted tracks").fill(ftpar[1]);
- double[] fitres = new double[5];
- for (int i = 0; i < 5; i++) {
- fitres[i] = ftpar[i] - ortpar[i];
- }
- if (fitres[1] > Math.PI)
- fitres[1] -= 2. * Math.PI;
- if (fitres[1] < -Math.PI)
- fitres[1] += 2. * Math.PI;
- double[] tmf = track.getMomentum();
- double Ptf = field / Math.abs(ftpar[2]);
- double bfom = ftpar[2];
- double ftP = Ptf * Math.sqrt(1. + ftpar[4] * ftpar[4]);
- double relPer = (ftP - mfmo) / mfmo;
-// aida.cloud1D(
-// "delta full mom over full mom for fitted trks, Tl "
-// + df.format(minTl) + " - " + df.format(maxTl))
-// .fill(relPer);
- double d0res = fitres[0];
-// aida.cloud1D("number of degrees of freedom for fitted tracks")
-// .fill(ndf);
-// aida.cloud1D(
-// "d0 residual for fitted tracks, Tl " + df.format(minTl)
-// + " - " + df.format(maxTl)).fill(d0res);
-// aida.cloud1D(
-// "phi0 residual for fitted tracks, Tl "
-// + df.format(minTl) + " - " + df.format(maxTl))
-// .fill(fitres[1]);
- double domovom = (fitres[2]) / Math.abs(ortpar[2]);
-// aida.cloud1D(
-// "domega over omega for fitted tracks, Tl "
-// + df.format(minTl) + " - " + df.format(maxTl))
-// .fill(domovom);
- domovom /= Pt;
-// aida.cloud1D(
-// "dPt over Pt**2 for fitted traks, Tl "
-// + df.format(minTl) + " - " + df.format(maxTl))
-// .fill(domovom);
- double noromr = fitres[0] / dd0exp;
-// aida.cloud1D(
-// "normalised to cov.mat d0 residual, Tl "
-// + df.format(minTl) + " - " + df.format(maxTl))
-// .fill(noromr);
-// aida.cloud2D("Track cov.mat. exp d0 error vs tanlam").fill(
-// Math.abs(ftl), dd0exp);
- noromr = fitres[1] / dphi0exp;
-// if (Math.abs(noromr) < 5.) {
-// aida.cloud1D(
-// "normalised to cov.mat Phi0 residual, Tl "
-// + df.format(minTl) + " - "
-// + df.format(maxTl)).fill(noromr);
-// }
-// aida.cloud2D("Track cov.mat. exp phi0 error vs tanlam").fill(
-// Math.abs(ftl), dphi0exp);
- noromr = fitres[2] / domexp;
-// if (Math.abs(noromr) < 5.) {
-// aida.cloud1D(
-// "normalised to cov.mat omega residual, Tl "
-// + df.format(minTl) + " - "
-// + df.format(maxTl)).fill(noromr);
-// aida.cloud1D("normalised to cov.mat omega residual for all EC tracks")
-// .fill(noromr);
-// }
-// aida.cloud2D("Track cov.mat. rel. omega error vs tanlam").fill(
-// Math.abs(ftl), domexp / Math.abs(ortpar[2]));
- noromr = fitres[3] / dz0exp;
-// if (Math.abs(noromr) < 5.) {
-// aida.cloud1D(
-// "normalised to cov.mat z0 residual, Tl "
-// + df.format(minTl) + " - "
-// + df.format(maxTl)).fill(noromr);
-// aida.cloud1D(
-// "z0 residual for fitted tracks, Tl "
-// + df.format(minTl) + " - "
-// + df.format(maxTl)).fill(fitres[3]);
-// }
-// aida.cloud2D("Track cov.mat. exp z0 error vs tanlam").fill(
-// Math.abs(ftl), dz0exp);
- noromr = fitres[4] / dtlexp;
-// if (Math.abs(noromr) < 5.) {
-// aida.cloud1D(
-// "normalised to cov.mat Tanlam residual, Tl "
-// + df.format(minTl) + " - "
-// + df.format(maxTl)).fill(noromr);
-// aida.cloud1D(
-// "tanl residual for fitted tracks, Tl "
-// + df.format(minTl) + " - "
-// + df.format(maxTl)).fill(fitres[4]);
-// }
-// aida.cloud2D("Track cov.mat. exp tanlam error vs tanlam").fill(
-// Math.abs(ftl), dtlexp);
- noromr = chi2 / ndf;
-// if (noromr < 10.) {
-// aida.cloud1D(
-// "Chi2 per dof for fitted tracks, Tl "
-// + df.format(minTl) + " - "
-// + df.format(maxTl)).fill(chi2 / ndf);
-// }
-// aida.cloud1D(
-// "Pt of fitted tracks, Tl " + df.format(minTl) + " - "
-// + df.format(maxTl)).fill(mPt);
-// aida.cloud1D(
-// "tan lambda of fitted tracks, Tl " + df.format(minTl)
-// + " - " + df.format(maxTl)).fill(tl);
-// if ((chi2 / ndf > 50.) || (chi2 < 0.))
-// System.out.println(" Event " + evnum + " track " + tn
-// + " has chi2: " + chi2);
- tn++;
- vtxCabdidate.add(new NewTrack(track));
-
- }
- }
+ int nprts = tracks.size();
+ System.out.println("Tracking cheater has found " + nprts
+ + " BaseTracks in event " + evnum);
+ int tn = 0;
+ double[] ortpar = new double[5];
+ List<Track> vtxCandidate = new ArrayList<Track>();
+ for (CheatedBaseTrack track : tracks) {// 4
+ double[] opar = track.getTrackParameters();
+ for (int i = 0; i < 5; i++) {// 5
+ ortpar[i] = opar[i];
+ }// 4
+ MCParticle mcp = track.getMCParticle();
+ Hep3Vector po = new BasicHep3Vector(1., 1., 1.);
+ int pdgid = -1;
+ double mfmo = 0.;
+ double mcpmas = 0.;
+ if (mcp == null) {
+ System.out.println("No Particle");
+ continue;
+ }
+ po = mcp.getOrigin();
+ Hep3Vector pm = mcp.getMomentum();
+ double mphi0 = Math.atan2(pm.y(), pm.x());
+ double mPt = Math.sqrt(pm.x() * pm.x() + pm.y() * pm.y());
+ mfmo = Math.sqrt(mPt * mPt + pm.z() * pm.z());
+ double mtlam = pm.z() / mPt;
+ pdgid = mcp.getType().getPDGID();
+ mcpmas = mcp.getMass();
+ double mcpip = Math.sqrt(po.x() * po.x() + po.y() * po.y() + po.z()
+ * po.z());
+ int iq = (int) mcp.getCharge();
+ double[] tm = track.getMomentum();
+ double Pt = Math.sqrt(tm[0] * tm[0] + tm[1] * tm[1]);
+ double d0 = Math.abs(track.getTrackParameter(0));
+ double phi0 = track.getTrackParameter(1);
+ double z0 = track.getTrackParameter(3);
+ double tl = track.getTrackParameter(4);
+ double momega = track.getTrackParameter(2);
+ System.out.printf("PT: %.3f mcpip: %.3f tl: %.3f ip: %d", Pt,
+ mcpip, tl, iq);
+ // if ((Pt > minPt) && (mcpip < 0.001) && (Math.abs(tl) < maxTl)
+ // && (iq != 0) && (Pt < maxPt) && (Math.abs(tl) >= minTl)) {
+ if (Pt > 1.0) {
+ ttn++;
+ fitter.setTrack(track);
+ fitter.fit(); // **************************************************
+ System.out.println(track.fitSuccess());
+ if (!track.fitSuccess()) {
+ System.out.println("No fit success");
+ continue;
+ }
+
+ // hist
+
+ double[] ftpar = track.getTrackParameters();
+ double ftl = ftpar[4];
+ aida.cloud2D("am. of mat vs TanLam").fill(Math.abs(ftl),
+ trtrlpa * 100.);
+ double chi2 = track.getChi2();
+ int ndf = track.getNDF();
+ SymmetricMatrix trcovm = track.getErrorMatrix();
+ double dd0exp = Math.sqrt(trcovm.diagonal(0));
+ double dphi0exp = Math.sqrt(trcovm.diagonal(1));
+ double domexp = Math.sqrt(trcovm.diagonal(2));
+ boolean badFit = false;
+ if (domexp < 0.00000001) {
+ System.out.println("domexp is too small");
+ badFit = true;
+ }
+ if (Double.toString(domexp).equals("NaN"))
+ badFit = true;
+ double dz0exp = Math.sqrt(trcovm.diagonal(3));
+ double dtlexp = Math.sqrt(trcovm.diagonal(4));
+ double chi2aft = chi2 / ndf;
+ if (badFit) {
+ System.out.println("bad fit");
+ continue;
+ }
+ System.out.println("Event:" + evnum
+ + "covariance matrix omega error is: " + domexp
+ + " divaded by omega: " + Math.abs(domexp / ortpar[2]));
+ ngdft++;
+ aida.cloud1D("Tan Lambda of fitted tracks").fill(ftpar[4]);
+ aida.cloud1D("Phi of fitted tracks").fill(ftpar[1]);
+ double[] fitres = new double[5];
+ for (int i = 0; i < 5; i++) {
+ fitres[i] = ftpar[i] - ortpar[i];
+ }
+ if (fitres[1] > Math.PI)
+ fitres[1] -= 2. * Math.PI;
+ if (fitres[1] < -Math.PI)
+ fitres[1] += 2. * Math.PI;
+ double[] tmf = track.getMomentum();
+ double Ptf = field / Math.abs(ftpar[2]);
+ double bfom = ftpar[2];
+ double ftP = Ptf * Math.sqrt(1. + ftpar[4] * ftpar[4]);
+ double relPer = (ftP - mfmo) / mfmo;
+ aida.cloud1D(
+ "delta full mom over full mom for fitted trks, Tl "
+ + df.format(minTl) + " - " + df.format(maxTl))
+ .fill(relPer);
+ double d0res = fitres[0];
+ aida.cloud1D("number of degrees of freedom for fitted tracks")
+ .fill(ndf);
+ aida.cloud1D(
+ "d0 residual for fitted tracks, Tl " + df.format(minTl)
+ + " - " + df.format(maxTl)).fill(d0res);
+ aida.cloud1D(
+ "phi0 residual for fitted tracks, Tl "
+ + df.format(minTl) + " - " + df.format(maxTl))
+ .fill(fitres[1]);
+ double domovom = (fitres[2]) / Math.abs(ortpar[2]);
+ aida.cloud1D(
+ "domega over omega for fitted tracks, Tl "
+ + df.format(minTl) + " - " + df.format(maxTl))
+ .fill(domovom);
+ domovom /= Pt;
+ aida.cloud1D(
+ "dPt over Pt**2 for fitted traks, Tl "
+ + df.format(minTl) + " - " + df.format(maxTl))
+ .fill(domovom);
+ double noromr = fitres[0] / dd0exp;
+ aida.cloud1D(
+ "normalised to cov.mat d0 residual, Tl "
+ + df.format(minTl) + " - " + df.format(maxTl))
+ .fill(noromr);
+ aida.cloud2D("Track cov.mat. exp d0 error vs tanlam").fill(
+ Math.abs(ftl), dd0exp);
+ noromr = fitres[1] / dphi0exp;
+ if (Math.abs(noromr) < 5.) {
+ aida.cloud1D(
+ "normalised to cov.mat Phi0 residual, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(noromr);
+ }
+ aida.cloud2D("Track cov.mat. exp phi0 error vs tanlam").fill(
+ Math.abs(ftl), dphi0exp);
+ noromr = fitres[2] / domexp;
+ if (Math.abs(noromr) < 5.) {
+ aida.cloud1D(
+ "normalised to cov.mat omega residual, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(noromr);
+ aida
+ .cloud1D(
+ "normalised to cov.mat omega residual for all EC tracks")
+ .fill(noromr);
+ }
+ aida.cloud2D("Track cov.mat. rel. omega error vs tanlam").fill(
+ Math.abs(ftl), domexp / Math.abs(ortpar[2]));
+ noromr = fitres[3] / dz0exp;
+ if (Math.abs(noromr) < 5.) {
+ aida.cloud1D(
+ "normalised to cov.mat z0 residual, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(noromr);
+ aida.cloud1D(
+ "z0 residual for fitted tracks, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(fitres[3]);
+ }
+ aida.cloud2D("Track cov.mat. exp z0 error vs tanlam").fill(
+ Math.abs(ftl), dz0exp);
+ noromr = fitres[4] / dtlexp;
+ if (Math.abs(noromr) < 5.) {
+ aida.cloud1D(
+ "normalised to cov.mat Tanlam residual, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(noromr);
+ aida.cloud1D(
+ "tanl residual for fitted tracks, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(fitres[4]);
+ }
+ aida.cloud2D("Track cov.mat. exp tanlam error vs tanlam").fill(
+ Math.abs(ftl), dtlexp);
+ noromr = chi2 / ndf;
+ if (noromr < 10.) {
+ aida.cloud1D(
+ "Chi2 per dof for fitted tracks, Tl "
+ + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(chi2 / ndf);
+ }
+ aida.cloud1D(
+ "Pt of fitted tracks, Tl " + df.format(minTl) + " - "
+ + df.format(maxTl)).fill(mPt);
+ aida.cloud1D(
+ "tan lambda of fitted tracks, Tl " + df.format(minTl)
+ + " - " + df.format(maxTl)).fill(tl);
+ if ((chi2 / ndf > 50.) || (chi2 < 0.))
+ System.out.println(" Event " + evnum + " track " + tn
+ + " has chi2: " + chi2);
+ tn++;
+ vtxCandidate.add(new NewTrack(track));
+ }
+ }
for (MCParticle part : event.getMCParticles()) {
- // K0: 310, J/Psi: 443
+ // K0: 310, J/Psi: 443
if (abs(part.getPDGID()) != 443)
continue;
-
- Fitter fitter = new Fitter(event);
- Vertex newVtx2 = fitter.fit(vtxCabdidate, new SpacePoint());
- double sigmaX = sqrt(newVtx2.getSpatialCovarianceMatrix().get(0, 0));
- double sigmaY = sqrt(newVtx2.getSpatialCovarianceMatrix().get(1, 1));
- double sigmaZ = sqrt(newVtx2.getSpatialCovarianceMatrix().get(2, 2));
- double deltaX = newVtx2.location().x() - part.getEndPoint().x();
- double deltaY = newVtx2.location().y() - part.getEndPoint().y();
- double deltaZ = newVtx2.location().z() - part.getEndPoint().z();
- aida.cloud1D("decayLength").fill(new SpacePoint(part.getEndPoint()).magnitude());
- aida.cloud1D("vtx_smeared_deltaX").fill(deltaX);
- aida.cloud1D("vtx_smeared_deltaY").fill(deltaY);
- aida.cloud1D("vtx_smeared_deltaZ").fill(deltaZ);
- aida.cloud1D("vtx_sigma_x").fill(sigmaX);
- aida.cloud1D("vtx_sigma_y").fill(sigmaY);
- aida.cloud1D("vtx_sigma_z").fill(sigmaZ);
- aida.cloud1D("vtx_pull_x").fill(deltaX/sigmaX);
- aida.cloud1D("vtx_pull_y").fill(deltaY/sigmaY);
- aida.cloud1D("vtx_pull_z").fill(deltaZ/sigmaZ);
+
+ // Fitter fitter = new Fitter(event);
+ // Vertex newVtx2 = fitter.fit(vtxCandidate, new SpacePoint());
+ // double sigmaX = sqrt(newVtx2.getSpatialCovarianceMatrix().get(0,
+ // 0));
+ // double sigmaY = sqrt(newVtx2.getSpatialCovarianceMatrix().get(1,
+ // 1));
+ // double sigmaZ = sqrt(newVtx2.getSpatialCovarianceMatrix().get(2,
+ // 2));
+ // double deltaX = newVtx2.location().x() - part.getEndPoint().x();
+ // double deltaY = newVtx2.location().y() - part.getEndPoint().y();
+ // double deltaZ = newVtx2.location().z() - part.getEndPoint().z();
+ // aida.cloud1D("decayLength").fill(
+ // new SpacePoint(part.getEndPoint()).magnitude());
+ // aida.cloud1D("vtx_smeared_deltaX").fill(deltaX);
+ // aida.cloud1D("vtx_smeared_deltaY").fill(deltaY);
+ // aida.cloud1D("vtx_smeared_deltaZ").fill(deltaZ);
+ // aida.cloud1D("vtx_sigma_x").fill(sigmaX);
+ // aida.cloud1D("vtx_sigma_y").fill(sigmaY);
+ // aida.cloud1D("vtx_sigma_z").fill(sigmaZ);
+ // aida.cloud1D("vtx_pull_x").fill(deltaX / sigmaX);
+ // aida.cloud1D("vtx_pull_y").fill(deltaY / sigmaY);
+ // aida.cloud1D("vtx_pull_z").fill(deltaZ / sigmaZ);
}
- evnum++;
- }
+ evnum++;
+ }
}
lcsim/src/org/lcsim/contrib/JanStrube/standalone
diff -u -r1.3 -r1.4
--- Tracks.py 29 Nov 2006 02:32:14 -0000 1.3
+++ Tracks.py 7 Dec 2006 06:11:24 -0000 1.4
@@ -7,6 +7,7 @@
from org.lcsim.contrib.JanStrube.tracking import HelixSwimmer as Swimmer_new
from math import sqrt
from java.lang import Boolean
+from java.util import Random
false = Boolean('false')
true = Boolean('true')
@@ -43,6 +44,7 @@
p = momentum
c = int(particle.getCharge())
newT = trackFac.getTrack(CartesianVector(p.x(), p.y(), p.z()), CartesianPoint(o.x(), o.y(), o.z()), CartesianPoint(0, 0, 0), c)
+ newT_unsmeared = trackFac.getTrack(CartesianVector(p.x(), p.y(), p.z()), CartesianPoint(o.x(), o.y(), o.z()), CartesianPoint(0, 0, 0), c, Random(), false)
emap_new = newT.getParameters()
vals_old = emap_old.getValues()
self.aida.cloud1D('OldTrack/d0').fill(vals_old[0])
@@ -104,10 +106,10 @@
self.aida.cloud1D('OldTrack/diffOrigin y').fill(point_old.y() - origin.y())
self.aida.cloud1D('OldTrack/diffOrigin z').fill(point_old.z() - origin.z())
momentum_old = swimmer_old.getTrajectory().getTangentAtDistance(dist_old)
- #m = momentum.magnitude()
- #self.aida.cloud1D('OldTrack/diffMomentum x').fill(momentum.x() - m*momentum_old.x())
- #self.aida.cloud1D('OldTrack/diffMomentum y').fill(momentum.y() - m*momentum_old.y())
- #self.aida.cloud1D('OldTrack/diffMomentum z').fill(momentum.z() - m*momentum_old.z())
+ m = momentum.magnitude()
+ self.aida.cloud1D('OldTrack/diffMomentum x').fill(momentum.x() - m*momentum_old.x())
+ self.aida.cloud1D('OldTrack/diffMomentum y').fill(momentum.y() - m*momentum_old.y())
+ self.aida.cloud1D('OldTrack/diffMomentum z').fill(momentum.z() - m*momentum_old.z())
self.aida.cloud1D('OldTrack distance to Origin').fill(dist_old)
swimmer_new.setTrack(newT)
@@ -118,10 +120,17 @@
self.aida.cloud1D('NewTrack/diffOrigin x').fill(point_new.x() - origin.x())
self.aida.cloud1D('NewTrack/diffOrigin y').fill(point_new.y() - origin.y())
self.aida.cloud1D('NewTrack/diffOrigin z').fill(point_new.z() - origin.z())
- #momentum_new = swimmer_new.getMomentumAtLength(dist_new)
- #self.aida.cloud1D('NewTrack/diffMomentum x').fill(momentum.x() - momentum_new.x())
- #self.aida.cloud1D('NewTrack/diffMomentum y').fill(momentum.y() - momentum_new.y())
- #self.aida.cloud1D('NewTrack/diffMomentum z').fill(momentum.z() - momentum_new.z())
+ momentum_new = swimmer_new.getMomentumAtLength(dist_new)
+ swimmer_new.setTrack(newT_unsmeared)
+ dist_newUnsmeared = swimmer_new.getTrackLengthToPoint(CartesianPoint(origin.v()))
+ momentum_newUnsmeared = swimmer_new.getMomentumAtLength(dist_newUnsmeared)
+ self.aida.cloud1D('NewTrack/diffMomentum x').fill(momentum.x() - momentum_new.x())
+ self.aida.cloud1D('NewTrack/diffMomentum y').fill(momentum.y() - momentum_new.y())
+ self.aida.cloud1D('NewTrack/diffMomentum z').fill(momentum.z() - momentum_new.z())
+ self.aida.cloud1D('NewTrack/diffMomentum unsmeared x').fill(momentum.x() - momentum_newUnsmeared.x())
+ self.aida.cloud1D('NewTrack/diffMomentum unsmeared y').fill(momentum.y() - momentum_newUnsmeared.y())
+ self.aida.cloud1D('NewTrack/diffMomentum unsmeared z').fill(momentum.z() - momentum_newUnsmeared.z())
+ self.aida.cloud1D('NewTrack/smeared - unsmeared tracklength to point').fill(dist_new - dist_newUnsmeared)
#print emap_old
#print emap_new
lcsim/src/org/lcsim/contrib/JanStrube/tracking
diff -u -r1.7 -r1.8
--- Helix.java 28 Oct 2006 00:48:42 -0000 1.7
+++ Helix.java 7 Dec 2006 06:11:24 -0000 1.8
@@ -25,236 +25,260 @@
* by Paul Avery.
*
* @author tonyj
- * @version $Id: Helix.java,v 1.7 2006/10/28 00:48:42 jstrube Exp $
+ * @version $Id: Helix.java,v 1.8 2006/12/07 06:11:24 jstrube Exp $
*/
public class Helix implements Trajectory {
- Hep3Vector origin;
+ Hep3Vector origin;
- double xCenter;
+ double xCenter;
- double yCenter;
+ double yCenter;
- private double radius;
+ private double radius;
- private double sinLambda;
-
- private double cosLambda;
-
- private double sinPhi;
-
- private double cosPhi;
-
- private double phi;
-
- // parameterization in terms of 'momentum'
- // A helix is a mathematical object and doesn't have "momentum",
- // but some of the used algorithms are expressed in terms of it.
- // That's OK, it's a private variable.
- SpaceVector momentum;
-
- private double abs_r;
-
- private double rho;
-
- /**
- * Creates a new instance of Helix. Parameters according to <a
- * href="doc-files/L3_helix.pdf">the L3 conventions</a><br />
- * Please also have a look at <img src="doc-files/Helix-1.png"
- * alt="Helix-1.png"> <img src="doc-files/Helix-2.png" alt="Helix-2.png">
- *
- * @param origin
- * A point on the helix
- * @param radius
- * The <em>signed</em> radius of curvature of the helix. The
- * conventions is such that for <em>positive</em> radii, the
- * direction is <em>clockwise</em>.
- * @param phi
- * The polar angle of the helix <em>momentum</em> in x-y plane
- * w/rt positive x-axis at the origin
- * @param lambda
- * The dip angle w/rt positive part of the x-y plane
- */
- public Helix(SpacePoint org, double r, double Phi, double lambda) {
- // if (abs(lambda) > PI/2.)
- // throw new IllegalArgumentException("lambda = " + lambda + " is
- // outside of -PI/2<lambda<PI/2");
- origin = org;
- radius = r;
- phi = Phi;
-
- // Calculate some useful quantities
- cosPhi = cos(phi);
- sinPhi = sin(phi);
- cosLambda = cos(lambda);
- sinLambda = sin(lambda);
- xCenter = origin.x() + radius * sinPhi;
- yCenter = origin.y() - radius * cosPhi;
- setSpatialParameters();
- }
-
- /**
- * returns a point on the Helix at a distance alpha from the origin along
- * the trajectory. alpha == 2*PI*radius/cosLambda corresponds to one
- * rotation in the x-y plane
- */
- public SpacePoint getPointAtLength(double alpha) {
- double angle = phi + alpha * rho;
- double x = xCenter - radius * sin(angle);
- double y = yCenter + radius * cos(angle);
- double z = origin.z() + alpha * sinLambda;
- return new CartesianPoint(x, y, z);
- }
-
- public SpacePoint getCenterXY() {
- return new CartesianPoint(xCenter, yCenter, 0);
- }
-
- /**
- * Calculates the distance along the Helix until it hits a cylinder centered
- * along z
- *
- * @param r
- * the radius of the cylinder
- * @return the distance along the trajectory or Double.NAN if the cylinder
- * can never be hit
- */
- public double getDistanceToInfiniteCylinder(double r) {
- double phiToCenter = atan2(yCenter, xCenter);
- double radiusOfCenter = sqrt(xCenter * xCenter + yCenter * yCenter);
- // Negative radius doesn't make sense
- if (r < 0)
- throw new IllegalArgumentException("radius " + r + "<0");
- double darg = r * r / (2. * radius * radiusOfCenter) - radiusOfCenter
- / (2. * radius) - radius / (2. * radiusOfCenter);
- double deltaphi = phi-phiToCenter;
- if (deltaphi < - Math.PI)
- deltaphi += 2.* Math.PI;
- if (deltaphi > Math.PI)
- deltaphi -= 2.*Math.PI;
- double diff = asin(darg) + deltaphi;
- double result = (radius / cosLambda) * diff;
- while (result < 0)
- result += Math.abs(radius / cosLambda) * 2 * Math.PI;
- return result;
- }
-
- public double getDistanceToZPlane(double z) {
- return (z - origin.z()) / sinLambda;
- }
-
- /**
- * @param alpha
- * The distance along the trajectory in the x-y plane
- * @return The unit vector of the momentum
- */
- public SpaceVector getUnitTangentAtLength(double alpha) {
- double angle = phi + alpha * rho;
- return new CartesianVector(cos(angle), sin(angle), sinLambda
- / cosLambda);
- }
-
- public double getRadius() {
- return radius;
- }
-
- /**
- * Calculates the <em>signed</em> distance in mm between the Helix and an
- * arbitrary point in Space
- *
- * @param point
- * the point in space to calculate the distance to
- * @return the distance in mm between the point and the helix at the point
- * of closest approach
- */
- // FIXME The translation in Z should be replaced, because it can lead to an
- // infinite loop and the formula is only valid in 2D
- public double getSignedClosestDifferenceToPoint(SpacePoint point) {
- double tanLambda = sinLambda / cosLambda;
- Hep3Vector pCrossB = new CartesianVector(momentum.y(), -momentum.x(), 0);
- Hep3Vector xDiff = VecOp.sub(origin, point);
- double pMag = momentum.magnitude();
-
- // translate along Z because algorithm can handle only numbers in the
- // first quadrant
- double zPos = xDiff.z();
- // these are two mutually exclusive cases and two while loops
- // may not be the best way to express this
- while (zPos > abs(radius * tanLambda * Math.PI / 4)) {
- zPos -= abs(radius * tanLambda * Math.PI / 4);
- }
- while (zPos < -abs(radius * tanLambda * Math.PI / 4)) {
- zPos += abs(radius * tanLambda * Math.PI / 4);
- }
- xDiff = new CartesianVector(xDiff.x(), xDiff.y(), zPos);
-
- double numerator = (-2 * VecOp.dot(xDiff, pCrossB) + pMag * rho
- * (xDiff.magnitudeSquared() - xDiff.z() * xDiff.z()))
- / radius;
- double denominator = 1 + sqrt(1 - 2 * rho * pMag
- * VecOp.dot(xDiff, pCrossB) / radius / radius + pMag * pMag
- * rho * rho
- * (xDiff.magnitudeSquared() - xDiff.z() * xDiff.z()) / radius
- / radius);
- return numerator / denominator;
- }
-
- /**
- * Swims the helix along its trajectory to the point of closest approach to
- * the given SpacePoint. For more info on swimming see <a
- * href=doc-files/fitting/transport.pdf> Paul Avery's excellent text</a>
- *
- * @param point
- * Point in Space to swim to.
- * @return the length along the trajectory
- */
- public double getTrackLengthToPoint(SpacePoint point) {
- double tanLambda = sinLambda / cosLambda;
- double pMag = momentum.magnitude();
- Hep3Vector pCrossB = new CartesianVector(momentum.y(), -momentum.x(), 0);
-
- Hep3Vector xDiff = VecOp.sub(origin, point);
- int addedQuarterPeriods = 0;
- if (abs(tanLambda) > 1e-10) {
- double zPos = xDiff.z();
- while (abs(zPos) > abs(radius*tanLambda*Math.PI)) {
- zPos -= signum(zPos)*abs(radius*tanLambda*Math.PI);
- ++addedQuarterPeriods;
- }
- // Make sure the helix is in the right quadrant for the atan
- if (zPos > 0 && addedQuarterPeriods > 0)
- addedQuarterPeriods *= -1;
- if (addedQuarterPeriods % 2 != 0)
- addedQuarterPeriods += signum(addedQuarterPeriods);
- xDiff = new CartesianVector(xDiff.x(), xDiff.y(), zPos);
- }//int addedHalfPeriods = 0;
- //while (abs(xDiff.z()) > abs(tanLambda * radius * 0.5*Math.PI)) {
- // addedHalfPeriods++;
- // xDiff = VecOp.sub(getPointAtLength(addedHalfPeriods*0.5*Math.PI/rho), point);
- //}
- double pz = momentum.z();
- double factorA1 = pMag - pz * pz / pMag - (VecOp.dot(xDiff, pCrossB))
- * rho;
- double factorA2 = (VecOp.dot(xDiff, momentum) - xDiff.z() * pz) * rho;
- return addedQuarterPeriods*abs(radius/cosLambda*Math.PI) -Math.atan2(factorA2, factorA1) / rho;
- }
-
- /**
- * Sets the parametrization in terms of "momentum" and charge
- *
- */
- private void setSpatialParameters() {
- abs_r = abs(radius);
- momentum = new CartesianVector(abs_r * cosPhi, abs_r * sinPhi, abs_r
- * sinLambda / cosLambda);
- rho = -cosLambda / radius;
- }
-
- public String toString() {
- String helix = String.format("Helix:\n");
- helix += String.format("radius: %f\n", radius);
- helix += String.format("phi: %f\n", phi);
- helix += String.format("lambda: %f\n", acos(cosLambda));
- helix += String.format("rho: %f\n", rho);
- return helix;
- }
+ private double sinLambda;
+
+ private double cosLambda;
+
+ private double sinPhi;
+
+ private double cosPhi;
+
+ private double phi;
+
+ // parameterization in terms of 'momentum'
+ // A helix is a mathematical object and doesn't have "momentum",
+ // but some of the used algorithms are expressed in terms of it.
+ // That's OK, it's a private variable.
+ SpaceVector momentum;
+
+ private double abs_r;
+
+ private double rho;
+
+ /**
+ * Creates a new instance of Helix. Parameters according to <a
+ * href="doc-files/L3_helix.pdf">the L3 conventions</a><br />
+ * Please also have a look at <img src="doc-files/Helix-1.png"
+ * alt="Helix-1.png"> <img src="doc-files/Helix-2.png" alt="Helix-2.png">
+ *
+ * @param origin
+ * A point on the helix
+ * @param radius
+ * The <em>signed</em> radius of curvature of the helix. The
+ * conventions is such that for <em>positive</em> radii, the
+ * direction is <em>clockwise</em>.
+ * @param phi
+ * The polar angle of the helix <em>momentum</em> in x-y plane
+ * w/rt positive x-axis at the origin
+ * @param lambda
+ * The dip angle w/rt positive part of the x-y plane
+ */
+ public Helix(SpacePoint org, double r, double Phi, double lambda) {
+ // if (abs(lambda) > PI/2.)
+ // throw new IllegalArgumentException("lambda = " + lambda + " is
+ // outside of -PI/2<lambda<PI/2");
+ origin = org;
+ radius = r;
+ phi = Phi;
+
+ // Calculate some useful quantities
+ cosPhi = cos(phi);
+ sinPhi = sin(phi);
+ cosLambda = cos(lambda);
+ sinLambda = sin(lambda);
+ xCenter = origin.x() + radius * sinPhi;
+ yCenter = origin.y() - radius * cosPhi;
+ setSpatialParameters();
+ }
+
+ /**
+ * returns a point on the Helix at a distance alpha from the origin along
+ * the trajectory. alpha == 2*PI*radius/cosLambda corresponds to one
+ * rotation in the x-y plane
+ */
+ public SpacePoint getPointAtLength(double alpha) {
+ double angle = phi + alpha * rho;
+ double x = xCenter - radius * sin(angle);
+ double y = yCenter + radius * cos(angle);
+ double z = origin.z() + alpha * sinLambda;
+ return new CartesianPoint(x, y, z);
+ }
+
+ public SpacePoint getCenterXY() {
+ return new CartesianPoint(xCenter, yCenter, 0);
+ }
+
+ /**
+ * Calculates the distance along the Helix until it hits a cylinder centered
+ * along z
+ *
+ * @param r
+ * the radius of the cylinder
+ * @return the distance along the trajectory or Double.NAN if the cylinder
+ * can never be hit
+ */
+ public double getDistanceToInfiniteCylinder(double r) {
+ double phiToCenter = atan2(yCenter, xCenter);
+ double radiusOfCenter = sqrt(xCenter * xCenter + yCenter * yCenter);
+ // Negative radius doesn't make sense
+ if (r < 0)
+ throw new IllegalArgumentException("radius " + r + "<0");
+ double darg = r * r / (2. * radius * radiusOfCenter) - radiusOfCenter
+ / (2. * radius) - radius / (2. * radiusOfCenter);
+ double deltaphi = phi - phiToCenter;
+ if (deltaphi < -Math.PI)
+ deltaphi += 2. * Math.PI;
+ if (deltaphi > Math.PI)
+ deltaphi -= 2. * Math.PI;
+ double diff = asin(darg) + deltaphi;
+ double result = (radius / cosLambda) * diff;
+ while (result < 0)
+ result += Math.abs(radius / cosLambda) * 2 * Math.PI;
+ return result;
+ }
+
+ public double getDistanceToZPlane(double z) {
+ return (z - origin.z()) / sinLambda;
+ }
+
+ /**
+ * @param alpha
+ * The distance along the trajectory in the x-y plane
+ * @return The unit vector of the momentum
+ */
+ public SpaceVector getUnitTangentAtLength(double alpha) {
+ double angle = phi + alpha * rho;
+ return new CartesianVector(cos(angle), sin(angle), sinLambda
+ / cosLambda);
+ }
+
+ public double getRadius() {
+ return radius;
+ }
+
+ /**
+ * Calculates the <em>signed</em> distance in mm between the Helix and an
+ * arbitrary point in Space
+ *
+ * @param point
+ * the point in space to calculate the distance to
+ * @return the distance in mm between the point and the helix at the point
+ * of closest approach
+ */
+ // FIXME The translation in Z should be replaced, because it can lead to an
+ // infinite loop and the formula is only valid in 2D
+ public double getSignedClosestDifferenceToPoint(SpacePoint point) {
+ double tanLambda = sinLambda / cosLambda;
+ Hep3Vector pCrossB = new CartesianVector(momentum.y(), -momentum.x(), 0);
+ Hep3Vector xDiff = VecOp.sub(origin, point);
+ double pMag = momentum.magnitude();
+
+ // translate along Z because algorithm can handle only numbers in the
+ // first quadrant
+ double zPos = xDiff.z();
+ // these are two mutually exclusive cases and two while loops
+ // may not be the best way to express this
+ while (zPos > abs(radius * tanLambda * Math.PI / 4)) {
+ zPos -= abs(radius * tanLambda * Math.PI / 4);
+ }
+ while (zPos < -abs(radius * tanLambda * Math.PI / 4)) {
+ zPos += abs(radius * tanLambda * Math.PI / 4);
+ }
+ xDiff = new CartesianVector(xDiff.x(), xDiff.y(), zPos);
+
+ double numerator = (-2 * VecOp.dot(xDiff, pCrossB) + pMag * rho
+ * (xDiff.magnitudeSquared() - xDiff.z() * xDiff.z()))
+ / radius;
+ double denominator = 1 + sqrt(1 - 2 * rho * pMag
+ * VecOp.dot(xDiff, pCrossB) / radius / radius + pMag * pMag
+ * rho * rho
+ * (xDiff.magnitudeSquared() - xDiff.z() * xDiff.z()) / radius
+ / radius);
+ return numerator / denominator;
+ }
+
+ // added by Nick Sinev
+ public double getSecondDistanceToInfiniteCylinder(double r) {
+ double result = getDistanceToInfiniteCylinder(r);
+ SpacePoint first = getPointAtLength(result);
+ double angto0 = Math.atan2(-yCenter, -xCenter);
+ double angtofirst = Math
+ .atan2(first.y() - yCenter, first.x() - xCenter);
+ double dang = angtofirst - angto0;
+ if (dang < -Math.PI)
+ dang += 2. * Math.PI;
+ if (dang > Math.PI)
+ dang -= 2. * Math.PI;
+ double angofarc = 2. * (Math.PI - Math.abs(dang));
+ double arc = Math.abs(radius) * angofarc / cosLambda;
+ return result + arc;
+ }
+
+
+ public double getZPeriod() {
+ return Math.PI * radius * 2. * sinLambda / cosLambda;
+ }
+
+ /**
+ * Swims the helix along its trajectory to the point of closest approach to
+ * the given SpacePoint. For more info on swimming see <a
+ * href=doc-files/fitting/transport.pdf> Paul Avery's excellent text</a>
+ *
+ * @param point
+ * Point in Space to swim to.
+ * @return the length along the trajectory
+ */
+ public double getTrackLengthToPoint(SpacePoint point) {
+ double tanLambda = sinLambda / cosLambda;
+ double pMag = momentum.magnitude();
+ Hep3Vector pCrossB = new CartesianVector(momentum.y(), -momentum.x(), 0);
+
+ Hep3Vector xDiff = VecOp.sub(origin, point);
+ int addedQuarterPeriods = 0;
+ if (abs(tanLambda) > 1e-10) {
+ double zPos = xDiff.z();
+ while (abs(zPos) > abs(radius * tanLambda * Math.PI)) {
+ zPos -= signum(zPos) * abs(radius * tanLambda * Math.PI);
+ ++addedQuarterPeriods;
+ }
+ // Make sure the helix is in the right quadrant for the atan
+ if (zPos > 0 && addedQuarterPeriods > 0)
+ addedQuarterPeriods *= -1;
+ if (addedQuarterPeriods % 2 != 0)
+ addedQuarterPeriods += signum(addedQuarterPeriods);
+ xDiff = new CartesianVector(xDiff.x(), xDiff.y(), zPos);
+ }// int addedHalfPeriods = 0;
+ // while (abs(xDiff.z()) > abs(tanLambda * radius * 0.5*Math.PI)) {
+ // addedHalfPeriods++;
+ // xDiff = VecOp.sub(getPointAtLength(addedHalfPeriods*0.5*Math.PI/rho),
+ // point);
+ // }
+ double pz = momentum.z();
+ double factorA1 = pMag - pz * pz / pMag - (VecOp.dot(xDiff, pCrossB))
+ * rho;
+ double factorA2 = (VecOp.dot(xDiff, momentum) - xDiff.z() * pz) * rho;
+ return addedQuarterPeriods * abs(radius / cosLambda * Math.PI)
+ - Math.atan2(factorA2, factorA1) / rho;
+ }
+
+ /**
+ * Sets the parametrization in terms of "momentum" and charge
+ *
+ */
+ private void setSpatialParameters() {
+ abs_r = abs(radius);
+ momentum = new CartesianVector(abs_r * cosPhi, abs_r * sinPhi, abs_r
+ * sinLambda / cosLambda);
+ rho = -cosLambda / radius;
+ }
+
+ public String toString() {
+ String helix = String.format("Helix:\n");
+ helix += String.format("radius: %f\n", radius);
+ helix += String.format("phi: %f\n", phi);
+ helix += String.format("lambda: %f\n", acos(cosLambda));
+ helix += String.format("rho: %f\n", rho);
+ return helix;
+ }
}
CVSspam 0.2.8
