 |
 |
| Commit in hps-java/src/main/java/org/lcsim/hps/recon/ecal on MAIN | |||
| GainCalibrationDriver.java | +705 | added 1.1 | |
| GainCalib.java | +591 | added 1.1 | |
| +1296 | |||
A start for gain calibration.
diff -N GainCalibrationDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ GainCalibrationDriver.java 24 Jul 2012 23:28:36 -0000 1.1 @@ -0,0 +1,705 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.hps.recon.ecal;
+
+import hep.aida.*;
+import hep.aida.ref.histogram.Histogram1D;
+import hep.aida.ref.plotter.PlotterRegion;
+import java.awt.Button;
+import java.awt.FlowLayout;
+import java.awt.Frame;
+import java.awt.event.ActionEvent;
+import java.awt.event.ActionListener;
+import java.awt.event.WindowEvent;
+import java.awt.event.WindowListener;
+import java.io.*;
+import java.util.ArrayList;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import javax.swing.JButton;
+import javax.swing.JComboBox;
+import javax.swing.JLabel;
+import org.lcsim.detector.identifier.ExpandedIdentifier;
+import org.lcsim.detector.identifier.IExpandedIdentifier;
+import org.lcsim.detector.identifier.IIdentifier;
+import org.lcsim.detector.identifier.IIdentifierHelper;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.Track;
+import org.lcsim.event.base.BaseRawCalorimeterHit;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.IDDecoder;
+import org.lcsim.geometry.Subdetector;
+import org.lcsim.hps.monitoring.AIDAFrame;
+import org.lcsim.hps.monitoring.Redrawable;
+import org.lcsim.hps.monitoring.Resettable;
+import org.lcsim.hps.recon.ecal.HPSEcalRawConverterDriver;
+import org.lcsim.hps.recon.ecal.HPSEcalCluster;
+import org.lcsim.hps.recon.ecal.HPSEcalConditions;
+import org.lcsim.hps.recon.tracking.EcalTrackMatch;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.util.lcio.LCIOConstants;
+
+/**
+ *
+ * @author phansson+
+ */
+public class GainCalibrationDriver extends Driver implements Resettable, ActionListener, Redrawable {
+
+ int nevents = 0;
+ boolean debug = true;
+ protected IDDecoder dec = null;
+ protected Subdetector ecal;
+ private String ecalName = "Ecal";
+ String rawCollectionName = "EcalReadoutHits";
+ String ecalReadoutName = "EcalHits";
+ String ecalTmpCollectionName = "EcalCalHitsTmp";
+ private String trackCollectionName = "MatchedTracks";
+ private String ecalClusterCollectionName = "EcalClusters";
+ private String outputPlotFileName = "test.aida";
+ private String ecalGainFileName = "";
+ private String ecalGainCalibCorrFileName = "";
+
+ //I would like to steel these from the converter -> FIX THIS!
+ private boolean applyBadCrystalMap = true;
+ private boolean dropBadFADC = false;
+ double threshold = Double.NEGATIVE_INFINITY;
+ private int integralWindow = 30;
+ private int _numIterations = 0;
+
+
+ HPSEcalClusterer clusterer = null;
+ EcalTrackMatch trkMatchTool = null;
+ HPSEcalRawConverter converter = null;
+
+ GainCalib gainCalib = new GainCalib();
+
+ private boolean hideFrame = false;
+ private int refreshRate = 100;
+ private AIDA aida = AIDA.defaultInstance();
+ private IAnalysisFactory af = aida.analysisFactory();
+ IHistogramFactory hf = aida.histogramFactory();
+ private AIDAFrame plotterFrame;
+ private AIDAFrame pePlotterFrame;
+ IPlotter plotter;
+ JComboBox xCombo;
+ JLabel xLabel;
+ JComboBox yCombo;
+ JLabel yLabel;
+ Integer xList[];
+ Integer yList[];
+ JButton blankButton;
+
+ boolean alive = true;
+
+
+
+ IHistogram1D pePlots[][][] = new IHistogram1D[47][11][5];
+ IHistogram2D mpePlot;
+ IHistogram2D spePlot;
+ IHistogram2D hitmap;
+ IHistogram1D[] h_PE_t = new IHistogram1D[5];
+ IHistogram1D[] h_PE_b = new IHistogram1D[5];
+
+
+ class AL extends Frame implements ActionListener,WindowListener {
+ Button b;
+
+ public AL(String title) {
+ super(title);
+ setLayout(new FlowLayout());
+ b= new Button("Click to stop program");
+ add(b);
+ b.addActionListener(this);
+ }
+ @Override
+ public void actionPerformed(ActionEvent ae) {
+ if(ae.getSource() == b) {
+ System.out.println("Clicked to exit!");
+ dispose();
+ alive=false;//System.exit(0);
+ }
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ @Override
+ public void windowOpened(WindowEvent we) {
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ @Override
+ public void windowClosing(WindowEvent we) {
+ System.out.println("windowClosing to exit!");
+ dispose();
+ alive=false;//System.exit(0);
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ @Override
+ public void windowClosed(WindowEvent we) {
+ System.out.println("windowClosing to exit!");
+ dispose();
+ alive=false;//System.exit(0);
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ @Override
+ public void windowIconified(WindowEvent we) {
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ @Override
+ public void windowDeiconified(WindowEvent we) {
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ @Override
+ public void windowActivated(WindowEvent we) {
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ @Override
+ public void windowDeactivated(WindowEvent we) {
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ }
+
+
+
+ public void startOfData() {
+ }
+
+ public void detectorChanged(Detector detector) {
+ // Get the Subdetector.
+ ecal = detector.getSubdetector(ecalName);
+
+ clusterer.setEcalName("Ecal");
+ clusterer.detectorChanged(detector);
+
+ // Cache ref to decoder.
+ dec = ecal.getIDDecoder();
+
+
+
+
+ gainCalib.setClusterer(clusterer);
+ gainCalib.setConverter(converter);
+ gainCalib.setDetector(detector);
+ gainCalib.setEcalReadoutName(ecalReadoutName);
+ gainCalib.setEcalName(ecalName);
+ gainCalib.setThreshold(threshold);
+ gainCalib.setIntegralWindow(integralWindow);
+ gainCalib.setApplyBadCrystalMap(this.applyBadCrystalMap);
+
+ //Use gain calibration file
+ System.out.println("Using gain calibration file \""+ecalGainFileName+"\"");
+ if(!"".equals(ecalGainFileName)) {
+ gainCalib.setAllGain(ecalGainFileName);
+ } else {
+ gainCalib.setAllGain(converter.gain);
+ }
+
+
+ pePlotterFrame = new AIDAFrame();
+ pePlotterFrame.setTitle("Gain Frame");
+
+ //plotterFrame = new AIDAFrame();
+ //plotterFrame.setTitle("Gain General");
+
+ IPlotterStyle style;
+
+
+ IPlotter plotter_hitmap_gr = af.createPlotterFactory().create();
+ plotter_hitmap_gr.createRegions(1,3,0);
+ plotter_hitmap_gr.setTitle("Cluster hit map");
+ plotter_hitmap_gr.style().statisticsBoxStyle().setVisible(false);
+ pePlotterFrame.addPlotter(plotter_hitmap_gr);
+
+ hitmap = aida.histogram2D("Cluster hit map", 47, -23.5, 23.5, 11, -5.5, 5.5);
+ plotter_hitmap_gr.region(0).plot(hitmap);
+
+ style = plotter_hitmap_gr.region(0).style();
+ style.setParameter("hist2DStyle", "colorMap");
+ style.dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ ((PlotterRegion) plotter_hitmap_gr.region(0)).getPlot().setAllowUserInteraction(true);
+ ((PlotterRegion) plotter_hitmap_gr.region(0)).getPlot().setAllowPopupMenus(true);
+
+ IPlotter[] plotter_PoverE = new IPlotter[5];
+
+
+ for(int iE=0;iE<=4;++iE) {
+
+ String str = iE==0?"":(" iE="+iE);
+
+ h_PE_t[iE] = aida.histogram1D("E over p top"+str,50,0,2);
+ h_PE_b[iE] = aida.histogram1D("E over p bottom"+str,50,0,2);
+
+ plotter_PoverE[iE] = af.createPlotterFactory().create();
+ plotter_PoverE[iE].createRegions(1,2,0);
+ plotter_PoverE[iE].setTitle("E over P"+str);
+ plotter_PoverE[iE].style().statisticsBoxStyle().setVisible(true);
+ pePlotterFrame.addPlotter(plotter_PoverE[iE]);
+
+ plotter_PoverE[iE].region(0).plot(h_PE_t[iE]);
+ plotter_PoverE[iE].region(1).plot(h_PE_b[iE]);
+ }
+
+ plotter = af.createPlotterFactory().create();
+ plotter.createRegions(1,3,0);
+ plotter.setTitle("Gain Plots");
+
+ pePlotterFrame.addPlotter(plotter);
+
+ mpePlot = aida.histogram2D("<E over p>",47,-23.5,23.5,11,-5.5,5.5);
+ plotter.region(0).plot(mpePlot);
+ spePlot = aida.histogram2D("RMS(E over p)",47,-23.5,23.5,11,-5.5,5.5);
+ plotter.region(1).plot(spePlot);
+ plotter.region(0).style().statisticsBoxStyle().setVisible(false);
+ plotter.region(0).style().setParameter("hist2DStyle", "colorMap");
+ plotter.region(0).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ plotter.region(1).style().statisticsBoxStyle().setVisible(false);
+ plotter.region(1).style().setParameter("hist2DStyle", "colorMap");
+ plotter.region(1).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+
+
+
+ for(int iE=0;iE<=4;++iE) {
+ for(int irow=-5;irow<=5;++irow) {
+ for(int icol=-23;icol<=23;++icol) {
+ if(iE==0) {
+ pePlots[icol+23][irow+5][iE] = aida.histogram1D("E over p x=" + icol + " y=" + irow, 50, 0,2);
+ } else {
+ pePlots[icol+23][irow+5][iE] = aida.histogram1D("E over p x=" + icol + " y=" + irow + " iE="+iE, 50, 0,2);
+ }
+ }
+ }
+ }
+
+
+ xList = new Integer[46];
+ yList = new Integer[10];
+ int in=0;
+ for(int i=-5;i<=5;++i) {
+ if(i!=0) {
+ yList[in] = i;
+ ++in;
+ }
+ }
+ in=0;
+ for(int i=-23;i<=23;++i) {
+ if(i!=0) {
+ xList[in] = i;
+ ++in;
+ }
+ }
+
+
+ xCombo = new JComboBox(xList);
+ xCombo.addActionListener(this);
+ xLabel = new JLabel("x");
+ xLabel.setLabelFor(xCombo);
+ pePlotterFrame.getControlsPanel().add(xLabel);
+ pePlotterFrame.getControlsPanel().add(xCombo);
+ yCombo = new JComboBox(yList);
+ yCombo.addActionListener(this);
+ yLabel = new JLabel("y");
+ yLabel.setLabelFor(xCombo);
+ pePlotterFrame.getControlsPanel().add(yLabel);
+ pePlotterFrame.getControlsPanel().add(yCombo);
+
+
+ plotter.region(2).plot(pePlots[-5 + 23][2 + 5 - 1][0]);
+ xCombo.setSelectedIndex(-5 + 23);
+ yCombo.setSelectedIndex(2 + 5 - 1);
+
+ blankButton = new JButton("Hide histogram");
+ pePlotterFrame.getControlsPanel().add(blankButton);
+ blankButton.addActionListener(this);
+
+
+
+
+
+ if(!hideFrame) {
+
+ //plotterFrame.pack();
+ //plotterFrame.setVisible(true);
+
+ pePlotterFrame.pack();
+ pePlotterFrame.setVisible(true);
+
+ }
+
+
+
+
+
+
+
+ }
+
+ public GainCalibrationDriver() {
+ converter = new HPSEcalRawConverter();
+ trkMatchTool = new EcalTrackMatch(false);
+ clusterer = new HPSEcalClusterer();
+ gainCalib = new GainCalib();
+
+
+
+ }
+
+
+
+ public void setNumIterations(int i) {
+ this._numIterations = i;
+ }
+
+ public void setDebug(boolean flag) {
+ this.debug = flag;
+ }
+
+ public void setOutputPlotFileName( String name ) {
+ this.outputPlotFileName = name;
+ }
+
+ public void setEcalGainFileName( String name ) {
+ this.ecalGainFileName = name;
+ }
+ public void setEcalGainCalibCorrFileName( String name ) {
+ this.ecalGainCalibCorrFileName = name;
+ }
+ public void setHideFrame( boolean val ) {
+ this.hideFrame = val;
+ }
+
+
+
+ public void process(EventHeader event) {
+ ++nevents;
+ if( debug ) {
+ System.out.println("Processing event " + nevents);
+ }
+
+ if( refreshRate > 0 && nevents % refreshRate == 0 ) {
+ redraw();
+ }
+
+ // Find the raw hits used to build calorimeter hits
+
+/*
+ List<RawTrackerHit> hits = null;
+ if (event.hasCollection(RawTrackerHit.class, rawCollectionName)) {
+ hits = event.get(RawTrackerHit.class, rawCollectionName);
+ }
+
+ if(hits==null) {
+ if(debug) System.out.println("No raw ecal hits");
+ return;
+ }
+*/
+
+ List<BaseRawCalorimeterHit> hits = null;
+ if (event.hasCollection(BaseRawCalorimeterHit.class, rawCollectionName)) {
+ hits = event.get(BaseRawCalorimeterHit.class, rawCollectionName);
+ }
+
+ if(hits==null) {
+ if(debug) System.out.println("No raw ecal hits ");
+ return;
+ }
+
+
+ if(debug) System.out.println(hits.size() + " raw ecal hits");
+
+ if(hits.isEmpty()) return;
+
+
+
+
+ ArrayList<CalorimeterHit> newHits = new ArrayList<CalorimeterHit>();
+
+
+
+
+
+ for (BaseRawCalorimeterHit hit : hits) {
+ double g = gainCalib.getGain(hit);
+ CalorimeterHit newHit = converter.HitDtoA(hit,integralWindow);
+ if(applyBadCrystalMap && HPSEcalRawConverterDriver.isBadCrystal(newHit)) continue;
+ if (dropBadFADC && HPSEcalRawConverterDriver.isBadFADC(newHit)) continue;
+ if (newHit.getRawEnergy() > threshold) {
+ newHits.add(newHit);
+ }
+ }
+
+ if(debug) System.out.println(newHits.size() + " calibrated calorimeter hits");
+
+ //To be able to use decoders I need to associate with readout. Use event put to do that?! -> FIX THIS
+ event.put(ecalTmpCollectionName, newHits, CalorimeterHit.class, 0, ecalReadoutName);
+
+ //Get the newly created hits
+
+ List<CalorimeterHit> newCalHits = null;
+ if (event.hasCollection(CalorimeterHit.class, ecalTmpCollectionName)) {
+ newCalHits = event.get(CalorimeterHit.class, ecalTmpCollectionName);
+ }
+
+ if(newCalHits==null) {
+ if(debug) System.out.println("No newCalHits ");
+ System.exit(1);
+
+ }
+
+
+ if(debug) System.out.println(newCalHits.size() + " newCalHits");
+
+ if(newCalHits.isEmpty()) return;
+
+
+ //Now we have a set of calibrated ecal hits
+ // Clusterize them to get a starting point -> copy the clusterer
+
+ // Make a hit map for quick lookup by ID.
+ Map<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
+ for (CalorimeterHit hit : newCalHits) {
+ hitMap.put(hit.getCellID(), hit);
+ }
+
+ List<HPSEcalCluster> clusters = clusterer.createClusters(hitMap);
+
+ if(debug) System.out.println(clusters.size() + " clusters found");
+
+
+ // Now we want to produce the corrections for each module: c = p/E
+
+
+
+ List<Track> tracks = null;
+ if(event.hasCollection(Track.class, trackCollectionName)) {
+ tracks = event.get(Track.class, trackCollectionName);
+ }
+
+ if(debug) System.out.println( tracks.size() + " tracks in this event");
+
+ if(tracks==null) return;
+
+
+
+ for(HPSEcalCluster cl : clusters) {
+
+
+
+
+ int[] pos = getCrystalPair(cl);
+ int side = pos[1]>0 ? 0 : 1; //top or bottom
+
+ if(debug) System.out.println("Looking at cluster at ix=" + pos[0] + " iy=" + pos[1]);
+
+ trkMatchTool.setCluster(cl);
+ trkMatchTool.match(tracks);
+
+ if(!trkMatchTool.isMatchedY(50)) {
+ if(debug) System.out.println("Cluster not matched to a track");
+ continue;
+ }
+
+
+ if(debug) System.out.println("Cluster matched to track at distance Y " + trkMatchTool.getDistanceToTrackInY() + " and X " + trkMatchTool.getDistanceToTrackInX());
+
+
+ double P = trkMatchTool.getMatchedTrack().getPX()*1000;
+ double E = cl.getEnergy();
+ double samplingFraction = 0.9;
+ double Ep = E*samplingFraction;
+ double Eoverp = Ep/P;
+
+ if(debug) System.out.println("P " + P + " E " + E + " Ep " + Ep);
+
+
+ double Eseed = cl.getSeedHit().getRawEnergy();
+ double ErawSum = 0;
+ for(CalorimeterHit hit : cl.getCalorimeterHits()) {
+ ErawSum += hit.getRawEnergy();
+ }
+
+ if(Eseed/ErawSum<0.6) continue;
+
+ int ebin = -1;
+ if(P>500&&P<=700) ebin=1;
+ else if(P>700&&P<=900) ebin=2;
+ else if(P>900&&P<=1100) ebin=3;
+ else ebin=4;
+
+ if(side==0) {
+ h_PE_t[0].fill(Eoverp);
+ h_PE_t[ebin].fill(Eoverp);
+
+ }
+ else {
+ h_PE_b[0].fill(Eoverp);
+ h_PE_b[ebin].fill(Eoverp);
+ }
+
+ hitmap.fill(pos[0], pos[1]);
+
+ pePlots[pos[0]+23][pos[1]+5][0].fill(Eoverp);
+
+ pePlots[pos[0]+23][pos[1]+5][ebin].fill(Eoverp);
+
+
+ gainCalib.setOriginalHits(event.getEventNumber(), hits, converter.gain);
+ for(BaseRawCalorimeterHit hit : hits) {
+ gainCalib.setTrack(event.getEventNumber(), hit.getCellID(), trkMatchTool.getMatchedTrack());
+ }
+
+
+
+ }
+
+
+
+
+ }
+
+
+
+
+
+ public int[] getCrystalPair(HPSEcalCluster cluster) {
+ int[] pos = new int[2];
+ pos[0] = cluster.getSeedHit().getIdentifierFieldValue("ix");
+ pos[1] = cluster.getSeedHit().getIdentifierFieldValue("iy");
+
+ //System.out.println("cluster ix,iy " + pos[0] + "," + pos[1] + " from pos " + cluster.getSeedHit().getPositionVec().toString());
+ return pos;
+ //getCrystalPair(cluster.getPosition());
+ }
+
+
+
+
+
+
+
+
+
+ public void endOfData() {
+
+ redraw();
+
+ if (outputPlotFileName != "")
+ try {
+ aida.saveAs(outputPlotFileName);
+ } catch (IOException ex) {
+ Logger.getLogger(GainCalibrationDriver.class.getName()).log(Level.SEVERE, "Couldn't save aida plots to file " + outputPlotFileName, ex);
+ }
+ //displayFastTrackingPlots();
+
+ PrintWriter out=null;
+ try {
+ out = new PrintWriter(ecalGainCalibCorrFileName);
+ } catch (FileNotFoundException ex) {
+ Logger.getLogger(GainCalibrationDriver.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ //Save a new file with corrected gains
+ for(int iE=0;iE<=4;++iE) {
+ for(int irow=-5;irow<=5;++irow) {
+ for(int icol=-23;icol<=23;++icol) {
+ IHistogram1D h = pePlots[icol+23][irow+5][iE];
+ double g = gainCalib.getGain(icol,irow);
+ int N = h.entries();
+ if(N>6) {
+ double m = h.mean();
+ double c = 1.0/m;
+ g = g*c;
+ out.println(iE + " " + icol + " " + irow + " " + g);
+ }
+
+ }
+ }
+ }
+ out.close();
+
+ //gainCalib.runOptimization(this._numIterations);
+
+ //AL alWindow = new AL("Click to quit");
+ //alWindow.setSize(350,100);
+ //alWindow.setVisible(true);
+ //while(alive) {
+ // sleep
+ //}
+
+
+ }
+
+ @Override
+ public void reset() {
+ if (plotter != null) {
+ plotter.hide();
+ plotter.destroyRegions();
+ for (int x = -23; x <= 23; x++) { // slot
+ for (int y = -5; y <= 5; y++) { // crate
+ for(int iE=0;iE<=4;++iE) {
+ pePlots[x + 23][y + 5][iE].reset();
+ }
+ }
+ }
+ }
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ @Override
+ public void actionPerformed(ActionEvent ae) {
+
+ if(ae.getSource()==blankButton) {
+ plotter.region(2).clear();
+ } else {
+ //get the col and row to display
+ Integer x = (Integer) xCombo.getSelectedItem();
+ Integer y = (Integer) yCombo.getSelectedItem();
+ plotter.region(2).clear();
+ plotter.region(2).plot(pePlots[x+23][y+5][0]);
+ }
+
+
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ @Override
+ public void redraw() {
+
+ //do something if needed
+ mpePlot.reset();
+ spePlot.reset();
+
+ for(int irow=-5;irow<=5;++irow) {
+ for(int icol=-23;icol<=23;++icol) {
+ //System.out.println(" redraw irow " + irow + " icol " + icol + " entries " + pePlots[icol+23][irow+5].entries());
+ if(pePlots[icol+23][irow+5][0].entries()>10) {
+ mpePlot.fill(icol,irow,pePlots[icol+23][irow+5][0].mean());
+ spePlot.fill(icol,irow,pePlots[icol+23][irow+5][0].rms());
+
+ }
+ }
+ }
+
+
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+ @Override
+ public void setEventRefreshRate(int eventRefreshRate) {
+ refreshRate = eventRefreshRate;
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+
+}
diff -N GainCalib.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ GainCalib.java 24 Jul 2012 23:28:36 -0000 1.1 @@ -0,0 +1,591 @@
+/*
+ * To change this template, choose Tools | Templates
+ * and open the template in the editor.
+ */
+package org.lcsim.hps.recon.ecal;
+
+import hep.aida.*;
+import java.awt.event.ActionEvent;
+import java.awt.event.ActionListener;
+import java.io.BufferedReader;
+import java.io.FileNotFoundException;
+import java.io.FileReader;
+import java.io.IOException;
+import java.util.*;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+import javax.swing.JComboBox;
+import javax.swing.JLabel;
+import org.lcsim.detector.identifier.ExpandedIdentifier;
+import org.lcsim.detector.identifier.IExpandedIdentifier;
+import org.lcsim.detector.identifier.IIdentifierHelper;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Track;
+import org.lcsim.event.base.BaseLCSimEvent;
+import org.lcsim.event.base.BaseRawCalorimeterHit;
+import org.lcsim.geometry.Detector;
+import org.lcsim.hps.monitoring.AIDAFrame;
+import org.lcsim.hps.recon.tracking.EcalTrackMatch;
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * @author phansson
+ */
+public class GainCalib implements ActionListener {
+
+ boolean _debug = false;
+
+ EventHeader _event = null;
+ HPSEcalClusterer _clusterer = null;
+ HPSEcalRawConverter _converter = null;
+ private int _integralWindow = 30;
+ private double _threshold = Double.NEGATIVE_INFINITY;
+ private boolean _applyBadCrystalMap = true;
+
+ private Map<Long,Double> _gains = new HashMap<Long,Double>();
+ //original list
+ private Map<Integer,List<BaseRawCalorimeterHit>> _hits = new HashMap<Integer,List<BaseRawCalorimeterHit>>();
+ private Map<Integer,List<MatchedTrack>> _tracks = new HashMap<Integer,List<MatchedTrack>>();
+
+ private Detector _detector;
+ private String _ecalReadoutName;
+ private String _ecalName;
+ private String _ecalTmpCollectionName = "temporaryEcalCalHits";
+ private static IIdentifierHelper helper = null;
+ private static IExpandedIdentifier expId =null;
+
+
+ private AIDA aida = AIDA.defaultInstance();
+ private IAnalysisFactory af = aida.analysisFactory();
+ IHistogramFactory hf = aida.histogramFactory();
+ IHistogram1D pe[][] = new IHistogram1D[47][11];
+ List<PEHistograms> iterPE = new ArrayList<PEHistograms>();
+ private AIDAFrame plotFrame = new AIDAFrame();
+ private List<IPlotter> plotterList = new ArrayList<IPlotter>();
+
+ JComboBox xCombo=null;
+ JLabel xLabel=null;
+ Integer xList[];
+ JComboBox yCombo=null;
+ JLabel yLabel=null;
+ Integer yList[];
+
+ @Override
+ public void actionPerformed(ActionEvent ae) {
+ //if(ae == blankButton) {
+ // //
+ //}
+ Integer x = (Integer)xCombo.getSelectedItem();
+ Integer y = (Integer)yCombo.getSelectedItem();
+ for(IPlotter p : plotterList) {
+ String title = p.title(); //("Gain iteration " + iter);
+ String v[] = title.split(" ");
+ if(v.length<3) {
+ System.out.println("ERROR this title " + title + " doesn't have the correct syntax");
+ System.exit(1);
+ }
+ Integer iter = Integer.parseInt(v[2]);
+ for(PEHistograms peh : iterPE) {
+ if(peh._i==iter) {
+ p.region(4).clear();
+ p.region(4).plot(peh.pe[x+23][y+5]);
+ }
+ }
+
+ }
+
+
+
+ //throw new UnsupportedOperationException("Not supported yet.");
+ }
+
+
+
+
+
+ class MatchedTrack {
+ long _cellid;
+ Track _track;
+ MatchedTrack(long id,Track trk) {
+ _cellid = id;
+ _track = trk;
+ }
+ }
+ class PEHistograms {
+ int _i;
+ IHistogram1D pe[][] = new IHistogram1D[47][11];
+ PEHistograms(int i) {
+ _i = i;
+ }
+ }
+
+ GainCalib() {
+
+ for(int irow=-5;irow<=5;++irow) {
+ for(int icol=-23;icol<=23;++icol) {
+ pe[icol+23][irow+5] = aida.histogram1D("(GainCalib) E over p x=" + icol + " y=" + irow, 50, 0,3);
+
+ }
+ }
+
+
+
+
+ xList = new Integer[46];
+ yList = new Integer[10];
+ int in=0;
+ for(int i=-5;i<=5;++i) {
+ if(i!=0) {
+ yList[in] = i;
+ ++in;
+ }
+ }
+ in=0;
+ for(int i=-23;i<=23;++i) {
+ if(i!=0) {
+ xList[in] = i;
+ ++in;
+ }
+ }
+
+
+
+
+
+
+
+
+ }
+
+ public void setApplyBadCrystalMap(boolean d) {
+ this._applyBadCrystalMap = d;
+ }
+
+ public void setDebug(boolean d) {
+ this._debug = d;
+ }
+
+ public void setClusterer(HPSEcalClusterer clust) {
+ _clusterer = clust;
+ }
+
+ public void setConverter(HPSEcalRawConverter conv) {
+ _converter = conv;
+ }
+
+ public void setOriginalHits(int evt, List<BaseRawCalorimeterHit> hits, double gain) {
+ this._hits.put(evt, hits);
+ for(BaseRawCalorimeterHit hit: hits) {
+ this._gains.put(hit.getCellID(), gain);
+ }
+ }
+
+ public void setTrack(int evt, long ecalCellID, Track track) {
+ MatchedTrack mt = new MatchedTrack(ecalCellID,track);
+
+
+ if(!this._tracks.containsKey(evt)) {
+ this._tracks.put(evt, new ArrayList<MatchedTrack>());
+ }
+
+ List<MatchedTrack> list = this._tracks.get(evt);
+ list.add(mt);
+
+
+
+ }
+
+ public void setIntegralWindow(int val) {
+ _integralWindow = val;
+ }
+
+ public void setThreshold(double val) {
+ _threshold = val;
+ }
+
+ public void setEcalReadoutName(String str) {
+ this._ecalReadoutName = str;
+ }
+
+ public void setEcalName(String str) {
+ this._ecalName = str;
+ }
+
+ public void setDetector(Detector det) {
+ this._detector = det;
+ }
+
+
+ public void runOptimization(int nIterations) {
+
+ System.out.println("Start gain calibration");
+
+
+
+
+ System.out.println("nIterations " + nIterations);
+
+
+ if(_hits==null) {
+ System.out.println("No original hits available");
+ System.exit(1);
+ }
+
+ int event;
+ int iter = 0;
+ while(iter<nIterations) {
+ System.out.println("Iteration " + iter);
+
+ Iterator it = _hits.entrySet().iterator();
+ while(it.hasNext()) {
+ event = ((Map.Entry<Integer,ArrayList<BaseRawCalorimeterHit>>)it.next()).getKey();
+ if(_debug) System.out.println("Event " + event);
+
+ List<HPSEcalCluster> clusters = getCalibratedClusters(event);
+
+ if(_debug) System.out.println("Found " + clusters.size() + " clusters");
+
+ for(HPSEcalCluster cluster : clusters) {
+
+ int[] pos = new int[2];
+ pos[0] = cluster.getSeedHit().getIdentifierFieldValue("ix");
+ pos[1] = cluster.getSeedHit().getIdentifierFieldValue("iy");
+
+ //find the tracks matched to those clusters
+
+ MatchedTrack track = getTrack(event,cluster.getSeedHit().getCellID());
+
+ if(_debug) System.out.println("Found track with with cell id " + cluster.getSeedHit().getCellID());
+
+
+ double P = track._track.getPX()*1000;
+ double E = cluster.getEnergy();
+ double samplingFraction = 0.9;
+ double Ep = E*samplingFraction;
+ double Eoverp = Ep/P;
+
+
+ if(_debug)
+ System.out.println("P " + P + " E " + E + " Ep " + Ep + " --> E/p=" + Eoverp);
+
+
+ pe[pos[0]+23][pos[1]+5].fill(Eoverp);
+
+ if(_debug)
+ System.out.println("Cluster icol " + pos[0] + " irow " + pos[1] + " entries " + pe[pos[0]+23][pos[1]+5].entries());
+
+
+ }
+
+
+ }
+
+ copyPE(iter);
+
+
+ ++iter;
+ }
+
+
+
+
+
+
+ }
+
+
+ void copyPE(int iter) {
+
+ //copy the histograms, extract gains and reset for new iteration
+
+ IHistogram2D mpePlot = aida.histogram2D("Iter " + iter + " <E over p>",47,-23.5,23.5,11,-5.5,5.5);
+ IHistogram2D spePlot = aida.histogram2D("Iter " + iter + " RMS(E over p)",47,-23.5,23.5,11,-5.5,5.5);
+ IHistogram1D pePlot = hf.createCopy("Iter " + iter + " E over p", pe[0][0]);
+ IHistogram2D gainPlot = aida.histogram2D("Iter " + iter + " gain",47,-23.5,23.5,11,-5.5,5.5);
+
+
+ PEHistograms hists = new PEHistograms(iter);
+ for(int irow=-5;irow<=5;++irow) {
+ for(int icol=-23;icol<=23;++icol) {
+ hists.pe[icol+23][irow+5] = hf.createCopy("Iter " + iter + " E over p x=" + icol + " y=" + irow, pe[icol+23][irow+5]);
+ pePlot = hf.add("Iter " + iter + " E over p", pePlot, pe[icol+23][irow+5]);
+ if(pe[icol+23][irow+5].entries()>10) {
+ mpePlot.fill(icol,irow,pe[icol+23][irow+5].mean());
+ spePlot.fill(icol,irow,pe[icol+23][irow+5].rms());
+ double corr = 1.0/pe[icol+23][irow+5].mean(); //Correction to get p/E==1
+ double gain = this.getGain(icol, irow)*corr;
+ setGain(icol,irow,gain);
+ gain = this.getGain(icol, irow);
+ gainPlot.fill(icol,irow,gain);
+ }
+
+ }
+ }
+ hists._i = iter;
+ iterPE.add(hists);
+
+
+
+
+
+/*
+ for(int irow=-5;irow<=5;++irow) {
+ for(int icol=-23;icol<=23;++icol) {
+ for(PEHistograms peh : iterPE) {
+ if(peh._i==iter) {
+ if(peh.pe[icol+23][irow+5].entries()>10) {
+ mpePlot.fill(icol,irow,peh.pe[icol+23][irow+5].mean());
+ spePlot.fill(icol,irow,peh.pe[icol+23][irow+5].rms());
+ }
+ }
+ }
+ }
+ }
+
+ */
+
+
+ //Plot
+ if(xCombo==null) {
+ xCombo = new JComboBox(xList);
+ xCombo.addActionListener(this);
+ xLabel = new JLabel("x");
+ xLabel.setLabelFor(xCombo);
+ plotFrame.getControlsPanel().add(xLabel);
+ plotFrame.getControlsPanel().add(xCombo);
+ yCombo = new JComboBox(yList);
+ yCombo.addActionListener(this);
+ yLabel = new JLabel("y");
+ yLabel.setLabelFor(xCombo);
+ plotFrame.getControlsPanel().add(yLabel);
+ plotFrame.getControlsPanel().add(yCombo);
+ xCombo.setSelectedIndex(-15+23);
+ yCombo.setSelectedIndex(1+5-1);
+ }
+
+
+
+
+ IPlotter plotter = af.createPlotterFactory().create();
+ plotter.createRegions(1,5,0);
+ plotter.setTitle("Gain iteration " + iter);
+ //plotter.style().statisticsBoxStyle().setVisible(false);
+ plotFrame.addPlotter(plotter);
+ plotterList.add(plotter);
+
+
+ plotter.region(0).plot(mpePlot);
+ plotter.region(1).plot(spePlot);
+ plotter.region(2).plot(gainPlot);
+ plotter.region(3).plot(pePlot);
+ plotter.region(4).plot(hists.pe[-15+23][1+5-1+1]);
+
+ plotter.region(0).style().statisticsBoxStyle().setVisible(false);
+ plotter.region(0).style().setParameter("hist2DStyle", "colorMap");
+ plotter.region(0).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ plotter.region(1).style().statisticsBoxStyle().setVisible(false);
+ plotter.region(1).style().setParameter("hist2DStyle", "colorMap");
+ plotter.region(1).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ plotter.region(2).style().statisticsBoxStyle().setVisible(false);
+ plotter.region(2).style().setParameter("hist2DStyle", "colorMap");
+ plotter.region(2).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+
+ //Extract gain corrections
+
+ plotFrame.pack();
+ plotFrame.setVisible(false);
+
+ //Reset basic histograms
+
+ for(int irow=-5;irow<=5;++irow) {
+ for(int icol=-23;icol<=23;++icol) {
+ pe[icol+23][irow+5].reset();
+ }
+ }
+
+
+ }
+
+
+ public void setAllGain(double g) {
+
+ for(int irow=-5;irow<=5;++irow) {
+ for(int icol=-23;icol<=23;++icol) {
+ this.setGain(icol,irow,g);
+ }
+ }
+ }
+
+ public void setAllGain(String fileName) {
+
+ FileReader fReader = null;
+ BufferedReader bReader = null;
+ try {
+ fReader = new FileReader(fileName);
+ } catch (FileNotFoundException ex) {
+ Logger.getLogger(GainCalib.class.getName()).log(Level.SEVERE, null, ex);
+ }
+ bReader = new BufferedReader(fReader);
+ String line;
+ try {
+ while((line = bReader.readLine()) != null) {
+ if(line.contains("#")) continue;
+ String[] vec = line.split("\\s+");
+ if(vec.length!=4) {
+ System.out.println("Wrong format for line: " + line);
+ throw new RuntimeException("");
+ }
+ int E = Integer.valueOf(vec[0]);
+ int x = Integer.valueOf(vec[1]);
+ int y = Integer.valueOf(vec[2]);
+ double g = Double.valueOf(vec[3]);
+ if(E==0) this.setGain(x,y,g);
+ }
+ } catch (IOException ex) {
+ Logger.getLogger(GainCalib.class.getName()).log(Level.SEVERE, null, ex);
+ }
+
+
+ }
+
+ public void setGain(int icol,int irow,double corr) {
+
+ if(helper==null) {
+ helper = _detector.getSubdetector(_ecalName).getDetectorElement().getIdentifierHelper();
+ expId = new ExpandedIdentifier(helper.getIdentifierDictionary().getNumberOfFields());
+ expId.setValue(helper.getFieldIndex("system"), _detector.getSubdetector(_ecalName).getSystemID());
+ }
+
+
+ //int x = Integer.valueOf(lineTok.nextToken());
+ //int y = Integer.valueOf(lineTok.nextToken());
+ expId.setValue(helper.getFieldIndex("ix"), icol);
+ expId.setValue(helper.getFieldIndex("iy"), irow);
+ long cellid = helper.pack(expId).getValue();
+
+ //if(!this._gains.containsKey(cellid)) {
+ // System.out.println("ERROR updating gain from icol " + icol + " irow " + irow + " this cellid " + cellid + " doesn't have a previous gain?");
+ // System.exit(1);
+ //}
+ double gain = corr;//this._gains.get(cellid);
+ this._gains.put(cellid, corr);
+
+
+ }
+
+ public double getGain(int icol,int irow) {
+ if(helper==null) {
+ helper = _detector.getSubdetector(_ecalName).getDetectorElement().getIdentifierHelper();
+ expId = new ExpandedIdentifier(helper.getIdentifierDictionary().getNumberOfFields());
+ expId.setValue(helper.getFieldIndex("system"), _detector.getSubdetector(_ecalName).getSystemID());
+ }
+ expId.setValue(helper.getFieldIndex("ix"), icol);
+ expId.setValue(helper.getFieldIndex("iy"), irow);
+ long cellid = helper.pack(expId).getValue();
+ if(!this._gains.containsKey(cellid)) {
+ System.out.println("ERROR getting gain from icol " + icol + " irow " + irow + " this cellid " + cellid + " doesn't have a previous gain?");
+ System.exit(1);
+ }
+ return this._gains.get(cellid);
+
+ }
+
+
+ public MatchedTrack getTrack(int evt,long id) {
+ MatchedTrack mt = null;
+
+ if(this._tracks.containsKey(evt)) {
+ for(MatchedTrack trk : this._tracks.get(evt)) {
+ if(trk._cellid == id) {
+ mt = trk;
+ }
+ }
+ }
+ return mt;
+
+ }
+
+ public double getGain(BaseRawCalorimeterHit hit) {
+ long id = hit.getCellID();
+ if(!this._gains.containsKey(id)) {
+ System.out.println("ERROR hit with id " + id + " doesn't exist in gain map!!!");
+ }
+ return this._gains.get(id);
+ }
+
+ public List<HPSEcalCluster> getCalibratedClusters(int eventNumber) {
+
+ if(!_hits.containsKey(eventNumber)) {
+ System.out.println("ERROR no raw hits for event " + eventNumber);
+ System.exit(1);
+ }
+
+
+ if(_hits.get(eventNumber).isEmpty()) {
+ System.out.println("empy hits for " + eventNumber);
+ System.exit(1);
+ }
+
+ if(_debug) System.out.println("Converting " + _hits.get(eventNumber).size() + " raw ecal hits");
+
+
+
+
+ ArrayList<CalorimeterHit> newHits = new ArrayList<CalorimeterHit>();
+
+ for (BaseRawCalorimeterHit hit : _hits.get(eventNumber)) {
+
+ if(_debug) System.out.println("Converting hit id " + hit.getCellID());
+
+ double gain = getGain(hit);
+ CalorimeterHit newHit = _converter.HitDtoA(hit,_integralWindow,gain);
+ if(_applyBadCrystalMap && HPSEcalRawConverterDriver.isBadCrystal(newHit)) continue;
+ //if(applyBadCrystalMap && HPSEcalRawConverterDriver.isBadCrystal(newHit)) continue;
+ //if (dropBadFADC && HPSEcalRawConverterDriver.isBadFADC(newHit)) continue;
+ if (newHit.getRawEnergy() > _threshold) {
+ newHits.add(newHit);
+ }
+
+ if(_debug)
+ System.out.println("Converted hit id " + hit.getCellID() + " with gain " + gain);
+
+
+ }
+
+ if(_debug) System.out.println(newHits.size() + " calibrated calorimeter hits for event " + eventNumber);
+
+ //To be able to use decoders I need to associate with readout. Use event put to do that?! -> FIX THIS
+ int run =999999;
+ BaseLCSimEvent event = new BaseLCSimEvent(999999,eventNumber,_detector.getDetectorName());
+ event.put(_ecalTmpCollectionName, newHits, CalorimeterHit.class, 0, _ecalReadoutName);
+
+
+ //Now we have a set of calibrated ecal hits
+ // Clusterize them
+
+ // Make a hit map for quick lookup by ID.
+ Map<Long, CalorimeterHit> hitMap = new HashMap<Long, CalorimeterHit>();
+ for (CalorimeterHit hit : newHits) {
+ hitMap.put(hit.getCellID(), hit);
+ }
+
+ List<HPSEcalCluster> clusters = _clusterer.createClusters(hitMap);
+
+ if(_debug) System.out.println(clusters.size() + " clusters found from event " + eventNumber);
+
+
+ //clean up for this event!
+ //LCSimEvent should clear by itself!?
+ newHits.clear();
+
+ return clusters;
+
+
+ }
+
+
+
+
+}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1