HBTAN/AliHBTFits.cxx

   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>
  18 #include <TH2.h>
  19 #include <TH3.h>
  20 #include <TF1.h>
  21 #include <TF2.h>
  22 #include <TF3.h>
  23 #include <TError.h>
  24 #include <TDirectory.h>
  25 #include <TMath.h>
  26
  27 ClassImp(AliHBTFits)
  28
  29 TF1*     AliHBTFits::fgF1 = 0x0;
  30 TF1*     AliHBTFits::fgF2 = 0x0;
  31
  32 /******************************************************************************************/
  33 AliHBTFits::~AliHBTFits()
  34 {
  35  //dtor
  36  delete fgF1;
  37  delete fgF2;
  38 }
  39 /******************************************************************************************/
  40
  41 /******************************************************************************************/
  42 /********    Qout-Qside-Qlong  Cyl Surf   *************************************************/
  43 /******************************************************************************************/
  44
  45 void 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 /******************************************************************************************/
  93 Double_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;
 170 }
 171 /******************************************************************************************/
 172 /********    Qout Qside Cyl Surf   ********************************************************/
 173 /******************************************************************************************/
 174
 175 void 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
 220 Double_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
 279 }
 280 /******************************************************************************************/
 281 /********    Qout  Cyl Surf   *************************************************************/
 282 /******************************************************************************************/
 283
 284 void 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
 317 Double_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
 369 void 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
 401 Double_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