Skip to content

Commit 05ab92f

Browse files
JeremyMcCormickJeremyMcCormick
authored
Merge pull request #220 from JeffersonLab/iss83
Iss83: Implement reconstruction ITs
2 parents aa5692b + b9deed2 commit 05ab92f

14 files changed

Lines changed: 2886 additions & 104 deletions

‎integration-tests/src/test/java/org/hps/test/it/EngRun2015FeeRecon.java‎

Lines changed: 311 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,311 @@
1+
package org.hps.test.it;
2+
3+
import hep.aida.IHistogram1D;
4+
import hep.physics.vec.BasicHep3Matrix;
5+
import hep.physics.vec.Hep3Vector;
6+
import hep.physics.vec.VecOp;
7+
import java.io.File;
8+
import java.io.IOException;
9+
import java.util.List;
10+
import org.hps.recon.ecal.cluster.ClusterUtilities;
11+
import org.hps.recon.tracking.TrackData;
12+
import org.hps.recon.tracking.TrackType;
13+
import org.hps.recon.tracking.TrackUtils;
14+
import org.hps.record.triggerbank.AbstractIntData;
15+
import org.hps.record.triggerbank.TIData;
16+
import org.lcsim.detector.DetectorElementStore;
17+
import org.lcsim.detector.IDetectorElement;
18+
import org.lcsim.detector.identifier.IExpandedIdentifier;
19+
import org.lcsim.detector.identifier.IIdentifier;
20+
import org.lcsim.detector.identifier.IIdentifierDictionary;
21+
import org.lcsim.detector.tracker.silicon.HpsSiSensor;
22+
import org.lcsim.detector.tracker.silicon.SiSensor;
23+
import org.lcsim.event.Cluster;
24+
import org.lcsim.event.EventHeader;
25+
import org.lcsim.event.GenericObject;
26+
import org.lcsim.event.RawTrackerHit;
27+
import org.lcsim.event.ReconstructedParticle;
28+
import org.lcsim.event.Track;
29+
import org.lcsim.event.TrackerHit;
30+
import org.lcsim.geometry.Detector;
31+
import org.lcsim.math.chisq.ChisqProb;
32+
import org.lcsim.util.Driver;
33+
import org.lcsim.util.aida.AIDA;
34+
35+
/**
36+
* Analysis Driver to test reconstruction of 2015 Full Energy Electron candidates.
37+
* All selection cuts should have already been made prior to calling this Driver.
38+
*
39+
* @author Norman A. Graf
40+
*/
41+
public class EngRun2015FeeRecon extends Driver {
42+
43+
boolean debug = true;
44+
private AIDA aida = AIDA.defaultInstance();
45+
private IHistogram1D trkChisqNdfTop = aida.histogram1D("Top Track Chisq per DoF", 100, 0., 100.);
46+
private IHistogram1D trkChisqProbTop = aida.histogram1D("Top Track Chisq Prob", 100, 0., 1.);
47+
private IHistogram1D trkNhitsTop = aida.histogram1D("Top Track Number of Hits", 7, -0.5, 6.5);
48+
private IHistogram1D trkMomentumTop = aida.histogram1D("Top Track Momentum", 200, 0., 3.);
49+
private IHistogram1D trkMomentumTop5 = aida.histogram1D("Top 5 Hit Track Momentum", 200, 0., 3.);
50+
private IHistogram1D trkMomentumTop6 = aida.histogram1D("Top 6 Hit Track Momentum", 200, 0., 3.);
51+
private IHistogram1D trkdEdXTop5 = aida.histogram1D("Top 5 Track dEdx", 100, 0., .00015);
52+
private IHistogram1D trkdEdXTop6 = aida.histogram1D("Top 6 Track dEdx", 100, 0., .00015);
53+
private IHistogram1D trkthetaTop = aida.histogram1D("Top Track theta", 100, 0.01, 0.05);
54+
private IHistogram1D trkX0Top = aida.histogram1D("Top Track X0", 100, -0.5, 0.5);
55+
private IHistogram1D trkY0Top = aida.histogram1D("Top Track Y0", 100, -5.0, 5.0);
56+
private IHistogram1D trkZ0Top = aida.histogram1D("Top Track Z0", 100, -1.0, 1.0);
57+
58+
private IHistogram1D trkChisqNdfBottom = aida.histogram1D("Bottom Track Chisq per DoF", 100, 0., 100.);
59+
private IHistogram1D trkChisqProbBottom = aida.histogram1D("Bottom Track Chisq Prob", 100, 0., 1.);
60+
private IHistogram1D trkNhitsBottom = aida.histogram1D("Bottom Track Number of Hits", 7, -0.5, 6.5);
61+
private IHistogram1D trkMomentumBottom = aida.histogram1D("Bottom Track Momentum", 200, 0., 3.);
62+
private IHistogram1D trkMomentumBottom5 = aida.histogram1D("Bottom 5 Hit Track Momentum", 200, 0., 3.);
63+
private IHistogram1D trkMomentumBottom6 = aida.histogram1D("Bottom 6 Hit Track Momentum", 200, 0., 3.);
64+
private IHistogram1D trkdEdXBottom5 = aida.histogram1D("Bottom 5 Track dEdx", 100, 0., .00015);
65+
private IHistogram1D trkdEdXBottom6 = aida.histogram1D("Bottom 6 Track dEdx", 100, 0., .00015);
66+
private IHistogram1D trkthetaBottom = aida.histogram1D("Bottom Track theta", 100, 0.01, 0.05);
67+
private IHistogram1D trkX0Bottom = aida.histogram1D("Bottom Track X0", 100, -0.5, 0.5);
68+
private IHistogram1D trkY0Bottom = aida.histogram1D("Bottom Track Y0", 100, -5.0, 5.0);
69+
private IHistogram1D trkZ0Bottom = aida.histogram1D("Bottom Track Z0", 100, -1.0, 1.0);
70+
71+
private String _aidaFileName = "EngRun2015FeeRecon.aida";
72+
private final String finalStateParticlesColName = "FinalStateParticles";
73+
74+
private Double _beamEnergy = 1.056;
75+
private double _percentFeeCut = 0.8;
76+
private final BasicHep3Matrix beamAxisRotation = new BasicHep3Matrix();
77+
78+
//Set min seed energy value, default to 2015 run
79+
private double seedCut = 0.4; //= 0.4
80+
81+
//set min cluster energy value, default to 2015 run
82+
private double clusterCut = 0.6;
83+
84+
//minimum number of hits per cluster
85+
private int minHits = 3; // = 3;
86+
87+
double ctMin = 40.;
88+
double ctMax = 49.;
89+
90+
private boolean _dumpRunAndEventNumber = false;
91+
92+
protected void detectorChanged(Detector detector) {
93+
beamAxisRotation.setActiveEuler(Math.PI / 2, -0.0305, -Math.PI / 2);
94+
}
95+
96+
protected void process(EventHeader event) {
97+
if (event.getRunNumber() > 7000) {
98+
_beamEnergy = 2.306;
99+
seedCut = 1.2;
100+
clusterCut = 1.6;
101+
ctMin = 55.;
102+
ctMax = 61.;
103+
}
104+
// only keep singles triggers:
105+
if (!event.hasCollection(GenericObject.class, "TriggerBank")) {
106+
return;
107+
}
108+
boolean isSingles = false;
109+
for (GenericObject gob : event.get(GenericObject.class, "TriggerBank")) {
110+
if (!(AbstractIntData.getTag(gob) == TIData.BANK_TAG)) {
111+
continue;
112+
}
113+
TIData tid = new TIData(gob);
114+
if (tid.isSingle0Trigger() || tid.isSingle1Trigger()) {
115+
isSingles = true;
116+
break;
117+
}
118+
}
119+
if (!isSingles) {
120+
return;
121+
}
122+
if (!event.hasCollection(ReconstructedParticle.class, finalStateParticlesColName)) {
123+
return;
124+
}
125+
List<ReconstructedParticle> rpList = event.get(ReconstructedParticle.class, finalStateParticlesColName);
126+
setupSensors(event);
127+
for (ReconstructedParticle rp : rpList) {
128+
if (!TrackType.isGBL(rp.getType())) {
129+
continue;
130+
}
131+
if (rp.getMomentum().magnitude() > 1.5 * _beamEnergy) {
132+
continue;
133+
}
134+
// require both track and cluster
135+
if (rp.getClusters().size() != 1) {
136+
continue;
137+
}
138+
139+
if (rp.getTracks().size() != 1) {
140+
continue;
141+
}
142+
143+
Track t = rp.getTracks().get(0);
144+
double p = rp.getMomentum().magnitude();
145+
Cluster c = rp.getClusters().get(0);
146+
// debug diagnostics to set cuts
147+
if (debug) {
148+
aida.cloud1D("clusterSeedHit energy").fill(ClusterUtilities.findSeedHit(c).getCorrectedEnergy());
149+
aida.cloud1D("cluster nHits").fill(c.getCalorimeterHits().size());
150+
aida.cloud2D("clusterSeedHit energy vs p").fill(p, ClusterUtilities.findSeedHit(c).getCorrectedEnergy());
151+
aida.cloud2D("cluster nHits vs p").fill(p, c.getCalorimeterHits().size());
152+
aida.cloud2D("cluster time vs p").fill(p, ClusterUtilities.getSeedHitTime(c));
153+
}
154+
double ct = ClusterUtilities.getSeedHitTime(c);
155+
156+
if (c.getEnergy() > clusterCut
157+
&& ClusterUtilities.findSeedHit(c).getCorrectedEnergy() > seedCut
158+
&& c.getCalorimeterHits().size() >= minHits
159+
&& ct > ctMin
160+
&& ct < ctMax) {
161+
double chiSquared = t.getChi2();
162+
int ndf = t.getNDF();
163+
double chi2Ndf = t.getChi2() / t.getNDF();
164+
double chisqProb = ChisqProb.gammp(ndf, chiSquared);
165+
int nHits = t.getTrackerHits().size();
166+
double dEdx = t.getdEdx();
167+
//rotate into physiscs frame of reference
168+
Hep3Vector rprot = VecOp.mult(beamAxisRotation, rp.getMomentum());
169+
double theta = Math.acos(rprot.z() / rprot.magnitude());
170+
171+
// debug diagnostics to set cuts
172+
if (debug) {
173+
aida.cloud1D("Track chisq per df").fill(chiSquared / ndf);
174+
aida.cloud1D("Track chisq prob").fill(chisqProb);
175+
aida.cloud1D("Track nHits").fill(t.getTrackerHits().size());
176+
aida.cloud1D("Track momentum").fill(p);
177+
aida.cloud1D("Track deDx").fill(t.getdEdx());
178+
aida.cloud1D("Track theta").fill(theta);
179+
aida.cloud2D("Track theta vs p").fill(theta, p);
180+
aida.cloud1D("rp x0").fill(TrackUtils.getX0(t));
181+
aida.cloud1D("rp y0").fill(TrackUtils.getY0(t));
182+
aida.cloud1D("rp z0").fill(TrackUtils.getZ0(t));
183+
}
184+
185+
double trackDataTime = TrackData.getTrackTime(TrackData.getTrackData(event, t));
186+
aida.cloud1D("track data time").fill(trackDataTime);
187+
if (isTopTrack(t)) {
188+
if(_dumpRunAndEventNumber) System.out.println(event.getRunNumber() + " " + event.getEventNumber() + " top");
189+
trkChisqNdfTop.fill(chi2Ndf);
190+
trkChisqProbTop.fill(chisqProb);
191+
trkNhitsTop.fill(nHits);
192+
trkMomentumTop.fill(p);
193+
trkthetaTop.fill(theta);
194+
trkX0Top.fill(TrackUtils.getX0(t));
195+
trkY0Top.fill(TrackUtils.getY0(t));
196+
trkZ0Top.fill(TrackUtils.getZ0(t));
197+
198+
if (nHits == 5) {
199+
trkMomentumTop5.fill(p);
200+
trkdEdXTop5.fill(dEdx);
201+
} else if (nHits == 6) {
202+
trkMomentumTop6.fill(p);
203+
trkdEdXTop6.fill(dEdx);
204+
if (_dumpRunAndEventNumber) {
205+
System.out.println(event.getRunNumber() + " " + event.getEventNumber());
206+
}
207+
}
208+
} else {
209+
if(_dumpRunAndEventNumber) System.out.println(event.getRunNumber() + " " + event.getEventNumber() + " bottom");
210+
trkChisqNdfBottom.fill(chi2Ndf);
211+
trkChisqProbBottom.fill(chisqProb);
212+
trkNhitsBottom.fill(nHits);
213+
trkMomentumBottom.fill(p);
214+
trkthetaBottom.fill(theta);
215+
trkX0Bottom.fill(TrackUtils.getX0(t));
216+
trkY0Bottom.fill(TrackUtils.getY0(t));
217+
trkZ0Bottom.fill(TrackUtils.getZ0(t));
218+
219+
if (nHits == 5) {
220+
trkMomentumBottom5.fill(p);
221+
trkdEdXBottom5.fill(dEdx);
222+
} else if (nHits == 6) {
223+
trkMomentumBottom6.fill(p);
224+
trkdEdXBottom6.fill(dEdx);
225+
if (_dumpRunAndEventNumber) {
226+
System.out.println(event.getRunNumber() + " " + event.getEventNumber());
227+
}
228+
}
229+
}
230+
}// end of cluster cuts
231+
}
232+
}
233+
234+
public void setDumpRunAndEventNumber(boolean b) {
235+
_dumpRunAndEventNumber = b;
236+
}
237+
238+
public void setAidaFileName(String s)
239+
{
240+
_aidaFileName = s;
241+
}
242+
243+
@Override
244+
protected void endOfData() {
245+
try {
246+
AIDA.defaultInstance().saveAs(_aidaFileName);
247+
//AIDA.defaultInstance().saveAs(testOutputDir.getPath() + File.separator + this.getClass().getSimpleName() + ".root");
248+
} catch (IOException e) {
249+
throw new RuntimeException(e);
250+
}
251+
}
252+
253+
254+
private boolean isTopTrack(Track t) {
255+
List<TrackerHit> hits = t.getTrackerHits();
256+
int n[] = {0, 0};
257+
int nHits = hits.size();
258+
for (TrackerHit h : hits) {
259+
HpsSiSensor sensor = ((HpsSiSensor) ((RawTrackerHit) h.getRawHits().get(0)).getDetectorElement());
260+
if (sensor.isTopLayer()) {
261+
n[0] += 1;
262+
} else {
263+
n[1] += 1;
264+
}
265+
}
266+
if (n[0] == nHits && n[1] == 0) {
267+
return true;
268+
}
269+
if (n[1] == nHits && n[0] == 0) {
270+
return false;
271+
}
272+
throw new RuntimeException("mixed top and bottom hits on same track");
273+
274+
}
275+
276+
private void setupSensors(EventHeader event) {
277+
List<RawTrackerHit> rawTrackerHits = event.get(RawTrackerHit.class, "SVTRawTrackerHits");
278+
EventHeader.LCMetaData meta = event.getMetaData(rawTrackerHits);
279+
// Get the ID dictionary and field information.
280+
IIdentifierDictionary dict = meta.getIDDecoder().getSubdetector().getDetectorElement().getIdentifierHelper().getIdentifierDictionary();
281+
int fieldIdx = dict.getFieldIndex("side");
282+
int sideIdx = dict.getFieldIndex("strip");
283+
for (RawTrackerHit hit : rawTrackerHits) {
284+
// The "side" and "strip" fields needs to be stripped from the ID for sensor lookup.
285+
IExpandedIdentifier expId = dict.unpack(hit.getIdentifier());
286+
expId.setValue(fieldIdx, 0);
287+
expId.setValue(sideIdx, 0);
288+
IIdentifier strippedId = dict.pack(expId);
289+
// Find the sensor DetectorElement.
290+
List<IDetectorElement> des = DetectorElementStore.getInstance().find(strippedId);
291+
if (des == null || des.size() == 0) {
292+
throw new RuntimeException("Failed to find any DetectorElements with stripped ID <0x" + Long.toHexString(strippedId.getValue()) + ">.");
293+
} else if (des.size() == 1) {
294+
hit.setDetectorElement((SiSensor) des.get(0));
295+
} else {
296+
// Use first sensor found, which should work unless there are sensors with duplicate IDs.
297+
for (IDetectorElement de : des) {
298+
if (de instanceof SiSensor) {
299+
hit.setDetectorElement((SiSensor) de);
300+
break;
301+
}
302+
}
303+
}
304+
// No sensor was found.
305+
if (hit.getDetectorElement() == null) {
306+
throw new RuntimeException("No sensor was found for hit with stripped ID <0x" + Long.toHexString(strippedId.getValue()) + ">.");
307+
}
308+
}
309+
}
310+
311+
}

0 commit comments

Comments
 (0)