1 /*
2 * ClusterEnergyAnalysis.java
3 *
4 * Created on July 14, 2008, 6:18 PM
5 *
6 * $Id: ClusterEnergyAnalysis.java,v 1.11 2008/07/23 16:16:37 ngraf Exp $
7 */
8
9 package org.lcsim.cal.calib;
10
11 import hep.aida.IAnalysisFactory;
12 import hep.aida.ICloud1D;
13 import hep.aida.IDataPoint;
14 import hep.aida.IDataPointSet;
15 import hep.aida.IDataPointSetFactory;
16 import hep.aida.IFitFactory;
17 import hep.aida.IFitResult;
18 import hep.aida.IFitter;
19 import hep.aida.IFunction;
20 import hep.aida.IFunctionFactory;
21 import hep.aida.IHistogram1D;
22 import hep.aida.IPlotter;
23 import hep.aida.IPlotterStyle;
24 import hep.aida.ITree;
25 import hep.physics.vec.Hep3Vector;
26 import java.io.IOException;
27 import java.util.List;
28 import org.lcsim.conditions.ConditionsManager;
29 import org.lcsim.conditions.ConditionsManager.ConditionsSetNotFoundException;
30 import org.lcsim.conditions.ConditionsSet;
31 import org.lcsim.event.CalorimeterHit;
32 import org.lcsim.event.Cluster;
33 import org.lcsim.event.EventHeader;
34 import org.lcsim.geometry.IDDecoder;
35 import org.lcsim.geometry.Subdetector;
36 import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer;
37 import org.lcsim.recon.cluster.fixedcone.FixedConeClusterer.FixedConeDistanceMetric;
38 import org.lcsim.spacegeom.CartesianPoint;
39 import org.lcsim.spacegeom.SpacePoint;
40 import org.lcsim.util.Driver;
41 import org.lcsim.util.aida.AIDA;
42
43 import static java.lang.Math.sin;
44 import static java.lang.Math.PI;
45 import static java.lang.Math.abs;
46 import static java.lang.Math.sqrt;
47 import java.util.Arrays;
48 import java.util.HashMap;
49 import java.util.Map;
50 import java.util.StringTokenizer;
51 import org.lcsim.event.MCParticle;
52
53 /**
54 *
55 * @author Norman Graf
56 */
57 public class ClusterEnergyAnalysis extends Driver
58 {
59 private ConditionsSet _cond;
60 private FixedConeClusterer _fcc;
61 private CollectionManager _collectionmanager = CollectionManager.defaultInstance();
62
63 private double[] _ecalLayering;
64 boolean _useFirstLayer;
65
66 private double emCalInnerRadius = 0.;
67 private double emCalInnerZ = 0.;
68
69 private Map<String, Double> _fitParameters = new HashMap<String, Double>();
70
71
72 private AIDA aida = AIDA.defaultInstance();
73 private ITree _tree;
74
75 private boolean _initialized;
76 /** Creates a new instance of ClusterEnergyAnalysis */
77 public ClusterEnergyAnalysis()
78 {
79 _tree = aida.tree();
80 }
81
82 protected void process(EventHeader event)
83 {
84 if(!_initialized)
85 {
86 ConditionsManager mgr = ConditionsManager.defaultInstance();
87 try
88 {
89 _cond = mgr.getConditions("CalorimeterCalibration");
90 }
91 catch(ConditionsSetNotFoundException e)
92 {
93 System.out.println("ConditionSet CalorimeterCalibration not found for detector "+mgr.getDetector());
94 System.out.println("Please check that this properties file exists for this detector ");
95 }
96 double radius = .5;
97 double seed = 0.;//.1;
98 double minE = .05; //.25;
99 _fcc = new FixedConeClusterer(radius, seed, minE, FixedConeDistanceMetric.DPHIDTHETA);
100
101 _ecalLayering = _cond.getDoubleArray("ECalLayering");
102 _useFirstLayer = _cond.getDouble("IsFirstEmLayerSampling")==1.;
103
104 // photons
105 String photonFitParametersList = _cond.getString("PhotonFitParameters");
106 String[] photonFitParameters = photonFitParametersList.split(",\\s");
107 for(int i=0; i<photonFitParameters.length; ++i)
108 {
109 _fitParameters.put(photonFitParameters[i], _cond.getDouble(photonFitParameters[i]));
110 }
111 // neutral hadrons
112 String hadronFitParametersList = _cond.getString("NeutralHadronFitParameters");
113 String[] hadronFitParameters = hadronFitParametersList.split(",\\s");
114 for(int i=0; i<hadronFitParameters.length; ++i)
115 {
116 _fitParameters.put(hadronFitParameters[i], _cond.getDouble(hadronFitParameters[i]));
117 }
118
119 _initialized = true;
120 }
121
122 if(event.getEventNumber()%1000==0) System.out.println("Event "+event.getEventNumber());
123
124 //start processing the event.
125 // organize the histogram tree by species and energy
126 List<MCParticle> mcparts = event.getMCParticles();
127 MCParticle mcpart = mcparts.get(mcparts.size()-1);
128 String particleType = mcpart.getType().getName();
129 double mcEnergy = mcpart.getEnergy();
130 long mcIntegerEnergy = Math.round(mcEnergy);
131 boolean meV = false;
132 if(mcEnergy<.99)
133 {
134 mcIntegerEnergy = Math.round(mcEnergy*1000);
135 meV = true;
136 }
137
138 // TODO: make this cluster type selection more robust
139 String type = "gamma";
140 if(!particleType.equals("gamma")) type = "neutralHadron";
141
142
143 // _tree.mkdirs(particleType);
144 // _tree.cd(particleType);
145 _tree.mkdirs(type);
146 _tree.cd(type);
147 _tree.mkdirs(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
148 _tree.cd(mcIntegerEnergy+(meV ? "_MeV": "_GeV"));
149
150
151 // this analysis is intended for single particle calorimeter response.
152 // let's make sure that the primary particle did not interact in the
153 // tracker...
154 Hep3Vector endpoint = mcpart.getEndPoint();
155 // this is just crap. Why not use SpacePoint?
156 double radius = sqrt(endpoint.x()*endpoint.x()+endpoint.y()*endpoint.y());
157 double z = endpoint.z();
158
159 boolean doit = true;
160 if(radius<emCalInnerRadius && abs(z) < emCalInnerZ) doit = false;
161 if(doit)
162 {
163 String processedHitsName = _cond.getString("ProcessedHitsCollectionName");
164 List<CalorimeterHit> hitsToCluster = _collectionmanager.getList(processedHitsName);
165 // cluster the hits
166 List<Cluster> clusters = _fcc.createClusters(hitsToCluster);
167
168
169 // simple sanity check to make sure things don't go awry
170 for(Cluster c : clusters)
171 {
172 aida.cloud1D("uncorrected cluster energy for all clusters").fill(c.getEnergy());
173 double e = correctClusterEnergy(c, type);
174 aida.cloud1D("corrected cluster energy for all clusters").fill(e);
175 }
176
177 // for simplicity proceed only with the highest energy cluster...
178 if(clusters.size()>0)
179 {
180 Cluster c = clusters.get(0);
181 aida.cloud1D("uncorrected cluster energy for highest energy cluster").fill(c.getEnergy());
182 double e = correctClusterEnergy(c, type);
183 aida.cloud1D("corrected cluster energy for highest energy cluster").fill(e);
184 // SpacePoint p = new CartesianPoint(c.getPosition());
185 // double clusterTheta = p.theta();
186 // aida.cloud2D("uncorrected cluster energy for all clusters vs theta").fill(clusterTheta, c.getEnergy());
187 // aida.cloud2D("corrected cluster energy for all clusters vs theta").fill(clusterTheta, e);
188 // aida.cloud2D("raw-corrected cluster energy for all clusters vs theta").fill(clusterTheta, c.getEnergy()-e);
189 }
190
191 }// end of check on decays outside tracker volume
192 //reset our histogram tree
193 _tree.cd("/");
194
195 }
196
197
198
199
200 // protected void endOfData()
201 // {
202 // boolean showPlots = false;
203 // String fileType = "png";
204 // String[] pars = {"amplitude", "mean","sigma"};
205 // // int[] intEnergy = {1, 2, 5, 10, 20, 50, 100, 500 };
206 //// String fileFullPath = "C:/orglcsimAnalyses/SamplingFractionAnalysis_gamma_Theta90_acme0605.aida";
207 //// if(args.length>0) fileFullPath = args[0];
208 //// if(args.length>1)
209 //// {
210 //// fileType = args[1];
211 //// showPlots = false;
212 //// }
213 // IAnalysisFactory af = IAnalysisFactory.create();
214 //// ITree tree = af.createTreeFactory().create(fileFullPath,"xml", true, false);
215 // System.out.println("calling ls on tree");
216 // // TODO find out if I can capture output of ls command to get a list
217 // // of available directories...
218 // _tree.ls("./gamma");
219 // String[] onames = _tree.listObjectNames("./gamma");
220 //
221 // sortDirectoriesByEnergy(onames);
222 //
223 // int numberOfPoints = onames.length;
224 // double[] energies = new double[onames.length];
225 // for(int j=0; j<onames.length; ++j)
226 // {
227 //// System.out.println(onames[j]);
228 // String subDir = onames[j].substring(8); // length of "./gamma/" is 8, so remove leading directory
229 // StringTokenizer st = new StringTokenizer(subDir,"_");
230 // String e = st.nextToken();
231 // String unit = st.nextToken();
232 //// System.out.println(e+" "+unit);
233 // energies[j] = Double.parseDouble(e);
234 // if(unit.contains("MeV")) energies[j]/=1000.;
235 //// System.out.println("energy: "+energies[j]);
236 // }
237 // IFunctionFactory functionFactory = af.createFunctionFactory(_tree);
238 // IFitFactory fitFactory = af.createFitFactory();
239 // IFitter jminuit = fitFactory.createFitter("Chi2","jminuit");
240 // IFunction gauss = functionFactory.createFunctionByName("gauss","G");
241 // IFunction line = functionFactory.createFunctionByName("line","P1");
242 // IDataPointSetFactory dpsf = af.createDataPointSetFactory(_tree);
243 //
244 // IPlotter plotter = af.createPlotterFactory().create("sampling fraction plot");
245 // plotter.createRegions(3, 4, 0);
246 // IPlotterStyle style2 = plotter.region(7).style();
247 // style2.legendBoxStyle().setVisible(false);
248 // style2.statisticsBoxStyle().setVisible(false);
249 //
250 //
251 // IPlotterStyle style;
252 //
253 // double[] fitMeans = new double[numberOfPoints];
254 // double[] fitSigmas = new double[numberOfPoints];
255 // IDataPointSet energyMeans = dpsf.create("energy means vs E",2);
256 // IDataPointSet energySigmas = dpsf.create("sigma \\/ E vs E",2);
257 // IDataPointSet resolutionFit = dpsf.create("sigma \\/ E vs 1 \\/ \u221a E",2);
258 // IDataPointSet energyResiduals = dpsf.create("energy residuals (%) vs E",2);
259 // double eMax = 0;
260 // for(int i=0; i< numberOfPoints; ++i)
261 // {
262 // if(energies[i] > .1) // do not analyze 100MeV and below...
263 // {
264 // System.out.println("Energy "+energies[i]);
265 //
266 // ICloud1D e = (ICloud1D) _tree.find(onames[i]+"corrected cluster energy for highest energy cluster");
267 // e.convertToHistogram();
268 // IHistogram1D eHist = e.histogram();
269 // gauss.setParameter("amplitude",eHist.maxBinHeight());
270 // gauss.setParameter("mean",eHist.mean());
271 // gauss.setParameter("sigma",eHist.rms());
272 // style = plotter.region(i).style();
273 // style.legendBoxStyle().setVisible(false);
274 // style.statisticsBoxStyle().setVisible(false);
275 // double loElimit = energies[i] - .6*sqrt(energies[i]); // expect ~20% resolution, and go out 3 sigma
276 // double hiElimit = energies[i] + .6*sqrt(energies[i]);;
277 // plotter.region(i).setXLimits(loElimit, hiElimit);
278 // plotter.region(i).plot(eHist);
279 // IFitResult jminuitResult = jminuit.fit(eHist,gauss);
280 // double[] fitErrors = jminuitResult.errors();
281 // IFunction fit = jminuitResult.fittedFunction();
282 // for(int j=0; j<pars.length; ++j)
283 // {
284 // System.out.println(" "+pars[j]+": "+ fit.parameter(pars[j])+" +/- "+fitErrors[j]);
285 // }
286 // fitMeans[i] = fit.parameter("mean");
287 // fitSigmas[i] = fit.parameter("sigma");
288 // plotter.region(i).plot(fit);
289 //// plotter.region(7).plot(eHist);
290 //
291 // // the means
292 // IDataPoint point = energyMeans.addPoint();
293 // point.coordinate(0).setValue(energies[i]);
294 // point.coordinate(1).setValue(fitMeans[i]);
295 // point.coordinate(1).setErrorPlus(fitErrors[1]);
296 // point.coordinate(1).setErrorMinus(fitErrors[1]);
297 //
298 // // sigma
299 // IDataPoint point1 = energySigmas.addPoint();
300 // point1.coordinate(0).setValue(energies[i]);
301 // point1.coordinate(1).setValue(fitSigmas[i]/energies[i]);
302 // point1.coordinate(1).setErrorPlus(fitErrors[2]/energies[i]);
303 // point1.coordinate(1).setErrorMinus(fitErrors[2]/energies[i]);
304 //
305 // // sigma/E vs 1/sqrt(E)
306 //
307 // IDataPoint point3 = resolutionFit.addPoint();
308 // point3.coordinate(0).setValue(1./sqrt(energies[i]));
309 // point3.coordinate(1).setValue(fitSigmas[i]/energies[i]);
310 // point3.coordinate(1).setErrorPlus(fitErrors[2]/energies[i]);
311 // point3.coordinate(1).setErrorMinus(fitErrors[2]/energies[i]);
312 //
313 // // residuals
314 // IDataPoint point2 = energyResiduals.addPoint();
315 // point2.coordinate(0).setValue(energies[i]);
316 // point2.coordinate(1).setValue(100.*(fitMeans[i]-energies[i])/energies[i]);
317 //
318 // // axis bookeeping...
319 // if(energies[i] > eMax) eMax = energies[i];
320 // } // end of 100 MeV cut
321 // }
322 //
323 //
324 //
325 // IPlotter results = af.createPlotterFactory().create("linearity");
326 // style = results.region(0).style();
327 // style.xAxisStyle().setLabel("MC Energy [GeV]");
328 // style.yAxisStyle().setLabel("Cluster Energy [GeV]");
329 // style.titleStyle().setVisible(false);
330 // style.statisticsBoxStyle().setVisibileStatistics("011");
331 // style.legendBoxStyle().setVisible(false);
332 // IFitResult fitLine = jminuit.fit(energyMeans, line);
333 // System.out.println(" fit status: "+fitLine.fitStatus());
334 // double eMaxBin = eMax+10.;
335 // results.region(0).setXLimits(0., eMaxBin);
336 // results.region(0).setYLimits(0., eMaxBin);
337 // results.region(0).plot(energyMeans);
338 // results.region(0).plot(fitLine.fittedFunction());
339 //
340 //
341 // IPlotter resolution = af.createPlotterFactory().create("resolution");
342 // style = resolution.region(0).style();
343 // style.xAxisStyle().setLabel("Energy [GeV]");
344 // style.yAxisStyle().setLabel("sigma/E");
345 // style.titleStyle().setVisible(false);
346 // style.statisticsBoxStyle().setVisible(false);
347 // style.legendBoxStyle().setVisible(false);
348 // resolution.region(0).setXLimits(0., eMaxBin);
349 // resolution.region(0).setYLimits(0., .2);
350 // resolution.region(0).plot(energySigmas);
351 //
352 //
353 // IPlotter resolution2 = af.createPlotterFactory().create("sigma/E vs 1/E");
354 // style = resolution2.region(0).style();
355 // style.xAxisStyle().setLabel("1/ \u221a Energy [1/GeV]");
356 // style.yAxisStyle().setLabel("sigma/E");
357 //// style.statisticsBoxStyle().setVisibileStatistics("011");
358 // style.legendBoxStyle().setVisible(false);
359 // IFitResult resFitLine = jminuit.fit(resolutionFit, line);
360 // System.out.println(" fit status: "+resFitLine.fitStatus());
361 //// resolution2.region(0).setXLimits(0., 1.05);
362 //// resolution2.region(0).setYLimits(0., .2);
363 // resolution2.region(0).plot(resolutionFit);
364 // resolution2.region(0).plot(resFitLine.fittedFunction());
365 //
366 // IPlotter residuals = af.createPlotterFactory().create("residuals (%)");
367 // style = residuals.region(0).style();
368 // style.xAxisStyle().setLabel("Energy [GeV]");
369 // style.yAxisStyle().setLabel("Residuals [%]");
370 // style.statisticsBoxStyle().setVisible(false);
371 // style.titleStyle().setVisible(false);
372 //
373 // residuals.region(0).setXLimits(0., eMaxBin);
374 //
375 // residuals.region(0).plot(energyResiduals);
376 //
377 // if(showPlots)
378 // {
379 // plotter.show();
380 // results.show();
381 // resolution.show();
382 // resolution2.show();
383 // residuals.show();
384 // }
385 // else
386 // {
387 // try
388 // {
389 // plotter.writeToFile("energyPlots."+fileType,fileType);
390 // results.writeToFile("linearity."+fileType,fileType);
391 // resolution.writeToFile("resolution."+fileType,fileType);
392 // resolution2.writeToFile("resolutionLinear."+fileType,fileType);
393 // residuals.writeToFile("residuals."+fileType,fileType);
394 // }
395 // catch(IOException e)
396 // {
397 // System.out.println("problem writing out hardcopy in "+fileType+" format");
398 // e.printStackTrace();
399 // }
400 // }
401 // }
402
403
404
405 private double correctClusterEnergy(Cluster c, String type)
406 {
407 double e = 0.;
408 // brute force for the time being
409 // in the future we will either move this to the sampling fraction corrections
410 // or use the cluster direction for calculating the angle effect.
411 List<CalorimeterHit> hits = c.getCalorimeterHits();
412 for (CalorimeterHit hit : hits)
413 {
414 Subdetector det = hit.getSubdetector();
415 String detName = det.getName();
416 // System.out.println(detName);
417 boolean isEM = detName.startsWith("EM");
418 boolean isEndcap = det.isEndcap();
419 IDDecoder decoder = hit.getIDDecoder();
420 decoder.setID(hit.getCellID());
421 int layer = decoder.getLayer();
422 int caltype = 0;
423 if(isEM)
424 {
425 for(int i=1; i<_ecalLayering.length+1; ++i)
426 {
427 if(layer >= _ecalLayering[i-1]) caltype = i-1;
428 }
429 }
430 //TODO change this when had cal layering changes...
431 else
432 {
433 caltype = 3;
434 }
435 // System.out.println("layer= "+layer+" caltype= "+caltype);
436 String name = type + "_" +(isEM ? "em"+caltype : "had") + (isEndcap ? "e" : "b");
437 // System.out.println("fit parameter name "+name);
438 SpacePoint p = new CartesianPoint(hit.getPosition());
439 double hitTheta = p.theta();
440 double normalTheta = hitTheta;
441 // now calculate normal to the sampling plane
442 if(isEndcap)
443 {
444 normalTheta -= PI/2.;
445 normalTheta = abs(normalTheta);
446 }
447 else
448 {
449 normalTheta -= PI;
450 }
451 double a = 0.;
452 double b = 0.;
453 if(caltype==0 && !_useFirstLayer)
454 {
455 a = 0.;
456 b = 1.;
457 }
458 else
459 {
460 a = _fitParameters.get(name+"_0");
461 b = _fitParameters.get(name+"_1");
462 }
463
464 double correctionFactor = a + b*sin(normalTheta);
465 // aida.cloud2D(name+" "+layer+"hit Theta vs correction factor ").fill(hitTheta, correctionFactor);
466 // aida.cloud2D(name+" hit Theta vs correction factor ").fill(hitTheta, correctionFactor);
467 // if(isEM) aida.cloud2D("EM layer vs caltype").fill(layer, caltype);
468
469 // now apply the correction
470
471 e += hit.getRawEnergy()/correctionFactor;
472 }
473 return e;
474 }
475
476 // private static void sortDirectoriesByEnergy(String[] s)
477 // {
478 // Map<Double, String> map = new HashMap<Double, String>();
479 // double[] energies = new double[s.length];
480 // for(int j=0; j<s.length; ++j)
481 // {
482 //// System.out.println(onames[j]);
483 // String subDir = s[j].substring(8); // length of "./gamma/" is 8, so remove leading directory
484 // StringTokenizer st = new StringTokenizer(subDir,"_");
485 // String e = st.nextToken();
486 // String unit = st.nextToken();
487 //// System.out.println(e+" "+unit);
488 // energies[j] = Double.parseDouble(e);
489 // if(unit.contains("MeV")) energies[j]/=1000.;
490 // map.put(energies[j], s[j]);
491 //// System.out.println("energy: "+energies[j]);
492 // }
493 // Arrays.sort(energies);
494 // for(int j=0; j<s.length; ++j)
495 // {
496 // s[j] = map.get(energies[j]);
497 // }
498 // for(int j=0; j<s.length; ++j)
499 // {
500 // System.out.println(s[j]);
501 // }
502 //
503 //
504 //
505 // }
506
507 protected void endOfData()
508 {
509 boolean showPlots = false;
510 boolean doit = false;
511 if(doit)
512 {
513 String fileType = "png";
514 String[] pars = {"amplitude", "mean","sigma"};
515
516 IAnalysisFactory af = IAnalysisFactory.create();
517 String[] dirs = _tree.listObjectNames(".");
518 for (int ii=0; ii<dirs.length; ++ii)
519 {
520 // System.out.println("dirs["+i+"]= "+dirs[i]);
521 String[] parts = dirs[ii].split("/");
522 // for(int k=0; k<parts.length; ++k)
523 // {
524 // System.out.println("parts["+k+"]= "+parts[k]);
525 // }
526 _tree.cd(dirs[ii]);
527 String[] objects = _tree.listObjectNames(".");
528
529 // for(int j=0; j<objects.length;++j)
530 // {
531 // System.out.println("obj["+j+"]= "+objects[i]);
532 // }
533
534 sortDirectoriesByEnergy(objects);
535
536 int numberOfPoints = objects.length;
537 double[] energies = new double[objects.length];
538 for(int j=0; j<objects.length; ++j)
539 {
540 // System.out.println(objects[j]);
541
542 String subDir =parts[1];
543 String[] st = objects[j].split("/")[1].split("_");
544 String e = st[0];
545 String unit = st[1];
546 //// System.out.println(e+" "+unit);
547 energies[j] = Double.parseDouble(e);
548 if(unit.contains("MeV")) energies[j]/=1000.;
549 // System.out.println("energy: "+energies[j]);
550 }
551 IFunctionFactory functionFactory = af.createFunctionFactory(_tree);
552 IFitFactory fitFactory = af.createFitFactory();
553 IFitter jminuit = fitFactory.createFitter("Chi2","jminuit");
554 IFunction gauss = functionFactory.createFunctionByName("gauss","G");
555 IFunction line = functionFactory.createFunctionByName("line","P1");
556 IDataPointSetFactory dpsf = af.createDataPointSetFactory(_tree);
557
558 IPlotter plotter = af.createPlotterFactory().create("sampling fraction plot");
559 plotter.createRegions(3, 4, 0);
560 IPlotterStyle style2 = plotter.region(7).style();
561 style2.legendBoxStyle().setVisible(false);
562 style2.statisticsBoxStyle().setVisible(false);
563
564
565 IPlotterStyle style;
566
567 double[] fitMeans = new double[numberOfPoints];
568 double[] fitSigmas = new double[numberOfPoints];
569 IDataPointSet energyMeans = dpsf.create("energy means vs E",2);
570 IDataPointSet energySigmas = dpsf.create("sigma \\/ E vs E",2);
571 IDataPointSet resolutionFit = dpsf.create("sigma \\/ E vs 1 \\/ \u221a E",2);
572 IDataPointSet energyResiduals = dpsf.create("energy residuals (%) vs E",2);
573 double eMax = 0;
574 for(int i=0; i< numberOfPoints; ++i)
575 {
576 if(energies[i] > .1) // do not analyze 100MeV and below...
577 {
578 System.out.println("Energy "+energies[i]);
579
580 ICloud1D e = (ICloud1D) _tree.find(objects[i]+"corrected cluster energy for all clusters");
581 if(!e.isConverted()) e.convertToHistogram();
582 IHistogram1D eHist = e.histogram();
583 gauss.setParameter("amplitude",eHist.maxBinHeight());
584 gauss.setParameter("mean",eHist.mean());
585 gauss.setParameter("sigma",eHist.rms());
586 style = plotter.region(i).style();
587 style.legendBoxStyle().setVisible(false);
588 style.statisticsBoxStyle().setVisible(false);
589 double loElimit = energies[i] - .6*sqrt(energies[i]); // expect ~20% resolution, and go out 3 sigma
590 double hiElimit = energies[i] + .6*sqrt(energies[i]);;
591 plotter.region(i).setXLimits(loElimit, hiElimit);
592 plotter.region(i).plot(eHist);
593 IFitResult jminuitResult = jminuit.fit(eHist,gauss);
594 double[] fitErrors = jminuitResult.errors();
595 IFunction fit = jminuitResult.fittedFunction();
596 for(int j=0; j<pars.length; ++j)
597 {
598 System.out.println(" "+pars[j]+": "+ fit.parameter(pars[j])+" +/- "+fitErrors[j]);
599 }
600 fitMeans[i] = fit.parameter("mean");
601 fitSigmas[i] = fit.parameter("sigma");
602 plotter.region(i).plot(fit);
603 // plotter.region(7).plot(eHist);
604
605 // the means
606 IDataPoint point = energyMeans.addPoint();
607 point.coordinate(0).setValue(energies[i]);
608 point.coordinate(1).setValue(fitMeans[i]);
609 point.coordinate(1).setErrorPlus(fitErrors[1]);
610 point.coordinate(1).setErrorMinus(fitErrors[1]);
611
612 // sigma
613 IDataPoint point1 = energySigmas.addPoint();
614 point1.coordinate(0).setValue(energies[i]);
615 point1.coordinate(1).setValue(fitSigmas[i]/energies[i]);
616 point1.coordinate(1).setErrorPlus(fitErrors[2]/energies[i]);
617 point1.coordinate(1).setErrorMinus(fitErrors[2]/energies[i]);
618
619 // sigma/E vs 1/sqrt(E)
620
621 IDataPoint point3 = resolutionFit.addPoint();
622 point3.coordinate(0).setValue(1./sqrt(energies[i]));
623 point3.coordinate(1).setValue(fitSigmas[i]/energies[i]);
624 point3.coordinate(1).setErrorPlus(fitErrors[2]/energies[i]);
625 point3.coordinate(1).setErrorMinus(fitErrors[2]/energies[i]);
626
627 // residuals
628 IDataPoint point2 = energyResiduals.addPoint();
629 point2.coordinate(0).setValue(energies[i]);
630 point2.coordinate(1).setValue(100.*(fitMeans[i]-energies[i])/energies[i]);
631
632 // axis bookeeping...
633 if(energies[i] > eMax) eMax = energies[i];
634 } // end of 100 MeV cut
635 }
636
637 IPlotter results = af.createPlotterFactory().create("linearity");
638 style = results.region(0).style();
639 style.xAxisStyle().setLabel("MC Energy [GeV]");
640 style.yAxisStyle().setLabel("Cluster Energy [GeV]");
641 style.titleStyle().setVisible(false);
642 style.statisticsBoxStyle().setVisibileStatistics("011");
643 style.legendBoxStyle().setVisible(true);
644 IFitResult fitLine = jminuit.fit(energyMeans, line);
645 System.out.println(" fit status: "+fitLine.fitStatus());
646 double eMaxBin = eMax+10.;
647 results.region(0).setXLimits(0., eMaxBin);
648 results.region(0).setYLimits(0., eMaxBin);
649 results.region(0).plot(energyMeans);
650 results.region(0).plot(fitLine.fittedFunction());
651
652
653 IPlotter resolution = af.createPlotterFactory().create("resolution");
654 style = resolution.region(0).style();
655 style.xAxisStyle().setLabel("Energy [GeV]");
656 style.yAxisStyle().setLabel("sigma/E");
657 style.titleStyle().setVisible(false);
658 style.statisticsBoxStyle().setVisible(false);
659 style.legendBoxStyle().setVisible(false);
660 resolution.region(0).setXLimits(0., eMaxBin);
661 resolution.region(0).setYLimits(0., .2);
662 resolution.region(0).plot(energySigmas);
663
664
665 IPlotter resolution2 = af.createPlotterFactory().create("sigma/E vs 1/E");
666 style = resolution2.region(0).style();
667 style.xAxisStyle().setLabel("1/ \u221a Energy [1/GeV]");
668 style.yAxisStyle().setLabel("sigma/E");
669 // style.statisticsBoxStyle().setVisibileStatistics("011");
670 style.legendBoxStyle().setVisible(false);
671 IFitResult resFitLine = jminuit.fit(resolutionFit, line);
672 System.out.println(" fit status: "+resFitLine.fitStatus());
673 // resolution2.region(0).setXLimits(0., 1.05);
674 // resolution2.region(0).setYLimits(0., .2);
675 resolution2.region(0).plot(resolutionFit);
676 resolution2.region(0).plot(resFitLine.fittedFunction());
677
678 IPlotter residuals = af.createPlotterFactory().create("residuals (%)");
679 style = residuals.region(0).style();
680 style.xAxisStyle().setLabel("Energy [GeV]");
681 style.yAxisStyle().setLabel("Residuals [%]");
682 style.statisticsBoxStyle().setVisible(false);
683 style.titleStyle().setVisible(false);
684
685 residuals.region(0).setXLimits(0., eMaxBin);
686
687 residuals.region(0).plot(energyResiduals);
688
689 if(showPlots)
690 {
691 plotter.show();
692 results.show();
693 resolution.show();
694 resolution2.show();
695 residuals.show();
696 }
697 else
698 {
699 try
700 {
701 // hardcopy
702 plotter.writeToFile("energyPlots."+fileType,fileType);
703 results.writeToFile("linearity."+fileType,fileType);
704 resolution.writeToFile("resolution."+fileType,fileType);
705 resolution2.writeToFile("resolutionLinear."+fileType,fileType);
706 residuals.writeToFile("residuals."+fileType,fileType);
707 }
708 catch(IOException e)
709 {
710 System.out.println("problem writing out hardcopy in "+fileType+" format");
711 e.printStackTrace();
712 }
713
714 }
715 }// end of loop over directories
716 }
717 }
718
719 private static void sortDirectoriesByEnergy(String[] s)
720 {
721 Map<Double, String> map = new HashMap<Double, String>();
722 double[] energies = new double[s.length];
723 for(int j=0; j<s.length; ++j)
724 {
725 // System.out.println(s[j]);
726 String subDir = s[j].split("/")[1]; // first token should be "."
727 // System.out.println(subDir);
728 String[] st = subDir.split("_");
729 String e = st[0];
730 String unit = st[1];
731 // System.out.println(e+" "+unit);
732 energies[j] = Double.parseDouble(e);
733 if(unit.contains("MeV")) energies[j]/=1000.;
734 map.put(energies[j], s[j]);
735 // System.out.println("energy: "+energies[j]);
736 }
737 Arrays.sort(energies);
738 for(int j=0; j<s.length; ++j)
739 {
740 s[j] = map.get(energies[j]);
741 }
742 // for(int j=0; j<s.length; ++j)
743 // {
744 // System.out.println(s[j]);
745 // }
746 }
747
748 }