7 added files
java/trunk/analysis/src/main/java/org/hps/analysis/dataquality
--- java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/DataQualityMonitor.java (rev 0)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/DataQualityMonitor.java 2014-04-30 02:40:51 UTC (rev 533)
@@ -0,0 +1,89 @@
+package org.hps.analysis.dataquality;
+
+import java.sql.ResultSet;
+import java.sql.SQLException;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ * sort of an interface for DQM analysis drivers
+ * creates the DQM database manager, checks whether row exists in db etc
+ * @author mgraham on Apr 15, 2014
+ */
+public class DataQualityMonitor extends Driver {
+
+ public AIDA aida = AIDA.defaultInstance();
+ public DQMDatabaseManager manager;
+ public String recoVersion = "v0.0";
+ public static int runNumber = 1350;
+ public boolean overwriteDB = false;
+
+ public void setRecoVersion(String recoVersion) {
+ this.recoVersion = recoVersion;
+ }
+
+ public void setOverwriteDB(boolean overwrite) {
+ this.overwriteDB = overwrite;
+ }
+
+ public void DataQualityMonitor() {
+
+ }
+
+ public void endOfData() {
+ manager = DQMDatabaseManager.getInstance();
+ //fill any plots that only get filled at end of data...e.g. SVT occupancy plots
+ fillEndOfRunPlots();
+
+ //check to see if I need to make a new db entry
+ boolean entryExists = false;
+ try {
+ entryExists = checkRowExists();
+ if(entryExists)System.out.println("Found an existing run/reco entry in the dqm database");
+ } catch (SQLException ex) {
+ Logger.getLogger(DataQualityMonitor.class.getName()).log(Level.SEVERE, null, ex);
+ }
+
+ if (!entryExists)
+ makeNewRow();
+ else
+ if (!overwriteDB)
+ return; //entry exists and I don't want to overwrite
+
+ dumpDQMData();
+ }
+
+ private void makeNewRow() {
+ System.out.println("is the data base connected? " + manager.isConnected);
+ if (manager.isConnected) {
+ String ins = "insert into dqm SET run=" + runNumber;
+// System.out.println(ins);
+ manager.updateQuery(ins);
+ ins = "update dqm SET recoversion='" + recoVersion + "' where run=" + runNumber;
+ manager.updateQuery(ins);
+ }
+ System.out.println("Made a new row for run=" + runNumber + "; recon version=" + recoVersion);
+ }
+
+ private boolean checkRowExists() throws SQLException {
+ String ins = "select * from dqm where " + getRunRecoString();
+ ResultSet res = manager.selectQuery(ins);
+ if (res.next()) //this is a funny way of determining if the ResultSet has any entries
+ return true;
+ return false;
+ }
+
+ public String getRunRecoString() {
+ return "run=" + runNumber + " and recoversion='"+recoVersion+"'";
+ }
+
+ //override this method to do something interesting
+ public void fillEndOfRunPlots() {
+ }
+ //override this method to do something interesting
+ public void dumpDQMData() {
+ }
+
+}
java/trunk/analysis/src/main/java/org/hps/analysis/dataquality
--- java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/EcalMonitoring.java (rev 0)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/EcalMonitoring.java 2014-04-30 02:40:51 UTC (rev 533)
@@ -0,0 +1,11 @@
+package org.hps.analysis.dataquality;
+
+
+/**
+ *
+ * @author mgraham on Mar 28, 2014
+ */
+
+public class EcalMonitoring {
+
+}
java/trunk/analysis/src/main/java/org/hps/analysis/dataquality
--- java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/HitMCEfficiency.java (rev 0)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/HitMCEfficiency.java 2014-04-30 02:40:51 UTC (rev 533)
@@ -0,0 +1,163 @@
+package org.hps.analysis.dataquality;
+
+import hep.aida.IHistogramFactory;
+import hep.aida.IProfile1D;
+import java.util.List;
+import java.util.Set;
+import org.hps.recon.tracking.FittedRawTrackerHit;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+
+/**
+ * DQM driver for the monte carlo SVT hit efficiency
+ * April 29 -- first pass, makes the SimTrackerHits-->SiClusters efficiency vs position (with a settable t0 cut)
+ * @author mgraham on April 29, 2014
+ */
+// TODO: Add HelicalTrackHit efficiency...this should include the fitted track hit cuts (t0 & chi^2) automatically since that where the cut is applied
+// TODO: Add some quantities for DQM monitoring: e.g. <efficiency>, probably within first 1 cm or so.
+public class HitMCEfficiency extends DataQualityMonitor {
+
+ private String rawTrackerHitCollectionName = "SVTRawTrackerHits";
+ private String helicalTrackHitCollectionName = "HelicalTrackHits";
+ private String rotatedTrackHitCollectionName = "RotatedHelicalTrackHits";
+ private String fittedTrackerHitCollectionName = "SVTFittedRawTrackerHits";
+ private String trackerHitCollectionName = "TrackerHits";
+ private String siClusterCollectionName = "StripClusterer_SiTrackerHitStrip1D";
+ private String svtTrueHitRelationName = "SVTTrueHitRelations";
+ private String trackerName = "Tracker";
+ private Detector detector = null;
+ private double t0Cut=16.0;
+ private static final String nameStrip = "Tracker_TestRunModule_";
+ private List<SiSensor> sensors;
+
+ public void setHelicalTrackHitCollectionName(String helicalTrackHitCollectionName) {
+ this.helicalTrackHitCollectionName = helicalTrackHitCollectionName;
+ }
+
+ public void setT0Cut(double cut){
+ this.t0Cut=cut;
+ }
+
+ @Override
+ protected void detectorChanged(Detector detector) {
+ this.detector = detector;
+ aida.tree().cd("/");
+ IHistogramFactory hf = aida.histogramFactory();
+
+
+ // Make a list of SiSensors in the SVT.
+ sensors = this.detector.getSubdetector(trackerName).getDetectorElement().findDescendants(SiSensor.class);
+
+ // Setup the efficiency plots.
+ //currently, just the Si cluster efficiency
+ aida.tree().cd("/");
+ for (int kk = 1; kk < 13; kk++) {
+ IProfile1D clEffic = createLayerPlot("clusterEfficiency", kk, 50, 0, 25.);
+ }
+ }
+
+ @Override
+ public void process(EventHeader event) {
+
+ aida.tree().cd("/");
+
+ //make sure the required collections exist
+ if (!event.hasCollection(RawTrackerHit.class, rawTrackerHitCollectionName))
+ return;
+ if (!event.hasCollection(FittedRawTrackerHit.class, fittedTrackerHitCollectionName))
+ return;
+
+ if (!event.hasCollection(SiTrackerHitStrip1D.class, siClusterCollectionName))
+ return;
+
+ if (!event.hasCollection(SimTrackerHit.class, trackerHitCollectionName))
+ return;
+
+ if (!event.hasCollection(LCRelation.class, svtTrueHitRelationName))
+ return;
+
+ RelationalTable mcHittomcP = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ // Get the collections of SimTrackerHits
+ List<List<SimTrackerHit>> simcols = event.get(SimTrackerHit.class);
+ // Loop over the SimTrackerHits and fill in the relational table
+ for (List<SimTrackerHit> simlist : simcols) {
+ for (SimTrackerHit simhit : simlist) {
+ if (simhit.getMCParticle() != null)
+ mcHittomcP.add(simhit, simhit.getMCParticle());
+ }
+ }
+ RelationalTable rawtomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ if (event.hasCollection(LCRelation.class,svtTrueHitRelationName)) {
+ List<LCRelation> trueHitRelations = event.get(LCRelation.class,svtTrueHitRelationName);
+ for (LCRelation relation : trueHitRelations) {
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+ rawtomc.add(relation.getFrom(), relation.getTo());
+ }
+ }
+ List<SimTrackerHit> simHits = event.get(SimTrackerHit.class, trackerHitCollectionName);
+ // make relational table for strip clusters to mc particle
+ List<SiTrackerHitStrip1D> siClusters = event.get(SiTrackerHitStrip1D.class, siClusterCollectionName);
+ RelationalTable clustertosimhit = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ for (SiTrackerHit cluster : siClusters) {
+ List<RawTrackerHit> rawHits = cluster.getRawHits();
+ for (RawTrackerHit rth : rawHits) {
+ Set<SimTrackerHit> simTrackerHits = rawtomc.allFrom(rth);
+ if (simTrackerHits != null)
+ for (SimTrackerHit simhit : simTrackerHits) {
+ clustertosimhit.add(cluster, simhit);
+ }
+ }
+ }
+//relational tables from mc particle to raw and fitted tracker hits
+ RelationalTable fittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ List<FittedRawTrackerHit> fittedTrackerHits = event.get(FittedRawTrackerHit.class, fittedTrackerHitCollectionName);
+ for (FittedRawTrackerHit hit : fittedTrackerHits) {
+ RawTrackerHit rth = hit.getRawTrackerHit();
+ Set<SimTrackerHit> simTrackerHits = rawtomc.allFrom(rth);
+ if (simTrackerHits != null)
+ for (SimTrackerHit simhit : simTrackerHits) {
+ if (simhit.getMCParticle() != null)
+ fittomc.add(hit, simhit.getMCParticle());
+ }
+ }
+
+ for (SimTrackerHit simhit : simHits) {
+ double wgt = 0.0;
+ Set<SiTrackerHitStrip1D> clusters = clustertosimhit.allTo(simhit);
+ if (clusters != null) {
+ for (SiTrackerHitStrip1D clust : clusters) {
+ if (Math.abs(clust.getTime()) < t0Cut)
+ wgt = 1.0;
+ }
+ }
+ getLayerPlot("clusterEfficiency", simhit.getLayer()).fill(Math.abs(simhit.getPoint()[1]), wgt);
+ }
+ }
+
+ @Override
+ public void fillEndOfRunPlots() {
+ }
+
+ @Override
+ public void dumpDQMData() {
+ }
+
+ private IProfile1D getLayerPlot(String prefix, int layer) {
+ return aida.profile1D(prefix + "_layer" + layer);
+ }
+
+ private IProfile1D createLayerPlot(String prefix, int layer, int nchan, double min, double max) {
+ IProfile1D hist = aida.profile1D(prefix + "_layer" + layer, nchan, min, max);
+ return hist;
+ }
+
+
+}
java/trunk/analysis/src/main/java/org/hps/analysis/dataquality
--- java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/ReconMonitoring.java (rev 0)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/ReconMonitoring.java 2014-04-30 02:40:51 UTC (rev 533)
@@ -0,0 +1,11 @@
+package org.hps.analysis.dataquality;
+
+
+/**
+ *
+ * @author mgraham on Mar 28, 2014
+ */
+
+public class ReconMonitoring {
+
+}
java/trunk/analysis/src/main/java/org/hps/analysis/dataquality
--- java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/SvtMonitoring.java (rev 0)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/SvtMonitoring.java 2014-04-30 02:40:51 UTC (rev 533)
@@ -0,0 +1,187 @@
+package org.hps.analysis.dataquality;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IFitFactory;
+import hep.aida.IHistogram1D;
+import hep.aida.IPlotter;
+import hep.aida.IPlotterStyle;
+import org.lcsim.geometry.Detector;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import org.hps.recon.tracking.FittedRawTrackerHit;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.RawTrackerHit;
+
+/**
+ * DQM driver for the monte carlo for reconstructed track quantities
+ * plots things like occupancy, t0, amplitude, chi^2 (from APV25 sampling fit); each on a per/sensor basis
+ * saves to DQM database: <occupancy>
+ * @author mgraham on Mar 28, 2014
+ */
+//TODO: add some more quantities to DQM database: <t0> or <sigma>_t0 for intime events; <chi^2>, <amplitude> etc
+public class SvtMonitoring extends DataQualityMonitor {
+
+ private String rawTrackerHitCollectionName = "SVTRawTrackerHits";
+ private String fittedTrackerHitCollectionName = "SVTFittedRawTrackerHits";
+ private String trackerHitCollectionName = "StripClusterer_SiTrackerHitStrip1D";
+ private Detector detector = null;
+ private IPlotter plotter;
+ private String trackerName = "Tracker";
+ private List<SiSensor> sensors;
+ private Map<String, int[]> occupancyMap;
+ private Map<String, Double> avgOccupancyMap;
+ private Map<String, String> avgOccupancyNames;
+ private int eventCountRaw = 0;
+ private int eventCountFit = 0;
+ private int eventCountCluster = 0;
+ private static final String nameStrip = "Tracker_TestRunModule_";
+ private static final int maxChannels = 640;
+
+ public void setRawTrackerHitCollectionName(String inputCollection) {
+ this.rawTrackerHitCollectionName = inputCollection;
+ }
+
+ public void setFittedTrackerHitCollectionName(String inputCollection) {
+ this.fittedTrackerHitCollectionName = inputCollection;
+ }
+
+ public void setTrackerHitCollectionName(String inputCollection) {
+ this.trackerHitCollectionName = inputCollection;
+ }
+
+ protected void detectorChanged(Detector detector) {
+ System.out.println("SvtMonitoring::detectorChanged Setting up the plotter");
+ this.detector = detector;
+ aida.tree().cd("/");
+
+ // Make a list of SiSensors in the SVT.
+ sensors = this.detector.getSubdetector(trackerName).getDetectorElement().findDescendants(SiSensor.class);
+
+ // Reset the data structure that keeps track of strip occupancies.
+ resetOccupancyMap();
+
+ // Setup the occupancy plots.
+ aida.tree().cd("/");
+ for (SiSensor sensor : sensors) {
+ //IHistogram1D occupancyPlot = aida.histogram1D(sensor.getName().replaceAll("Tracker_TestRunModule_", ""), 640, 0, 639);
+ IHistogram1D occupancyPlot = createSensorPlot("occupancy_",sensor, maxChannels, 0, maxChannels - 1);
+ IHistogram1D t0Plot = createSensorPlot("t0_",sensor,50,-50.,50.);
+ IHistogram1D amplitudePlot = createSensorPlot("amplitude_",sensor,50,0,2000);
+ IHistogram1D chi2Plot = createSensorPlot("chi2_",sensor,50,0,25);
+ occupancyPlot.reset();
+ }
+
+ }
+
+ public SvtMonitoring() {
+ }
+
+ public void process(EventHeader event) {
+ /* increment the strip occupancy arrays */
+ if (event.hasCollection(RawTrackerHit.class, rawTrackerHitCollectionName)) {
+ List<RawTrackerHit> rawTrackerHits = event.get(RawTrackerHit.class, rawTrackerHitCollectionName);
+ for (RawTrackerHit hit : rawTrackerHits) {
+ int[] strips = occupancyMap.get(hit.getDetectorElement().getName());
+ strips[hit.getIdentifierFieldValue("strip")] += 1;
+ }
+ ++eventCountRaw;
+ } else
+ return; /* kick out of this if the even has none of these...*/
+ /* fill the FittedTrackerHit related histograms */
+ if (event.hasCollection(FittedRawTrackerHit.class, fittedTrackerHitCollectionName)){
+ List<FittedRawTrackerHit> fittedTrackerHits = event.get(FittedRawTrackerHit.class, fittedTrackerHitCollectionName);
+ for(FittedRawTrackerHit hit: fittedTrackerHits){
+ String sensorName= hit.getRawTrackerHit().getDetectorElement().getName();
+ double t0=hit.getT0();
+ double amp=hit.getAmp();
+ double chi2=hit.getShapeFitParameters().getChiSq();
+ getSensorPlot("t0_",sensorName).fill(t0);
+ getSensorPlot("amplitude_",sensorName).fill(amp);
+ getSensorPlot("chi2_",sensorName).fill(chi2);
+ }
+ ++eventCountFit;
+ } else
+ return;
+ }
+
+ private IHistogram1D getSensorPlot( String prefix, SiSensor sensor) {
+ return aida.histogram1D(prefix+sensor.getName());
+ }
+
+ private IHistogram1D getSensorPlot( String prefix, String sensorName) {
+ return aida.histogram1D(prefix+sensorName);
+ }
+
+ private IHistogram1D createSensorPlot( String prefix,SiSensor sensor, int nchan, double min, double max) {
+ IHistogram1D hist = aida.histogram1D(prefix+sensor.getName(),nchan,min,max);
+ hist.setTitle(sensor.getName().replaceAll(nameStrip, "")
+ .replace("module", "mod")
+ .replace("layer", "lyr")
+ .replace("sensor", "sens"));
+ return hist;
+ }
+
+ private void resetOccupancyMap() {
+ occupancyMap = new HashMap<>();
+ avgOccupancyMap = new HashMap<>();
+ avgOccupancyNames = new HashMap<>();
+ for (SiSensor sensor : sensors) {
+ occupancyMap.put(sensor.getName(), new int[640]);
+ avgOccupancyMap.put(sensor.getName(), -999.);
+ String occName = "avgOcc_" + getNiceSensorName(sensor);
+ avgOccupancyNames.put(sensor.getName(), occName);
+ }
+ }
+
+ private String getNiceSensorName(SiSensor sensor) {
+ return sensor.getName().replaceAll(nameStrip, "")
+ .replace("module", "mod")
+ .replace("layer", "lyr")
+ .replace("sensor", "sens");
+ }
+
+ public void reset() {
+ eventCountRaw = 0;
+ eventCountFit = 0;
+ eventCountCluster = 0;
+ resetOccupancyMap();
+ }
+
+ @Override
+ public void fillEndOfRunPlots() {
+ // Plot strip occupancies.
+ System.out.println("SvtMonitoring::endOfData filling occupancy plots");
+ for (SiSensor sensor : sensors) {
+ Double avg = 0.0;
+ //IHistogram1D sensorHist = aida.histogram1D(sensor.getName());
+ IHistogram1D sensorHist = getSensorPlot("occupancy_",sensor);
+ sensorHist.reset();
+ int[] strips = occupancyMap.get(sensor.getName());
+ for (int i = 0; i < strips.length; i++) {
+ double stripOccupancy = (double) strips[i] / (double) (eventCountRaw);
+ if (stripOccupancy != 0)
+ sensorHist.fill(i, stripOccupancy);
+ avg += stripOccupancy;
+ }
+ avg /= strips.length;
+ avgOccupancyMap.put(sensor.getName(), avg);
+ }
+ }
+
+ @Override
+ public void dumpDQMData() {
+ System.out.println("SvtMonitoring::endOfData filling DQM database");
+ double s1occ = 0.99;
+ String put = "update dqm SET avgOcc_T1=" + s1occ + " WHERE " + getRunRecoString();
+// manager.updateQuery(put);
+ printDQMData();
+ }
+
+ public void printDQMData() {
+ for (SiSensor sensor : sensors) {
+ System.out.println(avgOccupancyNames.get(sensor.getName()) + ": " + avgOccupancyMap.get(sensor.getName()));
+ }
+ }
+}
java/trunk/analysis/src/main/java/org/hps/analysis/dataquality
--- java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/TrackMCEfficiency.java (rev 0)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/TrackMCEfficiency.java 2014-04-30 02:40:51 UTC (rev 533)
@@ -0,0 +1,362 @@
+package org.hps.analysis.dataquality;
+
+import hep.aida.IHistogramFactory;
+import hep.aida.IProfile1D;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import org.hps.analysis.examples.TrackAnalysis;
+import org.hps.recon.tracking.FindableTrack;
+import org.hps.recon.tracking.FittedRawTrackerHit;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCRelation;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.RelationalTable;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.event.base.BaseRelationalTable;
+import org.lcsim.fit.helicaltrack.HelicalTrackCross;
+import org.lcsim.fit.helicaltrack.HelixParamCalculator;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHit;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+
+/**
+ * DQM driver for the monte carlo track efficiency; makes a bunch of efficiency vs variable plots
+ * for all tracks and just electrons from trident/A' event, as well as "findable" tracks
+ * use the debugTrackEfficiency flag to print out info regarding individual failed events
+ * @author mgraham on Mar 28, 2014
+ */
+// TODO: Add some quantities for DQM monitoring: e.g. <efficiency>, <eff>_findable
+public class TrackMCEfficiency extends DataQualityMonitor {
+
+ private String rawTrackerHitCollectionName = "SVTRawTrackerHits";
+ private String helicalTrackHitCollectionName = "HelicalTrackHits";
+ private String fittedTrackerHitCollectionName = "SVTFittedRawTrackerHits";
+ private String trackerHitCollectionName = "TrackerHits";
+ private String siClusterCollectionName = "StripClusterer_SiTrackerHitStrip1D";
+ private String rotatedMCRelationsCollectionName = "RotatedHelicalTrackMCRelations";
+ private String trackCollectionName = "MatchedTracks";
+ private String trackerName = "Tracker";
+ private Detector detector = null;
+ IDDecoder dec;
+ private IProfile1D peffFindable;
+ private IProfile1D thetaeffFindable;
+ private IProfile1D phieffFindable;
+ private IProfile1D ctheffFindable;
+ private IProfile1D d0effFindable;
+ private IProfile1D z0effFindable;
+ private IProfile1D peffElectrons;
+ private IProfile1D thetaeffElectrons;
+ private IProfile1D phieffElectrons;
+ private IProfile1D ctheffElectrons;
+ private IProfile1D d0effElectrons;
+ private IProfile1D z0effElectrons;
+ double beamP = 2.2;
+ int nlayers = 12;
+ int totelectrons = 0;
+ double foundelectrons = 0;
+ int findableelectrons = 0;
+ int findableTracks = 0;
+ double foundTracks = 0;
+ private static final String nameStrip = "Tracker_TestRunModule_";
+ private List<SiSensor> sensors;
+ private boolean debugTrackEfficiency = false;
+
+ public void setHelicalTrackHitCollectionName(String helicalTrackHitCollectionName) {
+ this.helicalTrackHitCollectionName = helicalTrackHitCollectionName;
+ }
+
+ public void setTrackCollectionName(String trackCollectionName) {
+ this.trackCollectionName = trackCollectionName;
+ }
+
+ public void setDebugTrackEfficiency(boolean debug) {
+ this.debugTrackEfficiency = debug;
+ }
+
+ @Override
+ protected void detectorChanged(Detector detector) {
+ this.detector = detector;
+ aida.tree().cd("/");
+ IHistogramFactory hf = aida.histogramFactory();
+
+ peffFindable = hf.createProfile1D("Findable Efficiency vs p", "", 20, 0., beamP);
+ thetaeffFindable = hf.createProfile1D("Findable Efficiency vs theta", "", 20, 80, 100);
+ phieffFindable = hf.createProfile1D("Findable Efficiency vs phi", "", 25, -0.25, 0.25);
+ ctheffFindable = hf.createProfile1D("Findable Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
+ d0effFindable = hf.createProfile1D("Findable Efficiency vs d0", "", 50, -2., 2.);
+ z0effFindable = hf.createProfile1D("Findable Efficiency vs z0", "", 50, -2., 2.);
+
+ peffElectrons = hf.createProfile1D("Electrons Efficiency vs p", "", 20, 0., beamP);
+ thetaeffElectrons = hf.createProfile1D("Electrons Efficiency vs theta", "", 20, 80, 100);
+ phieffElectrons = hf.createProfile1D("Electrons Efficiency vs phi", "", 25, -0.25, 0.25);
+ ctheffElectrons = hf.createProfile1D("Electrons Efficiency vs cos(theta)", "", 25, -0.25, 0.25);
+ d0effElectrons = hf.createProfile1D("Electrons Efficiency vs d0", "", 20, -1., 1.);
+ z0effElectrons = hf.createProfile1D("Electrons Efficiency vs z0", "", 20, -1., 1.);
+
+ // Make a list of SiSensors in the SVT.
+ sensors = this.detector.getSubdetector(trackerName).getDetectorElement().findDescendants(SiSensor.class);
+
+ // Setup the occupancy plots.
+ aida.tree().cd("/");
+ for (int kk = 1; kk < 13; kk++) {
+ IProfile1D clEffic = createLayerPlot("clusterEfficiency", kk, 50, 0, 25.);
+ }
+ }
+
+ @Override
+ public void process(EventHeader event) {
+
+ aida.tree().cd("/");
+
+ //make sure the required collections exist
+ if (!event.hasCollection(RawTrackerHit.class, rawTrackerHitCollectionName))
+ return;
+ if (!event.hasCollection(FittedRawTrackerHit.class, fittedTrackerHitCollectionName))
+ return;
+ if (!event.hasCollection(Track.class, trackCollectionName))
+ return;
+ if (!event.hasCollection(LCRelation.class, rotatedMCRelationsCollectionName))
+ return;
+ if (!event.hasCollection(SiTrackerHitStrip1D.class, siClusterCollectionName))
+ return;
+
+ if (!event.hasCollection(SimTrackerHit.class, trackerHitCollectionName))
+ return;
+ //
+ //get the b-field
+ Hep3Vector IP = new BasicHep3Vector(0., 0., 1.);
+ double bfield = event.getDetector().getFieldMap().getField(IP).y();
+ //make some maps and relation tables
+ Map<Track, TrackAnalysis> tkanalMap = new HashMap<Track, TrackAnalysis>();
+ RelationalTable hittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ List<LCRelation> mcrelations = event.get(LCRelation.class, rotatedMCRelationsCollectionName);
+ for (LCRelation relation : mcrelations) {
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+ hittomc.add(relation.getFrom(), relation.getTo());
+ }
+ RelationalTable mcHittomcP = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ // Get the collections of SimTrackerHits
+ List<List<SimTrackerHit>> simcols = event.get(SimTrackerHit.class);
+ // Loop over the SimTrackerHits and fill in the relational table
+ for (List<SimTrackerHit> simlist : simcols) {
+ for (SimTrackerHit simhit : simlist) {
+ if (simhit.getMCParticle() != null)
+ mcHittomcP.add(simhit, simhit.getMCParticle());
+ }
+ }
+ RelationalTable trktomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ RelationalTable rawtomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ if (event.hasCollection(LCRelation.class, "SVTTrueHitRelations")) {
+ List<LCRelation> trueHitRelations = event.get(LCRelation.class, "SVTTrueHitRelations");
+ for (LCRelation relation : trueHitRelations) {
+ if (relation != null && relation.getFrom() != null && relation.getTo() != null)
+ rawtomc.add(relation.getFrom(), relation.getTo());
+ }
+ }
+ // make relational table for strip clusters to mc particle
+ List<SiTrackerHitStrip1D> siClusters = event.get(SiTrackerHitStrip1D.class, siClusterCollectionName);
+ RelationalTable clustertosimhit = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ for (SiTrackerHit cluster : siClusters) {
+ List<RawTrackerHit> rawHits = cluster.getRawHits();
+ for (RawTrackerHit rth : rawHits) {
+ Set<SimTrackerHit> simTrackerHits = rawtomc.allFrom(rth);
+ if (simTrackerHits != null)
+ for (SimTrackerHit simhit : simTrackerHits) {
+ clustertosimhit.add(cluster, simhit);
+ }
+ }
+ }
+ //relational tables from mc particle to raw and fitted tracker hits
+ RelationalTable fittomc = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
+ List<FittedRawTrackerHit> fittedTrackerHits = event.get(FittedRawTrackerHit.class, fittedTrackerHitCollectionName);
+ for (FittedRawTrackerHit hit : fittedTrackerHits) {
+ RawTrackerHit rth = hit.getRawTrackerHit();
+ Set<SimTrackerHit> simTrackerHits = rawtomc.allFrom(rth);
+ if (simTrackerHits != null)
+ for (SimTrackerHit simhit : simTrackerHits) {
+ if (simhit.getMCParticle() != null)
+ fittomc.add(hit, simhit.getMCParticle());
+ }
+ }
+
+ // Instantiate the class that determines if a track is "findable"
+ FindableTrack findable = new FindableTrack(event);
+
+ List<Track> tracks = event.get(Track.class, trackCollectionName);
+ aida.histogram1D("Tracks per Event").fill(tracks.size());
+ for (Track trk : tracks) {
+ TrackAnalysis tkanal = new TrackAnalysis(trk, hittomc);
+ tkanalMap.put(trk, tkanal);
+ MCParticle mcp = tkanal.getMCParticleNew();
+ if (mcp != null)
+ // Create a map between the tracks found and the assigned MC particle
+ trktomc.add(trk, tkanal.getMCParticleNew());
+ }
+
+ // Now loop over all MC Particles
+ List<MCParticle> mclist = event.getMCParticles();
+ int _nchMCP = 0;
+ int _nchMCPBar = 0;
+ 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));
+ double phi = Math.atan2(py, px);
+ // Find the number of layers hit by this mc particle
+// System.out.println("MC pt=" + pt);
+ int nhits = findable.LayersHit(mcp);
+ boolean isFindable = findable.InnerTrackerIsFindable(mcp, nlayers - 2);
+
+ // 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();
+
+// Set<Track> trklistAxial = trktomcAxial.allTo(mcp);
+// int ntrkAxial = trklistAxial.size();
+ if (mcp.getPDGID() == 622) {
+ boolean bothreco = true;
+ boolean bothfindable = true;
+ //it's the A'...let's see if we found both tracks.
+ List<MCParticle> daughters = mcp.getDaughters();
+ for (MCParticle d : daughters) {
+ if (trktomc.allTo(d).isEmpty())
+ bothreco = false;
+ if (!findable.InnerTrackerIsFindable(d, nlayers - 2))
+ bothfindable = false;
+ }
+ double vtxWgt = 0;
+ if (bothreco)
+ vtxWgt = 1.0;
+// VxEff.fill(mcp.getOriginX(), vtxWgt);
+// VyEff.fill(mcp.getOriginY(), vtxWgt);
+// VzEff.fill(mcp.getOriginZ(), vtxWgt);
+ if (bothfindable) {
+// VxEffFindable.fill(mcp.getOriginX(), vtxWgt);
+// VyEffFindable.fill(mcp.getOriginY(), vtxWgt);
+// VzEffFindable.fill(mcp.getOriginZ(), vtxWgt);
+ }
+ }
+
+// if (nhits == nlayers[0]) {
+ if (isFindable) {
+ _nchMCP++;
+ findableTracks++;
+ double wgt = 0.;
+ if (ntrk > 0)
+ wgt = 1.;
+ foundTracks += wgt;
+ peffFindable.fill(p, wgt);
+ phieffFindable.fill(phi, wgt);
+ thetaeffFindable.fill(theta, wgt);
+ ctheffFindable.fill(cth, wgt);
+ d0effFindable.fill(d0, wgt);
+ z0effFindable.fill(z0, wgt);
+
+ if (wgt == 0) {
+ Set<SimTrackerHit> mchitlist = mcHittomcP.allTo(mcp);
+ Set<HelicalTrackCross> hitlist = hittomc.allTo(mcp);
+ Set<FittedRawTrackerHit> fitlist = fittomc.allTo(mcp);
+ if (debugTrackEfficiency) {
+ System.out.println("TrackMCEfficiencyMonitoring:: Missed a findable track with MC p = " + p);
+ if (!hasHTHInEachLayer(hitlist, fitlist))
+ System.out.println("This track failed becasue it's missing a helical track hit");
+ }
+ }
+
+ }
+ if (mcp.getParents().size() == 1 && mcp.getParents().get(0).getPDGID() == 622) {
+ totelectrons++;
+// findableelectrons++;
+ double wgt = 0.;
+ if (ntrk > 0)
+ wgt = 1.;
+ foundelectrons += wgt;
+ peffElectrons.fill(p, wgt);
+ phieffElectrons.fill(phi, wgt);
+ thetaeffElectrons.fill(theta, wgt);
+ ctheffElectrons.fill(cth, wgt);
+ d0effElectrons.fill(d0, wgt);
+ z0effElectrons.fill(z0, wgt);
+
+ // }
+ }
+ }
+ }
+
+ @Override
+ public void fillEndOfRunPlots() {
+ }
+
+ @Override
+ public void dumpDQMData() {
+ }
+
+ private IProfile1D getLayerPlot(String prefix, int layer) {
+ return aida.profile1D(prefix + "_layer" + layer);
+ }
+
+ private IProfile1D createLayerPlot(String prefix, int layer, int nchan, double min, double max) {
+ IProfile1D hist = aida.profile1D(prefix + "_layer" + layer, nchan, min, max);
+ return hist;
+ }
+
+ private boolean hasHTHInEachLayer(Set<HelicalTrackCross> list, Set<FittedRawTrackerHit> fitlist) {
+ for (int layer = 1; layer < nlayers - 2; layer += 2) {
+ boolean hasThisLayer = false;
+ for (HelicalTrackCross hit : list) {
+ if (hit.Layer() == layer)
+ hasThisLayer = true;
+ }
+ if (!hasThisLayer) {
+ System.out.println("Missing reconstructed hit in layer = " + layer);
+ boolean hasFitHitSL1 = false;
+ boolean hasFitHitSL2 = false;
+ FittedRawTrackerHit fitSL1 = null;
+ FittedRawTrackerHit fitSL2 = null;
+ System.out.println("fitted hit list size = " + fitlist.size());
+ for (FittedRawTrackerHit fit : fitlist) {
+ System.out.println("fitted hit layer number = " + fit.getRawTrackerHit().getLayerNumber());
+ if (fit.getRawTrackerHit().getLayerNumber() == layer) {
+ hasFitHitSL1 = true;
+ fitSL1 = fit;
+ System.out.println("Found a hit in SL1 with t0 = " + fitSL1.getT0() + "; amp = " + fitSL1.getAmp() + "; chi^2 = " + fitSL1.getShapeFitParameters().getChiSq() + "; strip = " + fitSL1.getRawTrackerHit().getCellID());
+ }
+ if (fit.getRawTrackerHit().getLayerNumber() == layer + 1) {
+ hasFitHitSL2 = true;
+ fitSL2 = fit;
+ System.out.println("Found a hit in SL2 with t0 = " + fitSL2.getT0() + "; amp = " + fitSL2.getAmp() + "; chi^2 = " + fitSL2.getShapeFitParameters().getChiSq() + "; strip = " + fitSL2.getRawTrackerHit().getCellID());
+
+ }
+ }
+ if (!hasFitHitSL1)
+ System.out.println("MISSING a hit in SL1!!!");
+ if (!hasFitHitSL2)
+ System.out.println("MISSING a hit in SL2!!!");
+
+ return false;
+ }
+ }
+ return true;
+ }
+
+}
java/trunk/analysis/src/main/java/org/hps/analysis/dataquality
--- java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/TrackingMonitoring.java (rev 0)
+++ java/trunk/analysis/src/main/java/org/hps/analysis/dataquality/TrackingMonitoring.java 2014-04-30 02:40:51 UTC (rev 533)
@@ -0,0 +1,123 @@
+package org.hps.analysis.dataquality;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IProfile;
+import java.util.List;
+import org.hps.conditions.deprecated.SvtUtils;
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.LCIOParameters;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.fit.helicaltrack.HelicalTrackCross;
+import org.lcsim.fit.helicaltrack.HelicalTrackHit;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.IDDecoder;
+
+/**
+ * DQM driver for the monte carlo for reconstructed track quantities
+ * plots things like number of tracks/event, momentum, chi^2, track parameters (d0/z0/theta/phi/curvature)
+ * @author mgraham on Mar 28, 2014
+ */
+// TODO: Add some quantities for DQM monitoring: e.g. <tracks>, <hits/track>, etc
+public class TrackingMonitoring extends DataQualityMonitor {
+
+ private String helicalTrackHitCollectionName = "HelicalTrackHits";
+ private String rotatedTrackHitCollectionName = "RotatedHelicalTrackHits";
+ private String helicalTrackHitRelationsCollectionName = "HelicalTrackHitRelations";
+ private String trackCollectionName = "MatchedTracks";
+ private String trackerName = "Tracker";
+ String ecalSubdetectorName = "Ecal";
+ String ecalCollectionName = "EcalClusters";
+ private Detector detector = null;
+ IDDecoder dec;
+
+ public void setHelicalTrackHitCollectionName(String helicalTrackHitCollectionName) {
+ this.helicalTrackHitCollectionName = helicalTrackHitCollectionName;
+ }
+
+ public void setTrackCollectionName(String trackCollectionName) {
+ this.trackCollectionName = trackCollectionName;
+ }
+
+ @Override
+ protected void detectorChanged(Detector detector) {
+ this.detector = detector;
+ aida.tree().cd("/");
+
+ IProfile avgLayersTopPlot = aida.profile1D("Number of Stereo Hits per layer in Top Half", 6, 1, 13);
+ IProfile avgLayersBottomPlot = aida.profile1D("Number of Stereo Hits per layer in Bottom Half", 6, 1, 13);
+
+ IHistogram1D trkPx = aida.histogram1D("Track Momentum (Px)", 25, -0.1, 0.200);
+ IHistogram1D trkPy = aida.histogram1D("Track Momentum (Py)", 25, -0.1, 0.1);
+ IHistogram1D trkPz = aida.histogram1D("Track Momentum (Pz)", 25, 0, 2.4);
+ IHistogram1D trkChi2 = aida.histogram1D("Track Chi2", 25, 0, 25.0);
+ IHistogram1D nTracks = aida.histogram1D("Tracks per Event", 6, 0, 6);
+ IHistogram1D trkd0 = aida.histogram1D("d0 ", 25, -5.0, 5.0);
+ IHistogram1D trkphi = aida.histogram1D("sinphi ", 25, -0.2, 0.2);
+ IHistogram1D trkomega = aida.histogram1D("omega ", 25, -0.00025, 0.00025);
+ IHistogram1D trklam = aida.histogram1D("tan(lambda) ", 25, -0.1, 0.1);
+ IHistogram1D trkz0 = aida.histogram1D("z0 ", 25, -1.0, 1.0);
+ IHistogram1D nHits = aida.histogram1D("Hits per Track", 2, 5, 7);
+
+ }
+
+ @Override
+ public void process(EventHeader event) {
+
+ aida.tree().cd("/");
+ if (!event.hasCollection(HelicalTrackHit.class, helicalTrackHitCollectionName))
+ return;
+
+ List<HelicalTrackHit> hthList = event.get(HelicalTrackHit.class, helicalTrackHitCollectionName);
+ int[] layersTop = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0,0};
+ int[] layersBot = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0,0};
+ for (HelicalTrackHit hth : hthList) {
+ HelicalTrackCross htc = (HelicalTrackCross) hth;
+ double x = htc.getPosition()[0];
+ double y = htc.getPosition()[1];
+ SiSensor sensor = ((SiSensor) ((RawTrackerHit) htc.getRawHits().get(0)).getDetectorElement());
+ if (SvtUtils.getInstance().isTopLayer(sensor))
+ layersTop[htc.Layer() - 1]++;
+ else
+ layersBot[htc.Layer() - 1]++;
+ }
+ for (int i = 0; i < 12; i++) {
+ aida.profile1D("Number of Stereo Hits per layer in Top Half").fill(i + 1, layersTop[i]);
+ aida.profile1D("Number of Stereo Hits per layer in Bottom Half").fill(i + 1, layersBot[i]);
+ }
+
+ if (!event.hasCollection(Track.class, trackCollectionName)) {
+// System.out.println(trackCollectionName + " does not exist; skipping event");
+ aida.histogram1D("Tracks per Event").fill(0);
+ return;
+ }
+
+ List<Track> tracks = event.get(Track.class, trackCollectionName);
+ aida.histogram1D("Tracks per Event").fill(tracks.size());
+ for (Track trk : tracks) {
+
+ aida.histogram1D("Track Momentum (Px)").fill(trk.getPY());
+ aida.histogram1D("Track Momentum (Py)").fill(trk.getPZ());
+ aida.histogram1D("Track Momentum (Pz)").fill(trk.getPX());
+ aida.histogram1D("Track Chi2").fill(trk.getChi2());
+
+ aida.histogram1D("Hits per Track").fill(trk.getTrackerHits().size());
+
+ aida.histogram1D("d0 ").fill(trk.getTrackParameter(LCIOParameters.ParameterName.d0.ordinal()));
+ aida.histogram1D("sinphi ").fill(Math.sin(trk.getTrackParameter(LCIOParameters.ParameterName.phi0.ordinal())));
+ aida.histogram1D("omega ").fill(trk.getTrackParameter(LCIOParameters.ParameterName.omega.ordinal()));
+ aida.histogram1D("tan(lambda) ").fill(trk.getTrackParameter(LCIOParameters.ParameterName.tanLambda.ordinal()));
+ aida.histogram1D("z0 ").fill(trk.getTrackParameter(LCIOParameters.ParameterName.z0.ordinal()));
+ }
+ }
+
+ @Override
+ public void fillEndOfRunPlots() {
+ }
+
+ @Override
+ public void dumpDQMData() {
+ }
+
+}
SVNspam 0.1