StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructFitter.cxx
1 #include "StEStructFitter.h"
2 #include "TMath.h"
3 #include "TFile.h"
4 #include "Stiostream.h"
5 
6 
7 ClassImp(StEStructFitter)
8 
9 StEStructFitter* StEStructFitter::mInstance=0;
10 
11 StEStructFitter* StEStructFitter::Instance(){
12 
13  if(!mInstance)mInstance=new StEStructFitter();
14  return mInstance;
15 
16 }
17 
18 
19 StEStructFitter::StEStructFitter(): mhists(NULL) {
20  mpi=3.141593; m2pi=2.0*mpi;
21  mmean=0.;//3.3935;
22  msigma=0.;//1.38;//0.9745;
23 
24 };
25 
26 StEStructFitter::~StEStructFitter(){};
27 
28 //
29 //-----------------------------------------------------------
30 double StEStructFitter::detadphiFit(double* x, double* par){
31 
32  double arg1=0;
33  double arg2=0;
34  double arg3=0;
35  double arg4=0;
36  double yval=x[1];
37  double xval=x[0];
38 
39 
40  if(par[0]!=0)arg1=(xval)/par[0];
41  if(par[1]!=0)arg2=(yval)/par[1];
42  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
43  fitval*=TMath::Exp(-0.5*arg2*arg2);
44 
45  if(par[3]!=0)arg3=(yval-mpi)/par[3];
46  fitval+=par[4]*TMath::Exp(-0.5*arg3*arg3);
47  fitval+=par[5];
48  if(par[6]!=0)arg4=(xval)/par[6];
49  fitval+=par[7]*TMath::Exp(-0.5*arg4*arg4);
50 
51 
52  // added a small 2d-gaussian about phi_delta=pi,eta_delta=0
53  // mainly for US low pt
54  arg1=0;
55  arg2=0;
56  if(par[8]!=0)arg1=(xval)/par[8];
57  if(par[9]!=0)arg2=(yval-mpi)/par[9];
58  double fitval2=par[10]*TMath::Exp(-0.5*arg1*arg1);
59  fitval2*=TMath::Exp(-0.5*arg2*arg2);
60  fitval+=fitval2;
61 
62  if(yval<0. || yval>mpi){
63  double xnew[2];
64  xnew[0]=xval;
65  if(yval<0){
66  yval+=m2pi;
67  } else {
68  yval-=m2pi;
69  }
70  xnew[1]=yval;
71  fitval+=seconddetadphiFit(xnew,par);
72  }
73  return fitval;
74 }
75 
76 
77 double StEStructFitter::seconddetadphiFit(double* x, double* par){
78 
79  double arg1=0;
80  double arg2=0;
81  double arg3=0;
82  // double arg4=0;
83  double yval=x[1];
84  double xval=x[0];
85 
86 
87  if(par[0]!=0)arg1=(xval)/par[0];
88  if(par[1]!=0)arg2=(yval)/par[1];
89  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
90  fitval*=TMath::Exp(-0.5*arg2*arg2);
91 
92  if(par[3]!=0)arg3=(yval-mpi)/par[3];
93  fitval+=par[4]*TMath::Exp(-0.5*arg3*arg3);
94 
95  // added a small 2d-gaussian about phi_delta=pi,eta_delta=0
96  // mainly for US low pt
97  arg1=0;
98  arg2=0;
99  if(par[8]!=0)arg1=(xval)/par[8];
100  if(par[9]!=0)arg2=(yval-mpi)/par[9];
101  double fitval2=par[10]*TMath::Exp(-0.5*arg1*arg1);
102  fitval2*=TMath::Exp(-0.5*arg2*arg2);
103  fitval+=fitval2;
104 
105 
106  // fitval+=par[5];
107  // if(par[6]!=0)arg4=(xval)/par[6];
108  // fitval+=par[7]*TMath::Exp(-0.5*arg4*arg4);
109 
110 
111  return fitval;
112 }
113 
114 double StEStructFitter::softLS(double* x,double* par){
115 
116  double arg1=0;
117  double arg2=0;
118  double yval=x[1];
119  double xval=x[0];
120 
121  if(par[0]!=0)arg1=(xval)/par[0];
122  if(par[1]!=0)arg2=(yval)/par[1];
123  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
124  fitval*=TMath::Exp(-0.5*arg2*arg2);
125  fitval+=par[3];
126  if(yval<0. || yval>mpi){
127  double xnew[2];
128  xnew[0]=xval;
129  if(yval<0){
130  yval+=m2pi;
131  } else {
132  yval-=m2pi;
133  }
134  xnew[1]=yval;
135  fitval+=secondSoftLS(xnew,par);
136  }
137 
138  return fitval;
139 }
140 
141 double StEStructFitter::secondSoftLS(double* x,double* par){
142 
143  double arg1=0;
144  double arg2=0;
145  double yval=x[1];
146  double xval=x[0];
147 
148  if(par[0]!=0)arg1=(xval)/par[0];
149  if(par[1]!=0)arg2=(yval)/par[1];
150  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
151  fitval*=TMath::Exp(-0.5*arg2*arg2);
152 
153  return fitval;
154 }
155 
156 double StEStructFitter::softUS(double* x,double* par){
157 
158  double arg1=0;
159  double xval=x[0];
160 
161  if(par[0]!=0)arg1=(xval)/par[0];
162  double fitval=par[1]*TMath::Exp(-0.5*arg1*arg1);
163  fitval+=par[2];
164 
165  return fitval;
166 }
167 
168 double StEStructFitter::syt(double* x,double* par){
169 
170  double arg1=0;
171  double arg2=0;
172  double xval=x[0];
173 
174  if(par[2]!=0)arg1=(xval-par[1])/par[2];
175  double fitval=par[3]*TMath::Exp(-0.5*arg1*arg1);
176  arg2=(xval-mmean)/msigma;
177  fitval+=(par[0]*TMath::Exp(-0.5*arg2*arg2));
178 
179  return fitval;
180 }
181 
182 
183 //------------------------------------------------------------
184 double StEStructFitter::softCD(double* x,double* par){
185 
186  double arg1=0;
187  double arg2=0;
188  double arg3=0;
189  double yval=x[1];
190  double xval=x[0];
191 
192  if(par[0]!=0)arg1=(xval)/par[0];
193  if(par[1]!=0)arg2=(yval)/par[1];
194  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
195  fitval*=TMath::Exp(-0.5*arg2*arg2);
196  if(par[3]!=0)arg3=(xval)/par[3];
197  fitval+=par[4]*TMath::Exp(-0.5*arg3*arg3);
198  fitval+=par[5];
199 
200  if(yval<0. || yval>mpi){
201  double xnew[2];
202  xnew[0]=xval;
203  if(yval<0){
204  yval+=m2pi;
205  } else {
206  yval-=m2pi;
207  }
208  xnew[1]=yval;
209  fitval+=secondSoftCD(xnew,par);
210  }
211 
212  return fitval;
213 }
214 
215 
216 double StEStructFitter::secondSoftCD(double* x,double* par){
217 
218  double arg1=0;
219  double arg2=0;
220  double yval=x[1];
221  double xval=x[0];
222 
223  if(par[0]!=0)arg1=(xval)/par[0];
224  if(par[1]!=0)arg2=(yval)/par[1];
225  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
226  fitval*=TMath::Exp(-0.5*arg2*arg2);
227 
228  return fitval;
229 }
230 
231 double StEStructFitter::hardCI(double* x,double* par){
232 
233  double arg1=0;
234  double arg2=0;
235  double arg3=0;
236  double yval=x[1];
237  double xval=x[0];
238 
239  if(par[0]!=0)arg1=(xval)/par[0];
240  if(par[1]!=0)arg2=(yval)/par[1];
241  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
242  fitval*=TMath::Exp(-0.5*arg2*arg2);
243  if(par[3]!=0)arg3=(yval-mpi)/par[3];
244  fitval+=par[4]*TMath::Exp(-0.5*arg3*arg3);
245  fitval+=(par[5]-mmean*xval*xval);
246 
247  if(yval<0. || yval>mpi){
248  double xnew[2];
249  xnew[0]=xval;
250  if(yval<0){
251  yval+=m2pi;
252  } else {
253  yval-=m2pi;
254  }
255  xnew[1]=yval;
256  fitval+=secondHardCI(xnew,par);
257  }
258 
259  return fitval;
260 }
261 
262 
263 double StEStructFitter::secondHardCI(double* x,double* par){
264 
265  double arg1=0;
266  double arg2=0;
267  double arg3=0;
268  double yval=x[1];
269  double xval=x[0];
270 
271  if(par[0]!=0)arg1=(xval)/par[0];
272  if(par[1]!=0)arg2=(yval)/par[1];
273  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
274  fitval*=TMath::Exp(-0.5*arg2*arg2);
275  if(par[3]!=0)arg3=(yval-mpi)/par[3];
276  fitval+=par[4]*TMath::Exp(-0.5*arg3*arg3);
277 
278  return fitval;
279 }
280 
281 double StEStructFitter::hardCICosine(double* x,double* par){
282 
283  double arg1=0;
284  double arg2=0;
285  double arg3=0;
286  double yval=x[1];
287  double xval=x[0];
288 
289  if(par[0]!=0)arg1=(xval)/par[0];
290  if(par[1]!=0)arg2=(yval)/par[1];
291  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
292  fitval*=TMath::Exp(-0.5*arg2*arg2);
293  if(par[3]!=0)arg3=(yval-mpi)/par[3];
294  fitval+=par[4]*TMath::Exp(-0.5*arg3*arg3);
295  fitval+=(par[5]-mmean*xval*xval);
296 
297  if(yval<0. || yval>mpi){
298  double xnew[2];
299  xnew[0]=xval;
300  if(yval<0){
301  yval+=m2pi;
302  } else {
303  yval-=m2pi;
304  }
305  xnew[1]=yval;
306  fitval+=secondHardCI(xnew,par);
307  }
308 
309  return fitval;
310 }
311 
312 
313 double StEStructFitter::secondHardCICosine(double* x,double* par){
314 
315  double arg1=0;
316  double arg2=0;
317  double arg3=0;
318  double yval=x[1];
319  double xval=x[0];
320 
321  if(par[0]!=0)arg1=(xval)/par[0];
322  if(par[1]!=0)arg2=(yval)/par[1];
323  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
324  fitval*=TMath::Exp(-0.5*arg2*arg2);
325  if(par[3]!=0)arg3=(yval-mpi)/par[3];
326  fitval+=par[4]*TMath::Exp(-0.5*arg3*arg3);
327 
328  return fitval;
329 }
330 
331 
332 
333 double StEStructFitter::detadphiSS(double* x, double* par){
334 
335  double arg1=0;
336  double arg2=0;
337  double yval=x[1];
338  double xval=x[0];
339 
340 
341  if(par[0]!=0)arg1=(xval)/par[0];
342  if(par[1]!=0)arg2=(yval)/par[1];
343  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
344  fitval*=TMath::Exp(-0.5*arg2*arg2);
345 
346  fitval+=par[3];
347 
348  /*
349  if(yval<0. || yval>mpi){
350  double xnew[2];
351  xnew[0]=xval;
352  if(yval<0){
353  yval+=m2pi;
354  } else {
355  yval-=m2pi;
356  }
357  xnew[1]=yval;
358  fitval+=seconddetadphiSS(xnew,par);
359  }
360  */
361 
362  return fitval;
363 }
364 
365 
366 double StEStructFitter::seconddetadphiSS(double* x, double* par){
367 
368  double arg1=0;
369  double arg2=0;
370  double yval=x[1];
371  double xval=x[0];
372 
373  if(par[0]!=0)arg1=(xval)/par[0];
374  if(par[1]!=0)arg2=(yval)/par[1];
375  double fitval=par[2]*TMath::Exp(-0.5*arg1*arg1);
376  fitval*=TMath::Exp(-0.5*arg2*arg2);
377 
378  return fitval;
379 }
380 
381 
382 //------------------------------------------------------
383 
384 double StEStructFitter::mcComponents(double* x, double* par){
385 
386  if(!mhists){
387  cout<<"opening file rmxTest2.root"<<endl;
388  TFile* tf=new TFile("rmxTest2.root");
389  mhists=new TH2F*[6];
390 
391  const char* snames[]={"S001","S010","S100","S011","S101","S110"};
392 
393  double xsum=0;
394  for(int i=0;i<6;i++){
395  mhists[i]=(TH2F*)(tf->Get(snames[i]))->Clone();
396  float xval=mhists[i]->Integral();
397  cout<<"scale for i="<<i<<" is "<<xval<<endl;
398  xsum+=mhists[i]->Integral();
399  }
400  for(int i=0;i<6;i++)mhists[i]->Scale(1.0/xsum);
401  }
402 
403  int ibin=mhists[0]->FindBin(x[0],x[1]);
404 
405  double retVal=par[0]*mhists[0]->GetBinContent(ibin);
406  retVal+=par[1]*mhists[1]->GetBinContent(ibin);
407  retVal+=par[2]*mhists[2]->GetBinContent(ibin);
408  retVal+=par[3]*mhists[3]->GetBinContent(ibin);
409  retVal+=par[4]*mhists[4]->GetBinContent(ibin);
410  retVal+=par[5]*mhists[5]->GetBinContent(ibin);
411  retVal*=par[6];
412 
413  return retVal;
414 
415 }
416 
417 //-------------------------------------------------------
418 double StEStructFitter::dytGsytG(double* x, double* par){
419 
420 
421  double arg1=0;
422  double arg2=0;
423 
424  double yval=x[1];
425  double xval=x[0];
426 
427  if(par[0]!=0)arg1=(xval-par[1])/par[0];
428  if(par[1]!=0)arg2=(yval)/par[2];
429 
430  double fitval=par[3]*TMath::Exp(-0.5*arg1*arg1);
431  fitval*=TMath::Exp(-0.5*arg2*arg2);
432  fitval+=par[4];
433 
434  return fitval;
435 }
436 
437 
438 //---------------------------------------------------------
439 double StEStructFitter::DoubleE(double* x, double* par){
440 
441  double arg1=0;
442  double arg2=0;
443  double yval=x[1];
444  double xval=x[0]; // assume x independent....
445 
446  if(par[0]!=0)arg1=(yval)/par[0];
447  if(par[2]!=0)arg2=(yval-mpi)/par[2];
448  double fitval=par[1]*TMath::Exp(-0.5*arg1*arg1);
449  fitval+=par[3]*TMath::Exp(-0.5*arg2*arg2);
450  fitval+=par[4];
451 
452  if(yval<0. || yval>mpi){
453  double xnew[2];
454  xnew[0]=xval;
455  if(yval<0){
456  yval+=m2pi;
457  } else {
458  yval-=m2pi;
459  }
460  xnew[1]=yval;
461  fitval+=secondDoubleE(xnew,par);
462  }
463 
464  return fitval;
465 }
466 
467 
468 double StEStructFitter::secondDoubleE(double* x,double* par){
469 
470  double arg1=0;
471  double arg2=0;
472  double yval=x[1];
473  //double xval=x[0];
474 
475  if(par[0]!=0)arg1=(yval)/par[0];
476  if(par[2]!=0)arg2=(yval-mpi)/par[2];
477  double fitval=par[1]*TMath::Exp(-0.5*arg1*arg1);
478  fitval+=par[3]*TMath::Exp(-0.5*arg2*arg2);
479 
480  return fitval;
481 }