lcsim-cal-calib/src/org/lcsim/cal/calib
diff -u -r1.7 -r1.8
--- ClusterEnergyAnalysis.java 15 Jul 2008 03:36:38 -0000 1.7
+++ ClusterEnergyAnalysis.java 22 Jul 2008 23:01:45 -0000 1.8
@@ -3,7 +3,7 @@
*
* Created on July 14, 2008, 6:18 PM
*
- * $Id: ClusterEnergyAnalysis.java,v 1.7 2008/07/15 03:36:38 ngraf Exp $
+ * $Id: ClusterEnergyAnalysis.java,v 1.8 2008/07/22 23:01:45 ngraf Exp $
*/
package org.lcsim.cal.calib;
@@ -119,6 +119,8 @@
_initialized = true;
}
+ if(event.getEventNumber()%1000==0) System.out.println("Event "+event.getEventNumber());
+
//start processing the event.
// organize the histogram tree by species and energy
List<MCParticle> mcparts = event.getMCParticles();
@@ -189,208 +191,210 @@
}
- protected void endOfData()
- {
- boolean showPlots = false;
- String fileType = "png";
- String[] pars = {"amplitude", "mean","sigma"};
- // int[] intEnergy = {1, 2, 5, 10, 20, 50, 100, 500 };
-// String fileFullPath = "C:/orglcsimAnalyses/SamplingFractionAnalysis_gamma_Theta90_acme0605.aida";
-// if(args.length>0) fileFullPath = args[0];
-// if(args.length>1)
+
+
+// protected void endOfData()
+// {
+// boolean showPlots = false;
+// String fileType = "png";
+// String[] pars = {"amplitude", "mean","sigma"};
+// // int[] intEnergy = {1, 2, 5, 10, 20, 50, 100, 500 };
+//// String fileFullPath = "C:/orglcsimAnalyses/SamplingFractionAnalysis_gamma_Theta90_acme0605.aida";
+//// if(args.length>0) fileFullPath = args[0];
+//// if(args.length>1)
+//// {
+//// fileType = args[1];
+//// showPlots = false;
+//// }
+// IAnalysisFactory af = IAnalysisFactory.create();
+//// ITree tree = af.createTreeFactory().create(fileFullPath,"xml", true, false);
+// System.out.println("calling ls on tree");
+// // TODO find out if I can capture output of ls command to get a list
+// // of available directories...
+// _tree.ls("./gamma");
+// String[] onames = _tree.listObjectNames("./gamma");
+//
+// sortDirectoriesByEnergy(onames);
+//
+// int numberOfPoints = onames.length;
+// double[] energies = new double[onames.length];
+// for(int j=0; j<onames.length; ++j)
// {
-// fileType = args[1];
-// showPlots = false;
+//// System.out.println(onames[j]);
+// String subDir = onames[j].substring(8); // length of "./gamma/" is 8, so remove leading directory
+// StringTokenizer st = new StringTokenizer(subDir,"_");
+// String e = st.nextToken();
+// String unit = st.nextToken();
+//// System.out.println(e+" "+unit);
+// energies[j] = Double.parseDouble(e);
+// if(unit.contains("MeV")) energies[j]/=1000.;
+//// System.out.println("energy: "+energies[j]);
// }
- IAnalysisFactory af = IAnalysisFactory.create();
-// ITree tree = af.createTreeFactory().create(fileFullPath,"xml", true, false);
- System.out.println("calling ls on tree");
- // TODO find out if I can capture output of ls command to get a list
- // of available directories...
- _tree.ls("./gamma");
- String[] onames = _tree.listObjectNames("./gamma");
-
- sortDirectoriesByEnergy(onames);
-
- int numberOfPoints = onames.length;
- double[] energies = new double[onames.length];
- for(int j=0; j<onames.length; ++j)
- {
-// System.out.println(onames[j]);
- String subDir = onames[j].substring(8); // length of "./gamma/" is 8, so remove leading directory
- StringTokenizer st = new StringTokenizer(subDir,"_");
- String e = st.nextToken();
- String unit = st.nextToken();
-// System.out.println(e+" "+unit);
- 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...
- {
- System.out.println("Energy "+energies[i]);
-
- ICloud1D e = (ICloud1D) _tree.find(onames[i]+"corrected cluster energy for highest energy cluster");
- 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(false);
- 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");
+// 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...
+// {
+// System.out.println("Energy "+energies[i]);
+//
+// ICloud1D e = (ICloud1D) _tree.find(onames[i]+"corrected cluster energy for highest energy cluster");
+// 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(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
- {
- 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();
- }
- }
- }
+// style.legendBoxStyle().setVisible(false);
+// 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());
+//// 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
+// {
+// 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();
+// }
+// }
+// }
@@ -465,17 +469,259 @@
return e;
}
+// private static void sortDirectoriesByEnergy(String[] s)
+// {
+// Map<Double, String> map = new HashMap<Double, String>();
+// double[] energies = new double[s.length];
+// for(int j=0; j<s.length; ++j)
+// {
+//// System.out.println(onames[j]);
+// String subDir = s[j].substring(8); // length of "./gamma/" is 8, so remove leading directory
+// StringTokenizer st = new StringTokenizer(subDir,"_");
+// String e = st.nextToken();
+// String unit = st.nextToken();
+//// System.out.println(e+" "+unit);
+// energies[j] = Double.parseDouble(e);
+// if(unit.contains("MeV")) energies[j]/=1000.;
+// map.put(energies[j], s[j]);
+//// System.out.println("energy: "+energies[j]);
+// }
+// Arrays.sort(energies);
+// for(int j=0; j<s.length; ++j)
+// {
+// s[j] = map.get(energies[j]);
+// }
+// for(int j=0; j<s.length; ++j)
+// {
+// System.out.println(s[j]);
+// }
+//
+//
+//
+// }
+
+ protected void endOfData()
+ {
+ boolean showPlots = false;
+ 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("/");
+// for(int k=0; k<parts.length; ++k)
+// {
+// System.out.println("parts["+k+"]= "+parts[k]);
+// }
+ _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];
+//// System.out.println(e+" "+unit);
+ 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...
+ {
+ System.out.println("Energy "+energies[i]);
+
+ ICloud1D e = (ICloud1D) _tree.find(objects[i]+"corrected cluster energy");
+ 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");
+// style.statisticsBoxStyle().setVisibileStatistics("011");
+ 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
+ {
+ // 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
+ }
+
private static void sortDirectoriesByEnergy(String[] s)
{
Map<Double, String> map = new HashMap<Double, String>();
double[] energies = new double[s.length];
for(int j=0; j<s.length; ++j)
{
-// System.out.println(onames[j]);
- String subDir = s[j].substring(8); // length of "./gamma/" is 8, so remove leading directory
- StringTokenizer st = new StringTokenizer(subDir,"_");
- String e = st.nextToken();
- String unit = st.nextToken();
+// System.out.println(s[j]);
+ String subDir = s[j].split("/")[1]; // first token should be "."
+// System.out.println(subDir);
+ String[] st = subDir.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.;
@@ -487,15 +733,10 @@
{
s[j] = map.get(energies[j]);
}
- for(int j=0; j<s.length; ++j)
- {
- System.out.println(s[j]);
- }
-
-
-
+// for(int j=0; j<s.length; ++j)
+// {
+// System.out.println(s[j]);
+// }
}
-
-
}