1 package org.lcsim.analysis;
2
3 import hep.aida.IAnalysisFactory;
4 import hep.aida.ICloud1D;
5 import hep.aida.IDataPoint;
6 import hep.aida.IDataPointSet;
7 import hep.aida.IDataPointSetFactory;
8 import hep.aida.IFitFactory;
9 import hep.aida.IFitResult;
10 import hep.aida.IFitter;
11 import hep.aida.IFunction;
12 import hep.aida.IFunctionFactory;
13 import hep.aida.IHistogram1D;
14 import hep.aida.IPlotter;
15 import hep.aida.IPlotterStyle;
16 import hep.aida.ITree;
17 import java.io.IOException;
18 import java.text.DecimalFormat;
19 import java.util.Arrays;
20 import java.util.Calendar;
21 import java.util.Date;
22 import java.util.GregorianCalendar;
23 import java.util.HashMap;
24 import java.util.List;
25 import java.util.Map;
26 import org.lcsim.event.EventHeader;
27 import org.lcsim.event.MCParticle;
28 import org.lcsim.event.ReconstructedParticle;
29 import org.lcsim.util.Driver;
30 import org.lcsim.util.aida.AIDA;
31
32 import static java.lang.Math.sqrt;
33
34
35
36
37
38
39
40 public class SinglePhotonAnalysis extends Driver
41 {
42
43 private double _expectedResolution;
44 private AIDA aida = AIDA.defaultInstance();
45 private ITree _tree = aida.tree();
46 private boolean _showPlots = true;
47 private String _fileType = "png";
48 private boolean _writeOutAidaFile = false;
49 private String _defaultAidaFileName = "test";
50 private String _detectorName;
51 private String _particleType;
52
53 public SinglePhotonAnalysis()
54 {
55 }
56
57 protected void process(EventHeader event)
58 {
59 _detectorName = event.getDetectorName();
60 _tree.mkdirs(_detectorName);
61 _tree.cd(_detectorName);
62
63 List<MCParticle> mcparts = event.getMCParticles();
64 MCParticle mcpart = mcparts.get(mcparts.size() - 1);
65 _particleType = mcpart.getType().getName();
66 double mcEnergy = mcpart.getEnergy();
67 long mcIntegerEnergy = Math.round(mcEnergy);
68 boolean meV = false;
69 if (mcEnergy < .99)
70 {
71 mcIntegerEnergy = Math.round(mcEnergy * 1000);
72 meV = true;
73 }
74
75 _tree.mkdirs(_particleType);
76 _tree.cd(_particleType);
77 _tree.mkdirs(mcIntegerEnergy + (meV ? "_MeV" : "_GeV"));
78 _tree.cd(mcIntegerEnergy + (meV ? "_MeV" : "_GeV"));
79
80
81
82
83 if (mcpart.getSimulatorStatus().isDecayedInCalorimeter())
84 {
85 List<ReconstructedParticle> rpCollection = event.get(ReconstructedParticle.class, "PandoraPFOCollection");
86
87 aida.cloud1D("Number of ReconstructedParticles found").fill(rpCollection.size());
88 for (ReconstructedParticle p : rpCollection)
89 {
90 aida.cloud1D("Cluster Energy", 99999).fill(p.getEnergy());
91 int id = p.getType();
92 aida.cloud1D("Cluster Energy pid= " + id).fill(p.getEnergy());
93 }
94 }
95
96 _tree.cd("/");
97 }
98
99 @Override
100 protected void endOfData()
101 {
102 String AidaFileName = _defaultAidaFileName + "_" + _detectorName + "_" + this._particleType+"_"+this.getClass().getSimpleName() + "_" + date() + ".aida";
103
104 if (_writeOutAidaFile)
105 {
106 System.out.println("writing aida file to " + AidaFileName);
107
108
109 try
110 {
111 AIDA.defaultInstance().saveAs(AidaFileName);
112
113
114 } catch (IOException x)
115 {
116 System.out.println("Experienced an IOException during");
117 x.printStackTrace();
118 return;
119 }
120 }
121 fitPlots();
122 }
123
124 private void fitPlots()
125 {
126 boolean doit = true;
127 if (doit)
128 {
129 _tree.cd(_detectorName);
130 String[] pars =
131 {
132 "amplitude", "mean", "sigma"
133 };
134
135 IAnalysisFactory af = IAnalysisFactory.create();
136 String[] dirs = _tree.listObjectNames(".");
137 for (int ii = 0; ii < dirs.length; ++ii)
138 {
139
140 String[] parts = dirs[ii].split("/");
141
142
143
144
145 _tree.cd(dirs[ii]);
146 String[] objects = _tree.listObjectNames(".");
147
148
149
150
151
152
153 sortDirectoriesByEnergy(objects);
154
155 int numberOfPoints = objects.length;
156 double[] energies = new double[objects.length];
157 for (int j = 0; j < objects.length; ++j)
158 {
159
160
161 String subDir = parts[1];
162 String[] st = objects[j].split("/")[1].split("_");
163 String e = st[0];
164 String unit = st[1];
165
166 energies[j] = Double.parseDouble(e);
167 if (unit.contains("MeV"))
168 energies[j] /= 1000.;
169
170 }
171 IFunctionFactory functionFactory = af.createFunctionFactory(_tree);
172 IFitFactory fitFactory = af.createFitFactory();
173 IFitter jminuit = fitFactory.createFitter("Chi2", "jminuit");
174 IFunction gauss = functionFactory.createFunctionByName("gauss", "G");
175 IFunction line = functionFactory.createFunctionByName("line", "P1");
176 IDataPointSetFactory dpsf = af.createDataPointSetFactory(_tree);
177
178 IPlotter plotter = af.createPlotterFactory().create("sampling fraction plot");
179 plotter.createRegions(3, 4, 0);
180 IPlotterStyle style2 = plotter.region(7).style();
181 style2.legendBoxStyle().setVisible(false);
182 style2.statisticsBoxStyle().setVisible(false);
183
184
185 IPlotterStyle style;
186
187 double[] fitMeans = new double[numberOfPoints];
188 double[] fitSigmas = new double[numberOfPoints];
189 IDataPointSet energyMeans = dpsf.create("energy means vs E", 2);
190 IDataPointSet energySigmas = dpsf.create("sigma \\/ E vs E", 2);
191 IDataPointSet resolutionFit = dpsf.create("sigma \\/ E vs 1 \\/ \u221a E", 2);
192 IDataPointSet energyResiduals = dpsf.create("energy residuals (%) vs E", 2);
193 double eMax = 0;
194 for (int i = 0; i < numberOfPoints; ++i)
195 {
196 if (energies[i] > .1)
197 {
198 System.out.println("Energy " + energies[i]);
199
200 int nSigma = 5;
201 double expectedSigma = _expectedResolution * sqrt(energies[i]);
202 double lowE = energies[i] - (expectedSigma * nSigma);
203 double hiE = energies[i] + (expectedSigma * nSigma);
204
205 ICloud1D e = (ICloud1D) _tree.find(objects[i] + "Cluster Energy");
206 if (!e.isConverted())
207 {
208 int nbins = 100;
209
210 e.convert(nbins, lowE, hiE);
211 }
212 IHistogram1D eHist = e.histogram();
213 gauss.setParameter("amplitude", eHist.maxBinHeight());
214
215
216
217 gauss.setParameter("mean", energies[i]);
218 gauss.setParameter("sigma", expectedSigma);
219 style = plotter.region(i).style();
220 style.legendBoxStyle().setVisible(false);
221 style.statisticsBoxStyle().setVisible(false);
222
223 style.xAxisStyle().setLabel(_particleType+" "+energies[i]+" GeV");
224 style.titleStyle().setVisible(false);
225
226 double loElimit = lowE;
227 double hiElimit = hiE;
228 plotter.region(i).setXLimits(loElimit, hiElimit);
229 plotter.region(i).plot(eHist);
230 IFitResult jminuitResult = jminuit.fit(eHist, gauss);
231 double[] fitErrors = jminuitResult.errors();
232 IFunction fit = jminuitResult.fittedFunction();
233 for (int j = 0; j < pars.length; ++j)
234 {
235 System.out.println(" " + pars[j] + ": " + fit.parameter(pars[j]) + " +/- " + fitErrors[j]);
236 }
237 fitMeans[i] = fit.parameter("mean");
238 fitSigmas[i] = fit.parameter("sigma");
239 plotter.region(i).plot(fit);
240
241
242
243 IDataPoint point = energyMeans.addPoint();
244 point.coordinate(0).setValue(energies[i]);
245 point.coordinate(1).setValue(fitMeans[i]);
246 point.coordinate(1).setErrorPlus(fitErrors[1]);
247 point.coordinate(1).setErrorMinus(fitErrors[1]);
248
249
250 IDataPoint point1 = energySigmas.addPoint();
251 point1.coordinate(0).setValue(energies[i]);
252 point1.coordinate(1).setValue(fitSigmas[i] / energies[i]);
253 point1.coordinate(1).setErrorPlus(fitErrors[2] / energies[i]);
254 point1.coordinate(1).setErrorMinus(fitErrors[2] / energies[i]);
255
256
257
258 IDataPoint point3 = resolutionFit.addPoint();
259 point3.coordinate(0).setValue(1. / sqrt(energies[i]));
260 point3.coordinate(1).setValue(fitSigmas[i] / energies[i]);
261 point3.coordinate(1).setErrorPlus(fitErrors[2] / energies[i]);
262 point3.coordinate(1).setErrorMinus(fitErrors[2] / energies[i]);
263
264
265 IDataPoint point2 = energyResiduals.addPoint();
266 point2.coordinate(0).setValue(energies[i]);
267 point2.coordinate(1).setValue(100. * (fitMeans[i] - energies[i]) / energies[i]);
268
269
270 if (energies[i] > eMax)
271 eMax = energies[i];
272 }
273 }
274
275 IPlotter results = af.createPlotterFactory().create("linearity");
276 style = results.region(0).style();
277 style.xAxisStyle().setLabel("MC Energy [GeV]");
278 style.yAxisStyle().setLabel("Cluster Energy [GeV]");
279 style.titleStyle().setVisible(false);
280 style.statisticsBoxStyle().setVisibileStatistics("011");
281 style.legendBoxStyle().setVisible(true);
282 IFitResult fitLine = jminuit.fit(energyMeans, line);
283 System.out.println("Linearity fit:");
284 printLineFitResult(fitLine);
285 double eMaxBin = eMax + 10.;
286 results.region(0).setXLimits(0., eMaxBin);
287
288 results.region(0).plot(energyMeans);
289 results.region(0).plot(fitLine.fittedFunction());
290
291
292 IPlotter resolution = af.createPlotterFactory().create("resolution");
293 style = resolution.region(0).style();
294 style.xAxisStyle().setLabel("Energy [GeV]");
295 style.yAxisStyle().setLabel("sigma/E");
296 style.titleStyle().setVisible(false);
297 style.statisticsBoxStyle().setVisible(false);
298 style.legendBoxStyle().setVisible(false);
299 resolution.region(0).setXLimits(0., eMaxBin);
300
301 resolution.region(0).plot(energySigmas);
302
303
304 IPlotter resolution2 = af.createPlotterFactory().create("sigma/E vs 1/E");
305 style = resolution2.region(0).style();
306 style.xAxisStyle().setLabel("1/ \u221a Energy [1/GeV]");
307 style.yAxisStyle().setLabel("sigma/E");
308 style.statisticsBoxStyle().setVisibileStatistics("011");
309 style.legendBoxStyle().setVisible(false);
310 IFitResult resFitLine = jminuit.fit(resolutionFit, line);
311 System.out.println("Resolution fit:");
312 printLineFitResult(resFitLine);
313
314
315
316 resolution2.region(0).plot(resolutionFit);
317 resolution2.region(0).plot(resFitLine.fittedFunction());
318
319 IPlotter residuals = af.createPlotterFactory().create("residuals (%)");
320 style = residuals.region(0).style();
321 style.xAxisStyle().setLabel("Energy [GeV]");
322 style.yAxisStyle().setLabel("Residuals [%]");
323 style.statisticsBoxStyle().setVisible(false);
324 style.titleStyle().setVisible(false);
325
326 residuals.region(0).setXLimits(0., eMaxBin);
327
328 residuals.region(0).plot(energyResiduals);
329
330 if (_showPlots)
331 {
332 plotter.show();
333 results.show();
334 resolution.show();
335 resolution2.show();
336 residuals.show();
337 }
338 else
339 {
340 try
341 {
342
343 plotter.writeToFile("energyPlots." + _fileType, _fileType);
344 results.writeToFile("linearity." + _fileType, _fileType);
345 resolution.writeToFile("resolution." + _fileType, _fileType);
346 resolution2.writeToFile("resolutionLinear." + _fileType, _fileType);
347 residuals.writeToFile("residuals." + _fileType, _fileType);
348 } catch (IOException e)
349 {
350 System.out.println("problem writing out hardcopy in " + _fileType + " format");
351 e.printStackTrace();
352 }
353
354 }
355 }
356 }
357 }
358
359 private static void sortDirectoriesByEnergy(String[] s)
360 {
361 Map<Double, String> map = new HashMap<Double, String>();
362 double[] energies = new double[s.length];
363 for (int j = 0; j < s.length; ++j)
364 {
365
366 String subDir = s[j].split("/")[1];
367
368 String[] st = subDir.split("_");
369 String e = st[0];
370 String unit = st[1];
371
372 energies[j] = Double.parseDouble(e);
373 if (unit.contains("MeV"))
374 energies[j] /= 1000.;
375 map.put(energies[j], s[j]);
376
377 }
378 Arrays.sort(energies);
379 for (int j = 0; j < s.length; ++j)
380 {
381 s[j] = map.get(energies[j]);
382 }
383
384
385
386
387 }
388
389 public void setExpectedResolution(double resolution)
390 {
391 System.out.println("setting expected resolution to " + resolution);
392 _expectedResolution = resolution;
393 }
394
395 public void setShowPlots(boolean showPlots)
396 {
397 _showPlots = showPlots;
398 }
399
400 public void setFileType(String fileType)
401 {
402 _fileType = fileType;
403 }
404
405 public void setWriteAidaFile(boolean writeAidaFile)
406 {
407 _writeOutAidaFile = writeAidaFile;
408 }
409
410 public void setAidaFileName(String aidaFileName)
411 {
412 _defaultAidaFileName = aidaFileName;
413 }
414
415 private void printLineFitResult(IFitResult fit)
416 {
417 System.out.println(" fit status: " + fit.fitStatus());
418
419 String[] parNames =
420 {
421 "intercept", "slope "
422 };
423 double[] parVals = fit.fittedParameters();
424 double[] parErrs = fit.errors();
425 for (int i = 0; i < parNames.length; ++i)
426 {
427 System.out.println(parNames[i] + " : " + parVals[i] + " +/- " + parErrs[i]);
428 }
429 }
430
431 private static String date()
432 {
433 Calendar cal = new GregorianCalendar();
434 Date date = new Date();
435 cal.setTime(date);
436 DecimalFormat formatter = new DecimalFormat("00");
437 String day = formatter.format(cal.get(Calendar.DAY_OF_MONTH));
438 String month = formatter.format(cal.get(Calendar.MONTH) + 1);
439 return cal.get(Calendar.YEAR) + month + day;
440 }
441 }
442