lcsim-cal-calib/src/org/lcsim/cal/calib
diff -u -r1.10 -r1.11
--- ClusterEnergyAnalysis.java 23 Jul 2008 00:19:33 -0000 1.10
+++ ClusterEnergyAnalysis.java 23 Jul 2008 16:16:37 -0000 1.11
@@ -3,7 +3,7 @@
*
* Created on July 14, 2008, 6:18 PM
*
- * $Id: ClusterEnergyAnalysis.java,v 1.10 2008/07/23 00:19:33 ngraf Exp $
+ * $Id: ClusterEnergyAnalysis.java,v 1.11 2008/07/23 16:16:37 ngraf Exp $
*/
package org.lcsim.cal.calib;
@@ -135,15 +135,19 @@
meV = true;
}
- _tree.mkdirs(particleType);
- _tree.cd(particleType);
- _tree.mkdirs(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
- _tree.cd(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
-
// TODO: make this cluster type selection more robust
String type = "gamma";
if(!particleType.equals("gamma")) type = "neutralHadron";
+
+// _tree.mkdirs(particleType);
+// _tree.cd(particleType);
+ _tree.mkdirs(type);
+ _tree.cd(type);
+ _tree.mkdirs(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
+ _tree.cd(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
+
+
// this analysis is intended for single particle calorimeter response.
// let's make sure that the primary particle did not interact in the
// tracker...
@@ -505,213 +509,211 @@
boolean showPlots = false;
boolean doit = false;
if(doit)
- {
- String fileType = "png";
- String[] pars = {"amplitude", "mean","sigma"};
-
- IAnalysisFactory af = IAnalysisFactory.create();
- String[] dirs = _tree.listObjectNames(".");
- for (int ii=0; ii<dirs.length; ++ii)
{
+ String fileType = "png";
+ String[] pars = {"amplitude", "mean","sigma"};
+
+ IAnalysisFactory af = IAnalysisFactory.create();
+ String[] dirs = _tree.listObjectNames(".");
+ for (int ii=0; ii<dirs.length; ++ii)
+ {
// System.out.println("dirs["+i+"]= "+dirs[i]);
- String[] parts = dirs[ii].split("/");
+ String[] parts = dirs[ii].split("/");
// for(int k=0; k<parts.length; ++k)
// {
// System.out.println("parts["+k+"]= "+parts[k]);
// }
- _tree.cd(dirs[ii]);
- String[] objects = _tree.listObjectNames(".");
-
+ _tree.cd(dirs[ii]);
+ String[] objects = _tree.listObjectNames(".");
+
// for(int j=0; j<objects.length;++j)
// {
// System.out.println("obj["+j+"]= "+objects[i]);
// }
-
- sortDirectoriesByEnergy(objects);
-
- int numberOfPoints = objects.length;
- double[] energies = new double[objects.length];
- for(int j=0; j<objects.length; ++j)
- {
-// System.out.println(objects[j]);
- String subDir =parts[1];
- String[] st = objects[j].split("/")[1].split("_");
- String e = st[0];
- String unit = st[1];
+ sortDirectoriesByEnergy(objects);
+
+ int numberOfPoints = objects.length;
+ double[] energies = new double[objects.length];
+ for(int j=0; j<objects.length; ++j)
+ {
+// System.out.println(objects[j]);
+
+ String subDir =parts[1];
+ String[] st = objects[j].split("/")[1].split("_");
+ String e = st[0];
+ String unit = st[1];
//// System.out.println(e+" "+unit);
- energies[j] = Double.parseDouble(e);
- if(unit.contains("MeV")) energies[j]/=1000.;
+ energies[j] = Double.parseDouble(e);
+ if(unit.contains("MeV")) energies[j]/=1000.;
// System.out.println("energy: "+energies[j]);
- }
- IFunctionFactory functionFactory = af.createFunctionFactory(_tree);
- IFitFactory fitFactory = af.createFitFactory();
- IFitter jminuit = fitFactory.createFitter("Chi2","jminuit");
- IFunction gauss = functionFactory.createFunctionByName("gauss","G");
- IFunction line = functionFactory.createFunctionByName("line","P1");
- IDataPointSetFactory dpsf = af.createDataPointSetFactory(_tree);
-
- IPlotter plotter = af.createPlotterFactory().create("sampling fraction plot");
- plotter.createRegions(3, 4, 0);
- IPlotterStyle style2 = plotter.region(7).style();
- style2.legendBoxStyle().setVisible(false);
- style2.statisticsBoxStyle().setVisible(false);
-
-
- IPlotterStyle style;
-
- double[] fitMeans = new double[numberOfPoints];
- double[] fitSigmas = new double[numberOfPoints];
- IDataPointSet energyMeans = dpsf.create("energy means vs E",2);
- IDataPointSet energySigmas = dpsf.create("sigma \\/ E vs E",2);
- IDataPointSet resolutionFit = dpsf.create("sigma \\/ E vs 1 \\/ \u221a E",2);
- IDataPointSet energyResiduals = dpsf.create("energy residuals (%) vs E",2);
- double eMax = 0;
- for(int i=0; i< numberOfPoints; ++i)
- {
- if(energies[i] > .1) // do not analyze 100MeV and below...
+ }
+ IFunctionFactory functionFactory = af.createFunctionFactory(_tree);
+ IFitFactory fitFactory = af.createFitFactory();
+ IFitter jminuit = fitFactory.createFitter("Chi2","jminuit");
+ IFunction gauss = functionFactory.createFunctionByName("gauss","G");
+ IFunction line = functionFactory.createFunctionByName("line","P1");
+ IDataPointSetFactory dpsf = af.createDataPointSetFactory(_tree);
+
+ IPlotter plotter = af.createPlotterFactory().create("sampling fraction plot");
+ plotter.createRegions(3, 4, 0);
+ IPlotterStyle style2 = plotter.region(7).style();
+ style2.legendBoxStyle().setVisible(false);
+ style2.statisticsBoxStyle().setVisible(false);
+
+
+ IPlotterStyle style;
+
+ double[] fitMeans = new double[numberOfPoints];
+ double[] fitSigmas = new double[numberOfPoints];
+ IDataPointSet energyMeans = dpsf.create("energy means vs E",2);
+ IDataPointSet energySigmas = dpsf.create("sigma \\/ E vs E",2);
+ IDataPointSet resolutionFit = dpsf.create("sigma \\/ E vs 1 \\/ \u221a E",2);
+ IDataPointSet energyResiduals = dpsf.create("energy residuals (%) vs E",2);
+ double eMax = 0;
+ for(int i=0; i< numberOfPoints; ++i)
{
- System.out.println("Energy "+energies[i]);
-
- ICloud1D e = (ICloud1D) _tree.find(objects[i]+"corrected cluster energy for all clusters");
- if(!e.isConverted()) e.convertToHistogram();
- IHistogram1D eHist = e.histogram();
- gauss.setParameter("amplitude",eHist.maxBinHeight());
- gauss.setParameter("mean",eHist.mean());
- gauss.setParameter("sigma",eHist.rms());
- style = plotter.region(i).style();
- style.legendBoxStyle().setVisible(false);
- style.statisticsBoxStyle().setVisible(false);
- double loElimit = energies[i] - .6*sqrt(energies[i]); // expect ~20% resolution, and go out 3 sigma
- double hiElimit = energies[i] + .6*sqrt(energies[i]);;
- plotter.region(i).setXLimits(loElimit, hiElimit);
- plotter.region(i).plot(eHist);
- IFitResult jminuitResult = jminuit.fit(eHist,gauss);
- double[] fitErrors = jminuitResult.errors();
- IFunction fit = jminuitResult.fittedFunction();
- for(int j=0; j<pars.length; ++j)
+ if(energies[i] > .1) // do not analyze 100MeV and below...
{
- System.out.println(" "+pars[j]+": "+ fit.parameter(pars[j])+" +/- "+fitErrors[j]);
- }
- fitMeans[i] = fit.parameter("mean");
- fitSigmas[i] = fit.parameter("sigma");
- plotter.region(i).plot(fit);
+ System.out.println("Energy "+energies[i]);
+
+ ICloud1D e = (ICloud1D) _tree.find(objects[i]+"corrected cluster energy for all clusters");
+ if(!e.isConverted()) e.convertToHistogram();
+ IHistogram1D eHist = e.histogram();
+ gauss.setParameter("amplitude",eHist.maxBinHeight());
+ gauss.setParameter("mean",eHist.mean());
+ gauss.setParameter("sigma",eHist.rms());
+ style = plotter.region(i).style();
+ style.legendBoxStyle().setVisible(false);
+ style.statisticsBoxStyle().setVisible(false);
+ double loElimit = energies[i] - .6*sqrt(energies[i]); // expect ~20% resolution, and go out 3 sigma
+ double hiElimit = energies[i] + .6*sqrt(energies[i]);;
+ plotter.region(i).setXLimits(loElimit, hiElimit);
+ plotter.region(i).plot(eHist);
+ IFitResult jminuitResult = jminuit.fit(eHist,gauss);
+ double[] fitErrors = jminuitResult.errors();
+ IFunction fit = jminuitResult.fittedFunction();
+ for(int j=0; j<pars.length; ++j)
+ {
+ System.out.println(" "+pars[j]+": "+ fit.parameter(pars[j])+" +/- "+fitErrors[j]);
+ }
+ fitMeans[i] = fit.parameter("mean");
+ fitSigmas[i] = fit.parameter("sigma");
+ plotter.region(i).plot(fit);
// plotter.region(7).plot(eHist);
-
- // the means
- IDataPoint point = energyMeans.addPoint();
- point.coordinate(0).setValue(energies[i]);
- point.coordinate(1).setValue(fitMeans[i]);
- point.coordinate(1).setErrorPlus(fitErrors[1]);
- point.coordinate(1).setErrorMinus(fitErrors[1]);
-
- // sigma
- IDataPoint point1 = energySigmas.addPoint();
- point1.coordinate(0).setValue(energies[i]);
- point1.coordinate(1).setValue(fitSigmas[i]/energies[i]);
- point1.coordinate(1).setErrorPlus(fitErrors[2]/energies[i]);
- point1.coordinate(1).setErrorMinus(fitErrors[2]/energies[i]);
-
- // sigma/E vs 1/sqrt(E)
-
- IDataPoint point3 = resolutionFit.addPoint();
- point3.coordinate(0).setValue(1./sqrt(energies[i]));
- point3.coordinate(1).setValue(fitSigmas[i]/energies[i]);
- point3.coordinate(1).setErrorPlus(fitErrors[2]/energies[i]);
- point3.coordinate(1).setErrorMinus(fitErrors[2]/energies[i]);
-
- // residuals
- IDataPoint point2 = energyResiduals.addPoint();
- point2.coordinate(0).setValue(energies[i]);
- point2.coordinate(1).setValue(100.*(fitMeans[i]-energies[i])/energies[i]);
-
- // axis bookeeping...
- if(energies[i] > eMax) eMax = energies[i];
- } // end of 100 MeV cut
- }
-
-
-
- IPlotter results = af.createPlotterFactory().create("linearity");
- style = results.region(0).style();
- style.xAxisStyle().setLabel("MC Energy [GeV]");
- style.yAxisStyle().setLabel("Cluster Energy [GeV]");
- style.titleStyle().setVisible(false);
- style.statisticsBoxStyle().setVisibileStatistics("011");
- style.legendBoxStyle().setVisible(true);
- IFitResult fitLine = jminuit.fit(energyMeans, line);
- System.out.println(" fit status: "+fitLine.fitStatus());
- double eMaxBin = eMax+10.;
- results.region(0).setXLimits(0., eMaxBin);
- results.region(0).setYLimits(0., eMaxBin);
- results.region(0).plot(energyMeans);
- results.region(0).plot(fitLine.fittedFunction());
-
-
- IPlotter resolution = af.createPlotterFactory().create("resolution");
- style = resolution.region(0).style();
- style.xAxisStyle().setLabel("Energy [GeV]");
- style.yAxisStyle().setLabel("sigma/E");
- style.titleStyle().setVisible(false);
- style.statisticsBoxStyle().setVisible(false);
- style.legendBoxStyle().setVisible(false);
- resolution.region(0).setXLimits(0., eMaxBin);
- resolution.region(0).setYLimits(0., .2);
- resolution.region(0).plot(energySigmas);
-
-
- IPlotter resolution2 = af.createPlotterFactory().create("sigma/E vs 1/E");
- style = resolution2.region(0).style();
- style.xAxisStyle().setLabel("1/ \u221a Energy [1/GeV]");
- style.yAxisStyle().setLabel("sigma/E");
+
+ // the means
+ IDataPoint point = energyMeans.addPoint();
+ point.coordinate(0).setValue(energies[i]);
+ point.coordinate(1).setValue(fitMeans[i]);
+ point.coordinate(1).setErrorPlus(fitErrors[1]);
+ point.coordinate(1).setErrorMinus(fitErrors[1]);
+
+ // sigma
+ IDataPoint point1 = energySigmas.addPoint();
+ point1.coordinate(0).setValue(energies[i]);
+ point1.coordinate(1).setValue(fitSigmas[i]/energies[i]);
+ point1.coordinate(1).setErrorPlus(fitErrors[2]/energies[i]);
+ point1.coordinate(1).setErrorMinus(fitErrors[2]/energies[i]);
+
+ // sigma/E vs 1/sqrt(E)
+
+ IDataPoint point3 = resolutionFit.addPoint();
+ point3.coordinate(0).setValue(1./sqrt(energies[i]));
+ point3.coordinate(1).setValue(fitSigmas[i]/energies[i]);
+ point3.coordinate(1).setErrorPlus(fitErrors[2]/energies[i]);
+ point3.coordinate(1).setErrorMinus(fitErrors[2]/energies[i]);
+
+ // residuals
+ IDataPoint point2 = energyResiduals.addPoint();
+ point2.coordinate(0).setValue(energies[i]);
+ point2.coordinate(1).setValue(100.*(fitMeans[i]-energies[i])/energies[i]);
+
+ // axis bookeeping...
+ if(energies[i] > eMax) eMax = energies[i];
+ } // end of 100 MeV cut
+ }
+
+ IPlotter results = af.createPlotterFactory().create("linearity");
+ style = results.region(0).style();
+ style.xAxisStyle().setLabel("MC Energy [GeV]");
+ style.yAxisStyle().setLabel("Cluster Energy [GeV]");
+ style.titleStyle().setVisible(false);
+ style.statisticsBoxStyle().setVisibileStatistics("011");
+ style.legendBoxStyle().setVisible(true);
+ IFitResult fitLine = jminuit.fit(energyMeans, line);
+ System.out.println(" fit status: "+fitLine.fitStatus());
+ double eMaxBin = eMax+10.;
+ results.region(0).setXLimits(0., eMaxBin);
+ results.region(0).setYLimits(0., eMaxBin);
+ results.region(0).plot(energyMeans);
+ results.region(0).plot(fitLine.fittedFunction());
+
+
+ IPlotter resolution = af.createPlotterFactory().create("resolution");
+ style = resolution.region(0).style();
+ style.xAxisStyle().setLabel("Energy [GeV]");
+ style.yAxisStyle().setLabel("sigma/E");
+ style.titleStyle().setVisible(false);
+ style.statisticsBoxStyle().setVisible(false);
+ style.legendBoxStyle().setVisible(false);
+ resolution.region(0).setXLimits(0., eMaxBin);
+ resolution.region(0).setYLimits(0., .2);
+ resolution.region(0).plot(energySigmas);
+
+
+ IPlotter resolution2 = af.createPlotterFactory().create("sigma/E vs 1/E");
+ style = resolution2.region(0).style();
+ style.xAxisStyle().setLabel("1/ \u221a Energy [1/GeV]");
+ style.yAxisStyle().setLabel("sigma/E");
// style.statisticsBoxStyle().setVisibileStatistics("011");
- style.legendBoxStyle().setVisible(false);
- IFitResult resFitLine = jminuit.fit(resolutionFit, line);
- System.out.println(" fit status: "+resFitLine.fitStatus());
+ style.legendBoxStyle().setVisible(false);
+ IFitResult resFitLine = jminuit.fit(resolutionFit, line);
+ System.out.println(" fit status: "+resFitLine.fitStatus());
// resolution2.region(0).setXLimits(0., 1.05);
// resolution2.region(0).setYLimits(0., .2);
- resolution2.region(0).plot(resolutionFit);
- resolution2.region(0).plot(resFitLine.fittedFunction());
-
- IPlotter residuals = af.createPlotterFactory().create("residuals (%)");
- style = residuals.region(0).style();
- style.xAxisStyle().setLabel("Energy [GeV]");
- style.yAxisStyle().setLabel("Residuals [%]");
- style.statisticsBoxStyle().setVisible(false);
- style.titleStyle().setVisible(false);
-
- residuals.region(0).setXLimits(0., eMaxBin);
-
- residuals.region(0).plot(energyResiduals);
-
- if(showPlots)
- {
- plotter.show();
- results.show();
- resolution.show();
- resolution2.show();
- residuals.show();
- }
- else
- {
- try
+ resolution2.region(0).plot(resolutionFit);
+ resolution2.region(0).plot(resFitLine.fittedFunction());
+
+ IPlotter residuals = af.createPlotterFactory().create("residuals (%)");
+ style = residuals.region(0).style();
+ style.xAxisStyle().setLabel("Energy [GeV]");
+ style.yAxisStyle().setLabel("Residuals [%]");
+ style.statisticsBoxStyle().setVisible(false);
+ style.titleStyle().setVisible(false);
+
+ residuals.region(0).setXLimits(0., eMaxBin);
+
+ residuals.region(0).plot(energyResiduals);
+
+ if(showPlots)
{
- // hardcopy
- plotter.writeToFile("energyPlots."+fileType,fileType);
- results.writeToFile("linearity."+fileType,fileType);
- resolution.writeToFile("resolution."+fileType,fileType);
- resolution2.writeToFile("resolutionLinear."+fileType,fileType);
- residuals.writeToFile("residuals."+fileType,fileType);
+ plotter.show();
+ results.show();
+ resolution.show();
+ resolution2.show();
+ residuals.show();
}
- catch(IOException e)
+ else
{
- System.out.println("problem writing out hardcopy in "+fileType+" format");
- e.printStackTrace();
+ try
+ {
+ // hardcopy
+ plotter.writeToFile("energyPlots."+fileType,fileType);
+ results.writeToFile("linearity."+fileType,fileType);
+ resolution.writeToFile("resolution."+fileType,fileType);
+ resolution2.writeToFile("resolutionLinear."+fileType,fileType);
+ residuals.writeToFile("residuals."+fileType,fileType);
+ }
+ catch(IOException e)
+ {
+ System.out.println("problem writing out hardcopy in "+fileType+" format");
+ e.printStackTrace();
+ }
+
}
-
- }
- }// end of loop over directories
- }
+ }// end of loop over directories
+ }
}
private static void sortDirectoriesByEnergy(String[] s)