StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructAutoFit.cxx
1 #include "StEStructAutoFit.h"
2 #include "TMath.h"
3 #include "TH2D.h"
4 #include "TVector.h"
5 #include "TMatrixD.h"
6 #include "TDecompSVD.h"
7 #include "TF1.h"
8 #include "TF2.h"
9 #include "TRandom2.h"
10 #include "Stiostream.h"
11 
12 
13 using namespace TMath;
14 
15 
16 const int numIterations=500;
17 
18 // Electron peak, dipole
19 Double_t fun1(Double_t *x, Double_t *par) {
20  //Double_t v1 = par[0] * Cos(x[1]);
21  Double_t v1 = par[0] * (1- Cos(x[1]))/2.;
22  Double_t v2 = par[1] * Cos(2.*x[1]);
23  Double_t geta = par[2] * Exp(-pow((x[0]/par[3]), 2)/2.);
24  Double_t d1 = Exp(-pow(x[1]/par[4], 2)/2.);
25  Double_t d2 = Exp(-pow((x[1]-2.*Pi())/par[4], 2)/2.);
26  Double_t g1 = geta*(d1+d2);
27  Double_t g2 = par[5] * Exp(-pow((x[0]/par[6]), 2)/2.);
28  Double_t eg = par[8] * Exp(-Sqrt(pow(x[0]/par[9], 2)+pow(x[1]/par[10], 2))/2.);
29 
30  return v1+v2+g1+g2+par[7]+eg;
31 }
32 
33 // No electron peak, dipole
34 Double_t fun2(Double_t *x, Double_t *par) {
35  //Double_t v1 = par[0] * Cos(x[1]);
36  Double_t v1 = par[0] * (1- Cos(x[1]))/2.;
37  Double_t v2 = par[1] * Cos(2.*x[1]);
38  Double_t geta = par[2] * Exp(-pow((x[0]/par[3]), 2)/2.);
39  Double_t d1 = Exp(-pow(x[1]/par[4], 2)/2.);
40  Double_t d2 = Exp(-pow((x[1]-2.*Pi())/par[4], 2)/2.);
41  Double_t g1 = geta*(d1+d2);
42  Double_t g2 = par[5] * Exp(-pow((x[0]/par[6]), 2)/2.);
43 
44  return v1+v2+g1+g2+par[7];
45 }
46 // Same as above but with v3 term
47 Double_t fun2v3(Double_t *x, Double_t *par) {
48  //Double_t v1 = par[0] * Cos(x[1]);
49  Double_t v1 = par[0] * (1- Cos(x[1]))/2.;
50  Double_t v2 = par[1] * Cos(2.*x[1]);
51  Double_t v3 = par[8] * Cos(3.*x[1]);
52  Double_t geta = par[2] * Exp(-pow((x[0]/par[3]), 2)/2.);
53  Double_t d1 = Exp(-pow(x[1]/par[4], 2)/2.);
54  Double_t d2 = Exp(-pow((x[1]-2.*Pi())/par[4], 2)/2.);
55  Double_t g1 = geta*(d1+d2);
56  Double_t g2 = par[5] * Exp(-pow((x[0]/par[6]), 2)/2.);
57 
58  return v1+v2+g1+g2+par[7]+v3;
59 }
60 Double_t fun3(Double_t *x, Double_t *par) {
61  //Double_t v1 = par[0] * Cos(x[1]);
62  Double_t v1 = par[0] * (1- Cos(x[1]))/2.;
63  Double_t v2 = par[1] * Cos(2.*x[1]);
64  Double_t geta = par[2] * Exp(-pow((x[0]/par[3]), 2)/2.);
65  Double_t d1 = Exp(-pow(x[1]/par[4], 2)/2.);
66  Double_t d2 = Exp(-pow((x[1]-2.*Pi())/par[4], 2)/2.);
67  Double_t g1 = geta*(d1+d2);
68  Double_t g2a = par[5] * Exp(-pow((x[0]/par[6]), 2)/2.);
69  Double_t g2b = par[7] * Exp(-pow((x[0]/par[8]), 2)/2.);
70 
71  return v1+v2+g1+g2a+g2b+par[9];
72 }
73 Double_t fun4(Double_t *x, Double_t *par) {
74  //Double_t v1 = par[0] * Cos(x[1]);
75  Double_t v1 = par[0] * (1- Cos(x[1]))/2.;
76  Double_t v2 = par[1] * Cos(2.*x[1]);
77  Double_t geta = par[2] * Exp(-pow((x[0]/par[3]), 2)/2.);
78  Double_t d1 = Exp(-pow(x[1]/par[4], 2)/2.);
79  Double_t d2 = Exp(-pow((x[1]-2.*Pi())/par[4], 2)/2.);
80  Double_t g1 = geta*(d1+d2);
81  Double_t g2 = par[5] * Exp(-pow((x[0]/par[6]), 10)/10.);
82 
83  return v1+v2+g1+g2+par[7];
84 }
85 
86 // No electron peak, AS Gauss
87 Double_t fun5(Double_t *x, Double_t *par) {
88  Double_t v1 = par[0] * (Exp(-pow((x[1]-Pi())/par[8], 2)/2.)+Exp(-pow((x[1]+Pi())/par[8], 2)/2.));
89  Double_t v2 = par[1] * Cos(2.*x[1]);
90  Double_t geta = par[2] * Exp(-pow((x[0]/par[3]), 2)/2.);
91  Double_t d1 = Exp(-pow(x[1]/par[4], 2)/2.);
92  Double_t d2 = Exp(-pow((x[1]-2.*Pi())/par[4], 2)/2.);
93  Double_t g1 = geta*(d1+d2);
94  Double_t g2 = par[5] * Exp(-pow((x[0]/par[6]), 2)/2.);
95 
96  return v1+v2+g1+g2+par[7];
97 }
98 
99 // Electron peak, AS Gauss
100 Double_t fun6(Double_t *x, Double_t *par) {
101  Double_t v1 = par[0] * (Exp(-pow((x[1]-Pi())/par[8], 2)/2.)+Exp(-pow((x[1]+Pi())/par[8], 2)/2.));
102  Double_t v2 = par[1] * Cos(2.*x[1]);
103  Double_t geta = par[2] * Exp(-pow((x[0]/par[3]), 2)/2.);
104  Double_t d1 = Exp(-pow(x[1]/par[4], 2)/2.);
105  Double_t d2 = Exp(-pow((x[1]-2.*Pi())/par[4], 2)/2.);
106  Double_t g1 = geta*(d1+d2);
107  Double_t g2 = par[5] * Exp(-pow((x[0]/par[6]), 2)/2.);
108  Double_t eg = par[9] * Exp(-Sqrt(pow(x[0]/par[10], 2)+pow(x[1]/par[11], 2))/2.);
109 
110  return v1+v2+g1+g2+par[7]+eg;
111 }
112 // No electron peak, dipole, two SS Gaussians
113 Double_t fun7(Double_t *x, Double_t *par) {
114  //Double_t v1 = par[0] * Cos(x[1]);
115  Double_t v1 = par[0] * (1- Cos(x[1]))/2.;
116  Double_t v2 = par[1] * Cos(2.*x[1]);
117  Double_t geta = par[2] * Exp(-pow((x[0]/par[3]), 2)/2.);
118  Double_t d1 = Exp(-pow(x[1]/par[4], 2)/2.);
119  Double_t d2 = Exp(-pow((x[1]-2.*Pi())/par[4], 2)/2.);
120  Double_t g1 = geta*(d1+d2);
121  Double_t g2 = par[5] * Exp(-pow((x[0]/par[6]), 2)/2.);
122  Double_t getab = par[8] * Exp(-pow((x[0]/par[9]), 2)/2.);
123  Double_t d1b = Exp(-pow(x[1]/par[10], 2)/2.);
124  Double_t d2b = Exp(-pow((x[1]-2.*Pi())/par[10], 2)/2.);
125  Double_t g1b = getab*(d1b+d2b);
126 
127  return v1+v2+g1+g2+par[7]+g1b;
128 }
129 
130 // Electron peak, dipole
131 Double_t fun8(Double_t *x, Double_t *par) {
132  //Double_t v1 = par[0] * Cos(x[1]);
133  Double_t v1 = par[0] * (1- Cos(x[1]))/2.;
134  Double_t v2 = par[1] * Cos(2.*x[1]);
135  Double_t v3 = par[11] * Cos(3.*x[1]);
136  Double_t geta = par[2] * Exp(-pow((x[0]/par[3]), 2)/2.);
137  Double_t d1 = Exp(-pow(x[1]/par[4], 2)/2.);
138  Double_t d2 = Exp(-pow((x[1]-2.*Pi())/par[4], 2)/2.);
139  Double_t g1 = geta*(d1+d2);
140  Double_t g2 = par[5] * Exp(-pow((x[0]/par[6]), 2)/2.);
141  Double_t eg = par[8] * Exp(-Sqrt(pow(x[0]/par[9], 2)+pow(x[1]/par[10], 2))/2.);
142 
143  return v1+v2+g1+g2+par[7]+eg+v3;
144 }
145 
146 
147 // New algorithm just based on random sampling
148 // Less sophisticated, but perhaps more reliable
149 
150 double* StEStructAutoFit::autofit8Par(double* best, TH2D* plot, int type, double* allchisq) {
151  // First set the error in the (0,0) bin high so it isn't part of the fit
152  plot->SetBinError(13, 7, 1000.);
153  plot->SetBinError(12, 7, 1000.);
154  plot->SetBinError(11, 7, 1000.);
155  plot->SetBinError(14, 7, 1000.);
156  plot->SetBinError(15, 7, 1000.);
157  plot->SetBinError(13, 6, 1000.);
158  plot->SetBinError(13, 8, 1000.);
159 
160  double amprange = plot->GetMaximum() - plot->GetMinimum();
161 
162  // The range of values will just be based on RMS and means of the histogramsa
163  double mean=0.;
164  for(int i=1; i<=25; i++) {
165  for(int k=1; k<=24; k++) {
166  mean += plot->GetBinContent(i, k);
167  }
168  }
169  mean *= 1./(25.*24.);
170 
171  double rms2=0.;
172  for(int i=1; i<=25; i++) {
173  for(int k=1; k<=24; k++) {
174  rms2 += pow(plot->GetBinContent(i, k)-mean, 2);
175  //rms += (plot->GetBinContent(i, k)-mean) * (plot->GetBinContent(i, k)-mean);
176  //rms2 += plot->GetBinContent(i, k);
177  }
178  }
179  rms2 *= 1./(25.*24.);
180  rms2 = sqrt(rms2);
181 
182  cout << "\nmean=" << mean << "\trms=" << rms2 << endl;
183 
184  double minphi=.35;
185  double maxphi=.9;
186 
187  double mineta=.4;
188  double maxeta=3.0;
189 
190  double mincos=0.;
191  //double maxcos=2.0*rms2;
192  double maxcos=amprange;
193 
194  double minv2=0.;
195  //double maxv2=2.0*rms2;
196  double maxv2=amprange;
197 
198  double ming2=0.;
199  //double maxg2=5.*rms2;
200  double maxg2=amprange;
201 
202  double ming1=-2.*rms2;
203  if(type == 1)
204  ming1=0.;
205  double maxg1=2.*rms2;
206 
207  double mingw=.1;
208  double maxgw=5;
209 
210  double minoffset=mean-3.*rms2;
211  double maxoffset=mean+rms2;
212 
213  TRandom2 rand2;
214  rand2.SetSeed(0);
215 
216  double chisqsum=0;
217  double chisqsqsum=0;
218  double bestchisq=999999;
219  double worstchisq=0;
220  int numConverge=0;
221  for(int i=0; i<numIterations; i++) {
222  double phi = rand2.Rndm()*(maxphi-minphi)+minphi;
223  double eta = rand2.Rndm()*(maxeta-mineta)+mineta;
224  double cosv = rand2.Rndm()*(maxcos-mincos)+mincos;
225  double v2 = rand2.Rndm()*(maxv2-minv2)+minv2;
226  double g2 = rand2.Rndm()*(maxg2-ming2)+ming2;
227  double g1 = rand2.Rndm()*(maxg1-ming1)+ming1;
228  double gw = rand2.Rndm()*(maxgw-mingw)+mingw;
229  double offset = rand2.Rndm()*(maxoffset-minoffset)+minoffset;
230 
231  cout << "\nphi=" << phi << "\teta=" << eta << "\tcos=" << cosv << "\tv2=" << v2 << "\tg2=" << g2 << "\tg1=" << g1 << "\tgw=" << gw << "\toffset=" << offset << endl;
232 
233 
234  const Int_t npar = 8;
235 
236  Double_t f2params[npar] = {cosv, v2, g2, eta, phi, g1, gw, offset};
237 
238  TF2 *f2 = new TF2("f2",fun2,-2,2,-Pi()/2.,3*Pi()/2, npar);
239  f2->SetParameters(f2params);
240  //f2->SetParLimits(0, -5, 0);
241  //f2->FixParameter(1, 0.011704);
242  //f2->FixParameter(1, 0.00364);
243  if(type == 2)
244  f2->SetParLimits(1, 0, 5);
245  f2->SetParLimits(2, 0, 5);
246  f2->SetParLimits(3, 0.1, 5);
247  f2->SetParLimits(4, 0.1, 2);
248  if(type == 1 || type == 2)
249  f2->SetParLimits(5, 0, 3);
250  f2->SetParLimits(6, 0.12, 7);
251 
252  int fitStatus=plot->Fit("f2","N");
253  fitStatus=plot->Fit("f2","EN");
254 
255  if(fitStatus==0) {
256  double chisq=f2->GetChisquare();
257  cout << "\n" << f2->GetParameter(2) << endl;
258  chisqsum+=chisq;
259  chisqsqsum+=chisq*chisq;
260  if(chisq<bestchisq) {
261  for(int i=0; i<npar; i++) {
262  best[i*2]=f2->GetParameter(i);
263  best[i*2+1]=f2->GetParError(i);
264  }
265  bestchisq=chisq;
266  }
267  if(chisq>worstchisq) {
268  worstchisq=chisq;
269  }
270  numConverge++;
271  if(allchisq!=NULL) {allchisq[i]=chisq;}
272  } else {
273  if(allchisq!=NULL) {allchisq[i]=0.;}
274  }
275  delete f2;
276  }
277 
278  cout << "\nOut of " << numIterations << " attempts, " << numConverge << " actually converged." << endl;
279  cout << "\nThe best chisq is: " << bestchisq << endl;
280  cout << "\nThe worst chisq is: " << worstchisq << endl;
281 
282 // if(storeChi) {
283 // chiPlot = new TH1D("chiPlot","Chi Squared",100,bestchisq,worstchisq);
284 // //chiPlot->SetBins(100,bestchisq,worstchisq);
285 // for(int i=0; i<numIterations; i++) {
286 // cerr << "\tIteration: " << i;
287 // if(allchisq[i]!=0.) {
288 // chiPlot->Fill(allchisq[i]);
289 // }
290 // }
291 // }
292 // cerr << "\nIn function plot is: " << chiPlot;
293 // cerr << "\nIn function RMS is: " << chiPlot->GetRMS();
294 
295 
296  return best;
297 }
298 
299 // Now with v3
300 
301 double* StEStructAutoFit::autofit8Parv3(double* best, TH2D* plot, int type, double* allchisq) {
302  // First set the error in the (0,0) bin high so it isn't part of the fit
303  plot->SetBinError(13, 7, 1000.);
304  plot->SetBinError(12, 7, 1000.);
305  plot->SetBinError(11, 7, 1000.);
306  plot->SetBinError(14, 7, 1000.);
307  plot->SetBinError(15, 7, 1000.);
308  plot->SetBinError(13, 6, 1000.);
309  plot->SetBinError(13, 8, 1000.);
310 
311  double amprange = plot->GetMaximum() - plot->GetMinimum();
312 
313  // The range of values will just be based on RMS and means of the histogramsa
314  double mean=0.;
315  for(int i=1; i<=25; i++) {
316  for(int k=1; k<=24; k++) {
317  mean += plot->GetBinContent(i, k);
318  }
319  }
320  mean *= 1./(25.*24.);
321 
322  double rms2=0.;
323  for(int i=1; i<=25; i++) {
324  for(int k=1; k<=24; k++) {
325  rms2 += pow(plot->GetBinContent(i, k)-mean, 2);
326  //rms += (plot->GetBinContent(i, k)-mean) * (plot->GetBinContent(i, k)-mean);
327  //rms2 += plot->GetBinContent(i, k);
328  }
329  }
330  rms2 *= 1./(25.*24.);
331  rms2 = sqrt(rms2);
332 
333  cout << "\nmean=" << mean << "\trms=" << rms2 << endl;
334 
335  double minphi=.35;
336  double maxphi=.9;
337 
338  double mineta=.4;
339  double maxeta=3.0;
340 
341  double mincos=0.;
342  //double maxcos=2.0*rms2;
343  double maxcos=amprange;
344 
345  double minv2=0.;
346  //double maxv2=2.0*rms2;
347  double maxv2=amprange;
348 
349  double minv3=0.;
350  //double maxv3=2.0*rms2;
351  double maxv3=amprange;
352 
353  double ming2=0.;
354  //double maxg2=5.*rms2;
355  double maxg2=amprange;
356 
357  double ming1=-2.*rms2;
358  if(type == 1)
359  ming1=0.;
360  double maxg1=2.*rms2;
361 
362  double mingw=.1;
363  double maxgw=5;
364 
365  double minoffset=mean-3.*rms2;
366  double maxoffset=mean+rms2;
367 
368  TRandom2 rand2;
369  rand2.SetSeed(0);
370 
371  double chisqsum=0;
372  double chisqsqsum=0;
373  double bestchisq=999999;
374  double worstchisq=0;
375  int numConverge=0;
376  for(int i=0; i<numIterations; i++) {
377  double phi = rand2.Rndm()*(maxphi-minphi)+minphi;
378  double eta = rand2.Rndm()*(maxeta-mineta)+mineta;
379  double cosv = rand2.Rndm()*(maxcos-mincos)+mincos;
380  double v2 = rand2.Rndm()*(maxv2-minv2)+minv2;
381  double v3 = rand2.Rndm()*(maxv3-minv3)+minv3;
382  double g2 = rand2.Rndm()*(maxg2-ming2)+ming2;
383  double g1 = rand2.Rndm()*(maxg1-ming1)+ming1;
384  double gw = rand2.Rndm()*(maxgw-mingw)+mingw;
385  double offset = rand2.Rndm()*(maxoffset-minoffset)+minoffset;
386 
387  cout << "\nphi=" << phi << "\teta=" << eta << "\tcos=" << cosv << "\tv2=" << v2 << "\tg2=" << g2 << "\tg1=" << g1 << "\tgw=" << gw << "\toffset=" << offset << endl;
388 
389 
390  const Int_t npar = 9;
391 
392  Double_t f2params[npar] = {cosv, v2, g2, eta, phi, g1, gw, offset, v3};
393 
394  TF2 *f2 = new TF2("f2",fun2v3,-2,2,-Pi()/2.,3*Pi()/2, npar);
395  f2->SetParameters(f2params);
396  //f2->SetParLimits(0, -5, 0);
397  //f2->FixParameter(1, 0.011704);
398  //f2->FixParameter(1, 0.00364);
399  if(type == 2)
400  f2->SetParLimits(1, 0, 5);
401  f2->SetParLimits(2, 0, 5);
402  f2->SetParLimits(3, 0.1, 5);
403  f2->SetParLimits(4, 0.1, 2);
404  if(type == 1 || type == 2)
405  f2->SetParLimits(5, 0, 3);
406  f2->SetParLimits(6, 0.01, 7);
407 
408  int fitStatus=plot->Fit("f2","N");
409  fitStatus=plot->Fit("f2","EN");
410 
411  if(fitStatus==0) {
412  double chisq=f2->GetChisquare();
413  cout << "\n" << f2->GetParameter(2) << endl;
414  chisqsum+=chisq;
415  chisqsqsum+=chisq*chisq;
416  if(chisq<bestchisq) {
417  for(int i=0; i<npar; i++) {
418  best[i*2]=f2->GetParameter(i);
419  best[i*2+1]=f2->GetParError(i);
420  }
421  bestchisq=chisq;
422  }
423  if(chisq>worstchisq) {
424  worstchisq=chisq;
425  }
426  numConverge++;
427  if(allchisq!=NULL) {allchisq[i]=chisq;}
428  } else {
429  if(allchisq!=NULL) {allchisq[i]=0.;}
430  }
431  delete f2;
432  }
433 
434  cout << "\nOut of " << numIterations << " attempts, " << numConverge << " actually converged." << endl;
435  cout << "\nThe best chisq is: " << bestchisq << endl;
436  cout << "\nThe worst chisq is: " << worstchisq << endl;
437 
438 // if(storeChi) {
439 // chiPlot = new TH1D("chiPlot","Chi Squared",100,bestchisq,worstchisq);
440 // //chiPlot->SetBins(100,bestchisq,worstchisq);
441 // for(int i=0; i<numIterations; i++) {
442 // cerr << "\tIteration: " << i;
443 // if(allchisq[i]!=0.) {
444 // chiPlot->Fill(allchisq[i]);
445 // }
446 // }
447 // }
448 // cerr << "\nIn function plot is: " << chiPlot;
449 // cerr << "\nIn function RMS is: " << chiPlot->GetRMS();
450 
451 
452  return best;
453 }
454 
455 // Exponential electron peak
456 double* StEStructAutoFit::autofit11Par(double* best, TH2D* plot, int type, double* allchisq) {
457 
458  double amprange = plot->GetMaximum() - plot->GetMinimum();
459 
460  // The range of values will just be based on RMS and means of the histogramsa
461  double mean=0.;
462  for(int i=1; i<=25; i++) {
463  for(int k=1; k<=24; k++) {
464  mean += plot->GetBinContent(i, k);
465  }
466  }
467  mean *= 1./(25.*24.);
468 
469  double rms2=0.;
470  for(int i=1; i<=25; i++) {
471  for(int k=1; k<=24; k++) {
472  rms2 += pow(plot->GetBinContent(i, k)-mean, 2);
473  //rms += (plot->GetBinContent(i, k)-mean) * (plot->GetBinContent(i, k)-mean);
474  //rms2 += plot->GetBinContent(i, k);
475  }
476  }
477  rms2 *= 1./(25.*24.);
478  rms2 = sqrt(rms2);
479 
480  cout << "\nmean=" << mean << "\trms=" << rms2 << endl;
481 
482  double minphi=.35;
483  double maxphi=.9;
484 
485  double mineta=.4;
486  double maxeta=2.8;
487 
488  double mincos=0.;
489  //double maxcos=2.0*rms2;
490  double maxcos=amprange;
491 
492  double minv2=0.;
493  //double maxv2=2.0*rms2;
494  double maxv2=amprange;
495 
496  double ming2=0.;
497  //double maxg2=5.*rms2;
498  double maxg2=amprange;
499 
500  double ming1=-2.*rms2;
501  if(type == 1)
502  ming1=0.;
503  double maxg1=2.*rms2;
504 
505  double mingw=.1;
506  double maxgw=5;
507 
508  double minoffset=mean-3.*rms2;
509  double maxoffset=mean+rms2;
510 
511  double minexp=0.;
512  double maxexp=10.*rms2;
513 
514  double minew=.01;
515  double maxew=.8;
516 
517  TRandom2 rand2;
518  rand2.SetSeed(0);
519 
520  double chisqsum=0;
521  double chisqsqsum=0;
522  double bestchisq=999999;
523  double worstchisq=0;
524  int numConverge=0;
525  for(int i=0; i<numIterations; i++) {
526  double phi = rand2.Rndm()*(maxphi-minphi)+minphi;
527  double eta = rand2.Rndm()*(maxeta-mineta)+mineta;
528  double cosv = rand2.Rndm()*(maxcos-mincos)+mincos;
529  double v2 = rand2.Rndm()*(maxv2-minv2)+minv2;
530  double g2 = rand2.Rndm()*(maxg2-ming2)+ming2;
531  double g1 = rand2.Rndm()*(maxg1-ming1)+ming1;
532  double gw = rand2.Rndm()*(maxgw-mingw)+mingw;
533  double offset = rand2.Rndm()*(maxoffset-minoffset)+minoffset;
534  double exp = rand2.Rndm()*(maxexp-minexp)+minexp;
535  double expwe = rand2.Rndm()*(maxew-minew)+minew;
536  double expwp = rand2.Rndm()*(maxew-minew)+minew;
537 
538  cout << "\nphi=" << phi << "\teta=" << eta << "\tcos=" << cosv << "\tv2=" << v2 << "\tg2=" << g2 << "\tg1=" << g1 << "\tgw=" << gw << "\toffset=" << offset << endl;
539 
540 
541  const Int_t npar = 11;
542 
543  Double_t f2params[npar] = {cosv, v2, g2, eta, phi, g1, gw, offset, exp, expwe, expwp};
544 
545  TF2 *f2 = new TF2("f2",fun1,-2,2,-Pi()/2.,3*Pi()/2, npar);
546  f2->SetParameters(f2params);
547  //f2->SetParLimits(0, -5, 0);
548  //f2->FixParameter(1, 0.011704);
549  //f2->FixParameter(1, 0.00364);
550  if(type == 2)
551  f2->SetParLimits(1, 0, 5);
552  f2->SetParLimits(2, 0, 5);
553  f2->SetParLimits(3, 0.1, 5);
554  f2->SetParLimits(4, 0.1, 1.2);
555  if(type == 1 || type == 2)
556  f2->SetParLimits(5, 0, 3);
557  f2->SetParLimits(6, 0.01, 7);
558  f2->SetParLimits(8, 0., 5);
559  f2->SetParLimits(9, 0.005, .9);
560  f2->SetParLimits(10, 0.005, .9);
561 
562  int fitStatus=plot->Fit("f2","N");
563  fitStatus=plot->Fit("f2","EN");
564 
565  if(fitStatus==0) {
566  double chisq=f2->GetChisquare();
567  cout << "\n" << f2->GetParameter(2) << endl;
568  chisqsum+=chisq;
569  chisqsqsum+=chisq*chisq;
570  if(chisq<bestchisq) {
571  for(int i=0; i<npar; i++) {
572  best[i*2]=f2->GetParameter(i);
573  best[i*2+1]=f2->GetParError(i);
574  }
575  bestchisq=chisq;
576  }
577  if(chisq>worstchisq) {
578  worstchisq=chisq;
579  }
580  numConverge++;
581  if(allchisq!=NULL) {allchisq[i]=chisq;}
582  } else {
583  if(allchisq!=NULL) {allchisq[i]=0.;}
584  }
585  delete f2;
586  }
587 
588  cout << "\nOut of " << numIterations << " attempts, " << numConverge << " actually converged." << endl;
589  cout << "\nThe best chisq is: " << bestchisq << endl;
590  cout << "\nThe worst chisq is: " << worstchisq << endl;
591 
592  return best;
593 }
594 
595 
596 // Exponential electron peak + v3
597 double* StEStructAutoFit::autofit11Parv3(double* best, TH2D* plot, int type, double* allchisq) {
598 
599  double amprange = plot->GetMaximum() - plot->GetMinimum();
600 
601  // The range of values will just be based on RMS and means of the histogramsa
602  double mean=0.;
603  for(int i=1; i<=25; i++) {
604  for(int k=1; k<=24; k++) {
605  mean += plot->GetBinContent(i, k);
606  }
607  }
608  mean *= 1./(25.*24.);
609 
610  double rms2=0.;
611  for(int i=1; i<=25; i++) {
612  for(int k=1; k<=24; k++) {
613  rms2 += pow(plot->GetBinContent(i, k)-mean, 2);
614  //rms += (plot->GetBinContent(i, k)-mean) * (plot->GetBinContent(i, k)-mean);
615  //rms2 += plot->GetBinContent(i, k);
616  }
617  }
618  rms2 *= 1./(25.*24.);
619  rms2 = sqrt(rms2);
620 
621  cout << "\nmean=" << mean << "\trms=" << rms2 << endl;
622 
623  double minphi=.35;
624  double maxphi=.9;
625 
626  double mineta=.4;
627  double maxeta=2.8;
628 
629  double mincos=0.;
630  //double maxcos=2.0*rms2;
631  double maxcos=amprange;
632 
633  double minv2=0.;
634  //double maxv2=2.0*rms2;
635  double maxv2=amprange;
636 
637  double minv3=0.;
638  //double maxv3=2.0*rms2;
639  double maxv3=amprange;
640 
641  double ming2=0.;
642  //double maxg2=5.*rms2;
643  double maxg2=amprange;
644 
645  double ming1=-2.*rms2;
646  if(type == 1)
647  ming1=0.;
648  double maxg1=2.*rms2;
649 
650  double mingw=.1;
651  double maxgw=5;
652 
653  double minoffset=mean-3.*rms2;
654  double maxoffset=mean+rms2;
655 
656  double minexp=0.;
657  double maxexp=10.*rms2;
658 
659  double minew=.01;
660  double maxew=1.;
661 
662  TRandom2 rand2;
663  rand2.SetSeed(0);
664 
665  double chisqsum=0;
666  double chisqsqsum=0;
667  double bestchisq=999999;
668  double worstchisq=0;
669  int numConverge=0;
670  for(int i=0; i<numIterations; i++) {
671  double phi = rand2.Rndm()*(maxphi-minphi)+minphi;
672  double eta = rand2.Rndm()*(maxeta-mineta)+mineta;
673  double cosv = rand2.Rndm()*(maxcos-mincos)+mincos;
674  double v2 = rand2.Rndm()*(maxv2-minv2)+minv2;
675  double v3 = rand2.Rndm()*(maxv3-minv3)+minv3;
676  double g2 = rand2.Rndm()*(maxg2-ming2)+ming2;
677  double g1 = rand2.Rndm()*(maxg1-ming1)+ming1;
678  double gw = rand2.Rndm()*(maxgw-mingw)+mingw;
679  double offset = rand2.Rndm()*(maxoffset-minoffset)+minoffset;
680  double exp = rand2.Rndm()*(maxexp-minexp)+minexp;
681  double expwe = rand2.Rndm()*(maxew-minew)+minew;
682  double expwp = rand2.Rndm()*(maxew-minew)+minew;
683 
684  cout << "\nphi=" << phi << "\teta=" << eta << "\tcos=" << cosv << "\tv2=" << v2 << "\tg2=" << g2 << "\tg1=" << g1 << "\tgw=" << gw << "\toffset=" << offset << endl;
685 
686 
687  const Int_t npar = 12;
688 
689  Double_t f2params[npar] = {cosv, v2, g2, eta, phi, g1, gw, offset, exp, expwe, expwp, v3};
690 
691  TF2 *f2 = new TF2("f2",fun1,-2,2,-Pi()/2.,3*Pi()/2, npar);
692  f2->SetParameters(f2params);
693  //f2->SetParLimits(0, -5, 0);
694  //f2->FixParameter(1, 0.011704);
695  //f2->FixParameter(1, 0.00364);
696  if(type == 2)
697  f2->SetParLimits(1, 0, 5);
698  f2->SetParLimits(2, 0, 5);
699  f2->SetParLimits(3, 0.1, 5);
700  f2->SetParLimits(4, 0.1, 2);
701  if(type == 1 || type == 2)
702  f2->SetParLimits(5, 0, 3);
703  f2->SetParLimits(6, 0.01, 7);
704  f2->SetParLimits(8, 0., 5);
705  f2->SetParLimits(9, 0.005, 2);
706  f2->SetParLimits(10, 0.005, 2);
707 
708  int fitStatus=plot->Fit("f2","N");
709  fitStatus=plot->Fit("f2","EN");
710 
711  if(fitStatus==0) {
712  double chisq=f2->GetChisquare();
713  cout << "\n" << f2->GetParameter(2) << endl;
714  chisqsum+=chisq;
715  chisqsqsum+=chisq*chisq;
716  if(chisq<bestchisq) {
717  for(int i=0; i<npar; i++) {
718  best[i*2]=f2->GetParameter(i);
719  best[i*2+1]=f2->GetParError(i);
720  }
721  bestchisq=chisq;
722  }
723  if(chisq>worstchisq) {
724  worstchisq=chisq;
725  }
726  numConverge++;
727  if(allchisq!=NULL) {allchisq[i]=chisq;}
728  } else {
729  if(allchisq!=NULL) {allchisq[i]=0.;}
730  }
731  delete f2;
732  }
733 
734  cout << "\nOut of " << numIterations << " attempts, " << numConverge << " actually converged." << endl;
735  cout << "\nThe best chisq is: " << bestchisq << endl;
736  cout << "\nThe worst chisq is: " << worstchisq << endl;
737 
738  return best;
739 }
740 
741 // New algorithm just based on random sampling
742 // Less sophisticated, but perhaps more reliable
743 // Away-side Gaussian verion
744 
745 double* StEStructAutoFit::autofit9Par(double* best, TH2D* plot, int type, double* allchisq) {
746  // First set the error in the (0,0) bin high so it isn't part of the fit
747  plot->SetBinError(13, 7, 1000.);
748  plot->SetBinError(12, 7, 1000.);
749  plot->SetBinError(11, 7, 1000.);
750  plot->SetBinError(14, 7, 1000.);
751  plot->SetBinError(15, 7, 1000.);
752  plot->SetBinError(13, 6, 1000.);
753  plot->SetBinError(13, 8, 1000.);
754 
755  double amprange = plot->GetMaximum() - plot->GetMinimum();
756 
757  // The range of values will just be based on RMS and means of the histogramsa
758  double mean=0.;
759  for(int i=1; i<=25; i++) {
760  for(int k=1; k<=24; k++) {
761  mean += plot->GetBinContent(i, k);
762  }
763  }
764  mean *= 1./(25.*24.);
765 
766  double rms2=0.;
767  for(int i=1; i<=25; i++) {
768  for(int k=1; k<=24; k++) {
769  rms2 += pow(plot->GetBinContent(i, k)-mean, 2);
770  //rms += (plot->GetBinContent(i, k)-mean) * (plot->GetBinContent(i, k)-mean);
771  //rms2 += plot->GetBinContent(i, k);
772  }
773  }
774  rms2 *= 1./(25.*24.);
775  rms2 = sqrt(rms2);
776 
777  cout << "\nmean=" << mean << "\trms=" << rms2 << endl;
778 
779  double minphi=.35;
780  double maxphi=.9;
781 
782  double mineta=.4;
783  double maxeta=2.8;
784 
785  double mincos=0.;
786  //double maxcos=2.0*rms2;
787  double maxcos=amprange;
788 
789  double minv2=0.;
790  //double maxv2=2.0*rms2;
791  double maxv2=amprange;
792 
793  double ming2=0.;
794  //double maxg2=5.*rms2;
795  double maxg2=amprange;
796 
797  double ming1=-2.*rms2;
798  if(type == 1)
799  ming1=0.;
800  double maxg1=2.*rms2;
801 
802  double mingw=.1;
803  double maxgw=5;
804 
805  double minoffset=mean-3.*rms2;
806  double maxoffset=mean+rms2;
807 
808  double minagw=.5;
809  double maxagw=2.5;
810 
811  TRandom2 rand2;
812  rand2.SetSeed(0);
813 
814  double chisqsum=0;
815  double chisqsqsum=0;
816  double bestchisq=999999;
817  double worstchisq=0;
818  int numConverge=0;
819  for(int i=0; i<numIterations; i++) {
820  double phi = rand2.Rndm()*(maxphi-minphi)+minphi;
821  double eta = rand2.Rndm()*(maxeta-mineta)+mineta;
822  double cosv = rand2.Rndm()*(maxcos-mincos)+mincos;
823  double v2 = rand2.Rndm()*(maxv2-minv2)+minv2;
824  double g2 = rand2.Rndm()*(maxg2-ming2)+ming2;
825  double g1 = rand2.Rndm()*(maxg1-ming1)+ming1;
826  double gw = rand2.Rndm()*(maxgw-mingw)+mingw;
827  double offset = rand2.Rndm()*(maxoffset-minoffset)+minoffset;
828  double agw = rand2.Rndm()*(maxagw-minagw)+minagw;
829 
830  cout << "\nphi=" << phi << "\teta=" << eta << "\tcos=" << cosv << "\tv2=" << v2 << "\tg2=" << g2 << "\tg1=" << g1 << "\tgw=" << gw << "\toffset=" << offset << endl;
831 
832 
833  const Int_t npar = 9;
834 
835  Double_t f2params[npar] = {cosv, v2, g2, eta, phi, g1, gw, offset, agw};
836 
837  TF2 *f2 = new TF2("f2",fun5,-2,2,-Pi()/2.,3*Pi()/2, npar);
838  f2->SetParameters(f2params);
839  //f2->SetParLimits(0, -5, 0);
840  //f2->FixParameter(1, 0.011704);
841  //f2->FixParameter(1, 0.00364);
842  if(type == 2)
843  f2->SetParLimits(1, 0, 5);
844  f2->SetParLimits(2, 0, 5);
845  f2->SetParLimits(3, 0.1, 5);
846  f2->SetParLimits(4, 0.1, 2);
847  if(type == 1 || type == 2)
848  f2->SetParLimits(5, 0, 3);
849  f2->SetParLimits(6, 0.01, 7);
850  f2->SetParLimits(8, 0.01, 5);
851 
852  int fitStatus=plot->Fit("f2","N");
853  fitStatus=plot->Fit("f2","EN");
854 
855  if(fitStatus==0) {
856  double chisq=f2->GetChisquare();
857  cout << "\n" << f2->GetParameter(2) << endl;
858  chisqsum+=chisq;
859  chisqsqsum+=chisq*chisq;
860  if(chisq<bestchisq) {
861  for(int i=0; i<npar; i++) {
862  best[i*2]=f2->GetParameter(i);
863  best[i*2+1]=f2->GetParError(i);
864  }
865  bestchisq=chisq;
866  }
867  if(chisq>worstchisq) {
868  worstchisq=chisq;
869  }
870  numConverge++;
871  if(allchisq!=NULL) {allchisq[i]=chisq;}
872  } else {
873  if(allchisq!=NULL) {allchisq[i]=0.;}
874  }
875  delete f2;
876  }
877 
878  cout << "\nOut of " << numIterations << " attempts, " << numConverge << " actually converged." << endl;
879  cout << "\nThe best chisq is: " << bestchisq << endl;
880  cout << "\nThe worst chisq is: " << worstchisq << endl;
881 
882  return best;
883 }
884 
885 // Away-side Gaussian verion
886 // Exponential peak
887 
888 double* StEStructAutoFit::autofit12Par(double* best, TH2D* plot, int type, double* allchisq) {
889 
890  double amprange = plot->GetMaximum() - plot->GetMinimum();
891 
892  // The range of values will just be based on RMS and means of the histogramsa
893  double mean=0.;
894  for(int i=1; i<=25; i++) {
895  for(int k=1; k<=24; k++) {
896  mean += plot->GetBinContent(i, k);
897  }
898  }
899  mean *= 1./(25.*24.);
900 
901  double rms2=0.;
902  for(int i=1; i<=25; i++) {
903  for(int k=1; k<=24; k++) {
904  rms2 += pow(plot->GetBinContent(i, k)-mean, 2);
905  //rms += (plot->GetBinContent(i, k)-mean) * (plot->GetBinContent(i, k)-mean);
906  //rms2 += plot->GetBinContent(i, k);
907  }
908  }
909  rms2 *= 1./(25.*24.);
910  rms2 = sqrt(rms2);
911 
912  cout << "\nmean=" << mean << "\trms=" << rms2 << endl;
913 
914  double minphi=.35;
915  double maxphi=.9;
916 
917  double mineta=.4;
918  double maxeta=2.8;
919 
920  double mincos=0.;
921  //double maxcos=2.0*rms2;
922  double maxcos=amprange;
923 
924  double minv2=0.;
925  //double maxv2=2.0*rms2;
926  double maxv2=amprange;
927 
928  double ming2=0.;
929  //double maxg2=5.*rms2;
930  double maxg2=amprange;
931 
932  double ming1=-2.*rms2;
933  if(type == 1)
934  ming1=0.;
935  double maxg1=2.*rms2;
936 
937  double mingw=.1;
938  double maxgw=5;
939 
940  double minoffset=mean-3.*rms2;
941  double maxoffset=mean+rms2;
942 
943  double minagw=.5;
944  double maxagw=2.5;
945 
946  double minexp=0.;
947  double maxexp=10.*rms2;
948 
949  double minew=.01;
950  double maxew=1.;
951 
952  TRandom2 rand2;
953  rand2.SetSeed(0);
954 
955  double chisqsum=0;
956  double chisqsqsum=0;
957  double bestchisq=999999;
958  double worstchisq=0;
959  int numConverge=0;
960  for(int i=0; i<numIterations; i++) {
961  double phi = rand2.Rndm()*(maxphi-minphi)+minphi;
962  double eta = rand2.Rndm()*(maxeta-mineta)+mineta;
963  double cosv = rand2.Rndm()*(maxcos-mincos)+mincos;
964  double v2 = rand2.Rndm()*(maxv2-minv2)+minv2;
965  double g2 = rand2.Rndm()*(maxg2-ming2)+ming2;
966  double g1 = rand2.Rndm()*(maxg1-ming1)+ming1;
967  double gw = rand2.Rndm()*(maxgw-mingw)+mingw;
968  double offset = rand2.Rndm()*(maxoffset-minoffset)+minoffset;
969  double agw = rand2.Rndm()*(maxagw-minagw)+minagw;
970  double exp = rand2.Rndm()*(maxexp-minexp)+minexp;
971  double expwe = rand2.Rndm()*(maxew-minew)+minew;
972  double expwp = rand2.Rndm()*(maxew-minew)+minew;
973 
974  cout << "\nphi=" << phi << "\teta=" << eta << "\tcos=" << cosv << "\tv2=" << v2 << "\tg2=" << g2 << "\tg1=" << g1 << "\tgw=" << gw << "\toffset=" << offset << endl;
975 
976 
977  const Int_t npar = 12;
978 
979  Double_t f2params[npar] = {cosv, v2, g2, eta, phi, g1, gw, offset, agw, exp, expwe, expwp};
980 
981  TF2 *f2 = new TF2("f2",fun6,-2,2,-Pi()/2.,3*Pi()/2, npar);
982  f2->SetParameters(f2params);
983  //f2->SetParLimits(0, -5, 0);
984  //f2->FixParameter(1, 0.011704);
985  //f2->FixParameter(1, 0.00364);
986  if(type == 2)
987  f2->SetParLimits(1, 0, 5);
988  f2->SetParLimits(2, 0, 5);
989  f2->SetParLimits(3, 0.1, 5);
990  f2->SetParLimits(4, 0.1, 2);
991  if(type == 1 || type == 2)
992  f2->SetParLimits(5, 0, 3);
993  f2->SetParLimits(6, 0.01, 7);
994  f2->SetParLimits(8, 0.01, 5);
995  f2->SetParLimits(9, 0., 5);
996  f2->SetParLimits(10, 0.005, 2);
997  f2->SetParLimits(11, 0.005, 2);
998 
999  int fitStatus=plot->Fit("f2","N");
1000  fitStatus=plot->Fit("f2","EN");
1001 
1002  if(fitStatus==0) {
1003  double chisq=f2->GetChisquare();
1004  cout << "\n" << f2->GetParameter(2) << endl;
1005  chisqsum+=chisq;
1006  chisqsqsum+=chisq*chisq;
1007  if(chisq<bestchisq) {
1008  for(int i=0; i<npar; i++) {
1009  best[i*2]=f2->GetParameter(i);
1010  best[i*2+1]=f2->GetParError(i);
1011  }
1012  bestchisq=chisq;
1013  }
1014  if(chisq>worstchisq) {
1015  worstchisq=chisq;
1016  }
1017  numConverge++;
1018  if(allchisq!=NULL) {allchisq[i]=chisq;}
1019  } else {
1020  if(allchisq!=NULL) {allchisq[i]=0.;}
1021  }
1022  delete f2;
1023  }
1024 
1025  cout << "\nOut of " << numIterations << " attempts, " << numConverge << " actually converged." << endl;
1026  cout << "\nThe best chisq is: " << bestchisq << endl;
1027  cout << "\nThe worst chisq is: " << worstchisq << endl;
1028 
1029  return best;
1030 }
1031 
1032 
1033 double* StEStructAutoFit::autofit11Par2G(double* best, TH2D* plot, int type, double yt, double* allchisq) {
1034  // First set the error in the (0,0) bin high so it isn't part of the fit
1035  plot->SetBinError(13, 7, 1000.);
1036  plot->SetBinError(12, 7, 1000.);
1037  plot->SetBinError(11, 7, 1000.);
1038  plot->SetBinError(14, 7, 1000.);
1039  plot->SetBinError(15, 7, 1000.);
1040  plot->SetBinError(13, 6, 1000.);
1041  plot->SetBinError(13, 8, 1000.);
1042 
1043  double amprange = plot->GetMaximum() - plot->GetMinimum();
1044 
1045  // The range of values will just be based on RMS and means of the histogramsa
1046  double mean=0.;
1047  for(int i=1; i<=25; i++) {
1048  for(int k=1; k<=24; k++) {
1049  mean += plot->GetBinContent(i, k);
1050  }
1051  }
1052  mean *= 1./(25.*24.);
1053 
1054  double rms2=0.;
1055  for(int i=1; i<=25; i++) {
1056  for(int k=1; k<=24; k++) {
1057  rms2 += pow(plot->GetBinContent(i, k)-mean, 2);
1058  //rms += (plot->GetBinContent(i, k)-mean) * (plot->GetBinContent(i, k)-mean);
1059  //rms2 += plot->GetBinContent(i, k);
1060  }
1061  }
1062  rms2 *= 1./(25.*24.);
1063  rms2 = sqrt(rms2);
1064 
1065  cout << "\nmean=" << mean << "\trms=" << rms2 << endl;
1066 
1067  double minphi=.35;
1068  double maxphi=.9;
1069 
1070  double mineta=.4;
1071  double maxeta=2.8;
1072 
1073  double mincos=0.;
1074  //double maxcos=2.0*rms2;
1075  double maxcos=amprange;
1076 
1077  double minv2=0.;
1078  //double maxv2=2.0*rms2;
1079  double maxv2=amprange;
1080 
1081  double ming2=0.;
1082  //double maxg2=5.*rms2;
1083  double maxg2=amprange;
1084 
1085  double ming2b=0.;
1086  //double maxg2b=5.*rms2;
1087  double maxg2b=amprange;
1088 
1089  double ming1=-2.*rms2;
1090  if(type == 1)
1091  ming1=0.;
1092  double maxg1=2.*rms2;
1093 
1094  double mingw=.1;
1095  double maxgw=5;
1096 
1097  double minoffset=mean-3.*rms2;
1098  double maxoffset=mean+rms2;
1099 
1100  double eta2 = .9/sqrt(yt);
1101  double phi2 = 1.8/sqrt(3.*yt-2.);
1102 
1103  TRandom2 rand2;
1104  rand2.SetSeed(0);
1105 
1106  double chisqsum=0;
1107  double chisqsqsum=0;
1108  double bestchisq=999999;
1109  double worstchisq=0;
1110  int numConverge=0;
1111  for(int i=0; i<numIterations; i++) {
1112  double phi = rand2.Rndm()*(maxphi-minphi)+minphi;
1113  double eta = rand2.Rndm()*(maxeta-mineta)+mineta;
1114  double cosv = rand2.Rndm()*(maxcos-mincos)+mincos;
1115  double v2 = rand2.Rndm()*(maxv2-minv2)+minv2;
1116  double g2 = rand2.Rndm()*(maxg2-ming2)+ming2;
1117  double g1 = rand2.Rndm()*(maxg1-ming1)+ming1;
1118  double gw = rand2.Rndm()*(maxgw-mingw)+mingw;
1119  double offset = rand2.Rndm()*(maxoffset-minoffset)+minoffset;
1120  double g2b = rand2.Rndm()*(maxg2-ming2)+ming2;
1121 
1122  cout << "\nphi=" << phi << "\teta=" << eta << "\tcos=" << cosv << "\tv2=" << v2 << "\tg2=" << g2 << "\tg1=" << g1 << "\tgw=" << gw << "\toffset=" << offset << endl;
1123 
1124 
1125  const Int_t npar = 11;
1126 
1127  Double_t f2params[npar] = {cosv, v2, g2, eta, phi, g1, gw, offset, g2b, eta2, phi2};
1128 
1129  TF2 *f2 = new TF2("f2",fun7,-2,2,-Pi()/2.,3*Pi()/2, npar);
1130  f2->SetParameters(f2params);
1131  //f2->SetParLimits(0, -5, 0);
1132  //f2->FixParameter(1, 0.011704);
1133  //f2->FixParameter(1, 0.00364);
1134  if(type == 2) {
1135  f2->SetParLimits(1, 0, 5);
1136  f2->SetParLimits(2, 0, 5);
1137  }
1138  f2->SetParLimits(3, 0.1, 5);
1139  f2->SetParLimits(4, 0.1, 2);
1140  if(type == 1 || type == 2)
1141  f2->SetParLimits(5, 0, 3);
1142  f2->SetParLimits(6, 0.01, 7);
1143  if(type == 1 || type == 2)
1144  f2->SetParLimits(8, 0, 3);
1145  f2->FixParameter(9, eta2);
1146  f2->FixParameter(10, phi2);
1147 
1148  int fitStatus=plot->Fit("f2","N");
1149  fitStatus=plot->Fit("f2","EN");
1150 
1151  if(fitStatus==0) {
1152  double chisq=f2->GetChisquare();
1153  cout << "\n" << f2->GetParameter(2) << endl;
1154  chisqsum+=chisq;
1155  chisqsqsum+=chisq*chisq;
1156  if(chisq<bestchisq) {
1157  for(int i=0; i<npar; i++) {
1158  best[i*2]=f2->GetParameter(i);
1159  best[i*2+1]=f2->GetParError(i);
1160  }
1161  bestchisq=chisq;
1162  }
1163  if(chisq>worstchisq) {
1164  worstchisq=chisq;
1165  }
1166  numConverge++;
1167  if(allchisq!=NULL) {allchisq[i]=chisq;}
1168  } else {
1169  if(allchisq!=NULL) {allchisq[i]=0.;}
1170  }
1171  delete f2;
1172  }
1173 
1174  cout << "\nOut of " << numIterations << " attempts, " << numConverge << " actually converged." << endl;
1175  cout << "\nThe best chisq is: " << bestchisq << endl;
1176  cout << "\nThe worst chisq is: " << worstchisq << endl;
1177 
1178  return best;
1179 }
Generated by   doxygen 1.8.5