java/trunk/analysis/src/main/java/org/hps/analysis/dataquality
--- java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/TrackMCEfficiency.java 2014-09-10 18:29:35 UTC (rev 983)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/TrackMCEfficiency.java 2014-09-10 20:31:54 UTC (rev 984)
@@ -8,6 +8,7 @@
import java.util.List;
import java.util.Map;
import java.util.Set;
+import org.hps.analysis.examples.LCIOTrackAnalysis;
import org.hps.analysis.examples.TrackAnalysis;
import org.hps.recon.tracking.FindableTrack;
import org.hps.recon.tracking.FittedRawTrackerHit;
@@ -19,6 +20,7 @@
import org.lcsim.event.RelationalTable;
import org.lcsim.event.SimTrackerHit;
import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
import org.lcsim.event.base.BaseRelationalTable;
import org.lcsim.fit.helicaltrack.HelicalTrackCross;
import org.lcsim.fit.helicaltrack.HelixParamCalculator;
@@ -28,9 +30,11 @@
import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
/**
- * DQM driver for the monte carlo track efficiency; makes a bunch of efficiency vs variable plots
- * for all tracks and just electrons from trident/A' event, as well as "findable" tracks
- * use the debugTrackEfficiency flag to print out info regarding individual failed events
+ * DQM driver for the monte carlo track efficiency; makes a bunch of efficiency
+ * vs variable plots for all tracks and just electrons from trident/A' event, as
+ * well as "findable" tracks use the debugTrackEfficiency flag to print out info
+ * regarding individual failed events
+ *
* @author mgraham on Mar 28, 2014
*/
// TODO: Add some quantities for DQM monitoring: e.g. <efficiency>, <eff>_findable
@@ -42,6 +46,8 @@
private String trackerHitCollectionName = "TrackerHits";
private String siClusterCollectionName = "StripClusterer_SiTrackerHitStrip1D";
private String rotatedMCRelationsCollectionName = "RotatedHelicalTrackMCRelations";
+ private final String helicalTrackHitRelationsCollectionName = "HelicalTrackHitRelations";
+ private final String rotatedHelicalTrackHitRelationsCollectionName = "RotatedHelicalTrackHitRelations";
private String trackCollectionName = "MatchedTracks";
private String trackerName = "Tracker";
private Detector detector = null;
@@ -62,7 +68,8 @@
private static final String nameStrip = "Tracker_TestRunModule_";
private List<SiSensor> sensors;
private boolean debugTrackEfficiency = false;
- private String plotDir = "TrackMCEfficiency/";
+ private String plotDir = "TrackMCEfficiency/";
+
public void setHelicalTrackHitCollectionName(String helicalTrackHitCollectionName) {
this.helicalTrackHitCollectionName = helicalTrackHitCollectionName;
}
@@ -82,14 +89,14 @@
aida.tree().cd("/");
IHistogramFactory hf = aida.histogramFactory();
- peffFindable = hf.createProfile1D(plotDir+"Findable Efficiency vs p", "", 20, 0., beamP);
- phieffFindable = hf.createProfile1D(plotDir+"Findable Efficiency vs phi", "", 25, -0.25, 0.25);
- ctheffFindable = hf.createProfile1D(plotDir+"Findable Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
+ peffFindable = hf.createProfile1D(plotDir + "Findable Efficiency vs p", "", 20, 0., beamP);
+ phieffFindable = hf.createProfile1D(plotDir + "Findable Efficiency vs phi", "", 25, -0.25, 0.25);
+ ctheffFindable = hf.createProfile1D(plotDir + "Findable Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
- peffElectrons = hf.createProfile1D(plotDir+"Electrons Efficiency vs p", "", 20, 0., beamP);
- phieffElectrons = hf.createProfile1D(plotDir+"Electrons Efficiency vs phi", "", 25, -0.25, 0.25);
- ctheffElectrons = hf.createProfile1D(plotDir+"Electrons Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
-
+ peffElectrons = hf.createProfile1D(plotDir + "Electrons Efficiency vs p", "", 20, 0., beamP);
+ phieffElectrons = hf.createProfile1D(plotDir + "Electrons Efficiency vs phi", "", 25, -0.25, 0.25);
+ ctheffElectrons = hf.createProfile1D(plotDir + "Electrons Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
+
}
@Override
@@ -98,19 +105,25 @@
aida.tree().cd("/");
//make sure the required collections exist
- if (!event.hasCollection(RawTrackerHit.class, rawTrackerHitCollectionName))
+ if (!event.hasCollection(RawTrackerHit.class, rawTrackerHitCollectionName)) {
return;
- if (!event.hasCollection(FittedRawTrackerHit.class, fittedTrackerHitCollectionName))
+ }
+ if (!event.hasCollection(LCRelation.class, fittedTrackerHitCollectionName)) {
return;
- if (!event.hasCollection(Track.class, trackCollectionName))
+ }
+ if (!event.hasCollection(Track.class, trackCollectionName)) {
return;
- if (!event.hasCollection(LCRelation.class, rotatedMCRelationsCollectionName))
+ }
+ if (!event.hasCollection(LCRelation.class, rotatedMCRelationsCollectionName)) {
return;
- if (!event.hasCollection(SiTrackerHitStrip1D.class, siClusterCollectionName))
+ }
+ if (!event.hasCollection(TrackerHit.class, siClusterCollectionName)) {
return;
+ }
- if (!event.hasCollection(SimTrackerHit.class, trackerHitCollectionName))
+ if (!event.hasCollection(SimTrackerHit.class, trackerHitCollectionName)) {
return;
+ }
//
//get the b-field
Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
@@ -120,17 +133,20 @@
RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
List<LCRelation> mcrelations = event.get(LCRelation.class, rotatedMCRelationsCollectionName);
for (LCRelation relation : mcrelations) {
- if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null) {
hittomc.add(relation.getFrom(), relation.getTo());
+ }
}
+
RelationalTable mcHittomcP = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
// Get the collections of SimTrackerHits
List<List<SimTrackerHit>> simcols = event.get(SimTrackerHit.class);
// Loop over the SimTrackerHits and fill in the relational table
for (List<SimTrackerHit> simlist : simcols) {
for (SimTrackerHit simhit : simlist) {
- if (simhit.getMCParticle() != null)
+ if (simhit.getMCParticle() != null) {
mcHittomcP.add(simhit, simhit.getMCParticle());
+ }
}
}
RelationalTable trktomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
@@ -138,47 +154,68 @@
if (event.hasCollection(LCRelation.class, "SVTTrueHitRelations")) {
List<LCRelation> trueHitRelations = event.get(LCRelation.class, "SVTTrueHitRelations");
for (LCRelation relation : trueHitRelations) {
- if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null) {
rawtomc.add(relation.getFrom(), relation.getTo());
+ }
}
}
// make relational table for strip clusters to mc particle
- List<SiTrackerHitStrip1D> siClusters = event.get(SiTrackerHitStrip1D.class, siClusterCollectionName);
+ List<TrackerHit> siClusters = event.get(TrackerHit.class, siClusterCollectionName);
RelationalTable clustertosimhit = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
- for (SiTrackerHit cluster : siClusters) {
+ for (TrackerHit cluster : siClusters) {
List<RawTrackerHit> rawHits = cluster.getRawHits();
for (RawTrackerHit rth : rawHits) {
Set<SimTrackerHit> simTrackerHits = rawtomc.allFrom(rth);
- if (simTrackerHits != null)
+ if (simTrackerHits != null) {
for (SimTrackerHit simhit : simTrackerHits) {
clustertosimhit.add(cluster, simhit);
}
+ }
}
}
//relational tables from mc particle to raw and fitted tracker hits
RelationalTable fittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
- List<FittedRawTrackerHit> fittedTrackerHits = event.get(FittedRawTrackerHit.class, fittedTrackerHitCollectionName);
- for (FittedRawTrackerHit hit : fittedTrackerHits) {
- RawTrackerHit rth = hit.getRawTrackerHit();
+ List<LCRelation> fittedTrackerHits = event.get(LCRelation.class, fittedTrackerHitCollectionName);
+ for (LCRelation hit : fittedTrackerHits) {
+ RawTrackerHit rth = FittedRawTrackerHit.getRawTrackerHit(hit);
Set<SimTrackerHit> simTrackerHits = rawtomc.allFrom(rth);
- if (simTrackerHits != null)
+ if (simTrackerHits != null) {
for (SimTrackerHit simhit : simTrackerHits) {
- if (simhit.getMCParticle() != null)
+ if (simhit.getMCParticle() != null) {
fittomc.add(hit, simhit.getMCParticle());
+ }
}
+ }
}
+ RelationalTable hittostrip = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> hitrelations = event.get(LCRelation.class, helicalTrackHitRelationsCollectionName);
+ for (LCRelation relation : hitrelations) {
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null) {
+ hittostrip.add(relation.getFrom(), relation.getTo());
+ }
+ }
+
+ RelationalTable hittorotated = new BaseRelationalTable(RelationalTable.Mode.ONE_TO_ONE, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> rotaterelations = event.get(LCRelation.class, rotatedHelicalTrackHitRelationsCollectionName);
+ for (LCRelation relation : rotaterelations) {
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null) {
+ hittorotated.add(relation.getFrom(), relation.getTo());
+ }
+ }
+
// Instantiate the class that determines if a track is "findable"
FindableTrack findable = new FindableTrack(event);
- List<Track> tracks = event.get(Track.class, trackCollectionName);
+ List<Track> tracks = event.get(Track.class, trackCollectionName);
for (Track trk : tracks) {
- TrackAnalysis tkanal = new TrackAnalysis(trk, hittomc);
+ TrackAnalysis tkanal = new TrackAnalysis(trk, hittomc, rawtomc, hittostrip, hittorotated);
tkanalMap.put(trk, tkanal);
MCParticle mcp = tkanal.getMCParticleNew();
- if (mcp != null)
- // Create a map between the tracks found and the assigned MC particle
+ if (mcp != null) // Create a map between the tracks found and the assigned MC particle
+ {
trktomc.add(trk, tkanal.getMCParticleNew());
+ }
}
// Now loop over all MC Particles
@@ -193,10 +230,10 @@
double pz = mcp.getPZ();
double pt = Math.sqrt(px * px + py * py);
double p = Math.sqrt(pt * pt + pz * pz);
- double cth = pz / p;
+ double cth = py / p;
double theta = 180. * Math.acos(cth) / Math.PI;
double eta = -Math.log(Math.tan(Math.atan2(pt, pz) / 2));
- double phi = Math.atan2(py, px);
+ double phi = Math.atan2(px, pz);
// Find the number of layers hit by this mc particle
// System.out.println("MC pt=" + pt);
int nhits = findable.LayersHit(mcp);
@@ -219,14 +256,17 @@
//it's the A'...let's see if we found both tracks.
List<MCParticle> daughters = mcp.getDaughters();
for (MCParticle d : daughters) {
- if (trktomc.allTo(d).isEmpty())
+ if (trktomc.allTo(d).isEmpty()) {
bothreco = false;
- if (!findable.InnerTrackerIsFindable(d, nlayers - 2))
+ }
+ if (!findable.InnerTrackerIsFindable(d, nlayers - 2)) {
bothfindable = false;
+ }
}
double vtxWgt = 0;
- if (bothreco)
+ if (bothreco) {
vtxWgt = 1.0;
+ }
// VxEff.fill(mcp.getOriginX(), vtxWgt);
// VyEff.fill(mcp.getOriginY(), vtxWgt);
// VzEff.fill(mcp.getOriginZ(), vtxWgt);
@@ -242,13 +282,13 @@
_nchMCP++;
findableTracks++;
double wgt = 0.;
- if (ntrk > 0)
+ if (ntrk > 0) {
wgt = 1.;
+ }
foundTracks += wgt;
peffFindable.fill(p, wgt);
phieffFindable.fill(phi, wgt);
ctheffFindable.fill(cth, wgt);
-
if (wgt == 0) {
Set<SimTrackerHit> mchitlist = mcHittomcP.allTo(mcp);
@@ -256,8 +296,9 @@
Set<FittedRawTrackerHit> fitlist = fittomc.allTo(mcp);
if (debugTrackEfficiency) {
System.out.println("TrackMCEfficiencyMonitoring:: Missed a findable track with MC p = " + p);
- if (!hasHTHInEachLayer(hitlist, fitlist))
+ if (!hasHTHInEachLayer(hitlist, fitlist)) {
System.out.println("This track failed becasue it's missing a helical track hit");
+ }
}
}
@@ -266,13 +307,13 @@
totelectrons++;
// findableelectrons++;
double wgt = 0.;
- if (ntrk > 0)
+ if (ntrk > 0) {
wgt = 1.;
+ }
foundelectrons += wgt;
peffElectrons.fill(p, wgt);
phieffElectrons.fill(phi, wgt);
ctheffElectrons.fill(cth, wgt);
-
// }
}
@@ -300,8 +341,9 @@
for (int layer = 1; layer < nlayers - 2; layer += 2) {
boolean hasThisLayer = false;
for (HelicalTrackCross hit : list) {
- if (hit.Layer() == layer)
+ if (hit.Layer() == layer) {
hasThisLayer = true;
+ }
}
if (!hasThisLayer) {
System.out.println("Missing reconstructed hit in layer = " + layer);
@@ -324,10 +366,12 @@
}
}
- if (!hasFitHitSL1)
+ if (!hasFitHitSL1) {
System.out.println("MISSING a hit in SL1!!!");
- if (!hasFitHitSL2)
+ }
+ if (!hasFitHitSL2) {
System.out.println("MISSING a hit in SL2!!!");
+ }
return false;
}
java/trunk/analysis/src/main/java/org/hps/analysis/examples
--- java/trunk/analysis/src/main/java/org/hps/analysis/examples/TrackAnalysis.java 2014-09-10 18:29:35 UTC (rev 983)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/examples/TrackAnalysis.java 2014-09-10 20:31:54 UTC (rev 984)
@@ -6,18 +6,32 @@
*/
package org.hps.analysis.examples;
+import hep.physics.matrix.BasicMatrix;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Matrix;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
import hep.physics.vec.VecOp;
-
import java.util.ArrayList;
import java.util.HashMap;
+import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Set;
-
+import java.util.TreeMap;
+import org.apache.commons.math3.linear.Array2DRowRealMatrix;
+import org.apache.commons.math3.linear.ArrayRealVector;
+import org.apache.commons.math3.linear.EigenDecomposition;
+import org.apache.commons.math3.linear.RealMatrix;
+import org.hps.conditions.deprecated.SvtUtils;
+import static org.hps.recon.tracking.CoordinateTransformations.transformVectorToTracking;
+import org.hps.recon.tracking.TrackerHitUtils;
+import org.lcsim.detector.identifier.IIdentifier;
+import org.lcsim.detector.identifier.Identifier;
import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
import org.lcsim.event.RelationalTable;
+import org.lcsim.event.SimTrackerHit;
import org.lcsim.event.Track;
import org.lcsim.event.TrackerHit;
import org.lcsim.fit.helicaltrack.HelicalTrack2DHit;
@@ -32,8 +46,11 @@
public class TrackAnalysis {
private enum HelixPar {
+
Curvature, Phi0, DCA, Z0, Slope
};
+ private static final Hep3Vector axial = new BasicHep3Vector(0, 1, 0);
+
private MCParticle _mcp = null;
private int _nhits;
private int _nbadhits;
@@ -47,25 +64,35 @@
private int _nbadAxialhits;
private int _nbadZhits;
private boolean _hasLayerOne;
- List<Integer> badHitList = new ArrayList();
- List<Integer> sharedHitList = new ArrayList();
- List<Integer> trackLayerList = new ArrayList();
- Map<MCParticle, HelicalTrackCross> badhits = new HashMap<MCParticle, HelicalTrackCross>();
- private int[] _nMCHitsPerLayer={0,0,0,0,0,0,0,0,0,0,0,0};
- private int[] _nStripHitsPerLayer={0,0,0,0,0,0,0,0,0,0,0,0};
- Map<Integer, Hep3Vector> _hitLocationPerLayer = new HashMap<Integer,Hep3Vector>();
+ private List<Integer> badHitList = new ArrayList();
+ private List<Integer> sharedHitList = new ArrayList();
+ private List<Integer> trackLayerList = new ArrayList();
+ private Map<MCParticle, HelicalTrackCross> badhits = new HashMap<MCParticle, HelicalTrackCross>();
+ // Create a map containing the number of hits for each MCParticle associated with the track
+ private Map<MCParticle, Integer> mcmap = new HashMap<MCParticle, Integer>();
+ private Map<MCParticle, Integer> mcmapAll = new HashMap<MCParticle, Integer>();
+ private Map<MCParticle, Integer> mcmapAxial = new HashMap<MCParticle, Integer>();
+ private Map<MCParticle, Integer> mcmapZ = new HashMap<MCParticle, Integer>();
+ private int[] _nMCHitsPerLayer = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ private int[] _nStripHitsPerLayer = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+ private Map<Integer, Hep3Vector> _hitLocationPerLayer = new HashMap<Integer, Hep3Vector>();
- /** Creates a new instance of TrackAnalysis */
+ /**
+ * Creates a new instance of TrackAnalysis
+ */
+ public TrackAnalysis(Track trk, RelationalTable hittomc, RelationalTable rthtosimhit, RelationalTable hittostrip, RelationalTable hittorotated) {
+ doAnalysis(trk, hittomc, rthtosimhit, hittostrip, hittorotated);
+ }
+
public TrackAnalysis(Track trk, RelationalTable hittomc) {
+ doAnalysis(trk, hittomc, null, null, null);
+ }
+ private void doAnalysis(Track trk, RelationalTable hittomc, RelationalTable rthtosimhit, RelationalTable hittostrip, RelationalTable hittorotated) {
+
// Get the number of hits on the track
_nhits = trk.getTrackerHits().size();
- // Create a map containing the number of hits for each MCParticle associated with the track
- Map<MCParticle, Integer> mcmap = new HashMap<MCParticle, Integer>();
- Map<MCParticle, Integer> mcmapAll = new HashMap<MCParticle, Integer>();
- Map<MCParticle, Integer> mcmapAxial = new HashMap<MCParticle, Integer>();
- Map<MCParticle, Integer> mcmapZ = new HashMap<MCParticle, Integer>();
_hasLayerOne = false;
// Loop over the hits on the track and make sure we have HelicalTrackHits (which contain the MC particle)
for (TrackerHit hit : trk.getTrackerHits()) {
@@ -73,79 +100,24 @@
Set<MCParticle> mclist = hittomc.allFrom(hit);
for (MCParticle mcp : mclist) {
Integer mchits = 0;
- if (mcmap.containsKey(mcp))
+ if (mcmap.containsKey(mcp)) {
mchits = mcmap.get(mcp);
+ }
mchits++;
mcmap.put(mcp, mchits);
}
- BasicHep3Vector axial = new BasicHep3Vector();
- axial.setV(0, 1, 0);
- HelicalTrackHit htc = (HelicalTrackHit) hit;
+// HelicalTrackHit htc = (HelicalTrackHit) hit;
if (hit instanceof HelicalTrackCross) {
- HelicalTrackCross cross = (HelicalTrackCross) hit;
- List<HelicalTrackStrip> clusterlist = cross.getStrips();
-
- for (HelicalTrackStrip cl : clusterlist) {
- int layer = cl.layer();
- if (layer == 1) _hasLayerOne = true;
-
- _nStripHitsPerLayer[layer - 1] = cl.rawhits().size();
- _hitLocationPerLayer.put(layer,clusterPosition(cl));
- _nhitsNew++;
- double axdotu = VecOp.dot(cl.u(), axial);
- boolean isAxial = false;
- if (axdotu > 0.5) {
- isAxial = true;
- _nAxialhits++;
- } else _nZhits++;
- List<MCParticle> mcPartList = cl.MCParticles();
- _nMCHitsPerLayer[layer-1] = mcPartList.size();
- for (MCParticle mcp : mcPartList) {
- Integer mchits = 0;
- if (mcmapAll.containsKey(mcp))
- mchits = mcmapAll.get(mcp);
- mchits++;
- mcmapAll.put(mcp, mchits);
- if (isAxial) {
- Integer mchitsAxial = 0;
- if (mcmapAxial.containsKey(mcp))
- mchitsAxial = mcmapAxial.get(mcp);
- mchitsAxial++;
- mcmapAxial.put(mcp, mchitsAxial);
- } else {
- Integer mchitsZ = 0;
- if (mcmapZ.containsKey(mcp))
- mchitsZ = mcmapZ.get(mcp);
- mchitsZ++;
- mcmapZ.put(mcp, mchitsZ);
- }
- }
- }
+ countHit((HelicalTrackCross) hit);
+ } else if (hit instanceof HelicalTrack2DHit) {
+ countHit((HelicalTrack2DHit) hit);
} else {
- _nhitsNew++;
- _nAxialhits++;
- HelicalTrack2DHit hit2d = (HelicalTrack2DHit) hit;
- List<MCParticle> mcPartList = hit2d.getMCParticles();
- //assume that lone hits are all axial
- boolean isAxial = true;
- for (MCParticle mcp : mcPartList) {
- Integer mchits = 0;
- if (mcmapAll.containsKey(mcp))
- mchits = mcmapAll.get(mcp);
- mchits++;
- mcmapAll.put(mcp, mchits);
- Integer mchitsAxial = 0;
- if (mcmapAxial.containsKey(mcp))
- mchitsAxial = mcmapAxial.get(mcp);
- mchitsAxial++;
- mcmapAxial.put(mcp, mchitsAxial);
- }
+ countHit(hit, rthtosimhit, hittostrip, hittorotated);
}
}
// Find the MCParticle that has the most hits on the track
-
int nbest = 0;
MCParticle mcbest = null;
for (MCParticle mcp : mcmap.keySet()) {
@@ -156,12 +128,12 @@
}
}
- if (nbest > 0)
+ if (nbest > 0) {
_mcp = mcbest;
+ }
_purity = (double) nbest / (double) _nhits;
_nbadhits = _nhits - nbest;
-
//single strip layer accounting.
int nbestAll = 0;
MCParticle mcbestAll = null;
@@ -173,42 +145,202 @@
}
}
- if (nbestAll > 0)
+ if (nbestAll > 0) {
_mcpNew = mcbestAll;
+ }
_purityNew = (double) nbestAll / (double) _nhitsNew;
_nbadhitsNew = _nhitsNew - nbestAll;
for (TrackerHit hit : trk.getTrackerHits()) {
- HelicalTrackHit htc = (HelicalTrackHit) hit;
if (hit instanceof HelicalTrackCross) {
- HelicalTrackCross cross = (HelicalTrackCross) hit;
- List<HelicalTrackStrip> clusterlist = cross.getStrips();
- for (HelicalTrackStrip cl : clusterlist){
- trackLayerList.add(cl.layer());
- if (!(cl.MCParticles().contains(_mcpNew))) {
- badHitList.add(cl.layer());
- badhits.put(_mcpNew, cross);
+ checkForBadHit((HelicalTrackCross) hit);
+ }
+ }
+
+ if (_nAxialhits > 0) {
+ if (mcmapAxial.containsKey(_mcpNew)) {
+ _nbadAxialhits = _nAxialhits - mcmapAxial.get(_mcpNew);
+ } else {
+ _nbadAxialhits = _nAxialhits;
+ }
+ }
+ if (_nZhits > 0) {
+ if (mcmapZ.containsKey(_mcpNew)) {
+ _nbadZhits = _nZhits - mcmapZ.get(_mcpNew);
+ } else {
+ _nbadZhits = _nZhits;
+ }
+ }
+ }
+
+ private void countHit(HelicalTrackCross cross) {
+ List<HelicalTrackStrip> clusterlist = cross.getStrips();
+
+ for (HelicalTrackStrip cl : clusterlist) {
+ int layer = cl.layer();
+ if (layer == 1) {
+ _hasLayerOne = true;
+ }
+
+ _nStripHitsPerLayer[layer - 1] = cl.rawhits().size();
+ _hitLocationPerLayer.put(layer, clusterPosition(cl));
+ _nhitsNew++;
+ double axdotu = VecOp.dot(cl.u(), axial);
+// System.out.println(new BasicHep3Vector(cross.getPosition()).toString() + cl.u());
+ boolean isAxial = false;
+ if (axdotu > 0.5) {
+ isAxial = true;
+ _nAxialhits++;
+ } else {
+ _nZhits++;
+ }
+ List<MCParticle> mcPartList = cl.MCParticles();
+ _nMCHitsPerLayer[layer - 1] = mcPartList.size();
+ for (MCParticle mcp : mcPartList) {
+ Integer mchits = 0;
+ if (mcmapAll.containsKey(mcp)) {
+ mchits = mcmapAll.get(mcp);
+ }
+ mchits++;
+ mcmapAll.put(mcp, mchits);
+ if (isAxial) {
+ Integer mchitsAxial = 0;
+ if (mcmapAxial.containsKey(mcp)) {
+ mchitsAxial = mcmapAxial.get(mcp);
}
- if(cl.MCParticles().size()>1)
- sharedHitList.add(cl.layer());
+ mchitsAxial++;
+ mcmapAxial.put(mcp, mchitsAxial);
+ } else {
+ Integer mchitsZ = 0;
+ if (mcmapZ.containsKey(mcp)) {
+ mchitsZ = mcmapZ.get(mcp);
+ }
+ mchitsZ++;
+ mcmapZ.put(mcp, mchitsZ);
}
}
}
+ }
+ private void countHit(TrackerHit hit, RelationalTable rthtosimhit, RelationalTable hittostrip, RelationalTable hittorotated) {
+ TrackerHit unrotatedHit = (TrackerHit) hittorotated.from(hit);
+// System.out.println("ID: " + unrotatedHit.getCellID());
+ Set<TrackerHit> hitlist = hittostrip.allFrom(unrotatedHit);
+// System.out.println("size: " + hitlist.size());
+ for (TrackerHit cl : hitlist) {
+ int layer = -1;
+ int module = -1;
+ List<RawTrackerHit> rawHits = cl.getRawHits();
+// System.out.println("RawHits: " + rawHits.size());
+ for (RawTrackerHit rawHit : rawHits) {
+// System.out.println(rawHit.getCellID());
+ IIdentifier id = new Identifier(rawHit.getCellID());
+ int newLayer = SvtUtils.getInstance().getHelper().getValue(id, "layer");
+ if (layer != -1 && layer != newLayer) {
+ System.out.format("TrackerHit has hits from multiple layers: %d and %d\n", layer, newLayer);
+ }
+ layer = newLayer;
+ int newModule = SvtUtils.getInstance().getHelper().getValue(id, "module");
+ if (module != -1 && module != newModule) {
+ System.out.format("TrackerHit has hits from multiple modules: %d and %d\n", module, newModule);
+ }
+ module = newModule;
+// System.out.println(SvtUtils.getInstance().getHelper().getValue(id, "strip"));
+ }
+ if (layer == 1) {
+ _hasLayerOne = true;
+ }
+ DiagonalizedCovarianceMatrix covariance = new DiagonalizedCovarianceMatrix(cl);
+ _nStripHitsPerLayer[layer - 1] = cl.getRawHits().size();
+ _hitLocationPerLayer.put(layer, new BasicHep3Vector(hit.getPosition()));
+ _nhitsNew++;
- if (_nAxialhits > 0)
- if (mcmapAxial.containsKey(_mcpNew))
- _nbadAxialhits = _nAxialhits - mcmapAxial.get(_mcpNew);
- else _nbadAxialhits = _nAxialhits;
- if (_nZhits > 0)
- if (mcmapZ.containsKey(_mcpNew))
- _nbadZhits = _nZhits - mcmapZ.get(_mcpNew);
- else _nbadZhits = _nZhits;
+ double axdotu = VecOp.dot(transformVectorToTracking(covariance.getMeasuredVector()), axial);
+// System.out.println(transformVectorToTracking(new BasicHep3Vector(cl.getPosition())).toString() + transformVectorToTracking(covariance.getMeasuredVector()));
+ boolean isAxial = false;
+ if (axdotu > 0.5) {
+ isAxial = true;
+ _nAxialhits++;
+ } else {
+ _nZhits++;
+ }
+ // get the set of MCParticles associated with this hit and update the hit count for each MCParticle
+ Set<MCParticle> mcPartList = new HashSet<MCParticle>();
+ for (RawTrackerHit rawHit : rawHits) {
+ Set<SimTrackerHit> simhits = (Set<SimTrackerHit>) rthtosimhit.allFrom(rawHit);
+ for (SimTrackerHit simhit : simhits) {
+ if (simhit != null && simhit.getMCParticle() != null) {
+ mcPartList.add(simhit.getMCParticle());
+ }
+ }
+ }
+// System.out.println("MCParticle count: " + mcPartList.size());
+ _nMCHitsPerLayer[layer - 1] = mcPartList.size();
+ for (MCParticle mcp : mcPartList) {
+ Integer mchits = 0;
+ if (mcmapAll.containsKey(mcp)) {
+ mchits = mcmapAll.get(mcp);
+ }
+ mchits++;
+ mcmapAll.put(mcp, mchits);
+ if (isAxial) {
+ Integer mchitsAxial = 0;
+ if (mcmapAxial.containsKey(mcp)) {
+ mchitsAxial = mcmapAxial.get(mcp);
+ }
+ mchitsAxial++;
+ mcmapAxial.put(mcp, mchitsAxial);
+ } else {
+ Integer mchitsZ = 0;
+ if (mcmapZ.containsKey(mcp)) {
+ mchitsZ = mcmapZ.get(mcp);
+ }
+ mchitsZ++;
+ mcmapZ.put(mcp, mchitsZ);
+ }
+ }
+ }
}
- public Hep3Vector clusterPosition(HelicalTrackStrip cl) {
+ private void countHit(HelicalTrack2DHit hit2d) {
+ _nhitsNew++;
+ _nAxialhits++;
+ List<MCParticle> mcPartList = hit2d.getMCParticles();
+ //assume that lone hits are all axial
+ boolean isAxial = true;
+ for (MCParticle mcp : mcPartList) {
+ Integer mchits = 0;
+ if (mcmapAll.containsKey(mcp)) {
+ mchits = mcmapAll.get(mcp);
+ }
+ mchits++;
+ mcmapAll.put(mcp, mchits);
+ Integer mchitsAxial = 0;
+ if (mcmapAxial.containsKey(mcp)) {
+ mchitsAxial = mcmapAxial.get(mcp);
+ }
+ mchitsAxial++;
+ mcmapAxial.put(mcp, mchitsAxial);
+ }
+ }
+
+ private void checkForBadHit(HelicalTrackCross cross) {
+ List<HelicalTrackStrip> clusterlist = cross.getStrips();
+ for (HelicalTrackStrip cl : clusterlist) {
+ trackLayerList.add(cl.layer());
+ if (!(cl.MCParticles().contains(_mcpNew))) {
+ badHitList.add(cl.layer());
+ badhits.put(_mcpNew, cross);
+ }
+ if (cl.MCParticles().size() > 1) {
+ sharedHitList.add(cl.layer());
+ }
+ }
+ }
+
+ public static Hep3Vector clusterPosition(HelicalTrackStrip cl) {
Hep3Vector corigin = cl.origin();
Hep3Vector u = cl.u();
double umeas = cl.umeas();
@@ -284,15 +416,97 @@
public List<Integer> getBadHitList() {
return badHitList;
}
- public List<Integer> getSharedHitList() {
+
+ public List<Integer> getSharedHitList() {
return sharedHitList;
}
-
- public List<Integer> getTrackLayerList() {
+
+ public List<Integer> getTrackLayerList() {
return trackLayerList;
}
public Map<MCParticle, HelicalTrackCross> getBadHits() {
return badhits;
}
+
+ public static class DiagonalizedCovarianceMatrix {
+
+ double[] measurement_errors = new double[3];
+ Hep3Vector[] measurement_vectors = new Hep3Vector[3];
+
+ public DiagonalizedCovarianceMatrix(TrackerHit hit) {
+ SymmetricMatrix cov = new SymmetricMatrix(3, hit.getCovMatrix(), true);
+ RealMatrix covMatrix = new Array2DRowRealMatrix(3, 3);
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ covMatrix.setEntry(i, j, cov.e(i, j));
+ }
+ }
+ EigenDecomposition decomposed = new EigenDecomposition(covMatrix);
+ BasicHep3Matrix localToGlobal = new BasicHep3Matrix();
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ localToGlobal.setElement(i, j, decomposed.getV().getEntry(i, j));
+ }
+ }
+// SymmetricMatrix localToGlobal = decomposed.getV().operate(new ArrayRealVector(3))
+ {
+ double eigenvalue = decomposed.getRealEigenvalue(0);
+// Hep3Vector eigenvector = VecOp.mult(localToGlobal, new BasicHep3Vector());
+ Hep3Vector eigenvector = VecOp.mult(Math.signum(eigenvalue), new BasicHep3Vector(decomposed.getVT().getRow(0)));
+ measurement_errors[0] = eigenvalue;
+ measurement_vectors[0] = eigenvector;
+ measurement_errors[2] = eigenvalue;
+ measurement_vectors[2] = eigenvector;
+ }
+ {
+ double eigenvalue = decomposed.getRealEigenvalue(1);
+ Hep3Vector eigenvector = VecOp.mult(Math.signum(eigenvalue), new BasicHep3Vector(decomposed.getVT().getRow(1)));
+ if (eigenvalue > measurement_errors[0]) {
+ measurement_errors[0] = eigenvalue;
+ measurement_vectors[0] = eigenvector;
+ }
+ if (eigenvalue < measurement_errors[2]) {
+ measurement_errors[2] = eigenvalue;
+ measurement_vectors[2] = eigenvector;
+ }
+ }
+ {
+ double eigenvalue = decomposed.getRealEigenvalue(2);
+ Hep3Vector eigenvector = VecOp.mult(Math.signum(eigenvalue), new BasicHep3Vector(decomposed.getVT().getRow(2)));
+ if (eigenvalue > measurement_errors[0]) {
+ measurement_errors[1] = measurement_errors[0];
+ measurement_vectors[1] = measurement_vectors[0];
+ measurement_errors[0] = eigenvalue;
+ measurement_vectors[0] = eigenvector;
+ }
+ if (eigenvalue < measurement_errors[2]) {
+ measurement_errors[1] = measurement_errors[2];
+ measurement_vectors[1] = measurement_vectors[2];
+ measurement_errors[2] = eigenvalue;
+ measurement_vectors[2] = eigenvector;
+ }
+ if (measurement_vectors[1] == null) {
+ measurement_errors[1] = eigenvalue;
+ measurement_vectors[1] = eigenvector;
+ }
+ }
+// for (int i = 0; i < 3; i++) {
+// System.out.format("%d: resolution %f, vector %s\n", i, measurement_errors[i], measurement_vectors[i].toString());
+// }
+ }
+
+ public Hep3Vector getUnmeasuredVector() {
+ return measurement_vectors[0];
+ }
+
+ public Hep3Vector getMeasuredVector() {
+ return measurement_vectors[1];
+ }
+
+ public Hep3Vector getNormalVector() {
+ return measurement_vectors[2];
+ }
+
+ }
}