Print

Print


Commit in java/trunk/analysis/src/main/java/org/hps/analysis/dataquality on MAIN
DataQualityMonitor.java+89added 533
EcalMonitoring.java+11added 533
HitMCEfficiency.java+163added 533
ReconMonitoring.java+11added 533
SvtMonitoring.java+187added 533
TrackMCEfficiency.java+362added 533
TrackingMonitoring.java+123added 533
+946
7 added files
some more data quality stuff

java/trunk/analysis/src/main/java/org/hps/analysis/dataquality
DataQualityMonitor.java added at 533
--- 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
EcalMonitoring.java added at 533
--- 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
HitMCEfficiency.java added at 533
--- 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
ReconMonitoring.java added at 533
--- 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
SvtMonitoring.java added at 533
--- 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
TrackMCEfficiency.java added at 533
--- 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
TrackingMonitoring.java added at 533
--- 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