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