1 #include "AliHBTFits.h"
2 //_________________________________________________
3 ///////////////////////////////////////////////////////////////////////////////////
7 // Sets of methods for fittig correlation functions
12 // Piotr.Skowronski@cern.ch
14 ///////////////////////////////////////////////////////////////////////////////////
24 #include <TDirectory.h>
29 TF1* AliHBTFits::fgF1 = 0x0;
30 TF1* AliHBTFits::fgF2 = 0x0;
32 /******************************************************************************************/
33 AliHBTFits::~AliHBTFits()
39 /******************************************************************************************/
41 /******************************************************************************************/
42 /******** Qout-Qside-Qlong Cyl Surf *************************************************/
43 /******************************************************************************************/
45 void AliHBTFits::FitQOutQSideQLongCylSurf(const TString& hname,Option_t* fopt, Float_t xmax, Float_t ymax,Float_t zmax)
47 //Fits 3D histo with QOutQSideQLongCylSurf
48 TH3* h = dynamic_cast<TH3*>(gDirectory->Get(hname));
51 ::Error("FitQOutCylSurf","Can not find histogram %s",hname.Data());
55 delete gROOT->GetListOfFunctions()->FindObject("BorisTomasikCylSurfQOutQSideQLong");
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());
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();
73 TF3 *theory = new TF3("BorisTomasikCylSurfQOutQSideQLong",&AliHBTFits::QOutQSideQLongCylSurf,0.0,xmax,0.0,ymax,0.0,zmax,5);//par :aLambda, aR, aJ,aPhi,
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);
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");
88 h->Fit(theory,fopt,"E");
92 /******************************************************************************************/
93 Double_t AliHBTFits::QOutQSideQLongCylSurf(Double_t *x, Double_t *par)
95 //QOutQSideQLongCylSurf Function
97 Double_t qside = x[1];
98 Double_t qlong = x[2];
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];
106 static Double_t oldlam = aLambda;
107 static Double_t oldPhi = aPhi;
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);
116 if (oldlam != aLambda)
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);
123 Double_t result=aLambda*TMath::Exp(-(qout*qout+qside*qside+qlong*qlong)*aJ*aJ/0.0389273);
125 if ( (qlong != 0.0) && (aL != 0.0) )
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);
136 fgF1->SetParameter(0,aR);
137 fgF2->SetParameter(0,aR);
139 fgF1->SetParameter(1,aPhi);
140 fgF2->SetParameter(1,aPhi);
142 fgF1->SetParameter(2,qout);
143 fgF2->SetParameter(2,qout);
145 fgF1->SetParameter(3,qside);
146 fgF2->SetParameter(3,qside);
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);
151 // ::Info("QSideCylSurf","sp=%f, cp=%f",cp,sp);
153 result *= cp*cp + sp*sp;
154 // ::Info("QSideCylSurf","2: result=%f",result);
156 Double_t tmp = TMath::BesselI0(1./aPhi);
157 tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
160 ::Error("QOutCylSurf","Div by 0");
168 // ::Info("QSideCylSurf","Final: result=%f",result);
171 /******************************************************************************************/
172 /******** Qout Qside Cyl Surf ********************************************************/
173 /******************************************************************************************/
175 void AliHBTFits::FitQOutQSideCylSurf(const TString& hname,Option_t* fopt,Float_t xmax, Float_t ymax)
177 //Fits 3D histo with QOutQSideQLongCylSurf
178 TH2* h = dynamic_cast<TH2*>(gDirectory->Get(hname));
181 ::Error("FitQOutCylSurf","Can not find histogram %s",hname.Data());
185 delete gROOT->GetListOfFunctions()->FindObject("BorisTomasikCylSurfQOutQSide");
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());
199 if (xmax <= 0.0) xmax = h->GetXaxis()->GetXmax();
200 if (ymax <= 0.0) ymax = h->GetYaxis()->GetXmax();
202 TF2 *theory = new TF2("BorisTomasikCylSurfQOutQSide",&AliHBTFits::QOutQSideCylSurf,0.0,xmax,0.0,ymax,4);//par :aLambda, aR, aJ,aPhi,
203 theory->SetNpx(1000);
205 theory->SetParameter(0,.75);
206 theory->SetParameter(1,6.0);
207 theory->SetParameter(2,1.0);
208 theory->SetParameter(3,0.9);
210 theory->SetParName(0,"\\lambda");
211 theory->SetParName(1,"R");
212 theory->SetParName(2,"J");
213 theory->SetParName(3,"\\Phi");
215 h->Fit(theory,fopt,"E");
218 /******************************************************************************************/
220 Double_t AliHBTFits::QOutQSideCylSurf(Double_t *x, Double_t *par)
222 //QOutQSideQLongCylSurf Function
223 Double_t qout = x[0];
224 Double_t qside = x[1];
226 Double_t aLambda = par[0];
227 Double_t aR = par[1];
228 Double_t aJ = par[2];
229 Double_t aPhi = par[3];
231 static Double_t oldlam = aLambda;
233 if (oldlam != aLambda)
235 ::Info("QOutQSideQLongCylSurf","out=%f, side=%f, lambda=%f, R=%f, J=%f, Phi=%f",
236 qout,qside,aLambda,aR,aJ,aPhi);
240 Double_t result=aLambda*TMath::Exp(-(qout*qout+qside*qside)*aJ*aJ/0.0389273);
244 fgF1->SetParameter(0,aR);
245 fgF2->SetParameter(0,aR);
247 fgF1->SetParameter(1,aPhi);
248 fgF2->SetParameter(1,aPhi);
250 fgF1->SetParameter(2,qout);
251 fgF2->SetParameter(2,qout);
253 fgF1->SetParameter(3,qside);
254 fgF2->SetParameter(3,qside);
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);
259 // ::Info("QSideCylSurf","sp=%f, cp=%f",cp,sp);
261 result *= cp*cp + sp*sp;
262 // ::Info("QSideCylSurf","2: result=%f",result);
264 Double_t tmp = TMath::BesselI0(1./aPhi);
265 tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
268 ::Error("QOutCylSurf","Div by 0");
276 // ::Info("QSideCylSurf","Final: result=%f",result);
280 /******************************************************************************************/
281 /******** Qout Cyl Surf *************************************************************/
282 /******************************************************************************************/
284 void AliHBTFits::FitQOutCylSurf(const TString& hname, Option_t* fopt, Float_t max)
286 TH1D* h = dynamic_cast<TH1D*>(gDirectory->Get(hname));
289 ::Error("FitQOutCylSurf","Can not find histogram %s",hname.Data());
293 delete gROOT->GetListOfFunctions()->FindObject("BorisTomasikCylSurfQOut");
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());
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");
312 h->Fit(theory,fopt,"E");
315 /******************************************************************************************/
317 Double_t AliHBTFits::QOutCylSurf(Double_t *x, Double_t *par)
320 Double_t qout = x[0];
321 Double_t aR = par[0];
322 Double_t aPhi = par[1];
323 Double_t aLambda = par[2];
325 // ::Info("QOutCylSurf","q=%f, lambda=%f, Rout=%f, Phi=%f",qout,aLambda,aR,aPhi);
327 Double_t result=aLambda*TMath::Exp(-qout*qout/0.0389273);
329 // ::Info("QOutCylSurf","1: result=%f",result);
333 fgF1->SetParameter(0,aR);
334 fgF2->SetParameter(0,aR);
336 fgF1->SetParameter(1,aPhi);
337 fgF2->SetParameter(1,aPhi);
339 fgF1->SetParameter(2,qout);
340 fgF2->SetParameter(2,qout);
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);
345 // ::Info("QOutCylSurf","sp=%f, cp=%f",cp,sp);
347 result = result * (cp*cp + sp*sp);
348 // ::Info("QOutCylSurf","2: result=%f",result);
350 Double_t tmp = TMath::BesselI0(1./aPhi);
351 tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
354 ::Error("QOutCylSurf","Div by 0");
362 // ::Info("QOutCylSurf","Final: result=%f",result);
365 /******************************************************************************************/
366 /******** Qside Cyl Surf ************************************************************/
367 /******************************************************************************************/
369 void AliHBTFits::FitQSideCylSurf(const TString& hname, Option_t* fopt, Float_t max)
371 //Fits QSide According to Boris Tomasik formula
372 TH1D* h = dynamic_cast<TH1D*>(gDirectory->Get(hname));
375 ::Error("FitQSideCylSurf","Can not find histogram %s",hname.Data());
378 delete gROOT->GetListOfFunctions()->FindObject("BorisTomasikCylSurfQSide");
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());
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");
396 h->Fit(theory,fopt,"E");
399 /******************************************************************************************/
401 Double_t AliHBTFits::QSideCylSurf(Double_t *x, Double_t *par)
404 Double_t qside = x[0];
405 Double_t aR = par[0];
406 Double_t aPhi = par[1];
407 Double_t aLambda = par[2];
409 // ::Info("QSideCylSurf","q=%f, lambda=%f, Rside=%f, Phi=%f",qside,aLambda,aR,aPhi);
411 Double_t result=aLambda*TMath::Exp(-qside*qside/0.0389273);
414 fgF1->SetParameter(0,aR);
415 fgF2->SetParameter(0,aR);
417 fgF1->SetParameter(1,aPhi);
418 fgF2->SetParameter(1,aPhi);
420 fgF1->SetParameter(2,qside);
421 fgF2->SetParameter(2,qside);
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);
426 result *= cp*cp + sp*sp;
428 Double_t tmp = TMath::BesselI0(1./aPhi);
429 tmp = tmp*tmp*TMath::TwoPi()*TMath::TwoPi();
432 ::Error("QSideCylSurf","Div by 0");
439 /******************************************************************************************/