 |
 |
| Commit in lcsim-contrib/src/main/java/org/lcsim/contrib/sATLAS/UTOPIA7 on MAIN | |||
| SiVertexEndcapSensorSetup.java | +138 | added 1.1 | |
| MakeSensorsDriver.java | +22 | added 1.1 | |
| SiTrackerEndcap2SensorSetup.java | +107 | added 1.1 | |
| TrackAnalysisDriver.java | +862 | added 1.1 | |
| TrackReconstructionDriver.java | +47 | added 1.1 | |
| StandAloneDriver.java | +28 | added 1.1 | |
| OccupancyDriver.java | +230 | added 1.1 | |
| SiVertexBarrelSensorSetup.java | +100 | added 1.1 | |
| sATLASMainDriver.java | +40 | added 1.1 | |
| SiTrackerBarrelSensorSetup.java | +100 | added 1.1 | |
| TrackerHitDriver_sATLAS.java | +243 | added 1.1 | |
| +1917 | |||
sATLAS UTOPIA7 drivers
diff -N SiVertexEndcapSensorSetup.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ SiVertexEndcapSensorSetup.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,138 @@
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+import hep.physics.matrix.BasicMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.util.List;
+
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IRotation3D;
+import org.lcsim.detector.ITranslation3D;
+import org.lcsim.detector.RotationPassiveXYZ;
+import org.lcsim.detector.Transform3D;
+import org.lcsim.detector.Translation3D;
+import org.lcsim.detector.solids.IPolyhedron;
+import org.lcsim.detector.solids.Polygon3D;
+import org.lcsim.detector.tracker.silicon.ChargeCarrier;
+import org.lcsim.detector.tracker.silicon.SiPixels;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.compact.Subdetector;
+import org.lcsim.geometry.subdetector.SiTrackerEndcap;
+import org.lcsim.geometry.subdetector.SiTrackerEndcap2;
+import org.lcsim.util.Driver;
+
+public class SiVertexEndcapSensorSetup extends Driver {
+
+ String subdetectorName;
+
+ public SiVertexEndcapSensorSetup() {
+ }
+
+ public SiVertexEndcapSensorSetup(String subdetectorName) {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void setSubdetectorName(String subdetectorName) {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void detectorChanged(Detector detector) {
+ if (subdetectorName == null) {
+ throw new RuntimeException("The subdetectorName was not set.");
+ }
+
+ Subdetector subdetector = detector.getSubdetector(subdetectorName);
+ if (subdetector instanceof SiTrackerEndcap || subdetector instanceof SiTrackerEndcap2) {
+ setupSensorDetectorElements(subdetector);
+ } else {
+ throw new RuntimeException("The subdetector " + subdetectorName + " is not an instance of SiTrackerBarrel.");
+ }
+ }
+
+ private void setupSensorDetectorElements(Subdetector subdet) {
+ System.out.println(this.getClass().getCanonicalName() + " - Setting up sensors for " + subdet.getName() + " ...");
+ for (IDetectorElement endcap : subdet.getDetectorElement().getChildren()) {
+ for (IDetectorElement layer : endcap.getChildren()) {
+ int nwedges = layer.getChildren().size();
+ for (IDetectorElement wedge : layer.getChildren()) {
+ for (IDetectorElement module : wedge.getChildren()) {
+ // find sensors on the module
+ List<SiSensor> sensors = module.findDescendants(SiSensor.class);
+
+ // require that sensors are found
+ if (sensors.size() == 0) {
+ throw new RuntimeException("No sensors found in module.");
+ }
+
+ // loop over sensors (can be double-sided)
+ for (SiSensor sensor : sensors) {
+ // get sensor field from id
+ int sensorId = sensor.getIdentifierHelper().getValue(sensor.getIdentifier(), "sensor");
+
+ // Get the sensor solid.
+ IPolyhedron sensor_solid = (IPolyhedron) sensor.getGeometry().getLogicalVolume().getSolid();
+
+ // Get solids for inner and outer surfaces.
+ Polygon3D inner_surface = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, -1, 0)).get(0);
+ Polygon3D outer_surface = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, 1, 0)).get(0);
+
+ //
+ // Determine p and n sides based on sensor id.
+ //
+/*
+ Polygon3D p_side;
+ Polygon3D n_side;
+ int side;
+
+ if (sensorId == 0) // inner sensor
+ {
+ p_side = inner_surface;
+ n_side = outer_surface;
+ side = 1;
+ } else // outer sensor
+ {
+ p_side = outer_surface;
+ n_side = inner_surface;
+ side = -1;
+ }
+ */
+ Polygon3D n_side = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, -1, 0)).get(0);
+ Polygon3D p_side = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, 1, 0)).get(0);
+ double strip_angle = Math.PI / nwedges;
+
+ // Compute the geometric propertes of the electrodes.
+ ITranslation3D electrodes_position = new Translation3D(VecOp.mult(-p_side.getDistance(), new BasicHep3Vector(0, 0, 1))); // translate to outside of polygon
+// IRotation3D electrodes_rotation = new RotationPassiveXYZ(side * (Math.PI / 2), 0, strip_angle); //
+ IRotation3D electrodes_rotation = new RotationPassiveXYZ(-(Math.PI / 2), 0, 0); // Do Not Rotate Pixels
+ Transform3D electrodes_transform = new Transform3D(electrodes_position, electrodes_rotation);
+
+ //
+ // Pixel-specific code starts here.
+ //
+
+ // Set the bias surfaces.
+ sensor.setBiasSurface(ChargeCarrier.ELECTRON, p_side);
+ sensor.setBiasSurface(ChargeCarrier.HOLE, n_side);
+
+ // Define the pixel electrodes.
+ SiSensorElectrodes readout_electrodes = new SiPixels(ChargeCarrier.ELECTRON, 0.05, 0.25, sensor, electrodes_transform);
+ SiSensorElectrodes sense_electrodes = new SiPixels(ChargeCarrier.ELECTRON, 0.05, 0.25, sensor, electrodes_transform);
+
+ // Tell the sensor about the electrodes.
+ sensor.setSenseElectrodes(sense_electrodes);
+ sensor.setReadoutElectrodes(readout_electrodes);
+
+ // Define the transfer efficiency from sense electrodes to readout electrodes.
+ // For pixels, we do a direct transfer of charge from sense electrodes to readout electrodes.
+ double[][] transfer_efficiencies = {{1.0}};
+ sensor.setTransferEfficiencies(ChargeCarrier.ELECTRON, new BasicMatrix(transfer_efficiencies));
+ }
+ }
+ }
+ }
+ }
+ }
+}
\ No newline at end of file
diff -N MakeSensorsDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ MakeSensorsDriver.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,22 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author mgraham
+ */
+public class MakeSensorsDriver extends Driver {
+
+ public MakeSensorsDriver() {
+ add(new SiVertexBarrelSensorSetup("VtxPixelBarrel"));
+ add(new SiTrackerBarrelSensorSetup("SCTShortBarrel"));
+ add(new SiTrackerBarrelSensorSetup("SCTLongBarrel"));
+ add(new SiVertexEndcapSensorSetup("VtxPixelEndcap"));
+ add(new SiTrackerEndcap2SensorSetup("SCTEndcap"));
+ }
+}
diff -N SiTrackerEndcap2SensorSetup.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ SiTrackerEndcap2SensorSetup.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,107 @@
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+import hep.physics.matrix.BasicMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.util.List;
+
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IRotation3D;
+import org.lcsim.detector.ITranslation3D;
+import org.lcsim.detector.RotationPassiveXYZ;
+import org.lcsim.detector.Transform3D;
+import org.lcsim.detector.Translation3D;
+import org.lcsim.detector.solids.Polygon3D;
+import org.lcsim.detector.solids.Trd;
+import org.lcsim.detector.tracker.silicon.ChargeCarrier;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
+import org.lcsim.detector.tracker.silicon.SiStrips;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.compact.Subdetector;
+import org.lcsim.geometry.subdetector.SiTrackerEndcap2;
+import org.lcsim.util.Driver;
+
+public class SiTrackerEndcap2SensorSetup extends Driver {
+
+ String subdetectorName;
+
+ public SiTrackerEndcap2SensorSetup() {
+ }
+
+ public SiTrackerEndcap2SensorSetup(String subdetectorName) {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void setSubdetectorName(String subdetectorName) {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void detectorChanged(Detector detector) {
+ if (subdetectorName == null) {
+ throw new RuntimeException("The subdetectorName was not set.");
+ }
+
+ Subdetector subdetector = detector.getSubdetector(subdetectorName);
+ if (subdetector instanceof SiTrackerEndcap2) {
+ setupSensorDetectorElements(subdetector);
+ } else {
+ throw new RuntimeException("The subdetector " + subdetectorName + " is not an instance of SiTrackerEndcap.");
+ }
+ }
+
+ private void setupSensorDetectorElements(Subdetector subdet) {
+ System.out.println(this.getClass().getCanonicalName() + " - Setting up sensors for " + subdet.getName() + " ...");
+
+ for (IDetectorElement endcap : subdet.getDetectorElement().getChildren()) {
+ for (IDetectorElement layer : endcap.getChildren()) {
+ //int nwedges = layer.getChildren().size();
+ for (IDetectorElement wedge : layer.getChildren()) {
+ for (IDetectorElement module : wedge.getChildren()) {
+ List<SiSensor> sensors = module.findDescendants(SiSensor.class);
+
+ if (sensors.size() == 0) {
+ throw new RuntimeException("No sensors found in module.");
+ }
+
+ //int sensorId = 0;
+ for (SiSensor sensor : sensors) {
+ Trd sensor_solid = (Trd) sensor.getGeometry().getLogicalVolume().getSolid();
+
+ Polygon3D n_side = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, -1, 0)).get(0);
+ Polygon3D p_side = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, 1, 0)).get(0);
+
+ // Bias the sensor
+ sensor.setBiasSurface(ChargeCarrier.HOLE, p_side);
+ sensor.setBiasSurface(ChargeCarrier.ELECTRON, n_side);
+
+ double strip_angle = Math.atan2(sensor_solid.getXHalfLength2() - sensor_solid.getXHalfLength1(), sensor_solid.getZHalfLength() * 2);
+
+ ITranslation3D electrodes_position = new Translation3D(VecOp.mult(-p_side.getDistance(), new BasicHep3Vector(0, 0, 1))); // translate to outside of polygon
+ IRotation3D electrodes_rotation = new RotationPassiveXYZ(-Math.PI / 2, 0, strip_angle);
+ Transform3D electrodes_transform = new Transform3D(electrodes_position, electrodes_rotation);
+
+ // Free calculation of readout electrodes, sense electrodes determined thereon
+ SiSensorElectrodes readout_electrodes = new SiStrips(ChargeCarrier.HOLE, 0.075, sensor, electrodes_transform);
+// SiSensorElectrodes sense_electrodes = new SiStrips(ChargeCarrier.HOLE,0.025,(readout_electrodes.getNCells()*2-1),sensor,electrodes_transform);
+ SiSensorElectrodes sense_electrodes = new SiStrips(ChargeCarrier.HOLE, 0.075, sensor, electrodes_transform);
+
+ sensor.setSenseElectrodes(sense_electrodes);
+ sensor.setReadoutElectrodes(readout_electrodes);
+
+ double[][] transfer_efficiencies = {{1.0}};
+//double[][] transfer_efficiencies = { {0.986,0.419} };
+ sensor.setTransferEfficiencies(ChargeCarrier.HOLE, new BasicMatrix(transfer_efficiencies));
+ // here
+
+ //++sensorId;
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+
diff -N TrackAnalysisDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TrackAnalysisDriver.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,862 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+import java.io.IOException;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.contrib.mgraham.sATLASDigi.*;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+
+import hep.physics.vec.VecOp;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+import hep.aida.*;
+
+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.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.digitization.sisim.SiTrackerHitPixel;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+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;
+
+/**
+ *
+ * @author partridge
+ */
+public class TrackAnalysisDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+ private IHistogram1D pTeff1;
+ private IHistogram1D pTeff2;
+ private IHistogram1D thetaeff;
+ private IHistogram1D ctheff;
+ private IHistogram1D etaeff;
+ private IHistogram1D d0eff1;
+ private IHistogram1D z0eff1;
+ private IHistogram1D z0eff2;
+ private IHistogram1D pTeff1Findable;
+ private IHistogram1D pTeff2Findable;
+ private IHistogram1D thetaeffFindable;
+ private IHistogram1D ctheffFindable;
+ private IHistogram1D etaeffFindable;
+ private IHistogram1D d0eff1Findable;
+ private IHistogram1D d0eff2Findable;
+ private IHistogram1D z0eff1Findable;
+ private IHistogram1D z0eff2Findable;
+ private IHistogram1D fakes;
+ private IHistogram1D nfakes;
+ private IHistogram1D etafake;
+ public String outputPlots = "myplots.aida";
+ Map<String, IProfile1D> clsizeMap = new HashMap<String, IProfile1D>();
+// String[] detNames = {"VtxPixelBarrel", "VtxPixelEndcap", "SCTShortBarrel", "SCTLongBarrel", "SCTShortEndcap", "SCTLongEndcap"};
+ String[] detNames = {"VtxPixelBarrel", "VtxPixelEndcap", "SCTShortBarrel", "SCTLongBarrel", "SCTEndcap", };
+// Integer[] nlayers = {4, 6, 3, 2, 5, 5};
+ Integer[] nlayers = {4, 6, 3, 2, 5};
+ int trk_count = 0;
+ int nevt = 0;
+ int _nmcTrk = 0;
+ double _nrecTrk = 0;
+ IAnalysisFactory af = IAnalysisFactory.create();
+ ITree tree = af.createTreeFactory().create();
+ ITupleFactory tf = af.createTupleFactory(tree);
+ String columnString = "int ievt=0 , nMCPart=0; ITuple MCInfo={double pXMC=0, pYMC=0, pXMC=0};";
+ ITuple tuple = tf.create("tuple", "MyNtpule", columnString);
+
+ public TrackAnalysisDriver() {
+
+ // 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.5, 2.5, "type=efficiency");
+ d0eff1 = hf.createHistogram1D("Efficiency vs d0", "", 50, -2., 2., "type=efficiency");
+
+ z0eff1 = hf.createHistogram1D("Efficiency vs z0", "", 50, -50., 50., "type=efficiency");
+ z0eff2 = hf.createHistogram1D("Efficiency vs z0 full", "", 50, -200., 200., "type=efficiency");
+
+ pTeff1Findable = hf.createHistogram1D("Findable Efficiency vs pT", "", 100, 0., 5., "type=efficiency");
+ pTeff2Findable = hf.createHistogram1D("Findable Efficiency vs pT full", "", 100, 0., 50., "type=efficiency");
+ thetaeffFindable = hf.createHistogram1D("Findable Efficiency vs theta", "", 72, 0., 180., "type=efficiency");
+ ctheffFindable = hf.createHistogram1D("Findable Efficiency vs cos(theta)", "", 50, -1., 1., "type=efficiency");
+ etaeffFindable = hf.createHistogram1D("Findable Efficiency vs eta", "", 50, -2.5, 2.5, "type=efficiency");
+ d0eff1Findable = hf.createHistogram1D("Findable Efficiency vs d0", "", 50, -0.5, 0.5, "type=efficiency");
+ d0eff2Findable = hf.createHistogram1D("Findable Efficiency vs d0 full", "", 50, -5., 5., "type=efficiency");
+ z0eff1Findable = hf.createHistogram1D("Findable Efficiency vs z0", "", 50, -50., 50., "type=efficiency");
+ z0eff2Findable = hf.createHistogram1D("Findable Efficiency vs z0 full", "", 50, -200., 200., "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.);
+ etafake = hf.createHistogram1D("Fake rate vs eta", "", 50, -2.5, 2.5, "type=efficiency");
+ int i, j;
+// for (i = 0; i < 6; i++) {
+ for (i = 0; i < 5; 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, 25, -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/";
+ // 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"};
+ String[] input_hit_collections = {"VtxBarrHits", "VtxEndcapHits", "SCTShortBarrHits", "SCTLongBarrHits", "SCTEndcapHits"};
+ 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) {
+
+ //System.out.println("Number of pixels: "+stripCluster.getReadoutElectrodes().getNCells());
+ stripCluster.getSensor().getReadout().getHits(RawTrackerHit.class); //gets all hits on sensor
+ 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();
+ //System.out.println("Number of rawhits: "+nhits);
+ String detlayer = "Foobar";
+ for (RawTrackerHit rth : rthList) {
+ IDetectorElement rhDetE = rth.getDetectorElement();
+ String rhDetName = rhDetE.getName();
+// System.out.println(rhDetName);
+ int rhLayer = rth.getLayerNumber();
+ String[] detLayerName = rhDetName.split("_module");
+// System.out.println(detLayerName[0]);
+ 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) {
+ }
+
+ // Get the list of strategies being used
+// String sfile = "autogen_ttbar_sid02_vs.xml";
+// List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(
+// StrategyXMLUtils.getDefaultStrategiesPrefix() + sfile);
+
+// String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "sATLASFull-JeffMarch26.xml";
+// String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "sATLASFull-Simple.xml";
+ String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "sATLASFull-UTOPIA7.xml";
+ List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(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(debugDir + "Track Chi2-Circle Fit").fill(chisq[0]);
+ aida.cloud1D(debugDir + "Track Chi2-RZ Fit").fill(chisq[1]);
+ aida.cloud1D(debugDir + "NH Track Chi2").fill(nhchisq);
+ if (nhchisq != 0) {
+ aida.cloud1D(debugDir + "NH!=0 Track Chi2-Circle Fit").fill(chisq[0]);
+ aida.cloud1D(debugDir + "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(debugDir + 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(debugDir + 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 umeasErr = cluster.du();
+
+ 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(debugDir + hit.Detector() + " associated ST hits").fill(nsthits);
+ aida.cloud1D(debugDir + 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);
+ }
+ }
+
+
+ // System.out.println("filling...");
+ if (umc != -999999) {
+ aida.cloud2D(debugDir + hit.Detector() + "cluster vs STHit dedx").fill(stenergy, charge);
+ aida.cloud2D(debugDir + 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(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u)").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u) Error").fill(umeasErr);
+ if (nstrips == 1) {
+ aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)--1 strip").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u)--1 strip").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u) Error--1 strip").fill(umeasErr);
+ }
+ if (nstrips == 2) {
+ aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)--2 strip").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u)--2 strip").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u) Error--2 strip").fill(umeasErr);
+ }
+ if (nstrips == 3) {
+ aida.cloud1D(debugDir + hit.Detector() + "layer=" + stripdir + " delta(u)--3 strip").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u)--3 strip").fill(umeas - umc);
+ aida.cloud1D(debugDir + hit.Detector() + " delta(u) Error--3 strip").fill(umeasErr);
+ }
+ }
+
+ }
+ aida.cloud2D(debugDir + 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(resDir + hit.Detector() + " dxdy").fill(dxdy);
+ aida.cloud1D(resDir + hit.Detector() + " dz").fill(dz);
+ aida.cloud1D(resDir + hit.Detector() + " dxdy Pull").fill(dxdy / dxdyErr);
+ aida.cloud1D(resDir + hit.Detector() + " dz Pull").fill(dz / dzErr);
+ aida.cloud2D(resDir + hit.Detector() + " dz vs hit z").fill(hitpos.z(), dz);
+ aida.cloud2D(resDir + hit.Detector() + " dz Pull vs hit z").fill(hitpos.z(), dz / dzErr);
+ if (Math.abs(dz) > 4) {
+ aida.cloud1D(debugDir + hit.Detector() + "Bad dz--nHits").fill(nhits);
+ }
+ aida.cloud1D(debugDir + "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.cloud1D("Mis-matched hits for all tracks").fill(nbad);
+ aida.cloud1D("Mis-matched hits " + nhits + " hit tracks").fill(nbad);
+
+ // Generate a normalized histogram after 1000 events
+ trk_count++;
+ if (nevt <= 1000) {
+ fakes.fill(nbad);
+ }
+
+ // Make plots for fake, non-fake, and all tracks
+ if (purity < 0.5) {
+ aida.cloud1D("Hits for fake tracks").fill(nhits);
+ aida.cloud1D("pT for fake tracks").fill(pt);
+ aida.cloud1D("cos(theta) for fake tracks").fill(cth);
+ aida.cloud1D("d0 for fake tracks").fill(d0);
+ aida.cloud1D("z0 for fake tracks").fill(z0);
+ aida.cloud1D("eta for fake tracks").fill(eta);
+ etafake.fill(eta, 1.0);
+ } else {
+ aida.cloud1D("Hits for non-fake tracks").fill(nhits);
+ aida.cloud1D("pT for non-fake tracks").fill(pt);
+ aida.cloud1D("cos(theta) for non-fake tracks").fill(cth);
+ aida.cloud1D("d0 for non-fake tracks").fill(d0);
+ aida.cloud1D("z0 for non-fake tracks").fill(z0);
+ aida.cloud1D("eta for non-fake tracks").fill(eta);
+ etafake.fill(eta, 0.0);
+ }
+ aida.cloud1D("Hits for all tracks").fill(nhits);
+ aida.cloud1D("pT for all tracks").fill(pt);
+ aida.cloud1D("cos(theta) for all tracks").fill(cth);
+ aida.cloud1D("d0 for all tracks").fill(d0);
+ aida.cloud1D("z0 for all tracks").fill(z0);
+ aida.cloud1D("eta for all tracks").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(resDir + "pT MC vs pT Reco for 0 Bad Hits",
+ 100, 0., 5., 100, 0., 5.).fill(ptmc, pttk);
+ aida.histogram2D(resDir + "d0 MC vs d0 Reco for 0 Bad Hits",
+ 100, -0.2, 0.2, 100, -0.2, 0.2).fill(d0mc, d0tk);
+ aida.cloud1D(resDir + "pT Pull for 0 Bad Hits").fill(ptpull);
+ aida.cloud1D(resDir + "d0 pull for 0 Bad Hits").fill(d0pull);
+ aida.cloud1D(resDir + "pT Residual for 0 Bad Hits").fill(ptresid);
+ aida.cloud1D(resDir + "d0 Residual for 0 Bad Hits").fill(d0resid);
+ aida.cloud1D(resDir + "1/pT for 0 Bad Hits").fill(1 / pttk);
+
+ } else if (purity > 0.5) {
+ aida.histogram2D(resDir + "pT MC vs pT Reco for 0.5 < purity < 1",
+ 100, 0., 5., 100, 0., 5.).fill(ptmc, pttk);
+ aida.histogram2D(resDir + "d0 MC vs d0 Reco for 0.5 < purity < 1",
+ 100, -0.2, 0.2, 100, -0.2, 0.2).fill(d0mc, d0tk);
+ aida.cloud1D(resDir + "pT Pull for 0.5 < purity < 1").fill(ptpull);
+ aida.cloud1D(resDir + "d0 pull for 0.5 < purity < 1").fill(d0pull);
+ aida.cloud1D(resDir + "pT Residual for 0.5 < purity < 1").fill(ptresid);
+ aida.cloud1D(resDir + "d0 Residual for 0.5 < purity < 1").fill(d0resid);
+ aida.cloud1D(resDir + "1/pT for 0.5 < purity < 1").fill(1 / pttk);
+ } else if (purity < 0.5) {
+ aida.histogram2D(resDir + "pT MC vs pT Reco for purity <= 0.5",
+ 100, 0., 5., 100, 0., 5.).fill(ptmc, pttk);
+ aida.histogram2D(resDir + "d0 MC vs d0 Reco for purity <= 0.5",
+ 100, -0.2, 0.2, 100, -0.2, 0.2).fill(d0mc, d0tk);
+ aida.cloud1D(resDir + "pT Pull for purity <= 0.5").fill(ptpull);
+ aida.cloud1D(resDir + "d0 pull for purity <= 0.5").fill(d0pull);
+ aida.cloud1D(resDir + "pT Residual for purity <= 0.5").fill(ptresid);
+ aida.cloud1D(resDir + "d0 Residial for purity <= 0.5").fill(d0resid);
+ aida.cloud1D(resDir + "1/pT for purity <= 0.5").fill(1 / pttk);
+ }
+ }
+ }
+
+ // Make the normalized fake plot after the specified number of events
+ if (nevt == 1000) {
+ double wgt = 1. / trk_count;
+ 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 = " + trk_count);
+ }
+
+ // 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.;
+ }
+ pTeff1Findable.fill(pt, wgt);
+ pTeff2Findable.fill(pt, wgt);
+ }
+
+ // Make angular efficiency plot
+ if (findable.isFindable(mcp, slist)) {
+ double wgt = 0.;
+ if (ntrk > 0) {
+ wgt = 1.;
+ } else {
+ System.out.println("Findable Track Not Found! eta=" + eta);
+ }
+ thetaeffFindable.fill(theta, wgt);
+ ctheffFindable.fill(cth, wgt);
+ etaeffFindable.fill(eta, wgt);
+
+ }
+
+ // Make d0 efficiency plot
+ if (findable.isFindable(mcp, slist, Ignore.NoDCACut)) {
+ double wgt = 0.;
+ if (ntrk > 0) {
+ wgt = 1.;
+ }
+ d0eff1Findable.fill(d0, wgt);
+ d0eff2Findable.fill(d0, wgt);
+ }
+
+ // Make z0 efficiency plot
+ if (findable.isFindable(mcp, slist, Ignore.NoZ0Cut)) {
+ double wgt = 0.;
+ if (ntrk > 0) {
+ wgt = 1.;
+ }
+ z0eff1Findable.fill(z0, wgt);
+ z0eff2Findable.fill(z0, wgt);
+ }
+ */
+ // Select charged MC particles
+ if (mcp.getCharge() == 0) {
+ continue;
+ }
+// make the true efficiency plots
+ double ptTrkCut = 1.0; //GeV
+ double d0TrkCut = 2.0; //mm
+ double z0TrkCut = 200.0; //mm
+ double etaTrkCut = 2.5;
+// System.out.println("Final Stat Part? "+mcp.FINAL_STATE+"; pt = "+pt+"; d0 = "+d0);
+ if (pt > ptTrkCut && mcp.getGeneratorStatus() == mcp.FINAL_STATE && Math.abs(d0) < d0TrkCut && Math.abs(eta) < etaTrkCut && Math.abs(z0) < z0TrkCut) {
+ double wgt = 0.0;
+ if (ntrk > 0) {
+ wgt = 1.0;
+// System.out.println("found track!");a
+ System.out.println("Found this track! eta = " + eta + "; pT = " + pt + "; z0 = " + z0 + "; d0 = " + d0);
+ } else {
+ System.out.println("Missed this track! eta = " + eta + "; pT = " + pt + "; z0 = " + z0 + "; d0 = " + d0);
+ }
+
+ pTeff1.fill(pt, wgt);
+ pTeff2.fill(pt, wgt);
+ ctheff.fill(cth, wgt);
+ thetaeff.fill(theta, wgt);
+ etaeff.fill(eta, wgt);
+ d0eff1.fill(d0, wgt);
+ z0eff1.fill(z0, wgt);
+ z0eff2.fill(z0, wgt);
+ if (eta < etaTrkCut) {
+ _nmcTrk++;
+ _nrecTrk += wgt;
+ }
+ }
+
+
+ if (mcp.getGeneratorStatus() == mcp.FINAL_STATE) {
+ aida.cloud1D("findable/eta for final state particles").fill(eta);
+ }
+ if (mcp.getGeneratorStatus() != mcp.FINAL_STATE && mcp.getGeneratorStatus() != mcp.INTERMEDIATE) {
+ aida.cloud1D("findable/eta for other particles").fill(eta);
+ }
+ // Select mcp that fail the final state requirement
+ if (mcp.getGeneratorStatus() != mcp.FINAL_STATE) {
+ aida.cloud1D("findable/Hits for non-final state particles").fill(nhits);
+ aida.cloud1D("findable/pT for non-final state particles").fill(pt);
+ aida.cloud1D("findable/cos(theta) for non-final state particles").fill(cth);
+ aida.cloud1D("findable/eta for non-final state particles").fill(eta);
+ aida.cloud1D("findable/d0 for non-final state particles").fill(d0);
+ aida.cloud1D("findable/z0 for non-final state particles").fill(z0);
+ aida.cloud2D("findable/Hits vs eta for non-final state particles").fill(eta, nhits);
+ double zOrig = mcp.getOriginZ();
+ double rOrig = Math.sqrt(mcp.getOriginX() * mcp.getOriginX() + mcp.getOriginY() * mcp.getOriginY());
+ // aida.histogram2D("r vs z for non-final state particles",3000,0,3000,1200,0,1200).fill(Math.abs(zOrig),rOrig);
+ continue;
+ }
+
+ // Make plots for the base sample
+ aida.cloud1D("findable/Hits for base MC selection").fill(nhits);
+ aida.cloud1D("findable/pT for base MC selection").fill(pt);
+ aida.cloud1D("findable/cos(theta) for base MC selection").fill(cth);
+ aida.cloud1D("findable/eta for base MC selection").fill(eta);
+ aida.cloud1D("findable/d0 for base MC selection").fill(d0);
+ aida.cloud1D("findable/z0 for base MC selection").fill(z0);
+ aida.cloud2D("findable/Hits vs eta for base MC selection").fill(eta, nhits);
+ /*
+ // Make plots for findable tracks
+ if (findable.isFindable(mcp, slist)) {
+ aida.cloud1D("findable/Hits for findable tracks").fill(nhits);
+ aida.cloud1D("findable/pT for findable tracks").fill(pt);
+ aida.cloud1D("findable/cos(theta) for findable tracks").fill(cth);
+ aida.cloud1D("findable/eta for findable tracks").fill(eta);
+ aida.cloud1D("findable/d0 for findable tracks").fill(d0);
+ aida.cloud1D("findable/z0 for findable tracks").fill(z0);
+ aida.cloud2D("findable/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.cloud1D("findable/Hits for z0 check failures").fill(nhits);
+ aida.cloud1D("findable/pT for z0 check failures").fill(pt);
+ aida.cloud1D("findable/cos(theta) for z0 check failures").fill(cth);
+ aida.cloud1D("findable/eta for z0 check failures").fill(eta);
+ aida.cloud1D("findable/d0 for z0 check failures").fill(d0);
+ aida.cloud1D("findable/z0 for z0 check failures").fill(z0);
+ continue;
+ }
+
+ // Select mc particles that fail on the d0 cut
+ ignores.add(Ignore.NoDCACut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.cloud1D("findable/Hits for d0 check failures").fill(nhits);
+ aida.cloud1D("findable/pT for d0 check failures").fill(pt);
+ aida.cloud1D("findable/cos(theta) for d0 check failures").fill(cth);
+ aida.cloud1D("findable/eta for d0 check failures").fill(eta);
+ aida.cloud1D("findable/d0 for d0 check failures").fill(d0);
+ aida.cloud1D("findable/z0 for d0 check failures").fill(z0);
+ continue;
+ }
+
+ // select mc particles that fail the confirm check
+ ignores.add(Ignore.NoConfirmCheck);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.cloud1D("findable/Hits for confirm check failures").fill(nhits);
+ aida.cloud1D("findable/pT for confir check failures").fill(pt);
+ aida.cloud1D("findable/cos(theta) for confirm check failures").fill(cth);
+ aida.cloud1D("findable/eta for confirm check failures").fill(eta);
+ aida.cloud1D("findable/d0 for seed confirm failures").fill(d0);
+ aida.cloud1D("findable/z0 for seed confirm failures").fill(z0);
+ continue;
+ }
+
+ // select mc particles that fail on the seed check
+ ignores.add(Ignore.NoSeedCheck);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.cloud1D("findable/Hits for seed check failures").fill(nhits);
+ aida.cloud1D("findable/pT for seed check failures").fill(pt);
+ aida.cloud1D("findable/cos(theta) for seed check failures").fill(cth);
+ aida.cloud1D("findable/eta for seed check failures").fill(eta);
+ aida.cloud1D("findable/d0 for seed check failures").fill(d0);
+ aida.cloud1D("findable/z0 for seed check failures").fill(z0);
+ continue;
+ }
+
+ // Select mc particles that fail the number of hit cut
+ ignores.add(Ignore.NoMinHitCut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.cloud1D("findable/Hits for nhit check failures").fill(nhits);
+ aida.cloud1D("findable/pT for nhit check failures").fill(pt);
+ aida.cloud1D("findable/cos(theta) for nhit check failures").fill(cth);
+ aida.cloud1D("findable/eta for nhit check failures").fill(eta);
+ aida.cloud1D("findable/d0 for nhit check failures").fill(d0);
+ aida.cloud1D("findable/z0 for nhit check failures").fill(z0);
+ continue;
+ }
+
+ // Select mc particles that fail on the pT cut
+ ignores.add(Ignore.NoPTCut);
+ if (findable.isFindable(mcp, slist, ignores)) {
+ aida.cloud1D("findable/Hits for pT check failures").fill(nhits);
+ aida.cloud1D("findable/pT for pT check failures").fill(pt);
+ aida.cloud1D("findable/cos(theta) for pT check failures").fill(cth);
+ aida.cloud1D("findable/eta for pT check failures").fill(eta);
+ aida.cloud1D("findable/d0 for pT check failures").fill(d0);
+ aida.cloud1D("findable/z0 for pT check failures").fill(z0);
+ } else {
+ System.out.println("MC Particle is not findable with all ignores set!!");
+ }
+ */
+ }
+ 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;
+ }
+}
+
+
diff -N TrackReconstructionDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TrackReconstructionDriver.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,47 @@
+/*
+ * TrackReconstructionDriver class
+ */
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+import java.util.List;
+
+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;
+
+/**
+ * Driver to perform hit digitization and track reconstruction for the sATLAS detector
+ *
+ * @author M. Graham and R. Partridge
+ */
+public class TrackReconstructionDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ public TrackReconstructionDriver() {
+
+ // Digitization and hit making driver for planar sensors
+ TrackerHitDriver_sATLAS thd = new TrackerHitDriver_sATLAS();
+ add(thd);
+
+ // Driver to make HelicalTrackHits for tracking
+ HelicalTrackHitDriver hitdriver = new HelicalTrackHitDriver();
+ hitdriver.addCollection(((TrackerHitDriver_sATLAS) thd).getStripHits1DName(), HitType.Digitized);
+ hitdriver.addCollection(((TrackerHitDriver_sATLAS) thd).getPixelHitsName(), HitType.Digitized);
+ hitdriver.OutputCollection("HelicalTrackHits");
+ add(hitdriver);
+
+ // Tracking code
+// String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "sATLASFull-JeffMarch26.xml";
+// String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "sATLASFull-Simple.xml";
+ String sfile = StrategyXMLUtils.getDefaultStrategiesPrefix() + "sATLASFull-UTOPIA7.xml";
+ List<SeedStrategy> slist = StrategyXMLUtils.getStrategyListFromResource(sfile);
+
+ SeedTracker st = new SeedTracker(slist);
+ add(st);
+ }
+}
diff -N StandAloneDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ StandAloneDriver.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,28 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+import java.io.File;
+import java.io.IOException;
+import org.freehep.record.loop.LoopException;
+import org.lcsim.util.Driver;
+import org.lcsim.util.loop.LCSimLoop;
+
+/**
+ *
+ * @author richp
+ */
+public class StandAloneDriver extends Driver {
+
+ public static void main(String[] args) throws IOException, LoopException {
+ LCSimLoop loop = new LCSimLoop();
+// loop.setLCIORecordSource(new File("/nfs/sulky21/g.ec.u12/users/mgraham/AtlasUpgrade/testStereo-muons-150events.slcio"));
+// loop.setLCIORecordSource(new File("/nfs/sulky21/g.ec.u12/users/mgraham/AtlasUpgrade/mbussonn/out/outdir-With-010-Muons-025-cat/sATLAS-With-10-Muons-025-IntPerEvent-160-EventPerFile-2.slcio"));
+ loop.setLCIORecordSource(new File("/nfs/sulky21/g.ec.u12/users/mgraham/AtlasUpgrade/sATLASFull-UTOPIA7-muons-1000events.slcio"));
+ loop.add(new sATLASMainDriver());
+ loop.loop(10); // 0 means loop forever
+ loop.dispose();
+ }
+}
diff -N OccupancyDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ OccupancyDriver.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,230 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+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.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.detector.IDetectorElement;
+import org.lcsim.detector.IReadout;
+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.MCParticle;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.tracking.digitization.sisim.GenericReadoutChip;
+import org.lcsim.recon.tracking.digitization.sisim.ReadoutChip;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author partridge
+ */
+public class OccupancyDriver 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> occMap = new HashMap<String, IProfile1D>();
+ int nevt = 0;
+
+ public OccupancyDriver() {
+
+ // 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();
+ }
+
+ @Override
+ public void process(
+ EventHeader event) {
+
+ // Increment the event counter
+ nevt++;
+
+ System.out.println("B: "+event.getDetector().getFieldMap().getField(new BasicHep3Vector(0., 0.,0.)));
+
+ ReadoutChip chip = new GenericReadoutChip();
+
+ int hittot = 0;
+ for (SiSensor sensor : _process_sensors) {
+
+ SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE);
+ if(sensor.hasPixels()){
+// System.out.println("Found some pixels");
+ electrodes = sensor.getReadoutElectrodes(ChargeCarrier.ELECTRON);
+ } else{
+// System.out.println("...got strips");
+ }
+ int nchan = electrodes.getNCells();
+
+ IReadout readout = sensor.getReadout();
+
+ List<RawTrackerHit> raw_hits = readout.getHits(RawTrackerHit.class);
+ int nhits = raw_hits.size();
+ hittot += nhits;
+ Set<SimTrackerHit> simhits = new HashSet<SimTrackerHit>();
+/*
+ for (RawTrackerHit hit : raw_hits) {
+ double charge = chip.decodeCharge(hit);
+ aida.histogram1D("Hit charge", 100, 0., 100000.).fill(charge);
+ aida.cloud1D("Number of SimTrackerHits per raw hit").fill(hit.getSimTrackerHits().size());
+ double energy = 0.;
+ for (SimTrackerHit shit : hit.getSimTrackerHits()) {
+ simhits.add(shit);
+ energy += 1000 * shit.getdEdx();
+ }
+ aida.histogram1D("SimTrackerHit energy", 100, 0., 10.).fill(energy);
+// aida.cloud2D("charge vs SimTrackerHit energy").fill(charge, energy);
+ }
+*/
+
+ double occ = 4. * ((double) nhits) / ((double) nchan);
+// double occ = ((double) nhits) / ((double) nchan);
+// if (nhits > 0) {
+// System.out.println("sensor "+sensor.getName()+" has nhits: " + nhits + " nchan: " + nchan + " occ: " + occ);
+// }
+ 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");
+ if (!barrel) continue;
+ String identifier;
+ String identifier2;
+ if (barrel) {
+ identifier = "Barrel at r = " + ri;
+ identifier2= "Barrel vs r";
+ } else {
+ identifier = "Endcap at z = " + zi;
+ identifier2 = "Endcap vs z";
+ }
+
+ if (!occMap.containsKey(identifier)) {
+ occMap.put(identifier, _hf.createProfile1D(identifier, 50, 0., 120.));
+ }
+ if (!occMap.containsKey(identifier2)) {
+ occMap.put(identifier2, _hf.createProfile1D(identifier2, 50, 0., 120.));
+ }
+ if (barrel) {
+ occMap.get(identifier).fill(zi, occ);
+ occMap.get(identifier2).fill(ri, occ);
+ } else {
+ occMap.get(identifier).fill(ri, occ);
+ occMap.get(identifier2).fill(zi, occ);
+ }
+ }
+
+ System.out.println("Total number of hit channels: " + hittot);
+/*
+ List<List<SimTrackerHit>> simall =
+ (List<List<SimTrackerHit>>) event.get(SimTrackerHit.class);
+
+ for (List<SimTrackerHit> simcol : simall) {
+ for (SimTrackerHit hit : simcol) {
+ IDetectorElement de = hit.getDetectorElement();
+ Hep3Vector pos = hit.getPositionVec();
+ double x = pos.x();
+ double y = pos.y();
+ double z = pos.z();
+ double r = Math.sqrt(x * x + y * y);
+ Hep3Vector p = hit.getMCParticle().getMomentum();
+ double px = p.x();
+ double py = p.y();
+ double pz = p.z();
+ double pr = Math.sqrt(px * px + py * py);
+ double eta = -Math.log(Math.tan(0.5 * Math.atan2(pr, pz)));
+ String detname = hit.getSubdetector().getName();
+ int layer = hit.getLayer();
+ aida.histogram1D("MC eta for Layer " + detname + layer, 50, -2.5, 2.5).fill(eta);
+ aida.histogram1D("SimTracker Hi z for Layer " + detname + layer, 50, -120., 120.).fill(z / 10);
+// de.
+ }
+ }
+
+ for (MCParticle mcp : event.getMCParticles()) {
+ if (mcp.getCharge() == 0) continue;
+ Hep3Vector p = mcp.getMomentum();
+ double px = p.x();
+ double py = p.y();
+ double pz = p.z();
+ double pr = Math.sqrt(px * px + py * py);
+ double eta = -Math.log(Math.tan(0.5 * Math.atan2(pr, pz)));
+ if (mcp.getGeneratorStatus() == mcp.FINAL_STATE) {
+ aida.histogram1D("MCP eta - final state",50,-2.5, 2.5).fill(eta);
+ }
+ if (mcp.getSimulatorStatus().isCreatedInSimulation()) {
+ aida.histogram1D("MCP eta - created in simulation",50, -2.5, 2.5).fill(eta);
+ }
+ }
+*/
+ return;
+ }
+
+ @Override
+ public void endOfData() {
+ try {
+ aida.saveAs(outputPlots);
+ } catch (IOException ex) {
+ Logger.getLogger(TrackAnalysisDriver.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ }
+
+ public void setOutputPlots(String output) {
+ this.outputPlots = output;
+
+ }
+
+ public void detectorChanged(Detector detector) {
+ System.out.println("In Occupancy Driver : "+detector.getName());
+ super.detectorChanged(detector);
+
+ // Process detectors specified by path, otherwise process entire detector
+ IDetectorElement detector_de = detector.getDetectorElement();
+ System.out.println("In Occupancy Driver : detector_de = "+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));
+ }
+
+ }
+}
\ No newline at end of file
diff -N SiVertexBarrelSensorSetup.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ SiVertexBarrelSensorSetup.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,100 @@
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+import hep.physics.matrix.BasicMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.util.List;
+
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IRotation3D;
+import org.lcsim.detector.ITranslation3D;
+import org.lcsim.detector.RotationPassiveXYZ;
+import org.lcsim.detector.Transform3D;
+import org.lcsim.detector.Translation3D;
+import org.lcsim.detector.solids.IPolyhedron;
+import org.lcsim.detector.solids.Polygon3D;
+import org.lcsim.detector.tracker.silicon.ChargeCarrier;
+import org.lcsim.detector.tracker.silicon.SiPixels;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.compact.Subdetector;
+import org.lcsim.geometry.subdetector.SiTrackerBarrel;
+import org.lcsim.util.Driver;
+
+public class SiVertexBarrelSensorSetup extends Driver
+{
+ String subdetectorName;
+
+ public SiVertexBarrelSensorSetup()
+ {}
+
+ public SiVertexBarrelSensorSetup(String subdetectorName)
+ {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void setSubdetectorName(String subdetectorName)
+ {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void detectorChanged(Detector detector)
+ {
+ if (subdetectorName == null)
+ throw new RuntimeException("The subdetectorName was not set.");
+
+ Subdetector subdetector = detector.getSubdetector(subdetectorName);
+ if (subdetector instanceof SiTrackerBarrel)
+ setupSensorDetectorElements(subdetector);
+ else
+ throw new RuntimeException("The subdetector " + subdetectorName + " is not an instance of SiTrackerBarrel.");
+ }
+
+ private void setupSensorDetectorElements(Subdetector subdet)
+ {
+ System.out.println(this.getClass().getCanonicalName() + " - Setting up sensors for " + subdet.getName() + " ...");
+
+ for ( IDetectorElement layer : subdet.getDetectorElement().getChildren() )
+ {
+ for ( IDetectorElement module : layer.getChildren() )
+ {
+ List<SiSensor> sensors = module.findDescendants(SiSensor.class);
+
+ if (sensors.size() == 0)
+ throw new RuntimeException("No sensors found in module.");
+
+ SiSensor sensor = sensors.get(0);
+ IPolyhedron sensor_solid = (IPolyhedron) sensor.getGeometry().getLogicalVolume().getSolid();
+
+ Polygon3D top_side = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, 0, 1)).get(0);
+ Polygon3D bot_side = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, 0, -1)).get(0);
+
+ // collect electrons on the top side
+ sensor.setBiasSurface(ChargeCarrier.HOLE, bot_side);
+ sensor.setBiasSurface(ChargeCarrier.ELECTRON, top_side);
+
+ // Add sense and readout electrodes
+ ITranslation3D electrodes_position = new Translation3D(VecOp.mult(-top_side.getDistance(), top_side.getNormal())); // translate to p_side
+ IRotation3D electrodes_rotation = new RotationPassiveXYZ(0.0, 0.0, 0.0);
+ // no rotation (global x-y = local x-y for axial strips)
+
+ Transform3D electrodes_transform = new Transform3D(electrodes_position, electrodes_rotation);
+
+ // Define the pixel electrodes...collecting holes;
+ SiSensorElectrodes readout_electrodes = new SiPixels(ChargeCarrier.ELECTRON, 0.05, 0.25, sensor, electrodes_transform);
+ SiSensorElectrodes sense_electrodes = new SiPixels(ChargeCarrier.ELECTRON, 0.05, 0.25, sensor, electrodes_transform);
+
+ // Tell the sensor about the electrodes
+ sensor.setSenseElectrodes(sense_electrodes);
+ sensor.setReadoutElectrodes(readout_electrodes);
+
+ // Define the transfer efficiency from sense electrodes to readout electrodes
+ // For pixels, we do a direct transfer of charge from sense electrodes to readout electrodes
+ double[][] transfer_efficiencies = {{1.0}};
+ sensor.setTransferEfficiencies(ChargeCarrier.ELECTRON, new BasicMatrix(transfer_efficiencies));
+ }
+ }
+ }
+}
diff -N sATLASMainDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ sATLASMainDriver.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,40 @@
+/*
+ * sATLASMainDriver class
+ */
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+import org.lcsim.util.Driver;
+import org.lcsim.util.loop.LCIODriver;
+
+/**
+ * Driver for track reconstruction and analysis of sATLAS detector
+ *
+ * @author M. Graham and R. Partridge
+ */
+public class sATLASMainDriver extends Driver {
+
+ public String outputFile = "foobar.slcio";
+ public String plotsFile = "myplots.aida";
+ TrackAnalysisDriver tad;
+// OccupancyDriver tad;
+
+ public sATLASMainDriver() {
+ add(new MakeSensorsDriver());
+ add(new TrackReconstructionDriver());
+ tad = new TrackAnalysisDriver();
+ add(tad);
+//
+// tad=new OccupancyDriver();
+// add(tad);
+ }
+
+ public void setOutputFile(String outputFile) {
+ System.out.println("Will output events to " + outputFile);
+ add(new LCIODriver(outputFile));
+ }
+
+ public void setPlotsFile(String plotsFile) {
+ System.out.println("Will output plots to " + plotsFile);
+ tad.setOutputPlots(plotsFile);
+ }
+}
diff -N SiTrackerBarrelSensorSetup.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ SiTrackerBarrelSensorSetup.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,100 @@
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+import hep.physics.matrix.BasicMatrix;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+
+import java.util.List;
+
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IPhysicalVolume;
+import org.lcsim.detector.IRotation3D;
+import org.lcsim.detector.ITranslation3D;
+import org.lcsim.detector.RotationPassiveXYZ;
+import org.lcsim.detector.Transform3D;
+import org.lcsim.detector.Translation3D;
+import org.lcsim.detector.solids.IPolyhedron;
+import org.lcsim.detector.solids.Polygon3D;
+import org.lcsim.detector.tracker.silicon.ChargeCarrier;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
+import org.lcsim.detector.tracker.silicon.SiStrips;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.compact.Subdetector;
+import org.lcsim.geometry.subdetector.SiTrackerBarrel;
+import org.lcsim.util.Driver;
+
+public class SiTrackerBarrelSensorSetup extends Driver {
+
+ String subdetectorName;
+
+ public SiTrackerBarrelSensorSetup() {
+ }
+
+ public SiTrackerBarrelSensorSetup(String subdetectorName) {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void setSubdetectorName(String subdetectorName) {
+ this.subdetectorName = subdetectorName;
+ }
+
+ public void detectorChanged(Detector detector) {
+ if (subdetectorName == null) {
+ throw new RuntimeException("The subdetectorName was not set.");
+ }
+
+ Subdetector subdetector = detector.getSubdetector(subdetectorName);
+ if (subdetector instanceof SiTrackerBarrel) {
+ setupSensorDetectorElements(subdetector);
+ } else {
+ throw new RuntimeException("The subdetector " + subdetectorName + " is not an instance of SiTrackerBarrel.");
+ }
+ }
+
+ private void setupSensorDetectorElements(Subdetector subdet) {
+ System.out.println(this.getClass().getCanonicalName() + " - Setting up sensors for " + subdet.getName() + " ...");
+
+ for (IDetectorElement layer : subdet.getDetectorElement().getChildren()) {
+ for (IDetectorElement module : layer.getChildren()) {
+ List<SiSensor> sensors = module.findDescendants(SiSensor.class);
+
+ if (sensors.size() == 0) {
+ throw new RuntimeException("No sensors found in module.");
+ }
+
+
+ for (SiSensor sensor : sensors) {
+ // Set up SiStrips for the sensors
+ IPolyhedron sensor_solid = (IPolyhedron) sensor.getGeometry().getLogicalVolume().getSolid();
+
+ // Bias the sensor
+ Polygon3D p_side = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, 0, 1)).get(0);
+
+ Polygon3D n_side = sensor_solid.getFacesNormalTo(new BasicHep3Vector(0, 0, -1)).get(0);
+
+ sensor.setBiasSurface(ChargeCarrier.HOLE, p_side);
+ sensor.setBiasSurface(ChargeCarrier.ELECTRON, n_side);
+
+ // Add sense and readout electrodes
+ ITranslation3D electrodes_position = new Translation3D(VecOp.mult(-p_side.getDistance(), p_side.getNormal())); // translate to p_side
+ IRotation3D electrodes_rotation = new RotationPassiveXYZ(0.0, 0.0, 0.0); // no rotation (global x-y = local x-y for axial strips)
+ Transform3D electrodes_transform = new Transform3D(electrodes_position, electrodes_rotation);
+
+ // Free calculation of readout electrodes, sense electrodes determined thereon
+ SiSensorElectrodes readout_electrodes = new SiStrips(ChargeCarrier.HOLE, 0.075, sensor, electrodes_transform);
+// SiSensorElectrodes sense_electrodes = new SiStrips(ChargeCarrier.HOLE,0.025,(readout_electrodes.getNCells()*2-1),sensor,electrodes_transform);
+ SiSensorElectrodes sense_electrodes = new SiStrips(ChargeCarrier.HOLE, 0.075, sensor, electrodes_transform);
+
+ sensor.setSenseElectrodes(sense_electrodes);
+ sensor.setReadoutElectrodes(readout_electrodes);
+ double[][] transfer_efficiencies = {{1.0}};
+
+// double[][] transfer_efficiencies = { {0.986,0.419} };
+ sensor.setTransferEfficiencies(ChargeCarrier.HOLE, new BasicMatrix(transfer_efficiencies));
+// System.out.println("Set up the sensor!"+sensor.getName());
+ }
+ }
+ }
+ }
+}
diff -N TrackerHitDriver_sATLAS.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ TrackerHitDriver_sATLAS.java 20 Aug 2009 05:06:59 -0000 1.1 @@ -0,0 +1,243 @@
+/*
+ * TrackerHitDriver Class
+ *
+ */
+package org.lcsim.contrib.sATLAS.UTOPIA7;
+
+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.ClusteringAlgorithm;
+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 _strip_digitizer;
+ SiDigitizer _pixel_digitizer;
+ StripHitMaker _strip_clusterer;
+ PixelHitMaker _pixel_clusterer;
+ String _digitizer_name;
+ 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.setNoiseIntercept(1000.);
+ // strip_readout.setNoiseSlope(100.);
+ strip_readout.setNoiseThreshold(9000.);
+ strip_readout.setNeighborThreshold(9000.);
+
+ GenericReadoutChip pixel_readout = new GenericReadoutChip();
+ pixel_readout.setNoiseIntercept(300.);
+ pixel_readout.setNoiseSlope(30.);
+ // Instantiate the digitizer that produces the raw hits
+ _strip_digitizer = new RawTrackerHitMaker(strip_simulation, strip_readout);
+ _pixel_digitizer = new RawTrackerHitMaker(pixel_simulation, pixel_readout);
+ _digitizer_name = _strip_digitizer.getName();
+// pixel_readout.setNoiseIntercept(1000.);
+// pixel_readout.setNoiseSlope(100.);
+ pixel_readout.setNoiseThreshold(7000.);
+ pixel_readout.setNeighborThreshold(7000.);
+
+ // Instantiate the digitizer that produces the raw hits
+ _strip_digitizer = new RawTrackerHitMaker(strip_simulation, strip_readout);
+ _pixel_digitizer = new RawTrackerHitMaker(pixel_simulation, pixel_readout);
+ _digitizer_name = _strip_digitizer.getName();
+
+ // Instantiate a nearest neighbor clustering algorithm for the pixels
+ NearestNeighbor strip_clustering = new NearestNeighbor();
+ strip_clustering.setSeedThreshold(9000.);
+ strip_clustering.setNeighborThreshold(9000.);
+
+ // Instantiate a nearest neighbor clustering algorithm for the pixels
+ NearestNeighbor pixel_clustering = new NearestNeighbor();
+ pixel_clustering.setSeedThreshold(7000.);
+ pixel_clustering.setNeighborThreshold(7000.);
+
+ // 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);
+ _strip_clusterer.SetOneClusterErr(1 / Math.sqrt(12.));
+ _strip_clusterer.SetTwoClusterErr(1 / 5.0);
+ _strip_clusterer.SetThreeClusterErr(1 / 3.0);
+ _strip_clusterer.SetFourClusterErr(1 / 2.0);
+ _strip_clusterer.SetFiveClusterErr(1 / 1.0);
+
+ _pixel_clusterer = new PixelHitMaker(pixel_simulation, pixel_readout, pixel_clustering);
+ _pixel_clusterer.SetOneClusterErr(1 / Math.sqrt(12.));
+ _pixel_clusterer.SetTwoClusterErr(1 / 5.0);
+ _pixel_clusterer.SetThreeClusterErr(1 / 3.0);
+ _pixel_clusterer.SetFourClusterErr(1 / 2.0);
+ _pixel_clusterer.SetFiveClusterErr(1 / 1.0);
+
+
+
+ // Specify the readouts to process
+ _readouts.add("VtxBarrHits");
+ _readouts.add("SCTShortBarrHits");
+ _readouts.add("SCTLongBarrHits");
+ _readouts.add("VtxEndcapHits");
+// _readouts.add("SCTShortEndcapHits");
+// _readouts.add("SCTLongEndcapHits");
+ _readouts.add("SCTEndcapHits");
+ // 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");
+ _process_paths.add("SCTEndcap");
+ }
+
+ /**
+ * 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 Name ="+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) {
+
+ if (sensor.hasStrips()) {
+ raw_hits.addAll(_strip_digitizer.makeHits(sensor));
+ hits_strip1D.addAll(_strip_clusterer.makeHits(sensor));
+
+ }
+
+
+ if (sensor.hasPixels()) {
+ raw_hits.addAll(_pixel_digitizer.makeHits(sensor));
+ 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_name + "_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