Print

Print


Author: [log in to unmask]
Date: Wed Apr 22 23:47:00 2015
New Revision: 2795

Log:
Analysis driver used to study SVT clusters.

Added:
    java/trunk/users/src/main/java/org/hps/users/omoreno/SvtClusterAnalysis.java   (with props)

Added: java/trunk/users/src/main/java/org/hps/users/omoreno/SvtClusterAnalysis.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/omoreno/SvtClusterAnalysis.java	(added)
+++ java/trunk/users/src/main/java/org/hps/users/omoreno/SvtClusterAnalysis.java	Wed Apr 22 23:47:00 2015
@@ -0,0 +1,295 @@
+package org.hps.users.omoreno;
+
+import java.io.IOException;
+import java.util.List;
+import java.util.Map;
+import java.util.HashMap;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IHistogramFactory;
+import hep.aida.IPlotterFactory;
+import hep.aida.IPlotter;
+import hep.aida.IHistogram1D;
+import hep.aida.IPlotterStyle;
+import hep.aida.ITree;
+import hep.aida.ref.rootwriter.RootFileStore;
+
+import org.hps.recon.tracking.FittedRawTrackerHit;
+import org.lcsim.recon.tracking.digitization.sisim.SiTrackerHitStrip1D;
+import org.lcsim.util.Driver; 
+import org.lcsim.geometry.Detector;
+import org.lcsim.detector.tracker.silicon.HpsSiSensor;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.RawTrackerHit;
+
+/**
+ * 
+ * @author Omar Moreno
+ *
+ */
+public class SvtClusterAnalysis extends Driver {
+   
+    
+    // Use JFreeChart as the default plotting backend
+    static { 
+        hep.aida.jfree.AnalysisFactory.register();
+    }
+
+    private List<HpsSiSensor> sensors;
+    private Map<RawTrackerHit, FittedRawTrackerHit> fittedRawTrackerHitMap 
+        = new HashMap<RawTrackerHit, FittedRawTrackerHit>();
+  
+    // Plotting
+    ITree tree; 
+    IHistogramFactory histogramFactory; 
+	IPlotterFactory plotterFactory = IAnalysisFactory.create().createPlotterFactory();
+	
+	protected Map<String, IPlotter> plotters = new HashMap<String, IPlotter>(); 
+	private Map<HpsSiSensor, IHistogram1D> clusterChargePlots = new HashMap<HpsSiSensor, IHistogram1D>();
+	private Map<HpsSiSensor, IHistogram1D> signalToNoisePlots = new HashMap<HpsSiSensor, IHistogram1D>();
+	private Map<HpsSiSensor, IHistogram1D> singleHitClusterChargePlots = new HashMap<HpsSiSensor, IHistogram1D>();
+	private Map<HpsSiSensor, IHistogram1D> singleHitSignalToNoisePlots = new HashMap<HpsSiSensor, IHistogram1D>();
+	private Map<HpsSiSensor, IHistogram1D> multHitClusterChargePlots = new HashMap<HpsSiSensor, IHistogram1D>();
+	private Map<HpsSiSensor, IHistogram1D> multHitSignalToNoisePlots = new HashMap<HpsSiSensor, IHistogram1D>();
+    
+    // Detector name
+    private static final String SUBDETECTOR_NAME = "Tracker";
+    
+    // Collections
+    private String clusterCollectionName = "StripClusterer_SiTrackerHitStrip1D";
+    private String fittedHitsCollectionName = "SVTFittedRawTrackerHits";
+    
+    private int runNumber = -1; 
+    
+    /**
+     * Default Ctor
+     */
+    public SvtClusterAnalysis() { }
+    
+    private int computePlotterRegion(HpsSiSensor sensor) {
+
+		if (sensor.getLayerNumber() < 7) {
+			if (sensor.isTopLayer()) {
+				return 6*(sensor.getLayerNumber() - 1); 
+			} else { 
+				return 6*(sensor.getLayerNumber() - 1) + 1;
+			} 
+		} else { 
+		
+			if (sensor.isTopLayer()) {
+				if (sensor.getSide() == HpsSiSensor.POSITRON_SIDE) {
+					return 6*(sensor.getLayerNumber() - 7) + 2;
+				} else { 
+					return 6*(sensor.getLayerNumber() - 7) + 3;
+				}
+			} else if (sensor.isBottomLayer()) {
+				if (sensor.getSide() == HpsSiSensor.POSITRON_SIDE) {
+					return 6*(sensor.getLayerNumber() - 7) + 4;
+				} else {
+					return 6*(sensor.getLayerNumber() - 7) + 5;
+				}
+			}
+		}
+		return -1; 
+    }
+    
+    protected void detectorChanged(Detector detector) {
+       
+        tree = IAnalysisFactory.create().createTreeFactory().create();
+        histogramFactory = IAnalysisFactory.create().createHistogramFactory(tree);
+        
+        // Get the HpsSiSensor objects from the tracker detector element
+        sensors = detector.getSubdetector(SUBDETECTOR_NAME)
+                          .getDetectorElement().findDescendants(HpsSiSensor.class);
+   
+        // If the detector element had no sensors associated with it, throw
+        // an exception
+        if (sensors.size() == 0) {
+            throw new RuntimeException("No sensors were found in this detector.");
+        }
+       
+        plotters.put("Cluster Amplitude", plotterFactory.create("Cluster Amplitude"));
+        plotters.get("Cluster Amplitude").createRegions(6, 6);
+       
+        plotters.put("Signal to Noise", plotterFactory.create("Signal to Noise"));
+        plotters.get("Signal to Noise").createRegions(6, 6);
+        
+        for (HpsSiSensor sensor : sensors) { 
+            
+            clusterChargePlots.put(sensor, 
+                    histogramFactory.createHistogram1D(sensor.getName() + " - Cluster Charge", 100, 0, 5000));
+            plotters.get("Cluster Amplitude").region(this.computePlotterRegion(sensor))
+                                             .plot(clusterChargePlots.get(sensor), this.createStyle(1, "Cluster Amplitude [ADC Counts]", ""));
+        
+            singleHitClusterChargePlots.put(sensor, 
+                    histogramFactory.createHistogram1D(sensor.getName() + " - Single Hit Cluster Charge", 100, 0, 5000));
+            plotters.get("Cluster Amplitude").region(this.computePlotterRegion(sensor))
+                                             .plot(singleHitClusterChargePlots.get(sensor), this.createStyle(2, "Cluster Amplitude [ADC Counts]", ""));
+
+            multHitClusterChargePlots.put(sensor, 
+                    histogramFactory.createHistogram1D(sensor.getName() + " - Multiple Hit Cluster Charge", 100, 0, 5000));
+            plotters.get("Cluster Amplitude").region(this.computePlotterRegion(sensor))
+                                             .plot(multHitClusterChargePlots.get(sensor), this.createStyle(2, "Cluster Amplitude [ADC Counts]", ""));
+            
+            signalToNoisePlots.put(sensor, 
+                    histogramFactory.createHistogram1D(sensor.getName() + " - Signal to Noise", 50, 0, 50));
+            plotters.get("Signal to Noise").region(this.computePlotterRegion(sensor))
+                                           .plot(signalToNoisePlots.get(sensor), this.createStyle(1, "Signal to Noise", ""));
+
+            singleHitSignalToNoisePlots.put(sensor, 
+                    histogramFactory.createHistogram1D(sensor.getName() + " - Single Hit Signal to Noise", 50, 0, 50));
+            plotters.get("Signal to Noise").region(this.computePlotterRegion(sensor))
+                                           .plot(singleHitSignalToNoisePlots.get(sensor), this.createStyle(2, "Signal to Noise", ""));
+        
+            multHitSignalToNoisePlots.put(sensor, 
+                    histogramFactory.createHistogram1D(sensor.getName() + " - Multiple Hit Signal to Noise", 50, 0, 50));
+            plotters.get("Signal to Noise").region(this.computePlotterRegion(sensor))
+                                           .plot(multHitSignalToNoisePlots.get(sensor), this.createStyle(2, "Signal to Noise", ""));
+        }
+        
+		for (IPlotter plotter : plotters.values()) { 
+			plotter.show();
+		}
+    }
+
+    public void process(EventHeader event) { 
+     
+        if (runNumber == -1) runNumber = event.getRunNumber();
+        
+        // If the event doesn't contain fitted raw hits, skip it
+        if (!event.hasCollection(FittedRawTrackerHit.class, fittedHitsCollectionName)) return;
+        
+        // Get the list of fitted hits from the event
+        List<FittedRawTrackerHit> fittedHits = event.get(FittedRawTrackerHit.class, fittedHitsCollectionName);
+        
+        // Map the fitted hits to their corresponding raw hits
+        this.mapFittedRawHits(fittedHits);
+        
+        // If the event doesn't contain any clusters, skip it
+        if (!event.hasCollection(SiTrackerHitStrip1D.class, clusterCollectionName)) return;
+        
+        // Get the list of clusters in the event
+        List<SiTrackerHitStrip1D> clusters = event.get(SiTrackerHitStrip1D.class, clusterCollectionName);
+       
+        for (SiTrackerHitStrip1D cluster : clusters) { 
+            
+            // Get the sensor associated with this cluster
+            HpsSiSensor sensor = (HpsSiSensor) cluster.getSensor();
+            
+            // Get the raw hits composing this cluster and use them to calculate the amplitude of the hit
+            double amplitude = 0;
+            double noise = 0;
+            for (RawTrackerHit rawHit : cluster.getRawHits()) {
+                
+                // Get the channel of the raw hit
+                int channel = rawHit.getIdentifierFieldValue("strip");
+                
+                // Add the amplitude of that channel to the total amplitude
+                amplitude += this.getFittedHit(rawHit).getAmp();
+                
+                // Calculate the mean noise for the channel
+                double channelNoise = 0;
+                for (int sampleN = 0; sampleN < 6; sampleN++) { 
+                    channelNoise += sensor.getNoise(channel, sampleN);
+                }
+                channelNoise = channelNoise/6;
+                
+                noise += channelNoise * this.getFittedHit(rawHit).getAmp();
+            }
+       
+            // Calculate the signal weighted noise
+            noise = noise/amplitude;
+            
+            // Fill all plots
+            clusterChargePlots.get(sensor).fill(amplitude);
+            signalToNoisePlots.get(sensor).fill(amplitude/noise);
+            
+            if (cluster.getRawHits().size() == 1) { 
+                singleHitClusterChargePlots.get(sensor).fill(amplitude);
+                singleHitSignalToNoisePlots.get(sensor).fill(amplitude/noise);
+            } else { 
+                multHitClusterChargePlots.get(sensor).fill(amplitude);
+                multHitSignalToNoisePlots.get(sensor).fill(amplitude/noise);
+            }
+        }
+    }
+    
+    public void endOfData() { 
+        
+        String rootFile = "run" + runNumber + "_cluster_analysis.root";
+        RootFileStore store = new RootFileStore(rootFile);
+        try {
+            store.open();
+            store.add(tree);
+            store.close(); 
+        } catch (IOException e) {
+            e.printStackTrace();
+        }
+    }
+    
+    /**
+     *  Method that creates a map between a fitted raw hit and it's corresponding raw fit
+     *  
+     * @param fittedHits : List of fitted hits to map
+     */
+    private void mapFittedRawHits(List<FittedRawTrackerHit> fittedHits) { 
+        
+        // Clear the fitted raw hit map of old values
+        fittedRawTrackerHitMap.clear();
+       
+        // Loop through all fitted hits and map them to their corresponding raw hits
+        for (FittedRawTrackerHit fittedHit : fittedHits) { 
+            fittedRawTrackerHitMap.put(fittedHit.getRawTrackerHit(), fittedHit);
+        }
+    }
+  
+    /**
+     * 
+     * @param rawHit
+     * @return
+     */
+    private FittedRawTrackerHit getFittedHit(RawTrackerHit rawHit) { 
+        return fittedRawTrackerHitMap.get(rawHit);
+    }
+    
+    IPlotterStyle createStyle(int color, String xAxisTitle, String yAxisTitle) { 
+       
+        // Create a default style
+        IPlotterStyle style = this.plotterFactory.createPlotterStyle();
+        
+        // Set the style of the X axis
+        style.xAxisStyle().setLabel(xAxisTitle);
+        style.xAxisStyle().labelStyle().setFontSize(14);
+        style.xAxisStyle().setVisible(true);
+        
+        // Set the style of the Y axis
+        style.yAxisStyle().setLabel(yAxisTitle);
+        style.yAxisStyle().labelStyle().setFontSize(14);
+        style.yAxisStyle().setVisible(true);
+        
+        // Turn off the histogram grid 
+        style.gridStyle().setVisible(false);
+        
+        // Set the style of the data
+        style.dataStyle().lineStyle().setVisible(false);
+        style.dataStyle().outlineStyle().setVisible(false);
+        style.dataStyle().outlineStyle().setThickness(4);
+        style.dataStyle().fillStyle().setVisible(true);
+        
+        if (color == 1) { 
+            style.dataStyle().fillStyle().setColor("31, 137, 229, 1");
+            style.dataStyle().outlineStyle().setColor("31, 137, 229, 1");
+            style.dataStyle().fillStyle().setOpacity(.30);
+        } else if (color == 2) { 
+            style.dataStyle().fillStyle().setColor("93, 228, 47, 1");
+            style.dataStyle().outlineStyle().setColor("93, 228, 47, 1");
+            style.dataStyle().fillStyle().setOpacity(.70);
+        }
+        style.dataStyle().errorBarStyle().setVisible(false);
+        
+        // Turn off the legend
+        style.legendBoxStyle().setVisible(false);
+       
+        return style;
+    }
+}