View Javadoc

1   package org.lcsim.cal.calib;
2   /*
3    * ClusterEnergyAnalysisHistogramFitter.java
4    *
5    * Created on October 19, 2006, 12:46 PM
6    *
7    * $Id: ClusterEnergyAnalysisHistogramFitter.java,v 1.1 2008/07/30 22:10:10 ngraf Exp $
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   * @author Norman Graf
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          //      int[] intEnergy = {1, 2, 5, 10, 20, 50, 100,  500 };
41  //        String fileFullPath   = "C:/orglcsimAnalyses/SamplingFractionAnalysis_gamma_Theta90_acme0605.aida";
42  //        String fileFullPath   = "C:/lcsim/gammaTheta90_Runem_20080722.aida";
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  //            System.out.println("dirs["+i+"]= "+dirs[i]);
63              String[] parts = dirs[ii].split("/");
64  //            for(int k=0; k<parts.length; ++k)
65  //            {
66  //                System.out.println("parts["+k+"]= "+parts[k]);
67  //            }
68              tree.cd(dirs[ii]);
69              String[] objects = tree.listObjectNames(".");
70              
71  //            for(int j=0; j<objects.length;++j)
72  //            {
73  //                System.out.println("obj["+j+"]= "+objects[i]);
74  //            }
75              
76              sortDirectoriesByEnergy(objects);
77              
78              // standard style for fitted functions
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              // standard style for histogram data
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              // standard style for point data
93              IPlotterStyle pointStyle = af.createPlotterFactory().createPlotterStyle();
94              pointStyle.dataStyle().markerStyle().setColor("RED");
95              pointStyle.dataStyle().errorBarStyle().setColor("RED");
96  //            pointStyle.dataStyle().lineStyle().setVisible(true);
97              pointStyle.dataStyle().outlineStyle().setColor("BLACK");
98  //            pointStyle.dataStyle().lineStyle().setColor("RED");
99  //            pointStyle.dataStyle().fillStyle().setVisible(false);
100             
101             int numberOfPoints = objects.length;
102             double[] energies = new double[objects.length];
103             for(int j=0; j<objects.length; ++j)
104             {
105 //                System.out.println(objects[j]);
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 ////            System.out.println(e+" "+unit);
112                 energies[j] = Double.parseDouble(e);
113                 if(unit.contains("MeV")) energies[j]/=1000.;
114 //                System.out.println("energy: "+energies[j]);
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) // do not analyze 100MeV and below...
146                 {
147                     System.out.println("Energy "+energies[i]);
148                     
149                     ICloud1D e = (ICloud1D) tree.find(objects[i]+"corrected cluster energy for highest energy cluster");//"corrected cluster energy");
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]); // expect ~20% resolution, and go out 3 sigma
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 //            plotter.region(7).plot(eHist);
174                     
175                     // the means
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                     // sigma
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                     // sigma/E vs 1/sqrt(E)
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                     // residualsPercent %
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                     // residualsPercent [GeV]
203                     IDataPoint point4 = energyResiduals.addPoint();
204                     point4.coordinate(0).setValue(energies[i]);
205                     point4.coordinate(1).setValue(fitMeans[i]-energies[i]);
206                     
207                     // axis bookeeping...
208                     if(energies[i] > eMax) eMax = energies[i];
209                 } // end of 100 MeV cut
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 //            style.statisticsBoxStyle().setVisibileStatistics("011");
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 //        resolution2.region(0).setXLimits(0., 1.05);
262 //        resolution2.region(0).setYLimits(0., .2);
263             resolution2.region(0).plot(resolutionFit, pointStyle);
264             resolution2.region(0).plot(resFitLine.fittedFunction(),fitStyle);
265             
266             // residuals
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             // residuals %
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                 // hardcopy
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         }// end of loop over directories
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 //           System.out.println(s[j]);
323             String subDir = s[j].split("/")[1]; // first token should be "."
324 //            System.out.println(subDir);
325             String[] st = subDir.split("_");
326             String e = st[0];
327             String unit = st[1];
328 //            System.out.println(e+" "+unit);
329             energies[j] = Double.parseDouble(e);
330             if(unit.contains("MeV")) energies[j]/=1000.;
331             map.put(energies[j], s[j]);
332 //            System.out.println("energy: "+energies[j]);
333         }
334         Arrays.sort(energies);
335         for(int j=0; j<s.length; ++j)
336         {
337             s[j] = map.get(energies[j]);
338         }
339 //        for(int j=0; j<s.length; ++j)
340 //        {
341 //            System.out.println(s[j]);
342 //        }
343     }
344     
345 }