1 package org.lcsim.cal.calib;
2
3
4
5
6
7
8
9 import hep.aida.IAnalysisFactory;
10 import hep.aida.ICloud1D;
11 import hep.aida.IDataPoint;
12 import hep.aida.IDataPointSet;
13 import hep.aida.IDataPointSetFactory;
14 import hep.aida.IFitFactory;
15 import hep.aida.IFitResult;
16 import hep.aida.IFitter;
17 import hep.aida.IFunction;
18 import hep.aida.IFunctionFactory;
19 import hep.aida.IHistogram1D;
20 import hep.aida.IPlotter;
21 import hep.aida.IPlotterFactory;
22 import hep.aida.IPlotterStyle;
23 import hep.aida.ITree;
24 import static java.lang.Math.sqrt;
25 import java.util.Arrays;
26 import java.util.HashMap;
27 import java.util.Map;
28
29
30
31
32 public class ClusterEnergyAnalysisHistogramFitter
33 {
34
35 public static void main(String[] args) throws Exception
36 {
37 boolean showPlots = true;
38 String fileType = "png";
39 String[] pars = {"amplitude", "mean","sigma"};
40
41
42
43
44 String fileFullPath = null;
45 if(args.length == 0)
46 {
47 System.out.println("Please provide name of aida file to be analyzed");
48 return;
49 }
50 if(args.length>0) fileFullPath = args[0];
51 if(args.length>1)
52 {
53 fileType = args[1];
54 showPlots = false;
55 }
56 IAnalysisFactory af = IAnalysisFactory.create();
57 ITree tree = af.createTreeFactory().create(fileFullPath,"xml", true, false);
58
59 String[] dirs = tree.listObjectNames(".");
60 for (int ii=0; ii<dirs.length; ++ii)
61 {
62
63 String[] parts = dirs[ii].split("/");
64
65
66
67
68 tree.cd(dirs[ii]);
69 String[] objects = tree.listObjectNames(".");
70
71
72
73
74
75
76 sortDirectoriesByEnergy(objects);
77
78
79 IPlotterFactory pf = af.createPlotterFactory();
80 IPlotterStyle fitStyle = pf.createPlotterStyle();
81 fitStyle.dataStyle().outlineStyle().setVisible(true);
82 fitStyle.dataStyle().outlineStyle().setColor("BLACK");
83 fitStyle.dataStyle().outlineStyle().setThickness(5);
84
85
86 IPlotterStyle dataStyle = af.createPlotterFactory().createPlotterStyle();
87 dataStyle.dataStyle().lineStyle().setVisible(true);
88 dataStyle.dataStyle().lineStyle().setColor("BLACK");
89 dataStyle.dataStyle().fillStyle().setColor("RED");
90 dataStyle.dataStyle().fillStyle().setVisible(true);
91
92
93 IPlotterStyle pointStyle = af.createPlotterFactory().createPlotterStyle();
94 pointStyle.dataStyle().markerStyle().setColor("RED");
95 pointStyle.dataStyle().errorBarStyle().setColor("RED");
96
97 pointStyle.dataStyle().outlineStyle().setColor("BLACK");
98
99
100
101 int numberOfPoints = objects.length;
102 double[] energies = new double[objects.length];
103 for(int j=0; j<objects.length; ++j)
104 {
105
106
107 String subDir =parts[1];
108 String[] st = objects[j].split("/")[1].split("_");
109 String e = st[0];
110 String unit = st[1];
111
112 energies[j] = Double.parseDouble(e);
113 if(unit.contains("MeV")) energies[j]/=1000.;
114
115 }
116 IFunctionFactory functionFactory = af.createFunctionFactory(tree);
117 IFitFactory fitFactory = af.createFitFactory();
118 IFitter jminuit = fitFactory.createFitter("Chi2","jminuit");
119 IFunction gauss = functionFactory.createFunctionByName("gauss","G");
120 IFunction line = functionFactory.createFunctionByName("line","P1");
121 IDataPointSetFactory dpsf = af.createDataPointSetFactory(tree);
122
123 IPlotter plotter = af.createPlotterFactory().create("sampling fraction plot");
124 plotter.createRegions(3, 4, 0);
125 IPlotterStyle style2 = plotter.region(7).style();
126 style2.legendBoxStyle().setVisible(false);
127 style2.statisticsBoxStyle().setVisible(false);
128
129 style2.dataStyle().fillStyle().setColor("RED");
130
131 IPlotterStyle style;
132
133 double[] fitMeans = new double[numberOfPoints];
134 double[] fitSigmas = new double[numberOfPoints];
135 IDataPointSet energyMeans = dpsf.create("energy means vs E",2);
136 IDataPointSet energySigmas = dpsf.create("sigma \\/ E vs E",2);
137 IDataPointSet resolutionFit = dpsf.create("sigma \\/ E vs 1 \\/ \u221a E",2);
138 IDataPointSet energyResiduals = dpsf.create("energy residuals vs E",2);
139 IDataPointSet energyResidualsPercent = dpsf.create("energy residuals (%) vs E",2);
140 double eMax = 0;
141 double eSigma = 0.8;
142 if(dirs[ii].contains("gamma")) eSigma = 0.2;
143 for(int i=0; i< numberOfPoints; ++i)
144 {
145 if(energies[i] > .1)
146 {
147 System.out.println("Energy "+energies[i]);
148
149 ICloud1D e = (ICloud1D) tree.find(objects[i]+"corrected cluster energy for highest energy cluster");
150 if(!e.isConverted()) e.convertToHistogram();
151 IHistogram1D eHist = e.histogram();
152 gauss.setParameter("amplitude",eHist.maxBinHeight());
153 gauss.setParameter("mean",eHist.mean());
154 gauss.setParameter("sigma",eHist.rms());
155 style = plotter.region(i).style();
156 style.legendBoxStyle().setVisible(false);
157 style.statisticsBoxStyle().setVisible(false);
158 style.dataStyle().fillStyle().setColor("RED");
159 double loElimit = energies[i] - 3*eSigma*sqrt(energies[i]);
160 double hiElimit = energies[i] + 3*eSigma*sqrt(energies[i]);;
161 plotter.region(i).setXLimits(loElimit, hiElimit);
162 plotter.region(i).plot(eHist, dataStyle);
163 IFitResult jminuitResult = jminuit.fit(eHist,gauss);
164 double[] fitErrors = jminuitResult.errors();
165 IFunction fit = jminuitResult.fittedFunction();
166 for(int j=0; j<pars.length; ++j)
167 {
168 System.out.println(" "+pars[j]+": "+ fit.parameter(pars[j])+" +/- "+fitErrors[j]);
169 }
170 fitMeans[i] = fit.parameter("mean");
171 fitSigmas[i] = fit.parameter("sigma");
172 plotter.region(i).plot(fit, fitStyle);
173
174
175
176 IDataPoint point = energyMeans.addPoint();
177 point.coordinate(0).setValue(energies[i]);
178 point.coordinate(1).setValue(fitMeans[i]);
179 point.coordinate(1).setErrorPlus(fitErrors[1]);
180 point.coordinate(1).setErrorMinus(fitErrors[1]);
181
182
183 IDataPoint point1 = energySigmas.addPoint();
184 point1.coordinate(0).setValue(energies[i]);
185 point1.coordinate(1).setValue(fitSigmas[i]/energies[i]);
186 point1.coordinate(1).setErrorPlus(fitErrors[2]/energies[i]);
187 point1.coordinate(1).setErrorMinus(fitErrors[2]/energies[i]);
188
189
190
191 IDataPoint point3 = resolutionFit.addPoint();
192 point3.coordinate(0).setValue(1./sqrt(energies[i]));
193 point3.coordinate(1).setValue(fitSigmas[i]/energies[i]);
194 point3.coordinate(1).setErrorPlus(fitErrors[2]/energies[i]);
195 point3.coordinate(1).setErrorMinus(fitErrors[2]/energies[i]);
196
197
198 IDataPoint point2 = energyResidualsPercent.addPoint();
199 point2.coordinate(0).setValue(energies[i]);
200 point2.coordinate(1).setValue(100.*(fitMeans[i]-energies[i])/energies[i]);
201
202
203 IDataPoint point4 = energyResiduals.addPoint();
204 point4.coordinate(0).setValue(energies[i]);
205 point4.coordinate(1).setValue(fitMeans[i]-energies[i]);
206
207
208 if(energies[i] > eMax) eMax = energies[i];
209 }
210 }
211
212
213
214 IPlotter results = af.createPlotterFactory().create("linearity");
215 style = results.region(0).style();
216 style.xAxisStyle().setLabel("MC Energy [GeV]");
217 style.yAxisStyle().setLabel("Cluster Energy [GeV]");
218 style.titleStyle().setVisible(true);
219
220 style.statisticsBoxStyle().setVisible(true);
221 style.legendBoxStyle().setVisible(true);
222 IFitResult fitLine = jminuit.fit(energyMeans, line);
223 System.out.println(" fit status: "+fitLine.fitStatus());
224 double eMaxBin = eMax+10.;
225 results.region(0).setXLimits(0., eMaxBin);
226 results.region(0).setYLimits(0., eMaxBin);
227 results.region(0).plot(energyMeans, pointStyle);
228 results.region(0).plot(fitLine.fittedFunction(),fitStyle);
229
230
231 IPlotter resolution = af.createPlotterFactory().create("resolution");
232 style = resolution.region(0).style();
233 style.xAxisStyle().setLabel("Energy [GeV]");
234 style.yAxisStyle().setLabel("sigma/E");
235 style.titleStyle().setVisible(false);
236 style.statisticsBoxStyle().setVisible(false);
237 style.legendBoxStyle().setVisible(false);
238 resolution.region(0).setXLimits(0., eMaxBin);
239 resolution.region(0).setYLimits(0., eSigma);
240 resolution.region(0).plot(energySigmas, pointStyle);
241
242
243 IPlotter resolution2 = af.createPlotterFactory().create("sigma/E vs 1/E");
244 style = resolution2.region(0).style();
245 style.xAxisStyle().setLabel("1/ \u221a Energy [1/GeV]");
246 style.yAxisStyle().setLabel("sigma/E");
247 style.statisticsBoxStyle().setVisibileStatistics("011");
248 style.statisticsBoxStyle().setVisible(true);
249 style.legendBoxStyle().setVisible(true);
250 IFitResult resFitLine = jminuit.fit(resolutionFit, line);
251 System.out.println(" fit status: "+resFitLine.fitStatus());
252 double[] resFitLinePars = resFitLine.fittedParameters();
253 double[] resFitLineParErrors = resFitLine.errors();
254 String[] resFitLineParNames = resFitLine.fittedParameterNames();
255
256 System.out.println(" Energy Resolution Fit: ");
257 for (int i=0; i< resFitLinePars.length; ++i)
258 {
259 System.out.println(resFitLineParNames[i]+" : "+resFitLinePars[i]+" +/- "+resFitLineParErrors[i]);
260 }
261
262
263 resolution2.region(0).plot(resolutionFit, pointStyle);
264 resolution2.region(0).plot(resFitLine.fittedFunction(),fitStyle);
265
266
267 IPlotter residuals = af.createPlotterFactory().create("residuals");
268 style = residuals.region(0).style();
269 style.xAxisStyle().setLabel("Energy [GeV]");
270 style.yAxisStyle().setLabel("Residuals [GeV]");
271 style.statisticsBoxStyle().setVisible(false);
272 style.titleStyle().setVisible(false);
273
274 residuals.region(0).setXLimits(0., eMaxBin);
275
276 residuals.region(0).plot(energyResiduals, pointStyle);
277
278
279
280
281 IPlotter residualsPercent = af.createPlotterFactory().create("residuals (%)");
282 style = residualsPercent.region(0).style();
283 style.xAxisStyle().setLabel("Energy [GeV]");
284 style.yAxisStyle().setLabel("Residuals [%]");
285 style.statisticsBoxStyle().setVisible(false);
286 style.titleStyle().setVisible(false);
287
288 residualsPercent.region(0).setXLimits(0., eMaxBin);
289
290
291 residualsPercent.region(0).plot(energyResidualsPercent, pointStyle);
292
293 if(showPlots)
294 {
295 plotter.show();
296 results.show();
297 resolution.show();
298 resolution2.show();
299 residuals.show();
300 residualsPercent.style().dataStyle().outlineStyle().setColor("BLACK");
301 residualsPercent.show();
302 }
303 else
304 {
305
306 plotter.writeToFile("energyPlots."+fileType,fileType);
307 results.writeToFile("linearity."+fileType,fileType);
308 resolution.writeToFile("resolution."+fileType,fileType);
309 resolution2.writeToFile("resolutionLinear."+fileType,fileType);
310 residuals.writeToFile("residuals."+fileType,fileType);
311 residualsPercent.writeToFile("residualsPercent."+fileType,fileType);
312 }
313 }
314 }
315
316 private static void sortDirectoriesByEnergy(String[] s)
317 {
318 Map<Double, String> map = new HashMap<Double, String>();
319 double[] energies = new double[s.length];
320 for(int j=0; j<s.length; ++j)
321 {
322
323 String subDir = s[j].split("/")[1];
324
325 String[] st = subDir.split("_");
326 String e = st[0];
327 String unit = st[1];
328
329 energies[j] = Double.parseDouble(e);
330 if(unit.contains("MeV")) energies[j]/=1000.;
331 map.put(energies[j], s[j]);
332
333 }
334 Arrays.sort(energies);
335 for(int j=0; j<s.length; ++j)
336 {
337 s[j] = map.get(energies[j]);
338 }
339
340
341
342
343 }
344
345 }