Skip to content

Commit 75655dd

Browse files
bloodyyugobloodyyugo
authored
Merge pull request #1175 from mgignac/mgignac-lowp-tracking
Kalman per-iteration hit timing cuts, plus SVT truth-relation diagnostics
2 parents ec09f1b + 4c45eec commit 75655dd

12 files changed

Lines changed: 1650 additions & 33 deletions

File tree

‎analysis/src/main/java/org/hps/analysis/MC/SvtEventForensicsDriver.java‎

Lines changed: 397 additions & 0 deletions
Large diffs are not rendered by default.

‎analysis/src/main/java/org/hps/analysis/MC/SvtHitProvenanceDriver.java‎

Lines changed: 297 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,297 @@
1+
package org.hps.analysis.MC;
2+
3+
import java.util.ArrayList;
4+
import java.util.HashMap;
5+
import java.util.HashSet;
6+
import java.util.List;
7+
import java.util.Map;
8+
import java.util.Set;
9+
import java.util.TreeSet;
10+
11+
import org.lcsim.event.EventHeader;
12+
import org.lcsim.event.LCRelation;
13+
import org.lcsim.event.RawTrackerHit;
14+
import org.lcsim.event.Track;
15+
import org.lcsim.event.TrackerHit;
16+
import org.lcsim.util.Driver;
17+
18+
/**
19+
* Attributes the hits on each track to their origin, using the provenance relations
20+
* written by SvtDigitizationWithPulserDataMergingReadoutDriver.
21+
*
22+
* Under pulser overlay a large fraction of tracks carry no truth relation at all. Two
23+
* explanations survive the aggregate hit counts:
24+
*
25+
* 1. The tracks are fakes assembled from pulser data hits.
26+
* 2. The tracks follow a real trajectory but were built from pulser hits on strips
27+
* neighbouring the ones the MC particle actually hit. Adjacent strips are ~55 um
28+
* apart, which over a metre of lever arm is ~55 urad, so such a track would still
29+
* point at the true particle to well inside a milliradian while carrying zero
30+
* truth relations.
31+
*
32+
* These differ observably in whether the untruthed hits on a track sit next to channels
33+
* that did receive MC charge. That is what this driver measures, against the baseline
34+
* rate at which any pulser hit in the event happens to be near an MC channel -- with
35+
* high occupancy, adjacency alone proves nothing, so the comparison to the baseline is
36+
* the whole point.
37+
*
38+
* Requires writeHitOriginCollections=true on the digitization driver.
39+
*/
40+
public class SvtHitProvenanceDriver extends Driver {
41+
42+
private String trackCollectionName = "KalmanFullTracks";
43+
private String rawHitCollectionName = "SVTRawTrackerHits";
44+
private String truthRelationCollectionName = "SVTTrueHitRelations";
45+
private String pulserOriginCollectionName = "SVTHitOriginPulser";
46+
private String mcContribCollectionName = "SVTHitOriginMCContrib";
47+
48+
/** Neighbour distances, in strips, at which adjacency is reported. */
49+
private static final int[] DISTANCES = { 1, 2, 3, 5, 10 };
50+
/** Bins for the minimum strip distance to an MC-contributing channel. */
51+
private static final int MAX_DIST_BIN = 12;
52+
53+
private boolean debug = false;
54+
private int debugMaxPrint = 20;
55+
private int debugPrinted = 0;
56+
57+
// Hit categories, matching the digitization driver.
58+
private static final int NOISE = 0;
59+
private static final int MC_PURE = 1;
60+
private static final int MC_PURE_SUBTHRESH = 2;
61+
private static final int PULSER_PURE = 3;
62+
private static final int MERGED = 4;
63+
private static final int MERGED_SUBTHRESH = 5;
64+
private static final String[] CAT_NAME = {
65+
"NOISE ", "MC_PURE ", "MC_PURE_SUBTHRESH",
66+
"PULSER_PURE ", "MERGED ", "MERGED_SUBTHRESH "
67+
};
68+
69+
private long nEvents = 0;
70+
private long nTracks = 0;
71+
private long nZeroTruthTracks = 0;
72+
private boolean warnedMissing = false;
73+
74+
// [0] = tracks with at least one truthed hit, [1] = tracks with none
75+
private final long[] nTracksByClass = new long[2];
76+
private final long[] nHitsByClass = new long[2];
77+
private final long[][] nHitsByClassAndCat = new long[2][6];
78+
79+
// Minimum strip distance from a PULSER_PURE hit to an MC-contributing channel on the
80+
// same sensor. Index MAX_DIST_BIN is the overflow, MAX_DIST_BIN+1 means the sensor had
81+
// no MC-contributing channel at all.
82+
private final long[][] distHistByClass = new long[2][MAX_DIST_BIN + 2];
83+
private final long[] distHistBaseline = new long[MAX_DIST_BIN + 2];
84+
private long nPulserPureBaseline = 0;
85+
86+
public void setTrackCollectionName(String val) { this.trackCollectionName = val; }
87+
public void setRawHitCollectionName(String val) { this.rawHitCollectionName = val; }
88+
public void setTruthRelationCollectionName(String val) { this.truthRelationCollectionName = val; }
89+
public void setPulserOriginCollectionName(String val) { this.pulserOriginCollectionName = val; }
90+
public void setMcContribCollectionName(String val) { this.mcContribCollectionName = val; }
91+
public void setDebug(boolean val) { this.debug = val; }
92+
public void setDebugMaxPrint(int val) { this.debugMaxPrint = val; }
93+
94+
/** Collects the "from" side of a relation collection, tolerating a missing collection. */
95+
private Set<RawTrackerHit> fromSide(EventHeader event, String name) {
96+
Set<RawTrackerHit> out = new HashSet<RawTrackerHit>();
97+
if(!event.hasCollection(LCRelation.class, name)) { return out; }
98+
for(LCRelation rel : event.get(LCRelation.class, name)) {
99+
if(rel.getFrom() instanceof RawTrackerHit) {
100+
out.add((RawTrackerHit) rel.getFrom());
101+
}
102+
}
103+
return out;
104+
}
105+
106+
private static String sensorOf(RawTrackerHit hit) {
107+
return hit.getDetectorElement().getName();
108+
}
109+
110+
private static int channelOf(RawTrackerHit hit) {
111+
return hit.getIdentifierFieldValue("strip");
112+
}
113+
114+
private static int categorise(RawTrackerHit hit, Set<RawTrackerHit> pulser,
115+
Set<RawTrackerHit> mcContrib, Set<RawTrackerHit> truthed) {
116+
boolean p = pulser.contains(hit);
117+
boolean m = mcContrib.contains(hit);
118+
boolean t = truthed.contains(hit);
119+
if(!m) { return p ? PULSER_PURE : NOISE; }
120+
if(p) { return t ? MERGED : MERGED_SUBTHRESH; }
121+
return t ? MC_PURE : MC_PURE_SUBTHRESH;
122+
}
123+
124+
/**
125+
* Minimum distance in strips from this hit to a channel that received MC charge on the
126+
* same sensor. Returns MAX_DIST_BIN+1 if the sensor had no MC contribution anywhere.
127+
*/
128+
private int minDistanceToMC(RawTrackerHit hit, Map<String, TreeSet<Integer>> mcChannels) {
129+
TreeSet<Integer> chans = mcChannels.get(sensorOf(hit));
130+
if(chans == null || chans.isEmpty()) { return MAX_DIST_BIN + 1; }
131+
int ch = channelOf(hit);
132+
Integer lo = chans.floor(ch);
133+
Integer hi = chans.ceiling(ch);
134+
int best = Integer.MAX_VALUE;
135+
if(lo != null) { best = Math.min(best, ch - lo); }
136+
if(hi != null) { best = Math.min(best, hi - ch); }
137+
if(best == Integer.MAX_VALUE) { return MAX_DIST_BIN + 1; }
138+
return Math.min(best, MAX_DIST_BIN);
139+
}
140+
141+
@Override
142+
public void process(EventHeader event) {
143+
if(!event.hasCollection(Track.class, trackCollectionName)) { return; }
144+
145+
Set<RawTrackerHit> truthed = fromSide(event, truthRelationCollectionName);
146+
Set<RawTrackerHit> pulser = fromSide(event, pulserOriginCollectionName);
147+
Set<RawTrackerHit> mcContrib = fromSide(event, mcContribCollectionName);
148+
149+
if(!warnedMissing && !event.hasCollection(LCRelation.class, mcContribCollectionName)) {
150+
warnedMissing = true;
151+
System.out.println("SvtHitProvenanceDriver: WARNING collection '" + mcContribCollectionName
152+
+ "' not found. Was writeHitOriginCollections set on the digitization driver?");
153+
}
154+
155+
nEvents++;
156+
157+
// Channels that received MC charge, by sensor. Built from the MC-contribution
158+
// relation so it is independent of whether the truth gate kept the relation.
159+
Map<String, TreeSet<Integer>> mcChannels = new HashMap<String, TreeSet<Integer>>();
160+
for(RawTrackerHit hit : mcContrib) {
161+
String s = sensorOf(hit);
162+
TreeSet<Integer> set = mcChannels.get(s);
163+
if(set == null) { set = new TreeSet<Integer>(); mcChannels.put(s, set); }
164+
set.add(channelOf(hit));
165+
}
166+
167+
// Baseline: how near an MC channel does an arbitrary pulser hit in this event sit?
168+
// Tracks are compared against this, not against zero.
169+
if(event.hasCollection(RawTrackerHit.class, rawHitCollectionName)) {
170+
for(RawTrackerHit hit : event.get(RawTrackerHit.class, rawHitCollectionName)) {
171+
if(categorise(hit, pulser, mcContrib, truthed) != PULSER_PURE) { continue; }
172+
nPulserPureBaseline++;
173+
distHistBaseline[minDistanceToMC(hit, mcChannels)]++;
174+
}
175+
}
176+
177+
for(Track track : event.get(Track.class, trackCollectionName)) {
178+
nTracks++;
179+
180+
List<RawTrackerHit> rawHits = new ArrayList<RawTrackerHit>();
181+
for(TrackerHit th : track.getTrackerHits()) {
182+
for(Object o : th.getRawHits()) {
183+
if(o instanceof RawTrackerHit) { rawHits.add((RawTrackerHit) o); }
184+
}
185+
}
186+
if(rawHits.isEmpty()) { continue; }
187+
188+
int nTruthedOnTrack = 0;
189+
for(RawTrackerHit hit : rawHits) {
190+
if(truthed.contains(hit)) { nTruthedOnTrack++; }
191+
}
192+
final int cls = (nTruthedOnTrack == 0) ? 1 : 0;
193+
if(cls == 1) { nZeroTruthTracks++; }
194+
nTracksByClass[cls]++;
195+
nHitsByClass[cls] += rawHits.size();
196+
197+
for(RawTrackerHit hit : rawHits) {
198+
int cat = categorise(hit, pulser, mcContrib, truthed);
199+
nHitsByClassAndCat[cls][cat]++;
200+
if(cat == PULSER_PURE) {
201+
distHistByClass[cls][minDistanceToMC(hit, mcChannels)]++;
202+
}
203+
}
204+
205+
if(debug && cls == 1 && debugPrinted < debugMaxPrint) {
206+
debugPrinted++;
207+
StringBuilder sb = new StringBuilder("[SvtProv] zero-truth track, nRawHits="
208+
+ rawHits.size() + " hits:");
209+
for(RawTrackerHit hit : rawHits) {
210+
int cat = categorise(hit, pulser, mcContrib, truthed);
211+
sb.append(" ").append(CAT_NAME[cat].trim())
212+
.append("(").append(sensorOf(hit)).append(":").append(channelOf(hit));
213+
if(cat == PULSER_PURE) {
214+
int d = minDistanceToMC(hit, mcChannels);
215+
sb.append(", dMC=").append(d > MAX_DIST_BIN ? "none" : Integer.toString(d));
216+
}
217+
sb.append(")");
218+
}
219+
System.out.println(sb.toString());
220+
}
221+
}
222+
}
223+
224+
/** Fraction of entries in a distance histogram at or below d strips. */
225+
private static double fracWithin(long[] hist, int d) {
226+
long num = 0, den = 0;
227+
for(int i = 0; i < hist.length; i++) {
228+
den += hist[i];
229+
if(i <= d) { num += hist[i]; }
230+
}
231+
return den > 0 ? (double) num / den : 0.0;
232+
}
233+
234+
@Override
235+
public void endOfData() {
236+
String[] clsName = { "tracks with truth ", "ZERO-truth tracks " };
237+
System.out.println();
238+
System.out.println("============== SVT hit provenance by track ==============");
239+
System.out.println(" events : " + nEvents);
240+
System.out.println(" tracks (" + trackCollectionName + ") : " + nTracks);
241+
System.out.println(" zero-truth tracks : " + nZeroTruthTracks
242+
+ (nTracks > 0 ? String.format(" (%.4f)", (double) nZeroTruthTracks / nTracks) : ""));
243+
System.out.println();
244+
245+
for(int cls = 0; cls < 2; cls++) {
246+
System.out.println(" ---- " + clsName[cls] + " ----");
247+
System.out.println(" tracks : " + nTracksByClass[cls]
248+
+ " raw hits/track : "
249+
+ (nTracksByClass[cls] > 0
250+
? String.format("%.2f", (double) nHitsByClass[cls] / nTracksByClass[cls]) : "-"));
251+
for(int cat = 0; cat < 6; cat++) {
252+
long n = nHitsByClassAndCat[cls][cat];
253+
System.out.println(" " + CAT_NAME[cat] + " : " + n
254+
+ (nHitsByClass[cls] > 0
255+
? String.format(" (%.4f)", (double) n / nHitsByClass[cls]) : ""));
256+
}
257+
System.out.println();
258+
}
259+
260+
System.out.println(" ---- are PULSER_PURE hits next to channels that saw MC charge? ----");
261+
System.out.println(" The baseline is every PULSER_PURE hit in the event, so it already");
262+
System.out.println(" folds in the occupancy. Only an excess over baseline is meaningful.");
263+
StringBuilder hdr = new StringBuilder(String.format(" %-26s %10s", "population", "nHits"));
264+
for(int d : DISTANCES) { hdr.append(String.format(" within%3d", d)); }
265+
System.out.println(hdr.toString());
266+
267+
long nZ = 0, nT = 0;
268+
for(int i = 0; i < distHistByClass[1].length; i++) { nZ += distHistByClass[1][i]; }
269+
for(int i = 0; i < distHistByClass[0].length; i++) { nT += distHistByClass[0][i]; }
270+
271+
Object[][] rows = {
272+
{ "baseline (all in event)", distHistBaseline, nPulserPureBaseline },
273+
{ "on tracks with truth", distHistByClass[0], nT },
274+
{ "on ZERO-truth tracks", distHistByClass[1], nZ },
275+
};
276+
for(Object[] row : rows) {
277+
StringBuilder sb = new StringBuilder(String.format(" %-26s %10d", row[0], (Long) row[2]));
278+
for(int d : DISTANCES) {
279+
sb.append(String.format(" %8.4f", fracWithin((long[]) row[1], d)));
280+
}
281+
System.out.println(sb.toString());
282+
}
283+
284+
System.out.println();
285+
System.out.println(" minimum strip distance to an MC channel, ZERO-truth tracks:");
286+
StringBuilder sb = new StringBuilder(" ");
287+
for(int i = 0; i <= MAX_DIST_BIN; i++) {
288+
sb.append(i == MAX_DIST_BIN ? ">=" + MAX_DIST_BIN : Integer.toString(i))
289+
.append("=").append(distHistByClass[1][i]).append(" ");
290+
}
291+
sb.append("noMCOnSensor=").append(distHistByClass[1][MAX_DIST_BIN + 1]);
292+
System.out.println(sb.toString());
293+
System.out.println("=========================================================");
294+
System.out.println();
295+
super.endOfData();
296+
}
297+
}

‎detector-model/dependency-reduced-pom.xml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@
5151
<dependency>
5252
<groupId>junit</groupId>
5353
<artifactId>junit</artifactId>
54-
<version>4.13.1</version>
54+
<version>4.13.2</version>
5555
<scope>test</scope>
5656
<exclusions>
5757
<exclusion>

0 commit comments

Comments
 (0)