| Commit in lcsim-contrib/src/main/java/org/lcsim/contrib/mgraham/sATLASDigi on MAIN | |||
| FindableTrack.java | +223 | added 1.1 | |
| LOIEffFakeMG.java | +645 | added 1.1 | |
| MyDigiHitMaker.java | +64 | added 1.1 | |
| PixelCheater.java | +59 | added 1.1 | |
| TrackAnalysis.java | +86 | added 1.1 | |
| TrackerHitCheater.java | +159 | added 1.1 | |
| TrackerHitDriver_sATLAS.java | +70 | added 1.1 | |
| TrackerHitSmearing.java | +114 | added 1.1 | |
| sATLASDigiDriver.java | +21 | added 1.1 | |
| +1441 | |||
added mgraham contrib area and sATLASDigi project
diff -N FindableTrack.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ FindableTrack.java 20 Apr 2009 21:19:43 -0000 1.1 @@ -0,0 +1,223 @@
+/*
+ * FindableTrack.java
+ *
+ * Created on October 24, 2008, 9:50 PM
+ *
+ */
+package org.lcsim.contrib.mgraham.sATLASDigi;
+
+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.recon.tracking.seedtracker.SeedLayer;
+import org.lcsim.recon.tracking.seedtracker.SeedLayer.SeedType;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+
+/**
+ *
+ * @author Richard Partridge
+ * @version 1.0
+ */
+public class FindableTrack {
+
+ public enum Ignore {
+
+ NoPTCut, NoDCACut, NoZ0Cut, NoSeedCheck, NoConfirmCheck, NoMinHitCut
+ };
+ private double _bfield;
+ private RelationalTable _hittomc;
+ private HitIdentifier _ID;
+
+ /** Creates a new instance of FindableTrack */
+ public FindableTrack(EventHeader event) {
+
+ // Get the magnetic field
+ Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
+ _bfield = event.getDetector().getFieldMap().getField(IP).z();
+
+ // Instantiate the hit identifier class
+ _ID = new HitIdentifier();
+
+ // Create a relational table that maps SimTrackerHits to MCParticles
+ _hittomc = 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) {
+ _hittomc.add(simhit, simhit.getMCParticle());
+ }
+ }
+ }
+
+ 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 = _ID.Identifier(getDetectorElement(simhit));
+ if (!idset.contains(identifier)) idset.add(identifier);
+ }
+
+ return idset.size();
+ }
+
+ 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);
+
+ // 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;
+ 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;
+ }
+}
\ No newline at end of file
diff -N LOIEffFakeMG.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ LOIEffFakeMG.java 20 Apr 2009 21:19:43 -0000 1.1 @@ -0,0 +1,645 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.mgraham.sATLASDigi;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogramFactory;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+
+import hep.physics.vec.VecOp;
+import java.io.File;
+import java.util.ArrayList;
+import java.util.List;
+
+import java.util.Set;
+//import org.lcsim.contrib.Partridge.TrackingTest.FindableTrack.Ignore;
+//import org.lcsim.contrib.Partridge.TrackingTest.TrackAnalysis;
+import org.lcsim.contrib.mgraham.sATLASDigi.FindableTrack.Ignore;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.Track;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.fit.helicaltrack.HelicalTrackCross;
+import org.lcsim.fit.helicaltrack.HelicalTrackFit;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
+import org.lcsim.fit.helicaltrack.HelixUtils;
+import org.lcsim.fit.helicaltrack.MultipleScatter;
+import org.lcsim.fit.helicaltrack.TrackDirection;
+import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+import org.lcsim.recon.tracking.seedtracker.SeedTrack;
+import org.lcsim.recon.tracking.seedtracker.StrategyXMLUtils;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import sun.java2d.loops.FillPath;
+
+/**
+ *
+ * @author partridge
+ */
+public class LOIEffFakeMG extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ private IHistogram1D pTeff1;
+ private IHistogram1D pTeff2;
+ private IHistogram1D thetaeff;
+ private IHistogram1D ctheff;
+ private IHistogram1D etaeff;
+ private IHistogram1D etafake;
+ private IHistogram1D d0eff1;
+ private IHistogram1D d0eff2;
+ private IHistogram1D z0eff1;
+ private IHistogram1D z0eff2;
+ private IHistogram1D fakes;
+ private IHistogram1D nfakes;
+ int ntrk = 0;
+ int nevt = 0;
+
+ public LOIEffFakeMG() {
+
+ // Define the efficiency histograms
+ IHistogramFactory hf = aida.histogramFactory();
+ pTeff1 = hf.createHistogram1D("Efficiency vs pT", "", 100, 0., 5., "type=efficiency");
+ pTeff2 = hf.createHistogram1D("Efficiency vs pT full", "", 100, 0., 50., "type=efficiency");
+ thetaeff = hf.createHistogram1D("Efficiency vs theta", "", 72, 0., 180., "type=efficiency");
+ ctheff = hf.createHistogram1D("Efficiency vs cos(theta)", "", 50, -1., 1., "type=efficiency");
+ etaeff = hf.createHistogram1D("Efficiency vs eta", "", 50, -2., 2., "type=efficiency");
+ etafake = hf.createHistogram1D("Fake rate vs eta", "", 50, -2., 2., "type=efficiency");
+ d0eff1 = hf.createHistogram1D("Efficiency vs d0", "", 50, -5., 5., "type=efficiency");
+ d0eff2 = hf.createHistogram1D("Efficiency vs d0 full", "", 24, -12., 12., "type=efficiency");
+ z0eff1 = hf.createHistogram1D("Efficiency vs z0", "", 50, -5., 5., "type=efficiency");
+ z0eff2 = hf.createHistogram1D("Efficiency vs z0 full", "", 24, -12., 12., "type=efficiency");
+ fakes = hf.createHistogram1D("Number of mis-matched hits (unnormalized)", "", 10, 0., 10.);
+ nfakes = hf.createHistogram1D("Number of mis-matched hits (normalized)", "", 10, 0., 10.);
+
+ }
+
+ @Override
+ public void process(EventHeader event) {
+
+ // Increment the event counter
+ nevt++;
+
+ // Get the magnetic field
+ Hep3Vector IP = new BasicHep3Vector(0., 0., 0.);
+ double bfield = event.getDetector().getFieldMap().getField(IP).z();
+// dump SThit information
+ String[] input_hit_collections = {"VtxBarrHits", "VtxEndcapHits", "SCTShortBarrHits", "SCTLongBarrHits", "SCTShortEndcapHits", "SCTLongEndcapHits"};
+ for (String input : input_hit_collections) {
+ List<SimTrackerHit> sthits = event.getSimTrackerHits(input);
+ int[] nhits = {0, 0, 0, 0, 0, 0};
+ for (SimTrackerHit st : sthits) {
+ String detector = st.getDetectorElement().getName();
+ int layer = st.getLayerNumber();
+ double[] hp = st.getPoint();
+ Hep3Vector hitPos = new BasicHep3Vector(hp[0], hp[1], hp[2]);
+ double r = Math.sqrt(hp[0] * hp[0] + hp[1] * hp[1]);
+ double theta = Math.atan2(r, hp[2]);
+// double eta = -Math.log(r / (2 * Math.abs(hp[2])));
+ double eta = -Math.log(Math.tan(theta / 2));
+ double phi = Math.atan2(hp[1], hp[0]);
+// System.out.println("r= " + r + " theta = "+theta+" eta = " + eta+ " phi=" + phi);
+ nhits[layer]++;
+ aida.cloud1D(input + " layer " + layer + " STHit eta").fill(eta);
+ aida.cloud1D(input + " layer " + layer + " STHit phi").fill(phi);
+ aida.cloud2D(input + " layer " + layer + " STHit phi vs eta").fill(eta, phi);
+ aida.histogram2D(input + " layer " + layer + " STHit phi vs eta occupancy", 100, -2.5, 2.5, 100, -3.2, 3.2).fill(eta, phi);
+ }
+ int i = 0;
+ while (i < 6) {
+ if (nhits[i] > 0) {
+ aida.cloud1D(input + "layer " + i + " number of ST hits").fill(nhits[i]);
+ }
+ i++;
+ }
+ }
+ List<HelicalTrackHit> hthits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
+ for (HelicalTrackHit HelTrHit : hthits) {
+ }
+
+ // Get the list of strategies being used
+// String sfile = "autogen_ttbar_sid02_vs.xml";
+// List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(
+// StrategyXMLUtils.getDefaultStrategiesPrefix() + sfile);
+
+ String sfile = "/nfs/sulky21/g.ec.u12/users/mgraham/AtlasUpgrade/lcsim/resources/org/lcsim/recon/tracking/seedtracker/strategies/sATLASBarrel-SM08.xml";
+ List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromFile(new File(sfile));
+
+ // Find the minimum pT among the strategies
+ double ptCut = 9999.;
+ for (SeedStrategy s : slist) {
+ if (s.getMinPT() < ptCut) {
+ ptCut = s.getMinPT();
+ }
+ }
+
+ // Create a relational table that maps TrackerHits to MCParticles
+ RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> mcrelations = event.get(LCRelation.class, "HelicalTrackMCRelations");
+ for (LCRelation relation : mcrelations) {
+ hittomc.add(relation.getFrom(), relation.getTo());
+ }
+
+ // Instantiate the class that determines if a track is "findable"
+ FindableTrack findable = new FindableTrack(event);
+
+ // Create a map between tracks and the associated MCParticle
+ List<Track> tracklist = event.getTracks();
+ RelationalTable trktomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+
+ // Analyze the tracks in the event
+ for (Track track : tracklist) {
+
+ // Calculate the track pT and cos(theta)
+ double px = track.getPX();
+ double py = track.getPY();
+ double pz = track.getPZ();
+ double pt = Math.sqrt(px * px + py * py);
+ double cth = pz / Math.sqrt(pt * pt + pz * pz);
+ double th = Math.atan2(pt, pz);
+ double eta = -Math.log(Math.tan(th / 2));
+ double d0 = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double z0 = track.getTrackParameter(HelicalTrackFit.z0Index);
+
+ SeedTrack st = (SeedTrack) track;
+
+ SeedCandidate seed = st.getSeedCandidate();
+
+ HelicalTrackFit helixTrack = seed.getHelix();
+ double[] chisq = helixTrack.chisq();
+ double nhchisq = helixTrack.nhchisq();
+ aida.cloud1D("Track Chi2-Circle Fit").fill(chisq[0]);
+ aida.cloud1D("Track Chi2-RZ Fit").fill(chisq[1]);
+ aida.cloud1D("NH Track Chi2").fill(nhchisq);
+ if (nhchisq != 0) {
+ aida.cloud1D("NH!=0 Track Chi2-Circle Fit").fill(chisq[0]);
+ aida.cloud1D("NH!=0 Track Chi2-RZ Fit").fill(chisq[1]);
+
+ }
+ List<HelicalTrackHit> hitlist = seed.getHits();
+ for (HelicalTrackHit hit : hitlist) {
+ int nhits = hit.getRawHits().size();
+ aida.cloud1D(hit.Detector() + " nHits").fill(nhits);
+ Hep3Vector HTHPos = hit.getCorrectedPosition();
+ double rHit = Math.sqrt(HTHPos.x() * HTHPos.x() + HTHPos.y() * HTHPos.y());
+ double zHit = HTHPos.z();
+ double etaHit = -Math.log(Math.tan(Math.atan2(rHit, zHit) / 2));
+// System.out.println("Looping over hit in " + hit.Detector());
+ double hitchisq = hit.chisq();
+ double s = helixTrack.PathMap().get(hit);
+ Hep3Vector posonhelix = HelixUtils.PointOnHelix(helixTrack, s);
+
+ if (hit instanceof HelicalTrackCross) {
+ HelicalTrackCross cross = (HelicalTrackCross) hit;
+ TrackDirection trkdir = HelixUtils.CalculateTrackDirection(helixTrack, s);
+ cross.setTrackDirection(trkdir, helixTrack.covariance());
+ List<HelicalTrackStrip> clusterlist = cross.getStrips();
+ double du_stereo = 0;
+ double du_axial = 0;
+ for (HelicalTrackStrip cluster : clusterlist) {
+ int nstrips = cluster.rawhits().size();
+ aida.cloud1D(hit.Detector() + " nStrips-per-layer").fill(nstrips);
+ Hep3Vector corigin = cluster.origin();
+ Hep3Vector u = cluster.u();
+ List<RawTrackerHit> rawhits = cluster.rawhits();
+ double umc = -999999;
+ double stenergy = -999999;
+ String stripdir = "axial";
+ double umeas = cluster.umeas();
+ double charge = cluster.dEdx();
+ double layer = cluster.layer();
+ for (RawTrackerHit rhit : rawhits) {
+ String deName = rhit.getDetectorElement().getName();
+ if (deName.contains("sensor1")) {
+ stripdir = "stereo";
+ }
+ // System.out.println("Layer number " + rhit.getLayerNumber() + " " + deName);
+ List<SimTrackerHit> sthits = rhit.getSimTrackerHits();
+ int nsthits = sthits.size();
+ aida.cloud1D(hit.Detector() + " associated ST hits").fill(nsthits);
+ aida.cloud1D(hit.Detector() + " layer" + stripdir + " associated ST hits").fill(nsthits);
+ if (nsthits == 1) {
+ double[] sthitD = sthits.get(0).getPoint();
+ BasicHep3Vector sthit = new BasicHep3Vector(sthitD);
+ stenergy = sthits.get(0).getdEdx();
+ Hep3Vector vdiff = VecOp.sub(sthit, corigin);
+ umc = VecOp.dot(vdiff, u);
+ }
+ }
+
+
+ aida.cloud2D(hit.Detector() + "clusterSize vs eta").fill(etaHit, nstrips);
+ // System.out.println("filling...");
+ if (umc != -999999) {
+ aida.cloud2D(hit.Detector() + "cluster vs STHit dedx").fill(stenergy, charge);
+ aida.cloud2D(hit.Detector() + "cluster dedx vs delte(u)").fill(umeas - umc, charge);
+ if (stripdir.contains("stereo")) {
+ du_stereo = umeas - umc;
+ }
+ if (stripdir.contains("axial")) {
+ du_axial = umeas - umc;
+ }
+ aida.cloud1D(hit.Detector() + "layer=" + stripdir + " delta(u)").fill(umeas - umc);
+ aida.cloud1D(hit.Detector() + " delta(u)").fill(umeas - umc);
+ if (nstrips == 1) {
+ aida.cloud1D(hit.Detector() + "layer=" + stripdir + " delta(u)--1 strip").fill(umeas - umc);
+ aida.cloud1D(hit.Detector() + " delta(u)--1 strip").fill(umeas - umc);
+ }
+ if (nstrips == 2) {
+ aida.cloud1D(hit.Detector() + "layer=" + stripdir + " delta(u)--2 strip").fill(umeas - umc);
+ aida.cloud1D(hit.Detector() + " delta(u)--2 strip").fill(umeas - umc);
+ }
+ if (nstrips == 3) {
+ aida.cloud1D(hit.Detector() + "layer=" + stripdir + " delta(u)--3 strip").fill(umeas - umc);
+ aida.cloud1D(hit.Detector() + " delta(u)--3 strip").fill(umeas - umc);
+ }
+ }
+
+ }
+ aida.cloud2D(hit.Detector() + " delta(u) stereo v axial").fill(du_stereo, du_axial);
+ }
+ MultipleScatter ms = seed.getMSMap().get(hit);
+ double msphi = ms.drphi();
+ double msz = ms.dz();
+ Hep3Vector hitpos = hit.getCorrectedPosition();
+ SymmetricMatrix cov = hit.getCorrectedCovMatrix();
+// double rHit = getr(hitpos.x(), hitpos.y());
+// double rHel = getr(posonhelix.x(), posonhelix.y());
+// double phiHet = getphi(hitpos.x(), hitpos.y());
+// double phiHel = getr(posonhelix.x(), posonhelix.y());
+ double dxdy = Math.sqrt(Math.pow(hitpos.x() - posonhelix.x(), 2) + Math.pow(hitpos.y() - posonhelix.y(), 2));
+ double dxdyErr = getdxdyErr(hitpos, posonhelix, cov);
+ double dz = posonhelix.z() - hitpos.z();
+
+ double dzErr = Math.sqrt(cov.e(2, 2));
+ aida.cloud1D(hit.Detector() + " dxdy").fill(dxdy);
+ aida.cloud1D(hit.Detector() + " dz").fill(dz);
+ aida.cloud1D(hit.Detector() + " dxdy Pull").fill(dxdy / dxdyErr);
+ aida.cloud1D(hit.Detector() + " dz Pull").fill(dz / dzErr);
+ if (Math.abs(dz) > 4) {
+ aida.cloud1D(hit.Detector() + "Bad dz--nHits").fill(nhits);
+ }
+ aida.cloud1D("NH Chi2 for Hits on Track").fill(hitchisq);
+ }
+
+ // Analyze the hits on the track
+ TrackAnalysis tkanal = new TrackAnalysis(track, hittomc);
+
+ // Calculate purity and make appropriate plots
+ int nbad = tkanal.getNBadHits();
+ int nhits = tkanal.getNHits();
+ double purity = tkanal.getPurity();
+ aida.histogram1D("Mis-matched hits for all tracks", 10, 0., 10.).fill(nbad);
+ aida.histogram1D("Mis-matched hits " + nhits + " hit tracks", 10, 0., 10.).fill(nbad);
+
+ // Generate a normalized histogram after 1000 events
+ ntrk++;
+ if (nevt <= 1000) {
+ fakes.fill(nbad);
+ }
+
+ // Make plots for fake, non-fake, and all tracks
+ if (purity < 0.5) {
+ aida.histogram1D("Hits for fake tracks", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for fake tracks", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for fake tracks", 100, -1., 1.).fill(cth);
+ aida.histogram1D("d0 for fake tracks", 50, -10., 10.).fill(d0);
+ aida.histogram1D("z0 for fake tracks", 50, -10., 10.).fill(z0);
+ aida.histogram1D("eta for fake tracks", 100, -2., 2.).fill(eta);
+ etafake.fill(eta, 1.0);
+ } else {
+ aida.histogram1D("Hits for non-fake tracks", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for non-fake tracks", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for non-fake tracks", 100, -1., 1.).fill(cth);
+ aida.histogram1D("d0 for non-fake tracks", 50, -10., 10.).fill(d0);
+ aida.histogram1D("z0 for non-fake tracks", 50, -10., 10.).fill(z0);
+ aida.histogram1D("eta for non-fake tracks", 100, -2., 2.).fill(eta);
+ etafake.fill(eta, 0.0);
+ }
+ aida.histogram1D("Hits for all tracks", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for all tracks", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for all tracks", 100, -1., 1.).fill(cth);
+ aida.histogram1D("d0 for all tracks", 50, -10., 10.).fill(d0);
+ aida.histogram1D("z0 for all tracks", 50, -10., 10.).fill(z0);
+ aida.histogram1D("eta for all tracks", 100, -2., 2.).fill(eta);
+ aida.cloud2D("Hits vs eta for all tracks").fill(eta, nhits);
+ // Now analyze MC Particles on this track
+ MCParticle mcp = tkanal.getMCParticle();
+ if (mcp != null) {
+
+ // Create a map between the tracks found and the assigned MC particle
+ trktomc.add(track, tkanal.getMCParticle());
+
+ // Calculate the MC momentum and polar angle
+ Hep3Vector pmc = mcp.getMomentum();
+ double pxmc = pmc.x();
+ double pymc = pmc.y();
+ double ptmc = Math.sqrt(pxmc * pxmc + pymc * pymc);
+ double pxtk = track.getPX();
+ double pytk = track.getPY();
+ double pttk = Math.sqrt(pxtk * pxtk + pytk * pytk);
+
+ // Calculate the helix parameters for this MC particle and pulls in pT, d0
+ HelixParamCalculator helix = new HelixParamCalculator(mcp, bfield);
+ double d0tk = track.getTrackParameter(HelicalTrackFit.dcaIndex);
+ double d0mc = helix.getDCA();
+ double d0err = Math.sqrt(track.getErrorMatrix().diagonal(HelicalTrackFit.dcaIndex));
+ double curv = track.getTrackParameter(HelicalTrackFit.curvatureIndex);
+ double curverr = Math.sqrt(track.getErrorMatrix().diagonal(HelicalTrackFit.curvatureIndex));
+ double pterr = pttk * curverr / curv;
+ double d0pull = (d0tk - d0mc) / d0err;
+ double ptpull = (pttk - ptmc) / pterr;
+ double ptresid = (pttk - ptmc);
+ double d0resid = (d0tk - d0mc);
+ // Plot the pt and d0 pulls for various purity intervals
+ if (nbad == 0) {
+ aida.histogram2D("pT MC vs pT Reco for 0 Bad Hits",
+ 100, 0., 5., 100, 0., 5.).fill(ptmc, pttk);
+ aida.histogram2D("d0 MC vs d0 Reco for 0 Bad Hits",
+ 100, -0.2, 0.2, 100, -0.2, 0.2).fill(d0mc, d0tk);
+ aida.histogram1D("pT Pull for 0 Bad Hits", 100, -10., 10.).fill(ptpull);
+ aida.histogram1D("d0 pull for 0 Bad Hits", 100, -10., 10.).fill(d0pull);
+ aida.cloud1D("pT Pull for 0 Bad Hits").fill(ptresid);
+ aida.cloud1D("d0 pull for 0 Bad Hits").fill(d0resid);
+
+ } else if (purity > 0.5) {
+ aida.histogram2D("pT MC vs pT Reco for 0.5 < purity < 1",
+ 100, 0., 5., 100, 0., 5.).fill(ptmc, pttk);
+ aida.histogram2D("d0 MC vs d0 Reco for 0.5 < purity < 1",
+ 100, -0.2, 0.2, 100, -0.2, 0.2).fill(d0mc, d0tk);
+ aida.histogram1D("pT Pull for 0.5 < purity < 1", 100, -10., 10.).fill(ptpull);
+ aida.histogram1D("d0 pull for 0.5 < purity < 1", 100, -10., 10.).fill(d0pull);
+ aida.cloud1D("pT Pull for 0.5 < purity < 1").fill(ptresid);
+ aida.cloud1D("d0 pull for 0.5 < purity < 1").fill(d0resid);
+ } else if (purity < 0.5) {
+ aida.histogram2D("pT MC vs pT Reco for purity <= 0.5",
+ 100, 0., 5., 100, 0., 5.).fill(ptmc, pttk);
+ aida.histogram2D("d0 MC vs d0 Reco for purity <= 0.5",
+ 100, -0.2, 0.2, 100, -0.2, 0.2).fill(d0mc, d0tk);
+ aida.histogram1D("pT Pull for purity <= 0.5", 100, -10., 10.).fill(ptpull);
+ aida.histogram1D("d0 pull for purity <= 0.5", 100, -10., 10.).fill(d0pull);
+ aida.cloud1D("pT Pull for purity <= 0.5").fill(ptresid);
+ aida.cloud1D("d0 pull for purity <= 0.5").fill(d0resid);
+ }
+ }
+ }
+
+ // Make the normalized fake plot after the specified number of events
+ if (nevt == 1000) {
+ double wgt = 1. / ntrk;
+ for (int i = 0; i < 10; i++) {
+ System.out.println(" Entries: " + fakes.binEntries(i) + " for mismatches: " + i);
+ for (int j = 0; j < fakes.binHeight(i); j++) {
+ nfakes.fill(i, wgt);
+ }
+ }
+ System.out.println("Normalization: " + nfakes.sumAllBinHeights() + " after ntrk = " + ntrk);
+ }
+
+ // Now loop over all MC Particles
+ List<MCParticle> mclist = event.getMCParticles();
+ for (MCParticle mcp : mclist) {
+
+ // Calculate the pT and polar angle of the MC particle
+ double px = mcp.getPX();
+ double py = mcp.getPY();
+ 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 theta = 180. * Math.acos(cth) / Math.PI;
+ double eta = -Math.log(Math.tan(Math.atan2(pt, pz) / 2));
+ // Find the number of layers hit by this mc particle
+// System.out.println("MC pt=" + pt);
+ int nhits = findable.LayersHit(mcp);
+
+ // Calculate the helix parameters for this MC particle
+ HelixParamCalculator helix = new HelixParamCalculator(mcp, bfield);
+ double d0 = helix.getDCA();
+ double z0 = helix.getZ0();
+
+ // Check cases where we have multiple tracks associated with this MC particle
+ Set<Track> trklist = trktomc.allTo(mcp);
+ int ntrk = trklist.size();
+// if (ntrk > 1) {
+ // Count tracks where the assigned MC particle has more than 1 hit
+// int nmulthits = 0;
+// for (Track trk : trklist) {
+// TrackAnalysis tkanal = new TrackAnalysis(trk, hittomc);
+// if (tkanal.getNBadHits() < tkanal.getNHits() - 1)
+// nmulthits++;
+// }
+ // Flag any anomalous cases that we find
+// if (nmulthits > 1) {
+// System.out.println("2 tracks associated with a single MC Particle");
+// for (Track trk : trklist) System.out.println(trk.toString());
+// }
+// }
+
+ // Make pT efficiency plot
+ if (findable.isFindable(mcp, slist, Ignore.NoPTCut)) {
+ double wgt = 0.;
+ if (ntrk > 0) {
+ wgt = 1.;
+ }
+ pTeff1.fill(pt, wgt);
+ pTeff2.fill(pt, wgt);
+ }
+
+ // Make angular efficiency plot
+ if (findable.isFindable(mcp, slist)) {
+ double wgt = 0.;
+ if (ntrk > 0) {
+ wgt = 1.;
+ }
+ thetaeff.fill(theta, wgt);
+ ctheff.fill(cth, wgt);
+ etaeff.fill(eta, wgt);
+ }
+
+ // Make d0 efficiency plot
+ if (findable.isFindable(mcp, slist, Ignore.NoDCACut)) {
+ double wgt = 0.;
+ if (ntrk > 0) {
+ wgt = 1.;
+ }
+ d0eff1.fill(d0, wgt);
+ d0eff2.fill(d0, wgt);
+ }
+
+ // Make z0 efficiency plot
+ if (findable.isFindable(mcp, slist, Ignore.NoZ0Cut)) {
+ double wgt = 0.;
+ if (ntrk > 0) {
+ wgt = 1.;
+ }
+ z0eff1.fill(z0, wgt);
+ z0eff2.fill(z0, wgt);
+ }
+
+ // Select charged MC particles
+ if (mcp.getCharge() == 0) {
+ continue;
+ }
+
+ // Select mcp that fail the final state requirement
+ if (mcp.getGeneratorStatus() != mcp.FINAL_STATE) {
+ aida.histogram1D("Hits for non-final state particles", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for non-final state particles", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for non-final state particles", 100, -1., 1.).fill(cth);
+ aida.histogram1D("eta for non-final state particles", 100, -2., 2.).fill(eta);
+ aida.histogram1D("d0 for non-final state particles", 100, -100., 100.).fill(d0);
+ aida.histogram1D("z0 for non-final state particles", 100, -100., 100.).fill(z0);
+ aida.cloud2D("Hits vs eta for non-final state particles").fill(eta, nhits);
+ continue;
+ }
+
+ // Make plots for the base sample
+ aida.histogram1D("Hits for base MC selection", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for base MC selection", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for base MC selection", 100, -1., 1.).fill(cth);
+ aida.histogram1D("eta for base MC selection", 100, -2., 2.).fill(eta);
+ aida.histogram1D("d0 for base MC selection", 100, -100., 100.).fill(d0);
+ aida.histogram1D("z0 for base MC selection", 100, -100., 100.).fill(z0);
+ aida.cloud2D("Hits vs eta for base MC selection").fill(eta, nhits);
+
+ // Make plots for findable tracks
+ if (findable.isFindable(mcp, slist)) {
+ aida.histogram1D("Hits for findable tracks", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for findable tracks", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for findable tracks", 100, -1., 1.).fill(cth);
+ aida.histogram1D("eta for findable tracks", 100, -2., 2.).fill(eta);
+ aida.histogram1D("d0 for findable tracks", 100, -100., 100.).fill(d0);
+ aida.histogram1D("z0 for findable tracks", 100, -100., 100.).fill(z0);
+ aida.cloud2D("Hits vs eta for findable tracks").fill(eta, nhits);
+ continue;
+ }
+
+ // Create the running list of conditions to ignore
+ List<Ignore> ignores = new ArrayList<Ignore>();
+
+ // select mc particles that fail on the z0 cut
+ ignores.add(Ignore.NoZ0Cut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for z0 check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for z0 check failures", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for z0 check failures", 100, -1., 1.).fill(cth);
+ aida.histogram1D("eta for z0 check failures", 100, -2., 2.).fill(eta);
+ aida.histogram1D("d0 for z0 check failures", 100, -100., 100.).fill(d0);
+ aida.histogram1D("z0 for z0 check failures", 100, -100., 100.).fill(z0);
+ continue;
+ }
+
+ // Select mc particles that fail on the d0 cut
+ ignores.add(Ignore.NoDCACut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for d0 check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for d0 check failures", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for d0 check failures", 100, -1., 1.).fill(cth);
+ aida.histogram1D("eta for d0 check failures", 100, -2., 2.).fill(eta);
+ aida.histogram1D("d0 for d0 check failures", 100, -100., 100.).fill(d0);
+ aida.histogram1D("z0 for d0 check failures", 100, -100., 100.).fill(z0);
+ continue;
+ }
+
+ // select mc particles that fail the confirm check
+ ignores.add(Ignore.NoConfirmCheck);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for confirm check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for confir check failures", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for confirm check failures", 100, -1., 1.).fill(cth);
+ aida.histogram1D("eta for confirm check failures", 100, -2., 2.).fill(eta);
+ aida.histogram1D("d0 for seed confirm failures", 100, -100., 100.).fill(d0);
+ aida.histogram1D("z0 for seed confirm failures", 100, -100., 100.).fill(z0);
+ continue;
+ }
+
+ // select mc particles that fail on the seed check
+ ignores.add(Ignore.NoSeedCheck);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for seed check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for seed check failures", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for seed check failures", 100, -1., 1.).fill(cth);
+ aida.histogram1D("eta for seed check failures", 100, -1., 1.).fill(eta);
+ aida.histogram1D("d0 for seed check failures", 100, -100., 100.).fill(d0);
+ aida.histogram1D("z0 for seed check failures", 100, -100., 100.).fill(z0);
+ continue;
+ }
+
+ // Select mc particles that fail the number of hit cut
+ ignores.add(Ignore.NoMinHitCut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for nhit check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for nhit check failures", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for nhit check failures", 100, -1., 1.).fill(cth);
+ aida.histogram1D("eta for nhit check failures", 100, -2., 2.).fill(eta);
+ aida.histogram1D("d0 for nhit check failures", 100, -100., 100.).fill(d0);
+ aida.histogram1D("z0 for nhit check failures", 100, -100., 100.).fill(z0);
+ continue;
+ }
+
+ // Select mc particles that fail on the pT cut
+ ignores.add(Ignore.NoPTCut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.histogram1D("Hits for pT check failures", 20, 0., 20.).fill(nhits);
+ aida.histogram1D("pT for pT check failures", 100, 0., 10.).fill(pt);
+ aida.histogram1D("cos(theta) for pT check failures", 100, -1., 1.).fill(cth);
+ aida.histogram1D("eta for pT check failures", 100, -2., 2.).fill(eta);
+ aida.histogram1D("d0 for pT check failures", 100, -100., 100.).fill(d0);
+ aida.histogram1D("z0 for pT check failures", 100, -100., 100.).fill(z0);
+ } else {
+ System.out.println("MC Particle is not findable with all ignores set!!");
+ }
+ }
+ return;
+ }
+
+ private double getr(double x, double y) {
+ return Math.sqrt(x * x + y * y);
+ }
+
+ protected double drcalc(Hep3Vector pos, SymmetricMatrix cov) {
+ double x = pos.x();
+ double y = pos.y();
+ double r2 = x * x + y * y;
+ return Math.sqrt((x * x * cov.e(0, 0) + y * y * cov.e(1, 1) + 2. * x * y * cov.e(0, 1)) / r2);
+ }
+
+ protected double drphicalc(Hep3Vector pos, SymmetricMatrix cov) {
+ double x = pos.x();
+ double y = pos.y();
+ double r2 = x * x + y * y;
+ return Math.sqrt((y * y * cov.e(0, 0) + x * x * cov.e(1, 1) - 2. * x * y * cov.e(0, 1)) / r2);
+ }
+
+ private double getphi(double x, double y) {
+ double phi = Math.atan2(y, x);
+ if (phi < 0.) {
+ phi += 2. * Math.PI;
+ }
+ return phi;
+ }
+
+ private double getdxdy(Hep3Vector hitpos, Hep3Vector posonhelix) {
+ return Math.sqrt(Math.pow(hitpos.x() - posonhelix.x(), 2) + Math.pow(hitpos.y() - posonhelix.y(), 2));
+ }
+
+ private double getdxdyErr(Hep3Vector hitpos, Hep3Vector posonhelix, SymmetricMatrix cov) {
+ double dxdySq = Math.pow(hitpos.x() - posonhelix.x(), 2) + Math.pow(hitpos.y() - posonhelix.y(), 2);
+ double ErrSqDxDySq = 4 * (cov.e(0, 0) * Math.pow(hitpos.x() - posonhelix.x(), 2) + cov.e(1, 1) * Math.pow(hitpos.y() - posonhelix.y(), 2));
+ double error = Math.sqrt(ErrSqDxDySq / dxdySq) / 2;
+ return error;
+ }
+}
+
+
diff -N MyDigiHitMaker.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MyDigiHitMaker.java 20 Apr 2009 21:19:43 -0000 1.1 @@ -0,0 +1,64 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.mgraham.sATLASDigi;
+
+import hep.aida.IHistogram2D;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import java.io.File;
+import java.util.List;
+//import org.lcsim.recon.tracking.digitization.sistripsim.TrackerHitDriver_User;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.recon.vertexing.pixsim.PixilatedSensorManager;
+import org.lcsim.recon.vertexing.pixsim.SensorOption;
+import org.lcsim.fit.helicaltrack.HelicalTrackHitDriver;
+import org.lcsim.fit.helicaltrack.HelicalTrackHitDriver.HitType;
+import org.lcsim.recon.tracking.seedtracker.SeedStrategy;
+import org.lcsim.recon.tracking.seedtracker.SeedTracker;
+import org.lcsim.recon.tracking.seedtracker.StrategyXMLUtils;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author cozzy
+ */
+public class MyDigiHitMaker extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ public MyDigiHitMaker() {
+
+ //Initialize Nick's Pixel hit making code
+// PixilatedSensorManager psm = new PixilatedSensorManager(SensorOption.ClassicCCD, true);
+// add(psm);
+ PixelCheater pcB = new PixelCheater("VtxBarrHits", 0.0075, 0.075, true);
+// PixelCheater pcB = new PixelCheater("VtxBarrHits",0.000000001,0.000000001,true);
+ add(pcB);
+ PixelCheater pcE = new PixelCheater("VtxEndcapHits", 0.0075, 0.075, false);
+ add(pcE);
+// sInitialize Tim's TrackerHit making code
+ TrackerHitDriver_sATLAS thd = new TrackerHitDriver_sATLAS();
+ add(thd);
+
+ HelicalTrackHitDriver hitdriver = new HelicalTrackHitDriver();
+ hitdriver.addCollection(((TrackerHitDriver_sATLAS) thd).getStripHits1DName(), HitType.Digitized);
+ hitdriver.addCollection("VtxBarrHits_PixelCheater", HitType.Base);
+ hitdriver.addCollection("VtxEndcapHits_PixelCheater", HitType.Base);
+// hitdriver.addCollection("RecVtxBarrHits", HitType.Base);
+// hitdriver.addCollection("RecVtxEndcapHits", HitType.Base);
+
+ hitdriver.OutputCollection("HelicalTrackHits"); // the default, but it might change?
+ add(hitdriver);
+ String sfile = "/nfs/sulky21/g.ec.u12/users/mgraham/AtlasUpgrade/lcsim/resources/org/lcsim/recon/tracking/seedtracker/strategies/sATLASBarrel-SM08.xml"; List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromFile(new File(sfile));
+
+ SeedTracker st = new SeedTracker(slist);
+ add(st);
+ }
+
+
+
+}
diff -N PixelCheater.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ PixelCheater.java 20 Apr 2009 21:19:43 -0000 1.1 @@ -0,0 +1,59 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.mgraham.sATLASDigi;
+
+import java.util.ArrayList;
+import java.util.List;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author mgraham
+ */
+public class PixelCheater extends Driver {
+
+ TrackerHitCheater _cheat;
+ ArrayList<String> _sim_event_input = new ArrayList<String>();
+ List<TrackerHit> _input_hits = null;
+ String _detName;
+
+ public PixelCheater() {
+ }
+
+ public PixelCheater(String detName, double drphi, double dz, boolean barrel) {
+ _detName = detName;
+ addSimTrackerInputCollection(detName);
+ _cheat = new TrackerHitCheater(drphi, dz, barrel);
+ }
+
+ public void addSimTrackerInputCollection(String sim_tracker_hit_collection_name) {
+ _sim_event_input.add(sim_tracker_hit_collection_name);
+ }
+
+ public String getPixelCheaterName() {
+ return _detName + "_PixelCheater";
+ }
+
+ public void process(EventHeader event) {
+ if (_sim_event_input.size() > 0) // get any simtracker input from the event
+ {
+ for (String s : _sim_event_input) {
+ ArrayList<SimTrackerHit> simeventhits = (ArrayList) event.get(SimTrackerHit.class, s);
+ if (simeventhits != null) {
+ System.out.println("PixelCheater::process getting " + _detName + "hits");
+ _input_hits=_cheat.makeTrackerHits(simeventhits);
+// _input_hits.addAll((List<TrackerHit>) _cheat.makeTrackerHits(simeventhits));
+ System.out.println("PixelCheater::process found " + _input_hits.size() + "hits");
+ }
+ }
+ }
+
+ event.put(getPixelCheaterName(), _input_hits, TrackerHit.class, 0, toString());
+
+ }
+}
diff -N TrackAnalysis.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TrackAnalysis.java 20 Apr 2009 21:19:43 -0000 1.1 @@ -0,0 +1,86 @@
+/*
+ * TrackAnalysis.java
+ *
+ * Created on October 16, 2008, 6:09 PM
+ *
+ */
+
+package org.lcsim.contrib.mgraham.sATLASDigi;
+
+import java.util.HashMap;
+import java.util.Map;
+import java.util.Set;
+
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.Track;
+import org.lcsim.event.TrackerHit;
+
+/**
+ *
+ * @author Richard Partridge
+ */
+public class TrackAnalysis {
+ private enum HelixPar {Curvature, Phi0, DCA, Z0, Slope};
+ private MCParticle _mcp = null;
+ private int _nhits;
+ private int _nbadhits;
+ private double _purity;
+
+ /** Creates a new instance of TrackAnalysis */
+ public TrackAnalysis(Track trk, RelationalTable hittomc) {
+
+ // 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>();
+
+ // Loop over the hits on the track and make sure we have HelicalTrackHits (which contain the MC particle)
+ for (TrackerHit hit : trk.getTrackerHits()) {
+
+ // get the set of MCParticles associated with this hit and update the hit count for each MCParticle
+ Set<MCParticle> mclist = hittomc.allFrom(hit);
+ for (MCParticle mcp : mclist) {
+ Integer mchits = 0;
+ if (mcmap.containsKey(mcp)) mchits = mcmap.get(mcp);
+ mchits++;
+ mcmap.put(mcp, mchits);
+ }
+ }
+
+ // Find the MCParticle that has the most hits on the track
+
+ int nbest = 0;
+ MCParticle mcbest = null;
+ for (MCParticle mcp : mcmap.keySet()) {
+ int count = mcmap.get(mcp);
+ if (count > nbest) {
+ nbest = count;
+ mcbest = mcp;
+ }
+ }
+
+ if (nbest > 0) _mcp = mcbest;
+ _purity = (double) nbest / (double) _nhits;
+ _nbadhits = _nhits - nbest;
+
+ }
+
+ public MCParticle getMCParticle() {
+ return _mcp;
+ }
+
+ public int getNHits() {
+ return _nhits;
+ }
+
+ public int getNBadHits() {
+ return _nbadhits;
+ }
+
+ public double getPurity() {
+ return _purity;
+ }
+
+}
diff -N TrackerHitCheater.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TrackerHitCheater.java 20 Apr 2009 21:19:43 -0000 1.1 @@ -0,0 +1,159 @@
+/*
+ * TrackerHitCheater.java
+ *
+ * Created on July 11, 2007, 10:54 AM
+ *
+ * To change this template, choose Tools | Template Manager
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.mgraham.sATLASDigi;
+
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.base.BaseTrackerHitMC;
+
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.util.List;
+import java.util.ArrayList;
+
+/**
+ *
+ * @author tknelson
+ */
+public class TrackerHitCheater
+{
+
+ double _max_dist = 0.1;
+ TrackerHitSmearing _smear;
+
+ /** Creates a new instance of TrackerHitCheater */
+ public TrackerHitCheater(double drphi,double dz,boolean barrel)
+ {
+ _smear= new TrackerHitSmearing(drphi,dz,barrel);
+ }
+
+ public List<TrackerHit> makeTrackerHits(List<SimTrackerHit> simulated_hits)
+ {
+
+ List<SimTrackerHit> unclustered_hits = new ArrayList<SimTrackerHit>(simulated_hits);
+
+ // Loop over hits, find clusters and create TrackerHits
+ List<TrackerHit> trackerhits = new ArrayList<TrackerHit>();
+
+ for (SimTrackerHit seed_hit : simulated_hits)
+ {
+ if (unclustered_hits.contains(seed_hit))
+ {
+ // Create a cluster with seed hit
+ List<SimTrackerHit> cluster = new ArrayList<SimTrackerHit>();
+ cluster.add(seed_hit);
+ unclustered_hits.remove(seed_hit);
+
+ // find nearby hits, move to cluster
+ while (nearbyHits(cluster,unclustered_hits).size() != 0)
+ {
+ List<SimTrackerHit> found_hits = nearbyHits(cluster,unclustered_hits);
+ cluster.addAll(found_hits);
+ unclustered_hits.removeAll(found_hits);
+ }
+ trackerhits.add(makeTrackerHit(cluster));
+ }
+ }
+ System.out.println("found this many hits "+ trackerhits.size());
+ return trackerhits;
+
+ }
+
+
+ private List<SimTrackerHit> nearbyHits(List<SimTrackerHit> cluster, List<SimTrackerHit> unclustered_hits)
+ {
+ List<SimTrackerHit> nearby_hits = new ArrayList<SimTrackerHit>();
+
+ for (SimTrackerHit seed_hit : cluster)
+ {
+ for (SimTrackerHit hit : unclustered_hits)
+ {
+ double max_dist = Math.max(_max_dist,(seed_hit.getPathLength()+hit.getPathLength())/2.0*1.01);
+ Hep3Vector seed_point = new BasicHep3Vector(seed_hit.getPoint());
+ Hep3Vector hit_point = new BasicHep3Vector(hit.getPoint());
+
+ if (VecOp.sub(hit_point,seed_point).magnitude() <= max_dist)
+ {
+ nearby_hits.add(hit);
+ }
+ }
+ }
+ return nearby_hits;
+ }
+
+ private TrackerHit makeTrackerHit(List<SimTrackerHit> cluster)
+ {
+ double[] position = getPosition(cluster);
+ SymmetricMatrix covariance = getCovariance(cluster,position);
+ double time = getTime(cluster);
+ double energy = getEnergy(cluster);
+ int type = 0;
+
+ TrackerHit hit = new BaseTrackerHitMC(position, covariance.asPackedArray(true), time, energy, type, cluster);
+ return hit;
+ }
+
+ private double[] getPosition(List<SimTrackerHit> cluster)
+ {
+ double[] position = new double[3];
+
+ for (SimTrackerHit hit : cluster)
+ {
+ position[0] += hit.getPoint()[0] * hit.getdEdx();
+ position[1] += hit.getPoint()[1] * hit.getdEdx();
+ position[2] += hit.getPoint()[2] * hit.getdEdx();
+ }
+ for (int i = 0; i < 3; i++)
+ {
+ position[i] /= getEnergy(cluster);
+ }
+ double[] smposition = _smear.SmearedPosition(new BasicHep3Vector(position)).v();
+
+ return smposition;
+ }
+
+ private double getTime(List<SimTrackerHit> cluster)
+ {
+ double mean_time = 0.0;
+
+ for (SimTrackerHit hit : cluster)
+ {
+ mean_time += hit.getTime() * hit.getdEdx();
+ }
+ mean_time /= getEnergy(cluster);
+
+ return mean_time;
+ }
+
+ private SymmetricMatrix getCovariance(List<SimTrackerHit> cluster, double[] position)
+ {
+// double[] covariance = new double[3];
+// covariance[0] = Math.pow(0.007,2);
+// covariance[1] = Math.pow(0.007,2);
+// covariance[2] = Math.pow(0.007,2);
+ SymmetricMatrix smcov = _smear.getCovariance(new BasicHep3Vector(position));
+ return smcov;
+ }
+
+ private double getEnergy(List<SimTrackerHit> cluster)
+ {
+ double total_energy = 0.0;
+ for (SimTrackerHit hit : cluster)
+ {
+ total_energy += hit.getdEdx();
+ }
+ return total_energy;
+ }
+
+
+}
diff -N TrackerHitDriver_sATLAS.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TrackerHitDriver_sATLAS.java 20 Apr 2009 21:19:43 -0000 1.1 @@ -0,0 +1,70 @@
+/*
+ * TrackerHitDriver_User.java
+ *
+ * Created on April 4, 2008, 4:13 PM
+ *
+ * To change this template, choose Tools | Template Manager
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.mgraham.sATLASDigi;
+import org.lcsim.recon.tracking.digitization.sistripsim.TrackerHitDriver;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author tknelson
+ */
+public class TrackerHitDriver_sATLAS extends Driver
+{
+
+ TrackerHitDriver _trackerhit_driver;
+
+ /** Creates a new instance of TrackerHitDriver_sATLAS */
+ public TrackerHitDriver_sATLAS()
+ {
+ _trackerhit_driver = new TrackerHitDriver();
+
+ // _trackerhit_driver.addReadout("VtxBarrHits");
+ _trackerhit_driver.addReadout("SCTShortBarrHits");
+ _trackerhit_driver.addReadout("SCTLongBarrHits");
+ //_trackerhit_driver.addReadout("VtxEndcapHits");
+ _trackerhit_driver.addReadout("SCTShortEndcapHits");
+ _trackerhit_driver.addReadout("SCTLongEndcapHits");
+
+ //_trackerhit_driver.addElementToProcess("VtxPixelBarrel");
+ _trackerhit_driver.addElementToProcess("SCTShortBarrel");
+ _trackerhit_driver.addElementToProcess("SCTLongBarrel");
+ //_trackerhit_driver.addElementToProcess("VtxPixelEndcap");
+ _trackerhit_driver.addElementToProcess("SCTShortEndcap");
+ _trackerhit_driver.addElementToProcess("SCTLongEndcap");
+
+
+ super.add( _trackerhit_driver );
+ }
+
+ // Collection names
+ //-----------------
+ public String getRawHitsName()
+ {
+ return _trackerhit_driver.getRawHitsName();
+ }
+
+ public String getStripHits1DName()
+ {
+ return _trackerhit_driver.getStripHits1DName();
+ }
+
+// public String getPixelHitsName()
+// {
+// return _trackerhit_driver.getPixelHitsName();
+// }
+
+ public String getStripHits2DName()
+ {
+ return _trackerhit_driver.getStripHits2DName();
+ }
+
+
+}
diff -N TrackerHitSmearing.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TrackerHitSmearing.java 20 Apr 2009 21:19:43 -0000 1.1 @@ -0,0 +1,114 @@
+/*
+ * TrackerHitSmearing.java
+ *
+ * Created on July 18, 2007, 11:51 AM
+ *
+ */
+package org.lcsim.contrib.mgraham.sATLASDigi;
+
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import java.util.Random;
+
+/**
+ * This class smears hit positions and returns the covariance matrix associated with the smearing.
+ * Gaussian smearing is performed using the java.util.Random class.
+ * Currently, only smearing in the azimuthal (r*phi) and z coordinates is performed.
+ * @author Richard Partridge
+ */
+public class TrackerHitSmearing {
+
+ Random rn = new Random();
+ double drphi;
+ double dz;
+ boolean _barrel;
+
+ /**
+ * Creates a new instance of TrackerHitSmearing
+ * @param drphi RMS of smearing in the azimuthal coordinate r*phi
+ * @param dz RMS of smearing in the z coordinate
+ */
+ public TrackerHitSmearing(double drphi, double dz, boolean barrel) {
+ this.drphi = drphi;
+ this.dz = dz;
+ _barrel = barrel;
+ }
+
+ /**
+ * Creates a new instance of TrackerHitSmearing
+ * @param drphi RMS of smearing in the azimuthal coordinate r*phi
+ * @param dz RMS of smearing in the z coordinate
+ * @param seed Initial random number seed
+ */
+ public TrackerHitSmearing(double drphi, double dz, boolean barrel, long seed) {
+ this.drphi = drphi;
+ this.dz = dz;
+ rn.setSeed(seed);
+ _barrel = barrel;
+ }
+
+ /**
+ * Smear the hit position
+ * @param pos Unsmeared hit position
+ * @return Smeared hit position
+ */
+ public Hep3Vector SmearedPosition(Hep3Vector pos) {
+ double x = pos.x();
+ double y = pos.y();
+ double z = pos.z();
+ double r = Math.sqrt(x * x + y * y);
+ if (_barrel) {
+ System.out.println("Smearing Barrel Hit");
+ // Smear the phi coordinate of the hit position and re-calculate the cartesian coordinates
+ if (r > 0) {
+ double phi = Math.atan2(y, x);
+ phi = phi + (drphi / r) * rn.nextGaussian();
+ x = r * Math.cos(phi);
+ y = r * Math.sin(phi);
+ }
+ // Smear the z coordinate of the hit position
+ z = z + dz * rn.nextGaussian();
+ } else {
+ System.out.println("Smearing Endcap Hit");
+ //smear in x and y
+ x = x + drphi * rn.nextGaussian();
+ y = y + dz * rn.nextGaussian();
+ }
+ // Return the smeared hit position as a new Hep3Vector
+ Hep3Vector smpos = new BasicHep3Vector(x, y, z);
+ return smpos;
+ }
+
+ /**
+ * Calculate the covariance matrix for the hit smearing
+ * @param pos Unsmeared hit position
+ * @return Covariance matrix for hit smearing
+ */
+ public SymmetricMatrix getCovariance(Hep3Vector pos) {
+ // covariance matrix
+ SymmetricMatrix cov = new SymmetricMatrix(3);
+ double tmp_drphi=0.0125;
+ double tmp_dz=0.125;
+ double x = pos.x();
+ double y = pos.y();
+ double r = Math.sqrt(x * x + y * y);
+ if (_barrel) {
+ cov.setElement(0, 0, Math.pow(y * tmp_drphi / r, 2));
+ cov.setElement(1, 0, -x * y * Math.pow(tmp_drphi / r, 2));
+ cov.setElement(1, 1, Math.pow(x * tmp_drphi / r, 2));
+ cov.setElement(2, 0, 0.);
+ cov.setElement(2, 1, 0.);
+ cov.setElement(2, 2, tmp_dz * tmp_dz);
+ return cov;
+ } else {
+ cov.setElement(0, 0, tmp_drphi*tmp_drphi);
+ cov.setElement(1, 0, 0);
+ cov.setElement(1, 1, tmp_dz*tmp_dz);
+ cov.setElement(2, 0, 0.);
+ cov.setElement(2, 1, 0.);
+ cov.setElement(2, 2, 0.0000001);
+ return cov;
+ }
+ }
+}
diff -N sATLASDigiDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ sATLASDigiDriver.java 20 Apr 2009 21:19:43 -0000 1.1 @@ -0,0 +1,21 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+
+package org.lcsim.contrib.mgraham.sATLASDigi;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author mgraham
+ */
+public class sATLASDigiDriver extends Driver {
+ public sATLASDigiDriver(){
+ MyDigiHitMaker hm=new MyDigiHitMaker();
+ add(hm);
+ LOIEffFakeMG ef=new LOIEffFakeMG();
+ add(ef);
+ }
+}