]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HBTAN/AliHBTFits.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
New functions added
[u/mrichter/AliRoot.git] / HBTAN / AliHBTFits.cxx
CommitLineData
59cf5e71 1#include "AliHBTFits.h"
2//_________________________________________________
3///////////////////////////////////////////////////////////////////////////////////
4//
5// class AliHBTFits
6//
7// Sets of methods for fittig correlation functions
8//
9//
10//
11//
12// Piotr.Skowronski@cern.ch
13//
14///////////////////////////////////////////////////////////////////////////////////
15
16#include <TROOT.h>
17#include <TH1.h>
d0510576 18#include <TH2.h>
59cf5e71 19#include <TH3.h>
20#include <TF1.h>
d0510576 21#include <TF2.h>
59cf5e71 22#include <TF3.h>
23#include <TError.h>
24#include <TDirectory.h>
25#include <TMath.h>
26
27ClassImp(AliHBTFits)
28
29TF1* AliHBTFits::fgF1 = 0x0;
30TF1* AliHBTFits::fgF2 = 0x0;
31
32/******************************************************************************************/
33AliHBTFits::~AliHBTFits()
34{
35 //dtor
36 delete fgF1;
37 delete fgF2;
38}
39/******************************************************************************************/
40
41/******************************************************************************************/
42/******** Qout-Qside-Qlong Cyl Surf *************************************************/
43/******************************************************************************************/
44
45void AliHBTFits::FitQOutQSideQLongCylSurf(const TString& hname,Option_t* fopt, Float_t xmax, Float_t ymax,Float_t zmax)
46{
47 //Fits 3D histo with QOutQSideQLongCylSurf
48 TH3* h = dynamic_cast<TH3*>(gDirectory->Get(hname));
49 if (h == 0x0)
50 {
51 ::Error("FitQOutCylSurf","Can not find histogram %s",hname.Data());
52 return;
53 }
54
55 delete gROOT->GetListOfFunctions()->FindObject("BorisTomasikCylSurfQOutQSideQLong");
56
57 delete fgF1;
58 delete fgF2;
59 // here x is phi
60 //par 0: R
61 //par 1: Phi
62 //par 2: qOut
63 //par 3: qSide
64 fgF1 = new TF1("ff1","sin([0]*([2]*cos(x)+[3]*sin(x))/0.197)*exp(cos(x)/[1])",-TMath::Pi(),TMath::Pi());
65 fgF2 = new TF1("ff2","cos([0]*([2]*cos(x)+[3]*sin(x))/0.197)*exp(cos(x)/[1])",-TMath::Pi(),TMath::Pi());
66 fgF1->SetNpx(1000);
67 fgF2->SetNpx(1000);
68
69 if (xmax <= 0.0) xmax = h->GetXaxis()->GetXmax();
70 if (ymax <= 0.0) ymax = h->GetYaxis()->GetXmax();
71 if (zmax <= 0.0) zmax = h->GetZaxis()->GetXmax();
72
73 TF3 *theory = new TF3("BorisTomasikCylSurfQOutQSideQLong",&AliHBTFits::QOutQSideQLongCylSurf,0.0,xmax,0.0,ymax,0.0,zmax,5);//par :aLambda, aR, aJ,aPhi,
74 theory->SetNpx(1000);
75
76 theory->SetParameter(0,.75);
77 theory->SetParameter(1,6.0);
78 theory->SetParameter(2,1.0);
79 theory->SetParameter(3,2.0);
80 theory->SetParameter(4,0.9);
81
82 theory->SetParName(0,"\\lambda");
83 theory->SetParName(1,"R");
84 theory->SetParName(2,"J");
85 theory->SetParName(3,"L");
86 theory->SetParName(4,"\\Phi");
87
88 h->Fit(theory,fopt,"E");
89
90}
91
92/******************************************************************************************/
93Double_t AliHBTFits::QOutQSideQLongCylSurf(Double_t *x, Double_t *par)
94{
95 //QOutQSideQLongCylSurf Function
96 Double_t qout = x[0];
97 Double_t qside = x[1];
98 Double_t qlong = x[2];
99
100 Double_t aLambda = par[0];
101 Double_t aR = par[1];
102 Double_t aJ = par[2];
103 Double_t aL = par[3];
104 Double_t aPhi = par[4];
105
106 static Double_t oldlam = aLambda;
107 static Double_t oldPhi = aPhi;
108
109 if (oldPhi != aPhi)
110 {
111 ::Info("QOutQSideQLongCylSurf","out=%f, side=%f, long=%f\n lambda=%f, R=%f, J=%f, L=%f Phi=%f",
112 qout,qside,qlong,aLambda,aR,aJ,aL,aPhi);
113 oldPhi = aPhi;
114 }
115
116 if (oldlam != aLambda)
117 {
118 ::Info("QOutQSideQLongCylSurf","out=%f, side=%f, long=%f\n lambda=%f, R=%f, J=%f, L=%f Phi=%f",
119 qout,qside,qlong,aLambda,aR,aJ,aL,aPhi);
120 oldlam = aLambda;
121 }
122
123 Double_t result=aLambda*TMath::Exp(-(qout*qout+qside*qside+qlong*qlong)*aJ*aJ/0.0389273);
124
125 if ( (qlong != 0.0) && (aL != 0.0) )
126 {
127 Double_t qlL = qlong*aL/2.0;
128 Double_t sinqlL = TMath::Sin(qlL);
129 Double_t sin2qlL = sinqlL*sinqlL;
130 Double_t longpart = sin2qlL/(qlL*qlL);
131 result*= longpart;
132 }
133
134 if (aPhi > 0.0)
135 {
136 fgF1->SetParameter(0,aR);
137 fgF2->SetParameter(0,aR);
138
139 fgF1->SetParameter(1,aPhi);
140 fgF2->SetParameter(1,aPhi);
141
142 fgF1->SetParameter(2,qout);
143 fgF2->SetParameter(2,qout);
144
145 fgF1->SetParameter(3,qside);
146 fgF2->SetParameter(3,qside);
147
148 Double_t sp = fgF1->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
149 Double_t cp = fgF2->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
150
151// ::Info("QSideCylSurf","sp=%f, cp=%f",cp,sp);
152
153 result *= cp*cp + sp*sp;
154// ::Info("QSideCylSurf","2: result=%f",result);
155
156 Double_t tmp = TMath::BesselI0(1./aPhi);
157 tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
158 if (tmp == 0.0)
159 {
160 ::Error("QOutCylSurf","Div by 0");
161 return 1.0;
162 }
163 result=result/tmp;
164 }
165
166 result+=1.;
167
168// ::Info("QSideCylSurf","Final: result=%f",result);
169 return result;
d0510576 170}
171/******************************************************************************************/
172/******** Qout Qside Cyl Surf ********************************************************/
173/******************************************************************************************/
174
175void AliHBTFits::FitQOutQSideCylSurf(const TString& hname,Option_t* fopt,Float_t xmax, Float_t ymax)
176{
177 //Fits 3D histo with QOutQSideQLongCylSurf
178 TH2* h = dynamic_cast<TH2*>(gDirectory->Get(hname));
179 if (h == 0x0)
180 {
181 ::Error("FitQOutCylSurf","Can not find histogram %s",hname.Data());
182 return;
183 }
184
185 delete gROOT->GetListOfFunctions()->FindObject("BorisTomasikCylSurfQOutQSide");
186
187 delete fgF1;
188 delete fgF2;
189 // here x is phi
190 //par 0: R
191 //par 1: Phi
192 //par 2: qOut
193 //par 3: qSide
194 fgF1 = new TF1("ff1","sin([0]*([2]*cos(x)+[3]*sin(x))/0.197)*exp(cos(x)/[1])",-TMath::Pi(),TMath::Pi());
195 fgF2 = new TF1("ff2","cos([0]*([2]*cos(x)+[3]*sin(x))/0.197)*exp(cos(x)/[1])",-TMath::Pi(),TMath::Pi());
196 fgF1->SetNpx(1000);
197 fgF2->SetNpx(1000);
198
199 if (xmax <= 0.0) xmax = h->GetXaxis()->GetXmax();
200 if (ymax <= 0.0) ymax = h->GetYaxis()->GetXmax();
201
202 TF2 *theory = new TF2("BorisTomasikCylSurfQOutQSide",&AliHBTFits::QOutQSideCylSurf,0.0,xmax,0.0,ymax,4);//par :aLambda, aR, aJ,aPhi,
203 theory->SetNpx(1000);
204
205 theory->SetParameter(0,.75);
206 theory->SetParameter(1,6.0);
207 theory->SetParameter(2,1.0);
208 theory->SetParameter(3,0.9);
209
210 theory->SetParName(0,"\\lambda");
211 theory->SetParName(1,"R");
212 theory->SetParName(2,"J");
213 theory->SetParName(3,"\\Phi");
214
215 h->Fit(theory,fopt,"E");
216
217}
218/******************************************************************************************/
219
220Double_t AliHBTFits::QOutQSideCylSurf(Double_t *x, Double_t *par)
221{
222 //QOutQSideQLongCylSurf Function
223 Double_t qout = x[0];
224 Double_t qside = x[1];
225
226 Double_t aLambda = par[0];
227 Double_t aR = par[1];
228 Double_t aJ = par[2];
229 Double_t aPhi = par[3];
230
231 static Double_t oldlam = aLambda;
232
233 if (oldlam != aLambda)
234 {
235 ::Info("QOutQSideQLongCylSurf","out=%f, side=%f, lambda=%f, R=%f, J=%f, Phi=%f",
236 qout,qside,aLambda,aR,aJ,aPhi);
237 oldlam = aLambda;
238 }
239
240 Double_t result=aLambda*TMath::Exp(-(qout*qout+qside*qside)*aJ*aJ/0.0389273);
241
242 if (aPhi > 0.0)
243 {
244 fgF1->SetParameter(0,aR);
245 fgF2->SetParameter(0,aR);
246
247 fgF1->SetParameter(1,aPhi);
248 fgF2->SetParameter(1,aPhi);
249
250 fgF1->SetParameter(2,qout);
251 fgF2->SetParameter(2,qout);
252
253 fgF1->SetParameter(3,qside);
254 fgF2->SetParameter(3,qside);
255
256 Double_t sp = fgF1->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
257 Double_t cp = fgF2->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
258
259// ::Info("QSideCylSurf","sp=%f, cp=%f",cp,sp);
260
261 result *= cp*cp + sp*sp;
262// ::Info("QSideCylSurf","2: result=%f",result);
263
264 Double_t tmp = TMath::BesselI0(1./aPhi);
265 tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
266 if (tmp == 0.0)
267 {
268 ::Error("QOutCylSurf","Div by 0");
269 return 1.0;
270 }
271 result=result/tmp;
272 }
273
274 result+=1.;
275
276// ::Info("QSideCylSurf","Final: result=%f",result);
277 return result;
278
59cf5e71 279}
280/******************************************************************************************/
281/******** Qout Cyl Surf *************************************************************/
282/******************************************************************************************/
283
284void AliHBTFits::FitQOutCylSurf(const TString& hname, Option_t* fopt, Float_t max)
285{
286 TH1D* h = dynamic_cast<TH1D*>(gDirectory->Get(hname));
287 if (h == 0x0)
288 {
289 ::Error("FitQOutCylSurf","Can not find histogram %s",hname.Data());
290 return;
291 }
292
293 delete gROOT->GetListOfFunctions()->FindObject("BorisTomasikCylSurfQOut");
294
295 delete fgF1;
296 delete fgF2;
297 fgF1 = new TF1("ff1","sin([0]*[2]*cos(x)/0.197)*exp(cos(x)/[1])",-TMath::Pi(),TMath::Pi());
298 fgF2 = new TF1("ff2","cos([0]*[2]*cos(x)/0.197)*exp(cos(x)/[1])",-TMath::Pi(),TMath::Pi());
299 fgF1->SetNpx(1000);
300 fgF2->SetNpx(1000);
301
302 if (max <= 0.0) max = h->GetXaxis()->GetXmax();
303 TF1 *theory = new TF1("BorisTomasikCylSurfQOut",&AliHBTFits::QOutCylSurf,0.0,max,3);//par : aR, aPhi
304 theory->SetNpx(1000);
305 theory->SetParameter(0,6.0);
306 theory->SetParameter(1,1.0);
307 theory->SetParameter(2,0.75);
308 theory->SetParName(0,"R_{out}");
309 theory->SetParName(1,"\\Phi");
310 theory->SetParName(2,"\\lambda");
311
312 h->Fit(theory,fopt,"E");
313
314}
315/******************************************************************************************/
316
317Double_t AliHBTFits::QOutCylSurf(Double_t *x, Double_t *par)
318{
319 //Qout
320 Double_t qout = x[0];
321 Double_t aR = par[0];
322 Double_t aPhi = par[1];
323 Double_t aLambda = par[2];
324
325// ::Info("QOutCylSurf","q=%f, lambda=%f, Rout=%f, Phi=%f",qout,aLambda,aR,aPhi);
326
327 Double_t result=aLambda*TMath::Exp(-qout*qout/0.0389273);
328
329// ::Info("QOutCylSurf","1: result=%f",result);
330
331 if (aPhi > 0.0)
332 {
333 fgF1->SetParameter(0,aR);
334 fgF2->SetParameter(0,aR);
335
336 fgF1->SetParameter(1,aPhi);
337 fgF2->SetParameter(1,aPhi);
338
339 fgF1->SetParameter(2,qout);
340 fgF2->SetParameter(2,qout);
341
342 Double_t sp = fgF1->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
343 Double_t cp = fgF2->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-8);
344
345// ::Info("QOutCylSurf","sp=%f, cp=%f",cp,sp);
346
347 result = result * (cp*cp + sp*sp);
348// ::Info("QOutCylSurf","2: result=%f",result);
349
350 Double_t tmp = TMath::BesselI0(1./aPhi);
351 tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
352 if (tmp == 0.0)
353 {
354 ::Error("QOutCylSurf","Div by 0");
355 return 1.0;
356 }
357 result=result/tmp;
358 }
359
360 result+=1.;
361
362// ::Info("QOutCylSurf","Final: result=%f",result);
363 return result;
364}
365/******************************************************************************************/
366/******** Qside Cyl Surf ************************************************************/
367/******************************************************************************************/
368
369void AliHBTFits::FitQSideCylSurf(const TString& hname, Option_t* fopt, Float_t max)
370{
371//Fits QSide According to Boris Tomasik formula
372 TH1D* h = dynamic_cast<TH1D*>(gDirectory->Get(hname));
373 if (h == 0x0)
374 {
375 ::Error("FitQSideCylSurf","Can not find histogram %s",hname.Data());
376 return;
377 }
378 delete gROOT->GetListOfFunctions()->FindObject("BorisTomasikCylSurfQSide");
379 delete fgF1;
380 delete fgF2;
381 fgF1 = new TF1("ff1","sin([0]*[2]*sin(x)/0.197)*exp(cos(x)/[1])",-TMath::Pi(),TMath::Pi());
382 fgF2 = new TF1("ff2","cos([0]*[2]*sin(x)/0.197)*exp(cos(x)/[1])",-TMath::Pi(),TMath::Pi());
383 fgF1->SetNpx(1000);
384 fgF2->SetNpx(1000);
385
386 if (max <= 0.0) max = h->GetXaxis()->GetXmax();
387 TF1 *theory = new TF1("BorisTomasikCylSurfQSide",&AliHBTFits::QSideCylSurf,0.0,max,3);//par : aR, aPhi
388 theory->SetNpx(1000);
389 theory->SetParameter(0,6.0);
390 theory->SetParameter(1,1.0);
391 theory->SetParameter(2,0.75);
392 theory->SetParName(0,"R_{out}");
393 theory->SetParName(1,"\\Phi");
394 theory->SetParName(2,"\\lambda");
395
396 h->Fit(theory,fopt,"E");
397
398}
399/******************************************************************************************/
400
401Double_t AliHBTFits::QSideCylSurf(Double_t *x, Double_t *par)
402{
403 //Qout
404 Double_t qside = x[0];
405 Double_t aR = par[0];
406 Double_t aPhi = par[1];
407 Double_t aLambda = par[2];
408
409// ::Info("QSideCylSurf","q=%f, lambda=%f, Rside=%f, Phi=%f",qside,aLambda,aR,aPhi);
410
411 Double_t result=aLambda*TMath::Exp(-qside*qside/0.0389273);
412 if (aPhi > 0.0)
413 {
414 fgF1->SetParameter(0,aR);
415 fgF2->SetParameter(0,aR);
416
417 fgF1->SetParameter(1,aPhi);
418 fgF2->SetParameter(1,aPhi);
419
420 fgF1->SetParameter(2,qside);
421 fgF2->SetParameter(2,qside);
422
423 Double_t sp = fgF1->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-7);
424 Double_t cp = fgF2->Integral(-TMath::Pi(),TMath::Pi(),(Double_t*)0x0,1.e-7);
425
426 result *= cp*cp + sp*sp;
427
428 Double_t tmp = TMath::BesselI0(1./aPhi);
429 tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
430 if (tmp == 0.0)
431 {
432 ::Error("QSideCylSurf","Div by 0");
433 return 1.0;
434 }
435 result/=tmp;
436 }
437 return result + 1.;
438}
439/******************************************************************************************/
440
441
442
443