Commit in hps-java on MAIN | |||
src/main/java/org/lcsim/hps/recon/ecal/HPSEcalFADCReadoutDriver.java | +21 | -10 | 1.20 -> 1.21 |
/HPSEcalRawConverterDriver.java | +16 | 1.15 -> 1.16 | |
/HPSEcalCTPClusterer.java | +3 | -1 | 1.7 -> 1.8 |
/HPSEcalConverterAtoDDriver.java | -70 | 1.3 removed | |
/HPSEcalConverter.java | -52 | 1.5 removed | |
/GainCalib.java | -508 | 1.2 removed | |
/GainCalibrationDriver.java | -636 | 1.2 removed | |
/HPSEcalConverterDriver.java | -68 | 1.4 removed | |
sandbox/GainCalibrationDriver.java | +636 | added 1.1 | |
/HPSEcalConverterDriver.java | +68 | added 1.1 | |
/HPSEcalConverter.java | +52 | added 1.1 | |
/GainCalib.java | +508 | added 1.1 | |
/HPSEcalConverterAtoDDriver.java | +70 | added 1.1 | |
src/main/resources/org/lcsim/hps/steering/MultScatAna.lcsim | +2 | -2 | 1.7 -> 1.8 |
+1376 | -1347 |
clean up ECal code, add some debug, fix a bug in FADCReadoutDriver
diff -u -r1.20 -r1.21 --- HPSEcalFADCReadoutDriver.java 6 Nov 2012 02:06:13 -0000 1.20 +++ HPSEcalFADCReadoutDriver.java 20 Nov 2012 23:25:09 -0000 1.21 @@ -27,7 +27,7 @@
* Simulates time evolution of preamp output pulse. * * @author Sho Uemura <[log in to unmask]>
- * @version $Id: HPSEcalFADCReadoutDriver.java,v 1.20 2012/11/06 02:06:13 meeg Exp $
+ * @version $Id: HPSEcalFADCReadoutDriver.java,v 1.21 2012/11/20 23:25:09 meeg Exp $
*/ public class HPSEcalFADCReadoutDriver extends HPSEcalReadoutDriver<RawCalorimeterHit> {
@@ -190,7 +190,7 @@
if (debug) { System.out.format("trigger %d, %d: %d\n", cellID, i, pipeline.getValue(numSamplesBefore - i - 1)); }
- sumBefore += pipeline.getValue(numSamplesBefore - i - 1) - pedestal;
+ sumBefore += pipeline.getValue(numSamplesBefore - i - 1);
} sumMap.put(cellID, sumBefore); } else {
@@ -203,9 +203,9 @@
if (debug) { System.out.format("trigger %d, %d: %d\n", cellID, readoutCounter - timeMap.get(cellID) + numSamplesBefore - 1, pipeline.getValue(0)); }
- sumMap.put(cellID, sum + pipeline.getValue(0) - pedestal);
+ sumMap.put(cellID, sum + pipeline.getValue(0));
} else if (timeMap.get(cellID) + delay0 <= readoutCounter) {
-// System.out.printf("sum = %f\n",sum);
+// System.out.printf("sum = %f\n", sum);
outputQueue.add(new HPSFADCCalorimeterHit(cellID, (int) Math.round(sum / scaleFactor), 64 * timeMap.get(cellID),
@@ -227,17 +227,22 @@
} eDepBuffer.step(); }
- while (outputQueue.peek() != null && outputQueue.peek().getTimeStamp() <= readoutCounter - delay0) { - if (outputQueue.peek().getTimeStamp() < readoutCounter - delay0) {
+ while (outputQueue.peek() != null && outputQueue.peek().getTimeStamp() / 64 <= readoutCounter - delay0) { + if (outputQueue.peek().getTimeStamp() / 64 < readoutCounter - delay0) {
System.out.println("Stale hit in output queue"); outputQueue.poll(); } else { buffer.add(outputQueue.poll()); } }
- while (!buffer.isEmpty() && buffer.peek().getTimeStamp() <= readoutCounter - delay0 - coincidenceWindow) {
+ while (!buffer.isEmpty() && buffer.peek().getTimeStamp() / 64 <= readoutCounter - delay0 - coincidenceWindow) {
buffer.remove(); }
+ if (debug) { + for (RawCalorimeterHit hit : buffer) { + System.out.format("new hit: energy %d\n", hit.getAmplitude()); + } + }
hits.addAll(buffer); }
@@ -254,15 +259,21 @@
protected void processTrigger(EventHeader event) { switch (mode) { case EventConstants.ECAL_WINDOW_MODE:
- System.out.println("Reading out ECal in window mode");
+ if (debug) { + System.out.println("Reading out ECal in window mode"); + }
event.put(ecalReadoutCollectionName, readWindow(), RawTrackerHit.class, 0, ecalReadoutName); break; case EventConstants.ECAL_PULSE_MODE:
- System.out.println("Reading out ECal in pulse mode");
+ if (debug) { + System.out.println("Reading out ECal in pulse mode"); + }
event.put(ecalReadoutCollectionName, readPulses(), RawTrackerHit.class, 0, ecalReadoutName); break; case EventConstants.ECAL_PULSE_INTEGRAL_MODE:
- System.out.println("Reading out ECal in integral mode");
+ if (debug) { + System.out.println("Reading out ECal in integral mode"); + }
event.put(ecalReadoutCollectionName, readIntegrals(), RawCalorimeterHit.class, 0, ecalReadoutName); break; }
diff -u -r1.15 -r1.16 --- HPSEcalRawConverterDriver.java 27 Sep 2012 01:08:12 -0000 1.15 +++ HPSEcalRawConverterDriver.java 20 Nov 2012 23:25:09 -0000 1.16 @@ -63,6 +63,10 @@
this.applyBadCrystalMap = apply; }
+ public void setDebug(boolean debug) { + this.debug = debug; + } +
@Override public void startOfData() { if (ecalCollectionName == null) {
@@ -109,6 +113,9 @@
List<RawCalorimeterHit> hits = event.get(RawCalorimeterHit.class, rawCollectionName); for (RawCalorimeterHit hit : hits) {
+ if (debug) { + System.out.format("old hit energy %d\n", hit.getAmplitude()); + }
CalorimeterHit newHit = converter.HitDtoA(hit, integralWindow); if (newHit.getRawEnergy() > threshold) { if (applyBadCrystalMap && isBadCrystal(newHit)) {
@@ -117,6 +124,9 @@
if (dropBadFADC && isBadFADC(newHit)) { continue; }
+ if (debug) { + System.out.format("new hit energy %f\n", newHit.getRawEnergy()); + }
newHits.add(newHit); } }
@@ -129,8 +139,14 @@
List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName); for (CalorimeterHit hit : hits) {
+ if (debug) { + System.out.format("old hit energy %f\n", hit.getRawEnergy()); + }
RawCalorimeterHit newHit = converter.HitAtoD(hit, integralWindow); if (newHit.getAmplitude() > 0) {
+ if (debug) { + System.out.format("new hit energy %d\n", newHit.getAmplitude()); + }
newHits.add(newHit); } }
diff -u -r1.7 -r1.8 --- HPSEcalCTPClusterer.java 27 Aug 2012 21:53:47 -0000 1.7 +++ HPSEcalCTPClusterer.java 20 Nov 2012 23:25:09 -0000 1.8 @@ -29,7 +29,7 @@
* @author Jeremy McCormick <[log in to unmask]> * @author Tim Nelson <[log in to unmask]> * @author Sho Uemura <[log in to unmask]>
- * @version $Id: HPSEcalCTPClusterer.java,v 1.7 2012/08/27 21:53:47 meeg Exp $
+ * @version $Id: HPSEcalCTPClusterer.java,v 1.8 2012/11/20 23:25:09 meeg Exp $
*/ public class HPSEcalCTPClusterer extends Driver {
@@ -171,6 +171,7 @@
hitSums = new HashMap<Long, Double>(); // Loop over ECal hits to compute energy sums for (CalorimeterHit hit : hits) {
+// System.out.format("hit: %f\n",hit.getRawEnergy());
// Make a hit map for quick lookup by ID. hitMap.put(hit.getCellID(), hit);
@@ -287,6 +288,7 @@
cluster.addHit(clusterHit); } clusters.add(cluster);
+// System.out.println(cluster.getEnergy());
} } return clusters;
diff -N HPSEcalConverterAtoDDriver.java --- HPSEcalConverterAtoDDriver.java 27 Aug 2012 21:53:47 -0000 1.3 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,70 +0,0 @@
-package org.lcsim.hps.recon.ecal; - -import java.util.ArrayList; -import java.util.List; -import org.lcsim.event.CalorimeterHit; -import org.lcsim.event.EventHeader; -import org.lcsim.event.RawCalorimeterHit; -import org.lcsim.util.Driver; -import org.lcsim.util.lcio.LCIOConstants; - -/** - * - * @author Sho Uemura <[log in to unmask]> - * @version $Id: HPSEcalConverterAtoDDriver.java,v 1.3 2012/08/27 21:53:47 meeg Exp $ - */ -public class HPSEcalConverterAtoDDriver extends Driver { - - HPSEcalConverter converter = null; - String rawCollectionName = "EcalDigitizedHits"; - String ecalReadoutName = "EcalHits"; - String ecalCollectionName; - int flags; - - public HPSEcalConverterAtoDDriver() { - flags = 0; - flags += 1 << LCIOConstants.CHBIT_LONG; //store position - flags += 1 << LCIOConstants.RCHBIT_ID1; //store cell ID - converter = new HPSEcalConverter(); - } - - public void setPedestal(double pedestal) { - converter.setPedestal(pedestal); - } - - public void setScale(double scale) { - converter.setScale(scale); - } - - public void setEcalCollectionName(String ecalCollectionName) { - this.ecalCollectionName = ecalCollectionName; - } - - public void setRawCollectionName(String rawCollectionName) { - this.rawCollectionName = rawCollectionName; - } - - @Override - public void startOfData() { - if (ecalCollectionName == null) { - throw new RuntimeException("The parameter ecalCollectionName was not set!"); - } - } - - @Override - public void process(EventHeader event) { - // Get the list of ECal hits. - List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName); - if (hits == null) { - throw new RuntimeException("Event is missing ECal hits collection!"); - } - - ArrayList<RawCalorimeterHit> newHits = new ArrayList<RawCalorimeterHit>(); - - for (CalorimeterHit hit : hits) { - newHits.add(converter.HitAtoD(hit)); - } - - event.put(rawCollectionName, newHits, RawCalorimeterHit.class, flags, ecalReadoutName); - } -}
diff -N HPSEcalConverter.java --- HPSEcalConverter.java 27 Aug 2012 21:53:47 -0000 1.5 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,52 +0,0 @@
-package org.lcsim.hps.recon.ecal; - -import org.lcsim.event.CalorimeterHit; -import org.lcsim.event.RawCalorimeterHit; - -/** - * - * @author Sho Uemura <[log in to unmask]> - * @version $Id: HPSEcalConverter.java,v 1.5 2012/08/27 21:53:47 meeg Exp $ - */ -public class HPSEcalConverter { - - double scale = 0.01; - double pedestal = 0.0; - double period = 4.0; - double dt = 0.0; - - public HPSEcalConverter() { - } - - public double getPedestal() { - return pedestal; - } - - public void setPedestal(double pedestal) { - this.pedestal = pedestal; - } - - public double getScale() { - return scale; - } - - public void setScale(double scale) { - this.scale = scale; - } - - public int AtoD(double amplitude, long cellID) { - return (int) Math.floor((amplitude + pedestal) / scale); - } - - public double DtoA(int amplitude, long cellID) { - return scale * (amplitude + 0.5) + pedestal; - } - - public CalorimeterHit HitDtoA(RawCalorimeterHit hit) { - return new HPSRawCalorimeterHit(DtoA(hit.getAmplitude(), hit.getCellID()), period * hit.getTimeStamp() + dt, hit.getCellID(), 0); - } - - public RawCalorimeterHit HitAtoD(CalorimeterHit hit) { - return new HPSFADCCalorimeterHit(hit.getCellID(), AtoD(hit.getRawEnergy(), hit.getCellID()), (int) Math.round(hit.getTime() / period), 0); - } -}
diff -N GainCalib.java --- GainCalib.java 9 Aug 2012 00:51:15 -0000 1.2 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,508 +0,0 @@
-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.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.util.aida.AIDA; - -/** - * - * @author phansson - */ -public class GainCalib implements ActionListener { - - boolean _debug = false; - EventHeader _event = null; - HPSEcalClusterer _clusterer = 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 _ecalTmpCollectionName = "temporaryEcalCalHits"; - 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]); - } - } - - } - } - - public 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 setOriginalHits(int evt, List<BaseRawCalorimeterHit> hits) { - this._hits.put(evt, hits); - for (BaseRawCalorimeterHit hit : hits) { - this._gains.put(hit.getCellID(), HPSEcalConditions.physicalToGain(hit.getCellID())); - } - } - - 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 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() { - for (int irow = -5; irow <= 5; ++irow) { - for (int icol = -23; icol <= 23; ++icol) { - long cellid = HPSEcalConditions.makePhysicalID(icol, irow); - Double gain = HPSEcalConditions.physicalToGain(cellid); - if (gain != null) { - this.setGain(icol, irow, gain); - } else { - this.setGain(icol, irow, 0); - } - } - } - } - - 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) { - long cellid = HPSEcalConditions.makePhysicalID(icol, irow); - - //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) { - long cellid = HPSEcalConditions.makePhysicalID(icol, irow); - 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 = HPSEcalRawConverter.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 - 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); - } - - return clusters; - } -}
diff -N GainCalibrationDriver.java --- GainCalibrationDriver.java 9 Aug 2012 00:51:15 -0000 1.2 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,636 +0,0 @@
-package org.lcsim.hps.recon.ecal; - -import hep.aida.*; -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.event.CalorimeterHit; -import org.lcsim.event.EventHeader; -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.tracking.EcalTrackMatch; -import org.lcsim.util.Driver; -import org.lcsim.util.aida.AIDA; - -/** - * - * @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 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 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."); - } - } - - @Override - 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.setDetector(detector); - gainCalib.setEcalReadoutName(ecalReadoutName); - 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(); - } - - 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); - 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()); - } - - @Override - public void endOfData() { - - redraw(); - - if (!"".equals(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(); - } - } - } - } - } - - @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]); - } - } - - @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()); - - } - } - } - } - - @Override - public void setEventRefreshRate(int eventRefreshRate) { - refreshRate = eventRefreshRate; - } -}
diff -N HPSEcalConverterDriver.java --- HPSEcalConverterDriver.java 27 Aug 2012 21:53:47 -0000 1.4 +++ /dev/null 1 Jan 1970 00:00:00 -0000 @@ -1,68 +0,0 @@
-package org.lcsim.hps.recon.ecal; - -import java.util.ArrayList; -import java.util.List; -import org.lcsim.event.CalorimeterHit; -import org.lcsim.event.EventHeader; -import org.lcsim.util.Driver; -import org.lcsim.util.lcio.LCIOConstants; - -/** - * - * @author Sho Uemura <[log in to unmask]> - * @version $Id: HPSEcalConverterDriver.java,v 1.4 2012/08/27 21:53:47 meeg Exp $ - */ -public class HPSEcalConverterDriver extends Driver { - - HPSEcalConverter converter = null; - String rawCollectionName; - String ecalReadoutName = "EcalHits"; - String ecalCollectionName = "EcalCorrectedHits"; - int flags; - - public HPSEcalConverterDriver() { - flags = 0; - flags += 1 << LCIOConstants.CHBIT_LONG; //store position - flags += 1 << LCIOConstants.RCHBIT_ID1; //store cell ID - converter = new HPSEcalConverter(); - } - - public void setPedestal(double pedestal) { - converter.setPedestal(pedestal); - } - - public void setScale(double scale) { - converter.setScale(scale); - } - - public void setEcalCollectionName(String ecalCollectionName) { - this.ecalCollectionName = ecalCollectionName; - } - - public void setRawCollectionName(String rawCollectionName) { - this.rawCollectionName = rawCollectionName; - } - - @Override - public void startOfData() { - if (ecalCollectionName == null) { - throw new RuntimeException("The parameter ecalCollectionName was not set!"); - } - } - - @Override - public void process(EventHeader event) { - if (event.hasCollection(HPSFADCCalorimeterHit.class, rawCollectionName)) { - // Get the list of ECal hits. - List<HPSFADCCalorimeterHit> hits = event.get(HPSFADCCalorimeterHit.class, rawCollectionName); - - ArrayList<CalorimeterHit> newHits = new ArrayList<CalorimeterHit>(); - - for (HPSFADCCalorimeterHit hit : hits) { - newHits.add(converter.HitDtoA(hit)); - } - - event.put(ecalCollectionName, newHits, CalorimeterHit.class, flags, ecalReadoutName); - } - } -}
diff -N GainCalibrationDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ GainCalibrationDriver.java 20 Nov 2012 23:25:09 -0000 1.1 @@ -0,0 +1,636 @@
+package org.lcsim.hps.recon.ecal; + +import hep.aida.*; +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.event.CalorimeterHit; +import org.lcsim.event.EventHeader; +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.tracking.EcalTrackMatch; +import org.lcsim.util.Driver; +import org.lcsim.util.aida.AIDA; + +/** + * + * @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 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 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."); + } + } + + @Override + 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.setDetector(detector); + gainCalib.setEcalReadoutName(ecalReadoutName); + 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(); + } + + 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); + 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()); + } + + @Override + public void endOfData() { + + redraw(); + + if (!"".equals(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(); + } + } + } + } + } + + @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]); + } + } + + @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()); + + } + } + } + } + + @Override + public void setEventRefreshRate(int eventRefreshRate) { + refreshRate = eventRefreshRate; + } +}
diff -N HPSEcalConverterDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ HPSEcalConverterDriver.java 20 Nov 2012 23:25:09 -0000 1.1 @@ -0,0 +1,68 @@
+package org.lcsim.hps.recon.ecal; + +import java.util.ArrayList; +import java.util.List; +import org.lcsim.event.CalorimeterHit; +import org.lcsim.event.EventHeader; +import org.lcsim.util.Driver; +import org.lcsim.util.lcio.LCIOConstants; + +/** + * + * @author Sho Uemura <[log in to unmask]> + * @version $Id: HPSEcalConverterDriver.java,v 1.1 2012/11/20 23:25:09 meeg Exp $ + */ +public class HPSEcalConverterDriver extends Driver { + + HPSEcalConverter converter = null; + String rawCollectionName; + String ecalReadoutName = "EcalHits"; + String ecalCollectionName = "EcalCorrectedHits"; + int flags; + + public HPSEcalConverterDriver() { + flags = 0; + flags += 1 << LCIOConstants.CHBIT_LONG; //store position + flags += 1 << LCIOConstants.RCHBIT_ID1; //store cell ID + converter = new HPSEcalConverter(); + } + + public void setPedestal(double pedestal) { + converter.setPedestal(pedestal); + } + + public void setScale(double scale) { + converter.setScale(scale); + } + + public void setEcalCollectionName(String ecalCollectionName) { + this.ecalCollectionName = ecalCollectionName; + } + + public void setRawCollectionName(String rawCollectionName) { + this.rawCollectionName = rawCollectionName; + } + + @Override + public void startOfData() { + if (ecalCollectionName == null) { + throw new RuntimeException("The parameter ecalCollectionName was not set!"); + } + } + + @Override + public void process(EventHeader event) { + if (event.hasCollection(HPSFADCCalorimeterHit.class, rawCollectionName)) { + // Get the list of ECal hits. + List<HPSFADCCalorimeterHit> hits = event.get(HPSFADCCalorimeterHit.class, rawCollectionName); + + ArrayList<CalorimeterHit> newHits = new ArrayList<CalorimeterHit>(); + + for (HPSFADCCalorimeterHit hit : hits) { + newHits.add(converter.HitDtoA(hit)); + } + + event.put(ecalCollectionName, newHits, CalorimeterHit.class, flags, ecalReadoutName); + } + } +}
diff -N HPSEcalConverter.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ HPSEcalConverter.java 20 Nov 2012 23:25:09 -0000 1.1 @@ -0,0 +1,52 @@
+package org.lcsim.hps.recon.ecal; + +import org.lcsim.event.CalorimeterHit; +import org.lcsim.event.RawCalorimeterHit; + +/** + * + * @author Sho Uemura <[log in to unmask]> + * @version $Id: HPSEcalConverter.java,v 1.1 2012/11/20 23:25:09 meeg Exp $ + */ +public class HPSEcalConverter { + + double scale = 0.01; + double pedestal = 0.0; + double period = 4.0; + double dt = 0.0; + + public HPSEcalConverter() { + } + + public double getPedestal() { + return pedestal; + } + + public void setPedestal(double pedestal) { + this.pedestal = pedestal; + } + + public double getScale() { + return scale; + } + + public void setScale(double scale) { + this.scale = scale; + } + + public int AtoD(double amplitude, long cellID) { + return (int) Math.floor((amplitude + pedestal) / scale); + } + + public double DtoA(int amplitude, long cellID) { + return scale * (amplitude + 0.5) + pedestal; + } + + public CalorimeterHit HitDtoA(RawCalorimeterHit hit) { + return new HPSRawCalorimeterHit(DtoA(hit.getAmplitude(), hit.getCellID()), period * hit.getTimeStamp() + dt, hit.getCellID(), 0); + } + + public RawCalorimeterHit HitAtoD(CalorimeterHit hit) { + return new HPSFADCCalorimeterHit(hit.getCellID(), AtoD(hit.getRawEnergy(), hit.getCellID()), (int) Math.round(hit.getTime() / period), 0); + } +}
diff -N GainCalib.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ GainCalib.java 20 Nov 2012 23:25:09 -0000 1.1 @@ -0,0 +1,508 @@
+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.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.util.aida.AIDA; + +/** + * + * @author phansson + */ +public class GainCalib implements ActionListener { + + boolean _debug = false; + EventHeader _event = null; + HPSEcalClusterer _clusterer = 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 _ecalTmpCollectionName = "temporaryEcalCalHits"; + 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]); + } + } + + } + } + + public 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 setOriginalHits(int evt, List<BaseRawCalorimeterHit> hits) { + this._hits.put(evt, hits); + for (BaseRawCalorimeterHit hit : hits) { + this._gains.put(hit.getCellID(), HPSEcalConditions.physicalToGain(hit.getCellID())); + } + } + + 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 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() { + for (int irow = -5; irow <= 5; ++irow) { + for (int icol = -23; icol <= 23; ++icol) { + long cellid = HPSEcalConditions.makePhysicalID(icol, irow); + Double gain = HPSEcalConditions.physicalToGain(cellid); + if (gain != null) { + this.setGain(icol, irow, gain); + } else { + this.setGain(icol, irow, 0); + } + } + } + } + + 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) { + long cellid = HPSEcalConditions.makePhysicalID(icol, irow); + + //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) { + long cellid = HPSEcalConditions.makePhysicalID(icol, irow); + 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 = HPSEcalRawConverter.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 + 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); + } + + return clusters; + } +}
diff -N HPSEcalConverterAtoDDriver.java --- /dev/null 1 Jan 1970 00:00:00 -0000 +++ HPSEcalConverterAtoDDriver.java 20 Nov 2012 23:25:09 -0000 1.1 @@ -0,0 +1,70 @@
+package org.lcsim.hps.recon.ecal; + +import java.util.ArrayList; +import java.util.List; +import org.lcsim.event.CalorimeterHit; +import org.lcsim.event.EventHeader; +import org.lcsim.event.RawCalorimeterHit; +import org.lcsim.util.Driver; +import org.lcsim.util.lcio.LCIOConstants; + +/** + * + * @author Sho Uemura <[log in to unmask]> + * @version $Id: HPSEcalConverterAtoDDriver.java,v 1.1 2012/11/20 23:25:09 meeg Exp $ + */ +public class HPSEcalConverterAtoDDriver extends Driver { + + HPSEcalConverter converter = null; + String rawCollectionName = "EcalDigitizedHits"; + String ecalReadoutName = "EcalHits"; + String ecalCollectionName; + int flags; + + public HPSEcalConverterAtoDDriver() { + flags = 0; + flags += 1 << LCIOConstants.CHBIT_LONG; //store position + flags += 1 << LCIOConstants.RCHBIT_ID1; //store cell ID + converter = new HPSEcalConverter(); + } + + public void setPedestal(double pedestal) { + converter.setPedestal(pedestal); + } + + public void setScale(double scale) { + converter.setScale(scale); + } + + public void setEcalCollectionName(String ecalCollectionName) { + this.ecalCollectionName = ecalCollectionName; + } + + public void setRawCollectionName(String rawCollectionName) { + this.rawCollectionName = rawCollectionName; + } + + @Override + public void startOfData() { + if (ecalCollectionName == null) { + throw new RuntimeException("The parameter ecalCollectionName was not set!"); + } + } + + @Override + public void process(EventHeader event) { + // Get the list of ECal hits. + List<CalorimeterHit> hits = event.get(CalorimeterHit.class, ecalCollectionName); + if (hits == null) { + throw new RuntimeException("Event is missing ECal hits collection!"); + } + + ArrayList<RawCalorimeterHit> newHits = new ArrayList<RawCalorimeterHit>(); + + for (CalorimeterHit hit : hits) { + newHits.add(converter.HitAtoD(hit)); + } + + event.put(rawCollectionName, newHits, RawCalorimeterHit.class, flags, ecalReadoutName); + } +}
diff -u -r1.7 -r1.8 --- MultScatAna.lcsim 27 Aug 2012 18:40:10 -0000 1.7 +++ MultScatAna.lcsim 20 Nov 2012 23:25:09 -0000 1.8 @@ -48,14 +48,14 @@
<outputPlotFileName>/Users/phansson/work/HPS/software/reco/run/trigratefile.aida</outputPlotFileName> </driver>
- <driver name="GainCalibrationDriver"
+<!-- <driver name="GainCalibrationDriver"
type="org.lcsim.hps.recon.ecal.GainCalibrationDriver"> <debug>false</debug> <hideFrame>true</hideFrame> <outputPlotFileName>gaincalib.aida</outputPlotFileName> <ecalGainFileName></ecalGainFileName> <ecalGainCalibCorrFileName>gaincalibout.txt</ecalGainCalibCorrFileName>
- </driver>
+ </driver>-->
<driver name="EcalDaqPlots" type="org.lcsim.hps.monitoring.ecal.EcalDaqPlots"> </driver> <driver name="CleanupDriver"
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