lcsim-cal-calib/src/org/lcsim/cal/calib
diff -u -r1.3 -r1.4
--- ClusterEnergyAnalysis.java 15 Jul 2008 01:58:25 -0000 1.3
+++ ClusterEnergyAnalysis.java 15 Jul 2008 02:59:50 -0000 1.4
@@ -3,11 +3,24 @@
*
* Created on July 14, 2008, 6:18 PM
*
- * $Id: ClusterEnergyAnalysis.java,v 1.3 2008/07/15 01:58:25 ngraf Exp $
+ * $Id: ClusterEnergyAnalysis.java,v 1.4 2008/07/15 02:59:50 ngraf Exp $
*/
package org.lcsim.cal.calib;
+import hep.aida.IAnalysisFactory;
+import hep.aida.ICloud1D;
+import hep.aida.IDataPoint;
+import hep.aida.IDataPointSet;
+import hep.aida.IDataPointSetFactory;
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+import hep.aida.IHistogram1D;
+import hep.aida.IPlotter;
+import hep.aida.IPlotterStyle;
import hep.aida.ITree;
import hep.physics.vec.Hep3Vector;
import java.util.List;
@@ -30,8 +43,10 @@
import static java.lang.Math.PI;
import static java.lang.Math.abs;
import static java.lang.Math.sqrt;
+import java.util.Arrays;
import java.util.HashMap;
import java.util.Map;
+import java.util.StringTokenizer;
import org.lcsim.event.MCParticle;
/**
@@ -122,6 +137,10 @@
_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";
+
// this analysis is intended for single particle calorimeter response.
// let's make sure that the primary particle did not interact in the
// tracker...
@@ -129,7 +148,6 @@
// this is just crap. Why not use SpacePoint?
double radius = sqrt(endpoint.x()*endpoint.x()+endpoint.y()*endpoint.y());
double z = endpoint.z();
-// System.out.println("Input MCParticle endpoint: r="+radius+" z= "+z);
boolean doit = true;
if(radius<emCalInnerRadius && abs(z) < emCalInnerZ) doit = false;
@@ -139,24 +157,235 @@
List<CalorimeterHit> hitsToCluster = _collectionmanager.getList(processedHitsName);
// cluster the hits
List<Cluster> clusters = _fcc.createClusters(hitsToCluster);
- String type = "gamma";
+
+
+ // simple sanity check to make sure things don't go awry
for(Cluster c : clusters)
{
aida.cloud1D("uncorrected cluster energy for all clusters").fill(c.getEnergy());
double e = correctClusterEnergy(c, type);
aida.cloud1D("corrected cluster energy for all clusters").fill(e);
- SpacePoint p = new CartesianPoint(c.getPosition());
- double clusterTheta = p.theta();
- aida.cloud2D("uncorrected cluster energy for all clusters vs theta").fill(clusterTheta, c.getEnergy());
- aida.cloud2D("corrected cluster energy for all clusters vs theta").fill(clusterTheta, e);
- aida.cloud2D("raw-corrected cluster energy for all clusters vs theta").fill(clusterTheta, c.getEnergy()-e);
}
+
+ // for simplicity proceed only with the highest energy cluster...
+ if(clusters.size()>0)
+ {
+ Cluster c = clusters.get(0);
+ aida.cloud1D("uncorrected cluster energy for highest energy cluster").fill(c.getEnergy());
+ double e = correctClusterEnergy(c, type);
+ aida.cloud1D("corrected cluster energy for highest energy cluster").fill(e);
+// SpacePoint p = new CartesianPoint(c.getPosition());
+// double clusterTheta = p.theta();
+// aida.cloud2D("uncorrected cluster energy for all clusters vs theta").fill(clusterTheta, c.getEnergy());
+// aida.cloud2D("corrected cluster energy for all clusters vs theta").fill(clusterTheta, e);
+// aida.cloud2D("raw-corrected cluster energy for all clusters vs theta").fill(clusterTheta, c.getEnergy()-e);
+ }
+
}// end of check on decays outside tracker volume
//reset our histogram tree
- _tree.cd("/");
+ _tree.cd("/");
}
+
+ protected void endOfData()
+ {
+ boolean showPlots = true;
+ 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)
+ {
+// 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");
+ 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");
+// 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
+ {
+ // hardcopy
+// plotter.writeToFile(fileFullPath+"_energyPlots."+fileType,fileType);
+// results.writeToFile(fileFullPath+"_linearity."+fileType,fileType);
+// resolution.writeToFile(fileFullPath+"_resolution."+fileType,fileType);
+// resolution2.writeToFile(fileFullPath+"_resolutionLinear."+fileType,fileType);
+// residuals.writeToFile(fileFullPath+"_residuals."+fileType,fileType);
+ }
+ }
+
+
+
private double correctClusterEnergy(Cluster c, String type)
{
double e = 0.;
@@ -228,5 +457,37 @@
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]);
+ }
+
+
+
+ }
+
+
}
