Author: [log in to unmask]
Date: Wed Nov 30 11:12:55 2016
New Revision: 4586
Log:
adding analysis code for layer0.
Added:
java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/
java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/FindableTrackLayer0.java
java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Acceptance.java
java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Driver.java
java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Plots.java
java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0ReconParticle.java
java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/TrackingReconstructionPlots.java
Added: java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/FindableTrackLayer0.java
=============================================================================
--- java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/FindableTrackLayer0.java (added)
+++ java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/FindableTrackLayer0.java Wed Nov 30 11:12:55 2016
@@ -0,0 +1,380 @@
+/*
+ * FindableTrack.java
+ *
+ * Created on October 24, 2008, 9:50 PM
+ *
+ */
+package org.hps.users.mrsolt;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+
+import java.util.ArrayList;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Set;
+
+import org.lcsim.detector.DetectorElementStore;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IDetectorElementContainer;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
+import org.lcsim.fit.helicaltrack.HitIdentifier;
+import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
+import org.lcsim.recon.tracking.seedtracker.SeedLayer;
+import org.lcsim.recon.tracking.seedtracker.SeedLayer.SeedType;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+
+/**
+ *
+ * @author Richard Partridge
+ * @version $Id: FindableTrack.java,v 1.4 2012/11/08 01:22:41 omoreno Exp $
+ */
+// FIXME: This class is used in some analyses but is not actually created by the recon.
+// Should it be put someplace else like in the analysis module? --JM
+public class FindableTrackLayer0 {
+
+ public enum Ignore {
+ NoPTCut, NoDCACut, NoZ0Cut, NoSeedCheck, NoConfirmCheck, NoMinHitCut
+ };
+
+ private double _bfield;
+ private RelationalTable<SimTrackerHit, MCParticle> _hittomc;
+ private HitIdentifier _ID;
+ private int _nlayersTot = 10;
+
+ public FindableTrackLayer0(EventHeader event, List<SimTrackerHit> simTrackerHits) {
+
+ // Get the magnetic field
+ Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
+ _bfield = event.getDetector().getFieldMap().getField(IP).y();
+
+ // Instantiate the hit identifier class
+ _ID = new HitIdentifier();
+
+ // Create a relational table that maps SimTrackerHits to MCParticles
+ _hittomc = new BaseRelationalTable<SimTrackerHit, MCParticle>(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+
+ List<List<SimTrackerHit>> simTrackerHitCollections = new ArrayList<List<SimTrackerHit>>();
+
+ // If the collection of SimTrackerHits is not specified get the collection from the event.
+ // Otherwise, add the collection to the list of collections to be processed.
+ if (simTrackerHits == null)
+ simTrackerHitCollections.addAll(event.get(SimTrackerHit.class));
+ else
+ simTrackerHitCollections.add(simTrackerHits);
+
+ // Loop over the SimTrackerHits and fill in the relational table
+ for (List<SimTrackerHit> simlist : simTrackerHitCollections) {
+ for (SimTrackerHit simhit : simlist)
+ if (simhit.getMCParticle() != null)
+ _hittomc.add(simhit, simhit.getMCParticle());
+ }
+ }
+
+ public FindableTrackLayer0(EventHeader event, List<SimTrackerHit> simTrackerHits, int nLayersTot) {
+ this(event, simTrackerHits);
+ this._nlayersTot = nLayersTot;
+ }
+
+ public FindableTrackLayer0(EventHeader event) {
+ this(event, null);
+ }
+
+ public FindableTrackLayer0(EventHeader event, int nLayersTot) {
+ this(event, null, nLayersTot);
+ }
+
+ public boolean isFindable(MCParticle mcp, List<SeedStrategy> slist, Ignore ignore) {
+ List<Ignore> ignores = new ArrayList<Ignore>();
+ ignores.add(ignore);
+ return isFindable(mcp, slist, ignores);
+ }
+
+ public boolean isFindable(MCParticle mcp, List<SeedStrategy> slist) {
+ return isFindable(mcp, slist, new ArrayList<Ignore>());
+ }
+
+ public boolean isFindable(MCParticle mcp, List<SeedStrategy> slist, List<Ignore> ignores) {
+
+ // We can't find neutral particles'
+ if (mcp.getCharge() == 0)
+ return false;
+
+ // Find the helix parameters in the L3 convention used by org.lcsim
+ HelixParamCalculator helix = new HelixParamCalculator(mcp, _bfield);
+
+ // We haven't yet determined the track is findable
+ boolean findable = false;
+
+ // Loop over strategies and check if the track is findable
+ for (SeedStrategy strat : slist) {
+
+ // Check the MC Particle's pT
+ if (!CheckPT(helix, ignores, strat))
+ continue;
+
+ // Check the MC Particle's DCA
+ if (!CheckDCA(helix, ignores, strat))
+ continue;
+
+ // Check the MC Particle's Z0
+ if (!CheckZ0(helix, ignores, strat))
+ continue;
+
+ // Check that we have hits on the seed layers
+ if (!CheckSeed(mcp, ignores, strat))
+ continue;
+
+ // Check that we have the required confirmation hits
+ if (!CheckConfirm(mcp, ignores, strat))
+ continue;
+
+ // Check for the minimum number of hits
+ if (!CheckMinHits(mcp, ignores, strat))
+ continue;
+
+ // Passed all the checks - track is findable
+ findable = true;
+ break;
+ }
+
+ return findable;
+ }
+
+ public int LayersHit(MCParticle mcp) {
+
+ // Get the list of hits associated with the MCParticle
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+
+ // Create a set of the identifiers for the hit layers
+ Set<String> idset = new HashSet<String>();
+
+ // Create the set of identifiers
+ for (SimTrackerHit simhit : hitlist) {
+
+ String identifier_old = _ID.Identifier(getDetectorElement(simhit));
+ String identifier = _ID.Identifier(simhit);
+ if (!idset.contains(identifier))
+ idset.add(identifier);
+ }
+
+ return idset.size();
+ }
+
+ public boolean isTrackFindable(MCParticle mcParticle, int nLayers) {
+
+ if (nLayers % 2 == 1)
+ throw new RuntimeException(this.getClass().getSimpleName() + ": The required number of layers hit must be even");
+
+ // A neutral particle can't be found
+ if (mcParticle.getCharge() == 0)
+ return false;
+
+ // Get the list of SimTrackerHits associated with the MC particle
+ Set<SimTrackerHit> simHits = _hittomc.allTo(mcParticle);
+
+ // Find the layers hit
+ boolean[] layerHit = new boolean[_nlayersTot];
+ for (SimTrackerHit simHit : simHits) {
+ layerHit[simHit.getLayer() - 1] = true;
+ }
+
+ int nLayersHit = 0;
+ // Check how many pairs of layers were hit
+ for (int index = 0; index < _nlayersTot; index += 2) {
+ if (layerHit[index] && layerHit[index + 1])
+ nLayersHit += 2;
+ }
+
+ return nLayersHit >= nLayers;
+ }
+
+ public Set<SimTrackerHit> getSimTrackerHits(MCParticle mcParticle) {
+ return _hittomc.allTo(mcParticle);
+ }
+
+ public boolean InnerTrackerIsFindable(MCParticle mcp, int nlayers, boolean printout) {
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+ boolean[] layerHit = { false, false, false, false, false, false, false, false, false, false, false, false };
+ for (SimTrackerHit simhit : hitlist) {
+ layerHit[simhit.getLayer() - 1] = true;
+ }
+ for (int i = 0; i < nlayers; i++) {
+ System.out.println(layerHit[i]);
+ if (layerHit[i] == false)
+ return false;
+ }
+ return true;
+ }
+
+ public boolean TrackIsFindable(MCParticle mcp, int nlayers) {
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+ boolean[] layerHit = { false, false, false, false, false, false, false, false, false, false, false, false, false, false };
+ for (SimTrackerHit simhit : hitlist) {
+ layerHit[simhit.getLayer() - 1] = true;
+ }
+ for (int i = 0; i < nlayers; i++) {
+ if (layerHit[i] == false)
+ return false;
+ }
+ return true;
+ }
+
+ public boolean TrackIsFindableRequireLayer0(MCParticle mcp, int nlayers) {
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+ boolean[] layerHit = { false, false, false, false, false, false, false, false, false, false, false, false, false, false };
+ for (SimTrackerHit simhit : hitlist) {
+ layerHit[simhit.getLayer() - 1] = true;
+ }
+ for (int i = 0; i < nlayers; i++) {
+ if (layerHit[i] == false)
+ return false;
+ }
+ return true;
+ }
+
+ public boolean TrackIsFindableExcludeLayer0(MCParticle mcp, int nlayers) {
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+ boolean[] layerHit = { false, false, false, false, false, false, false, false, false, false, false, false };
+ for (SimTrackerHit simhit : hitlist) {
+ layerHit[simhit.getLayer() - 1] = true;
+ }
+ for (int i = 0; i < nlayers; i++) {
+ if (layerHit[i] == false)
+ return false;
+ }
+ return true;
+ }
+
+ public boolean OuterTrackerIsFindable(MCParticle mcp, int start) {
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+ boolean[] layerHit = { false, false, false, false, false, false, false, false, false, false, false, false };
+ for (SimTrackerHit simhit : hitlist) {
+ layerHit[simhit.getLayer() - 1] = true;
+ }
+ for (int i = start; i < _nlayersTot; i++) {
+ if (layerHit[i] == false)
+ return false;
+ }
+ return true;
+ }
+
+ private boolean CheckPT(HelixParamCalculator helix, List<Ignore> ignores, SeedStrategy strat) {
+
+ // First see if we are skipping this check
+ if (ignores.contains(Ignore.NoPTCut))
+ return true;
+
+ return helix.getMCTransverseMomentum() >= strat.getMinPT();
+ }
+
+ private boolean CheckDCA(HelixParamCalculator helix, List<Ignore> ignores, SeedStrategy strat) {
+
+ // First see if we are skipping this check
+ if (ignores.contains(Ignore.NoDCACut))
+ return true;
+
+ return Math.abs(helix.getDCA()) <= strat.getMaxDCA();
+ }
+
+ private boolean CheckZ0(HelixParamCalculator helix, List<Ignore> ignores, SeedStrategy strat) {
+
+ // First see if we are skipping this check
+ if (ignores.contains(Ignore.NoZ0Cut))
+ return true;
+
+ return Math.abs(helix.getZ0()) <= strat.getMaxZ0();
+ }
+
+ private boolean CheckSeed(MCParticle mcp, List<Ignore> ignores, SeedStrategy strat) {
+
+ // First see if we are skipping this check
+ if (ignores.contains(Ignore.NoSeedCheck))
+ return true;
+
+ return HitCount(mcp, strat.getLayers(SeedType.Seed)) == 3;
+ }
+
+ private boolean CheckConfirm(MCParticle mcp, List<Ignore> ignores, SeedStrategy strat) {
+
+ // First see if we are skipping this check
+ if (ignores.contains(Ignore.NoConfirmCheck))
+ return true;
+
+ return HitCount(mcp, strat.getLayers(SeedType.Confirm)) >= strat.getMinConfirm();
+ }
+
+ private boolean CheckMinHits(MCParticle mcp, List<Ignore> ignores, SeedStrategy strat) {
+
+ // First see if we are skipping this check
+ if (ignores.contains(Ignore.NoMinHitCut))
+ return true;
+
+ return HitCount(mcp, strat.getLayerList()) >= strat.getMinHits();
+ }
+
+ private int HitCount(MCParticle mcp, List<SeedLayer> lyrlist) {
+
+ // Get the list of hits associated with the MCParticle
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+
+ // Count the number of layers with hits in them
+ int hitcount = 0;
+ for (SeedLayer lyr : lyrlist)
+ // Loop over the hits for this MCParticle
+ for (SimTrackerHit simhit : hitlist) {
+
+ // Get the detector element for this hit
+ // IDetectorElement de = getDetectorElement(simhit);
+
+ // Check names
+ // String detname_old = _ID.getName(de);
+ String detname_new = simhit.getSubdetector().getName();
+ // if (!detname_old.equals(detname_new)) {
+ // System.out.println("Detector name mismatch - old: "+detname_old+
+ // " new: "+detname_new);
+ // }
+ // int layer_old = _ID.getLayer(de);
+ int layer_new = simhit.getLayer();
+ // if (layer_old != layer_new) {
+ // System.out.println("Layer number mismatch - old: "+layer_old+" new: "+layer_new);
+ // }
+ // BarrelEndcapFlag be_old = _ID.getBarrelEndcapFlag(de);
+ BarrelEndcapFlag be_new = simhit.getBarrelEndcapFlag();
+ // if (!be_old.equals(be_new)) {
+ // System.out.println("BarrelEndcapFlag mismatch - old: "+be_old+" new: "+be_new);
+ // }
+
+ // See if this hit is on the layer we are checking
+ // if (!lyr.getDetName().equals(_ID.getName(de))) continue;
+ // if (lyr.getLayer() != _ID.getLayer(de)) continue;
+ // if (!lyr.getBarrelEndcapFlag().equals(_ID.getBarrelEndcapFlag(de)))
+ // continue;
+ if (!lyr.getDetName().equals(detname_new))
+ continue;
+ if (lyr.getLayer() != layer_new)
+ continue;
+ if (!lyr.getBarrelEndcapFlag().equals(be_new))
+ continue;
+ hitcount++;
+ break;
+ }
+
+ return hitcount;
+ }
+
+ private IDetectorElement getDetectorElement(SimTrackerHit hit) {
+ IDetectorElementContainer cont = DetectorElementStore.getInstance().find(hit.getIdentifier());
+ IDetectorElement de;
+ if (cont.isEmpty())
+ throw new RuntimeException("Detector Container is empty!");
+ else
+ de = cont.get(0);
+ return de;
+ }
+}
Added: java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Acceptance.java
=============================================================================
--- java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Acceptance.java (added)
+++ java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Acceptance.java Wed Nov 30 11:12:55 2016
@@ -0,0 +1,1145 @@
+/**
+ * Analysis driver to calculate hit efficiencies in the SVT
+ */
+/**
+ * @author mrsolt
+ *
+ */
+
+package org.hps.users.mrsolt;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogramFactory;
+import hep.aida.IPlotterFactory;
+import hep.aida.IPlotter;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.aida.ITree;
+import hep.physics.vec.BasicHep3Matrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import org.lcsim.detector.converter.compact.subdetector.SvtStereoLayer;
+import org.lcsim.detector.tracker.silicon.HpsSiSensor;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.Vertex;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.hps.recon.tracking.TrackType;
+import org.hps.recon.tracking.TrackUtils;
+import org.hps.recon.tracking.TrackerHitUtils;
+import org.hps.recon.vertexing.BilliorTrack;
+import org.hps.recon.vertexing.BilliorVertex;
+import org.hps.recon.vertexing.BilliorVertexer;
+import org.hps.users.mrsolt.FindableTrackLayer0;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HitIdentifier;
+
+
+public class Layer0Acceptance extends Driver {
+
+
+ // Use JFreeChart as the default plotting backend
+ static {
+ hep.aida.jfree.AnalysisFactory.register();
+ }
+
+ // Plotting
+ protected AIDA aida = AIDA.defaultInstance();
+ ITree tree;
+ IHistogramFactory histogramFactory;
+ IPlotterFactory plotterFactory = IAnalysisFactory.create().createPlotterFactory();
+ protected Map<String, IPlotter> plotters = new HashMap<String, IPlotter>();
+ private Map<Integer, List<SvtStereoLayer>> topStereoLayers = new HashMap<Integer, List<SvtStereoLayer>>();
+ private Map<Integer, List<SvtStereoLayer>> bottomStereoLayers = new HashMap<Integer, List<SvtStereoLayer>>();
+ private static final String SUBDETECTOR_NAME = "Tracker";
+ RelationalTable<SimTrackerHit, MCParticle> _hittomc = null;
+ private List<HpsSiSensor> sensors = null;
+ private double _bfield;
+ private HitIdentifier _ID;
+ private int _nlayersTot = 14;
+ double eBeam = 1.05;
+ double maxFactor = 1.25;
+
+
+ /*IHistogram1D Findable4Hit;
+ IHistogram1D Findable4HitLayer0;
+ IHistogram1D Findable4HitExcludeLayer0;
+ IHistogram1D Findable5Hit;
+ IHistogram1D Findable5HitLayer0;
+ IHistogram1D Findable5HitExcludeLayer0;
+ IHistogram1D Findable6Hit;
+ IHistogram1D Findable6HitLayer0;
+ IHistogram1D Findable6HitExcludeLayer0;
+ IHistogram1D Findable7Hit;*/
+
+ int Findable4Hit = 0;
+ int Findable4HitLayer0 = 0;
+ int Findable4HitExcludeLayer0 = 0;
+ int Findable5Hit = 0;
+ int Findable5HitLayer0 = 0;
+ int Findable5HitExcludeLayer0 = 0;
+ int Findable6Hit = 0;
+ int Findable6HitLayer0 = 0;
+ int Findable6HitExcludeLayer0 = 0;
+ int Findable7Hit = 0;
+
+ int Findable4HitRecoil = 0;
+ int Findable4HitLayer0Recoil = 0;
+ int Findable4HitExcludeLayer0Recoil = 0;
+ int Findable5HitRecoil = 0;
+ int Findable5HitLayer0Recoil = 0;
+ int Findable5HitExcludeLayer0Recoil = 0;
+ int Findable6HitRecoil = 0;
+ int Findable6HitLayer0Recoil = 0;
+ int Findable6HitExcludeLayer0Recoil = 0;
+ int Findable7HitRecoil = 0;
+
+ int TridentsFindable4Hit = 0;
+ int TridentsFindable4HitLayer0 = 0;
+ int TridentsFindable4HitExcludeLayer0 = 0;
+ int TridentsFindable5Hit = 0;
+ int TridentsFindable5HitLayer0 = 0;
+ int TridentsFindable5HitExcludeLayer0 = 0;
+ int TridentsFindable6Hit = 0;
+ int TridentsFindable6HitLayer0 = 0;
+ int TridentsFindable6HitExcludeLayer0 = 0;
+ int TridentsFindable7Hit = 0;
+
+ int ThreeProngFindable4Hit = 0;
+ int ThreeProngFindable4HitLayer0 = 0;
+ int ThreeProngFindable4HitExcludeLayer0 = 0;
+ int ThreeProngFindable5Hit = 0;
+ int ThreeProngFindable5HitLayer0 = 0;
+ int ThreeProngFindable5HitExcludeLayer0 = 0;
+ int ThreeProngFindable6Hit = 0;
+ int ThreeProngFindable6HitLayer0 = 0;
+ int ThreeProngFindable6HitExcludeLayer0 = 0;
+ int ThreeProngFindable7Hit = 0;
+
+ int TracksFound4Hit = 0;
+ int TracksFound4HitLayer0 = 0;
+ int TracksFound4HitExcludeLayer0 = 0;
+ int TracksFound5Hit = 0;
+ int TracksFound5HitLayer0 = 0;
+ int TracksFound5HitExcludeLayer0 = 0;
+ int TracksFound6Hit = 0;
+ int TracksFound6HitLayer0 = 0;
+ int TracksFound6HitExcludeLayer0 = 0;
+ int TracksFound7Hit = 0;
+
+ int events = 0;
+ int goodTrident = 0;
+ int goodTridentInL0 = 0;
+ int goodTridentExL0 = 0;
+
+ IHistogram1D Findable;
+ IHistogram1D FindableLayer0;
+ IHistogram1D FindableExcludeLayer0;
+
+ IHistogram1D FindableRecoil;
+ IHistogram1D FindableLayer0Recoil;
+ IHistogram1D FindableExcludeLayer0Recoil;
+
+ IHistogram1D TridentFindable;
+ IHistogram1D TridentFindableLayer0;
+ IHistogram1D TridentFindableExcludeLayer0;
+
+ IHistogram1D ThreeProngFindable;
+ IHistogram1D ThreeProngFindableLayer0;
+ IHistogram1D ThreeProngFindableExcludeLayer0;
+
+ IHistogram1D TracksFound;
+ IHistogram1D TracksFoundLayer0;
+ IHistogram1D TracksFoundExcludeLayer0;
+
+ IHistogram1D NumberofEvents;
+
+ IHistogram2D TrackHitTopAxial;
+ IHistogram2D TrackHitBotAxial;
+ IHistogram2D TrackHitTopStereo;
+ IHistogram2D TrackHitBotStereo;
+
+ IHistogram1D TrackXHitTopAxial;
+ IHistogram1D TrackXHitBotAxial;
+ IHistogram1D TrackXHitTopStereo;
+ IHistogram1D TrackXHitBotStereo;
+
+ IHistogram2D TrackHitEleTopAxial;
+ IHistogram2D TrackHitEleBotAxial;
+ IHistogram2D TrackHitEleTopStereo;
+ IHistogram2D TrackHitEleBotStereo;
+
+ IHistogram1D TrackXHitEleTopAxial;
+ IHistogram1D TrackXHitEleBotAxial;
+ IHistogram1D TrackXHitEleTopStereo;
+ IHistogram1D TrackXHitEleBotStereo;
+
+ IHistogram2D TrackHitPosTopAxial;
+ IHistogram2D TrackHitPosBotAxial;
+ IHistogram2D TrackHitPosTopStereo;
+ IHistogram2D TrackHitPosBotStereo;
+
+ IHistogram1D TrackXHitPosTopAxial;
+ IHistogram1D TrackXHitPosBotAxial;
+ IHistogram1D TrackXHitPosTopStereo;
+ IHistogram1D TrackXHitPosBotStereo;
+
+ IHistogram1D TrackYHitTopAxial;
+ IHistogram1D TrackYHitBotAxial;
+ IHistogram1D TrackYHitTopStereo;
+ IHistogram1D TrackYHitBotStereo;
+
+ IHistogram1D TrackZHitTopAxialHole;
+ IHistogram1D TrackZHitBotAxialHole;
+ IHistogram1D TrackZHitTopStereoHole;
+ IHistogram1D TrackZHitBotStereoHole;
+ IHistogram1D TrackZHitTopAxialSlot;
+ IHistogram1D TrackZHitBotAxialSlot;
+ IHistogram1D TrackZHitTopStereoSlot;
+ IHistogram1D TrackZHitBotStereoSlot;
+
+ IHistogram1D NumberofTracks;
+ IHistogram1D NumberofHits;
+
+ IHistogram1D Chi2;
+ IHistogram1D TanLambda;
+ IHistogram1D D0;
+ IHistogram1D Z0;
+ IHistogram1D Phi;
+ IHistogram1D P;
+ IHistogram1D Px;
+ IHistogram1D Py;
+ IHistogram1D Pz;
+ IHistogram1D Pt;
+
+
+ private String stereoHitCollectionName = "HelicalTrackHits";
+ private String trackCollectionName = "MatchedTracks";
+ private String MCparticleCollectionName = "MCParticle";
+
+ int num_lay = 7;
+
+ //double ebeam = 1.056;
+ double zPosTop = 59.3105;
+ double zPosBot = 40.843;
+
+
+ public void detectorChanged(Detector detector){
+
+ aida.tree().cd("/");
+ tree = aida.tree();
+ histogramFactory = IAnalysisFactory.create().createHistogramFactory(tree);
+
+ // Get the HpsSiSensor objects from the tracker detector element
+ sensors = detector.getSubdetector(SUBDETECTOR_NAME)
+ .getDetectorElement().findDescendants(HpsSiSensor.class);
+
+ // If the detector element had no sensors associated with it, throw
+ // an exception
+ if (sensors.size() == 0) {
+ throw new RuntimeException("No sensors were found in this detector.");
+ }
+
+ /*Findable4Hit = aida.histogram1D("Findable Tracks 4 Hits", num_lay, 0, num_lay);
+ Findable4HitLayer0 = aida.histogram1D("Findable Tracks 4 Hits Layer 0 Included", num_lay, 0, num_lay);
+ Findable4HitExcludeLayer0 = aida.histogram1D("Findable Tracks 4 Hits Layer 0 Excluded", num_lay, 0, num_lay);
+ Findable5Hit = aida.histogram1D("Findable Tracks 5 Hits", num_lay, 0, num_lay);
+ Findable5HitLayer0 = aida.histogram1D("Findable Tracks 5 Hits Layer 0 Included", num_lay, 0, num_lay);
+ Findable5HitExcludeLayer0 = aida.histogram1D("Findable Tracks 5 Hits Layer 0 Excluded", num_lay, 0, num_lay);
+ Findable6Hit = aida.histogram1D("Findable Tracks 6 Hits", num_lay, 0, num_lay);
+ Findable6HitLayer0 = aida.histogram1D("Findable Tracks 6 Hits Layer 0 Included", num_lay, 0, num_lay);
+ Findable6HitExcludeLayer0 = aida.histogram1D("Findable Tracks 6 Hits Layer 0 Excluded", num_lay, 0, num_lay);
+ Findable7Hit = aida.histogram1D("Findable Tracks 7 Hits", num_lay, 0, num_lay);*/
+
+ Findable = aida.histogram1D("Findable Tracks as Function of Hits", num_lay, 1, num_lay+1);
+ FindableLayer0 = aida.histogram1D("Findable Tracks as Function of Hits with Layer 0", num_lay, 1, num_lay+1);
+ FindableExcludeLayer0 = aida.histogram1D("Findable Tracks as Function of Hits Exclude Layer0", num_lay, 1, num_lay+1);
+
+ FindableRecoil = aida.histogram1D("Findable Recoils as Function of Hits", num_lay, 1, num_lay+1);
+ FindableLayer0Recoil = aida.histogram1D("Findable Recoils as Function of Hits with Layer 0", num_lay, 1, num_lay+1);
+ FindableExcludeLayer0Recoil = aida.histogram1D("Findable Recoils as Function of Hits Exclude Layer0", num_lay, 1, num_lay+1);
+
+ TridentFindable = aida.histogram1D("Findable Tridents as Function of Hits", num_lay, 1, num_lay+1);
+ TridentFindableLayer0 = aida.histogram1D("Findable Tridents as Function of Hits with Layer 0", num_lay, 1, num_lay+1);
+ TridentFindableExcludeLayer0 = aida.histogram1D("Findable Tridents as Function of Hits Exclude Layer0", num_lay, 1, num_lay+1);
+
+ ThreeProngFindable = aida.histogram1D("Findable Three Prong Events as Function of Hits", num_lay, 1, num_lay+1);
+ ThreeProngFindableLayer0 = aida.histogram1D("Findable Three Prong Events as Function of Hits with Layer 0", num_lay, 1, num_lay+1);
+ ThreeProngFindableExcludeLayer0 = aida.histogram1D("Findable Three Prong Events as Function of Hits Exclude Layer0", num_lay, 1, num_lay+1);
+
+ TracksFound = aida.histogram1D("Found Tracks as Function of Hits", num_lay, 1, num_lay+1);
+ TracksFoundLayer0 = aida.histogram1D("Found Tracks as Function of Hits with Layer 0", num_lay, 1, num_lay+1);
+ TracksFoundExcludeLayer0 = aida.histogram1D("Found Tracks as Function of Hits Exclude Layer0", num_lay, 1, num_lay+1);
+
+ NumberofEvents = aida.histogram1D("Number of Events", 2, 0, 2);
+
+ TrackHitTopAxial = aida.histogram2D("Particle Hit Layer 0 Top Axial",120,-15,15,80,0,20);
+ TrackHitBotAxial = aida.histogram2D("Particle Hit Layer 0 Bot Axial",120,-15,15,80,-20,0);
+ TrackHitTopStereo = aida.histogram2D("Particle Hit Layer 0 Top Stereo",120,-15,15,80,0,20);
+ TrackHitBotStereo = aida.histogram2D("Particle Hit Layer 0 Bot Stereo",120,-15,15,80,-20,0);
+
+ TrackXHitTopAxial = aida.histogram1D("Particle X Hit Layer 0 Top Axial",120,-15,15);
+ TrackXHitBotAxial = aida.histogram1D("Particle X Hit Layer 0 Bot Axial",120,-15,15);
+ TrackXHitTopStereo = aida.histogram1D("Particle X Hit Layer 0 Top Stereo",120,-15,15);
+ TrackXHitBotStereo = aida.histogram1D("Particle X Hit Layer 0 Bot Stereo",120,-15,15);
+
+ TrackHitEleTopAxial = aida.histogram2D("Particle Hit Layer 0 Top Axial Electrons",120,-15,15,80,0,20);
+ TrackHitEleBotAxial = aida.histogram2D("Particle Hit Layer 0 Bot Axial Electrons",120,-15,15,80,-20,0);
+ TrackHitEleTopStereo = aida.histogram2D("Particle Hit Layer 0 Top Stereo Electrons",120,-15,15,80,0,20);
+ TrackHitEleBotStereo = aida.histogram2D("Particle Hit Layer 0 Bot Stereo Electrons",120,-15,15,80,-20,0);
+
+ TrackXHitEleTopAxial = aida.histogram1D("Particle X Hit Layer 0 Top Axial Electrons",120,-15,15);
+ TrackXHitEleBotAxial = aida.histogram1D("Particle X Hit Layer 0 Bot Axial Electrons",120,-15,15);
+ TrackXHitEleTopStereo = aida.histogram1D("Particle X Hit Layer 0 Top Stereo Electrons",120,-15,15);
+ TrackXHitEleBotStereo = aida.histogram1D("Particle X Hit Layer 0 Bot Stereo Electrons",120,-15,15);
+
+ TrackHitPosTopAxial = aida.histogram2D("Particle Hit Layer 0 Top Axial Positrons",120,-15,15,80,0,20);
+ TrackHitPosBotAxial = aida.histogram2D("Particle Hit Layer 0 Bot Axial Positrons",120,-15,15,80,-20,0);
+ TrackHitPosTopStereo = aida.histogram2D("Particle Hit Layer 0 Top Stereo Positrons",120,-15,15,80,0,20);
+ TrackHitPosBotStereo = aida.histogram2D("Particle Hit Layer 0 Bot Stereo Positrons",120,-15,15,80,-20,0);
+
+ TrackXHitPosTopAxial = aida.histogram1D("Particle X Hit Layer 0 Top Axial Positrons",120,-15,15);
+ TrackXHitPosBotAxial = aida.histogram1D("Particle X Hit Layer 0 Bot Axial Positrons",120,-15,15);
+ TrackXHitPosTopStereo = aida.histogram1D("Particle X Hit Layer 0 Top Stereo Positrons",120,-15,15);
+ TrackXHitPosBotStereo = aida.histogram1D("Particle X Hit Layer 0 Bot Stereo Positrons",120,-15,15);
+
+ TrackYHitTopAxial = aida.histogram1D("Particle Y Hit Layer 0 Top Axial",80,0,20);
+ TrackYHitBotAxial = aida.histogram1D("Particle Y Hit Layer 0 Bot Axial",80,-20,0);
+ TrackYHitTopStereo = aida.histogram1D("Particle Y Hit Layer 0 Top Stereo",80,0,20);
+ TrackYHitBotStereo = aida.histogram1D("Particle Y Hit Layer 0 Bot Stereo",80,-20,0);
+
+ TrackZHitTopAxialHole = aida.histogram1D("Particle Z Hit Layer 0 Top Axial Hole",400,0,100);
+ TrackZHitBotAxialHole = aida.histogram1D("Particle Z Hit Layer 0 Bot Axial Hole",400,0,100);
+ TrackZHitTopStereoHole = aida.histogram1D("Particle Z Hit Layer 0 Top Stereo Hole",400,0,100);
+ TrackZHitBotStereoHole = aida.histogram1D("Particle Z Hit Layer 0 Bot Stereo Hole",400,0,100);
+
+ TrackZHitTopAxialSlot = aida.histogram1D("Particle Z Hit Layer 0 Top Axial Slot",400,0,100);
+ TrackZHitBotAxialSlot = aida.histogram1D("Particle Z Hit Layer 0 Bot Axial Slot",400,0,100);
+ TrackZHitTopStereoSlot = aida.histogram1D("Particle Z Hit Layer 0 Top Stereo Slot",400,0,100);
+ TrackZHitBotStereoSlot = aida.histogram1D("Particle Z Hit Layer 0 Bot Stereo Slot",400,0,100);
+
+ NumberofTracks = aida.histogram1D("Number of Tracks Per Event",10,0,10);
+ NumberofHits = aida.histogram1D("Number of Hits Per Track",20,0,20);
+
+ Chi2 = aida.histogram1D("Chi2",50,0,100);
+ TanLambda = aida.histogram1D("Tan Lambda",50,-0.1,0.1);
+ D0 = aida.histogram1D("D0",50,-4,4);
+ Z0 = aida.histogram1D("Z0",50,-4,4);
+ Phi = aida.histogram1D("Phi",50,0,2*Math.PI);
+ P = aida.histogram1D("P",50,0,1.3*eBeam);
+ Px = aida.histogram1D("Px",50,0,1.3*eBeam);
+ Py = aida.histogram1D("Py",50,0,1.3*eBeam);
+ Pz = aida.histogram1D("Pz",50,0,1.3*eBeam);
+ Pt = aida.histogram1D("Pt",50,0,1.3*eBeam);
+
+
+ for (IPlotter plotter : plotters.values()) {
+ //plotter.show();
+ }
+ }
+
+ private void mapMCtoSimTrackerHit(List<SimTrackerHit> simTrackerHits) {
+
+ // Create a relational table that maps SimTrackerHits to MCParticles
+ _hittomc = new BaseRelationalTable<SimTrackerHit, MCParticle>(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+
+ // Loop over the SimTrackerHits and fill in the relational table
+ for (SimTrackerHit simhit : simTrackerHits) {
+ if (simhit.getMCParticle() != null)
+ _hittomc.add(simhit, simhit.getMCParticle());
+ }
+ }
+
+ @Override
+ public void process(EventHeader event){
+ //System.out.println("Pass1");
+ aida.tree().cd("/");
+
+ events++;
+
+ List<SimTrackerHit> simHits = event.get(SimTrackerHit.class, "TrackerHits");
+ this.mapMCtoSimTrackerHit(simHits);
+
+ RelationalTable hitToStrips = TrackUtils.getHitToStripsTable(event);
+ RelationalTable hitToRotated = TrackUtils.getHitToRotatedTable(event);
+
+ // If the event does not have tracks, skip it
+ //if(!event.hasCollection(Track.class, trackCollectionName)) return;
+
+ List<MCParticle> MCParticlelist = event.get(MCParticle.class, MCparticleCollectionName);
+ //List<Track> tracks = event.get(Track.class, trackCollectionName);
+ List<TrackerHit> stereoHits = event.get(TrackerHit.class, stereoHitCollectionName);
+ //System.out.println("Pass2");
+
+ for(SimTrackerHit hits : simHits){
+ //if(hits.getLayer() != 1 || hits.getLayer() !=2) continue;
+ double[] hitPos = hits.getPosition();
+ if(hits.getMomentum()[1] > 0 && hits.getLayer() == 1){
+ TrackHitTopAxial.fill(hitPos[0],hitPos[1]);
+ TrackXHitTopAxial.fill(hitPos[0]);
+ if(hits.getMCParticle().getPDGID() == 11){
+ TrackHitEleTopAxial.fill(hitPos[0],hitPos[1]);
+ TrackXHitEleTopAxial.fill(hitPos[0]);
+ }
+ else if(hits.getMCParticle().getPDGID() == -11){
+ TrackHitPosTopAxial.fill(hitPos[0],hitPos[1]);
+ TrackXHitPosTopAxial.fill(hitPos[0]);
+ }
+ TrackYHitTopAxial.fill(hitPos[1]);
+ if(hitPos[0] < 0){
+ TrackZHitTopAxialHole.fill(hitPos[2]);
+ }
+ else{
+ TrackZHitTopAxialSlot.fill(hitPos[2]);
+ }
+ }
+ //System.out.println("Pass3");
+ if(hits.getMomentum()[1] < 0 && hits.getLayer() == 2){
+ TrackHitBotAxial.fill(hitPos[0],hitPos[1]);
+ TrackXHitBotAxial.fill(hitPos[0]);
+ if(hits.getMCParticle().getPDGID() == 11){
+ TrackHitEleBotAxial.fill(hitPos[0],hitPos[1]);
+ TrackXHitEleBotAxial.fill(hitPos[0]);
+ }
+ else if(hits.getMCParticle().getPDGID() == -11){
+ TrackHitPosBotAxial.fill(hitPos[0],hitPos[1]);
+ TrackXHitPosBotAxial.fill(hitPos[0]);
+ }
+ TrackYHitBotAxial.fill(hitPos[1]);
+ if(hitPos[0] < 0){
+ TrackZHitBotAxialHole.fill(hitPos[2]);
+ }
+ else{
+ TrackZHitBotAxialSlot.fill(hitPos[2]);
+ }
+ }
+ //System.out.println("Pass4");
+ if(hits.getMomentum()[1] > 0 && hits.getLayer() == 2){
+ TrackHitTopStereo.fill(hitPos[0],hitPos[1]);
+ TrackXHitTopStereo.fill(hitPos[0]);
+ if(hits.getMCParticle().getPDGID() == 11){
+ TrackHitEleTopStereo.fill(hitPos[0],hitPos[1]);
+ TrackXHitEleTopStereo.fill(hitPos[0]);
+ }
+ else if(hits.getMCParticle().getPDGID() == -11){
+ TrackHitPosTopStereo.fill(hitPos[0],hitPos[1]);
+ TrackXHitPosTopStereo.fill(hitPos[0]);
+ }
+ TrackYHitTopStereo.fill(hitPos[1]);
+ if(hitPos[0] < 0){
+ TrackZHitTopStereoHole.fill(hitPos[2]);
+ }
+ else{
+ TrackZHitTopStereoSlot.fill(hitPos[2]);
+ }
+ }
+ //System.out.println("Pass5");
+ if(hits.getMomentum()[1] < 0 && hits.getLayer() == 1){
+ TrackHitBotStereo.fill(hitPos[0],hitPos[1]);
+ TrackXHitBotStereo.fill(hitPos[0]);
+ if(hits.getMCParticle().getPDGID() == 11){
+ TrackHitEleBotStereo.fill(hitPos[0],hitPos[1]);
+ TrackXHitEleBotStereo.fill(hitPos[0]);
+ }
+ else if(hits.getMCParticle().getPDGID() == -11){
+ TrackHitPosBotStereo.fill(hitPos[0],hitPos[1]);
+ TrackXHitPosBotStereo.fill(hitPos[0]);
+ }
+ TrackYHitBotStereo.fill(hitPos[1]);
+ if(hitPos[0] < 0){
+ TrackZHitBotStereoHole.fill(hitPos[2]);
+ }
+ else{
+ TrackZHitBotStereoSlot.fill(hitPos[2]);
+ }
+ }
+
+
+ }
+ //System.out.println("Pass6");
+ boolean eleFound4Hit = false;
+ boolean posFound4Hit = false;
+ boolean eleFound5Hit = false;
+ boolean posFound5Hit = false;
+ boolean eleFound6Hit = false;
+ boolean posFound6Hit = false;
+ boolean eleFound7Hit = false;
+ boolean posFound7Hit = false;
+
+ boolean eleFound4HitExL0 = false;
+ boolean posFound4HitExL0 = false;
+ boolean eleFound5HitExL0 = false;
+ boolean posFound5HitExL0 = false;
+ boolean eleFound6HitExL0 = false;
+ boolean posFound6HitExL0 = false;
+
+ boolean eleFound4HitInL0 = false;
+ boolean posFound4HitInL0 = false;
+ boolean eleFound5HitInL0 = false;
+ boolean posFound5HitInL0 = false;
+ boolean eleFound6HitInL0 = false;
+ boolean posFound6HitInL0 = false;
+
+ boolean eleTop = false;
+ boolean eleBot = false;
+ boolean posTop = false;
+ boolean posBot = false;
+
+ boolean eleTopExL0 = false;
+ boolean eleBotExL0 = false;
+ boolean posTopExL0 = false;
+ boolean posBotExL0 = false;
+
+ boolean eleTopInL0 = false;
+ boolean eleBotInL0 = false;
+ boolean posTopInL0 = false;
+ boolean posBotInL0 = false;
+
+ boolean recoilFound4Hit = false;
+ boolean recoilFound5Hit = false;
+ boolean recoilFound6Hit = false;
+ boolean recoilFound7Hit = false;
+
+ boolean recoilFound4HitExL0 = false;
+ boolean recoilFound5HitExL0 = false;
+ boolean recoilFound6HitExL0 = false;
+
+ boolean recoilFound4HitInL0 = false;
+ boolean recoilFound5HitInL0 = false;
+ boolean recoilFound6HitInL0 = false;
+
+ for(MCParticle MCParticle : MCParticlelist){
+ List<MCParticle> Parents = MCParticle.getParents();
+ boolean trident = false;
+ boolean recoil = false;
+ for(MCParticle parent : Parents){
+ if(parent.getPDGID() == 622)
+ trident = true;
+ if(parent.getPDGID() == -623 && parent.getPDGID() != 622 && MCParticle.getPDGID() == 11)
+ recoil = true;
+ }
+ if(recoil){
+ if(TrackIsFindable(MCParticle,4)){
+ Findable4HitRecoil++;
+ recoilFound4Hit = true;
+ }
+ if(TrackIsFindable(MCParticle,5)){
+ Findable5HitRecoil++;
+ recoilFound5Hit = true;
+ }
+ if(TrackIsFindable(MCParticle,6)){
+ Findable6HitRecoil++;
+ recoilFound6Hit = true;
+ }
+ if(TrackIsFindable(MCParticle,7)){
+ Findable7HitRecoil++;
+ recoilFound7Hit = true;
+ }
+ if(TrackIsFindableRequireLayer0(MCParticle,4)){
+ Findable4HitLayer0Recoil++;
+ recoilFound4HitExL0 = true;
+ }
+ if(TrackIsFindableRequireLayer0(MCParticle,5)){
+ Findable5HitLayer0Recoil++;
+ recoilFound5HitExL0 = true;
+ }
+ if(TrackIsFindableRequireLayer0(MCParticle,6)){
+ Findable6HitLayer0Recoil++;
+ recoilFound6HitExL0 = true;
+ }
+ if(TrackIsFindableExcludeLayer0(MCParticle,4)){
+ Findable4HitExcludeLayer0Recoil++;
+ recoilFound4HitInL0 = true;
+ }
+ if(TrackIsFindableExcludeLayer0(MCParticle,5)){
+ Findable5HitExcludeLayer0Recoil++;
+ recoilFound5HitInL0 = true;
+ }
+ if(TrackIsFindableExcludeLayer0(MCParticle,6)){
+ Findable6HitExcludeLayer0Recoil++;
+ recoilFound6HitInL0 = true;
+ }
+ }
+ if(!trident) continue;
+ if(TrackIsFindable(MCParticle,4)) {
+ Findable4Hit++;
+ if(MCParticle.getPDGID() == 11){
+ eleFound4Hit = true;
+ if(MCParticle.getPY() > 0){
+ eleTop = true;
+ }
+ else{
+ eleBot = true;
+ }
+ }
+ else if(MCParticle.getPDGID() == -11) {
+ posFound4Hit = true;
+ if(MCParticle.getPY() > 0){
+ posTop = true;
+ }
+ else{
+ posBot = true;
+ }
+ }
+ }
+ if(TrackIsFindable(MCParticle,5)) {
+
+ Findable5Hit++;
+ if(MCParticle.getPDGID() == 11){
+ eleFound5Hit = true;
+ if(MCParticle.getPY() > 0){
+ eleTop = true;
+ }
+ else{
+ eleBot = true;
+ }
+ }
+ else if(MCParticle.getPDGID() == -11) {
+ posFound5Hit = true;
+ if(MCParticle.getPY() > 0){
+ posTop = true;
+ }
+ else{
+ posBot = true;
+ }
+ }
+ }
+ if(TrackIsFindable(MCParticle,6)) {
+ Findable6Hit++;
+ if(MCParticle.getPDGID() == 11){
+ eleFound6Hit = true;
+ if(MCParticle.getPY() > 0){
+ eleTop = true;
+ }
+ else{
+ eleBot = true;
+ }
+ }
+ else if(MCParticle.getPDGID() == -11) {
+ posFound6Hit = true;
+ if(MCParticle.getPY() > 0){
+ posTop = true;
+ }
+ else{
+ posBot = true;
+ }
+ }
+ }
+ if(TrackIsFindable(MCParticle,7)) {
+ Findable7Hit++;
+ if(MCParticle.getPDGID() == 11){
+ eleFound7Hit = true;
+ if(MCParticle.getPY() > 0){
+ eleTop = true;
+ }
+ else{
+ eleBot = true;
+ }
+ }
+ else if(MCParticle.getPDGID() == -11) {
+ posFound7Hit = true;
+ if(MCParticle.getPY() > 0){
+ posTop = true;
+ }
+ else{
+ posBot = true;
+ }
+ }
+ }
+ //System.out.println("Pass7");
+
+ if(TrackIsFindableRequireLayer0(MCParticle,4)) {
+ Findable4HitLayer0++;
+ if(MCParticle.getPDGID() == 11){
+ eleFound4HitInL0 = true;
+ if(MCParticle.getPY() > 0){
+ eleTopInL0 = true;
+ }
+ else{
+ eleBotInL0 = true;
+ }
+ }
+ else if(MCParticle.getPDGID() == -11) {
+ posFound4HitInL0 = true;
+ if(MCParticle.getPY() > 0){
+ posTopInL0 = true;
+ }
+ else{
+ posBotInL0 = true;
+ }
+ }
+ }
+ if(TrackIsFindableRequireLayer0(MCParticle,5)) {
+ Findable5HitLayer0++;
+ if(MCParticle.getPDGID() == 11){
+ eleFound5HitInL0 = true;
+ if(MCParticle.getPY() > 0){
+ eleTopInL0 = true;
+ }
+ else{
+ eleBotInL0 = true;
+ }
+ }
+ else if(MCParticle.getPDGID() == -11) {
+ posFound5HitInL0 = true;
+ if(MCParticle.getPY() > 0){
+ posTopInL0 = true;
+ }
+ else{
+ posBotInL0 = true;
+ }
+ }
+ }
+ if(TrackIsFindableRequireLayer0(MCParticle,6)) {
+ Findable6HitLayer0++;
+ if(MCParticle.getPDGID() == 11){
+ eleFound6HitInL0 = true;
+ if(MCParticle.getPY() > 0){
+ eleTopInL0 = true;
+ }
+ else{
+ eleBotInL0 = true;
+ }
+ }
+ else if(MCParticle.getPDGID() == -11) {
+ posFound6HitInL0 = true;
+ if(MCParticle.getPY() > 0){
+ posTopInL0 = true;
+ }
+ else{
+ posBotInL0 = true;
+ }
+ }
+ }
+ //System.out.println("Pass8");
+
+ if(TrackIsFindableExcludeLayer0(MCParticle,4)) {
+ Findable4HitExcludeLayer0++;
+ if(MCParticle.getPDGID() == 11){
+ eleFound4HitExL0 = true;
+ if(MCParticle.getPY() > 0){
+ eleTopExL0 = true;
+ }
+ else{
+ eleBotExL0 = true;
+ }
+ }
+ else if(MCParticle.getPDGID() == -11) {
+ posFound4HitExL0 = true;
+ if(MCParticle.getPY() > 0){
+ posTopExL0 = true;
+ }
+ else{
+ posBotExL0 = true;
+ }
+ }
+ }
+ if(TrackIsFindableExcludeLayer0(MCParticle,5)) {
+ Findable5HitExcludeLayer0++;
+ if(MCParticle.getPDGID() == 11){
+ eleFound5HitExL0 = true;
+ if(MCParticle.getPY() > 0){
+ eleTopExL0 = true;
+ }
+ else{
+ eleBotExL0 = true;
+ }
+ }
+ else if(MCParticle.getPDGID() == -11) {
+ posFound5HitExL0 = true;
+ if(MCParticle.getPY() > 0){
+ posTopExL0 = true;
+ }
+ else{
+ posBotExL0 = true;
+ }
+ }
+ }
+ if(TrackIsFindableExcludeLayer0(MCParticle,6)) {
+ Findable6HitExcludeLayer0++;
+ if(MCParticle.getPDGID() == 11){
+ eleFound6HitExL0 = true;
+ if(MCParticle.getPY() > 0){
+ eleTopExL0 = true;
+ }
+ else{
+ eleBotExL0 = true;
+ }
+ }
+ else if(MCParticle.getPDGID() == -11) {
+ posFound6HitExL0 = true;
+ if(MCParticle.getPY() > 0){
+ posTopExL0 = true;
+ }
+ else{
+ posBotExL0 = true;
+ }
+ }
+ }
+ //System.out.println("Pass9");
+ }
+
+ if((eleTop && posBot) || (eleBot && posTop)){
+ goodTrident++;
+ if(eleFound7Hit && posFound7Hit){
+ TridentsFindable7Hit++;
+ if(recoilFound7Hit) ThreeProngFindable7Hit++;
+ }
+ if((eleFound6Hit && posFound6Hit) || (eleFound7Hit && posFound6Hit) || (eleFound6Hit && posFound7Hit)){
+ TridentsFindable6Hit++;
+ if(recoilFound6Hit || recoilFound7Hit) ThreeProngFindable6Hit++;
+ }
+ if((eleFound5Hit && posFound5Hit) || (eleFound5Hit && posFound6Hit) || (eleFound5Hit && posFound7Hit) || (eleFound6Hit && posFound5Hit) || (eleFound7Hit && posFound5Hit)){
+ TridentsFindable5Hit++;
+ if(recoilFound5Hit || recoilFound6Hit || recoilFound7Hit) ThreeProngFindable5Hit++;
+ }
+ if((eleFound4Hit && posFound4Hit) || (eleFound4Hit && posFound5Hit) || (eleFound4Hit && posFound6Hit) || (eleFound4Hit && posFound7Hit) || (eleFound5Hit && posFound4Hit) || (eleFound6Hit && posFound4Hit) || (eleFound7Hit && posFound4Hit)){
+ TridentsFindable4Hit++;
+ if(recoilFound4Hit || recoilFound5Hit || recoilFound6Hit || recoilFound7Hit) ThreeProngFindable4Hit++;
+ }
+ }
+
+ if((eleTopInL0 && posBotInL0) || (eleBotInL0 && posTopInL0)){
+ goodTridentInL0++;
+ if((eleFound6HitInL0 && posFound6HitInL0) || (eleFound7Hit && posFound6HitInL0) || (eleFound6HitInL0 && posFound7Hit)){
+ TridentsFindable6HitLayer0++;
+ if(recoilFound6HitInL0 || recoilFound7Hit) ThreeProngFindable6HitLayer0++;
+ }
+ if((eleFound5HitInL0 && posFound5HitInL0) || (eleFound5HitInL0 && posFound6HitInL0) || (eleFound5HitInL0 && posFound7Hit) || (eleFound6HitInL0 && posFound5HitInL0) || (eleFound7Hit && posFound5HitInL0)){
+ TridentsFindable5HitLayer0++;
+ if(recoilFound5HitInL0 || recoilFound6HitInL0 || recoilFound7Hit) ThreeProngFindable5HitLayer0++;
+ }
+ if((eleFound4HitInL0 && posFound4HitInL0) || (eleFound4HitInL0 && posFound5HitInL0) || (eleFound4HitInL0 && posFound6HitInL0) || (eleFound4HitInL0 && posFound7Hit) || (eleFound5HitInL0 && posFound4HitInL0) || (eleFound6HitInL0 && posFound4HitInL0) || (eleFound7Hit && posFound4HitInL0)){
+ TridentsFindable4HitLayer0++;
+ if(recoilFound4HitInL0 || recoilFound5HitInL0 || recoilFound6HitInL0 || recoilFound7Hit) ThreeProngFindable4HitLayer0++;
+ }
+ }
+
+ if((eleTopExL0 && posBotExL0) || (eleBotExL0 && posTopExL0)){
+ goodTridentExL0++;
+ if((eleFound6HitExL0 && posFound6HitExL0)){
+ TridentsFindable6HitExcludeLayer0++;
+ if(recoilFound6HitExL0) ThreeProngFindable6HitExcludeLayer0++;
+ }
+ if((eleFound5HitExL0 && posFound5HitExL0) || (eleFound5HitExL0 && posFound6HitExL0) || (eleFound6HitExL0 && posFound5HitExL0)){
+ TridentsFindable5HitExcludeLayer0++;
+ if(recoilFound5HitExL0 || recoilFound6HitExL0) ThreeProngFindable5HitExcludeLayer0++;
+ }
+ if((eleFound4HitExL0 && posFound4HitExL0) || (eleFound4HitExL0 && posFound5HitExL0) || (eleFound4HitExL0 && posFound6HitExL0) || (eleFound5HitExL0 && posFound4HitExL0) || (eleFound6HitExL0 && posFound4HitExL0)){
+ TridentsFindable4HitExcludeLayer0++;
+ if(recoilFound4HitExL0 || recoilFound5HitExL0 || recoilFound6HitExL0) ThreeProngFindable4HitExcludeLayer0++;
+ }
+ }
+
+ /*for(MCParticle MCParticle : MCParticlelist){
+ List<MCParticle> Parents = MCParticle.getParents();
+ boolean trident = false;
+
+ for(MCParticle parent : Parents){
+ if(parent.getPDGID() == 622) trident = true;
+ }
+ if(!trident) continue;
+
+ }*/
+
+ List<List<Track>> trackCollections = event.get(Track.class);
+
+ for (List<Track> tracks : trackCollections) {
+
+ for(Track track : tracks){
+
+ HpsSiSensor trackSensor = (HpsSiSensor) ((RawTrackerHit)track.getTrackerHits().get(0).getRawHits().get(0)).getDetectorElement();
+
+ double d0 = track.getTrackStates().get(0).getD0();
+ //if(Math.abs(d0) > 10) continue;
+
+ NumberofHits.fill(track.getTrackerHits().get(0).getRawHits().size());
+ //System.out.println("Get Tracker Hits " + track.getTrackerHits());
+ //System.out.println("Get Raw Hits " + track.getTrackerHits().get(0).getRawHits());
+ //System.out.println("Get Raw Hits Get 0 " + track.getTrackerHits().get(0).getRawHits().get(0));
+ NumberofTracks.fill(tracks.size());
+ int nHits = track.getTrackerHits().size();
+ NumberofHits.fill(nHits);
+ boolean Layer0Hit =false;
+ // Calculate the momentum of the track
+ HelicalTrackFit hlc_trk_fit = TrackUtils.getHTF(track.getTrackStates().get(0));
+
+
+
+ double B_field = Math.abs(TrackUtils.getBField(event.getDetector()).y());
+
+ double p = hlc_trk_fit.p(B_field);
+ double pT = hlc_trk_fit.pT(B_field);
+ //double phi = hlc_trk_fit.phi0();
+ double tanLambda = track.getTrackStates().get(0).getTanLambda();
+ double phi = track.getTrackStates().get(0).getPhi();
+
+
+ P.fill(p);
+ Pt.fill(pT);
+ Px.fill(pT * Math.cos(phi));
+ Py.fill(pT * Math.sin(phi));
+ Pz.fill(pT * tanLambda);
+
+ Chi2.fill(track.getChi2());
+ TanLambda.fill(tanLambda);
+ D0.fill(d0);
+ Z0.fill(track.getTrackStates().get(0).getZ0());
+ Phi.fill(track.getTrackStates().get(0).getPhi());
+
+ for (TrackerHit stereoHit : stereoHits) {
+ HpsSiSensor hitSensor = (HpsSiSensor) ((RawTrackerHit) stereoHit.getRawHits().get(0)).getDetectorElement();
+
+ // Retrieve the layer number by using the sensor
+ int layer = (hitSensor.getLayerNumber() + 1)/2;
+ //nHits++;
+ if (layer == 1) Layer0Hit = true;
+ //System.out.println("Layer = " + layer + " nHits = " + nHits);
+ }
+ if(nHits >= 4) TracksFound4Hit++;
+ if(nHits >= 5) TracksFound5Hit++;
+ if(nHits >= 6) TracksFound6Hit++;
+ if(nHits >= 7) TracksFound7Hit++;
+
+ if(nHits >= 4 && Layer0Hit) TracksFound4HitLayer0++;
+ if(nHits >= 5 && Layer0Hit) TracksFound5HitLayer0++;
+ if(nHits >= 6 && Layer0Hit) TracksFound6HitLayer0++;
+
+ if(nHits >= 4 && !Layer0Hit) TracksFound4HitExcludeLayer0++;
+ if(nHits >= 5 && !Layer0Hit) TracksFound5HitExcludeLayer0++;
+ if(nHits >= 6 && !Layer0Hit) TracksFound6HitExcludeLayer0++;
+
+ NumberofHits.fill(nHits);
+ }
+ }
+ }
+ @Override
+ public void endOfData(){
+ System.out.println("Findable Tracks 4 Hits: " + Findable4Hit);
+ Findable.fill(4,Findable4Hit);
+ System.out.println("Findable Tracks 5 Hits: " + Findable5Hit);
+ Findable.fill(5,Findable5Hit);
+ System.out.println("Findable Tracks 6 Hits: " + Findable6Hit);
+ Findable.fill(6,Findable6Hit);
+ System.out.println("Findable Tracks 7 Hits: " + Findable7Hit);
+ Findable.fill(7,Findable7Hit);
+
+ System.out.println("Findable Tracks 4 Hits Require Layer0: " + Findable4HitLayer0);
+ FindableLayer0.fill(4,Findable4HitLayer0);
+ System.out.println("Findable Tracks 5 Hits Require Layer0: " + Findable5HitLayer0);
+ FindableLayer0.fill(5,Findable5HitLayer0);
+ System.out.println("Findable Tracks 6 Hits Require Layer0: " + Findable6HitLayer0);
+ FindableLayer0.fill(6,Findable6HitLayer0);
+
+ System.out.println("Findable Tracks 4 Hits Exclude Layer0: " + Findable4HitExcludeLayer0);
+ FindableExcludeLayer0.fill(4,Findable4HitExcludeLayer0);
+ System.out.println("Findable Tracks 5 Hits Exclude Layer0: " + Findable5HitExcludeLayer0);
+ FindableExcludeLayer0.fill(5,Findable5HitExcludeLayer0);
+ System.out.println("Findable Tracks 6 Hits Exclude Layer0: " + Findable6HitExcludeLayer0);
+ FindableExcludeLayer0.fill(6,Findable6HitExcludeLayer0);
+
+
+ System.out.println("Findable Tridents 4 Hits: " + TridentsFindable4Hit);
+ TridentFindable.fill(4,TridentsFindable4Hit);
+ System.out.println("Findable Tridents 5 Hits: " + TridentsFindable5Hit);
+ TridentFindable.fill(5,TridentsFindable5Hit);
+ System.out.println("Findable Tridents 6 Hits: " + TridentsFindable6Hit);
+ TridentFindable.fill(6,TridentsFindable6Hit);
+ System.out.println("Findable Tridents 7 Hits: " + TridentsFindable7Hit);
+ TridentFindable.fill(7,TridentsFindable7Hit);
+
+ System.out.println("Findable Tridents 4 Hits Require Layer0: " + TridentsFindable4HitLayer0);
+ TridentFindableLayer0.fill(4,TridentsFindable4HitLayer0);
+ System.out.println("Findable Tridents 5 Hits Require Layer0: " + TridentsFindable5HitLayer0);
+ TridentFindableLayer0.fill(5,TridentsFindable5HitLayer0);
+ System.out.println("Findable Tridents 6 Hits Require Layer0: " + TridentsFindable6HitLayer0);
+ TridentFindableLayer0.fill(6,TridentsFindable6HitLayer0);
+
+ System.out.println("Findable Tridents 4 Hits Exclude Layer0: " + TridentsFindable4HitExcludeLayer0);
+ TridentFindableExcludeLayer0.fill(4,TridentsFindable4HitExcludeLayer0);
+ System.out.println("Findable Tridents 5 Hits Exclude Layer0: " + TridentsFindable5HitExcludeLayer0);
+ TridentFindableExcludeLayer0.fill(5,TridentsFindable5HitExcludeLayer0);
+ System.out.println("Findable Tridents 6 Hits Exclude Layer0: " + TridentsFindable6HitExcludeLayer0);
+ TridentFindableExcludeLayer0.fill(6,TridentsFindable6HitExcludeLayer0);
+
+ System.out.println("Findable Recoils 4 Hits: " + Findable4HitRecoil);
+ FindableRecoil.fill(4,Findable4HitRecoil);
+ System.out.println("Findable Recoils 5 Hits: " + Findable5HitRecoil);
+ FindableRecoil.fill(5,Findable5HitRecoil);
+ System.out.println("Findable Recoils 6 Hits: " + Findable6HitRecoil);
+ FindableRecoil.fill(6,Findable6HitRecoil);
+ System.out.println("Findable Recoils 7 Hits: " + Findable7HitRecoil);
+ FindableRecoil.fill(7,Findable7HitRecoil);
+
+ System.out.println("Findable Recoils 4 Hits Require Layer0: " + Findable4HitLayer0Recoil);
+ FindableLayer0Recoil.fill(4,Findable4HitLayer0Recoil);
+ System.out.println("Findable Recoils 5 Hits Require Layer0: " + Findable5HitLayer0Recoil);
+ FindableLayer0Recoil.fill(5,Findable5HitLayer0Recoil);
+ System.out.println("Findable Recoils 6 Hits Require Layer0: " + Findable6HitLayer0Recoil);
+ FindableLayer0Recoil.fill(6,Findable6HitLayer0Recoil);
+
+ System.out.println("Findable Recoils 4 Hits Exclude Layer0: " + Findable4HitExcludeLayer0Recoil);
+ FindableExcludeLayer0Recoil.fill(4,Findable4HitExcludeLayer0Recoil);
+ System.out.println("Findable Recoils 5 Hits Exclude Layer0: " + Findable5HitExcludeLayer0Recoil);
+ FindableExcludeLayer0Recoil.fill(5,Findable5HitExcludeLayer0Recoil);
+ System.out.println("Findable Recoils 6 Hits Exclude Layer0: " + Findable6HitExcludeLayer0Recoil);
+ FindableExcludeLayer0Recoil.fill(6,Findable6HitExcludeLayer0Recoil);
+
+
+ System.out.println("Findable Three Prongs 4 Hits: " + ThreeProngFindable4Hit);
+ ThreeProngFindable.fill(4,ThreeProngFindable4Hit);
+ System.out.println("Findable Three Prongs 5 Hits: " + ThreeProngFindable5Hit);
+ ThreeProngFindable.fill(5,ThreeProngFindable5Hit);
+ System.out.println("Findable Three Prongs 6 Hits: " + ThreeProngFindable6Hit);
+ ThreeProngFindable.fill(6,ThreeProngFindable6Hit);
+ System.out.println("Findable Three Prongs 7 Hits: " + ThreeProngFindable7Hit);
+ ThreeProngFindable.fill(7,ThreeProngFindable7Hit);
+
+ System.out.println("Findable Three Prongs 4 Hits Require Layer0: " + ThreeProngFindable4HitLayer0);
+ ThreeProngFindableLayer0.fill(4,ThreeProngFindable4HitLayer0);
+ System.out.println("Findable Three Prongs 5 Hits Require Layer0: " + ThreeProngFindable5HitLayer0);
+ ThreeProngFindableLayer0.fill(5,ThreeProngFindable5HitLayer0);
+ System.out.println("Findable Three Prongs 6 Hits Require Layer0: " + ThreeProngFindable6HitLayer0);
+ ThreeProngFindableLayer0.fill(6,ThreeProngFindable6HitLayer0);
+
+ System.out.println("Findable Three Prongs 4 Hits Exclude Layer0: " + ThreeProngFindable4HitExcludeLayer0);
+ ThreeProngFindableExcludeLayer0.fill(4,ThreeProngFindable4HitExcludeLayer0);
+ System.out.println("Findable Three Prongs 5 Hits Exclude Layer0: " + ThreeProngFindable5HitExcludeLayer0);
+ ThreeProngFindableExcludeLayer0.fill(5,ThreeProngFindable5HitExcludeLayer0);
+ System.out.println("Findable Three Prongs 6 Hits Exclude Layer0: " + ThreeProngFindable6HitExcludeLayer0);
+ ThreeProngFindableExcludeLayer0.fill(6,ThreeProngFindable6HitExcludeLayer0);
+
+
+
+ System.out.println("Found Tracks 4 Hits: " + TracksFound4Hit);
+ TracksFound.fill(4,TracksFound4Hit);
+ System.out.println("Found Tracks 5 Hits: " + TracksFound5Hit);
+ TracksFound.fill(5,TracksFound5Hit);
+ System.out.println("Found Tracks 6 Hits: " + TracksFound6Hit);
+ TracksFound.fill(6,TracksFound6Hit);
+ System.out.println("Found Tracks 7 Hits: " + TracksFound7Hit);
+ TracksFound.fill(7,TracksFound7Hit);
+
+ System.out.println("Found Tracks 4 Hits Require Layer0: " + TracksFound4HitLayer0);
+ TracksFoundLayer0.fill(4,TracksFound4HitLayer0);
+ System.out.println("Found Tracks 5 Hits Require Layer0: " + TracksFound5HitLayer0);
+ TracksFoundLayer0.fill(5,TracksFound5HitLayer0);
+ System.out.println("Found Tracks 6 Hits Require Layer0: " + TracksFound6HitLayer0);
+ TracksFoundLayer0.fill(6,TracksFound6HitLayer0);
+
+ System.out.println("Found Tracks 4 Hits Exclude Layer0: " + TracksFound4HitExcludeLayer0);
+ TracksFoundExcludeLayer0.fill(4,TracksFound4HitExcludeLayer0);
+ System.out.println("Found Tracks 5 Hits Exclude Layer0: " + TracksFound5HitExcludeLayer0);
+ TracksFoundExcludeLayer0.fill(5,TracksFound5HitExcludeLayer0);
+ System.out.println("Found Tracks 6 Hits Exclude Layer0: " + TracksFound6HitExcludeLayer0);
+ TracksFoundExcludeLayer0.fill(6,TracksFound6HitExcludeLayer0);
+
+ NumberofEvents.fill(1,events);
+ System.out.println("Good Tridents = " + goodTrident);
+ System.out.println("Good Tridents Include L0 = " + goodTridentInL0);
+ System.out.println("Good Tridents Exclude L0 = " + goodTridentExL0);
+
+ System.out.println("%===================================================================%");
+ System.out.println("%====================== Layer 0 Acceptance Complete==========================%");
+ System.out.println("%===================================================================% \n%");
+ System.out.println("% \n%===================================================================%");
+ }
+
+
+ public boolean TrackIsFindable(MCParticle mcp, int nlayers) {
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+ boolean[] sensorHit = { false, false, false, false, false, false, false, false, false, false, false, false, false, false };
+ boolean[] layerHit = { false, false, false, false, false, false, false };
+ int numHit = 0;
+ for (SimTrackerHit simhit : hitlist) {
+ sensorHit[simhit.getLayer() - 1] = true;
+ }
+ for(int i = 0; i < 7; i++){
+ if(sensorHit[2*i] == true && sensorHit[2*i+1] == true){
+ layerHit[i] = true;
+ }
+ }
+ for(int i = 0; i < 7; i++){
+ if(layerHit[i] == true){
+ numHit++;
+ }
+ }
+ if(numHit < nlayers){
+ return false;
+ }
+ else{
+ return true;
+ }
+ }
+
+ public boolean TrackIsFindableRequireLayer0(MCParticle mcp, int nlayers) {
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+ boolean[] sensorHit = { false, false, false, false, false, false, false, false, false, false, false, false, false, false };
+ boolean[] layerHit = { false, false, false, false, false, false, false };
+ int numHit = 0;
+ for (SimTrackerHit simhit : hitlist) {
+ sensorHit[simhit.getLayer() - 1] = true;
+ }
+ for(int i = 0; i < 7; i++){
+ if(sensorHit[2*i] == true && sensorHit[2*i+1] == true){
+ layerHit[i] = true;
+ }
+ }
+ for(int i = 0; i < 7; i++){
+ if(layerHit[i] == true){
+ numHit++;
+ }
+ }
+ if(numHit < nlayers || layerHit[0] == false){
+ return false;
+ }
+ else{
+ return true;
+ }
+ }
+
+ public boolean TrackIsFindableExcludeLayer0(MCParticle mcp, int nlayers) {
+ Set<SimTrackerHit> hitlist = _hittomc.allTo(mcp);
+ boolean[] sensorHit = { false, false, false, false, false, false, false, false, false, false, false, false, false, false };
+ boolean[] layerHit = { false, false, false, false, false, false, false };
+ int numHit = 0;
+ for (SimTrackerHit simhit : hitlist) {
+ sensorHit[simhit.getLayer() - 1] = true;
+ }
+ for(int i = 0; i < 7; i++){
+ if(sensorHit[2*i] == true && sensorHit[2*i+1] == true){
+ layerHit[i] = true;
+ }
+ }
+ for(int i = 1; i < 7; i++){
+ if(layerHit[i] == true){
+ numHit++;
+ }
+ }
+ if(numHit < nlayers){
+ return false;
+ }
+ else{
+ return true;
+ }
+ }
+
+}
Added: java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Driver.java
=============================================================================
--- java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Driver.java (added)
+++ java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Driver.java Wed Nov 30 11:12:55 2016
@@ -0,0 +1,477 @@
+package org.hps.users.mrsolt;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogramFactory;
+import hep.aida.IPlotter;
+import hep.aida.IPlotterFactory;
+import hep.aida.IPlotterStyle;
+import hep.aida.ITree;
+
+import java.io.IOException;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+import org.hps.analysis.dataquality.SvtMonitoring;
+import org.hps.recon.tracking.gbl.GBLStripClusterData;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.GenericObject;
+import org.lcsim.geometry.Detector;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+
+/**
+ * DQM driver for track residuals in position and time
+ *
+ * @author mgraham on June 5, 2014
+ */
+// TODO: Add some quantities for DQM monitoring:
+public class Layer0Driver extends Driver {
+
+ static {
+ hep.aida.jfree.AnalysisFactory.register();
+ }
+
+ private static final Logger LOGGER = Logger.getLogger(Layer0Driver.class.getPackage().getName());
+
+ // Collection Names
+ String trackTimeDataCollectionName = "TrackTimeData";
+ String trackResidualsCollectionName = "TrackResiduals";
+ String gblStripClusterDataCollectionName = "GBLStripClusterData";
+
+ int nEvents = 0;
+
+ private AIDA aida = AIDA.defaultInstance();
+ ITree tree;
+ IHistogramFactory histogramFactory;
+ IPlotterFactory plotterFactory = IAnalysisFactory.create().createPlotterFactory();
+
+ private final String plotDir = "TrackResiduals/";
+ String[] trackingQuantNames = {};
+ private final int nmodules = 7;
+ private final int nsensors = 20;
+ private final String posresDir = "PostionResiduals/";
+ private final String uresDir = "UResiduals/";
+ private final String timeresDir = "TimeResiduals/";
+ private final String outputPlotDir = "ResidualsPlots/";
+ private boolean outputPlots = true;
+ private Map<String, Double> xposTopMeanResidMap;
+ private Map<String, Double> yposTopMeanResidMap;
+ private Map<String, Double> xposBotMeanResidMap;
+ private Map<String, Double> yposBotMeanResidMap;
+ private Map<String, Double> timeMeanResidMap;
+ private Map<String, Double> xposTopSigmaResidMap;
+ private Map<String, Double> yposTopSigmaResidMap;
+ private Map<String, Double> xposBotSigmaResidMap;
+ private Map<String, Double> yposBotSigmaResidMap;
+ private Map<String, Double> timeSigmaResidMap;
+ IHistogram1D[] xresidtop = new IHistogram1D[nmodules];
+ IHistogram1D[] yresidtop = new IHistogram1D[nmodules];
+ IHistogram1D[] xresidbot = new IHistogram1D[nmodules];
+ IHistogram1D[] yresidbot = new IHistogram1D[nmodules];
+ IHistogram1D[] utopresid = new IHistogram1D[nsensors];
+ IHistogram1D[] ubotresid = new IHistogram1D[nsensors];
+ IHistogram1D[] tresid = new IHistogram1D[nmodules * 2];
+ IHistogram1D xtopresidBS;
+ IHistogram1D ytopresidBS;
+ IHistogram1D xbotresidBS;
+ IHistogram1D ybotresidBS;
+
+// IHistogram2D[] utopresidVsU = new IHistogram2D[nsensors];
+// IHistogram2D[] ubotresidVsU = new IHistogram2D[nsensors];
+// IHistogram2D[] utopresidVsV = new IHistogram2D[nsensors];
+// IHistogram2D[] ubotresidVsV = new IHistogram2D[nsensors];
+ @Override
+ protected void detectorChanged(Detector detector) {
+
+ aida.tree().cd("/");
+ tree = aida.tree();
+ histogramFactory = IAnalysisFactory.create().createHistogramFactory(tree);
+ resetOccupancyMap();
+ for (int i = 1; i <= nmodules; i++) {
+ xresidtop[i - 1] = aida.histogram1D(plotDir + "/" + posresDir + "Module " + i + " Top x Residual", 50, -getRange(i, true), getRange(i, true));
+ yresidtop[i - 1] = aida.histogram1D(plotDir + "/" + posresDir + "Module " + i + " Top y Residual", 50, -getRange(i, false), getRange(i, false));
+ xresidbot[i - 1] = aida.histogram1D(plotDir + "/" + posresDir + "Module " + i + " Bot x Residual", 50, -getRange(i, true), getRange(i, true));
+ yresidbot[i - 1] = aida.histogram1D(plotDir + "/" + posresDir + "Module " + i + " Bot y Residual", 50, -getRange(i, false), getRange(i, false));
+ }
+
+ for (int i = 1; i <= nmodules * 2; i++) {
+ tresid[i - 1] = aida.histogram1D(plotDir + "/" + timeresDir + "HalfModule " + i + " t Residual", 50, -20, 20);
+ }
+ for (int i = 1; i <= nsensors; i++) {
+ utopresid[i - 1] = aida.histogram1D(plotDir + "/" + uresDir + "HalfModule " + i + " Top u Residual", 200, -Math.max(getRange((i + 1) / 2, false), 0.01), Math.max(getRange((i + 1) / 2, false), 0.01));
+ ubotresid[i - 1] = aida.histogram1D(plotDir + "/" + uresDir + "HalfModule " + i + " Bot u Residual", 200, -Math.max(getRange((i + 1) / 2, false), 0.01), Math.max(getRange((i + 1) / 2, false), 0.01));
+ }
+
+ ytopresidBS = aida.histogram1D(plotDir + "/" + uresDir + "BeamSpot Top y Residual", 200, -0.5, 0.5);
+ ybotresidBS = aida.histogram1D(plotDir + "/" + uresDir + "Beamspot Bot y Residual", 200, -0.5, 0.5);
+
+ xtopresidBS = aida.histogram1D(plotDir + "/" + uresDir + "BeamSpot Top x Residual", 200, -0.5, 0.5);
+ xbotresidBS = aida.histogram1D(plotDir + "/" + uresDir + "Beamspot Bot x Residual", 200, -0.5, 0.5);
+
+ }
+
+
+ @Override
+ public void process(EventHeader event) {
+ aida.tree().cd("/");
+ if (!event.hasCollection(GenericObject.class, trackResidualsCollectionName)) {
+ return;
+ }
+ nEvents++;
+ List<GenericObject> trdList = event.get(GenericObject.class, trackResidualsCollectionName);
+ for (GenericObject trd : trdList) {
+ int nResid = trd.getNDouble();
+ int isBot = trd.getIntVal(trd.getNInt() - 1);//last Int is the top/bottom flag
+ for (int i = 1; i <= nResid; i++) {
+ if (isBot == 1) {
+ xresidbot[i - 1].fill(trd.getDoubleVal(i - 1));//x is the double value in the generic object
+ yresidbot[i - 1].fill(trd.getFloatVal(i - 1));//y is the float value in the generic object
+ } else {
+ xresidtop[i - 1].fill(trd.getDoubleVal(i - 1));//x is the double value in the generic object
+ yresidtop[i - 1].fill(trd.getFloatVal(i - 1));//y is the float value in the generic object
+ }
+ }
+ }
+
+ if (event.hasCollection(GenericObject.class, trackTimeDataCollectionName)) {
+ List<GenericObject> ttdList = event.get(GenericObject.class, trackTimeDataCollectionName);
+ for (GenericObject ttd : ttdList) {
+ int nResid = ttd.getNDouble();
+ for (int i = 1; i <= nResid; i++) {
+ tresid[i - 1].fill(ttd.getDoubleVal(i - 1));//x is the double value in the generic object
+ }
+ }
+ }
+
+ if (!event.hasCollection(GenericObject.class, gblStripClusterDataCollectionName)) {
+ return;
+ }
+ List<GenericObject> gblSCDList = event.get(GenericObject.class, gblStripClusterDataCollectionName);
+ for (GenericObject gblSCD : gblSCDList) {
+ double umeas = gblSCD.getDoubleVal(GBLStripClusterData.GBLDOUBLE.UMEAS);//TODO: implement generic methods into GBLStripClusterData so this isn't hard coded
+ double utrk = gblSCD.getDoubleVal(GBLStripClusterData.GBLDOUBLE.TPOSU);//implement generic methods into GBLStripClusterData so this isn't hard coded
+ double vtrk = gblSCD.getDoubleVal(GBLStripClusterData.GBLDOUBLE.TPOSV);//implement generic methods into GBLStripClusterData so this isn't hard coded
+ double resid = umeas - utrk;
+ double tanlambda = gblSCD.getDoubleVal(GBLStripClusterData.GBLDOUBLE.TLAMBDA);//use the slope as a proxy for the top/bottom half of tracker
+
+ int i = gblSCD.getIntVal(GBLStripClusterData.GBLINT.ID);//implement generic methods into GBLStripClusterData so this isn't hard coded
+ if (i == 666) {
+ if (tanlambda > 0) {
+ xtopresidBS.fill(resid);
+ } else {
+ xbotresidBS.fill(resid);
+ }
+ } else if (i == 667) {
+ if (tanlambda > 0) {
+ ytopresidBS.fill(resid);
+ } else {
+ ybotresidBS.fill(resid);
+ }
+ } else if (tanlambda > 0) {
+ utopresid[i - 1].fill(resid);//x is the double value in the generic object
+ } // aida.histogram2D(plotDir + triggerType + "/"+uresDir + "HalfModule " + i + " Top u Residual vs. u").fill(utrk,resid);//x is the double value in the generic object
+ // aida.histogram2D(plotDir + triggerType + "/"+uresDir + "HalfModule " + i + " Top u Residual vs. v").fill(vtrk,resid);//x is the double value in the generic object
+ else {
+ ubotresid[i - 1].fill(resid);//x is the double value in the generic object
+ } // aida.histogram2D(plotDir + triggerType + "/"+uresDir + "HalfModule " + i + " Bot u Residual vs. u").fill(utrk,resid);//x is the double value in the generic object
+ // aida.histogram2D(plotDir + triggerType + "/"+uresDir + "HalfModule " + i + " Bot u Residual vs. v").fill(vtrk,resid);//x is the double value in the generic object
+ }
+ }
+
+ public void calculateEndOfRunQuantities() {
+ IAnalysisFactory analysisFactory = IAnalysisFactory.create();
+ IFitFactory fitFactory = analysisFactory.createFitFactory();
+ IFitter fitter = fitFactory.createFitter("chi2");
+ IPlotter plotterXTop = analysisFactory.createPlotterFactory().create("x-residual Top");
+ IPlotter plotterXBottom = analysisFactory.createPlotterFactory().create("x-residual Bottom");
+ IPlotter plotterYTop = analysisFactory.createPlotterFactory().create("y-residual Top");
+ IPlotter plotterYBottom = analysisFactory.createPlotterFactory().create("y-residual Bottom");
+
+ IPlotter plotterTime = analysisFactory.createPlotterFactory().create("Track Time");
+
+ plotterXTop.createRegions(2, 3);
+ IPlotterStyle pstyle = plotterXTop.style();
+ pstyle.legendBoxStyle().setVisible(false);
+ plotterXBottom.createRegions(2, 3);
+ IPlotterStyle pstyle2 = plotterXBottom.style();
+ pstyle2.legendBoxStyle().setVisible(false);
+
+ plotterYTop.createRegions(2, 3);
+ IPlotterStyle pstyle3 = plotterYTop.style();
+ pstyle3.legendBoxStyle().setVisible(false);
+ plotterYBottom.createRegions(2, 3);
+ IPlotterStyle pstyle4 = plotterYBottom.style();
+ pstyle4.legendBoxStyle().setVisible(false);
+
+ plotterTime.createRegions(3, 4);
+ IPlotterStyle pstyle5 = plotterTime.style();
+ pstyle5.legendBoxStyle().setVisible(false);
+
+ int irXTop = 0;
+ int irXBot = 0;
+ int irYTop = 0;
+ int irYBot = 0;
+ for (int i = 1; i <= nmodules; i++) {
+ IFitResult xresultTop = fitGaussian(xresidtop[i - 1], fitter, "range=\"(-1.0,1.0)\"");
+ IFitResult yresultTop = fitGaussian(yresidtop[i - 1], fitter, "range=\"(-0.5,0.5)\"");
+ IFitResult xresultBot = fitGaussian(xresidbot[i - 1], fitter, "range=\"(-1.0,1.0)\"");
+ IFitResult yresultBot = fitGaussian(yresidbot[i - 1], fitter, "range=\"(-8.0,8.0)\"");
+ if (xresultTop != null) {
+ double[] parsXTop = xresultTop.fittedParameters();
+ plotterXTop.region(irXTop).plot(xresidtop[i - 1]);
+ plotterXTop.region(irXTop).plot(xresultTop.fittedFunction());
+ irXTop++;
+ xposTopMeanResidMap.put(trackResidualsCollectionName + " " + getQuantityName(0, 0, 1, i) + "_x", parsXTop[1]);
+ xposTopSigmaResidMap.put(trackResidualsCollectionName + " " + getQuantityName(0, 1, 1, i) + "_x", parsXTop[2]);
+ }
+ if (yresultTop != null) {
+ double[] parsYTop = yresultTop.fittedParameters();
+
+ plotterYTop.region(irYTop).plot(yresidtop[i - 1]);
+ plotterYTop.region(irYTop).plot(yresultTop.fittedFunction());
+ irYTop++;
+ yposTopMeanResidMap.put(trackResidualsCollectionName + " " + getQuantityName(0, 0, 1, i) + "_y", parsYTop[1]);
+ yposTopSigmaResidMap.put(trackResidualsCollectionName + " " + getQuantityName(0, 1, 1, i) + "_y", parsYTop[2]);
+ }
+ if (xresultBot != null) {
+ double[] parsXBot = xresultBot.fittedParameters();
+ plotterXBottom.region(irXBot).plot(xresidbot[i - 1]);
+ plotterXBottom.region(irXBot).plot(xresultBot.fittedFunction());
+ irXBot++;
+ xposBotMeanResidMap.put(trackResidualsCollectionName + " " + getQuantityName(0, 0, 0, i) + "_x", parsXBot[1]);
+ xposBotSigmaResidMap.put(trackResidualsCollectionName + " " + getQuantityName(0, 1, 0, i) + "_x", parsXBot[2]);
+ }
+ if (yresultBot != null) {
+ double[] parsYBot = yresultBot.fittedParameters();
+ plotterYBottom.region(irYBot).plot(yresidbot[i - 1]);
+ plotterYBottom.region(irYBot).plot(yresultBot.fittedFunction());
+ irYBot++;
+ yposBotMeanResidMap.put(trackResidualsCollectionName + " " + getQuantityName(0, 0, 0, i) + "_y", parsYBot[1]);
+ yposBotSigmaResidMap.put(trackResidualsCollectionName + " " + getQuantityName(0, 1, 0, i) + "_y", parsYBot[2]);
+ }
+
+ }
+ int iTime = 0;
+ for (int i = 1; i <= nmodules * 2; i++) {
+ IFitResult tresult = fitGaussian(tresid[i - 1], fitter, "range=\"(-15.0,15.0)\"");
+ if (tresult != null) {
+ double[] parsTime = tresult.fittedParameters();
+ plotterTime.region(iTime).plot(tresid[i - 1]);
+ plotterTime.region(iTime).plot(tresult.fittedFunction());
+ iTime++;
+ timeMeanResidMap.put(trackTimeDataCollectionName + " " + getQuantityName(1, 0, 2, i) + "_dt", parsTime[1]);
+ timeSigmaResidMap.put(trackTimeDataCollectionName + " " + getQuantityName(1, 1, 2, i) + "_dt", parsTime[2]);
+ }
+
+ }
+
+ if (outputPlots) {
+ try {
+ plotterXTop.writeToFile(outputPlotDir + "X-Residuals-Top.png");
+ } catch (IOException ex) {
+ Logger.getLogger(SvtMonitoring.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ try {
+ plotterYTop.writeToFile(outputPlotDir + "Y-Residuals-Top.png");
+ } catch (IOException ex) {
+ Logger.getLogger(SvtMonitoring.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ try {
+ plotterXBottom.writeToFile(outputPlotDir + "X-Residuals-Bottom.png");
+ } catch (IOException ex) {
+ Logger.getLogger(SvtMonitoring.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ try {
+ plotterYBottom.writeToFile(outputPlotDir + "Y-Residuals-Bottom.png");
+ } catch (IOException ex) {
+ Logger.getLogger(SvtMonitoring.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ try {
+ plotterTime.writeToFile(outputPlotDir + "Time-Residuals.png");
+ } catch (IOException ex) {
+ Logger.getLogger(SvtMonitoring.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+ }
+
+ private String getQuantityName(int itype, int iquant, int top, int nlayer) {
+ String typeString = "position_resid";
+ String quantString = "mean_";
+ if (itype == 1) {
+ typeString = "time_resid";
+ }
+ if (iquant == 1) {
+ quantString = "sigma_";
+ }
+
+ String botString = "bot_";
+ if (top == 1) {
+ botString = "top_";
+ }
+ if (top == 2) {
+ botString = "";
+ }
+
+ String layerString = "module" + nlayer;
+ if (itype == 1) {
+ layerString = "halfmodule" + nlayer;
+ }
+
+ return typeString + quantString + botString + layerString;
+ }
+
+ public void printDQMData() {
+ LOGGER.info("TrackingResiduals::printDQMData");
+ for (Map.Entry<String, Double> entry : xposTopMeanResidMap.entrySet()) {
+ LOGGER.info(entry.getKey() + " = " + entry.getValue());
+ }
+ for (Map.Entry<String, Double> entry : xposBotMeanResidMap.entrySet()) {
+ LOGGER.info(entry.getKey() + " = " + entry.getValue());
+ }
+ for (Map.Entry<String, Double> entry : xposTopSigmaResidMap.entrySet()) {
+ LOGGER.info(entry.getKey() + " = " + entry.getValue());
+ }
+ for (Map.Entry<String, Double> entry : xposBotSigmaResidMap.entrySet()) {
+ LOGGER.info(entry.getKey() + " = " + entry.getValue());
+ }
+ for (Map.Entry<String, Double> entry : yposTopMeanResidMap.entrySet()) {
+ LOGGER.info(entry.getKey() + " = " + entry.getValue());
+ }
+ for (Map.Entry<String, Double> entry : yposBotMeanResidMap.entrySet()) {
+ LOGGER.info(entry.getKey() + " = " + entry.getValue());
+ }
+ for (Map.Entry<String, Double> entry : yposTopSigmaResidMap.entrySet()) {
+ LOGGER.info(entry.getKey() + " = " + entry.getValue());
+ }
+ for (Map.Entry<String, Double> entry : yposBotSigmaResidMap.entrySet()) {
+ LOGGER.info(entry.getKey() + " = " + entry.getValue());
+ }
+ for (Map.Entry<String, Double> entry : timeMeanResidMap.entrySet()) {
+ LOGGER.info(entry.getKey() + " = " + entry.getValue());
+ }
+ for (Map.Entry<String, Double> entry : timeSigmaResidMap.entrySet()) {
+ LOGGER.info(entry.getKey() + " = " + entry.getValue());
+ }
+ LOGGER.info("*******************************");
+ }
+
+ public void printDQMStrings() {
+ for (Map.Entry<String, Double> entry : xposTopMeanResidMap.entrySet()) {
+ LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;");
+ }
+ for (Map.Entry<String, Double> entry : xposBotMeanResidMap.entrySet()) {
+ LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;");
+ }
+ for (Map.Entry<String, Double> entry : xposTopSigmaResidMap.entrySet()) {
+ LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;");
+ }
+ for (Map.Entry<String, Double> entry : xposBotSigmaResidMap.entrySet()) {
+ LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;");
+ }
+ for (Map.Entry<String, Double> entry : yposTopMeanResidMap.entrySet()) {
+ LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;");
+ }
+ for (Map.Entry<String, Double> entry : yposBotMeanResidMap.entrySet()) {
+ LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;");
+ }
+ for (Map.Entry<String, Double> entry : yposTopSigmaResidMap.entrySet()) {
+ LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;");
+ }
+ for (Map.Entry<String, Double> entry : yposBotSigmaResidMap.entrySet()) {
+ LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;");
+ }
+ for (Map.Entry<String, Double> entry : timeMeanResidMap.entrySet()) {
+ LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;");
+ }
+ for (Map.Entry<String, Double> entry : timeSigmaResidMap.entrySet()) {
+ LOGGER.info("ALTER TABLE dqm ADD " + entry.getKey() + " double;");
+ }
+ }
+
+ private void resetOccupancyMap() {
+ xposBotMeanResidMap = new HashMap<>();
+ xposBotSigmaResidMap = new HashMap<>();
+ yposBotMeanResidMap = new HashMap<>();
+ yposBotSigmaResidMap = new HashMap<>();
+ xposTopMeanResidMap = new HashMap<>();
+ xposTopSigmaResidMap = new HashMap<>();
+ yposTopMeanResidMap = new HashMap<>();
+ yposTopSigmaResidMap = new HashMap<>();
+ timeMeanResidMap = new HashMap<>();
+ timeSigmaResidMap = new HashMap<>();
+ for (int i = 0; i < nmodules; i++) {
+ xposTopMeanResidMap.put(getQuantityName(0, 0, 1, i) + "_x", -999.);
+ yposTopMeanResidMap.put(getQuantityName(0, 0, 1, i) + "_y", -999.);
+ xposTopSigmaResidMap.put(getQuantityName(0, 1, 1, i) + "_x", -999.);
+ yposTopSigmaResidMap.put(getQuantityName(0, 1, 1, i) + "_y", -999.);
+ xposBotMeanResidMap.put(getQuantityName(0, 0, 0, i) + "_x", -999.);
+ yposBotMeanResidMap.put(getQuantityName(0, 0, 0, i) + "_y", -999.);
+ xposBotSigmaResidMap.put(getQuantityName(0, 1, 0, i) + "_x", -999.);
+ yposBotSigmaResidMap.put(getQuantityName(0, 1, 0, i) + "_y", -999.);
+ }
+ for (int i = 0; i < nmodules * 2; i++) {
+ timeMeanResidMap.put(getQuantityName(1, 0, 2, i) + "_dt", -999.);
+ timeSigmaResidMap.put(getQuantityName(1, 1, 2, i) + "_dt", -999.);
+ }
+
+ }
+
+ IFitResult fitGaussian(IHistogram1D h1d, IFitter fitter, String range) {
+ double[] init = {20.0, 0.0, 0.2};
+ IFitResult ifr = null;
+ try {
+ ifr = fitter.fit(h1d, "g", init, range);
+ } catch (RuntimeException ex) {
+ LOGGER.info(this.getClass().getSimpleName() + ": caught exception in fitGaussian");
+ }
+ return ifr;
+// double[] init = {20.0, 0.0, 1.0, 20, -1};
+// return fitter.fit(h1d, "g+p1", init, range);
+ }
+
+ private double getRange(int layer, boolean isX) {
+ double range = 2.5;
+ if (isX) {
+ switch (layer) {
+ case 1:
+ return 0.2;
+ case 2:
+ return 0.5;
+ case 3:
+ return 0.5;
+ case 4:
+ return 1.0;
+ case 5:
+ return 1.0;
+ case 6:
+ return 1.0;
+ }
+ } else {
+ switch (layer) {
+ case 1:
+ return 0.005;
+ case 2:
+ return 0.5;
+ case 3:
+ return 0.5;
+ case 4:
+ return 1.0;
+ case 5:
+ return 1.0;
+ case 6:
+ return 1.5;
+ }
+ }
+ return range;
+
+ }
+
+}
Added: java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Plots.java
=============================================================================
--- java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Plots.java (added)
+++ java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0Plots.java Wed Nov 30 11:12:55 2016
@@ -0,0 +1,638 @@
+/**
+ * Analysis driver to calculate hit efficiencies in the SVT
+ */
+/**
+ * @author mrsolt
+ *
+ */
+package org.hps.users.mrsolt;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogramFactory;
+import hep.aida.IPlotterFactory;
+import hep.aida.IPlotter;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.aida.ITree;
+import hep.physics.vec.BasicHep3Matrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import org.lcsim.detector.ITransform3D;
+import org.lcsim.detector.converter.compact.subdetector.SvtStereoLayer;
+import org.lcsim.detector.converter.compact.subdetector.HpsTracker2;
+import org.lcsim.detector.solids.Box;
+import org.lcsim.detector.solids.Point3D;
+import org.lcsim.detector.solids.Polygon3D;
+import org.lcsim.detector.tracker.silicon.ChargeCarrier;
+import org.lcsim.detector.tracker.silicon.HpsSiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiStrips;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackCross;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.geometry.Detector;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.hps.recon.tracking.CoordinateTransformations;
+import org.hps.recon.tracking.TrackUtils;
+import org.hps.recon.tracking.TrackerHitUtils;
+import org.hps.users.mrsolt.FindableTrackLayer0;
+
+public class Layer0Plots extends Driver {
+
+
+ // Use JFreeChart as the default plotting backend
+ static {
+ hep.aida.jfree.AnalysisFactory.register();
+ }
+
+ // Plotting
+ protected AIDA aida = AIDA.defaultInstance();
+ ITree tree;
+ IHistogramFactory histogramFactory;
+ IPlotterFactory plotterFactory = IAnalysisFactory.create().createPlotterFactory();
+ protected Map<String, IPlotter> plotters = new HashMap<String, IPlotter>();
+ private Map<Integer, List<SvtStereoLayer>> topStereoLayers = new HashMap<Integer, List<SvtStereoLayer>>();
+ private Map<Integer, List<SvtStereoLayer>> bottomStereoLayers = new HashMap<Integer, List<SvtStereoLayer>>();
+ private static final String SUBDETECTOR_NAME = "Tracker";
+ private List<HpsSiSensor> sensors = null;
+
+ Map<Integer, IHistogram1D> UnbiasedResidualX_top = new HashMap<Integer,IHistogram1D>();
+ Map<Integer, IHistogram1D> UnbiasedResidualX_bot = new HashMap<Integer,IHistogram1D>();
+ Map<Integer, IHistogram1D> UnbiasedResidualY_top = new HashMap<Integer,IHistogram1D>();
+ Map<Integer, IHistogram1D> UnbiasedResidualY_bot = new HashMap<Integer,IHistogram1D>();
+
+ IHistogram1D UnbiasedResidualXTop;
+ IHistogram1D UnbiasedResidualYTop;
+ IHistogram1D UnbiasedResidualXBot;
+ IHistogram1D UnbiasedResidualYBot;
+ IHistogram1D trackMomentum;
+ IHistogram1D trackMomentumPlots_top;
+ IHistogram1D trackMomentumPlots_bot;
+ IHistogram1D trackChi2;
+ IHistogram1D trackChi2Top;
+ IHistogram1D trackChi2Bot;
+ IHistogram1D LayerHitTop;
+ IHistogram1D LayerHitBot;
+ IHistogram1D ExtrapTrackTopX;
+ IHistogram1D ExtrapTrackBotX;
+ IHistogram1D ExtrapTrackTopY;
+ IHistogram1D ExtrapTrackBotY;
+ IHistogram1D Hits0;
+ IHistogram1D Hits1;
+ IHistogram1D Hits2;
+ IHistogram1D Hits3;
+
+ IHistogram1D Accept0;
+ IHistogram1D Accept1;
+ IHistogram1D Accept2;
+ IHistogram1D Accept3;
+
+ IHistogram1D Findable4Hit;
+ IHistogram1D Findable4HitLayer0;
+ IHistogram1D Findable4HitExcludeLayer0;
+ IHistogram1D Findable5Hit;
+ IHistogram1D Findable5HitLayer0;
+ IHistogram1D Findable5HitExcludeLayer0;
+ IHistogram1D Findable6Hit;
+ IHistogram1D Findable6HitLayer0;
+ IHistogram1D Findable6HitExcludeLayer0;
+ IHistogram1D Findable7Hit;
+
+
+ IHistogram1D ThreeProngP;
+ IHistogram1D InvMass;
+ IHistogram1D Psum;
+
+
+ IHistogram2D ExtrapTrackTop;
+ IHistogram2D ExtrapTrackBot;
+ IHistogram2D Layer0TopNoAcceptance;
+ IHistogram2D Layer0BotNoAcceptance;
+ IHistogram2D Layer0TopAcceptance;
+ IHistogram2D Layer0BotAcceptance;
+
+ IHistogram1D ExtrapTrackTopYHighResidual;
+ IHistogram1D HitTopYHighResidual;
+ IHistogram1D ElectronPx;
+ IHistogram1D ElectronPy;
+ IHistogram1D ElectronPz;
+ IHistogram1D ElectronP;
+ IHistogram1D PositronPx;
+ IHistogram1D PositronPy;
+ IHistogram1D PositronPz;
+ IHistogram1D PositronP;
+ IHistogram1D Phi;
+
+ private String stereoHitCollectionName = "HelicalTrackHits";
+ private String trackCollectionName = "MatchedTracks";
+
+ int num_lay = 7;
+
+ double ebeam = 1.056;
+ double zPosTop = 59.3105;
+ double zPosBot = 40.843;
+ double eMass = 0.000511;
+
+ int[] countHits = new int[num_lay];
+ int[] countAccept = new int[num_lay];
+
+ public void detectorChanged(Detector detector){
+
+ aida.tree().cd("/");
+ tree = aida.tree();
+ histogramFactory = IAnalysisFactory.create().createHistogramFactory(tree);
+
+ // Get the HpsSiSensor objects from the tracker detector element
+ sensors = detector.getSubdetector(SUBDETECTOR_NAME)
+ .getDetectorElement().findDescendants(HpsSiSensor.class);
+
+ // If the detector element had no sensors associated with it, throw
+ // an exception
+ if (sensors.size() == 0) {
+ throw new RuntimeException("No sensors were found in this detector.");
+ }
+
+ // Get the stereo layers from the geometry and build the stereo
+ // layer maps
+ List<SvtStereoLayer> stereoLayers
+ = ((HpsTracker2) detector.getSubdetector(SUBDETECTOR_NAME).getDetectorElement()).getStereoPairs();
+ for (SvtStereoLayer stereoLayer : stereoLayers) {
+ if (stereoLayer.getAxialSensor().isTopLayer()) {
+ //System.out.println("Adding stereo layer " + stereoLayer.getLayerNumber());
+ if (!topStereoLayers.containsKey(stereoLayer.getLayerNumber())) {
+ topStereoLayers.put(stereoLayer.getLayerNumber(), new ArrayList<SvtStereoLayer>());
+ }
+ topStereoLayers.get(stereoLayer.getLayerNumber()).add(stereoLayer);
+ } else {
+ if (!bottomStereoLayers.containsKey(stereoLayer.getLayerNumber())) {
+ bottomStereoLayers.put(stereoLayer.getLayerNumber(), new ArrayList<SvtStereoLayer>());
+ }
+ bottomStereoLayers.get(stereoLayer.getLayerNumber()).add(stereoLayer);
+ }
+ }
+
+ trackMomentum = aida.histogram1D("Track Momentum", 50, 0, ebeam*1.3);
+ trackMomentumPlots_top = aida.histogram1D("Track Momentum Top", 50, 0, ebeam*1.3);
+ trackMomentumPlots_bot = aida.histogram1D("Track Momentum Bot", 50, 0, ebeam*1.3);
+ UnbiasedResidualXTop = aida.histogram1D("Unbiased Residual X Top",100,-10,10);
+ UnbiasedResidualYTop = aida.histogram1D("Unbiased Residual Y Top",100,-10,10);
+ UnbiasedResidualXBot = aida.histogram1D("Unbiased Residual X Bot",100,-10,10);
+ UnbiasedResidualYBot = aida.histogram1D("Unbiased Residual Y Bot",100,-10,10);
+ trackChi2 = aida.histogram1D("Track Chi2",50,0,100);
+ trackChi2Top = aida.histogram1D("Track Chi2 Top",50,0,100);
+ trackChi2Bot = aida.histogram1D("Track Chi2 Bot",50,0,100);
+ LayerHitTop = aida.histogram1D("Layer Hit Top",8,0,8);
+ LayerHitBot = aida.histogram1D("Layer Hit Bot",8,0,8);
+
+ ExtrapTrackTopX = aida.histogram1D("Extrapolated Track X Top",50,-10,10);
+ ExtrapTrackBotX = aida.histogram1D("Extrapolated Track X Bot",50,-10,10);
+ ExtrapTrackTopY = aida.histogram1D("Extrapolated Track Y Top",50,-10,10);
+ ExtrapTrackBotY = aida.histogram1D("Extrapolated Track Y Bot",50,-10,10);
+
+ ExtrapTrackTop = aida.histogram2D("Extrapolated Track Position Top",50,-10,10,50,-10,10);
+ ExtrapTrackBot = aida.histogram2D("Extrapolated Track Position Bot",50,-10,10,50,-10,10);
+
+ Layer0TopNoAcceptance = aida.histogram2D("Extrapolated Track Position Top No Acceptance",50,-10,10,50,-10,10);
+ Layer0BotNoAcceptance = aida.histogram2D("Extrapolated Track Position Bot No Acceptance",50,-10,10,50,-10,10);
+ Layer0TopAcceptance = aida.histogram2D("Extrapolated Track Position Top Acceptance",50,-10,10,50,-10,10);
+ Layer0BotAcceptance = aida.histogram2D("Extrapolated Track Position Bot Acceptance",50,-10,10,50,-10,10);
+
+ Hits0 = aida.histogram1D("0 Hits",num_lay,0,num_lay);
+ Hits1 = aida.histogram1D("1 Hits",num_lay,0,num_lay);
+ Hits2 = aida.histogram1D("2 Hits",num_lay,0,num_lay);
+ Hits3 = aida.histogram1D("3 Hits",num_lay,0,num_lay);
+
+ Accept0 = aida.histogram1D("0 in Acceptance",num_lay,0,num_lay);
+ Accept1 = aida.histogram1D("1 in Acceptance",num_lay,0,num_lay);
+ Accept2 = aida.histogram1D("2 in Acceptance",num_lay,0,num_lay);
+ Accept3 = aida.histogram1D("3 in Acceptance",num_lay,0,num_lay);
+
+ ThreeProngP = aida.histogram1D("Total P of 3 Prong Event", 50, 0, ebeam*1.3);
+ Psum = aida.histogram1D("P Sum Trident Pair", 50, 0, ebeam*1.3);
+ InvMass = aida.histogram1D("Invariant Mass", 100, 0, 0.2);
+
+ ExtrapTrackTopYHighResidual = aida.histogram1D("Extrapolated Track Y Top High Residuals",50,-10,10);
+ HitTopYHighResidual = aida.histogram1D("Hit Position Y Top High Residuals",50,-10,10);
+
+ ElectronPx = aida.histogram1D("Electron Px", 50, -ebeam*0.5, ebeam*0.5);
+ ElectronPy = aida.histogram1D("Electron Py", 50, -ebeam*0.5, ebeam*0.5);
+ ElectronPz = aida.histogram1D("Electron Pz", 50, 0, ebeam*1.3);
+ ElectronP = aida.histogram1D("Electron P", 50, 0, ebeam*1.3);
+ PositronPx = aida.histogram1D("Positron Px", 50, -ebeam*0.5, ebeam*0.5);
+ PositronPy = aida.histogram1D("Positron Py", 50, -ebeam*0.5, ebeam*0.5);
+ PositronPz = aida.histogram1D("Positron Pz", 50, 0, ebeam*1.3);
+ PositronP = aida.histogram1D("Positron P", 50, 0, ebeam*1.3);
+
+ Findable4Hit = aida.histogram1D("Findable Tracks 4 Hits", num_lay, 0, num_lay);
+ Findable4HitLayer0 = aida.histogram1D("Findable Tracks 4 Hits Layer 0 Included", num_lay, 0, num_lay);
+ Findable4HitExcludeLayer0 = aida.histogram1D("Findable Tracks 4 Hits Layer 0 Excluded", num_lay, 0, num_lay);
+ Findable5Hit = aida.histogram1D("Findable Tracks 5 Hits", num_lay, 0, num_lay);
+ Findable5HitLayer0 = aida.histogram1D("Findable Tracks 5 Hits Layer 0 Included", num_lay, 0, num_lay);
+ Findable5HitExcludeLayer0 = aida.histogram1D("Findable Tracks 5 Hits Layer 0 Excluded", num_lay, 0, num_lay);
+ Findable6Hit = aida.histogram1D("Findable Tracks 6 Hits", num_lay, 0, num_lay);
+ Findable6HitLayer0 = aida.histogram1D("Findable Tracks 6 Hits Layer 0 Included", num_lay, 0, num_lay);
+ Findable6HitExcludeLayer0 = aida.histogram1D("Findable Tracks 6 Hits Layer 0 Excluded", num_lay, 0, num_lay);
+ Findable7Hit = aida.histogram1D("Findable Tracks 7 Hits", num_lay, 0, num_lay);
+
+
+ Phi = aida.histogram1D("Phi", 500, -2*Math.PI, 2*Math.PI);
+
+
+ for(int i = 0; i < num_lay; i++){
+ UnbiasedResidualX_top.put((i+1),histogramFactory.createHistogram1D("Unbiased Residual X Top Layer #" + (i+1),100,-10,10));
+ UnbiasedResidualX_bot.put((i+1),histogramFactory.createHistogram1D("Unbiased Residual X Bot Layer #" + (i+1),100,-10,10));
+ UnbiasedResidualY_top.put((i+1),histogramFactory.createHistogram1D("Unbiased Residual Y Top Layer #" + (i+1),100,-10,10));
+ UnbiasedResidualY_bot.put((i+1),histogramFactory.createHistogram1D("Unbiased Residual Y Bot Layer #" + (i+1),100,-10,10));
+ }
+
+ for (IPlotter plotter : plotters.values()) {
+ //plotter.show();
+ }
+ }
+
+
+ @Override
+ public void process(EventHeader event){
+ aida.tree().cd("/");
+
+ System.out.println("Pass1");
+
+ // If the event does not have tracks, skip it
+ //if(!event.hasCollection(Track.class, trackCollectionName)) return;
+
+ // Get the list of tracks from the event
+ List<List<Track>> trackCollections = event.get(Track.class);
+
+ List<Track> tracks = event.get(Track.class, trackCollectionName);
+
+ List<TrackerHit> stereoHits = event.get(TrackerHit.class, stereoHitCollectionName);
+
+ for(int i = 0; i<num_lay; i++){
+ countHits[i] = 0;
+ countAccept[i] = 0;
+ }
+
+ double totalPx = 0;
+ double totalPy = 0;
+ double totalPz = 0;
+ double EleP = 0;
+ double ElePx = 0;
+ double ElePy = 0;
+ double ElePz = 0;
+ double PosP = 0;
+ double PosPx = 0;
+ double PosPy = 0;
+ double PosPz = 0;
+ int electrons = 0;
+ int positrons = 0;
+
+ //for (List<Track> tracks : trackCollections) {
+
+ System.out.println("Pass2");
+
+ for(Track track : tracks){
+
+ HpsSiSensor trackSensor = (HpsSiSensor) ((RawTrackerHit)track.getTrackerHits().get(0).getRawHits().get(0)).getDetectorElement();
+
+
+ // Calculate the momentum of the track
+ HelicalTrackFit hlc_trk_fit = TrackUtils.getHTF(track.getTrackStates().get(0));
+
+ System.out.println("Pass3");
+
+ double B_field = Math.abs(TrackUtils.getBField(event.getDetector()).y());
+
+ double p = hlc_trk_fit.p(B_field);
+ double pT = hlc_trk_fit.pT(B_field);
+ double phi = hlc_trk_fit.phi0();
+ double tanLambda = track.getTrackStates().get(0).getTanLambda();
+
+ //totalPx = totalPx + pT * Math.cos(phi);
+ totalPy = totalPy + pT * Math.sin(phi);
+ //totalPz = totalPz + Math.sqrt(p*p-pT*pT);
+ totalPz = pT * tanLambda;
+ totalPx = totalPx + Math.sqrt(p*p-Math.pow(pT*Math.sin(phi), 2)-Math.pow(pT*tanLambda, 2));
+
+ System.out.println("Pass4");
+
+ if(track.getCharge() < 0){
+ EleP = p;
+ //ElePx = pT * Math.cos(phi);
+ ElePy = pT * Math.sin(phi);
+ //ElePz = Math.sqrt(p*p-pT*pT);
+ ElePz = pT * tanLambda;
+ ElePx = Math.sqrt(p*p-ElePy*ElePy-ElePz*ElePz);
+ electrons++;
+ }
+ else{
+ PosP = p;
+ //PosPx = pT * Math.cos(phi);
+ PosPy = pT * Math.sin(phi);
+ //PosPz = Math.sqrt(p*p-pT*pT);
+ PosPz = pT * tanLambda;
+ PosPx = Math.sqrt(p*p-PosPy*PosPy-PosPz*PosPz);
+ positrons++;
+ }
+
+ System.out.println("Pass5");
+
+ Phi.fill(phi);
+ trackMomentum.fill(p);
+ trackChi2.fill(track.getChi2());
+ if (trackSensor.isTopLayer()) {
+ Hep3Vector extrapTrackPos = TrackUtils.extrapolateTrack(track, zPosTop);
+ double extrapTrackX = extrapTrackPos.x();
+ double extrapTrackY = extrapTrackPos.y();
+ trackMomentumPlots_top.fill(p);
+ trackChi2Top.fill(track.getChi2());
+ ExtrapTrackTopX.fill(extrapTrackX);
+ ExtrapTrackTopY.fill(extrapTrackY);
+ ExtrapTrackTop.fill(extrapTrackX,extrapTrackY);
+ } else {
+ trackMomentumPlots_bot.fill(p);
+ Hep3Vector extrapTrackPos = TrackUtils.extrapolateTrack(track, zPosBot);
+ double extrapTrackX = extrapTrackPos.x();
+ double extrapTrackY = extrapTrackPos.y();
+ trackChi2Bot.fill(track.getChi2());
+ ExtrapTrackBotX.fill(extrapTrackX);
+ ExtrapTrackBotY.fill(extrapTrackY);
+ ExtrapTrackBot.fill(extrapTrackX,extrapTrackY);
+ }
+ System.out.println("Pass6");
+ for(int i = 0; i < num_lay; i++){
+ if(isWithinAcceptance(track, i + 1)){
+ countAccept[i]++;
+ if(i == 0){
+ if (trackSensor.isTopLayer()) {
+ Hep3Vector extrapTrackPos = TrackUtils.extrapolateTrack(track, zPosTop);
+ double extrapTrackX = extrapTrackPos.x();
+ double extrapTrackY = extrapTrackPos.y();
+ Layer0TopAcceptance.fill(extrapTrackX,extrapTrackY);
+ }
+ else{
+ Hep3Vector extrapTrackPos = TrackUtils.extrapolateTrack(track, zPosBot);
+ double extrapTrackX = extrapTrackPos.x();
+ double extrapTrackY = extrapTrackPos.y();
+ Layer0BotAcceptance.fill(extrapTrackX,extrapTrackY);
+ }
+ }
+ }
+ else{
+ if(i == 0){
+ if (trackSensor.isTopLayer()) {
+ Hep3Vector extrapTrackPos = TrackUtils.extrapolateTrack(track, zPosTop);
+ double extrapTrackX = extrapTrackPos.x();
+ double extrapTrackY = extrapTrackPos.y();
+ Layer0TopNoAcceptance.fill(extrapTrackX,extrapTrackY);
+ }
+ else{
+ Hep3Vector extrapTrackPos = TrackUtils.extrapolateTrack(track, zPosBot);
+ double extrapTrackX = extrapTrackPos.x();
+ double extrapTrackY = extrapTrackPos.y();
+ Layer0BotNoAcceptance.fill(extrapTrackX,extrapTrackY);
+ }
+ }
+ }
+ }
+ System.out.println("Pass7");
+
+ for (TrackerHit stereoHit : stereoHits) {
+ HpsSiSensor hitSensor = (HpsSiSensor) ((RawTrackerHit) stereoHit.getRawHits().get(0)).getDetectorElement();
+ System.out.println("Pass8");
+ // Retrieve the layer number by using the sensor
+ int layer = (hitSensor.getLayerNumber() + 1)/2;
+ countHits[layer-1]++;
+
+ Hep3Vector stereoHitPosition = new BasicHep3Vector(stereoHit.getPosition());
+ Hep3Vector trackPosition = TrackUtils.extrapolateTrack(track, stereoHitPosition.z());
+ double xResidual = trackPosition.x() - stereoHitPosition.x();
+ double yResidual = trackPosition.y() - stereoHitPosition.y();
+ System.out.println("Pass9");
+ if (hitSensor.isTopLayer()) {
+ UnbiasedResidualX_top.get(layer).fill(xResidual);
+ UnbiasedResidualY_top.get(layer).fill(yResidual);
+ UnbiasedResidualXTop.fill(xResidual);
+ UnbiasedResidualYTop.fill(yResidual);
+ LayerHitTop.fill(layer);
+ if(layer == 1){
+ if(yResidual<-2){
+ ExtrapTrackTopYHighResidual.fill(trackPosition.y());
+ HitTopYHighResidual.fill(stereoHitPosition.y());
+ }
+ }
+ }
+ else{
+ UnbiasedResidualX_bot.get(layer).fill(xResidual);
+ UnbiasedResidualY_bot.get(layer).fill(yResidual);
+ UnbiasedResidualXBot.fill(xResidual);
+ UnbiasedResidualYBot.fill(yResidual);
+ LayerHitBot.fill(layer);
+ }
+ }
+ System.out.println("Pass10");
+ }
+ double massSqr = EleP*EleP + PosP*PosP + 2*eMass*eMass + 2*Math.sqrt(EleP*EleP+eMass*eMass)*Math.sqrt(PosP*PosP+eMass*eMass)-(Math.pow(ElePx+PosPx,2)+Math.pow(ElePy+PosPy,2)+Math.pow(ElePz+PosPz,2));
+ //double massSqr = Math.pow(EleP+PosP,2) - (Math.pow(ElePx+PosPx,2)+Math.pow(ElePy+PosPy,2)+Math.pow(ElePz+PosPz,2));
+ //double EleE = Math.sqrt(EleP*EleP + eMass*eMass);
+ //double PosE = Math.sqrt(PosP*PosP + eMass*eMass);
+ //double BetaE =;
+ //double BetaP;
+ //double theta;
+ //double massSqr = 2*(eMass*eMass + EleE*PosE*(1 - BetaE*BetaP*Math.cos(theta)));
+ //InvMass.fill(Math.sqrt(massSqr));
+ System.out.println("Pass11");
+ for(int i = 0; i < num_lay; i++){
+ if(countHits[i] == 3){
+ Hits3.fill(i);
+ }
+ else if(countHits[i] == 2){
+ Hits2.fill(i);
+ }
+ else if(countHits[i] == 1){
+ Hits1.fill(i);
+ }
+ else if(countHits[i] == 0){
+ Hits0.fill(i);
+ }
+ }
+ System.out.println("Pass12");
+ for(int i = 0; i < num_lay; i++){
+ if(countAccept[i] == 3){
+ Accept3.fill(i);
+ if(i == 0){
+ double totalP = Math.sqrt(totalPx*totalPx + totalPy*totalPy + totalPz*totalPz);
+ ThreeProngP.fill(totalP);
+ }
+ }
+ else if(countAccept[i] == 2){
+ Accept2.fill(i);
+ }
+ else if(countAccept[i] == 1){
+ Accept1.fill(i);
+ }
+ else if(countAccept[i] == 0){
+ Accept0.fill(i);
+ }
+ }
+ System.out.println("Pass13");
+ if(electrons == 1 && positrons == 1 && massSqr > 0){
+ InvMass.fill(Math.sqrt(massSqr));
+ Psum.fill(EleP+PosP);
+ ElectronPx.fill(ElePx);
+ ElectronPy.fill(ElePy);
+ ElectronPz.fill(ElePz);
+ ElectronP.fill(EleP);
+ PositronPx.fill(PosPx);
+ PositronPy.fill(PosPy);
+ PositronPz.fill(PosPz);
+ PositronP.fill(PosP);
+ }
+ System.out.println("Pass14");
+ }
+ //}
+
+
+ @Override
+ public void endOfData(){
+ System.out.println("%===================================================================%");
+ System.out.println("%====================== Layer 0 Plots Complete==========================%");
+ System.out.println("%===================================================================% \n%");
+ System.out.println("% \n%===================================================================%");
+ }
+
+ private boolean isWithinAcceptance(Track track, int layer) {
+
+ // TODO: Move this to a utility class
+
+ //System.out.println("Retrieving sensors for layer: " + layer);
+
+ // Since TrackUtils.isTop/BottomTrack does not work when running off
+ // a recon file, get the detector volume that a track is associated
+ // with by using the sensor. This assumes that a track is always
+ // composed by stereo hits that lie within a single detector volume
+ HpsSiSensor sensor = (HpsSiSensor) ((RawTrackerHit)track.getTrackerHits().get(0).getRawHits().get(0)).getDetectorElement();
+
+ // Get the sensors associated with the layer that the track
+ // will be extrapolated to
+ List<SvtStereoLayer> stereoLayers = null;
+
+ // if (TrackUtils.isTopTrack(track, track.getTrackerHits().size())) {
+ if (sensor.isTopLayer()) {
+ //System.out.println("Top track found.");
+ stereoLayers = this.topStereoLayers.get(layer);
+ //} else if (TrackUtils.isBottomTrack(track, track.getTrackerHits().size())) {
+ } else {
+ //System.out.println("Bottom track found.");
+ stereoLayers = this.bottomStereoLayers.get(layer);
+ }
+
+ for (SvtStereoLayer stereoLayer : stereoLayers) {
+ Hep3Vector axialSensorPosition = stereoLayer.getAxialSensor().getGeometry().getPosition();
+ Hep3Vector stereoSensorPosition = stereoLayer.getStereoSensor().getGeometry().getPosition();
+
+ //System.out.println("Axial sensor position: " + axialSensorPosition.toString());
+ //System.out.println("Stereo sensor position: " + stereoSensorPosition.toString());
+
+ Hep3Vector axialTrackPos = TrackUtils.extrapolateTrack(track, axialSensorPosition.z());
+ Hep3Vector stereoTrackPos = TrackUtils.extrapolateTrack(track, stereoSensorPosition.z());
+
+ //System.out.println("Track position at axial sensor: " + axialTrackPos.toString());
+ //System.out.println("Track position at stereo sensor: " + stereoTrackPos.toString());
+
+ if(this.sensorContainsTrack(axialTrackPos, stereoLayer.getAxialSensor())
+ && this.sensorContainsTrack(stereoTrackPos, stereoLayer.getStereoSensor())){
+ //System.out.println("Track lies within layer acceptance.");
+ return true;
+ }
+ }
+
+ return false;
+
+ /*int layerNumber = (layer - 1)/2 + 1;
+ String title = "Track Position - Layer " + layerNumber + " - Tracks Within Acceptance";
+ //aida.histogram2D(title).fill(trackPos.y(), trackPos.z());
+ //aida.cloud2D(title).fill(frontTrackPos.y(), frontTrackPos.z()); */
+
+ }
+
+ public boolean sensorContainsTrack(Hep3Vector trackPosition, HpsSiSensor sensor){
+
+
+ ITransform3D localToGlobal = sensor.getGeometry().getLocalToGlobal();
+
+ Hep3Vector sensorPos = sensor.getGeometry().getPosition();
+ Box sensorSolid = (Box) sensor.getGeometry().getLogicalVolume().getSolid();
+ Polygon3D sensorFace = sensorSolid.getFacesNormalTo(new BasicHep3Vector(0, 0, 1)).get(0);
+
+ List<Point3D> vertices = new ArrayList<Point3D>();
+ for(int index = 0; index < 4; index++){
+ vertices.add(new Point3D());
+ }
+ for(Point3D vertex : sensorFace.getVertices()){
+ if(vertex.y() < 0 && vertex.x() > 0){
+ localToGlobal.transform(vertex);
+ //vertices.set(0, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(0, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ //System.out.println(this.getClass().getSimpleName() + ": Vertex 1 Position: " + vertices.get(0).toString());
+ //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 1 Position: " + localToGlobal.transformed(vertex).toString());
+ }
+ else if(vertex.y() > 0 && vertex.x() > 0){
+ localToGlobal.transform(vertex);
+ //vertices.set(1, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(1, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ //System.out.println(this.getClass().getSimpleName() + ": Vertex 2 Position: " + vertices.get(1).toString());
+ //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 2 Position: " + localToGlobal.transformed(vertex).toString());
+ }
+ else if(vertex.y() > 0 && vertex.x() < 0){
+ localToGlobal.transform(vertex);
+ //vertices.set(2, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(2, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ //System.out.println(this.getClass().getSimpleName() + ": Vertex 3 Position: " + vertices.get(2).toString());
+ //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 3 Position: " + localToGlobal.transformed(vertex).toString());
+ }
+ else if(vertex.y() < 0 && vertex.x() < 0){
+ localToGlobal.transform(vertex);
+ //vertices.set(3, new Point3D(vertex.y() + sensorPos.x(), vertex.x() + sensorPos.y(), vertex.z() + sensorPos.z()));
+ vertices.set(3, new Point3D(vertex.x(), vertex.y(), vertex.z()));
+ //System.out.println(this.getClass().getSimpleName() + ": Vertex 4 Position: " + vertices.get(3).toString());
+ //System.out.println(this.getClass().getSimpleName() + ": Transformed Vertex 4 Position: " + localToGlobal.transformed(vertex).toString());
+ }
+ }
+
+ /*
+ double area1 = this.findTriangleArea(vertices.get(0).x(), vertices.get(0).y(), vertices.get(1).x(), vertices.get(1).y(), trackPosition.y(), trackPosition.z());
+ double area2 = this.findTriangleArea(vertices.get(1).x(), vertices.get(1).y(), vertices.get(2).x(), vertices.get(2).y(), trackPosition.y(), trackPosition.z());
+ double area3 = this.findTriangleArea(vertices.get(2).x(), vertices.get(2).y(), vertices.get(3).x(), vertices.get(3).y(), trackPosition.y(), trackPosition.z());
+ double area4 = this.findTriangleArea(vertices.get(3).x(), vertices.get(3).y(), vertices.get(0).x(), vertices.get(0).y(), trackPosition.y(), trackPosition.z());
+ */
+
+ double area1 = this.findTriangleArea(vertices.get(0).x(), vertices.get(0).y(), vertices.get(1).x(), vertices.get(1).y(), trackPosition.x(), trackPosition.y());
+ double area2 = this.findTriangleArea(vertices.get(1).x(), vertices.get(1).y(), vertices.get(2).x(), vertices.get(2).y(), trackPosition.x(), trackPosition.y());
+ double area3 = this.findTriangleArea(vertices.get(2).x(), vertices.get(2).y(), vertices.get(3).x(), vertices.get(3).y(), trackPosition.x(), trackPosition.y());
+ double area4 = this.findTriangleArea(vertices.get(3).x(), vertices.get(3).y(), vertices.get(0).x(), vertices.get(0).y(), trackPosition.x(), trackPosition.y());
+
+ if((area1 > 0 && area2 > 0 && area3 > 0 && area4 > 0) || (area1 < 0 && area2 < 0 && area3 < 0 && area4 < 0)) return true;
+ return false;
+ }
+
+ public double findTriangleArea(double x0, double y0, double x1, double y1, double x2, double y2){
+ return .5*(x1*y2 - y1*x2 -x0*y2 + y0*x2 + x0*y1 - y0*x1);
+ }
+
+ }
Added: java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0ReconParticle.java
=============================================================================
--- java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0ReconParticle.java (added)
+++ java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/Layer0ReconParticle.java Wed Nov 30 11:12:55 2016
@@ -0,0 +1,557 @@
+/**
+ * Analysis driver to calculate hit efficiencies in the SVT
+ */
+/**
+ * @author mrsolt
+ *
+ */
+
+package org.hps.users.mrsolt;
+
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogramFactory;
+import hep.aida.IPlotterFactory;
+import hep.aida.IPlotter;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.aida.ITree;
+import hep.physics.vec.BasicHep3Matrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+
+import org.lcsim.detector.converter.compact.subdetector.SvtStereoLayer;
+import org.lcsim.detector.tracker.silicon.HpsSiSensor;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.ReconstructedParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.Vertex;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.hps.recon.tracking.TrackType;
+import org.hps.recon.tracking.TrackUtils;
+import org.hps.recon.tracking.TrackerHitUtils;
+import org.hps.recon.vertexing.BilliorTrack;
+import org.hps.recon.vertexing.BilliorVertex;
+import org.hps.recon.vertexing.BilliorVertexer;
+import org.hps.users.mrsolt.FindableTrackLayer0;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HitIdentifier;
+
+
+public class Layer0ReconParticle extends Driver {
+
+
+ // Use JFreeChart as the default plotting backend
+ static {
+ hep.aida.jfree.AnalysisFactory.register();
+ }
+
+ // Plotting
+ protected AIDA aida = AIDA.defaultInstance();
+ ITree tree;
+ IHistogramFactory histogramFactory;
+ IPlotterFactory plotterFactory = IAnalysisFactory.create().createPlotterFactory();
+ protected Map<String, IPlotter> plotters = new HashMap<String, IPlotter>();
+ private Map<Integer, List<SvtStereoLayer>> topStereoLayers = new HashMap<Integer, List<SvtStereoLayer>>();
+ private Map<Integer, List<SvtStereoLayer>> bottomStereoLayers = new HashMap<Integer, List<SvtStereoLayer>>();
+ private static final String SUBDETECTOR_NAME = "Tracker";
+ RelationalTable<SimTrackerHit, MCParticle> _hittomc = null;
+ private List<HpsSiSensor> sensors = null;
+ private double _bfield;
+ private HitIdentifier _ID;
+ private int _nlayersTot = 14;
+ double eBeam = 1.05;
+ double maxFactor = 1.25;
+
+
+ IHistogram1D unconMass;
+ IHistogram1D unconVx;
+ IHistogram1D unconVy;
+ IHistogram1D unconVz;
+ IHistogram1D unconChi2;
+ IHistogram2D unconVzVsChi2;
+ IHistogram2D unconChi2VsTrkChi2;
+ /* beamspot constrained */
+
+ IHistogram1D nV0;
+
+ IHistogram1D v0Time;
+ IHistogram1D v0Dt;
+ IHistogram2D trigTimeV0Time;
+ IHistogram1D trigTime;
+
+ IHistogram1D bsconMass;
+ IHistogram1D bsconVx;
+ IHistogram1D bsconVy;
+ IHistogram1D bsconVz;
+ IHistogram1D bsconChi2;
+ IHistogram2D bsconVzVsChi2;
+ IHistogram2D bsconChi2VsTrkChi2;
+ /* target constrained */
+ IHistogram1D tarconMass;
+ IHistogram1D tarconVx;
+ IHistogram1D tarconVy;
+ IHistogram1D tarconVz;
+ IHistogram1D tarconChi2;
+ IHistogram2D tarconVzVsChi2;
+ IHistogram2D tarconChi2VsTrkChi2;
+
+ IHistogram2D pEleVspPos;
+ IHistogram2D pEleVspPosWithCut;
+ IHistogram2D pyEleVspyPos;
+ IHistogram2D pxEleVspxPos;
+
+ IHistogram2D VtxZVsMass;
+ IHistogram2D VtxYVsVtxZ;
+ IHistogram2D VtxXVsVtxZ;
+ IHistogram2D VtxXVsVtxY;
+ IHistogram2D VtxXVsVtxPx;
+ IHistogram2D VtxYVsVtxPy;
+ IHistogram2D VtxZVsVtxPx;
+ IHistogram2D VtxZVsVtxPy;
+ IHistogram2D VtxZVsVtxPz;
+
+ IHistogram2D VtxZVsL1Iso;
+ IHistogram2D VtxZVsTrkChi2;
+
+ IHistogram2D pEleVspEle;
+ IHistogram2D phiEleVsphiEle;
+ IHistogram2D pyEleVspyEle;
+ IHistogram2D pxEleVspxEle;
+ IHistogram2D pEleVspEleNoBeam;
+ IHistogram2D pyEleVspyEleNoBeam;
+ IHistogram2D pxEleVspxEleNoBeam;
+ IHistogram2D pEleVspEleMoller;
+ IHistogram2D pEleVsthetaMoller;
+ IHistogram2D thetaEleVsthetaMoller;
+ IHistogram2D pEleVspEleBeamBeam;
+ IHistogram2D pEleVsthetaBeamBeam;
+ IHistogram2D thetaEleVsthetaBeamBeam;
+
+ IHistogram1D mollerMass;
+ IHistogram1D mollerMassVtxCut;
+ IHistogram1D mollerVx;
+ IHistogram1D mollerVy;
+ IHistogram1D mollerVz;
+ IHistogram1D mollerVzVtxCut;
+ IHistogram2D mollerXVsVtxZ;
+ IHistogram2D mollerYVsVtxZ;
+ IHistogram2D mollerXVsVtxY;
+
+ IHistogram1D sumChargeHisto;
+ IHistogram1D numChargeHisto;
+
+ double v0ESumMinCut = 0.8 * eBeam;
+ double v0ESumMaxCut = 1.25 * eBeam;
+ double v0MaxPCut = 1.1;//GeV
+ double molPSumMin = 0.85;
+ double molPSumMax = 1.3;
+ double beambeamCut = 0.85;
+ double thetaMax = 0.06;
+ double thetaMin = 0.015;
+
+ private String stereoHitCollectionName = "HelicalTrackHits";
+ private String trackCollectionName = "MatchedTracks";
+ private String MCparticleCollectionName = "MCParticle";
+ String finalStateParticlesColName = "FinalStateParticles";
+ String unconstrainedV0CandidatesColName = "UnconstrainedV0Candidates";
+ String beamConV0CandidatesColName = "BeamspotConstrainedV0Candidates";
+ String targetV0ConCandidatesColName = "TargetConstrainedV0Candidates";
+ private final BasicHep3Matrix beamAxisRotation = new BasicHep3Matrix();
+
+ //some counters
+ int nRecoEvents = 0;
+ int nTotV0 = 0;
+ int nTot2Ele = 0;
+ //some summers
+ double sumMass = 0.0;
+ double sumVx = 0.0;
+ double sumVy = 0.0;
+ double sumVz = 0.0;
+ double sumChi2 = 0.0;
+
+ int num_lay = 7;
+
+ //double ebeam = 1.056;
+ double zPosTop = 59.3105;
+ double zPosBot = 40.843;
+
+
+ public void detectorChanged(Detector detector){
+
+ beamAxisRotation.setActiveEuler(Math.PI / 2, -0.0305, -Math.PI / 2);
+ aida.tree().cd("/");
+ tree = aida.tree();
+ histogramFactory = IAnalysisFactory.create().createHistogramFactory(tree);
+
+ // Get the HpsSiSensor objects from the tracker detector element
+ sensors = detector.getSubdetector(SUBDETECTOR_NAME)
+ .getDetectorElement().findDescendants(HpsSiSensor.class);
+
+ // If the detector element had no sensors associated with it, throw
+ // an exception
+ if (sensors.size() == 0) {
+ throw new RuntimeException("No sensors were found in this detector.");
+ }
+
+ unconMass = aida.histogram1D(unconstrainedV0CandidatesColName + "/" + "Invariant Mass (GeV)", 100, 0, 0.200);
+ unconVx = aida.histogram1D(unconstrainedV0CandidatesColName + "/" + "Vx (mm)", 50, -10, 10);
+ unconVy = aida.histogram1D(unconstrainedV0CandidatesColName + "/" + "Vy (mm)", 50, -10, 10);
+ unconVz = aida.histogram1D(unconstrainedV0CandidatesColName + "/" + "Vz (mm)", 50, -50, 50);
+ unconChi2 = aida.histogram1D(unconstrainedV0CandidatesColName + "/" + "Chi2", 25, 0, 25);
+ unconVzVsChi2 = aida.histogram2D(unconstrainedV0CandidatesColName + "/" + "Vz vs. Chi2", 25, 0, 25, 50, -50, 50);
+ unconChi2VsTrkChi2 = aida.histogram2D(unconstrainedV0CandidatesColName + "/" + "Chi2 vs. total track chi2", 50, 0, 50, 50, 0, 25);
+ /* beamspot constrained */
+ bsconMass = aida.histogram1D(beamConV0CandidatesColName + "/" + "Mass (GeV)", 100, 0, 0.200);
+ bsconVx = aida.histogram1D(beamConV0CandidatesColName + "/" + "Vx (mm)", 50, -10, 10);
+ bsconVy = aida.histogram1D(beamConV0CandidatesColName + "/" + "Vy (mm)", 50, -10, 10);
+ bsconVz = aida.histogram1D(beamConV0CandidatesColName + "/" + "Vz (mm)", 50, -50, 50);
+ bsconChi2 = aida.histogram1D(beamConV0CandidatesColName + "/" + "Chi2", 25, 0, 25);
+ bsconVzVsChi2 = aida.histogram2D(beamConV0CandidatesColName + "/" + "Vz vs. Chi2", 25, 0, 25, 50, -50, 50);
+ bsconChi2VsTrkChi2 = aida.histogram2D(beamConV0CandidatesColName + "/" + "Chi2 vs. total track chi2", 50, 0, 50, 50, 0, 25);
+ /* target constrained */
+ tarconMass = aida.histogram1D(targetV0ConCandidatesColName + "/" + "Mass (GeV)", 100, 0, 0.200);
+ tarconVx = aida.histogram1D(targetV0ConCandidatesColName + "/" + "Vx (mm)", 50, -1, 1);
+ tarconVy = aida.histogram1D(targetV0ConCandidatesColName + "/" + "Vy (mm)", 50, -1, 1);
+ tarconVz = aida.histogram1D(targetV0ConCandidatesColName + "/" + "Vz (mm)", 50, -10, 10);
+ tarconChi2 = aida.histogram1D(targetV0ConCandidatesColName + "/" + "Chi2", 25, 0, 25);
+ tarconVzVsChi2 = aida.histogram2D(targetV0ConCandidatesColName + "/" + "Vz vs. Chi2", 25, 0, 25, 50, -50, 50);
+ tarconChi2VsTrkChi2 = aida.histogram2D(targetV0ConCandidatesColName + "/" + "Chi2 vs. total track chi2", 50, 0, 50, 50, 0, 25);
+
+ nV0 = aida.histogram1D("Number of V0 per event", 10, 0, 10);
+ v0Time = aida.histogram1D("V0 mean time", 100, -25, 25);
+ v0Dt = aida.histogram1D("V0 time difference", 100, -25, 25);
+ trigTimeV0Time = aida.histogram2D("Trigger phase vs. V0 mean time", 100, -25, 25, 6, 0, 24);
+ trigTime = aida.histogram1D("Trigger phase", 6, 0, 24);
+
+ pEleVspPos = aida.histogram2D("P(e) vs P(p)", 50, 0, eBeam * maxFactor, 50, 0, eBeam * maxFactor);
+ pEleVspPosWithCut = aida.histogram2D("P(e) vs P(p): Radiative", 50, 0, eBeam * maxFactor, 50, 0, eBeam * maxFactor);
+ pyEleVspyPos = aida.histogram2D("Py(e) vs Py(p)", 50, -0.04, 0.04, 50, -0.04, 0.04);
+ pxEleVspxPos = aida.histogram2D("Px(e) vs Px(p)", 50, -0.04, 0.04, 50, -0.04, 0.04);
+ VtxZVsMass = aida.histogram2D("Vz vs Mass", 50, 0, 0.15, 50, -50, 80);
+ VtxXVsVtxZ = aida.histogram2D("Vx vs Vz", 100, -10, 10, 100, -50, 80);
+ VtxYVsVtxZ = aida.histogram2D("Vy vs Vz", 100, -5, 5, 100, -50, 80);
+ VtxXVsVtxY = aida.histogram2D("Vx vs Vy", 100, -10, 10, 100, -5, 5);
+ VtxXVsVtxPx = aida.histogram2D("Vx vs Px", 100, -0.1, 0.1, 100, -10, 10);
+ VtxYVsVtxPy = aida.histogram2D("Vy vs Py", 100, -0.1, 0.1, 100, -5, 5);
+ VtxZVsVtxPx = aida.histogram2D("Vz vs Px", 100, -0.1, 0.1, 100, -50, 80);
+ VtxZVsVtxPy = aida.histogram2D("Vz vs Py", 100, -0.1, 0.1, 100, -50, 80);
+ VtxZVsVtxPz = aida.histogram2D("Vz vs Pz", 100, 0.0, eBeam * maxFactor, 100, -50, 80);
+ VtxZVsL1Iso = aida.histogram2D("Vz vs L1 Isolation", 100, 0.0, 5.0, 50, -50, 80);
+ VtxZVsTrkChi2 = aida.histogram2D("Vz vs Track Chi2", 50, 0, 50, 50, -50, 80);
+ pEleVspEle = aida.histogram2D("2 Electron/P(e) vs P(e)", 50, 0, eBeam * maxFactor, 50, 0, eBeam * maxFactor);
+ phiEleVsphiEle = aida.histogram2D("2 Electron/phi(e) vs phi(e)", 50, -Math.PI, Math.PI, 50, -Math.PI, Math.PI);
+ pyEleVspyEle = aida.histogram2D("2 Electron/Py(e) vs Py(e)", 50, -0.04, 0.04, 50, -0.04, 0.04);
+ pxEleVspxEle = aida.histogram2D("2 Electron/Px(e) vs Px(e)", 50, -0.02, 0.06, 50, -0.02, 0.06);
+ pEleVspEleNoBeam = aida.histogram2D("2 Electron/P(e) vs P(e) NoBeam", 50, 0, beambeamCut, 50, 0, beambeamCut);
+ pEleVspEleMoller = aida.histogram2D("2 Electron/P(e) vs P(e) Moller", 50, 0, beambeamCut, 50, 0, beambeamCut);
+ pEleVspEleBeamBeam = aida.histogram2D("2 Electron/P(e) vs P(e) BeamBeam", 50, beambeamCut, eBeam * maxFactor, 50, beambeamCut, eBeam * maxFactor);
+ pyEleVspyEleNoBeam = aida.histogram2D("2 Electron/Py(e) vs Py(e) NoBeam", 50, -0.04, 0.04, 50, -0.04, 0.04);
+ pxEleVspxEleNoBeam = aida.histogram2D("2 Electron/Px(e) vs Px(e) NoBeam", 50, -0.02, 0.06, 50, -0.02, 0.06);
+ sumChargeHisto = aida.histogram1D("Total Charge of Event", 5, -2, 3);
+ numChargeHisto = aida.histogram1D("Number of Charged Particles", 6, 0, 6);
+
+ pEleVsthetaMoller = aida.histogram2D("2 Electron/P(e) vs Theta Moller", 50, 0, beambeamCut, 50, thetaMin, thetaMax);
+ thetaEleVsthetaMoller = aida.histogram2D("2 Electron/Theta vs Theta Moller", 50, thetaMin, thetaMax, 50, thetaMin, thetaMax);
+ pEleVsthetaBeamBeam = aida.histogram2D("2 Electron/P(e) vs Theta BeamBeam", 50, beambeamCut, eBeam * maxFactor, 50, thetaMin, thetaMax);
+ thetaEleVsthetaBeamBeam = aida.histogram2D("2 Electron/Theta vs Theta BeamBeam", 50, thetaMin, thetaMax, 50, thetaMin, thetaMax);
+
+ mollerMass = aida.histogram1D("2 Electron/Moller Mass (GeV)", 100, 0, 0.100);
+ mollerMassVtxCut = aida.histogram1D("2 Electron/Moller Mass (GeV): VtxCut", 100, 0, 0.100);
+ mollerVx = aida.histogram1D("2 Electron/Moller Vx (mm)", 50, -10, 10);
+ mollerVy = aida.histogram1D("2 Electron/Moller Vy (mm)", 50, -2, 2);
+ mollerVz = aida.histogram1D("2 Electron/Moller Vz (mm)", 50, -50, 50);
+ mollerVzVtxCut = aida.histogram1D("2 Electron/Moller Vz (mm): VtxCut", 50, -50, 50);
+ mollerXVsVtxZ = aida.histogram2D("2 Electron/Moller Vx vs Vz", 100, -5, 5, 100, -50, 50);
+ mollerYVsVtxZ = aida.histogram2D("2 Electron/Moller Vy vs Vz", 100, -2, 2, 100, -50, 50);
+ mollerXVsVtxY = aida.histogram2D("2 Electron/Moller Vx vs Vy", 100, -5, 5, 100, -2, 2);
+
+
+ for (IPlotter plotter : plotters.values()) {
+ //plotter.show();
+ }
+ }
+
+ private void mapMCtoSimTrackerHit(List<SimTrackerHit> simTrackerHits) {
+
+ // Create a relational table that maps SimTrackerHits to MCParticles
+ _hittomc = new BaseRelationalTable<SimTrackerHit, MCParticle>(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+
+ // Loop over the SimTrackerHits and fill in the relational table
+ for (SimTrackerHit simhit : simTrackerHits) {
+ if (simhit.getMCParticle() != null)
+ _hittomc.add(simhit, simhit.getMCParticle());
+ }
+ }
+
+ @Override
+ public void process(EventHeader event){
+ //System.out.println("Pass1");
+ aida.tree().cd("/");
+
+ List<SimTrackerHit> simHits = event.get(SimTrackerHit.class, "TrackerHits");
+ this.mapMCtoSimTrackerHit(simHits);
+
+ RelationalTable hitToStrips = TrackUtils.getHitToStripsTable(event);
+ RelationalTable hitToRotated = TrackUtils.getHitToRotatedTable(event);
+
+ // If the event does not have tracks, skip it
+ //if(!event.hasCollection(Track.class, trackCollectionName)) return;
+
+
+
+ List<ReconstructedParticle> unonstrainedV0List = event.get(ReconstructedParticle.class, unconstrainedV0CandidatesColName);
+ for (ReconstructedParticle uncV0 : unonstrainedV0List) {
+ Vertex uncVert = uncV0.getStartVertex();
+ Hep3Vector pVtxRot = VecOp.mult(beamAxisRotation, uncV0.getMomentum());
+ Hep3Vector vtxPosRot = VecOp.mult(beamAxisRotation, uncVert.getPosition());
+ double theta = Math.acos(pVtxRot.z() / pVtxRot.magnitude());
+ double phi = Math.atan2(pVtxRot.y(), pVtxRot.x());
+ unconVx.fill(vtxPosRot.x());
+ unconVy.fill(vtxPosRot.y());
+ unconVz.fill(vtxPosRot.z());
+ unconMass.fill(uncV0.getMass());
+ unconChi2.fill(uncVert.getChi2());
+ unconVzVsChi2.fill(uncVert.getChi2(), vtxPosRot.z());
+ unconChi2VsTrkChi2.fill(Math.max(uncV0.getParticles().get(0).getTracks().get(0).getChi2(), uncV0.getParticles().get(1).getTracks().get(0).getChi2()), uncVert.getChi2());
+
+ VtxZVsMass.fill(uncV0.getMass(), vtxPosRot.z());
+ VtxXVsVtxZ.fill(vtxPosRot.x(), vtxPosRot.z());
+ VtxYVsVtxZ.fill(vtxPosRot.y(), vtxPosRot.z());
+ VtxXVsVtxY.fill(vtxPosRot.x(), vtxPosRot.y());
+ VtxXVsVtxPx.fill(pVtxRot.x(), vtxPosRot.x());
+ VtxYVsVtxPy.fill(pVtxRot.y(), vtxPosRot.y());
+ VtxZVsVtxPx.fill(pVtxRot.x(), vtxPosRot.z());
+ VtxZVsVtxPy.fill(pVtxRot.y(), vtxPosRot.z());
+ VtxZVsVtxPz.fill(pVtxRot.z(), vtxPosRot.z());
+
+ //this always has 2 tracks.
+ List<ReconstructedParticle> trks = uncV0.getParticles();
+// Track ele = trks.get(0).getTracks().get(0);
+// Track pos = trks.get(1).getTracks().get(0);
+// //if track #0 has charge>0 it's the electron! This seems mixed up, but remember the track
+// //charge is assigned assuming a positive B-field, while ours is negative
+// if (trks.get(0).getCharge() > 0) {
+// pos = trks.get(0).getTracks().get(0);
+// ele = trks.get(1).getTracks().get(0);
+// }
+// aida.histogram2D(plotDir + trkType + triggerType + "/" + "P(e) vs P(p)").fill(getMomentum(ele), getMomentum(pos));
+// aida.histogram2D(plotDir + trkType + triggerType + "/" + "Px(e) vs Px(p)").fill(ele.getTrackStates().get(0).getMomentum()[1], pos.getTrackStates().get(0).getMomentum()[1]);
+// aida.histogram2D(plotDir + trkType + triggerType + "/" + "Py(e) vs Py(p)").fill(ele.getTrackStates().get(0).getMomentum()[2], pos.getTrackStates().get(0).getMomentum()[2]);
+ ReconstructedParticle ele = trks.get(0);
+ ReconstructedParticle pos = trks.get(1);
+ //ReconParticles have the charge correct.
+ if (trks.get(0).getCharge() > 0) {
+ pos = trks.get(0);
+ ele = trks.get(1);
+ }
+ if (ele.getCharge() < 0 && pos.getCharge() > 0) {
+ VtxZVsTrkChi2.fill(Math.max(uncV0.getParticles().get(0).getTracks().get(0).getChi2(), uncV0.getParticles().get(1).getTracks().get(0).getChi2()), uncVert.getPosition().z());
+
+ /*Double[] eleIso = TrackUtils.getIsolations(ele.getTracks().get(0), hitToStrips, hitToRotated);
+ Double[] posIso = TrackUtils.getIsolations(pos.getTracks().get(0), hitToStrips, hitToRotated);
+ if (eleIso[0] != null && posIso[0] != null) {
+ double eleL1Iso = Math.min(Math.abs(eleIso[0]), Math.abs(eleIso[1]));
+ double posL1Iso = Math.min(Math.abs(posIso[0]), Math.abs(posIso[1]));
+ double minL1Iso = Math.min(eleL1Iso, posL1Iso);
+ VtxZVsL1Iso.fill(minL1Iso, uncVert.getPosition().z());
+ }*/
+
+ double pe = ele.getMomentum().magnitude();
+ double pp = pos.getMomentum().magnitude();
+ Hep3Vector pEleRot = VecOp.mult(beamAxisRotation, ele.getMomentum());
+ Hep3Vector pPosRot = VecOp.mult(beamAxisRotation, pos.getMomentum());
+
+ pEleVspPos.fill(pe, pp);
+ pxEleVspxPos.fill(pEleRot.x(), pPosRot.x());
+ pyEleVspyPos.fill(pEleRot.y(), pPosRot.y());
+ if (pe < v0MaxPCut && pp < v0MaxPCut && (pe + pp) > v0ESumMinCut && (pe + pp) < v0ESumMaxCut)//enrich radiative-like events
+
+ pEleVspPosWithCut.fill(pe, pp);
+ }
+
+ double eleT = TrackUtils.getTrackTime(ele.getTracks().get(0), hitToStrips, hitToRotated);
+ double posT = TrackUtils.getTrackTime(pos.getTracks().get(0), hitToStrips, hitToRotated);
+ double meanT = (eleT + posT) / 2.0;
+ v0Time.fill(meanT);
+ v0Dt.fill(eleT - posT);
+ trigTimeV0Time.fill(meanT, event.getTimeStamp() % 24);
+ trigTime.fill(event.getTimeStamp() % 24);
+ }
+
+ List<ReconstructedParticle> beamConstrainedV0List = event.get(ReconstructedParticle.class, beamConV0CandidatesColName);
+ nV0.fill(beamConstrainedV0List.size());
+ for (ReconstructedParticle bsV0 : beamConstrainedV0List) {
+
+ nTotV0++;
+ Vertex bsVert = bsV0.getStartVertex();
+ Hep3Vector vtxPosRot = VecOp.mult(beamAxisRotation, bsVert.getPosition());
+ bsconVx.fill(vtxPosRot.x());
+ bsconVy.fill(vtxPosRot.y());
+ bsconVz.fill(vtxPosRot.z());
+ bsconMass.fill(bsV0.getMass());
+ bsconChi2.fill(bsVert.getChi2());
+ bsconVzVsChi2.fill(bsVert.getChi2(), vtxPosRot.z());
+ bsconChi2VsTrkChi2.fill(Math.max(bsV0.getParticles().get(0).getTracks().get(0).getChi2(), bsV0.getParticles().get(1).getTracks().get(0).getChi2()), bsVert.getChi2());
+ sumMass += bsV0.getMass();
+ sumVx += vtxPosRot.x();
+ sumVy += vtxPosRot.y();
+ sumVz += vtxPosRot.z();
+ sumChi2 += bsVert.getChi2();
+ }
+
+ List<ReconstructedParticle> targetConstrainedV0List = event.get(ReconstructedParticle.class, targetV0ConCandidatesColName);
+ for (ReconstructedParticle tarV0 : targetConstrainedV0List) {
+
+ Vertex tarVert = tarV0.getStartVertex();
+ Hep3Vector vtxPosRot = VecOp.mult(beamAxisRotation, tarVert.getPosition());
+ tarconVx.fill(vtxPosRot.x());
+ tarconVy.fill(vtxPosRot.y());
+ tarconVz.fill(vtxPosRot.z());
+ tarconMass.fill(tarV0.getMass());
+ tarconChi2.fill(tarVert.getChi2());
+ tarconVzVsChi2.fill(tarVert.getChi2(), vtxPosRot.z());
+ tarconChi2VsTrkChi2.fill(Math.max(tarV0.getParticles().get(0).getTracks().get(0).getChi2(), tarV0.getParticles().get(1).getTracks().get(0).getChi2()), tarVert.getChi2());
+ }
+ List<ReconstructedParticle> finalStateParticles = event.get(ReconstructedParticle.class, finalStateParticlesColName);
+
+ ReconstructedParticle ele1 = null;
+ ReconstructedParticle ele2 = null;
+ int sumCharge = 0;
+ int numChargedParticles = 0;
+ for (ReconstructedParticle fsPart : finalStateParticles) {
+ double charge = fsPart.getCharge();
+ sumCharge += charge;
+ if (charge != 0) {
+ numChargedParticles++;
+ if (charge < 1)
+ if (ele1 == null)
+ ele1 = fsPart;
+ else if (!hasSharedStrips(ele1, fsPart, hitToStrips, hitToRotated))
+ ele2 = fsPart;
+ }
+ }
+ sumChargeHisto.fill(sumCharge);
+ numChargeHisto.fill(numChargedParticles);
+
+ if (ele1 != null && ele2 != null) {
+ Hep3Vector p1 = VecOp.mult(beamAxisRotation, ele1.getMomentum());
+ Hep3Vector p2 = VecOp.mult(beamAxisRotation, ele2.getMomentum());
+// Hep3Vector beamAxis = new BasicHep3Vector(Math.sin(0.0305), 0, Math.cos(0.0305));
+// LOGGER.info(p1);
+// LOGGER.info(VecOp.mult(rot, p1));
+
+ double theta1 = Math.acos(p1.z() / p1.magnitude());
+ double theta2 = Math.acos(p2.z() / p2.magnitude());
+ double phi1 = Math.atan2(p1.y(), p1.x());
+ double phi2 = Math.atan2(p2.y(), p2.x());
+ phiEleVsphiEle.fill(phi1, phi2);
+ pEleVspEle.fill(ele1.getMomentum().magnitude(), ele2.getMomentum().magnitude());
+ pyEleVspyEle.fill(ele1.getMomentum().y(), ele2.getMomentum().y());
+ pxEleVspxEle.fill(ele1.getMomentum().x(), ele2.getMomentum().x());
+ //remove beam electrons
+ if (ele1.getMomentum().magnitude() < beambeamCut && ele2.getMomentum().magnitude() < beambeamCut) {
+ pEleVspEleNoBeam.fill(ele1.getMomentum().magnitude(), ele2.getMomentum().magnitude());
+ pyEleVspyEleNoBeam.fill(ele1.getMomentum().y(), ele2.getMomentum().y());
+ pxEleVspxEleNoBeam.fill(ele1.getMomentum().x(), ele2.getMomentum().x());
+ }
+ //look at beam-beam events
+ if (ele1.getMomentum().magnitude() > beambeamCut && ele2.getMomentum().magnitude() > beambeamCut) {
+ pEleVspEleBeamBeam.fill(ele1.getMomentum().magnitude(), ele2.getMomentum().magnitude());
+ pEleVsthetaBeamBeam.fill(p1.magnitude(), theta1);
+ pEleVsthetaBeamBeam.fill(p2.magnitude(), theta2);
+ thetaEleVsthetaBeamBeam.fill(theta1, theta2);
+ }
+
+ //look at "Moller" events (if that's what they really are
+ if (ele1.getMomentum().magnitude() + ele2.getMomentum().magnitude() > molPSumMin
+ && ele1.getMomentum().magnitude() + ele2.getMomentum().magnitude() < molPSumMax
+ && (p1.magnitude() < beambeamCut && p2.magnitude() < beambeamCut)) {
+
+ Track ele1trk = ele1.getTracks().get(0);
+ Track ele2trk = ele2.getTracks().get(0);
+ SeedTrack stEle1 = TrackUtils.makeSeedTrackFromBaseTrack(ele1trk);
+ SeedTrack stEle2 = TrackUtils.makeSeedTrackFromBaseTrack(ele2trk);
+ BilliorTrack btEle1 = new BilliorTrack(stEle1.getSeedCandidate().getHelix());
+ BilliorTrack btEle2 = new BilliorTrack(stEle2.getSeedCandidate().getHelix());
+ BilliorVertex bv = fitVertex(btEle1, btEle2);
+// LOGGER.info("ee vertex: "+bv.toString());
+ mollerMass.fill(bv.getParameters().get("invMass"));
+ mollerVx.fill(bv.getPosition().x());
+ mollerVy.fill(bv.getPosition().y());
+ mollerVz.fill(bv.getPosition().z());
+ mollerXVsVtxZ.fill(bv.getPosition().x(), bv.getPosition().z());
+ mollerYVsVtxZ.fill(bv.getPosition().y(), bv.getPosition().z());
+ mollerXVsVtxY.fill(bv.getPosition().x(), bv.getPosition().y());
+ if (Math.abs(bv.getPosition().x()) < 2
+ && Math.abs(bv.getPosition().y()) < 0.5) {
+ mollerMassVtxCut.fill(bv.getParameters().get("invMass"));
+ mollerVzVtxCut.fill(bv.getPosition().z());
+ }
+ pEleVspEleMoller.fill(p1.magnitude(), p2.magnitude());
+ pEleVsthetaMoller.fill(p1.magnitude(), theta1);
+ pEleVsthetaMoller.fill(p2.magnitude(), theta2);
+ thetaEleVsthetaMoller.fill(theta1, theta2);
+ }
+ }
+
+
+ }
+ @Override
+ public void endOfData(){
+
+
+ System.out.println("%===================================================================%");
+ System.out.println("%====================== Layer 0 Reconctructed Particle Complete==========================%");
+ System.out.println("%===================================================================% \n%");
+ System.out.println("% \n%===================================================================%");
+ }
+
+
+ private BilliorVertex fitVertex(BilliorTrack electron, BilliorTrack positron) {
+ // Create a vertex fitter from the magnetic field.
+ double bField = 0.24;
+ double[] beamSize = {0.001, 0.2, 0.02};
+ BilliorVertexer vtxFitter = new BilliorVertexer(bField);
+ // TODO: The beam size should come from the conditions database.
+ vtxFitter.setBeamSize(beamSize);
+
+ // Perform the vertexing based on the specified constraint.
+ vtxFitter.doBeamSpotConstraint(false);
+
+ // Add the electron and positron tracks to a track list for
+ // the vertex fitter.
+ List<BilliorTrack> billiorTracks = new ArrayList<BilliorTrack>();
+
+ billiorTracks.add(electron);
+
+ billiorTracks.add(positron);
+
+ // Find and return a vertex based on the tracks.
+ return vtxFitter.fitVertex(billiorTracks);
+ }
+
+ private static boolean hasSharedStrips(ReconstructedParticle vertex, RelationalTable hittostrip, RelationalTable hittorotated) {
+ return hasSharedStrips(vertex.getParticles().get(0), vertex.getParticles().get(1), hittostrip, hittorotated);
+ }
+
+ private static boolean hasSharedStrips(ReconstructedParticle fs1, ReconstructedParticle fs2, RelationalTable hittostrip, RelationalTable hittorotated) {
+ return TrackUtils.hasSharedStrips(fs1.getTracks().get(0), fs2.getTracks().get(0), hittostrip, hittorotated);
+ }
+
+}
Added: java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/TrackingReconstructionPlots.java
=============================================================================
--- java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/TrackingReconstructionPlots.java (added)
+++ java/branches/layer0-thin-branch/users/src/main/java/org/hps/users/mrsolt/TrackingReconstructionPlots.java Wed Nov 30 11:12:55 2016
@@ -0,0 +1,1563 @@
+package org.hps.users.mrsolt;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.aida.IPlotter;
+import hep.aida.IPlotterStyle;
+import hep.aida.IProfile;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+
+import java.io.IOException;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+import org.hps.recon.tracking.BeamlineConstants;
+import org.hps.recon.tracking.DumbShaperFit;
+import org.hps.recon.tracking.FittedRawTrackerHit;
+import org.hps.recon.tracking.HelixConverter;
+import org.hps.recon.tracking.ShaperFitAlgorithm;
+import org.hps.recon.tracking.StraightLineTrack;
+import org.hps.recon.tracking.TrackUtils;
+import org.hps.recon.tracking.gbl.HelicalTrackStripGbl;
+import org.lcsim.detector.tracker.silicon.HpsSiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCIOParameters.ParameterName;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+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.HelixUtils;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.geometry.compact.converter.HPSTrackerBuilder;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+import org.lcsim.recon.tracking.digitization.sisim.TrackerHitType;
+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;
+
+/**
+ * Analysis class to check recon.
+ *
+ * @author phansson
+ */
+public class TrackingReconstructionPlots extends Driver {
+
+ static {
+ hep.aida.jfree.AnalysisFactory.register();
+ }
+
+
+ private AIDA aida = AIDA.defaultInstance();
+ private String helicalTrackHitCollectionName = "HelicalTrackHits";
+ private String stripClusterCollectionName = "StripClusterer_SiTrackerHitStrip1D";
+
+ private String trackCollectionName = "MatchedTracks";
+ String ecalSubdetectorName = "Ecal";
+ String ecalCollectionName = "EcalClusters";
+ IDDecoder dec;
+ private String outputPlots = null;
+ IPlotter plotter;
+ IPlotter plotter2;
+ IPlotter plotter22;
+ IPlotter plotter2221;
+ IPlotter plotter2222;
+ IPlotter plotter222;
+ IPlotter plotter22299;
+ IPlotter plotter22298;
+ IPlotter plotter2224;
+ IPlotter plotter2223;
+ IPlotter plotter3;
+ IPlotter plotter3_1;
+ IPlotter plotter3_11;
+ IPlotter plotter3_2;
+ IPlotter plotter4;
+ IPlotter plotter5;
+ IPlotter plotter5_1;
+ IPlotter plotter55;
+ IPlotter plotter6;
+ IPlotter plotter66;
+ IPlotter plotter8;
+ IPlotter plotter88;
+ IPlotter plotter888;
+ IPlotter plotter8888;
+ IPlotter top1;
+ IPlotter top2;
+ IPlotter top3;
+ IPlotter top4;
+ IPlotter top44;
+ IPlotter bot1;
+ IPlotter bot2;
+ IPlotter bot3;
+ IPlotter bot4;
+ double zAtColl = -1500;
+ IHistogram1D trkPx;
+ IHistogram1D nTracks;
+ IHistogram1D nTracksTop;
+ IHistogram1D nTracksBot;
+ ShaperFitAlgorithm _shaper = new DumbShaperFit();
+ HelixConverter converter = new HelixConverter(0);
+ private boolean showPlots = true;
+ private double _bfield;
+ private static Logger LOGGER = Logger.getLogger(TrackingReconstructionPlots.class.getName());
+ private List<HpsSiSensor> sensors = new ArrayList<HpsSiSensor>();
+ private int nLayers = 7;
+
+ @Override
+ protected void detectorChanged(Detector detector) {
+ aida.tree().cd("/");
+ for(HpsSiSensor s : detector.getDetectorElement().findDescendants(HpsSiSensor.class)) {
+ if(s.getName().startsWith("module_") && s.getName().endsWith("sensor0")) {
+ sensors.add(s);
+ }
+ }
+ LOGGER.info("Found " + sensors.size() + " SiSensors.");
+
+ Hep3Vector bfieldvec = detector.getFieldMap().getField(new BasicHep3Vector(0., 0., 1.));
+ _bfield = bfieldvec.y();
+
+ setupPlots();
+ }
+
+
+
+
+
+
+
+ public TrackingReconstructionPlots() {
+ LOGGER.setLevel(Level.WARNING);
+ }
+
+ public void setOutputPlots(String output) {
+ this.outputPlots = output;
+ }
+
+ public void setShowPlots(boolean show) {
+ this.showPlots = show;
+ }
+
+ public void setHelicalTrackHitCollectionName(String helicalTrackHitCollectionName) {
+ this.helicalTrackHitCollectionName = helicalTrackHitCollectionName;
+ }
+
+ public void setTrackCollectionName(String trackCollectionName) {
+ this.trackCollectionName = trackCollectionName;
+ }
+
+ @Override
+ public void process(EventHeader event) {
+ aida.tree().cd("/");
+ if (!event.hasCollection(HelicalTrackHit.class, helicalTrackHitCollectionName)) {
+// System.out.println(helicalTrackHitCollectionName + " does not exist; skipping event");
+ return;
+ }
+
+ List<SiTrackerHitStrip1D> stripClusters = event.get(SiTrackerHitStrip1D.class, stripClusterCollectionName);
+ //System.out.printf("%s: Got %d SiTrackerHitStrip1D in this event\n", stripHits.size());
+ Map<HpsSiSensor, Integer> stripHits = new HashMap<HpsSiSensor, Integer>();
+ Map<HpsSiSensor, Double> stripHitsIso = new HashMap<HpsSiSensor, Double>();
+ for (SiTrackerHitStrip1D stripHit : stripClusters) {
+ HpsSiSensor sensor = (HpsSiSensor) stripHit.getRawHits().get(0).getDetectorElement();
+ int n;
+ if(stripHits.containsKey(sensor)) {
+ n = stripHits.get(sensor);
+ } else {
+ n=0;
+ }
+ n++;
+ stripHits.put(sensor, n);
+
+ // calculate isolation to other strip clusters
+
+ SiTrackerHitStrip1D local = stripHit.getTransformedHit(TrackerHitType.CoordinateSystem.SENSOR);
+
+ double stripIsoMin = 9999.9;
+ for (SiTrackerHitStrip1D stripHitOther : stripClusters) {
+ LOGGER.fine(stripHit.getPositionAsVector().toString() + " c.f. " + stripHitOther.getPositionAsVector().toString());
+
+ if(stripHitOther.equals(stripHit)) {
+ continue;
+ }
+
+ HpsSiSensor sensorOther = (HpsSiSensor) stripHitOther.getRawHits().get(0).getDetectorElement();
+ //System.out.println(sensor.getName() + " c.f. " + sensorOther.getName());
+ if(sensorOther.equals(sensor)) {
+ SiTrackerHitStrip1D localOther = stripHitOther.getTransformedHit(TrackerHitType.CoordinateSystem.SENSOR);
+ double d = Math.abs(local.getPosition()[0] - localOther.getPosition()[0]);
+ //System.out.println(sensor.getName() + " d " + Double.toString(d));
+ if (d < stripIsoMin && d > 0) {
+ stripIsoMin = d;
+ }
+ }
+ }
+ stripHitsIso.put(sensor, stripIsoMin);
+ }
+
+ for(Map.Entry<HpsSiSensor,Integer> sensor : stripHits.entrySet()) {
+ aida.histogram1D(sensor.getKey().getName() + " strip hits").fill(sensor.getValue());
+ aida.histogram1D(sensor.getKey().getName() + " strip hits iso").fill(stripHitsIso.get(sensor.getKey()));
+ }
+
+
+ List<HelicalTrackHit> hthList = event.get(HelicalTrackHit.class, helicalTrackHitCollectionName);
+ Map<Integer, Integer> layersTop = new HashMap<Integer,Integer>();
+ Map<Integer, Integer> layersBot = new HashMap<Integer,Integer>();
+ Map<HpsSiSensor, Integer> stripHitsFromStereoHits = new HashMap<HpsSiSensor, Integer>();
+ for (HelicalTrackHit hth : hthList) {
+ HelicalTrackCross htc = (HelicalTrackCross) hth;
+ HpsSiSensor sensor = ((HpsSiSensor) ((RawTrackerHit) htc.getRawHits().get(0)).getDetectorElement());
+ for(HelicalTrackStrip strip : htc.getStrips()) {
+ HpsSiSensor stripsensor = (HpsSiSensor) ((RawTrackerHit)strip.rawhits().get(0)).getDetectorElement();
+ if(stripHitsFromStereoHits.containsKey(stripsensor)) {
+ stripHitsFromStereoHits.put(stripsensor, stripHitsFromStereoHits.get(stripsensor) + 1);
+ } else {
+ stripHitsFromStereoHits.put(stripsensor, 0);
+ }
+ }
+ int l = sensor.getLayerNumberFromDetectorElement();
+ int n = 0;
+ if(sensor.isTopLayer()) {
+ if (layersTop.containsKey(l) )
+ n = layersTop.get(l);
+ else
+ layersTop.put(l,n+1);
+ }
+ else {
+ if (layersBot.containsKey(l) )
+ n = layersTop.get(l);
+ else
+ layersTop.put(l,n+1);
+ }
+ }
+
+
+ for(Map.Entry<HpsSiSensor,Integer> sensor : stripHitsFromStereoHits.entrySet()) {
+ aida.histogram1D(sensor.getKey().getName() + " strip hits from stereo").fill(sensor.getValue());
+ }
+
+ for (Map.Entry<Integer,Integer> e : layersTop.entrySet())
+ aida.profile1D("Number of Stereo Hits per layer in Top Half").fill(e.getKey(), e.getValue());
+
+ for (Map.Entry<Integer,Integer> e : layersBot.entrySet())
+ aida.profile1D("Number of Stereo Hits per layer in Bottom Half").fill(e.getKey(), e.getValue());
+
+
+
+ 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());
+ int nBotClusters = 0;
+ int nTopClusters = 0;
+ if (event.hasCollection(Cluster.class, ecalCollectionName)) {
+ List<Cluster> clusters = event.get(Cluster.class, ecalCollectionName);
+ for (Cluster cluster : clusters) {
+ // Get the ix and iy indices for the seed.
+// final int ix = cluster.getCalorimeterHits().get(0).getIdentifierFieldValue("ix");
+// final int iy = cluster.getCalorimeterHits().get(0).getIdentifierFieldValue("iy");
+
+ //System.out.println("cluser position = ("+cluster.getPosition()[0]+","+cluster.getPosition()[1]+") with energy = "+cluster.getEnergy());
+ if (cluster.getPosition()[1] > 0) {
+ nTopClusters++;
+ //System.out.println("cl " + cluster.getPosition()[0] + " " + cluster.getPosition()[1] + " ix " + ix + " iy " + iy);
+ aida.histogram2D("Top ECal Cluster Position").fill(cluster.getPosition()[0], cluster.getPosition()[1]);
+ aida.histogram1D("Top ECal Cluster Energy").fill(cluster.getEnergy());
+ }
+ if (cluster.getPosition()[1] < 0) {
+ nBotClusters++;
+ aida.histogram2D("Bottom ECal Cluster Position").fill(cluster.getPosition()[0], cluster.getPosition()[1]);
+ aida.histogram1D("Bottom ECal Cluster Energy").fill(cluster.getEnergy());
+ }
+
+ if (tracks.size() > 0) {
+ if (cluster.getPosition()[1] > 0) {
+ aida.histogram2D("Top ECal Cluster Position (>0 tracks)").fill(cluster.getPosition()[0], cluster.getPosition()[1]);
+ }
+ if (cluster.getPosition()[1] < 0) {
+ aida.histogram2D("Bottom ECal Cluster Position (>0 tracks)").fill(cluster.getPosition()[0], cluster.getPosition()[1]);
+ }
+
+ if (cluster.getEnergy() > 0.1) {
+ if (cluster.getPosition()[1] > 0) {
+ aida.histogram2D("Top ECal Cluster Position (E>0.1,>0 tracks)").fill(cluster.getPosition()[0], cluster.getPosition()[1]);
+ aida.histogram2D("Top ECal Cluster Position w_E (E>0.1,>0 tracks)").fill(cluster.getPosition()[0], cluster.getPosition()[1], cluster.getEnergy());
+ }
+ if (cluster.getPosition()[1] < 0) {
+ aida.histogram2D("Bottom ECal Cluster Position (E>0.1,>0 tracks)").fill(cluster.getPosition()[0], cluster.getPosition()[1]);
+ aida.histogram2D("Bottom ECal Cluster Position w_E (E>0.1,>0 tracks)").fill(cluster.getPosition()[0], cluster.getPosition()[1], cluster.getEnergy());
+ }
+ }
+ }
+
+ }
+ }
+
+ aida.histogram1D("Number of Clusters Top").fill(nTopClusters);
+ aida.histogram1D("Number of Clusters Bot").fill(nBotClusters);
+
+
+
+
+
+ //List<SiTrackerHitStrip1D> stripHits = event.get(SiTrackerHitStrip1D.class, "StripClusterer_SiTrackerHitStrip1D");
+ //int stripClustersPerLayerTop[] = getStripClustersPerLayer(stripHits, "up");
+ //int stripClustersPerLayerBottom[] = getStripClustersPerLayer(stripHits,"down");
+
+ //boolean hasSingleStripClusterPerLayer = singleStripClusterPerLayer(stripClustersPerLayerTop);
+
+ Map<Track, Cluster> eCanditates = new HashMap<Track, Cluster>();
+ Map<Track, Cluster> pCanditates = new HashMap<Track, Cluster>();
+
+ int ntracksTop = 0;
+ int ntracksBot = 0;
+
+ for (Track trk : tracks) {
+
+ //boolean isSingleHitPerLayerTrack = singleTrackHitPerLayer(trk);
+
+ aida.histogram1D("Track Momentum (Px)").fill(trk.getTrackStates().get(0).getMomentum()[1]);
+ aida.histogram1D("Track Momentum (Py)").fill(trk.getTrackStates().get(0).getMomentum()[2]);
+ aida.histogram1D("Track Momentum (Pz)").fill(trk.getTrackStates().get(0).getMomentum()[0]);
+ aida.histogram1D("Track Chi2").fill(trk.getChi2());
+
+ aida.histogram1D("Hits per Track").fill(trk.getTrackerHits().size());
+ SeedTrack stEle = (SeedTrack) trk;
+ SeedCandidate seedCandidate = stEle.getSeedCandidate();
+ HelicalTrackFit helicalTrackFit = seedCandidate.getHelix();
+ StraightLineTrack slt = converter.Convert(helicalTrackFit);
+
+ Hep3Vector posAtEcal = TrackUtils.getTrackPositionAtEcal(trk);
+
+ aida.histogram1D("X (mm) @ Z=-60cm").fill(slt.getYZAtX(BeamlineConstants.HARP_POSITION_TESTRUN)[0]); //this is y in the tracker frame
+ aida.histogram1D("Y (mm) @ Z=-60cm").fill(slt.getYZAtX(BeamlineConstants.HARP_POSITION_TESTRUN)[1]); //this is z in the tracker frame
+ aida.histogram1D("X (mm) @ Z=-150cm").fill(slt.getYZAtX(zAtColl)[0]);
+ aida.histogram1D("Y (mm) @ Z=-150cm").fill(slt.getYZAtX(zAtColl)[1]);
+
+ aida.histogram1D("X (mm) @ ECAL").fill(posAtEcal.x());
+ aida.histogram1D("Y (mm) @ ECAL").fill(posAtEcal.y());
+ if (trk.getTrackStates().get(0).getMomentum()[0] > 1.0) {
+ aida.histogram1D("X (mm) @ ECAL (Pz>1)").fill(posAtEcal.x());
+ aida.histogram1D("Y (mm) @ ECAL (Pz>1)").fill(posAtEcal.y());
+ }
+ aida.histogram1D("d0 ").fill(trk.getTrackStates().get(0).getParameter(ParameterName.d0.ordinal()));
+ aida.histogram1D("sinphi ").fill(Math.sin(trk.getTrackStates().get(0).getParameter(ParameterName.phi0.ordinal())));
+ aida.histogram1D("omega ").fill(trk.getTrackStates().get(0).getParameter(ParameterName.omega.ordinal()));
+ aida.histogram1D("tan(lambda) ").fill(trk.getTrackStates().get(0).getParameter(ParameterName.tanLambda.ordinal()));
+ aida.histogram1D("z0 ").fill(trk.getTrackStates().get(0).getParameter(ParameterName.z0.ordinal()));
+
+ int isTop = -1;
+ if (trk.getTrackerHits().get(0).getPosition()[2] > 0) {
+ isTop = 0;
+ }
+ if (isTop == 0) {
+ aida.histogram1D("Top Track Momentum (Px)").fill(trk.getTrackStates().get(0).getMomentum()[1]);
+ aida.histogram1D("Top Track Momentum (Py)").fill(trk.getTrackStates().get(0).getMomentum()[2]);
+ aida.histogram1D("Top Track Momentum (Pz)").fill(trk.getTrackStates().get(0).getMomentum()[0]);
+ aida.histogram1D("Top Track Chi2").fill(trk.getChi2());
+
+ aida.histogram1D("d0 Top").fill(trk.getTrackStates().get(0).getParameter(ParameterName.d0.ordinal()));
+ aida.histogram1D("sinphi Top").fill(Math.sin(trk.getTrackStates().get(0).getParameter(ParameterName.phi0.ordinal())));
+ aida.histogram1D("omega Top").fill(trk.getTrackStates().get(0).getParameter(ParameterName.omega.ordinal()));
+ aida.histogram1D("tan(lambda) Top").fill(trk.getTrackStates().get(0).getParameter(ParameterName.tanLambda.ordinal()));
+ aida.histogram1D("z0 Top").fill(trk.getTrackStates().get(0).getParameter(ParameterName.z0.ordinal()));
+ ntracksTop++;
+ } else {
+ aida.histogram1D("Bottom Track Momentum (Px)").fill(trk.getTrackStates().get(0).getMomentum()[1]);
+ aida.histogram1D("Bottom Track Momentum (Py)").fill(trk.getTrackStates().get(0).getMomentum()[2]);
+ aida.histogram1D("Bottom Track Momentum (Pz)").fill(trk.getTrackStates().get(0).getMomentum()[0]);
+ aida.histogram1D("Bottom Track Chi2").fill(trk.getChi2());
+
+ aida.histogram1D("d0 Bottom").fill(trk.getTrackStates().get(0).getParameter(ParameterName.d0.ordinal()));
+ aida.histogram1D("sinphi Bottom").fill(Math.sin(trk.getTrackStates().get(0).getParameter(ParameterName.phi0.ordinal())));
+ aida.histogram1D("omega Bottom").fill(trk.getTrackStates().get(0).getParameter(ParameterName.omega.ordinal()));
+ aida.histogram1D("tan(lambda) Bottom").fill(trk.getTrackStates().get(0).getParameter(ParameterName.tanLambda.ordinal()));
+ aida.histogram1D("z0 Bottom").fill(trk.getTrackStates().get(0).getParameter(ParameterName.z0.ordinal()));
+ ntracksBot++;
+ }
+ List<TrackerHit> hitsOnTrack = trk.getTrackerHits();
+ Map<HpsSiSensor, Integer> stripHitsOnTrack = new HashMap<HpsSiSensor, Integer>();
+ Map<HpsSiSensor, Double> stripHitsIsoOnTrack = new HashMap<HpsSiSensor, Double>();
+
+ for (TrackerHit hit : hitsOnTrack) {
+
+ HelicalTrackHit htc = (HelicalTrackHit) hit;
+ HelicalTrackCross htcross = (HelicalTrackCross) htc;
+ double sHit = helicalTrackFit.PathMap().get(htc);
+ Hep3Vector posonhelix = HelixUtils.PointOnHelix(helicalTrackFit, sHit);
+ int layer = ((HpsSiSensor) ((RawTrackerHit) htcross.getStrips().get(0).rawhits().get(0)).getDetectorElement()).getLayerNumberFromDetectorElement();
+ boolean isTopLayer = !((HpsSiSensor) ((RawTrackerHit) htcross.getStrips().get(0).rawhits().get(0)).getDetectorElement()).isBottomLayer();
+
+
+
+ //HpsSiSensor sensor = ((HpsSiSensor) ((RawTrackerHit) htc.getRawHits().get(0)).getDetectorElement());
+ for(HelicalTrackStrip strip : htcross.getStrips()) {
+ HpsSiSensor sensor = ((HpsSiSensor) ((RawTrackerHit) strip.rawhits().get(0)).getDetectorElement());
+ if(sensor.isTopLayer()) isTopLayer = true;
+ else isTopLayer=false;
+ HelicalTrackStripGbl stripGbl = new HelicalTrackStripGbl(strip, true);
+ Map<String, Double> stripResiduals = TrackUtils.calculateLocalTrackHitResiduals(helicalTrackFit, stripGbl, 0.,0.,_bfield);
+ LOGGER.fine("Sensor " + sensor.getName() + " ures = " + stripResiduals.get("ures"));
+ aida.histogram1D(sensor.getName() + " strip residual (mm)").fill(stripResiduals.get("ures"));
+
+
+ // calculate isolation to other strip clusters
+ double stripIsoMin = 9999.9;
+ for (SiTrackerHitStrip1D stripHit : stripClusters) {
+ if (stripHit.getRawHits().get(0).getDetectorElement().getName().equals(sensor.getName())) {
+ SiTrackerHitStrip1D local = stripHit.getTransformedHit(TrackerHitType.CoordinateSystem.SENSOR);
+ double d = Math.abs(strip.umeas() - local.getPosition()[0]);
+ if (d < stripIsoMin && d > 0) {
+ stripIsoMin = d;
+ }
+ }
+ }
+
+ if(stripHitsOnTrack.containsKey(sensor)) {
+ stripHitsOnTrack.put(sensor, stripHitsOnTrack.get(sensor) + 1);
+ } else {
+ stripHitsOnTrack.put(sensor, 1);
+ }
+ stripHitsIsoOnTrack.put(sensor, stripIsoMin);
+ }
+
+
+
+
+
+ double yTr = posonhelix.y();
+ double zTr = posonhelix.z();
+ String modNum = "Layer " + String.valueOf(layer) + " ";
+
+ aida.histogram1D(modNum + "Residual X(mm)").fill(htcross.getCorrectedPosition().y() - yTr);
+ aida.histogram1D(modNum + "Residual Y(mm)").fill(htcross.getCorrectedPosition().z() - zTr);
+
+ if (isTopLayer) {
+ aida.histogram1D(modNum + "Residual X(mm) Top").fill(htcross.getCorrectedPosition().y() - yTr);
+ aida.histogram1D(modNum + "Residual Y(mm) Top").fill(htcross.getCorrectedPosition().z() - zTr);
+
+ } else {
+ aida.histogram1D(modNum + "Residual X(mm) Bottom").fill(htcross.getCorrectedPosition().y() - yTr);
+ aida.histogram1D(modNum + "Residual Y(mm) Bottom").fill(htcross.getCorrectedPosition().z() - zTr);
+ }
+
+
+ boolean doAmplitudePlots = true;
+ if(doAmplitudePlots) {
+ for (HelicalTrackStrip hts : htcross.getStrips()) {
+ double clusterSum = 0;
+ double clusterT0 = 0;
+ int nHitsCluster = 0;
+
+ for (RawTrackerHit rawHit : (List<RawTrackerHit>) hts.rawhits()) {
+ if(event.hasCollection(LCRelation.class, "SVTFittedRawTrackerHits")) {
+ List<LCRelation> fittedHits = event.get(LCRelation.class, "SVTFittedRawTrackerHits");
+ for(LCRelation fittedHit : fittedHits) {
+ if(rawHit.equals((RawTrackerHit)fittedHit.getFrom())) {
+ double amp = FittedRawTrackerHit.getAmp(fittedHit);
+ double t0 = FittedRawTrackerHit.getT0(fittedHit);
+ //System.out.println("to="+t0 + " amp=" + amp);
+ aida.histogram1D("Amp (HitOnTrack)").fill(amp);
+ if (trk.getTrackStates().get(0).getMomentum()[0] > 0.8) {
+ aida.histogram1D("Amp Pz>0.8 (HitOnTrack)").fill(amp);
+ }
+ aida.histogram1D("t0 (HitOnTrack)").fill(t0);
+ if (trk.getTrackStates().get(0).getMomentum()[0] > 0.8) {
+ aida.histogram1D("t0 Pz>0.8 (HitOnTrack)").fill(t0);
+ }
+ clusterSum += amp;
+ clusterT0 += t0;
+ nHitsCluster++;
+ }
+ }
+ }
+ }
+
+ aida.histogram1D("Hits in Cluster (HitOnTrack)").fill(nHitsCluster);
+ aida.histogram1D("Cluster Amp (HitOnTrack)").fill(clusterSum);
+ if (trk.getTrackStates().get(0).getMomentum()[0] > 0.8) {
+ aida.histogram1D("Cluster Amp Pz>0.8 (HitOnTrack)").fill(clusterSum);
+ }
+ if(nHitsCluster>0) {
+ aida.histogram1D("Cluster t0 (HitOnTrack)").fill(clusterT0/nHitsCluster);
+ if (trk.getTrackStates().get(0).getMomentum()[0] > 0.8) {
+ aida.histogram1D("Cluster t0 Pz>0.8 (HitOnTrack)").fill(clusterT0/nHitsCluster);
+ }
+ }
+
+ }
+ }
+ }
+
+ for(Map.Entry<HpsSiSensor,Integer> sensor : stripHitsOnTrack.entrySet())
+ aida.histogram1D(sensor.getKey().getName() + " strip hits iso on track").fill(stripHitsIsoOnTrack.get(sensor.getKey()));
+
+
+
+ Cluster clust = null;
+ if(event.hasCollection(Cluster.class,ecalCollectionName)) {
+ List<Cluster> clusters = event.get(Cluster.class, ecalCollectionName);
+
+ Cluster cand_clust = findClosestCluster(posAtEcal, clusters);
+
+ if (cand_clust != null) {
+
+ // track matching requirement
+ if(Math.abs( posAtEcal.x() - cand_clust.getPosition()[0])<30.0 &&
+ Math.abs( posAtEcal.y() - cand_clust.getPosition()[1])<30.0)
+ {
+ clust = cand_clust;
+ if(trk.getCharge()<0) pCanditates.put(trk, clust);
+ else eCanditates.put(trk, clust);
+
+ posAtEcal = TrackUtils.extrapolateTrack(trk, clust.getPosition()[2]);//.positionAtEcal();
+
+ aida.histogram2D("Energy Vs Momentum").fill(clust.getEnergy(), trk.getTrackStates().get(0).getMomentum()[0]);
+ aida.histogram1D("Energy Over Momentum").fill(clust.getEnergy() / (trk.getTrackStates().get(0).getMomentum()[0] ));
+ aida.histogram1D("deltaX").fill(clust.getPosition()[0] - posAtEcal.x());
+ aida.histogram1D("deltaY").fill(clust.getPosition()[1] - posAtEcal.y());
+ aida.histogram2D("X ECal Vs Track").fill(clust.getPosition()[0], posAtEcal.x());
+ aida.histogram2D("Y ECal Vs Track").fill(clust.getPosition()[1], posAtEcal.y());
+
+ if (isTop == 0) {
+ aida.histogram2D("Top Energy Vs Momentum").fill(clust.getEnergy(), trk.getTrackStates().get(0).getMomentum()[0] );
+ // aida.histogram2D("Top Energy Vs Momentum").fill(posAtEcal.y(), trk.getTrackStates().get(0).getMomentum()[0]);
+ aida.histogram1D("Top Energy Over Momentum").fill(clust.getEnergy() / (trk.getTrackStates().get(0).getMomentum()[0]));
+ aida.histogram1D("Top deltaX").fill(clust.getPosition()[0] - posAtEcal.x());
+ aida.histogram1D("Top deltaY").fill(clust.getPosition()[1] - posAtEcal.y());
+ aida.histogram2D("Top deltaX vs X").fill(clust.getPosition()[0], clust.getPosition()[0] - posAtEcal.x());
+ aida.histogram2D("Top deltaY vs Y").fill(clust.getPosition()[1], clust.getPosition()[1] - posAtEcal.y());
+ aida.histogram2D("Top X ECal Vs Track").fill(clust.getPosition()[0], posAtEcal.x());
+ aida.histogram2D("Top Y ECal Vs Track").fill(clust.getPosition()[1], posAtEcal.y());
+ } else {
+ aida.histogram2D("Bottom Energy Vs Momentum").fill(clust.getEnergy(), trk.getTrackStates().get(0).getMomentum()[0] );
+ aida.histogram1D("Bottom Energy Over Momentum").fill(clust.getEnergy() / (trk.getTrackStates().get(0).getMomentum()[0]));
+ aida.histogram1D("Bottom deltaX").fill(clust.getPosition()[0] - posAtEcal.x());
+ aida.histogram1D("Bottom deltaY").fill(clust.getPosition()[1] - posAtEcal.y());
+ aida.histogram2D("Bottom deltaX vs X").fill(clust.getPosition()[0], clust.getPosition()[0] - posAtEcal.x());
+ aida.histogram2D("Bottom deltaY vs Y").fill(clust.getPosition()[1], clust.getPosition()[1] - posAtEcal.y());
+ aida.histogram2D("Bottom X ECal Vs Track").fill(clust.getPosition()[0], posAtEcal.x());
+ aida.histogram2D("Bottom Y ECal Vs Track").fill(clust.getPosition()[1], posAtEcal.y());
+ }
+ }
+ }
+ }
+
+ if (clust == null) {
+ aida.histogram1D("Tracks matched").fill(0);
+ if(trk.getTrackStates().get(0).getMomentum()[0] > 0.8){
+ aida.histogram1D("Tracks matched (Pz>0.8)").fill(0);
+ }
+ if(isTop == 0) {
+ aida.histogram1D("Tracks matched Top").fill(0);
+ if(trk.getTrackStates().get(0).getMomentum()[0] > 0.8){
+ aida.histogram1D("Tracks matched Top (Pz>0.8)").fill(0);
+ }
+ } else {
+ aida.histogram1D("Tracks matched Bottom").fill(0);
+ if(trk.getTrackStates().get(0).getMomentum()[0] > 0.8){
+ aida.histogram1D("Tracks matched Bottom (Pz>0.8)").fill(0);
+ }
+ }
+ } else {
+ aida.histogram1D("Tracks matched").fill(1);
+ if(trk.getTrackStates().get(0).getMomentum()[0] > 0.8){
+ aida.histogram1D("Tracks matched (Pz>0.8)").fill(1);
+ }
+
+ if (isTop == 0) {
+ aida.histogram1D("Tracks matched Top").fill(1);
+ if(trk.getTrackStates().get(0).getMomentum()[0] > 0.8){
+ aida.histogram1D("Tracks matched Top (Pz>0.8)").fill(1);
+ }
+ } else {
+ aida.histogram1D("Tracks matched Bottom").fill(1);
+ if(trk.getTrackStates().get(0).getMomentum()[0] > 0.8){
+ aida.histogram1D("Tracks matched Bottom (Pz>0.8)").fill(1);
+ }
+ }
+
+ }
+
+ }
+
+ nTracksBot.fill(ntracksBot);
+ nTracksTop.fill(ntracksTop);
+
+ // e+/e-
+ Map.Entry<Track, Cluster> ecand_highestP = null;
+ double e_pmax = -1;
+ Map.Entry<Track, Cluster> pcand_highestP = null;
+ double p_pmax = -1;
+ for(Map.Entry<Track, Cluster> ecand : eCanditates.entrySet()) {
+ double p = getMomentum(ecand.getKey());
+ aida.histogram1D("p(e-)").fill(p);
+ if(ecand_highestP == null ) {
+ ecand_highestP = ecand;
+ e_pmax = getMomentum(ecand_highestP.getKey());
+ } else {
+ if(p > e_pmax) {
+ ecand_highestP = ecand;
+ e_pmax = getMomentum(ecand_highestP.getKey());
+ }
+ }
+ }
+
+ for(Map.Entry<Track, Cluster> pcand : pCanditates.entrySet()) {
+ double p = getMomentum(pcand.getKey());
+ aida.histogram1D("p(e+)").fill(p);
+ if(pcand_highestP == null ) {
+ pcand_highestP = pcand;
+ p_pmax = getMomentum(pcand_highestP.getKey());
+ } else {
+ if(p > p_pmax) {
+ pcand_highestP = pcand;
+ p_pmax = getMomentum(pcand_highestP.getKey());
+ }
+ }
+ }
+
+ aida.histogram1D("n(e-)").fill(eCanditates.size());
+ aida.histogram1D("n(e+)").fill(pCanditates.size());
+ if(ecand_highestP!=null) {
+ aida.histogram1D("p(e-) max").fill(e_pmax);
+ }
+ if(pcand_highestP!=null) {
+ aida.histogram1D("p(e+) max").fill(p_pmax);
+ }
+ if(ecand_highestP!=null && pcand_highestP!=null) {
+ aida.histogram2D("p(e-) vs p(e+) max").fill(e_pmax, p_pmax);
+ }
+
+
+ }
+
+ private double getMomentum(Track trk) {
+ double p = Math.sqrt(trk.getTrackStates().get(0).getMomentum()[0]*trk.getTrackStates().get(0).getMomentum()[0] +
+ trk.getTrackStates().get(0).getMomentum()[1]*trk.getTrackStates().get(0).getMomentum()[1] +
+ trk.getTrackStates().get(0).getMomentum()[2]*trk.getTrackStates().get(0).getMomentum()[2]);
+ return p;
+ }
+
+ public int[] getTrackHitsPerLayer(Track trk) {
+ int n[] = {0, 0, 0, 0, 0, 0};
+ List<TrackerHit> hitsOnTrack = trk.getTrackerHits();
+ int layer;
+ for (TrackerHit hit : hitsOnTrack) {
+ HelicalTrackHit htc = (HelicalTrackHit) hit;
+// if (htc.getPosition()[2] < 0) {
+ layer = htc.Layer();
+ layer = (layer - 1) / 2;
+ n[layer] = n[layer] + 1;
+// }
+ }
+
+ return n;
+ }
+
+ public boolean singleTrackHitPerLayer(Track track) {
+ int hitsPerLayer[] = getTrackHitsPerLayer(track);
+ for (int i = 0; i < 5; ++i) {
+ if (hitsPerLayer[i] != 1) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ public boolean singleStripClusterPerLayer(int hitsPerLayer[]) {
+ //This includes both axial and stereo separately
+ // so for a hit in each double layer we need 10 hits
+ for (int i = 0; i < 10; ++i) {
+ if (hitsPerLayer[i] != 1) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ public int[] getStripClustersPerLayer(List<SiTrackerHitStrip1D> trackerHits, String side) {
+ String name;
+ int l;
+ int n[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ boolean ddd = false;
+
+ if (ddd) {
+ System.out.println("Get # hits per layer on side \"" + side + "\"");
+ }
+
+ for (SiTrackerHitStrip1D stripCluster : trackerHits) {
+
+ if (ddd) {
+ System.out.println("Processing stripCluster " + stripCluster.toString());
+ }
+
+ if (!"".equals(side)) {
+ String s;
+ if (stripCluster.getPosition()[1] >= 0.0) {
+ s = "up";
+ } else {
+ s = "down";
+ }
+ if (!s.equals(side)) {
+ continue;
+ }
+ }
+
+ name = stripCluster.getSensor().getName();
+ if (name.length() < 14) {
+ System.err.println("This name is too short!!");
+ throw new RuntimeException("SiSensor name " + name + " is invalid?");
+ }
+
+ if (ddd) {
+ System.out.println("sensor name " + name);
+ }
+
+ if(name.contains("layer") && name.contains("_module")) {
+ //String str_l = name.substring(13);
+ String str_l = name.substring(name.indexOf("layer") + 5, name.indexOf("_module"));
+ l = Integer.parseInt(str_l);
+ }
+ else if(name.contains("module") && name.contains("_halfmodule")) {
+ int ll = HPSTrackerBuilder.getLayerFromVolumeName(name);
+ boolean isAxial = HPSTrackerBuilder.isAxialFromName(name);
+ boolean isTopLayer = HPSTrackerBuilder.getHalfFromName(name).equals("top") ? true : false;
+ if(isAxial) {
+ if(isTopLayer) {
+ l = 2*ll-1;
+ }
+ else {
+ l = 2*ll;
+ }
+ } else {
+ if(isTopLayer) {
+ l = 2*ll;
+ } else {
+ l = 2*ll-1;
+ }
+ }
+ } else {
+ throw new RuntimeException("Cannot get layer from name " + name);
+ }
+
+ if (ddd) {
+ System.out.println("sensor name " + name + " --> layer " + l);
+ }
+
+ if (l < 1 || l > 12) {
+ System.out.println("This layer doesn't exist?");
+ throw new RuntimeException("SiSensor name " + name + " is invalid?");
+ }
+
+ n[l - 1] = n[l - 1] + 1;
+
+ }
+
+ return n;
+ }
+
+ @Override
+ public void endOfData() {
+ if (outputPlots != null) {
+ try {
+ aida.saveAs(outputPlots);
+ } catch (IOException ex) {
+ Logger.getLogger(TrackingReconstructionPlots.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+ //plotterFrame.dispose();
+ //topFrame.dispose();
+ //bottomFrame.dispose();
+ }
+
+
+
+
+private void setupPlots() {
+
+
+ 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 Momentum (Px)", 25, -0.25, 0.25);
+ IHistogram1D trkPy = aida.histogram1D("Track Momentum (Py)", 25, -0.5, 0.5);
+ IHistogram1D trkPz = aida.histogram1D("Track Momentum (Pz)", 25, 0, 1.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);
+
+ if(showPlots) plotter.show();
+
+// ******************************************************************
+ top1 = fac.createPlotterFactory().create("Top Tracking Plots");
+ top1.setTitle("Top Momentum");
+ IPlotterStyle stop1 = top1.style();
+ stop1.dataStyle().fillStyle().setColor("green");
+ stop1.dataStyle().errorBarStyle().setVisible(false);
+ top1.createRegions(2, 2);
+ //topFrame.addPlotter(top1);
+
+ IHistogram1D toptrkPx = aida.histogram1D("Top Track Momentum (Px)", 25, -0.25, 0.25);
+ IHistogram1D toptrkPy = aida.histogram1D("Top Track Momentum (Py)", 25, -0.5, 0.5);
+ IHistogram1D toptrkPz = aida.histogram1D("Top Track Momentum (Pz)", 25, 0, 1.5);
+ IHistogram1D toptrkChi2 = aida.histogram1D("Top Track Chi2", 25, 0, 25.0);
+
+ top1.region(0).plot(toptrkPx);
+ top1.region(1).plot(toptrkPy);
+ top1.region(2).plot(toptrkPz);
+ top1.region(3).plot(toptrkChi2);
+
+ if(showPlots) top1.show();
+
+ bot1 = fac.createPlotterFactory().create("Bottom Tracking Plots");
+ bot1.setTitle("Bottom Momentum");
+ IPlotterStyle sbot1 = bot1.style();
+ sbot1.dataStyle().fillStyle().setColor("blue");
+ sbot1.dataStyle().errorBarStyle().setVisible(false);
+ bot1.createRegions(2, 2);
+ //bottomFrame.addPlotter(bot1);
+
+ IHistogram1D bottrkPx = aida.histogram1D("Bottom Track Momentum (Px)", 25, -0.25, 0.25);
+ IHistogram1D bottrkPy = aida.histogram1D("Bottom Track Momentum (Py)", 25, -0.5, 0.5);
+ IHistogram1D bottrkPz = aida.histogram1D("Bottom Track Momentum (Pz)", 25, 0, 1.5);
+ IHistogram1D bottrkChi2 = aida.histogram1D("Bottom Track Chi2", 25, 0, 25.0);
+
+ bot1.region(0).plot(bottrkPx);
+ bot1.region(1).plot(bottrkPy);
+ bot1.region(2).plot(bottrkPz);
+ bot1.region(3).plot(bottrkChi2);
+
+ if(showPlots) bot1.show();
+
+// ******************************************************************
+ IHistogram1D trkd0 = aida.histogram1D("d0 ", 25, -10.0, 10.0);
+ IHistogram1D trkphi = aida.histogram1D("sinphi ", 25, -0.2, 0.2);
+ IHistogram1D trkomega = aida.histogram1D("omega ", 25, -0.0025, 0.0025);
+ IHistogram1D trklam = aida.histogram1D("tan(lambda) ", 25, -0.1, 0.1);
+ IHistogram1D trkz0 = aida.histogram1D("z0 ", 25, -6.0, 6.0);
+
+ plotter22 = fac.createPlotterFactory().create("HPS Track Params");
+ plotter22.setTitle("Track parameters");
+ //plotterFrame.addPlotter(plotter22);
+ IPlotterStyle style22 = plotter22.style();
+ style22.dataStyle().fillStyle().setColor("yellow");
+ style22.dataStyle().errorBarStyle().setVisible(false);
+ plotter22.createRegions(2, 3);
+ plotter22.region(0).plot(trkd0);
+ plotter22.region(1).plot(trkphi);
+ plotter22.region(2).plot(trkomega);
+ plotter22.region(3).plot(trklam);
+ plotter22.region(4).plot(trkz0);
+
+ if(showPlots) plotter22.show();
+
+ // ******************************************************************
+
+
+ trkd0 = aida.histogram1D("d0 Top", 25, -10.0, 10.0);
+ trkphi = aida.histogram1D("sinphi Top", 25, -0.2, 0.2);
+ trkomega = aida.histogram1D("omega Top", 25, -0.0025, 0.0025);
+ trklam = aida.histogram1D("tan(lambda) Top", 25, -0.1, 0.1);
+ trkz0 = aida.histogram1D("z0 Top", 25, -6.0, 6.0);
+
+ plotter2221 = fac.createPlotterFactory().create("HPS Track Params");
+ plotter2221.setTitle("Track parameters");
+ //plotterFrame.addPlotter(plotter22);
+ IPlotterStyle style2221 = plotter2221.style();
+ style2221.dataStyle().fillStyle().setColor("yellow");
+ style2221.dataStyle().errorBarStyle().setVisible(false);
+ plotter2221.createRegions(2, 3);
+ plotter2221.region(0).plot(trkd0);
+ plotter2221.region(1).plot(trkphi);
+ plotter2221.region(2).plot(trkomega);
+ plotter2221.region(3).plot(trklam);
+ plotter2221.region(4).plot(trkz0);
+
+ if(showPlots) plotter2221.show();
+
+
+ // ******************************************************************
+
+
+ trkd0 = aida.histogram1D("d0 Bottom", 25, -10.0, 10.0);
+ trkphi = aida.histogram1D("sinphi Bottom", 25, -0.2, 0.2);
+ trkomega = aida.histogram1D("omega Bottom", 25, -0.0025, 0.0025);
+ trklam = aida.histogram1D("tan(lambda) Bottom", 25, -0.1, 0.1);
+ trkz0 = aida.histogram1D("z0 Bottom", 25, -6.0, 6.0);
+
+ plotter2222 = fac.createPlotterFactory().create("HPS Track Params");
+ plotter2222.setTitle("Track parameters");
+ //plotterFrame.addPlotter(plotter22);
+ IPlotterStyle style2222 = plotter2222.style();
+ style2222.dataStyle().fillStyle().setColor("yellow");
+ style2222.dataStyle().errorBarStyle().setVisible(false);
+ plotter2222.createRegions(2, 3);
+ plotter2222.region(0).plot(trkd0);
+ plotter2222.region(1).plot(trkphi);
+ plotter2222.region(2).plot(trkomega);
+ plotter2222.region(3).plot(trklam);
+ plotter2222.region(4).plot(trkz0);
+
+ if(showPlots) plotter2222.show();
+
+
+
+ // ******************************************************************
+
+
+ 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 xAtEcal2 = aida.histogram1D("X (mm) @ ECAL (Pz>1)", 50, -500, 500);
+ IHistogram1D yAtEcal2 = aida.histogram1D("Y (mm) @ ECAL (Pz>1)", 50, -100, 100);
+
+ 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(xAtEcal2);
+ plotter2.region(7).plot(yAtEcal2);
+
+ if(showPlots) plotter2.show();
+
+ // ******************************************************************
+
+ plotter222 = fac.createPlotterFactory().create("HPS Tracking Plots");
+ plotter222.setTitle("HPS Tracking Plots");
+ //plotterFrame.addPlotter(plotter222);
+ IPlotterStyle style222 = plotter222.style();
+ style222.dataStyle().fillStyle().setColor("yellow");
+ style222.dataStyle().errorBarStyle().setVisible(false);
+ plotter222.createRegions(2, 2);
+
+ IHistogram1D nHits = aida.histogram1D("Hits per Track", 5, 3, 8);
+ nTracks = aida.histogram1D("Tracks per Event", 3, 0, 3);
+ IHistogram1D nHitsCluster = aida.histogram1D("Hits in Cluster (HitOnTrack)", 4, 0, 4);
+
+
+ plotter222.region(0).plot(nHits);
+ plotter222.region(1).plot(nTracks);
+ plotter222.region(2).plot(nHitsCluster);
+
+ if(showPlots) plotter222.show();
+
+
+ // ******************************************************************
+
+ plotter22299 = fac.createPlotterFactory().create("HPS Tracking Plots Top");
+ plotter22299.setTitle("HPS Tracking Plots Top");
+ //plotterFrame.addPlotter(plotter22299);
+ IPlotterStyle style22299 = plotter22299.style();
+ style22299.dataStyle().fillStyle().setColor("yellow");
+ style22299.dataStyle().errorBarStyle().setVisible(false);
+ plotter22299.createRegions(2, 2);
+
+ IHistogram1D nHitsTop = aida.histogram1D("Hits per Track Top", 5, 3, 8);
+ nTracksTop = aida.histogram1D("Tracks per Event Top", 3, 0, 3);
+ IHistogram1D nHitsClusterTop = aida.histogram1D("Hits in Cluster (HitOnTrack) Top", 4, 0, 4);
+
+
+ plotter22299.region(0).plot(nHitsTop);
+ plotter22299.region(1).plot(nTracksTop);
+ plotter22299.region(2).plot(nHitsClusterTop);
+
+ if(showPlots) plotter22299.show();
+
+// ******************************************************************
+
+ plotter22298 = fac.createPlotterFactory().create("HPS Tracking Plots Bottom");
+ plotter22298.setTitle("HPS Tracking Plots Bottom");
+ //plotterFrame.addPlotter(plotter22298);
+ IPlotterStyle style22298 = plotter22298.style();
+ style22298.dataStyle().fillStyle().setColor("yellow");
+ style22298.dataStyle().errorBarStyle().setVisible(false);
+ plotter22298.createRegions(2, 2);
+
+ IHistogram1D nHitsBot = aida.histogram1D("Hits per Track Bot", 5, 3, 8);
+ nTracksBot = aida.histogram1D("Tracks per Event Bot", 3, 0, 3);
+ IHistogram1D nHitsClusterBot = aida.histogram1D("Hits in Cluster (HitOnTrack) Bot", 4, 0, 4);
+
+
+ plotter22298.region(0).plot(nHitsBot);
+ plotter22298.region(1).plot(nTracksBot);
+ plotter22298.region(2).plot(nHitsClusterBot);
+
+ if(showPlots) plotter22298.show();
+
+
+ // ******************************************************************
+
+
+ plotter2223 = fac.createPlotterFactory().create("Cluster Amp Plots");
+ plotter2223.setTitle("Other");
+ //plotterFrame.addPlotter(plotter222);
+ IPlotterStyle style2223 = plotter2223.style();
+ style2223.dataStyle().fillStyle().setColor("yellow");
+ style2223.dataStyle().errorBarStyle().setVisible(false);
+ plotter2223.createRegions(2, 2);
+
+
+
+ IHistogram1D amp = aida.histogram1D("Amp (HitOnTrack)", 50, 0, 5000);
+ IHistogram1D ampcl = aida.histogram1D("Cluster Amp (HitOnTrack)", 50, 0, 5000);
+ IHistogram1D amp2 = aida.histogram1D("Amp Pz>0.8 (HitOnTrack)", 50, 0, 5000);
+ IHistogram1D ampcl2 = aida.histogram1D("Cluster Amp Pz>0.8 (HitOnTrack)", 50, 0, 5000);
+
+
+ plotter2223.region(0).plot(amp);
+ plotter2223.region(1).plot(amp2);
+ plotter2223.region(2).plot(ampcl);
+ plotter2223.region(3).plot(ampcl2);
+
+ if(showPlots) plotter2223.show();
+
+// ******************************************************************
+
+
+ plotter2224 = fac.createPlotterFactory().create("t0 Plots");
+ plotter2224.setTitle("Other");
+ IPlotterStyle style2224 = plotter2224.style();
+ style2224.dataStyle().fillStyle().setColor("yellow");
+ style2224.dataStyle().errorBarStyle().setVisible(false);
+ plotter2224.createRegions(2, 2);
+
+ IHistogram1D t0 = aida.histogram1D("t0 (HitOnTrack)", 50, -100, 100);
+ IHistogram1D t0cl = aida.histogram1D("Cluster t0 (HitOnTrack)", 50, -100, 100);
+ IHistogram1D t02 = aida.histogram1D("t0 Pz>0.8 (HitOnTrack)", 50, -100, 100);
+ IHistogram1D t0cl2 = aida.histogram1D("Cluster t0 Pz>0.8 (HitOnTrack)", 50, -100, 100);
+
+ plotter2224.region(0).plot(t0);
+ plotter2224.region(1).plot(t0cl);
+ plotter2224.region(2).plot(t02);
+ plotter2224.region(3).plot(t0cl2);
+
+ if(showPlots) plotter2224.show();
+
+
+ // ******************************************************************
+
+ plotter3 = fac.createPlotterFactory().create("HPS Layer Residual Plots");
+ plotter3.setTitle("Layer Residuals");
+ //plotterFrame.addPlotter(plotter3);
+ IPlotterStyle style3 = plotter3.style();
+ style3.dataStyle().fillStyle().setColor("yellow");
+ style3.dataStyle().errorBarStyle().setVisible(false);
+ plotter3.createRegions(6, 2);
+
+ plotter3_1 = fac.createPlotterFactory().create("HPS Residual Plots (Single hit per layer)");
+ plotter3_1.setTitle("Residuals (Top)");
+ //plotterFrame.addPlotter(plotter3_1);
+ IPlotterStyle style3_1 = plotter3_1.style();
+ style3_1.dataStyle().fillStyle().setColor("yellow");
+ style3_1.dataStyle().errorBarStyle().setVisible(false);
+ plotter3_1.createRegions(6, 2);
+
+ plotter3_2 = fac.createPlotterFactory().create("HPS Residual Plots (Single strip cluster per layer)");
+ plotter3_2.setTitle("Residuals (Bottom)");
+ //plotterFrame.addPlotter(plotter3_2);
+ IPlotterStyle style3_2 = plotter3_2.style();
+ style3_2.dataStyle().fillStyle().setColor("yellow");
+ style3_2.dataStyle().errorBarStyle().setVisible(false);
+ plotter3_2.createRegions(6, 2);
+
+
+
+ for(int l=0; l < nLayers; ++l) {
+ for(int h=0; h <3;++h) {
+ String half;
+ if (h==0) half = "";
+ else if (h==1) half = " Top";
+ else half = " Bottom";
+ String name = "Layer " + String.valueOf(l+1) + " Residual X(mm)" + half;
+ IHistogram1D mod1ResX = aida.histogram1D(name, 25, -1, 1);
+ name = "Layer " + String.valueOf(l+1) + " Residual Y(mm)" + half;
+ IHistogram1D mod1ResY = aida.histogram1D(name, 25, -1, 1);
+
+ if(l>5) {
+ LOGGER.warning("cannot make plots for layer " + String.valueOf(l + 1));
+ continue;
+ }
+ if (h==0) {
+ plotter3.region(l*2).plot(mod1ResX);
+ plotter3.region(l*2 + 1).plot(mod1ResY);
+ }
+ else if (h==0) {
+ plotter3_1.region(l*2).plot(mod1ResX);
+ plotter3_1.region(l*2 + 1).plot(mod1ResY);
+ }
+ else {
+ plotter3_2.region(l*2).plot(mod1ResX);
+ plotter3_2.region(l*2 + 1).plot(mod1ResY);
+ }
+ }
+ }
+ if(showPlots) {
+ plotter3.show();
+ plotter3_1.show();
+ plotter3_2.show();
+ }
+
+
+
+ plotter3_11 = fac.createPlotterFactory().create("HPS Strip Residual Plots");
+ plotter3_11.setTitle("Strip Residuals (Top)");
+ //plotterFrame.addPlotter(plotter3_11);
+ IPlotterStyle style3_11 = plotter3_11.style();
+ style3_11.dataStyle().fillStyle().setColor("yellow");
+ style3_11.dataStyle().errorBarStyle().setVisible(false);
+ plotter3_11.createRegions(6, 6);
+ int i=0;
+ double[] limits = {1.,1.5,3.,4.,5.,5.,5.0};
+ for(HpsSiSensor sensor : sensors) {
+ double min = 0.0;
+ double max = 0.0;
+ int l = sensor.getLayerNumberFromDetectorElement();
+ max = limits[l-1];
+ min = -1.0*limits[l-1];
+ IHistogram1D resX = aida.histogram1D(sensor.getName() + " strip residual (mm)", 50, min, max);
+ if(l>6) {
+ LOGGER.warning("cannot make plots for this sensor " + sensor.getName());
+ continue;
+ }
+ plotter3_11.region(i).plot(resX);
+ i++;
+ }
+
+ if(showPlots) plotter3_11.show();
+
+
+
+
+
+
+ plotter4 = fac.createPlotterFactory().create("HPS Track and ECal Plots");
+ plotter4.setTitle("Track and ECal Correlations");
+ //plotterFrame.addPlotter(plotter4);
+ IPlotterStyle style4 = plotter4.style();
+ style4.setParameter("hist2DStyle", "colorMap");
+ style4.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ style4.dataStyle().fillStyle().setColor("yellow");
+ style4.dataStyle().errorBarStyle().setVisible(false);
+ plotter4.createRegions(2, 3);
+
+ IHistogram2D eVsP = aida.histogram2D("Energy Vs Momentum", 50, 0, 0.50, 50, 0, 1.5);
+ IHistogram1D eOverP = aida.histogram1D("Energy Over Momentum", 50, 0, 2);
+
+ IHistogram1D distX = aida.histogram1D("deltaX", 50, -100, 100);
+ IHistogram1D distY = aida.histogram1D("deltaY", 50, -40, 40);
+
+ IHistogram2D xEcalVsTrk = aida.histogram2D("X ECal Vs Track", 100, -400, 400, 100, -400, 400);
+ IHistogram2D yEcalVsTrk = aida.histogram2D("Y ECal Vs Track", 100, -100, 100, 100, -100, 100);
+
+ plotter4.region(0).plot(eVsP);
+ plotter4.region(3).plot(eOverP);
+ plotter4.region(1).plot(distX);
+ plotter4.region(4).plot(distY);
+ plotter4.region(2).plot(xEcalVsTrk);
+ plotter4.region(5).plot(yEcalVsTrk);
+
+ if(showPlots) plotter4.show();
+
+ // ******************************************************************
+ top2 = fac.createPlotterFactory().create("Top ECal Plots");
+ top2.setTitle("Top ECal Correlations");
+ IPlotterStyle stop2 = top2.style();
+ stop2.dataStyle().fillStyle().setColor("green");
+ stop2.dataStyle().errorBarStyle().setVisible(false);
+ stop2.setParameter("hist2DStyle", "colorMap");
+ stop2.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ top2.createRegions(2, 3);
+ //topFrame.addPlotter(top2);
+
+ IHistogram2D topeVsP = aida.histogram2D("Top Energy Vs Momentum", 50, 0, 0.500, 50, 0, 1.5);
+ IHistogram1D topeOverP = aida.histogram1D("Top Energy Over Momentum", 50, 0, 2);
+
+ IHistogram1D topdistX = aida.histogram1D("Top deltaX", 50, -100, 100);
+ IHistogram1D topdistY = aida.histogram1D("Top deltaY", 50, -40, 40);
+
+ IHistogram2D topxEcalVsTrk = aida.histogram2D("Top X ECal Vs Track", 100, -400, 400, 100, -100, 100);
+ IHistogram2D topyEcalVsTrk = aida.histogram2D("Top Y ECal Vs Track", 100, 0, 100, 100, 0, 100);
+
+ top2.region(0).plot(topeVsP);
+ top2.region(3).plot(topeOverP);
+ top2.region(1).plot(topdistX);
+ top2.region(4).plot(topdistY);
+ top2.region(2).plot(topxEcalVsTrk);
+ top2.region(5).plot(topyEcalVsTrk);
+
+ if(showPlots) top2.show();
+
+ bot2 = fac.createPlotterFactory().create("Bottom ECal Plots");
+ bot2.setTitle("Bottom ECal Correlations");
+ IPlotterStyle sbot2 = bot2.style();
+ sbot2.dataStyle().fillStyle().setColor("green");
+ sbot2.dataStyle().errorBarStyle().setVisible(false);
+ sbot2.setParameter("hist2DStyle", "colorMap");
+ sbot2.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ bot2.createRegions(2, 3);
+ //bottomFrame.addPlotter(bot2);
+
+ IHistogram2D BottomeVsP = aida.histogram2D("Bottom Energy Vs Momentum", 50, 0, 0.500, 50, 0, 1.5);
+ IHistogram1D BottomeOverP = aida.histogram1D("Bottom Energy Over Momentum", 50, 0, 2);
+
+ IHistogram1D BottomdistX = aida.histogram1D("Bottom deltaX", 50, -100, 100);
+ IHistogram1D BottomdistY = aida.histogram1D("Bottom deltaY", 50, -40, 40);
+
+ IHistogram2D BottomxEcalVsTrk = aida.histogram2D("Bottom X ECal Vs Track", 100, -400, 400, 100, -400, 400);
+ IHistogram2D BottomyEcalVsTrk = aida.histogram2D("Bottom Y ECal Vs Track", 100, -100, 0, 100, -100, 0);
+
+ bot2.region(0).plot(BottomeVsP);
+ bot2.region(3).plot(BottomeOverP);
+ bot2.region(1).plot(BottomdistX);
+ bot2.region(4).plot(BottomdistY);
+ bot2.region(2).plot(BottomxEcalVsTrk);
+ bot2.region(5).plot(BottomyEcalVsTrk);
+
+ if(showPlots) bot2.show();
+
+
+ // ******************************************************************
+ top3 = fac.createPlotterFactory().create("Top ECal Plots");
+ top3.setTitle("Top ECal More Correlations");
+ IPlotterStyle stop3 = top3.style();
+ stop3.dataStyle().fillStyle().setColor("green");
+ stop3.dataStyle().errorBarStyle().setVisible(false);
+ stop3.setParameter("hist2DStyle", "colorMap");
+ stop3.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ top3.createRegions(1, 2);
+ //topFrame.addPlotter(top3);
+
+ IHistogram2D topdistXvsX = aida.histogram2D("Top deltaX vs X", 51, -400, 400, 25, -100, 100);
+ IHistogram2D topdistYvsY = aida.histogram2D("Top deltaY vs Y", 51, 0, 100, 25, -40, 40);
+
+ top3.region(0).plot(topdistXvsX);
+ top3.region(1).plot(topdistYvsY);
+
+ if(showPlots) top3.show();
+
+ bot3 = fac.createPlotterFactory().create("Bottom ECal Plots");
+ bot3.setTitle("Bottom ECal More Correlations");
+ IPlotterStyle sbot3 = bot3.style();
+ sbot3.dataStyle().fillStyle().setColor("green");
+ sbot3.dataStyle().errorBarStyle().setVisible(false);
+ sbot3.setParameter("hist2DStyle", "colorMap");
+ sbot3.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ bot3.createRegions(1, 2);
+ //bottomFrame.addPlotter(bot3);
+
+ IHistogram2D botdistXvsX = aida.histogram2D("Bottom deltaX vs X", 51, -400, 400, 25, -100, 100);
+ IHistogram2D botdistYvsY = aida.histogram2D("Bottom deltaY vs Y", 51, -100, 0, 25, -40, 40);
+
+ bot3.region(0).plot(botdistXvsX);
+ bot3.region(1).plot(botdistYvsY);
+
+ if(showPlots) bot3.show();
+
+ // ******************************************************************
+ top4 = fac.createPlotterFactory().create("Track Matching Plots");
+ top4.setTitle("Track Matching Plots");
+ IPlotterStyle stop4 = top4.style();
+ stop4.dataStyle().fillStyle().setColor("green");
+ stop4.dataStyle().errorBarStyle().setVisible(false);
+ stop4.setParameter("hist2DStyle", "colorMap");
+ stop4.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ top4.createRegions(2, 3);
+ //topFrame.addPlotter(top4);
+
+ IHistogram1D trackmatchN = aida.histogram1D("Tracks matched", 3, 0, 3);
+ IHistogram1D toptrackmatchN = aida.histogram1D("Tracks matched Top", 3, 0, 3);
+ IHistogram1D bottrackmatchN = aida.histogram1D("Tracks matched Bottom", 3, 0, 3);
+ IHistogram1D trackmatchN2 = aida.histogram1D("Tracks matched (Pz>0.8)", 3, 0, 3);
+ IHistogram1D toptrackmatchN2 = aida.histogram1D("Tracks matched Top (Pz>0.8)", 3, 0, 3);
+ IHistogram1D bottrackmatchN2 = aida.histogram1D("Tracks matched Bottom (Pz>0.8)", 3, 0, 3);
+
+ top4.region(0).plot(trackmatchN);
+ top4.region(1).plot(toptrackmatchN);
+ top4.region(2).plot(bottrackmatchN);
+ top4.region(3).plot(trackmatchN2);
+ top4.region(4).plot(toptrackmatchN2);
+ top4.region(5).plot(bottrackmatchN2);
+
+ if(showPlots) top4.show();
+
+ // ******************************************************************
+ top44 = fac.createPlotterFactory().create("e+e- Plots");
+ top44.setTitle("e+e- Plots");
+ IPlotterStyle stop44 = top44.style();
+ stop44.dataStyle().fillStyle().setColor("green");
+ stop44.dataStyle().errorBarStyle().setVisible(false);
+ stop44.setParameter("hist2DStyle", "colorMap");
+ stop44.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ top44.createRegions(2,4);
+ //topFrame.addPlotter(top44);
+
+ IHistogram2D trackPCorr = aida.histogram2D("p(e-) vs p(e+) max", 25, 0, 1.2, 25, 0, 1.2);
+ IHistogram1D ne = aida.histogram1D("n(e-)", 3, 0, 3);
+ IHistogram1D np = aida.histogram1D("n(e+)", 3, 0, 3);
+ IHistogram1D pem = aida.histogram1D("p(e-) max", 25, 0, 1.5);
+ IHistogram1D pe = aida.histogram1D("p(e-)", 25, 0, 1.5);
+ IHistogram1D ppm = aida.histogram1D("p(e+) max", 25, 0, 1.5);
+ IHistogram1D pp = aida.histogram1D("p(e+)", 25, 0, 1.5);
+
+ top44.region(0).plot(trackPCorr);
+ top44.region(1).plot(ne);
+ top44.region(2).plot(np);
+ top44.region(3).plot(pe);
+ top44.region(4).plot(pp);
+ top44.region(5).plot(pem);
+ top44.region(6).plot(ppm);
+
+ if(showPlots) top44.show();
+
+
+
+// ******************************************************************
+ plotter5 = fac.createPlotterFactory().create("HPS Hit Positions");
+ plotter5.setTitle("Hit Positions: Top");
+ //plotterFrame.addPlotter(plotter5);
+ IPlotterStyle style5 = plotter5.style();
+ style5.setParameter("hist2DStyle", "colorMap");
+ style5.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ style5.dataStyle().fillStyle().setColor("yellow");
+ style5.dataStyle().errorBarStyle().setVisible(false);
+ plotter5.createRegions(1, 2);
+
+ IHistogram2D l1Pos = aida.histogram2D("Layer 1 HTH Position: Top", 50, -55, 55, 55, -25, 25);
+ IHistogram2D l7Pos = aida.histogram2D("Layer 7 HTH Position: Top", 50, -55, 55, 55, -25, 25);
+
+ plotter5.region(0).plot(l1Pos);
+ plotter5.region(1).plot(l7Pos);
+
+ if(showPlots) plotter5.show();
+
+ plotter5_1 = fac.createPlotterFactory().create("HPS Hit Positions");
+ plotter5_1.setTitle("Hit Positions: Bottom");
+ //plotterFrame.addPlotter(plotter5_1);
+ IPlotterStyle style5_1 = plotter5_1.style();
+ style5_1.setParameter("hist2DStyle", "colorMap");
+ style5_1.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ style5_1.dataStyle().fillStyle().setColor("yellow");
+ style5_1.dataStyle().errorBarStyle().setVisible(false);
+ plotter5_1.createRegions(1, 2);
+
+
+ IHistogram2D l1PosBot = aida.histogram2D("Layer 1 HTH Position: Bottom", 50, -55, 55, 55, -25, 25);
+ IHistogram2D l7PosBot = aida.histogram2D("Layer 7 HTH Position: Bottom", 50, -55, 55, 55, -25, 25);
+ plotter5_1.region(0).plot(l1PosBot);
+ plotter5_1.region(1).plot(l7PosBot);
+
+ if(showPlots) plotter5_1.show();
+
+ plotter55 = fac.createPlotterFactory().create("HPS Hit Positions");
+ plotter55.setTitle("Helical Track Hits");
+ //plotterFrame.addPlotter(plotter55);
+ IPlotterStyle style55 = plotter55.style();
+ style55.dataStyle().fillStyle().setColor("Green");
+ style55.dataStyle().errorBarStyle().setVisible(false);
+ style55.dataStyle().markerStyle().setSize(20);
+ plotter55.createRegions(1, 2);
+
+ IProfile avgLayersTopPlot = aida.profile1D("Number of Stereo Hits per layer in Top Half", 10, 0, 10);
+ IProfile avgLayersBottomPlot = aida.profile1D("Number of Stereo Hits per layer in Bottom Half", 10, 0, 10);
+
+ plotter55.region(0).plot(avgLayersTopPlot);
+ plotter55.region(1).plot(avgLayersBottomPlot);
+
+ if(showPlots) plotter55.show();
+
+ plotter6 = fac.createPlotterFactory().create("HPS ECAL Hit Positions");
+ plotter6.setTitle("ECAL Positions");
+ //plotterFrame.addPlotter(plotter6);
+ IPlotterStyle style6 = plotter6.style();
+ style6.setParameter("hist2DStyle", "colorMap");
+ style6.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ style6.dataStyle().fillStyle().setColor("yellow");
+ style6.dataStyle().errorBarStyle().setVisible(false);
+ plotter6.createRegions(4, 2);
+
+ IHistogram2D topECal = aida.histogram2D("Top ECal Cluster Position", 50, -400, 400, 10, 0, 100);
+ IHistogram2D botECal = aida.histogram2D("Bottom ECal Cluster Position", 50, -400, 400, 10, -100, 0);
+ IHistogram2D topECal1 = aida.histogram2D("Top ECal Cluster Position (>0 tracks)", 50, -400, 400, 10, 0, 100);
+ IHistogram2D botECal1 = aida.histogram2D("Bottom ECal Cluster Position (>0 tracks)", 50, -400, 400, 10, -100, 0);
+ IHistogram2D topECal2 = aida.histogram2D("Top ECal Cluster Position (E>0.1,>0 tracks)", 50, -400, 400, 10, 0, 100);
+ IHistogram2D botECal2 = aida.histogram2D("Bottom ECal Cluster Position (E>0.1,>0 tracks)", 50, -400, 400, 10, -100, 0);
+ IHistogram2D topECal3 = aida.histogram2D("Top ECal Cluster Position w_E (E>0.1,>0 tracks)", 50, -400, 400, 10, 0, 100);
+ IHistogram2D botECal3 = aida.histogram2D("Bottom ECal Cluster Position w_E (E>0.1,>0 tracks)", 50, -400, 400, 10, -100, 0);
+
+ plotter6.region(0).plot(topECal);
+ plotter6.region(1).plot(botECal);
+ plotter6.region(2).plot(topECal1);
+ plotter6.region(3).plot(botECal1);
+ plotter6.region(4).plot(topECal2);
+ plotter6.region(5).plot(botECal2);
+ plotter6.region(6).plot(topECal3);
+ plotter6.region(7).plot(botECal3);
+
+ if(showPlots) plotter6.show();
+
+
+ plotter66 = fac.createPlotterFactory().create("HPS ECAL Basic Plots");
+ plotter66.setTitle("ECAL Basic Plots");
+ //plotterFrame.addPlotter(plotter6);
+ IPlotterStyle style66 = plotter66.style();
+ style66.dataStyle().fillStyle().setColor("yellow");
+ style66.dataStyle().errorBarStyle().setVisible(false);
+ plotter66.createRegions(2, 2);
+
+ IHistogram1D topECalE = aida.histogram1D("Top ECal Cluster Energy", 50, 0, 2);
+ IHistogram1D botECalE = aida.histogram1D("Bottom ECal Cluster Energy", 50, 0, 2);
+ IHistogram1D topECalN = aida.histogram1D("Number of Clusters Top", 6, 0, 6);
+ IHistogram1D botECalN = aida.histogram1D("Number of Clusters Bot", 6, 0, 6);
+
+ plotter66.region(0).plot(topECalE);
+ plotter66.region(1).plot(botECalE);
+ plotter66.region(2).plot(botECalN);
+ plotter66.region(3).plot(topECalN);
+
+ if(showPlots) plotter66.show();
+
+
+
+
+ plotter8 = fac.createPlotterFactory().create("HPS Strip Hit From Stereo Multiplicity");
+ plotter8.setTitle("Strip Hit Multiplicity");
+ //plotterFrame.addPlotter(plotter8);
+ IPlotterStyle style8 = plotter8.style();
+ style8.dataStyle().fillStyle().setColor("yellow");
+ style8.dataStyle().errorBarStyle().setVisible(false);
+ plotter8.createRegions(6, 6);
+ i=0;
+ for(SiSensor sensor : sensors) {
+
+ IHistogram1D resX = aida.histogram1D(sensor.getName() + " strip hits from stereo", 10, 0, 10);
+ if(sensor.getName().contains("L7")) {
+ LOGGER.warning("cannot setup this plot. Fix.");
+ continue;
+ }
+ plotter8.region(i).plot(resX);
+ i++;
+ }
+
+ if(showPlots) plotter8.show();
+
+ plotter88 = fac.createPlotterFactory().create("HPS Strip Hit Multiplicity");
+ plotter88.setTitle("Strip Hit Multiplicity");
+ //plotterFrame.addPlotter(plotter88);
+ plotter88.setStyle(style8);
+ plotter88.createRegions(6, 6);
+ i=0;
+ for(SiSensor sensor : sensors) {
+
+ IHistogram1D resX = aida.histogram1D(sensor.getName() + " strip hits", 10, 0, 10);
+ if(sensor.getName().contains("L7")) {
+ LOGGER.warning("cannot setup this plot. Fix.");
+ continue;
+ }
+ plotter88.region(i).plot(resX);
+ i++;
+ }
+
+ if(showPlots) plotter88.show();
+
+
+
+
+
+ plotter888 = fac.createPlotterFactory().create("HPS Strip Hit Isolation");
+ plotter888.setTitle("Strip Hit Isolation");
+ //plotterFrame.addPlotter(plotter88);
+ plotter888.setStyle(style8);
+ plotter888.createRegions(6, 6);
+ i=0;
+ for(SiSensor sensor : sensors) {
+
+ IHistogram1D resX = aida.histogram1D(sensor.getName() + " strip hits iso", 50, 0, 5);
+ if(sensor.getName().contains("L7")) {
+ LOGGER.warning("cannot setup this plot. Fix.");
+ continue;
+ }
+ plotter888.region(i).plot(resX);
+ i++;
+ }
+
+ if(showPlots) plotter888.show();
+
+ plotter8888 = fac.createPlotterFactory().create("HPS Strip Hit On Track Isolation");
+ plotter8888.setTitle("Strip Hit On Track Isolation");
+ //plotterFrame.addPlotter(plotter88);
+ plotter8888.setStyle(style8);
+ plotter8888.createRegions(6, 6);
+ i=0;
+ for(SiSensor sensor : sensors) {
+
+ IHistogram1D resX = aida.histogram1D(sensor.getName() + " strip hits iso on track", 50, 0, 5);
+ if(sensor.getName().contains("L7")) {
+ LOGGER.warning("cannot setup this plot. Fix.");
+ continue;
+ }
+ plotter8888.region(i).plot(resX);
+ i++;
+ }
+
+ if(showPlots) plotter8888.show();
+
+
+ }
+
+
+
+
+ private Cluster findClosestCluster(Hep3Vector posonhelix, List<Cluster> clusters) {
+ Cluster closest = null;
+ double minDist = 9999;
+ for (Cluster cluster : clusters) {
+ double[] clPos = cluster.getPosition();
+ double clEne = cluster.getEnergy();
+ double dist = Math.sqrt(Math.pow(clPos[0] - posonhelix.x(), 2) + Math.pow(clPos[1] - posonhelix.y(), 2)); //coordinates!!!
+ if (dist < minDist && clEne > 0.4) {
+ closest = cluster;
+ minDist = dist;
+ }
+ }
+ return closest;
+ }
+}
