lcsim-contrib/src/main/java/org/lcsim/contrib/Partridge/Occupancy
diff -N TrackAnalysisDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackAnalysisDriver.java 7 May 2009 23:59:14 -0000 1.1
@@ -0,0 +1,344 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.Partridge.Occupancy;
+
+import java.io.IOException;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import hep.aida.IHistogramFactory;
+import hep.aida.IProfile1D;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+
+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 org.lcsim.contrib.Partridge.TrackingTest.FindableTrack.Ignore;
+//import org.lcsim.contrib.Partridge.TrackingTest.TrackAnalysis;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.tracker.silicon.ChargeCarrier;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitPixel;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author partridge
+ */
+public class TrackAnalysisDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ private IHistogramFactory _hf;
+ List<String> _process_paths = new ArrayList<String>();
+ List<IDetectorElement> _process_de = new ArrayList<IDetectorElement>();
+ Set<SiSensor> _process_sensors = new HashSet<SiSensor>();
+ public String outputPlots = "myplots.aida";
+ Map<String, IProfile1D> clsizeMap = new HashMap<String, IProfile1D>();
+ Map<String, IProfile1D> occMap = new HashMap<String, IProfile1D>();
+ String[] detNames = {"VtxPixelBarrel", "VtxPixelEndcap", "SCTShortBarrel", "SCTLongBarrel", "SCTShortEndcap", "SCTLongEndcap"};
+ Integer[] nlayers = {4, 6, 3, 2, 5, 5};
+ int trk_count = 0;
+ int nevt = 0;
+ int _nmcTrk = 0;
+ double _nrecTrk = 0;
+
+ public TrackAnalysisDriver() {
+
+ // Specify the detectors to process
+ _process_paths.add("VtxPixelBarrel");
+ _process_paths.add("SCTShortBarrel");
+ _process_paths.add("SCTLongBarrel");
+ _process_paths.add("VtxPixelEndcap");
+ _process_paths.add("SCTShortEndcap");
+ _process_paths.add("SCTLongEndcap");
+
+ // Define the efficiency histograms
+ _hf = aida.histogramFactory();
+
+
+ int i, j;
+ for (i = 0; i < 6; i++) {
+ for (j = 0; j < nlayers[i]; j++) {
+ int laynum = j + 1;
+ String profname = detNames[i] + "_layer" + laynum + " cluster size vs eta";
+ String key = detNames[i] + "_layer" + laynum;
+ clsizeMap.put(key, _hf.createProfile1D(profname, 50, -2.5, 2.5));
+ }
+ }
+
+ }
+
+ @Override
+ public void process(
+ EventHeader event) {
+
+ // Increment the event counter
+ nevt++;
+// String resDir = "residualsPlots/";
+ String simDir = "STHitPlots/";
+ String debugDir = "debugPlots/";
+ String occDir = "occupancyPlots/";
+ String occ2Dir = "RPOccupancy";
+
+ // Get the magnetic field
+ Hep3Vector IP = new BasicHep3Vector(0., 0., 0.);
+ double bfield = event.getDetector().getFieldMap().getField(IP).z();
+
+ for (SiSensor sensor : _process_sensors) {
+ SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE);
+ int nchan = electrodes.getNCells();
+ int nhits = sensor.getReadout().getHits(RawTrackerHit.class).size();
+ double occ = (double) 40 * nhits / nchan;
+ Hep3Vector pos = sensor.getGeometry().getPosition();
+ double x = pos.x();
+ double y = pos.y();
+ double z = pos.z();
+ double r = Math.sqrt(x * x + y * y);
+ int ri = 2 * ((int) Math.round(r / 20));
+ int zi = 2 * ((int) Math.round(Math.abs(z) / 20));
+ boolean barrel = !sensor.getName().toUpperCase().contains("ENDCAP");
+ String identifier;
+ if (barrel) {
+ identifier = "Barrel at r = " + ri;
+ } else {
+ identifier = "Endcap at z = " + zi;
+ }
+
+ if (!occMap.containsKey(identifier)) {
+ occMap.put(identifier, _hf.createProfile1D(identifier, 50, 0., 120.));
+ }
+
+ if (barrel) {
+ occMap.get(identifier).fill(zi, occ);
+ } else {
+ occMap.get(identifier).fill(ri, occ);
+ }
+ }
+
+ // 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, 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(simDir + input + " layer " + layer + " STHit eta").fill(eta);
+ aida.cloud1D(simDir + input + " layer " + layer + " STHit phi").fill(phi);
+ aida.cloud2D(simDir + input + " layer " + layer + " STHit phi vs eta").fill(eta, phi);
+ aida.histogram2D(simDir + 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 < 7) {
+ if (nhits[i] > 0) {
+ aida.cloud1D(simDir + input + "layer " + i + " number of ST hits").fill(nhits[i]);
+ }
+ i++;
+ }
+ }
+ List<SiTrackerHitStrip1D> stripHits = event.get(SiTrackerHitStrip1D.class, "StripClusterer_SiTrackerHitStrip1D");
+ List<SiTrackerHitPixel> pixelHits = event.get(SiTrackerHitPixel.class, "PixelClusterer_SiTrackerHitPixel");
+ List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, "RawTrackerHitMaker_RawTrackerHits");
+ List<HelicalTrackHit> hthits = event.get(HelicalTrackHit.class, "HelicalTrackHits");
+ Map<String, Integer> occupancyMap = new HashMap<String, Integer>();
+ for (RawTrackerHit rh : rawHits) {
+ IDetectorElement rhDetE = rh.getDetectorElement();
+
+ String rhDetName = rhDetE.getName();
+// System.out.println(rhDetName);
+ int rhLayer = rh.getLayerNumber();
+// String[] shortrhDetName=rhDetName.split("^[A-Z]+_layer[0-9]");
+
+
+ for (String myname : detNames) {
+ if (rhDetName.contains(myname)) {
+ String detlayer = myname + "_" + rhLayer;
+ Integer myint = occupancyMap.get(detlayer);
+ if (myint == null) {
+ myint = 1;
+ }
+ myint++;
+ occupancyMap.put(detlayer, myint);
+ }
+ }
+ }
+ Set<String> mykeyset = (Set<String>) occupancyMap.keySet();
+ for (String keys : mykeyset) {
+ aida.cloud1D(occDir + keys + " # of hits").fill(occupancyMap.get(keys));
+ }
+ int detNum = 0;
+ int layerNum = 0;
+ for (SiTrackerHitPixel stripCluster : pixelHits) {
+ Hep3Vector strCluPos = stripCluster.getPositionAsVector();
+ double rHit = Math.sqrt(strCluPos.x() * strCluPos.x() + strCluPos.y() * strCluPos.y());
+ double zHit = strCluPos.z();
+ double etaHit = -Math.log(Math.tan(Math.atan2(rHit, zHit) / 2));
+ List<RawTrackerHit> rthList = stripCluster.getRawHits();
+ int nhits = rthList.size();
+ String detlayer = "Foobar";
+ for (RawTrackerHit rth : rthList) {
+ IDetectorElement rhDetE = rth.getDetectorElement();
+ String rhDetName = rhDetE.getName();
+ int rhLayer = rth.getLayerNumber();
+ for (String myname : detNames) {
+ if (rhDetName.contains(myname)) {
+ detlayer = myname + "_layer" + rhLayer;
+ }
+ }
+ }
+// System.out.println(detlayer);
+ clsizeMap.get(detlayer).fill(etaHit, nhits);
+ aida.cloud1D(occDir + detlayer + " cluster size").fill(nhits);
+// ClusterSize[detNum][layerNum].fill(etaHit, nhits);
+
+ }
+ for (SiTrackerHitStrip1D stripCluster : stripHits) {
+ Hep3Vector strCluPos = stripCluster.getPositionAsVector();
+ double rHit = Math.sqrt(strCluPos.x() * strCluPos.x() + strCluPos.y() * strCluPos.y());
+ double zHit = strCluPos.z();
+ double etaHit = -Math.log(Math.tan(Math.atan2(rHit, zHit) / 2));
+ List<RawTrackerHit> rthList = stripCluster.getRawHits();
+ int nhits = rthList.size();
+ String detlayer = "Foobar";
+ for (RawTrackerHit rth : rthList) {
+ IDetectorElement rhDetE = rth.getDetectorElement();
+ String rhDetName = rhDetE.getName();
+ int rhLayer = rth.getLayerNumber();
+ for (String myname : detNames) {
+ if (rhDetName.contains(myname)) {
+ detlayer = myname + "_layer" + rhLayer;
+ }
+ }
+ }
+// System.out.println(detlayer);
+ clsizeMap.get(detlayer).fill(etaHit, nhits);
+ aida.cloud1D(occDir + detlayer + " cluster size").fill(nhits);
+ }
+ for (HelicalTrackHit HelTrHit : hthits) {
+ }
+
+
+ // 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());
+ }
+
+ // Now loop over all MC Particles
+
+
+
+
+
+ return;
+ }
+
+ @Override
+ public void endOfData() {
+ try {
+ aida.saveAs(outputPlots);
+ } catch (IOException ex) {
+ Logger.getLogger(TrackAnalysisDriver.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ System.out.println("# of reco tracks = " + _nrecTrk + "; # of MC tracks = " + _nmcTrk + "; Efficiency = " + _nrecTrk / _nmcTrk);
+ }
+
+ public void setOutputPlots(String output) {
+ this.outputPlots = output;
+
+ }
+
+ 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;
+ }
+
+ public void detectorChanged(Detector detector) {
+ System.out.println(detector.getName());
+ super.detectorChanged(detector);
+
+ // Process detectors specified by path, otherwise process entire detector
+ IDetectorElement detector_de = detector.getDetectorElement();
+ System.out.println(detector_de.getName());
+ for (String de_path : _process_paths) {
+ _process_de.add(detector_de.findDetectorElement(de_path));
+ }
+
+ if (_process_de.size() == 0) {
+ _process_de.add(detector_de);
+ }
+
+
+
+
+ for (IDetectorElement detector_element : _process_de) {
+ _process_sensors.addAll(detector_element.findDescendants(SiSensor.class));
+ }
+
+ }
+}
+
+
lcsim-contrib/src/main/java/org/lcsim/contrib/Partridge/Occupancy
diff -N TrackerHitDriver_sATLAS.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackerHitDriver_sATLAS.java 7 May 2009 23:59:14 -0000 1.1
@@ -0,0 +1,216 @@
+/*
+ * TrackerHitDriver Class
+ *
+ */
+package org.lcsim.contrib.Partridge.Occupancy;
+
+import org.lcsim.contrib.sATLAS.*;
+import java.util.ArrayList;
+import java.util.HashSet;
+import java.util.List;
+import java.util.Set;
+
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiTrackerModule;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.tracking.digitization.sisim.CDFSiSensorSim;
+import org.lcsim.recon.tracking.digitization.sisim.GenericReadoutChip;
+import org.lcsim.recon.tracking.digitization.sisim.NearestNeighbor;
+import org.lcsim.recon.tracking.digitization.sisim.PixelHitMaker;
+import org.lcsim.recon.tracking.digitization.sisim.RawTrackerHitMaker;
+import org.lcsim.recon.tracking.digitization.sisim.SiDigitizer;
+import org.lcsim.recon.tracking.digitization.sisim.SiSensorSim;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitPixel;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+import org.lcsim.recon.tracking.digitization.sisim.SimTrackerHitReadoutDriver;
+import org.lcsim.recon.tracking.digitization.sisim.StripHitMaker;
+import org.lcsim.util.Driver;
+import org.lcsim.util.lcio.LCIOConstants;
+
+/**
+ *
+ * @author tknelson
+ */
+public class TrackerHitDriver_sATLAS extends Driver {
+
+ List<String> _readouts = new ArrayList<String>();
+ List<String> _process_paths = new ArrayList<String>();
+ List<IDetectorElement> _process_de = new ArrayList<IDetectorElement>();
+ Set<SiSensor> _process_sensors = new HashSet<SiSensor>();
+ Set<SiTrackerModule> _process_modules = new HashSet<SiTrackerModule>();
+ SiDigitizer _digitizer;
+ StripHitMaker _strip_clusterer;
+ PixelHitMaker _pixel_clusterer;
+ int _nev = 0;
+
+ /**
+ * Creates a new instance of TrackerHitDriver
+ */
+ public TrackerHitDriver_sATLAS() {
+
+ // Instantiate the sensor simulation classes and set the thresholds
+ SiSensorSim strip_simulation = new CDFSiSensorSim();
+ SiSensorSim pixel_simulation = new CDFSiSensorSim();
+
+ // Instantiate the readout chips and set the noise parameters
+ GenericReadoutChip strip_readout = new GenericReadoutChip();
+ strip_readout.setNoiseIntercept(300.);
+ strip_readout.setNoiseSlope(30.);
+ strip_readout.setNoiseThreshold(4000.);
+ strip_readout.setNeighborThreshold(4000.);
+ GenericReadoutChip pixel_readout = new GenericReadoutChip();
+ pixel_readout.setNoiseIntercept(300.);
+ pixel_readout.setNoiseSlope(30.);
+ pixel_readout.setNoiseThreshold(4000.);
+ pixel_readout.setNeighborThreshold(4000.);
+
+ // Instantiate the digitizer that produces the raw hits
+ _digitizer = new RawTrackerHitMaker(strip_simulation, strip_readout);
+
+ // Instantiate a nearest neighbor clustering algorithm for the pixels
+ NearestNeighbor strip_clustering = new NearestNeighbor();
+ strip_clustering.setSeedThreshold(4000.);
+ strip_clustering.setNeighborThreshold(4000.);
+
+ // Instantiate a nearest neighbor clustering algorithm for the pixels
+ NearestNeighbor pixel_clustering = new NearestNeighbor();
+ pixel_clustering.setSeedThreshold(4000.);
+ pixel_clustering.setNeighborThreshold(4000.);
+
+ // Instantiate the clusterers and set hit-making parameters
+ _strip_clusterer = new StripHitMaker(strip_simulation, strip_readout, strip_clustering);
+ _strip_clusterer.setMaxClusterSize(10);
+ _strip_clusterer.setCentralStripAveragingThreshold(4);
+ _pixel_clusterer = new PixelHitMaker(pixel_simulation, pixel_readout, pixel_clustering);
+
+ // Specify the readouts to process
+ _readouts.add("VtxBarrHits");
+ _readouts.add("SCTShortBarrHits");
+ _readouts.add("SCTLongBarrHits");
+ _readouts.add("VtxEndcapHits");
+ _readouts.add("SCTShortEndcapHits");
+ _readouts.add("SCTLongEndcapHits");
+
+ // Specify the detectors to process
+ _process_paths.add("VtxPixelBarrel");
+ _process_paths.add("SCTShortBarrel");
+ _process_paths.add("SCTLongBarrel");
+ _process_paths.add("VtxPixelEndcap");
+ _process_paths.add("SCTShortEndcap");
+ _process_paths.add("SCTLongEndcap");
+
+ }
+
+ /**
+ * Initialize whenever we have a new detector
+ *
+ * @param detector
+ */
+ public void detectorChanged(Detector detector) {
+ System.out.println(detector.getName());
+ super.detectorChanged(detector);
+
+ // Process detectors specified by path, otherwise process entire detector
+ IDetectorElement detector_de = detector.getDetectorElement();
+ System.out.println(detector_de.getName());
+ for (String de_path : _process_paths) {
+ _process_de.add(detector_de.findDetectorElement(de_path));
+ }
+
+ if (_process_de.size() == 0) {
+ _process_de.add(detector_de);
+ }
+
+ for (IDetectorElement detector_element : _process_de) {
+ _process_sensors.addAll(detector_element.findDescendants(SiSensor.class));
+ _process_modules.addAll(detector_element.findDescendants(SiTrackerModule.class));
+ }
+
+ }
+
+ /**
+ * Setup readouts
+ */
+ public void startOfData() {
+ // If readouts not already set, set them up
+ if (_readouts.size() != 0) {
+ System.out.println("Adding SimTrackerHitIdentifierReadoutDriver with readouts: " + _readouts);
+ super.add(new SimTrackerHitReadoutDriver(_readouts));
+ }
+
+ super.startOfData();
+ _readouts.clear();
+ _nev = 0;
+ }
+
+ /**
+ * Main digitization driver. Creates raw hits, forms clusters, and makes
+ * tracker hits using the sisim package.
+ *
+ * @param event
+ */
+ public void process(EventHeader event) {
+ super.process(event);
+
+ // Print out the event number
+ System.out.println("TrackerHitDriver processing event "+_nev);
+ _nev++;
+
+ // Lists of hits
+ List<RawTrackerHit> raw_hits = new ArrayList<RawTrackerHit>();
+ List<SiTrackerHit> hits_strip1D = new ArrayList<SiTrackerHit>();
+ List<SiTrackerHit> hits_pixel = new ArrayList<SiTrackerHit>();
+
+ for (SiSensor sensor : _process_sensors) {
+ raw_hits.addAll(_digitizer.makeHits(sensor));
+
+ if (sensor.hasStrips()) {
+ hits_strip1D.addAll(_strip_clusterer.makeHits(sensor));
+ }
+
+
+ if (sensor.hasPixels()) {
+ hits_pixel.addAll(_pixel_clusterer.makeHits(sensor));
+ }
+
+ }
+
+ int flag = (1 << LCIOConstants.RTHBIT_HITS | 1 << LCIOConstants.TRAWBIT_ID1); //correct flag for persistence
+ event.put(getRawHitsName(), raw_hits, RawTrackerHit.class, flag, toString());
+ event.put(getStripHits1DName(), hits_strip1D, SiTrackerHitStrip1D.class, 0, toString());
+ event.put(getPixelHitsName(), hits_pixel, SiTrackerHitPixel.class, 0, toString());
+
+ }
+
+ /**
+ * Return the name of the raw hits collection
+ *
+ * @return name of raw hits collection
+ */
+ public String getRawHitsName() {
+ return _digitizer.getName() + "_RawTrackerHits";
+ }
+
+ /**
+ * Return the name of the strip hits collection
+ *
+ * @return name of strip hits collection
+ */
+ public String getStripHits1DName() {
+ return _strip_clusterer.getName() + "_SiTrackerHitStrip1D";
+ }
+
+ /**
+ * Return the name of the pixel hits collection
+ *
+ * @return name of pixel hits collection
+ */
+ public String getPixelHitsName() {
+ return _pixel_clusterer.getName() + "_SiTrackerHitPixel";
+ }
+
+}
\ No newline at end of file
