Commit in hps-java/src/main/java/org/lcsim/hps on MAIN | |||
users/phansson/ecalGainAna.java | +106 | -12 | 1.6 -> 1.7 |
/ECalGainDriver.java | +15 | -2 | 1.8 -> 1.9 |
recon/ecal/HPSEcalEdepToTriggerConverterDriver.java | +5 | -4 | 1.3 -> 1.4 |
+126 | -18 |
ECal calibration with crystal weights
diff -u -r1.6 -r1.7 --- ecalGainAna.java 29 Aug 2012 20:49:23 -0000 1.6 +++ ecalGainAna.java 31 Aug 2012 01:26:33 -0000 1.7 @@ -33,6 +33,7 @@
option.addOption("f", true, "Input file"); option.addOption("s", false, "Save to file"); option.addOption("n", true, "Name added to plots");
+ option.addOption("m", true, "Minimum weight for calibration");
return option; }
@@ -83,6 +84,11 @@
System.exit(1); }
+ double minCount = 20.0; + if (cmd.hasOption("m")) { + minCount = Double.valueOf(cmd.getOptionValue("m")); + } +
ITree sim = null; ITree real = null; try {
@@ -109,6 +115,11 @@
} }
+ IHistogram2D weightPlotSim = (IHistogram2D) sim.find("Weights for correction factor"); + IHistogram2D sumPlotSim = (IHistogram2D) sim.find("Sum for correction factor"); + IHistogram2D weightPlotReal = (IHistogram2D) real.find("Weights for correction factor"); + IHistogram2D sumPlotReal = (IHistogram2D) real.find("Sum for correction factor"); +
// IHistogram2D mpePlotSim = (IHistogram2D) sim.find("<E over p>"); // IHistogram2D spePlotSim = (IHistogram2D) sim.find("RMS(E over p)"); // IHistogram2D mpePlotReal = (IHistogram2D) real.find("<E over p>");
@@ -119,19 +130,90 @@
IHistogramFactory hf = af.createHistogramFactory(null);
+ IHistogram2D corrPlotSim = hf.createHistogram2D("Correction factor - simulation", 47, -23.5, 23.5, 11, -5.5, 5.5); + IHistogram2D corrPlotReal = hf.createHistogram2D("Correction factor - data", 47, -23.5, 23.5, 11, -5.5, 5.5); + + for (int x = -23; x <= 23; x++) { + for (int y = -5; y <= 5; y++) { +// if (weightPlotSim.binHeight(x + 23, y + 5) > minCount) { + corrPlotSim.fill(x, y, sumPlotSim.binHeight(x + 23, y + 5) / weightPlotSim.binHeight(x + 23, y + 5)); +// } +// if (weightPlotReal.binHeight(x + 23, y + 5) > minCount) { + corrPlotReal.fill(x, y, sumPlotReal.binHeight(x + 23, y + 5) / weightPlotReal.binHeight(x + 23, y + 5)); +// } + } + } + +
IPlotter plotter = af.createPlotterFactory().create(); plotter.setParameter("plotterWidth", "600"); plotter.setParameter("plotterHeight", "400"); plotter.createRegion();
- IHistogram2D ecalPlot = hf.createHistogram2D("ECal map", 47, -23.5, 23.5, 11, -5.5, 5.5); - plotter.region(0).plot(ecalPlot);
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(0).plot(weightPlotSim); + plotter.region(0).setTitle(weightPlotSim.title() + " - simulation"); + try { + plotter.writeToFile(outName + "_weights_sim.png", "png"); + } catch (IOException e) { + throw new RuntimeException(e); + } + plotter.region(0).remove(weightPlotSim); + plotter.region(0).plot(weightPlotReal); + plotter.region(0).setTitle(weightPlotReal.title() + " - data"); + plotter.region(0).refresh(); + try { + plotter.writeToFile(outName + "_weights_data.png", "png"); + } catch (IOException e) { + throw new RuntimeException(e); + } + plotter.region(0).remove(weightPlotReal); + +// plotter.region(0).style().zAxisStyle().setParameter("scale", "log"); + plotter.region(0).plot(corrPlotSim); +// plotter.region(0).setZLimits(0, 10.0); +plotter.region(0).style().zAxisStyle().setParameter("upperLimit", "10.0"); + plotter.region(0).setTitle(corrPlotSim.title()); + plotter.region(0).refresh(); + plotter.region(0).style().zAxisStyle().setParameter("upperLimit", "10.0"); + try { + plotter.writeToFile(outName + "_corr_sim.png", "png"); + } catch (IOException e) { + throw new RuntimeException(e); + } + plotter.region(0).remove(corrPlotSim); + plotter.region(0).plot(corrPlotReal); +// plotter.region(0).setZLimits(0, 10.0); +//for (String opt :plotter.region(0).style().zAxisStyle().availableParameters()) { +// System.out.println(opt); +// } + plotter.region(0).setTitle(corrPlotReal.title()); + plotter.region(0).refresh(); + try { + plotter.writeToFile(outName + "_corr_data.png", "png"); + } catch (IOException e) { + throw new RuntimeException(e); + } + plotter.region(0).remove(corrPlotReal); +// plotter.region(0).style().zAxisStyle().setParameter("scale", "lin"); +// plotter.region(0).setZLimits(); + + + + IHistogram2D ecalPlot = hf.createHistogram2D("ECal map", 47, -23.5, 23.5, 11, -5.5, 5.5); + plotter.region(0).plot(ecalPlot); + +
IHistogram1D pePlotsHVSim[][] = new IHistogram1D[4][12]; IHistogram1D pePlotsHVReal[][] = new IHistogram1D[4][12];
+ double[][] weightsHVSim = new double[4][12]; + double[][] sumsHVSim = new double[4][12]; + double[][] weightsHVReal = new double[4][12]; + double[][] sumsHVReal = new double[4][12]; +
for (int quad = 1; quad <= 4; ++quad) { for (int module = 1; module <= 12; ++module) { pePlotsHVSim[quad - 1][module - 1] = hf.createHistogram1D("E over p quadrant=" + quad + " HV group=" + module, 50, 0, 2);
@@ -142,6 +224,10 @@
for (int y = -5; y <= 5; y++) { pePlotsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1].add(pePlotsSim[x + 23][y + 5][0]); pePlotsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1].add(pePlotsReal[x + 23][y + 5][0]);
+ weightsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1] += weightPlotSim.binHeight(x + 23, y + 5); + sumsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1] += sumPlotSim.binHeight(x + 23, y + 5); + weightsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1] += weightPlotReal.binHeight(x + 23, y + 5); + sumsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1] += sumPlotReal.binHeight(x + 23, y + 5);
} }
@@ -158,7 +244,7 @@
plotter.region(0).setTitle("Crystal gains before iterating"); plotter.region(0).refresh(); try {
- plotter.writeToFile(outName + "_initial_gains.png", "png");
+ plotter.writeToFile(outName + "_gains_initial.png", "png");
} catch (IOException e) { throw new RuntimeException(e); }
@@ -173,7 +259,7 @@
plotter.region(0).setTitle("Matched cluster counts in simulation"); plotter.region(0).refresh(); try {
- plotter.writeToFile(outName + "_sim_counts.png", "png");
+ plotter.writeToFile(outName + "_counts_sim.png", "png");
} catch (IOException e) { throw new RuntimeException(e); }
@@ -188,7 +274,7 @@
plotter.region(0).setTitle("Matched cluster counts in data"); plotter.region(0).refresh(); try {
- plotter.writeToFile(outName + "_data_counts.png", "png");
+ plotter.writeToFile(outName + "_counts_data.png", "png");
} catch (IOException e) { throw new RuntimeException(e); }
@@ -201,8 +287,7 @@
throw new RuntimeException(e); }
- int minCount = 50; -
+
gainStream.format("#x\ty\tgain\n"); for (int x = -23; x <= 23; x++) { for (int y = -5; y <= 5; y++) {
@@ -211,12 +296,21 @@
IHistogram1D peReal = pePlotsReal[x + 23][y + 5][0]; IHistogram1D peHVSim = pePlotsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1]; IHistogram1D peHVReal = pePlotsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1];
+ double weightHVReal = weightsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1]; + double sumHVReal = sumsHVReal[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1]; + double weightHVSim = weightsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1]; + double sumHVSim = sumsHVSim[HPSECalUtils.getQuadrant(x, y) - 1][HPSECalUtils.getHVGroup(x, y) - 1];
if (gain != null) {
- if (peSim.allEntries() >= minCount && peReal.allEntries() >= minCount) { - gain *= (peSim.mean() / peReal.mean()); - } else if (peHVSim.allEntries() >= minCount && peHVReal.allEntries() >= minCount) { - gain *= (peHVSim.mean() / peHVReal.mean());
+ if (weightPlotSim.binHeight(x + 23, y + 5) > minCount && weightPlotReal.binHeight(x + 23, y + 5) > minCount) { + gain *= (corrPlotReal.binHeight(x + 23, y + 5) / corrPlotSim.binHeight(x + 23, y + 5)); + } else if (weightHVReal > minCount && weightHVSim > minCount) { +// gain *= (sumHVReal / weightHVReal) / (sumHVSim / weightHVSim);
}
+// if (peSim.allEntries() >= minCount && peReal.allEntries() >= minCount) { +// gain *= (peSim.mean() / peReal.mean()); +// } else if (peHVSim.allEntries() >= minCount && peHVReal.allEntries() >= minCount) { +//// gain *= (peHVSim.mean() / peHVReal.mean()); +// }
gainStream.format("%d\t%d\t%f\n", x, y, gain); ecalPlot.fill(x, y, gain); // System.out.printf("%d\t%d\t%d\t%f\t%f\n", x, y, pePlotsSim[x + 23][y + 5][0].allEntries(), pePlotsSim[x + 23][y + 5][0].mean(), pePlotsSim[x + 23][y + 5][0].rms());
@@ -228,7 +322,7 @@
plotter.region(0).setTitle("Crystal gains after iterating"); plotter.region(0).refresh(); try {
- plotter.writeToFile(outName + "_final_gains.png", "png");
+ plotter.writeToFile(outName + "_gains_final.png", "png");
} catch (IOException e) { throw new RuntimeException(e); }
diff -u -r1.8 -r1.9 --- ECalGainDriver.java 29 Aug 2012 20:49:23 -0000 1.8 +++ ECalGainDriver.java 31 Aug 2012 01:26:33 -0000 1.9 @@ -62,6 +62,8 @@
private IHistogram2D hitmap; private IHistogram1D[] h_PE_t = new IHistogram1D[5]; private IHistogram1D[] h_PE_b = new IHistogram1D[5];
+ private IHistogram2D weightPlot; + private IHistogram2D sumPlot;
int clTop = 0; int clBot = 0; int trTop = 0;
@@ -136,6 +138,9 @@
plotter.region(1).style().setParameter("hist2DStyle", "colorMap"); plotter.region(1).style().dataStyle().fillStyle().setParameter("colorMapScheme", "rainbow");
+ weightPlot = aida.histogram2D("Weights for correction factor", 47, -23.5, 23.5, 11, -5.5, 5.5); + sumPlot = aida.histogram2D("Sum for correction factor", 47, -23.5, 23.5, 11, -5.5, 5.5); +
for (int iE = 0; iE <= 4; ++iE) { for (int irow = -5; irow <= 5; ++irow) { for (int icol = -23; icol <= 23; ++icol) {
@@ -373,8 +378,16 @@
double P = bestTrk.getPX() * 1000; double E = bestCl.getEnergy();
- double Ep = E; - double Eoverp = Ep / P;
+ double Eoverp = E / P; + double PoverE = P / E; + + for (CalorimeterHit hit : bestCl.getCalorimeterHits()) { + double weight = hit.getRawEnergy() / bestCl.getEnergy(); + int ix = hit.getIdentifierFieldValue("ix"); + int iy = hit.getIdentifierFieldValue("iy"); + weightPlot.fill(ix, iy, weight); + sumPlot.fill(ix, iy, weight * PoverE); + }
if (debug) { System.out.println("P " + P + " E " + E);
diff -u -r1.3 -r1.4 --- HPSEcalEdepToTriggerConverterDriver.java 27 Aug 2012 21:53:47 -0000 1.3 +++ HPSEcalEdepToTriggerConverterDriver.java 31 Aug 2012 01:26:33 -0000 1.4 @@ -27,6 +27,7 @@
private double pulseIntegral = tp * Math.E / readoutPeriod; //normalization constant from cal gain (MeV/integral bit) to amplitude gain (amplitude bit/GeV) private double gainNorm = 1000.0 / pulseIntegral;
+ private double gainScale = 1.0;
public HPSEcalEdepToTriggerConverterDriver() { }
@@ -99,7 +100,7 @@
public CalorimeterHit makeTriggerHit(CalorimeterHit hit) { long id = hit.getCellID();
- double amplitude = hit.getRawEnergy() * gainNorm / HPSEcalConditions.physicalToGain(id);
+ double amplitude = hit.getRawEnergy() * gainNorm / HPSEcalConditions.physicalToGain(id) * gainScale;
// double time = readoutPeriod * (Math.random() - 1); double time = 0 - hit.getTime();
@@ -134,12 +135,12 @@
public CalorimeterHit makeReadoutHit(CalorimeterHit hit) { long id = hit.getCellID();
- double amplitude = hit.getRawEnergy() * gainNorm / HPSEcalConditions.physicalToGain(id);
+ double amplitude = hit.getRawEnergy() * gainNorm / HPSEcalConditions.physicalToGain(id) * gainScale;
if (amplitude < readoutThreshold) { return null; }
- double integral = hit.getRawEnergy() * 1000.0; -
+ double integral = hit.getRawEnergy() * 1000.0 * gainScale; +
// double thresholdCrossingTime = 0 - hit.getTime(); // while (true) { // double currentValue = amplitude * pulseAmplitude(thresholdCrossingTime);
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