hps-java/src/main/java/org/lcsim/HPSDedicatedv2
diff -N TestConstrainedVertex.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TestConstrainedVertex.java 24 Jun 2010 20:22:53 -0000 1.1
@@ -0,0 +1,1736 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.HPSDedicatedv2;
+
+import Utilities.FindableTrack;
+import Utilities.TrackAnalysis;
+import java.io.IOException;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+import hep.aida.*;
+
+import hep.physics.matrix.BasicMatrix;
+import hep.physics.matrix.Matrix;
+import hep.physics.matrix.MutableMatrix;
+import hep.physics.vec.VecOp;
+import java.io.File;
+import java.io.FileWriter;
+import java.io.PrintWriter;
+import java.util.ArrayList;
+import java.util.Set;
+
+//import org.lcsim.contrib.sATLAS.FindableTrack;
+//import org.lcsim.contrib.sATLAS.TrackAnalysis;
+import org.lcsim.HPSVertexing.HelixConverter;
+import org.lcsim.HPSVertexing.StraightLineTrack;
+import org.lcsim.HPSVertexing.VertexFit;
+import org.lcsim.HPSVertexing.VertexFitter;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.Track;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.fit.helicaltrack.HelicalTrackCross;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
+import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.fit.helicaltrack.TrackDirection;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+import org.lcsim.recon.tracking.seedtracker.StrategyXMLUtils;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author partridge
+ */
+public class TestConstrainedVertex extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ private IAnalysisFactory af = aida.analysisFactory();
+ private ITupleFactory tf = af.createTupleFactory(aida.tree());
+ private ITuple tuple;
+ private IProfile1D phifake;
+ private IProfile1D pfake;
+ private IProfile1D cthfake;
+ private IProfile1D peffFindable;
+ private IProfile1D thetaeffFindable;
+ private IProfile1D phieffFindable;
+ private IProfile1D ctheffFindable;
+ private IProfile1D d0effFindable;
+ private IProfile1D z0effFindable;
+ private IProfile1D peffElectrons;
+ private IProfile1D thetaeffElectrons;
+ private IProfile1D phieffElectrons;
+ private IProfile1D ctheffElectrons;
+ private IProfile1D d0effElectrons;
+ private IProfile1D z0effElectrons;
+ private IProfile1D peffAxial;
+ private IProfile1D thetaeffAxial;
+ private IProfile1D phieffAxial;
+ private IProfile1D ctheffAxial;
+ private IProfile1D d0effAxial;
+ private IProfile1D z0effAxial;
+ private IHistogram1D fakes;
+ private IHistogram1D nfakes;
+ private IProfile1D HitEffdEdX;
+ private IProfile1D ClHitEffdEdX;
+ private IProfile1D ClHitEffY;
+ private IProfile1D ClHitEffZ;
+ private IProfile1D STdEdXY;
+ private IProfile1D STdEdXZ;
+ private IProfile1D frdEdXY;
+ private IProfile1D VxEff;
+ private IProfile1D VyEff;
+ private IProfile1D VzEff;
+ private IProfile1D VxEffFindable;
+ private IProfile1D VyEffFindable;
+ private IProfile1D VzEffFindable;
+ public String outputPlots = "myplots.aida";
+ Map<String, IProfile1D> clsizeMap = new HashMap<String, IProfile1D>();
+ String[] detNames = {"Tracker"};
+ Integer[] nlayers = {8};
+ int trk_count = 0;
+ int nevt = 0;
+ int _nmcTrk = 0;
+ double _nrecTrk = 0;
+ double phiTrkCut = 0.3;
+ double cosThTrkCutMax = 0.2;
+ double cosThTrkCutMin = 0.05;
+ double pTrkCut = 0.5; //GeV
+ double d0TrkCut = 2.0; //mm
+ double z0TrkCut = 2.0; //mm
+ double etaTrkCut = 2.5;
+ int totelectrons = 0;
+ double foundelectrons = 0;
+ int findableelectrons = 0;
+ int findableTracks = 0;
+ double foundTracks = 0;
+ double xref = 50.0; //mm
+ public String outputTextName = "myevents.txt";
+ FileWriter fw;
+ PrintWriter pw;
+
+ public TestConstrainedVertex(int layers) {
+ nlayers[0] = layers;
+ String names = new String(
+ "float chi2; float Xoca; float Yoca; float Zoca; float Px; float Py; float Pz" +
+ "; float mcXoca; float mcYoca; float mcZoca; float mcPx; float mcPy; float mcPz" +
+ "; int charge; int nbad; int nhits; ITuple layerTuple={int layer, int nstrips, int bad" +
+ ", boolean isAxial, float htsY, float htsZ, float sthY, float sthZ}" +
+ "; ITuple planeTuple={int plane, int bad" +
+ ", float hthY, float hthZ, float hthYErr, float hthZErr, float trkY, float trkZ}");
+ tuple = tf.create("trackFitterResults", "properties of track", names, "");
+ // Define the efficiency histograms
+ IHistogramFactory hf = aida.histogramFactory();
+
+
+ peffFindable = hf.createProfile1D("Findable Efficiency vs p", "", 20, 0., 6.);
+ thetaeffFindable = hf.createProfile1D("Findable Efficiency vs theta", "", 20, 80, 100);
+ phieffFindable = hf.createProfile1D("Findable Efficiency vs phi", "", 25, -0.25, 0.25);
+ ctheffFindable = hf.createProfile1D("Findable Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
+ d0effFindable = hf.createProfile1D("Findable Efficiency vs d0", "", 50, -2., 2.);
+ z0effFindable = hf.createProfile1D("Findable Efficiency vs z0", "", 50, -2., 2.);
+
+ peffElectrons = hf.createProfile1D("Electrons Efficiency vs p", "", 20, 0., 6.);
+ thetaeffElectrons = hf.createProfile1D("Electrons Efficiency vs theta", "", 20, 80, 100);
+ phieffElectrons = hf.createProfile1D("Electrons Efficiency vs phi", "", 25, -0.25, 0.25);
+ ctheffElectrons = hf.createProfile1D("Electrons Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
+ d0effElectrons = hf.createProfile1D("Electrons Efficiency vs d0", "", 20, -1., 1.);
+ z0effElectrons = hf.createProfile1D("Electrons Efficiency vs z0", "", 20, -1., 1.);
+
+ peffAxial = hf.createProfile1D("Axial Efficiency vs p", "", 20, 0., 6.);
+ thetaeffAxial = hf.createProfile1D("Axial Efficiency vs theta", "", 20, 80, 100);
+ phieffAxial = hf.createProfile1D("Axial Efficiency vs phi", "", 25, -0.25, 0.25);
+ ctheffAxial = hf.createProfile1D("Axial Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
+ d0effAxial = hf.createProfile1D("Axial Efficiency vs d0", "", 20, -1., 1.);
+ z0effAxial = hf.createProfile1D("Axial Efficiency vs z0", "", 20, -1., 1.);
+
+ cthfake = hf.createProfile1D("Fake rate vs cos(theta)", "", 25, -0.25, 0.25);
+ phifake = hf.createProfile1D("Fake rate vs phi", "", 25, -0.25, 0.25);
+ pfake = hf.createProfile1D("Fake rate vs p", "", 20, 0, 6);
+
+ fakes = hf.createHistogram1D("Number of mis-matched hits (unnormalized)", "", 10, 0., 10.);
+ nfakes = hf.createHistogram1D("Number of mis-matched hits (normalized)", "", 10, 0., 10.);
+
+ HitEffdEdX = hf.createProfile1D("Strip Hit Efficiency vs dEdX", "", 50, 0, 0.3);
+ ClHitEffdEdX = hf.createProfile1D("Cluster Hit Efficiency vs dEdX", "", 50, 0, 0.3);
+ ClHitEffY = hf.createProfile1D("Cluster Hit Efficiency vs y", "", 50, -100, 100);
+ ClHitEffZ = hf.createProfile1D("Cluster Hit Efficiency vs z", "", 50, -100, 100);
+ STdEdXY = hf.createProfile1D("SimTHit dEdX vs y", "", 50, -100, 100);
+ frdEdXY = hf.createProfile1D("fractional dEdX vs y", "", 50, -100, 100);
+ STdEdXZ = hf.createProfile1D("SimTHit dEdX vs z", "", 50, -100, 100);
+
+ VxEff = hf.createProfile1D("Aprime Efficiency vs Vx", "", 25, 0., 50.);
+ VyEff = hf.createProfile1D("Aprime Efficiency vs Vy", "", 40, -0.2, 0.2);
+ VzEff = hf.createProfile1D("Aprime Efficiency vs Vz", "", 40, -0.2, 0.2);
+
+ VxEffFindable = hf.createProfile1D("Aprime Efficiency vs Vx: Findable", "", 25, 0., 50.);
+ VyEffFindable = hf.createProfile1D("Aprime Efficiency vs Vy: Findable", "", 40, -0.2, 0.2);
+ VzEffFindable = hf.createProfile1D("Aprime Efficiency vs Vz: Findable", "", 40, -0.2, 0.2);
+
+ int i, j;
+ for (i = 0; i < 1; i++)
+ for (j = 0; j < nlayers[i]; j++) {
+ int laynum = j + 1;
+ String profname = detNames[i] + "_layer" + laynum + " cluster size vs y";
+ String key = detNames[i] + "_layer" + laynum;
+ clsizeMap.put(key, hf.createProfile1D(profname, 20, -15, 15));
+ }
+ }
+
+ @Override
+ public void process(
+ EventHeader event) {
+ if (nevt == 0)
+ try {
+//open things up
+ fw = new FileWriter(outputTextName);
+ pw = new PrintWriter(fw);
+ } catch (IOException ex) {
+ Logger.getLogger(TestConstrainedVertex.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ // Increment the event counter
+ nevt++;
+ String resDir = "residualsPlots/";
+ String resDirBar = "residualsBarrelPlots/";
+ String resDirEC = "residualsEndcapPlots/";
+ String simDir = "STHitPlots/";
+ String debugDir = "debugPlots/";
+ String occDir = "occupancyPlots/";
+ // Get the magnetic field
+ Hep3Vector IP = new BasicHep3Vector(200., 0., 0.);
+ double bfield = event.getDetector().getFieldMap().getField(IP).z();
+
+ List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, "RawTrackerHitMaker_RawTrackerHits");
+ List<SiTrackerHitStrip1D> stripHits = event.get(SiTrackerHitStrip1D.class, "StripClusterer_SiTrackerHitStrip1D");
+ // dump SThit information
+ String[] input_hit_collections = {"TrackerHits"};
+ for (String input : input_hit_collections) {
+ List<SimTrackerHit> sthits = event.getSimTrackerHits(input);
+ int[] nhits = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ for (SimTrackerHit st : sthits) {
+ String detector = st.getDetectorElement().getName();
+ int layer = st.getLayerNumber();
+ double[] hp = st.getPoint();
+ Hep3Vector hitPos = new BasicHep3Vector(hp[0], hp[1], hp[2]);
+ double r = Math.sqrt(hp[0] * hp[0] + hp[1] * hp[1]);
+ double theta = Math.atan2(r, hp[2]);
+// double eta = -Math.log(r / (2 * Math.abs(hp[2])));
+ double eta = -Math.log(Math.tan(theta / 2));
+ double phi = Math.atan2(hp[1], hp[0]);
+// System.out.println("r= " + r + " theta = "+theta+" eta = " + eta+ " phi=" + phi);
+ nhits[layer]++;
+ //double mom=0,dedx=0;
+ //if(st.getMCParticle()!=null){
+ // mom = st.getMCParticle().getMomentum().magnitude();
+ // dedx = st.getdEdx() * 1000.0;
+ // }
+ double hitwgt = 0;
+ double clhitwgt = 0;
+ for (RawTrackerHit rth : rawHits) {
+ List<SimTrackerHit> SthFromRth = rth.getSimTrackerHits();
+ if (SthFromRth.contains(st))
+ hitwgt = 1.0;
+ }
+ for (SiTrackerHitStrip1D cluster : stripHits) {
+ double measdedx = cluster.getdEdx() * 1000.0;
+
+ List<RawTrackerHit> RthFromSith = cluster.getRawHits();
+
+ for (RawTrackerHit rth : RthFromSith) {
+ List<SimTrackerHit> SthFromRth = rth.getSimTrackerHits();
+ if (SthFromRth.contains(st)) {
+ clhitwgt = 1.0;
+ double totdedx = 0;
+ for (SimTrackerHit sthtemp : SthFromRth)
+ totdedx = totdedx + sthtemp.getdEdx() * 1000.0;
+ aida.histogram1D(simDir + "delta dEdX", 50, -0.2, 0.2).fill(measdedx - totdedx);
+ aida.histogram1D(simDir + "fractional dEdX", 50, -1, 1.).fill((measdedx - totdedx) / totdedx);
+ aida.cloud1D(simDir + "fractional dEdX Cloud").fill((measdedx - totdedx) / totdedx);
+ // if (Math.abs((measdedx - dedx) / dedx) < 1)
+ frdEdXY.fill(hp[1], (measdedx - totdedx) / totdedx);
+ // if (dedx == 0)
+ // System.out.println("***************** dedx==0 ********");
+ }
+ }
+ }
+ //HitEffdEdX.fill(dedx, hitwgt);
+ //ClHitEffdEdX.fill(dedx, clhitwgt);
+ ClHitEffY.fill(hp[1], clhitwgt);
+ ClHitEffZ.fill(hp[2], clhitwgt);
+ //STdEdXY.fill(hp[1], dedx);
+ //STdEdXZ.fill(hp[2], dedx);
+ //aida.histogram1D(simDir + " dedx", 50, 0, 0.3).fill(dedx);
+// if (hitwgt == 0) {
+// System.out.println("TrackAnalysis: found an inefficiency hit: " + dedx);
+// }
+
+ //aida.cloud1D(simDir + input + " layer " + layer + " STHit p").fill(mom);
+ aida.cloud1D(simDir + input + " layer " + layer + " STHit y").fill(hp[1]);
+ aida.cloud1D(simDir + input + " layer " + layer + " STHit z").fill(hp[2]);
+ aida.cloud2D(simDir + input + " layer " + layer + " STHit y vs z").fill(hp[2], hp[1]);
+ aida.histogram2D(simDir + input + " layer " + layer + " STHit y vs z occupancy", 100, -15, 15, 500, -15, 15).fill(hp[2], hp[1]);
+
+ }
+ int i = 0;
+ while (i < nlayers[0]) {
+ if (nhits[i] > 0)
+ aida.cloud1D(simDir + input + "layer " + i + " number of ST hits").fill(nhits[i]);
+ i++;
+ }
+ }
+
+
+// List<HelicalTrackHit> hthits = event.get(HelicalTrackHit.class, "MatchedHTHits");
+// String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "DarkPhoton-Final.xml";
+ String strategyPrefix = "/nfs/sulky21/g.ec.u12/users/mgraham/AtlasUpgrade/hps-java/src/main/resources/";
+ String sfile = "DarkPhoton-Final.xml";
+// List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(sfile);
+ List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromFile(new File(strategyPrefix + sfile));
+ List<HelicalTrackHit> toththits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
+ List<HelicalTrackHit> axialhits = event.get(HelicalTrackHit.class, "AxialTrackHits");
+
+
+ Map<String, Integer> occupancyMap = new HashMap<String, Integer>();
+ for (RawTrackerHit rh : rawHits) {
+ IDetectorElement rhDetE = rh.getDetectorElement();
+
+ String rhDetName = rhDetE.getName();
+// System.out.println(rhDetName);
+ int rhLayer = rh.getLayerNumber();
+// String[] shortrhDetName=rhDetName.split("^[A-Z]+_layer[0-9]");
+
+ for (String myname : detNames)
+ if (rhDetName.contains(myname)) {
+ String detlayer = myname + "_" + rhLayer;
+ Integer myint = occupancyMap.get(detlayer);
+ if (myint == null)
+ myint = 1;
+ myint++;
+ occupancyMap.put(detlayer, myint);
+ }
+ }
+ Set<String> mykeyset = (Set<String>) occupancyMap.keySet();
+ for (String keys : mykeyset)
+ aida.cloud1D(occDir + keys + " # of hits").fill(occupancyMap.get(keys));
+ int detNum = 0;
+ int layerNum = 0;
+
+ for (SiTrackerHitStrip1D stripCluster : stripHits) {
+ Hep3Vector strCluPos = stripCluster.getPositionAsVector();
+ double xHit = strCluPos.x();
+ double yHit = strCluPos.y();
+
+ double rHit = Math.sqrt(strCluPos.x() * strCluPos.x() + strCluPos.y() * strCluPos.y());
+ double zHit = strCluPos.z();
+ double etaHit = -Math.log(Math.tan(Math.atan2(rHit, zHit) / 2));
+ Set<MCParticle> mcparts = stripCluster.getMCParticles();
+ aida.cloud1D(occDir + "associated MC Particles").fill(mcparts.size());
+
+
+ List<RawTrackerHit> rthList = stripCluster.getRawHits();
+ int nhits = rthList.size();
+ String detlayer = "Foobar";
+ for (RawTrackerHit rth : rthList) {
+ IDetectorElement rhDetE = rth.getDetectorElement();
+ String rhDetName = rhDetE.getName();
+ int rhLayer = rth.getLayerNumber();
+ for (String myname : detNames)
+ if (rhDetName.contains(myname))
+ detlayer = myname + "_layer" + rhLayer;
+ }
+// System.out.println(detlayer);
+ clsizeMap.get(detlayer).fill(yHit, nhits);
+ aida.cloud1D(occDir + detlayer + "associated MC Particles").fill(mcparts.size());
+ aida.cloud1D(occDir + detlayer + " cluster size").fill(nhits);
+ }
+
+ // Find the minimum pT among the strategies
+ double ptCut = 9999.;
+ for (SeedStrategy s : slist)
+ if (s.getMinPT() < ptCut)
+ ptCut = s.getMinPT();
+
+ // Create a relational table that maps TrackerHits to MCParticles
+ RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
+
+ for (LCRelation relation : mcrelations)
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+ hittomc.add(relation.getFrom(), relation.getTo());
+
+ RelationalTable hittomcAxial = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+// List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
+ List<LCRelation> mcrelationsAxial = event.get(LCRelation.class, "AxialTrackMCRelations");
+ for (LCRelation relation : mcrelationsAxial)
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+ hittomcAxial.add(relation.getFrom(), relation.getTo());
+
+ // Instantiate the class that determines if a track is "findable"
+ FindableTrack findable = new FindableTrack(event);
+
+ // Create a map between tracks and the associated MCParticle
+// List<Track> tracklist = event.getTracks();
+ List<Track> tracklist = event.get(Track.class, "MatchedTracks");
+ List<Track> axialtracklist = event.get(Track.class, "AxialTracks");
+
+ aida.cloud1D("Matched Tracks per Event").fill(tracklist.size());
+ aida.cloud1D("Axial Tracks per Event").fill(axialtracklist.size());
+ aida.cloud1D("HelicalTrackHits per Event").fill(toththits.size());
+ RelationalTable trktomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ RelationalTable trktomcAxial = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+
+ int _nchRec = 0;
+ int _nchRecBar = 0;
+ int _neleRec = 0;
+ int _nposRec = 0;
+ int _neleTru = 0;
+ int _nposTru = 0;
+ int _neleFake = 0;
+ int _nposFake = 0;
+ int totcharge = 0;
+
+ MutableMatrix BS = new BasicMatrix(4, 4);
+ BS.setElement(StraightLineTrack.y0Index, StraightLineTrack.y0Index, 0.01 * 0.01);
+ BS.setElement(StraightLineTrack.z0Index, StraightLineTrack.z0Index, 0.01 * 0.01);
+ BS.setElement(StraightLineTrack.dydxIndex, StraightLineTrack.dydxIndex, 1e-20);
+ BS.setElement(StraightLineTrack.dzdxIndex, StraightLineTrack.dzdxIndex, 1e-20);
+
+ SymmetricMatrix beamSpotCov = new SymmetricMatrix(BS);
+//make the beamspot centered at (0,0,0) with a 10um width in y,z
+ StraightLineTrack beamSpot = new StraightLineTrack(0, 0, 0, 0, 0, beamSpotCov);
+
+ RelationalTable mcHittomcP = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+
+ // Get the collections of SimTrackerHits
+ List<List<SimTrackerHit>> simcols = event.get(SimTrackerHit.class);
+
+ // Loop over the SimTrackerHits and fill in the relational table
+ for (List<SimTrackerHit> simlist : simcols)
+ for (SimTrackerHit simhit : simlist)
+ if (simhit.getMCParticle() != null)
+ mcHittomcP.add(simhit, simhit.getMCParticle());
+
+ Map<Track, TrackAnalysis> tkanalMap = new HashMap<Track, TrackAnalysis>();
+ Map<Track, StraightLineTrack> sltMap = new HashMap<Track, StraightLineTrack>();
+ RelationalTable nearestHitToTrack = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ Map<Track, Double> l1Isolation = new HashMap<Track, Double>();
+ String atrackdir = "TrackInfoAxial/";
+ // Analyze the tracks in the event
+ for (Track atrack : axialtracklist) {
+ double apx = atrack.getPX();
+ aida.cloud1D(atrackdir + "pX").fill(apx);
+ }
+ String trackdir = "TrackInfo/";
+ // Analyze the tracks in the event
+ for (Track track : tracklist) {
+ // Calculate the track pT and cos(theta)
+ double d0 = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
+ double phi0 = track.getTrackParameter(HelicalTrackFit.phi0Index);
+ double slope = track.getTrackParameter(HelicalTrackFit.slopeIndex);
+ double curve = track.getTrackParameter(HelicalTrackFit.curvatureIndex);
+ double d0Err = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.dcaIndex, HelicalTrackFit.dcaIndex));
+ double z0Err = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.z0Index, HelicalTrackFit.z0Index));
+ double phi0Err = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.phi0Index, HelicalTrackFit.phi0Index));
+ double slopeErr = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.slopeIndex, HelicalTrackFit.slopeIndex));
+ double curveErr = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.curvatureIndex, HelicalTrackFit.curvatureIndex));
+ _nchRec++;
+ if (track.getCharge() < 0)
+ _neleRec++;
+ if (track.getCharge() > 0)
+ _nposRec++;
+ double mom[] = track.getMomentum();
+
+ //extrapolate straight back...
+
+ SeedTrack stEle = (SeedTrack) track;
+ SeedCandidate seedEle = stEle.getSeedCandidate();
+ HelicalTrackFit ht = seedEle.getHelix();
+ HelixConverter converter = new HelixConverter(xref);
+ StraightLineTrack slt = converter.Convert(ht);
+
+
+ TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
+
+ tkanalMap.put(track, tkanal);
+ sltMap.put(track, slt);
+
+ double truedoca = slt.Doca();
+ double xoca = slt.xPoca();
+ double[] poca = slt.Poca();
+ double y0new = slt.TargetYZ()[0];
+ double z0new = slt.TargetYZ()[1];
+
+ Hep3Vector trueP = getTrueMomentum(track, slt);
+ double px = trueP.x();
+ double py = trueP.y();
+ double pz = trueP.z();
+ double pperp = Math.sqrt(py * py + pz * pz);
+ double pt = Math.sqrt(px * px + py * py);
+ double p = Math.sqrt(pt * pt + pz * pz);
+ double phi = Math.atan2(py, px);
+ double cth = pz / Math.sqrt(pt * pt + pz * pz);
+ double sth = pt / Math.sqrt(pt * pt + pz * pz);
+ double th = Math.atan2(pt, pz);
+ double eta = -Math.log(Math.tan(th / 2));
+
+
+
+ double s = HelixUtils.PathLength(ht, (HelicalTrackHit) track.getTrackerHits().get(0));
+ double y1 = HelixUtils.PointOnHelix(ht, s).y();
+ double z1 = HelixUtils.PointOnHelix(ht, s).z();
+
+ int nhits = tkanal.getNHitsNew();
+ double purity = tkanal.getPurityNew();
+ int nbad = tkanal.getNBadHitsNew();
+ int nbadAxial = tkanal.getNBadAxialHits();
+ int nbadZ = tkanal.getNBadZHits();
+ int nAxial = tkanal.getNAxialHits();
+ int nZ = tkanal.getNZHits();
+ List<Integer> badLayers = tkanal.getBadHitList();
+ Integer badLayerEle = encodeBadHitList(badLayers);
+ if (badLayers.size() > 0) {
+ System.out.println(badLayers.toString());
+ System.out.println("Bad Layer code: " + badLayerEle);
+ }
+ aida.cloud1D(trackdir + "Mis-matched hits for all tracks").fill(nbad);
+ aida.cloud1D(trackdir + "purityNew for all tracks").fill(purity);
+ aida.cloud1D(trackdir + "Bad Axial hits for all tracks").fill(nbadAxial);
+ aida.cloud1D(trackdir + "Bad Z hits for all tracks").fill(nbadZ);
+ aida.cloud1D(trackdir + "Number of Axial hits for all tracks").fill(nAxial);
+ aida.cloud1D(trackdir + "Number of Z hits for all tracks").fill(nZ);
+
+ for (Integer bhit : badLayers)
+ aida.histogram1D(trackdir + "Layer of Bad Hit", nlayers[0], 1, nlayers[0] + 1).fill(bhit);
+
+ // Generate a normalized histogram after 1000 events
+ trk_count++;
+ if (nevt <= 1000)
+ fakes.fill(nbad);
+
+ // Make plots for fake, non-fake, and all tracks
+ if (purity < 0.5) {
+ if (track.getCharge() < 0)
+ _neleFake++;
+ if (track.getCharge() > 0)
+ _nposFake++;
+ cthfake.fill(cth, 1.0);
+ phifake.fill(phi, 1.0);
+ pfake.fill(p, 1.0);
+
+ fillTrackInfo(trackdir, "fake tracks", track.getChi2(), nhits, p, pperp, px, py, pz, phi, cth, truedoca, xoca, poca[1], poca[2]);
+
+ } else {
+ if (track.getCharge() < 0)
+ _neleTru++;
+ if (track.getCharge() > 0)
+ _nposTru++;
+ cthfake.fill(cth, 0.0);
+ phifake.fill(phi, 0.0);
+ pfake.fill(p, 0.0);
+
+ fillTrackInfo(trackdir, "non-fake tracks", track.getChi2(), nhits, p, pperp, px, py, pz, phi, cth, truedoca, xoca, poca[1], poca[2]);
+
+
+
+ }
+ fillTrackInfo(trackdir, "all tracks", track.getChi2(), nhits, p, pperp, px, py, pz, phi, cth, truedoca, xoca, poca[1], poca[2]);
+ if (nbadZ == 3)
+ fillTrackInfo(trackdir, "3 Bad Z-hits", track.getChi2(), nhits, p, pperp, px, py, pz, phi, cth, truedoca, xoca, poca[1], poca[2]);
+
+
+ // Now analyze MC Particles on this track
+ MCParticle mcp = tkanal.getMCParticle();
+ if (mcp != null) {
+
+ // Create a map between the tracks found and the assigned MC particle
+ trktomc.add(track, tkanal.getMCParticle());
+
+ // Calculate the MC momentum and polar angle
+ Hep3Vector Pmc = mcp.getMomentum();
+ double pxmc = Pmc.x();
+ double pymc = Pmc.y();
+ double pzmc = Pmc.z();
+ double ptmc = Math.sqrt(pxmc * pxmc + pymc * pymc);
+ double pmc = Math.sqrt(ptmc * ptmc + pzmc * pzmc);
+ double pxtk = track.getPX();
+ double pytk = track.getPY();
+ double pztk = track.getPZ();
+ double pttk = Math.sqrt(pxtk * pxtk + pytk * pytk);
+ double ptk = Math.sqrt(pttk * pttk + pztk * pztk);
+
+ // Calculate the helix parameters for this MC particle and pulls in pT, d0
+ HelixParamCalculator helix = new HelixParamCalculator(mcp, bfield);
+ double d0mc = helix.getDCA();
+ double z0mc = helix.getZ0();
+ double phi0mc = helix.getPhi0();
+ double slopemc = helix.getSlopeSZPlane();
+ double curvemc = 1 / helix.getRadius();
+ double pinvresid = (1 / ptk - 1 / pmc);
+ double presid = (ptk - pmc);
+ double z0newMC = z0mc;
+ double y0newMC = d0mc;
+ aida.histogram1D(trackdir + "d0 Pull", 200, -8, 8).fill((d0 - d0mc) / d0Err);
+ aida.histogram1D(trackdir + "z0 Pull", 200, -8, 8).fill((z0 - z0mc) / z0Err);
+ aida.histogram1D(trackdir + "phi0 Pull", 200, -8, 8).fill((phi0 - phi0mc) / phi0Err);
+ aida.histogram1D(trackdir + "slope Pull", 200, -8, 8).fill((slope - slopemc) / slopeErr);
+ aida.histogram1D(trackdir + "curvature Pull", 200, -8, 8).fill((curve - curvemc) / curveErr);
+
+ double truedocaMC = findDoca(y0newMC, z0newMC, pxmc, pymc, pzmc);
+ double[] pocaMC = findPoca(y0newMC, z0newMC, pxmc, pymc, pzmc);
+ double z0Newresid = z0new - z0newMC;
+ double y0Newresid = y0new - y0newMC;
+ double docaresid = truedoca - truedocaMC;
+
+ double[] pocaresid = {0, 0, 0};
+ pocaresid[0] = poca[0] - pocaMC[0];
+ pocaresid[1] = poca[1] - pocaMC[1];
+ pocaresid[2] = poca[2] - pocaMC[2];
+ // Plot the pt and d0 pulls for various purity intervals
+ if (nbad == 0)
+ fillClouds("residuals/", "0 Bad Hits", nhits, p, cth, y0new, z0new, truedoca, poca, pinvresid, presid, y0Newresid, z0Newresid, docaresid, pocaresid);
+ else if (purity > 0.5)
+ fillClouds("residuals/", "purity.gt.0.5", nhits, p, cth, y0new, z0new, truedoca, poca, pinvresid, presid, y0Newresid, z0Newresid, docaresid, pocaresid);
+ else if (purity < 0.5)
+ fillClouds("residuals/", "purity.lt.0.5", nhits, p, cth, y0new, z0new, truedoca, poca, pinvresid, presid, y0Newresid, z0Newresid, docaresid, pocaresid);
+ if (nbadAxial == 0)
+ fillClouds("residuals/Axial/", "0 Bad Axial Hits", nhits, p, cth, y0new, z0new, truedoca, poca, pinvresid, presid, y0Newresid, z0Newresid, docaresid, pocaresid);
+ else if (nbadAxial == 1)
+ fillClouds("residuals/Axial/", "1 Bad Axial Hits", nhits, p, cth, y0new, z0new, truedoca, poca, pinvresid, presid, y0Newresid, z0Newresid, docaresid, pocaresid);
+ else
+ fillClouds("residuals/Axial/", "2 or more Bad Axial Hits", nhits, p, cth, y0new, z0new, truedoca, poca, pinvresid, presid, y0Newresid, z0Newresid, docaresid, pocaresid);
+ if (nbadZ == 0)
+ fillClouds("residuals/Z/", "0 Bad Z Hits", nhits, p, cth, y0new, z0new, truedoca, poca, pinvresid, presid, y0Newresid, z0Newresid, docaresid, pocaresid);
+ else if (nbadZ == 1)
+ fillClouds("residuals/Z/", "1 Bad Z Hits", nhits, p, cth, y0new, z0new, truedoca, poca, pinvresid, presid, y0Newresid, z0Newresid, docaresid, pocaresid);
+ else
+ fillClouds("residuals/Z/", "2 or more Bad Z Hits", nhits, p, cth, y0new, z0new, truedoca, poca, pinvresid, presid, y0Newresid, z0Newresid, docaresid, pocaresid);
+
+
+ BasicHep3Vector axial = new BasicHep3Vector();
+ axial.setV(0, 1, 0);
+ String hitdir = "HitsOnTrack/";
+ List<TrackerHit> hitsOnTrack = track.getTrackerHits();
+ MCParticle bestmcp = tkanal.getMCParticleNew();
+
+ String tkresid = "TrackResiduals/";
+
+ int ndaug = 0;
+ if (bestmcp != null) ndaug = bestmcp.getDaughters().size();
+ int imain = 0;
+ tuple.fill(imain++, (float) track.getChi2());
+ tuple.fill(imain++, (float) poca[0]);
+ tuple.fill(imain++, (float) poca[1]);
+ tuple.fill(imain++, (float) poca[2]);
+ tuple.fill(imain++, (float) px);
+ tuple.fill(imain++, (float) py);
+ tuple.fill(imain++, (float) pz);
+ tuple.fill(imain++, (float) pocaMC[0]);
+ tuple.fill(imain++, (float) pocaMC[1]);
+ tuple.fill(imain++, (float) pocaMC[2]);
+ tuple.fill(imain++, (float) pxmc);
+ tuple.fill(imain++, (float) pymc);
+ tuple.fill(imain++, (float) pzmc);
+ tuple.fill(imain++, (int) track.getCharge());
+ tuple.fill(imain++, (int) nbad);
+ tuple.fill(imain++, (int) nhits);
+ double mcmom = 0;
+ double prevmom = 0;
+ ITuple layerTuple = tuple.getTuple(imain);
+ ITuple planeTuple = tuple.getTuple(imain + 1);
+ for (TrackerHit hit : hitsOnTrack) {
+
+ int iplane = 0;
+ HelicalTrackHit htc = (HelicalTrackHit) hit;
+ List<MCParticle> mcpsHTH = htc.getMCParticles();
+ int isbad = 0;
+ if (mcpsHTH.size() == 0 || mcpsHTH.size() > 1 || !mcpsHTH.contains(bestmcp))
+ isbad = 1;
+ double sHit = ht.PathMap().get(htc);
+ Hep3Vector posonhelix = HelixUtils.PointOnHelix(ht, sHit);
+ double yTr = posonhelix.y();
+ double zTr = posonhelix.z();
+ HelicalTrackCross cross = (HelicalTrackCross) htc;
+ List<HelicalTrackStrip> clusterlist = cross.getStrips();
+ TrackDirection trkdir = HelixUtils.CalculateTrackDirection(ht, sHit);
+ cross.setTrackDirection(trkdir, ht.covariance());
+ double y = cross.y();
+ double z = cross.z();
+ double yerr = Math.sqrt(cross.getCorrectedCovMatrix().e(1, 1));
+ double zerr = Math.sqrt(cross.getCorrectedCovMatrix().e(2, 2));
+
+ int htlayer = htc.Layer();
+ planeTuple.fill(iplane++, (int) htlayer);
+ planeTuple.fill(iplane++, (int) isbad);
+ planeTuple.fill(iplane++, (float) y);
+ planeTuple.fill(iplane++, (float) z);
+ planeTuple.fill(iplane++, (float) yerr);
+ planeTuple.fill(iplane++, (float) zerr);
+ planeTuple.fill(iplane++, (float) yTr);
+ planeTuple.fill(iplane++, (float) zTr);
+ planeTuple.addRow();
+ if (purity == 1 && track.getCharge() > 0 && nhits == 10) {
+ if (clusterlist.get(0).rawhits().size() == 1 && clusterlist.get(1).rawhits().size() == 1) {
+ aida.cloud1D(hitdir + tkresid + "SingleStrip--Track delta y: Layer " + htlayer).fill(y - yTr);
+ aida.cloud1D(hitdir + tkresid + "SingleStrip--Track delta z: Layer " + htlayer).fill(z - zTr);
+ }
+ aida.cloud1D(hitdir + tkresid + " Measured y: Layer " + htlayer).fill(y);
+ aida.cloud1D(hitdir + tkresid + " Track y: Layer " + htlayer).fill(yTr);
+ aida.cloud1D(hitdir + tkresid + " Measured z: Layer " + htlayer).fill(z);
+ aida.cloud1D(hitdir + tkresid + " Track z: Layer " + htlayer).fill(zTr);
+ aida.cloud1D(hitdir + tkresid + " Measured y ").fill(y);
+ aida.cloud1D(hitdir + tkresid + " Track delta y: Layer " + htlayer + "; ndaug=" + ndaug).fill(y - yTr);
+ aida.cloud1D(hitdir + tkresid + " Track delta z: Layer " + htlayer + "; ndaug=" + ndaug).fill(z - zTr);
+ aida.cloud2D(hitdir + tkresid + " Track deltay vs delta z: Layer " + htlayer).fill(z - zTr, y - yTr);
+ if (htlayer == 1) {
+ aida.cloud2D(hitdir + tkresid + " Layer 1 deltay vs xoca").fill(poca[0], y - yTr);
+ aida.cloud2D(hitdir + tkresid + " Layer 1 deltay vs yoca").fill(poca[1], y - yTr);
+ aida.cloud2D(hitdir + tkresid + " Layer 1 deltay vs zoca").fill(poca[2], y - yTr);
+ aida.cloud2D(hitdir + tkresid + " Layer 1 deltaz vs xoca").fill(poca[0], z - zTr);
+ aida.cloud2D(hitdir + tkresid + " Layer 1 deltaz vs yoca").fill(poca[1], z - zTr);
+ aida.cloud2D(hitdir + tkresid + " Layer 1 deltaz vs zoca").fill(poca[2], z - zTr);
+ }
+ aida.cloud2D(hitdir + tkresid + " Track vs measured y: Layer " + htlayer).fill(y, yTr);
+ aida.cloud2D(hitdir + tkresid + " Track vs measured z: Layer " + htlayer).fill(z, zTr);
+ aida.cloud2D(hitdir + tkresid + " Track deltay vs S ").fill(sHit, y - yTr);
+ aida.cloud2D(hitdir + tkresid + " Track deltaz vs S ").fill(sHit, z - zTr);
+ aida.histogram1D(hitdir + tkresid + " Track pull y: Layer " + htlayer, 200, -8, 8).fill((y - yTr) / yerr);
+ aida.histogram1D(hitdir + tkresid + " Track pull z: Layer " + htlayer, 200, -8, 8).fill((z - zTr) / zerr);
+ aida.histogram2D(hitdir + tkresid + " Track pull y vs p: Layer " + htlayer, 200, -8, 8, 200, 0, 5).fill((y - yTr) / yerr, pmc);
+ aida.histogram2D(hitdir + tkresid + " Track pull z vs p: Layer " + htlayer, 200, -8, 8, 200, 0, 5).fill((z - zTr) / zerr, pmc);
+
+
+ }
+ for (HelicalTrackStrip cl : clusterlist) {
+ int ilayer = 0;
+ List<MCParticle> mcps = cl.MCParticles();
+
+ Hep3Vector corigin = cl.origin();
+ Hep3Vector u = cl.u();
+ double umeas = cl.umeas();
+ Hep3Vector uvec = VecOp.mult(umeas, u);
+ Hep3Vector clvec = VecOp.add(corigin, uvec);
+ int layer = cl.layer();
+ HelicalTrackStrip nearest = getNearestHit(cl, toththits);
+ if (layer == 1) {
+ Double l1Dist = getNearestDistance(cl, toththits);
+ if (l1Dist != null)
+ l1Isolation.put(track, l1Dist);
+ }
+ if (nearest != null)
+ nearestHitToTrack.add(track, nearest);
+
+ int badCl = 0;
+ if (mcps.size() == 0 || mcps.size() > 1 || !mcps.contains(bestmcp))
+ badCl = 1;
+ if (badCl == 1)
+ if (mcps.size() > 0 && mcps.get(0) != null) {
+ MCParticle tmpmc = mcps.get(0);
+ aida.cloud1D(hitdir + layer + " Momentum of bad hit ").fill(tmpmc.getMomentum().magnitude());
+ aida.cloud1D(hitdir + layer + " PDGID of bad hit ").fill(tmpmc.getPDGID());
+ for (MCParticle mymc : tmpmc.getParents())
+ aida.cloud1D(hitdir + layer + " PDGID of bad hit mother ").fill(mymc.getPDGID());
+ }
+ String label = "False hit";
+ if (badCl == 0)
+ label = "True Hit ";
+
+ SimTrackerHit mcbesthit;
+ Set<SimTrackerHit> mchitlist = mcHittomcP.allTo(bestmcp);
+
+ double ymc = 0, zmc = 0;
+ for (SimTrackerHit sthbest : mchitlist) {
+ int slayer = sthbest.getLayer();
+ if (layer == slayer) {
+ mcbesthit = sthbest;
+ ymc = mcbesthit.getPoint()[1];
+ zmc = mcbesthit.getPoint()[2];
+ mcmom = getMag(mcbesthit.getMomentum());
+ if (prevmom > 0 && badCl == 0) {
+ aida.histogram1D(hitdir + layer + " MC energy difference ", 100, -0.005, 0.0).fill(mcmom - prevmom);
+ aida.histogram1D(hitdir + " MC energy difference ", 100, -0.005, 0.0).fill(mcmom - prevmom);
+ }
+ prevmom = mcmom;
+
+ }
+ }
+
+ double axdotu = VecOp.dot(cl.u(), axial);
+ boolean isAxial = false;
+ if (axdotu > 0.5)
+ isAxial = true;
+// aida.cloud2D(hitdir + layer + " y vs z " + label).fill(z, y);
+ if (isAxial) {
+ aida.cloud1D(hitdir + layer + " y " + label).fill(clvec.y());
+ aida.cloud1D(hitdir + layer + " deltay " + label).fill(clvec.y() - ymc);
+ aida.cloud2D(hitdir + layer + " y vs yMC " + label).fill(ymc, clvec.y());
+ } else {
+ aida.cloud1D(hitdir + layer + " z " + label).fill(clvec.z());
+ aida.cloud1D(hitdir + layer + " deltaz " + label).fill(clvec.z() - zmc);
+ aida.cloud2D(hitdir + layer + " z vs zMC " + label).fill(zmc, clvec.z());
+ }
+ Set<MCParticle> mclist = hittomc.allFrom(hit);
+ aida.cloud1D(hitdir + layer + " Associated MC particles").fill(mclist.size());
+ layerTuple.fill(ilayer++, (int) layer);
+ layerTuple.fill(ilayer++, (int) cl.rawhits().size());
+ layerTuple.fill(ilayer++, (int) badCl);
+ layerTuple.fill(ilayer++, (boolean) isAxial);
+ layerTuple.fill(ilayer++, (float) clvec.y());
+ layerTuple.fill(ilayer++, (float) clvec.z());
+ layerTuple.fill(ilayer++, (float) ymc);
+ layerTuple.fill(ilayer++, (float) zmc);
+ layerTuple.addRow();
+
+ }
+ }
+ tuple.addRow();
+ }
+ }
+
+ // Make the normalized fake plot after the specified number of events
+ if (nevt == 1000) {
+ double wgt = 1. / trk_count;
+ for (int i = 0; i < 10; i++) {
+ System.out.println(" Entries: " + fakes.binEntries(i) + " for mismatches: " + i);
+ for (int j = 0; j < fakes.binHeight(i); j++)
+ nfakes.fill(i, wgt);
+ }
+ System.out.println("Normalization: " + nfakes.sumAllBinHeights() + " after ntrk = " + trk_count);
+ }
+
+//analyze the axial only tracks
+ for (Track track : axialtracklist) {
+ TrackAnalysis tkanal = new TrackAnalysis(track, hittomcAxial);
+ MCParticle mcp = tkanal.getMCParticle();
+ if (mcp != null)
+ // Create a map between the tracks found and the assigned MC particle
+ trktomcAxial.add(track, tkanal.getMCParticle());
+
+ }
+ for (HelicalTrackHit hit : toththits) {
+
+ int nAssHits = hit.getRawHits().size();
+ aida.cloud1D(debugDir + hit.Detector() + " nAssHits").fill(nAssHits);
+ Hep3Vector HTHPos = hit.getCorrectedPosition();
+ double rHit = Math.sqrt(HTHPos.x() * HTHPos.x() + HTHPos.y() * HTHPos.y());
+ double zHit = HTHPos.z();
+ double etaHit = -Math.log(Math.tan(Math.atan2(rHit, zHit) / 2));
+ double hitchisq = hit.chisq();
+// double s = helixTrack.PathMap().get(hit);
+// Hep3Vector posonhelix = HelixUtils.PointOnHelix(helixTrack, s);
+
+ if (hit instanceof HelicalTrackCross) {
+ HelicalTrackCross cross = (HelicalTrackCross) hit;
+// TrackDirection trkdir = HelixUtils.CalculateTrackDirection(helixTrack, s);
+// cross.setTrackDirection(trkdir, helixTrack.covariance());
+ List<HelicalTrackStrip> clusterlist = cross.getStrips();
+ double du_stereo = 0;
+ double du_axial = 0;
+ for (HelicalTrackStrip cluster : clusterlist) {
+
+ int nstrips = cluster.rawhits().size();
+ aida.cloud1D(debugDir + hit.Detector() + " nStrips-per-layer").fill(nstrips);
+ Hep3Vector corigin = cluster.origin();
+ Hep3Vector u = cluster.u();
+ List<RawTrackerHit> rawhits = cluster.rawhits();
+ double umc = -999999;
+ double stenergy = -999999;
+ String stripdir = "axial";
+ double umeas = cluster.umeas();
+ double charge = cluster.dEdx() * 1000.0;
+ double layer = cluster.layer();
+ for (RawTrackerHit rhit : rawhits) {
+
+ String deName = rhit.getDetectorElement().getName();
+ if (deName.contains("sensor1"))
+ stripdir = "stereo";
+ // System.out.println("Layer number " + rhit.getLayerNumber() + " " + deName);
+ List<SimTrackerHit> sthits = rhit.getSimTrackerHits();
+ int nsthits = sthits.size();
+ aida.cloud1D(debugDir + hit.Detector() + " associated ST hits").fill(nsthits);
+ aida.cloud1D(debugDir + hit.Detector() + " layer" + stripdir + " associated ST hits").fill(nsthits);
+ if (nsthits == 1) {
+ double[] sthitD = sthits.get(0).getPoint();
+ BasicHep3Vector sthit = new BasicHep3Vector(sthitD);
+ stenergy = sthits.get(0).getdEdx();
+ Hep3Vector vdiff = VecOp.sub(sthit, corigin);
+ umc = VecOp.dot(vdiff, u);
+ }
+ }
+
+
+ // System.out.println("filling...");
+ if (umc != -999999) {
+// aida.cloud1D(debugDir + hit.Detector() + "dedx").fill(charge);
+ aida.histogram1D(debugDir + hit.Detector() + "dedx", 50, 0, 0.3).fill(charge);
+ if (umc < 1 && umc > -1)
+ aida.cloud2D(debugDir + hit.Detector() + "cluster reco vs cluster mc").fill(umeas - umc, umc);
+ aida.cloud2D(debugDir + hit.Detector() + "cluster vs STHit dedx").fill(stenergy, charge);
+ aida.cloud2D(debugDir + hit.Detector() + "cluster dedx vs delte(u)").fill(umeas - umc, charge);
+ if (stripdir.contains("stereo"))
+ du_stereo = umeas - umc;
+ if (stripdir.contains("axial"))
+ du_axial = umeas - umc;
+ aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u)").fill(umeas - umc);
+ if (nstrips == 1) {
+ aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)--1 strip").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u)--1 strip").fill(umeas - umc);
+ }
+ if (nstrips == 2) {
+ aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)--2 strip").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u)--2 strip").fill(umeas - umc);
+ }
+ if (nstrips == 3) {
+ aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)--3 strip").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u)--3 strip").fill(umeas - umc);
+ }
+ }
+
+ }
+ aida.cloud2D(debugDir + hit.Detector() + " delta(u) stereo v axial").fill(du_stereo, du_axial);
+ }
+ }
+
+ //analyze the event
+ int ApCand = 0;
+ List<MCParticle> mcplist = event.getMCParticles();
+ String apdir = "Aprime/";
+ Track eleID = null;
+ Track posID = null;
+ MCParticle eleMC = null;
+ MCParticle posMC = null;
+ for (Track track : tracklist) {
+
+ TrackAnalysis tkanal = tkanalMap.get(track);
+ StraightLineTrack slt = sltMap.get(track);
+ // Calculate purity and make appropriate plots
+ MCParticle mcp = tkanal.getMCParticle();
+ if (mcp == null)
+ continue;
+ if (mcp.getParents().size() == 1 && mcp.getParents().get(0).getPDGID() == 622) {
+ int nbad = tkanal.getNBadHitsNew();
+ int nhits = tkanal.getNHitsNew();
+ double purity = tkanal.getPurityNew();
+ int nbadAxial = tkanal.getNBadAxialHits();
+ int nbadZ = tkanal.getNBadZHits();
+ double px = track.getPX();
+ double py = track.getPY();
+ double pz = track.getPZ();
+ double pt = Math.sqrt(px * px + py * py);
+ double pperp = Math.sqrt(py * py + pz * pz);
+ double p = Math.sqrt(pt * pt + pz * pz);
+ double phi = Math.atan2(py, px);
+ double cth = pz / Math.sqrt(pt * pt + pz * pz);
+ double sth = pt / Math.sqrt(pt * pt + pz * pz);
+ double th = Math.atan2(pt, pz);
+
+
+
+ double d0 = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
+ double phi0 = track.getTrackParameter(HelicalTrackFit.phi0Index);
+
+ double z0sinth = sth * z0;
+
+
+ double mom[] = track.getMomentum();
+
+
+ double doca = slt.Doca();
+ double xoca = slt.xPoca();
+ double[] poca = slt.Poca();
+ double y0new = slt.TargetYZ()[0];
+ double z0new = slt.TargetYZ()[1];
+
+ Hep3Vector Pmc = mcp.getMomentum();
+ double pxmc = Pmc.x();
+ double pymc = Pmc.y();
+ double pzmc = Pmc.z();
+ double ptmc = Math.sqrt(pxmc * pxmc + pymc * pymc);
+ double pmc = Math.sqrt(ptmc * ptmc + pzmc * pzmc);
+
+
+
+ SeedTrack st = (SeedTrack) track;
+ SeedCandidate seed = st.getSeedCandidate();
+
+ HelicalTrackFit tr = seed.getHelix();
+
+ if (mcp.getCharge() > 0) {
+ posID = track;
+ posMC = mcp;
+ fillTrackInfo(apdir, "positron", track.getChi2(), nhits, p, pperp, px, py, pz, phi, cth, doca, poca[0], poca[1], poca[2]);
+
+ } else {
+ eleID = track;
+ eleMC = mcp;
+ fillTrackInfo(apdir, "electron", track.getChi2(), nhits, p, pperp, px, py, pz, phi, cth, doca, poca[0], poca[1], poca[2]);
+ }
+ }
+
+ }
+ String vertex = "Vertexing/";
+ String selected = "Selection/";
[truncated at 1000 lines; 740 more skipped]
