hps-java/src/main/java/org/lcsim/hps/users/mgraham
diff -N ElwinsTrackingRecon.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ElwinsTrackingRecon.java 19 Nov 2012 17:12:29 -0000 1.1
@@ -0,0 +1,1558 @@
+package org.lcsim.hps.users.mgraham;
+
+import Jama.*;
+import hep.aida.*;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.Hep3Vector;
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.List;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.event.*;
+import org.lcsim.fit.helicaltrack.HelicalTrackCross;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.hps.monitoring.AIDAFrame;
+import org.lcsim.hps.monitoring.Resettable;
+import org.lcsim.hps.recon.ecal.HPSEcalCluster;
+import org.lcsim.hps.recon.tracking.*;
+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;
+
+/**
+ *
+ * @author elwinm
+ */
+public class ElwinsTrackingRecon extends Driver implements Resettable {
+
+ private AIDAFrame plotterFrame;
+ private AIDAFrame topFrame;
+ private AIDAFrame bottomFrame;
+ private AIDAFrame chargeFrame;
+ private AIDAFrame twotrackFrame;
+ private AIDA aida = AIDA.defaultInstance();
+ private String rawTrackerHitCollectionName = "SVTRawTrackerHits";
+ private String fittedTrackerHitCollectionName = "SVTFittedRawTrackerHits";
+ private String trackerHitCollectionName = "StripClusterer_SiTrackerHitStrip1D";
+ private String helicalTrackHitCollectionName = "HelicalTrackHits";
+ private String rotatedTrackHitCollectionName = "RotatedHelicalTrackHits";
+ private String helicalTrackHitRelationsCollectionName = "HelicalTrackHitRelations";
+ private String trackCollectionName = "MatchedTracks";
+ private String trackerName = "Tracker";
+ String ecalSubdetectorName = "Ecal";
+ String ecalCollectionName = "EcalClusters";
+ private Detector detector = null;
+ IDDecoder dec;
+ private int eventCount;
+ private List<SiSensor> sensors;
+ private double zAtConverter = -674.062;//mm
+ private String outputPlots = null;
+ IPlotter plotter;
+ IPlotter plotter2;
+ IPlotter plotter22;
+ IPlotter plotter222;
+ IPlotter plotter3;
+ IPlotter plotter3_1;
+ IPlotter plotter3_2;
+ IPlotter plotter4;
+ IPlotter plotter5;
+ IPlotter plotter5_1;
+ IPlotter plotter55;
+ IPlotter plotter6;
+ IPlotter plotter7;
+ IPlotter plotter9000;
+ IPlotter plotter9001;
+ IPlotter plotter9002;
+ IPlotter plotter9003;
+ IPlotter plotter9004;
+ IPlotter plotter9005;
+ IPlotter plotter9006;
+ IPlotter plotter9007;
+ IPlotter plotter9008;
+ IPlotter plotter9009;
+ IPlotter plotter9010;
+ IPlotter plotter9011;
+ IPlotter plotter9012;
+ IPlotter plotter9013;
+ IPlotter plotter9014;
+ IPlotter plotter9015;
+ IPlotter plotter9016;
+ IPlotter plotter9017;
+ IPlotter twotrkextra;
+ IPlotter twotrkextra2;
+ IPlotter threetrack;
+ IPlotter top1;
+ IPlotter top2;
+ IPlotter top3;
+ IPlotter top4;
+ IPlotter bot1;
+ IPlotter bot2;
+ IPlotter bot3;
+ IPlotter charge;
+ IPlotter bot4;
+ double zEcal = 1500;
+ double zAtDownStrPairSpec = 914.0; //mm
+ double zAtColl = -1500;
+ IHistogram1D trkPx;
+ IHistogram1D nTracks;
+ HPSShaperFitAlgorithm _shaper = new DumbShaperFit();
+
+ protected void detectorChanged(Detector detector) {
+ this.detector = detector;
+ aida.tree().cd("/");
+ plotterFrame = new AIDAFrame();
+ plotterFrame.setTitle("HPS Tracking Plots");
+
+ twotrackFrame = new AIDAFrame();
+ twotrackFrame.setTitle("Two Track Plots");
+
+ sensors = detector.getSubdetector(trackerName).getDetectorElement().findDescendants(SiSensor.class);
+
+ IAnalysisFactory fac = aida.analysisFactory();
+ plotter = fac.createPlotterFactory().create("HPS Tracking Plots");
+ plotter.setTitle("Momentum");
+ IPlotterStyle style = plotter.style();
+ style.dataStyle().fillStyle().setColor("yellow");
+ style.dataStyle().errorBarStyle().setVisible(false);
+ plotter.createRegions(2, 2);
+ plotterFrame.addPlotter(plotter);
+
+ trkPx = aida.histogram1D("Track X Momentum", 25, -0.25, 0.25);
+ IHistogram1D trkPy = aida.histogram1D("Track Y Momentum", 25, -0.1, 0.1);
+ IHistogram1D trkPz = aida.histogram1D("Track Z Momentum", 25, 0, 3.5);
+ IHistogram1D trkChi2 = aida.histogram1D("Track Chi2", 25, 0, 25.0);
+
+ plotter.region(0).plot(trkPx);
+ plotter.region(1).plot(trkPy);
+ plotter.region(2).plot(trkPz);
+ plotter.region(3).plot(trkChi2);
+
+
+ plotter2 = fac.createPlotterFactory().create("HPS Tracking Plots");
+ plotter2.setTitle("Track extrapolation");
+ plotterFrame.addPlotter(plotter2);
+ IPlotterStyle style2 = plotter2.style();
+ style2.dataStyle().fillStyle().setColor("yellow");
+ style2.dataStyle().errorBarStyle().setVisible(false);
+ plotter2.createRegions(2, 4);
+ IHistogram1D xAtConverter = aida.histogram1D("X (mm) @ Z=-60cm", 50, -50, 50);
+ IHistogram1D yAtConverter = aida.histogram1D("Y (mm) @ Z=-60cm", 50, -20, 20);
+ IHistogram1D xAtColl = aida.histogram1D("X (mm) @ Z=-150cm", 50, -200, 200);
+ IHistogram1D yAtColl = aida.histogram1D("Y (mm) @ Z=-150cm", 50, -200, 200);
+ IHistogram1D xAtEcal = aida.histogram1D("X (mm) @ ECAL", 50, -500, 500);
+ IHistogram1D yAtEcal = aida.histogram1D("Y (mm) @ ECAL", 50, -100, 100);
+ IHistogram1D xAtConvert = aida.histogram1D("X (mm) @ Converter", 50, -50, 50);
+ IHistogram1D yAtConvert = aida.histogram1D("Y (mm) @ Converter", 50, -20, 20);
+
+ plotter2.region(0).plot(xAtConverter);
+ plotter2.region(4).plot(yAtConverter);
+ plotter2.region(1).plot(xAtColl);
+ plotter2.region(5).plot(yAtColl);
+ plotter2.region(2).plot(xAtEcal);
+ plotter2.region(6).plot(yAtEcal);
+ plotter2.region(3).plot(xAtConvert);
+ plotter2.region(7).plot(yAtConvert);
+
+ twotrkextra = fac.createPlotterFactory().create("Two Trk Extrapolation");
+ twotrkextra.setTitle("Stuff");
+ plotterFrame.addPlotter(twotrkextra);
+ IPlotterStyle styletwo = twotrkextra.style();
+ styletwo.dataStyle().fillStyle().setColor("blue");
+ styletwo.dataStyle().errorBarStyle().setVisible(false);
+ twotrkextra.createRegions(3, 2);
+ IHistogram1D x1AtTarget = aida.histogram1D("Trk1 X @ Target", 50, 0, 50);
+ IHistogram1D y1AtTarget = aida.histogram1D("Trk1 Y @ Target", 50, -5, 5);
+ IHistogram1D x2AtTarget = aida.histogram1D("Trk2 X @ Target", 50, 0, 50);
+ IHistogram1D y2AtTarget = aida.histogram1D("Trk2 Y @ Target", 50, -5, 5);
+ IHistogram1D distatt = aida.histogram1D("Distance btwn Trks @ Target", 40, 0, 40);
+ IHistogram1D zdiff = aida.histogram1D("Z Diff", 40, -.1, .1);
+
+ twotrkextra.region(0).plot(x1AtTarget);
+ twotrkextra.region(1).plot(y1AtTarget);
+ twotrkextra.region(2).plot(x2AtTarget);
+ twotrkextra.region(3).plot(y2AtTarget);
+ twotrkextra.region(4).plot(distatt);
+ twotrkextra.region(5).plot(zdiff);
+
+
+ plotter222 = fac.createPlotterFactory().create("HPS Tracking Plots");
+ plotter222.setTitle("Other");
+ plotterFrame.addPlotter(plotter222);
+ IPlotterStyle style222 = plotter222.style();
+ style222.dataStyle().fillStyle().setColor("yellow");
+ style222.dataStyle().errorBarStyle().setVisible(false);
+ plotter222.createRegions(2, 3);
+
+ IHistogram1D nHits = aida.histogram1D("Hits per Track", 2, 4, 6);
+ IHistogram1D amp = aida.histogram1D("Amp (HitOnTrack)", 50, 0, 5000);
+ IHistogram1D ampcl = aida.histogram1D("Amp (CluOnTrack)", 50, 0, 5000);
+ IHistogram1D amp2 = aida.histogram1D("Amp Pz>1000 (HitOnTrack)", 50, 0, 5000);
+ IHistogram1D ampcl2 = aida.histogram1D("Amp Pz>1000 (CluOnTrack)", 50, 0, 5000);
+ nTracks = aida.histogram1D("Tracks per Event", 3, 0, 3);
+
+ plotter222.region(0).plot(nHits);
+ plotter222.region(3).plot(nTracks);
+ plotter222.region(1).plot(amp);
+ plotter222.region(4).plot(amp2);
+ plotter222.region(2).plot(ampcl);
+ plotter222.region(5).plot(ampcl2);
+
+
+ plotterFrame.pack();
+ plotterFrame.setVisible(true);
+
+
+ twotrkextra2 = fac.createPlotterFactory().create("Two Trk Uncertainties");
+ twotrkextra2.setTitle("Uncertainties");
+ plotter9000 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9000.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks = plotter9000.style();
+ plotter9001 = fac.createPlotterFactory().create("Two Track Plots 2");
+ plotter9001.setTitle("Two Track Plots Test 2");
+ IPlotterStyle TwoTracks1 = plotter9001.style();
+ plotter9002 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9002.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks2 = plotter9002.style();
+ plotter9003 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9003.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks3 = plotter9003.style();
+ plotter9004 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9004.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks4 = plotter9004.style();
+ plotter9005 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9005.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks5 = plotter9005.style();
+ plotter9006 = fac.createPlotterFactory().create("Two Track Versus");
+ plotter9006.setTitle("Two Track Versus");
+ IPlotterStyle TwoTracks6 = plotter9006.style();
+ plotter9007 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9007.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks7 = plotter9007.style();
+ plotter9008 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9008.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks8 = plotter9000.style();
+ plotter9009 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9009.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks9 = plotter9000.style();
+ plotter9010 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9010.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks10 = plotter9010.style();
+ plotter9011 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9011.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks11 = plotter9011.style();
+ plotter9012 = fac.createPlotterFactory().create("Two Track Plotz");
+ plotter9012.setTitle("Two Track Plotz Test");
+ IPlotterStyle TwoTracks12 = plotter9012.style();
+ plotter9013 = fac.createPlotterFactory().create("Two Track Plotz");
+ plotter9013.setTitle("Two Track Plotz Test");
+ IPlotterStyle TwoTracks13 = plotter9013.style();
+ plotter9014 = fac.createPlotterFactory().create("Two Track Plotz");
+ plotter9014.setTitle("Two Track Plotz Test");
+ IPlotterStyle TwoTracks14 = plotter9014.style();
+ plotter9015 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9015.setTitle("Two Track Plots Test");
+ IPlotterStyle TwoTracks15 = plotter9015.style();
+ plotter9016 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9016.setTitle("Test");
+ IPlotterStyle TwoTracks16 = plotter9016.style();
+ plotter9017 = fac.createPlotterFactory().create("Two Track Plots");
+ plotter9017.setTitle("Residuals");
+ IPlotterStyle TwoTracks17 = plotter9017.style();
+ threetrack = fac.createPlotterFactory().create("Three Track Plots");
+ threetrack.setTitle("Invariant Mass");
+
+
+ TwoTracks.dataStyle().fillStyle().setColor("green");
+ TwoTracks.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks.setParameter("hist2DStyle", "colorMap");
+ TwoTracks.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks1.dataStyle().fillStyle().setColor("green");
+ TwoTracks1.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks1.setParameter("hist2DStyle", "colorMap");
+ TwoTracks1.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks2.dataStyle().fillStyle().setColor("green");
+ TwoTracks2.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks2.setParameter("hist2DStyle", "colorMap");
+ TwoTracks2.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks3.dataStyle().fillStyle().setColor("green");
+ TwoTracks3.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks3.setParameter("hist2DStyle", "colorMap");
+ TwoTracks3.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks4.dataStyle().fillStyle().setColor("green");
+ TwoTracks4.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks4.setParameter("hist2DStyle", "colorMap");
+ TwoTracks4.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks5.dataStyle().fillStyle().setColor("green");
+ TwoTracks5.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks5.setParameter("hist2DStyle", "colorMap");
+ TwoTracks5.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks6.dataStyle().fillStyle().setColor("green");
+ TwoTracks6.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks6.setParameter("hist2DStyle", "colorMap");
+ TwoTracks6.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks7.dataStyle().fillStyle().setColor("green");
+ TwoTracks7.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks7.setParameter("hist2DStyle", "colorMap");
+ TwoTracks7.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks8.dataStyle().fillStyle().setColor("green");
+ TwoTracks8.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks8.setParameter("hist2DStyle", "colorMap");
+ TwoTracks8.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks9.dataStyle().fillStyle().setColor("green");
+ TwoTracks9.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks9.setParameter("hist2DStyle", "colorMap");
+ TwoTracks9.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks10.dataStyle().fillStyle().setColor("green");
+ TwoTracks10.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks10.setParameter("hist2DStyle", "colorMap");
+ TwoTracks10.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks11.dataStyle().fillStyle().setColor("green");
+ TwoTracks11.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks11.setParameter("hist2DStyle", "colorMap");
+ TwoTracks11.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks12.dataStyle().fillStyle().setColor("green");
+ TwoTracks12.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks13.dataStyle().fillStyle().setColor("green");
+ TwoTracks13.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks14.dataStyle().fillStyle().setColor("green");
+ TwoTracks14.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks15.dataStyle().fillStyle().setColor("green");
+ TwoTracks15.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks15.setParameter("hist2DStyle", "colorMap");
+ TwoTracks15.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ TwoTracks16.dataStyle().fillStyle().setColor("green");
+ TwoTracks16.dataStyle().errorBarStyle().setVisible(false);
+ TwoTracks16.setParameter("hist2DStyle", "colorMap");
+ TwoTracks16.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ IPlotterStyle styletwotwo = twotrkextra2.style();
+ styletwotwo.dataStyle().fillStyle().setColor("blue");
+ styletwotwo.dataStyle().errorBarStyle().setVisible(false);
+
+ IPlotterStyle threesty = threetrack.style();
+ threesty.dataStyle().fillStyle().setColor("blue");
+ threesty.dataStyle().errorBarStyle().setVisible(false);
+
+
+
+ twotrkextra2.createRegions(3, 2);
+ plotter9000.createRegion();
+ plotter9001.createRegion();
+ plotter9002.createRegions(2, 1);
+ plotter9003.createRegion();
+ plotter9004.createRegion();
+ plotter9005.createRegion();
+ plotter9006.createRegions(1, 3);
+ plotter9007.createRegions(2, 3);
+ plotter9008.createRegions(2, 2);
+ plotter9009.createRegions(2, 2);
+ plotter9010.createRegions(2, 2);
+ plotter9011.createRegions(2, 2);
+ plotter9012.createRegions(2, 2);
+ // plotter9013.createRegions(2, 2);
+ // plotter9014.createRegions(2, 2);
+ plotter9015.createRegions(2, 2);
+ plotter9016.createRegions(2, 2);
+ plotter9017.createRegions(2, 2);
+ threetrack.createRegion();
+
+ twotrackFrame.addPlotter(plotter9000);
+ twotrackFrame.addPlotter(plotter9001);
+ twotrackFrame.addPlotter(plotter9002);
+ twotrackFrame.addPlotter(plotter9003);
+ twotrackFrame.addPlotter(plotter9004);
+ twotrackFrame.addPlotter(plotter9005);
+ twotrackFrame.addPlotter(plotter9006);
+ twotrackFrame.addPlotter(plotter9007);
+ twotrackFrame.addPlotter(plotter9008);
+ twotrackFrame.addPlotter(plotter9009);
+ twotrackFrame.addPlotter(plotter9010);
+ twotrackFrame.addPlotter(plotter9011);
+ twotrackFrame.addPlotter(plotter9012);
+ // twotrackFrame.addPlotter(plotter9013);
+ // twotrackFrame.addPlotter(plotter9014);
+ twotrackFrame.addPlotter(plotter9015);
+ twotrackFrame.addPlotter(plotter9016);
+ twotrackFrame.addPlotter(plotter9017);
+ twotrackFrame.addPlotter(twotrkextra2);
+ twotrackFrame.addPlotter(threetrack);
+
+ IHistogram1D trkbins = aida.histogram1D("Track Distributions", 5, -2, 3);
+ IHistogram2D twtrkptot = aida.histogram2D("Total P+ vs. P-", 60, 0, 4, 60, 0, 4);
+ IHistogram1D sumtrks = aida.histogram1D("Sum of Track's Momentums", 100, -1, 7);
+ IHistogram1D invarmass = aida.histogram1D("Invariant Mass", 50, 0, .2);
+ IHistogram1D perptrks = aida.histogram1D("Perpendicular Momentum", 100, 0, .1);
+ IHistogram2D pyppm = aida.histogram2D("Py+ vs. Py-", 60, -.1, .1, 60, -.1, .1);
+ IHistogram2D pzppm = aida.histogram2D("Pz+ vs. Pz-", 60, -.1, .1, 60, -.1, .1);
+ IHistogram1D px = aida.histogram1D("Two Track X Momentum", 40, 0, 4);
+ IHistogram1D py = aida.histogram1D("Two Track Y Momentum", 40, -.1, .1);
+ IHistogram1D pz = aida.histogram1D("Two Track Z Momentum", 40, -.1, .1);
+ IHistogram1D chi2 = aida.histogram1D("Tracks Chi2", 25, 0, 25.0);
+ IHistogram1D bbpx = aida.histogram1D("Big Bump Track Momenta (Px)", 40, 0, 4);
+ IHistogram1D bbpy = aida.histogram1D("Big Bump Track Momenta (Py)", 40, -.1, .1);
+ IHistogram1D bbpz = aida.histogram1D("Big Bump Track Momenta (Pz)", 40, -.1, .1);
+ IHistogram1D bbchi2 = aida.histogram1D("Big Bump Tracks Chi2", 25, 0, 25.0);
+ IHistogram1D spx = aida.histogram1D("Split Track Momenta (Px)", 40, 0, 4);
+ IHistogram1D spy = aida.histogram1D("Split Track Momenta (Py)", 40, -.1, .1);
+ IHistogram1D spz = aida.histogram1D("Split Track Momenta (Pz)", 40, -.1, .1);
+ IHistogram1D schi2 = aida.histogram1D("Split Tracks Chi2", 25, 0, 25.0);
+ IHistogram1D bbsumtrks = aida.histogram1D("Big Bump Sum of Track's Momentums", 50, -1, 7);
+ IHistogram2D bbpppm = aida.histogram2D("Big Bump P+ vs. P-", 50, 0, 4, 50, 0, 4);
+ IHistogram2D lbpppm = aida.histogram2D("Little Bump P+ vs. P-", 50, 0, 4, 50, 0, 4);
+ IHistogram1D lbsumtrks = aida.histogram1D("Little Bump Sum of Track's Momentums", 50, -1, 7);
+ IHistogram1D lbpx = aida.histogram1D("Little Bump Track Momenta (Px)", 40, 0, 4);
+ IHistogram1D lbpy = aida.histogram1D("Little Bump Track Momenta (Py)", 40, -.1, .1);
+ IHistogram1D lbpz = aida.histogram1D("Little Bump Track Momenta (Pz)", 40, -.1, .1);
+ IHistogram1D lbchi2 = aida.histogram1D("Little Bump Tracks Chi2", 25, 0, 25.0);
+ // IHistogram1D q0spx = aida.histogram1D("Net Charge 0 Split Track Momenta (Px)", 40, 0, 4);
+ // IHistogram1D q0spy = aida.histogram1D("Net Charge 0 Split Track Momenta (Py)", 40, -.1, .1);
+ // IHistogram1D q0spz = aida.histogram1D("Net Charge 0 Split Track Momenta (Pz)", 40, -.1, .1);
+ // IHistogram1D q0schi2 = aida.histogram1D("Net Charge 0 Split Tracks Chi2", 25, 0, 25.0);
+ IHistogram2D xyemt = aida.histogram2D("X v Y - e- Top", 50, -30, 50, 50, -35, 30);
+ IHistogram2D xzemt = aida.histogram2D("X v Z - e- Top", 50, -30, 50, 50, -800, -450);
+ IHistogram2D yzemt = aida.histogram2D("Y v Z - e- Top", 50, -35, 30, 50, -800, -450);
+ IHistogram1D qbins = aida.histogram1D("Charge Distributions", 5, -2, 3);
+ IHistogram1D lbtp = aida.histogram1D("Little Bump Track Parity", 7, 0, 7);
+ IHistogram1D bbtp = aida.histogram1D("Big Bump Track Parity", 7, 0, 7);
+ IHistogram1D xvert = aida.histogram1D("XVertex", 40, -30, 50);
+ IHistogram1D yvert = aida.histogram1D("YVertex", 40, -35, 30);
+ IHistogram1D zvert = aida.histogram1D("ZVertex", 40, -800, -450);
+ IHistogram1D dist = aida.histogram1D("Distance btwn Trks @ Solution", 40, 0, 20);
+ IHistogram1D xres = aida.histogram1D("X Res Trk1", 40, -0.25, 0.25);
+ IHistogram1D yres = aida.histogram1D("Y Res Trk1", 40, -0.25, 0.25);
+ IHistogram1D xres2 = aida.histogram1D("X Res Trk2", 40, -0.25, 0.25);
+ IHistogram1D yres2 = aida.histogram1D("Y Res Trk2", 40, -0.25, 0.25);
+ IHistogram1D unx1 = aida.histogram1D("Uncert X Trk 1", 50, 0, 10);
+ IHistogram1D uny1 = aida.histogram1D("Uncert Y Trk 1", 50, 0, 10);
+ IHistogram1D unz1 = aida.histogram1D("Uncert Z Trk 1", 50, 0, 40);
+ IHistogram1D unx2 = aida.histogram1D("Uncert X Trk 2", 50, 0, 10);
+ IHistogram1D uny2 = aida.histogram1D("Uncert Y Trk 2", 50, 0, 10);
+ IHistogram1D unz2 = aida.histogram1D("Uncert Z Trk 2", 50, 0, 40);
+ IHistogram2D xy = aida.histogram2D("X v Y", 50, -30, 50, 50, -35, 30);
+ IHistogram2D xz = aida.histogram2D("X v Z", 50, -30, 50, 50, -800, -450);
+ IHistogram2D yz = aida.histogram2D("Y v Z", 50, -35, 30, 50, -800, -450);
+ IHistogram2D xyept = aida.histogram2D("X v Y - e+ Top", 50, -30, 50, 50, -35, 30);
+ IHistogram2D xzept = aida.histogram2D("X v Z - e+ Top", 50, -30, 50, 50, -800, -450);
+ IHistogram2D yzept = aida.histogram2D("Y v Z - e+ Top", 50, -35, 30, 50, -800, -450);
+ IHistogram1D three = aida.histogram1D("Three Track Invariant Mass", 50, 0, .4);
+
+ twotrackFrame.pack();
+ twotrackFrame.setVisible(true);
+
+ plotter9000.region(0).plot(trkbins);
+ plotter9001.region(0).plot(twtrkptot);
+ plotter9002.region(0).plot(sumtrks);
+ plotter9002.region(1).plot(invarmass);
+ plotter9003.region(0).plot(perptrks);
+ plotter9004.region(0).plot(pyppm);
+ plotter9005.region(0).plot(pzppm);
+ plotter9006.region(0).plot(xy);
+ plotter9006.region(1).plot(xz);
+ plotter9006.region(2).plot(yz);
+ plotter9007.region(0).plot(xyemt);
+ plotter9007.region(1).plot(xzemt);
+ plotter9007.region(2).plot(yzemt);
+ plotter9007.region(3).plot(xyept);
+ plotter9007.region(4).plot(xzept);
+ plotter9007.region(5).plot(yzept);
+ plotter9008.region(0).plot(px);
+ plotter9008.region(1).plot(py);
+ plotter9008.region(2).plot(pz);
+ plotter9008.region(3).plot(chi2);
+ plotter9009.region(0).plot(bbpx);
+ plotter9009.region(1).plot(bbpy);
+ plotter9009.region(2).plot(bbpz);
+ plotter9009.region(3).plot(bbchi2);
+ plotter9010.region(0).plot(spx);
+ plotter9010.region(1).plot(spy);
+ plotter9010.region(2).plot(spz);
+ plotter9010.region(3).plot(schi2);
+ plotter9011.region(0).plot(bbsumtrks);
+ plotter9011.region(1).plot(bbpppm);
+ plotter9011.region(2).plot(lbpppm);
+ plotter9011.region(3).plot(lbsumtrks);
+ plotter9012.region(0).plot(lbpx);
+ plotter9012.region(1).plot(lbpy);
+ plotter9012.region(2).plot(lbpz);
+ plotter9012.region(3).plot(lbchi2);
+ // plotter9013.region(0).plot(q0spx);
+ // plotter9013.region(1).plot(q0spy);
+ // plotter9013.region(2).plot(q0spz);
+ // plotter9013.region(3).plot(q0schi2);
+ plotter9015.region(0).plot(qbins);
+ plotter9015.region(1).plot(lbtp);
+ plotter9015.region(2).plot(bbtp);
+ plotter9016.region(0).plot(xvert);
+ plotter9016.region(1).plot(yvert);
+ plotter9016.region(2).plot(zvert);
+ plotter9016.region(3).plot(dist);
+ plotter9017.region(0).plot(xres);
+ plotter9017.region(1).plot(yres);
+ plotter9017.region(2).plot(xres2);
+ plotter9017.region(3).plot(yres2);
+
+ twotrkextra2.region(0).plot(unx1);
+ twotrkextra2.region(1).plot(uny1);
+ twotrkextra2.region(2).plot(unz1);
+ twotrkextra2.region(3).plot(unx2);
+ twotrkextra2.region(4).plot(uny2);
+ twotrkextra2.region(5).plot(unz2);
+ threetrack.region(0).plot(three);
+
+
+
+ }
+
+ public ElwinsTrackingRecon() {
+ }
+
+ public void setOutputPlots(String output) {
+ this.outputPlots = output;
+ }
+
+ public void setRawTrackerHitCollectionName(String rawTrackerHitCollectionName) {
+ this.rawTrackerHitCollectionName = rawTrackerHitCollectionName;
+ }
+
+ public void setFittedTrackerHitCollectionName(String fittedTrackerHitCollectionName) {
+ this.fittedTrackerHitCollectionName = fittedTrackerHitCollectionName;
+ }
+
+ public void setTrackerHitCollectionName(String trackerHitCollectionName) {
+ this.trackerHitCollectionName = trackerHitCollectionName;
+ }
+
+ public void setHelicalTrackHitCollectionName(String helicalTrackHitCollectionName) {
+ this.helicalTrackHitCollectionName = helicalTrackHitCollectionName;
+ }
+
+ public void setTrackCollectionName(String trackCollectionName) {
+ this.trackCollectionName = trackCollectionName;
+ }
+
+ public void process(EventHeader event) {
+ aida.tree().cd("/");
+ if (!event.hasCollection(HelicalTrackHit.class, helicalTrackHitCollectionName)) {
+ // System.out.println(helicalTrackHitCollectionName + " does not exist; skipping event");
+ return;
+ }
+ if (event.get(Track.class, trackCollectionName).size() < 2) {
+ // System.out.println(trackCollectionName + " has less than two tracks; skipping event");
+ return;
+ }
+
+ List<HelicalTrackHit> rotList = event.get(HelicalTrackHit.class, rotatedTrackHitCollectionName);
+ for (HelicalTrackHit hth : rotList) {
+ HelicalTrackCross htc = (HelicalTrackCross) hth;
+// System.out.println("TrackingReconstructionPlots::original helical track position = "+hth.getPosition()[0]+","+hth.getPosition()[1]+","+hth.getPosition()[2]);
+// System.out.println("TrackingReconstructionPlots::corrected helical track position = "+htc.getCorrectedPosition().toString());
+ }
+
+ List<HelicalTrackHit> hthList = event.get(HelicalTrackHit.class, helicalTrackHitCollectionName);
+ int[] layersTop = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ int[] layersBot = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ for (HelicalTrackHit hth : hthList) {
+ HelicalTrackCross htc = (HelicalTrackCross) hth;
+// System.out.println("TrackingReconstructionPlots::original helical track position = "+hth.getPosition()[0]+","+hth.getPosition()[1]+","+hth.getPosition()[2]);
+// System.out.println("TrackingReconstructionPlots::corrected helical track position = "+htc.getCorrectedPosition().toString());
+ //These Helical Track Hits are in the JLAB frame
+// htc.resetTrackDirection();
+ double x = htc.getPosition()[0];
+ double y = htc.getPosition()[1];
+ SiSensor sensor = ((SiSensor) ((RawTrackerHit) htc.getRawHits().get(0)).getDetectorElement());
+ if (SvtUtils.getInstance().isTopLayer(sensor)) {
+ layersTop[htc.Layer() - 1]++;
+ Hep3Vector sensorPos = ((SiSensor) ((RawTrackerHit) htc.getRawHits().get(0)).getDetectorElement()).getGeometry().getPosition();
+ if (htc.Layer() == 1) {
+// System.out.println(sensorPos.toString());
+// System.out.println("Hit X = " + x + "; Hit Y = " + y);
+// aida.histogram2D("Layer 1 HTH Position: Top").fill(x - sensorPos.x(), y - sensorPos.y());
+ }
+// if (htc.Layer() == 7)
+// aida.histogram2D("Layer 7 HTH Position: Top").fill(x - sensorPos.x(), y - sensorPos.y());
+ } else {
+ layersBot[htc.Layer() - 1]++;
+ Hep3Vector sensorPos = ((SiSensor) ((RawTrackerHit) htc.getRawHits().get(0)).getDetectorElement()).getGeometry().getPosition();
+ if (htc.Layer() == 1) {
+// System.out.println(sensorPos.toString());
+// System.out.println("Hit X = " + x + "; Hit Y = " + y);
+// aida.histogram2D("Layer 1 HTH Position: Bottom").fill(x - sensorPos.x(), y - sensorPos.y());
+ }
+// if (htc.Layer() == 7)
+// aida.histogram2D("Layer 7 HTH Position: Bottom").fill(x - sensorPos.x(), y - sensorPos.y());
+ }
+ }
+
+ if (!event.hasCollection(Track.class, trackCollectionName)) {
+// System.out.println(trackCollectionName + " does not exist; skipping event");
+ // aida.histogram1D("Number Tracks/Event").fill(0);
+ return;
+ }
+
+
+ List<Track> tracks = event.get(Track.class, trackCollectionName);
+ nTracks.fill(tracks.size());
+
+
+ if (tracks.size() == 2) { //uncert can be used here && (Ytrue || Ytrue2) && (Xtrue || Xtrue2)
+
+ Track trk1 = tracks.get(0);
+ Track trk2 = tracks.get(1);
+ int isTrk1Top = -1;
+ if (trk1.getTrackerHits().get(0).getPosition()[2] > 0) {
+ isTrk1Top = 1;
+ }
+ int isTrk2Top = -1;
+ if (trk2.getTrackerHits().get(0).getPosition()[2] > 0) {
+ isTrk2Top = 1;
+ }
+ boolean topbot = false;
+ if ((isTrk1Top + isTrk2Top) == 0) {
+ topbot = true;
+ }
+
+ SeedTrack stEle1 = (SeedTrack) trk1;
+ SeedCandidate seedEle1 = stEle1.getSeedCandidate();
+ HelicalTrackFit ht1 = seedEle1.getHelix();
+ HelixConverter converter1 = new HelixConverter(0);
+ StraightLineTrack slt1 = converter1.Convert(ht1);
+
+ SeedTrack stEle2 = (SeedTrack) trk2;
+ SeedCandidate seedEle2 = stEle2.getSeedCandidate();
+ HelicalTrackFit ht2 = seedEle2.getHelix();
+ HelixConverter converter2 = new HelixConverter(0);
+ StraightLineTrack slt2 = converter2.Convert(ht2);
+
+ HPSTrack hpstrack1 = new HPSTrack(ht1);
+ Hep3Vector[] trkatconver1 = hpstrack1.getPositionAtZMap(100, zAtConverter, 1);
+ HPSTrack hpstrack2 = new HPSTrack(ht2);
+ Hep3Vector[] trkatconver2 = hpstrack2.getPositionAtZMap(100, zAtConverter, 1);
+
+
+
+ List<TrackerHit> hitsOnTrack1 = trk1.getTrackerHits();
+ int layer1;
+ double y1 = 0;
+ double y2 = 0;
+ double z1 = 0;
+ double z2 = 0;
+ double dely1 = 0;
+ double dely2 = 0;
+ for (TrackerHit hit : hitsOnTrack1) {
+ HelicalTrackHit htc1 = (HelicalTrackHit) hit;
+ layer1 = htc1.Layer();
+ int y1layer = 0;
+ int y2layer = 0;
+ if (y1 == 0) {
+ y1 = htc1.getPosition()[2]; //
+ z1 = htc1.getPosition()[0]; // z1 is jlab but the get position refers to hps-tracking
+ y1layer = layer1;
+ SymmetricMatrix ErrorHitOne = htc1.getCorrectedCovMatrix();
+ dely1 = Math.sqrt(ErrorHitOne.diagonal(2)); //y in jlab is z in hps
+ } else {
+
+ if ((layer1 > y1layer) && (y2layer == 0)) {
+ y2 = htc1.getPosition()[2]; //
+ z2 = htc1.getPosition()[0]; // see above comments!
+ y2layer = layer1;
+ SymmetricMatrix ErrorHitTwo = htc1.getCorrectedCovMatrix();
+ dely2 = Math.sqrt(ErrorHitTwo.diagonal(2));
+ }
+
+ }
+
+ }
+ List<TrackerHit> hitsOnTrack2 = trk2.getTrackerHits();
+
+ double my1 = 0;
+ double my2 = 0;
+ double mz1 = 0;
+ double mz2 = 0;
+ double delymy1 = 0;
+ double delymy2 = 0;
+ int layer2;
+ for (TrackerHit hit : hitsOnTrack2) {
+ HelicalTrackHit htc2 = (HelicalTrackHit) hit;
+// if (htc.getPosition()[2] < 0) {
+
+ layer2 = htc2.Layer();
+ int my1layer = 0;
+ int my2layer = 0;
+ if (my1 == 0) {
+ my1 = htc2.getPosition()[2]; //see above comments
+ mz1 = htc2.getPosition()[0];
+ my1layer = layer2;
+ SymmetricMatrix ErrorHitOne = htc2.getCorrectedCovMatrix();
+ delymy1 = Math.sqrt(ErrorHitOne.diagonal(2));
+ } else {
+ if ((layer2 > my1layer) && (my2layer == 0)) {
+ my2 = htc2.getPosition()[2];
+ mz2 = htc2.getPosition()[0];
+ my2layer = layer2;
+ SymmetricMatrix ErrorHitTwo = htc2.getCorrectedCovMatrix();
+ delymy2 = Math.sqrt(ErrorHitTwo.diagonal(2));
+ }
+ }
+ }
+ // double dely = .00001; //mm
+ double b1;
+ double m1;
+ double b2;
+ double m2;
+ boolean check1 = true;
+ if (y1 == 0) {
+ check1 = false;
+ }
+ boolean check2 = true;
+ if (my1 == 0) {
+ check2 = false;
+ }
+ boolean check3 = true;
+ if (my2 == 0) {
+ check3 = false;
+ }
+
+
+
+
+ //-674.062;//mm
+ double X1 = slt1.getYZAtX(zAtConverter)[0];
+ double Y1 = slt1.getYZAtX(zAtConverter)[1];
+
+ // boolean Y1top = false;
+ // boolean X1plus = false;
+ // boolean Y1bot = false;
+ // boolean X1minus = false;
+ boolean X1cent = false;
+ boolean Y1cent = false; //for simulation
+
+ if (11 < X1 && X1 < 29) {
+ X1cent = true;
+ }
+ if (-3.5 < Y1 && Y1 < 3.5) {
+ Y1cent = true;
+ }
+
+ // if (1 < Y1 && Y1 < 6) { //1 < Y1 && Y1 < 6 +-2.5
+ // Y1top = true;
+ // }
+ // if (11 < X1 && X1 < 29) { // 4 < X1 && X1 < 16 +-6
+ // X1minus = true;
+ // }
+ // if (-5 < Y1 && Y1 < 0) { // -5 < Y1 && Y1 < 0 +-2.5
+ // Y1bot = true;
+ // }
+ // if (11 < X1 && X1 < 29) { // 24 < X1 && X1 < 36 +-6
+ // X1plus = true;
+ // }
+ double X2 = slt2.getYZAtX(zAtConverter)[0];
+ double Y2 = slt2.getYZAtX(zAtConverter)[1];
+
+ // boolean Y2top = false; //for data
+ // boolean X2plus = false;
+ // boolean Y2bot = false; //in general
+ // boolean X2minus = false;
+ boolean X2cent = false;
+ boolean Y2cent = false; //for simulation
+ if (11 < X2 && X2 < 29) {
+ X2cent = true;
+ }
+ if (-3.5 < Y2 && Y2 < 3.5) {
+ Y2cent = true;
+ }
+ // if (1 < Y2 && Y2 < 6) {
+ // Y2top = true;
+ // }
+ // if (11 < X2 && X2 < 29) {
+ // X2minus = true;
+ // }
+ // if (-5 < Y2 && Y2 < 0) {
+ // Y2bot = true;
+ // }
+ // if (11 < X2 && X2 < 29) {
+ // X2plus = true;
+ // }
+
+
+ // boolean Trk1Top = false;
+ // boolean Trk2Top = false;
+ // boolean Trk1Bot = false;
+ // boolean Trk2Bot = false;
+ // if (isTrk1Top == 1) {
+ // Trk1Top = true;
+ // }
+ // if (isTrk2Top == 1) {
+ // Trk2Top = true;
+ // }
+ // if (isTrk1Top == -1) {
+ // Trk1Bot = true;
+ // }
+ // if (isTrk2Top == -1) {
+ // Trk2Bot = true;
+ // }
+ // boolean Trk1goodTop = false;
+ // boolean Trk2goodTop = false;
+ // boolean Trk1goodBot = false;
+ // boolean Trk2goodBot = false;
+ // if (Trk1Top && Y1top) {
+ // Trk1goodTop = true;
+ // }
+ // if (Trk2Top && Y2top) {
+ // Trk2goodTop = true;
+ // }
+ // if (Trk1Bot && Y1bot) {
+ // Trk1goodBot = true;
+ // }
+ // if (Trk2Bot && Y2bot) {
+ // Trk2goodBot = true;
+ // }
+
+ int qtrk1 = trk1.getCharge();
+ int qtrk2 = trk2.getCharge();
+ boolean pm = false;
+ if ((qtrk1 + qtrk2) == 0) {
+ pm = true;
+ }
+
+ // boolean Trk1Plus = false;
+ // boolean Trk2Plus = false;
+ // boolean Trk1Minus = false;
+ // boolean Trk2Minus = false;
+ // if (qtrk1 > 0) {
+ // Trk1Plus = true;
+ // } else {
+ // Trk1Minus = true;
+ // }
+ // if (qtrk2 > 0) {
+ // Trk2Plus = true;
+ // } else {
+ // Trk2Minus = true;
+ // }
+
+ // boolean Trk1goodPlus = false;
+ // boolean Trk2goodPlus = false;
+ // boolean Trk1goodMinus = false;
+ // boolean Trk2goodMinus = false;
+
+ // if (Trk1Plus && X1plus) {
+ // Trk1goodPlus = true;
+ // }
+ // if (Trk2Plus && X2plus) {
+ // Trk2goodPlus = true;
+ // }
+ // if (Trk1Minus && X1minus) {
+ // Trk1goodMinus = true;
+ // }
+ // if (Trk2Minus && X2minus) {
+ // Trk2goodMinus = true;
+ // }
+ if (topbot && pm) {
+
+ double b1p;
+ double b2p;
+ double m1p;
+ double m2p;
+
+ if (check1 && check2 && check3) {
+ if (isTrk1Top == 1) {
+ double zc = -1 * z1 / (z2 - z1);
+ b1 = (zc * (y2 - y1 + .5 * (dely2 + dely1))) + y1 - (.5 * dely1);
+ m1 = (y2 - y1 + .5 * (dely2 + dely1)) / (z2 - z1);
+ m1p = (y2 - y1 - .5 * (dely2 + dely1)) / (z2 - z1);
+ b1p = y1 - (m1p * z1) + (.5 * dely1);
+ } else {
+ double zc = -1 * z1 / (z2 - z1);
+ b1 = (zc * (y2 - y1 - .5 * (dely2 + dely1))) + y1 + (.5 * dely1);
+ m1 = (y2 - y1 - .5 * (dely2 + dely1)) / (z2 - z1);
+ m1p = (y2 - y1 + .5 * (dely2 + dely1)) / (z2 - z1);
+ b1p = y1 - (m1p * z1) - (.5 * dely1);
+ }
+
+ if (isTrk2Top == 1) {
+ double zc = -1 * mz1 / (mz2 - mz1);
+ b2 = (zc * (my2 - my1 + .5 * (delymy2 + delymy1))) + my1 - (.5 * delymy1);
+ m2 = (my2 - my1 + .5 * (delymy2 + delymy1)) / (mz2 - mz1);
+ m2p = (my2 - my1 - .5 * (delymy2 + delymy1)) / (mz2 - mz1);
+ b2p = my1 - (m2p * mz1) + (.5 * delymy1);
+ } else {
+ double zc = -1 * mz1 / (mz2 - mz1);
+ b2 = (zc * (my2 - my1 - .5 * (delymy2 + delymy1))) + my1 + (.5 * delymy1);
+ m2 = (my2 - my1 - .5 * (delymy2 + delymy1)) / (mz2 - mz1);
+ m2p = (my2 - my1 + .5 * (delymy2 + delymy1)) / (mz2 - mz1);
+ b2p = my1 - (m2p * mz1) - (.5 * delymy1);
+ }
+ // System.out.println("y1 = " + y1);
+ // System.out.println("y2 = " + y2);
+ // System.out.println("y'1 = " + my1);
+ // System.out.println("y'2 = " + my2);
+ double zi = (b2 - b1) / (m1 - m2);
+ double zr = Math.abs(zi - zAtConverter);
+ double zs = 2 * zr / 100;
+ // System.out.println("Closest Possible Z to Tracker");
+ // System.out.println(zi);
+
+
+ List<double[]> Trk1 = new ArrayList<double[]>();
+ for (int i = 0; i < 100; i++) {
+ double z = zAtConverter - zr + (zs * i);
+ double[] posvec = new double[3];
+ Hep3Vector[] trk1atz = hpstrack1.getPositionAtZMap(100, z, 1);
+ posvec[0] = trk1atz[0].x();
+ posvec[1] = trk1atz[0].y();
+ posvec[2] = z;
+
+ Trk1.add(posvec);
+ }
+ // System.out.println("Vectors ");
+
+ // System.out.println(Trk1);
+
+ double xbar = 0;
+ double ybar = 0;
+ double zbar = 0;
+ double xsqbar = 0;
+ double ysqbar = 0;
+ double zsqbar = 0;
+ int n = 0;
+ for (double[] inttrk : Trk1) {
+ // System.out.println(inttrk[0]);
+ // System.out.println(inttrk[1]);
+ // System.out.println(inttrk[2]);
+ xbar = xbar + inttrk[0];
+ ybar = ybar + inttrk[1];
+ zbar = zbar + inttrk[2];
+ n = n + 1;
+ }
+ // System.out.println("n " + n);
+ xbar = xbar / n;
+ ybar = ybar / n;
+ zbar = zbar / n;
+ // System.out.println("Xbar is " + xbar);
+ // System.out.println("Ybar is " + ybar);
+ // System.out.println("Zbar is " + zbar);
+ Matrix d;
+ Matrix A = Matrix.random(n, 3);
+ int j1 = 0;
+ for (double[] inttrk : Trk1) {
+ A.set(j1, 0, inttrk[0] - xbar);
+ A.set(j1, 1, inttrk[1] - ybar);
+ A.set(j1, 2, inttrk[2] - zbar);
+ j1++;
+ }
+
+ // System.out.println("Matrix A");
+ // A.print(9, 6);
+ A.svd();
+ SingularValueDecomposition s = A.svd();
+ Matrix S = s.getS();
+ // System.out.println("S Matrix");
+ // S.print(9, 6);
+ Matrix V = s.getV();
+ // System.out.println("V Matrix");
+ // V.print(9, 6);
+ d = V.getMatrix(0, 2, 0, 0);
+ double[] dd;
+ dd = new double[3];
+
+ dd[0] = d.get(0, 0);
+ dd[1] = d.get(1, 0);
+ dd[2] = d.get(2, 0);
+ double nd = Math.sqrt((Math.pow(dd[0], 2)) + (Math.pow(dd[1], 2)) + (Math.pow(dd[2], 2)));
+
+ for (double[] inttrk : Trk1) {
+ double t1 = (inttrk[2] - zbar) / dd[2];
+ double restrk1[];
+ restrk1 = new double[3];
+ restrk1[0] = xbar + (t1) * dd[0] - inttrk[0];
+ restrk1[1] = ybar + (t1) * dd[1] - inttrk[1];
+ restrk1[2] = zbar + (t1) * dd[2] - inttrk[2];
+ aida.histogram1D("X Res Trk1").fill(restrk1[0]);
+ aida.histogram1D("Y Res Trk1").fill(restrk1[1]);
+ }
+
+ List<double[]> Trk2 = new ArrayList<double[]>();
+ for (int i = 0; i < 100; i++) {
+ double z = zAtConverter - zr + (zs * i);
+ double[] posvec2 = new double[3];
+ Hep3Vector[] trk2atz = hpstrack2.getPositionAtZMap(100, z, 1);
+ posvec2[0] = trk2atz[0].x();
+ posvec2[1] = trk2atz[0].y();
+ posvec2[2] = z;
+ Trk2.add(posvec2);
+ // System.out.println("Components");
+ // System.out.println(posvec2[0]);
+ // System.out.println(posvec2[1]);
+ // System.out.println(posvec2[2]);
+ }
+ double xbar2 = 0;
+ double ybar2 = 0;
+ double zbar2 = 0;
[truncated at 1000 lines; 562 more skipped]