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