PWGCF/FLOW/blastwave/AliBlastwaveFit2D.cxx

   1 #include<stdio.h>
   2
   3 #include"AliBlastwaveFit2D.h"
   4
   5 #include"TMath.h"
   6 #include"TH1.h"
   7 #include"TGraphErrors.h"
   8
   9 ClassImp(AliBlastwaveFit2D);
  10
  11
  12 TF2 *AliBlastwaveFit2D::fgFuncIntYield = NULL;
  13 TF2 *AliBlastwaveFit2D::fgFuncIntV2 = NULL;
  14 const char *AliBlastwaveFit2D::fgParName[7] = {"T_{FO}","s_{2}","mean #rho_{0}","#rho_{2}","#gamma","mass","norm"};
  15 Float_t AliBlastwaveFit2D::fgStartValues[5] = {0.1,0.057,1.2,0.025,1.1};
  16 const Float_t AliBlastwaveFit2D::fgStepValues[5] = {0.001,0.001,0.001,0.001,0.001};
  17 const Float_t AliBlastwaveFit2D::fgMinValues[5] = {0.0001,-0.45,0.01,-0.45,0.1};
  18 const Float_t AliBlastwaveFit2D::fgMaxValues[5] = {1.0,0.45,10.0,0.45,5};
  19
  20 AliBlastwaveFit2D::AliBlastwaveFit2D(const char *name,Double_t mass) :
  21   AliBlastwaveFit(name,mass)
  22 {
  23   if(!fgFuncIntYield){
  24      fgFuncIntYield = new TF2("fFuncInt2Yield",AliBlastwaveFit2D::FunctionIntYield,0,2*TMath::Pi(),0,1,7);
  25      fgFuncIntYield->SetNpx(100);
  26      fgFuncIntYield->SetNpy(100);
  27   }
  28   if(!fgFuncIntV2){
  29      fgFuncIntV2 = new TF2("fFuncInt2V2",AliBlastwaveFit2D::FunctionIntV2,0,2*TMath::Pi(),0,1,7);
  30      fgFuncIntV2->SetNpx(100);
  31      fgFuncIntV2->SetNpy(100);
  32   }
  33
  34   Initialize();
  35 }
  36 //------------------------------------------------------------------------------
  37 AliBlastwaveFit2D::AliBlastwaveFit2D() :
  38   AliBlastwaveFit()
  39 {
  40   if(!fgFuncIntYield){
  41      fgFuncIntYield = new TF2("fFuncInt2Yield",AliBlastwaveFit2D::FunctionIntYield,0,2*TMath::Pi(),0,1,7);
  42      fgFuncIntYield->SetNpx(100);
  43      fgFuncIntYield->SetNpy(100);
  44   }
  45   if(!fgFuncIntV2){
  46      fgFuncIntV2 = new TF2("fFuncInt2V2",AliBlastwaveFit2D::FunctionIntV2,0,2*TMath::Pi(),0,1,7);
  47      fgFuncIntV2->SetNpx(100);
  48      fgFuncIntV2->SetNpy(100);
  49   }
  50   Initialize();
  51 }
  52 //------------------------------------------------------------------------------
  53 AliBlastwaveFit2D::~AliBlastwaveFit2D(){
  54 }
  55 //------------------------------------------------------------------------------
  56 Double_t AliBlastwaveFit2D::FunctionIntYield(Double_t x[],Double_t par[]){
  57   // x[0] = phi
  58   // x[1] = r/R
  59   // par[0] = pt
  60   // par[1] = T
  61   // par[2] = s2
  62   // par[3] = rh0_0
  63   // par[4] = rho_a
  64   // par[5] = R_power
  65   // par[6] = mass
  66
  67   Double_t mt = TMath::Sqrt(par[0]*par[0]+par[6]*par[6]);
  68   Double_t rho = (par[3]+par[4]*TMath::Cos(2*x[0]))*TMath::Power(x[1],par[5]);
  69   Double_t alfat = (par[0]/par[1])*TMath::SinH(rho);
  70   Double_t betat = (mt/par[1])*TMath::CosH(rho);
  71   Double_t results;
  72   if(betat < 200) results = TMath::BesselI(0,alfat)*TMath::BesselK(1,betat)*(1+2*par[2]*TMath::Cos(2*x[0])) *x[1];
  73   else results = 0.5*TMath::Exp(alfat-betat)/sqrt(alfat*betat)*(1+2*par[2]*TMath::Cos(2*x[0])) *x[1];
  74
  75   return results*par[0]*mt;
  76 }
  77 //------------------------------------------------------------------------------
  78 Double_t AliBlastwaveFit2D::FunctionIntV2(Double_t x[],Double_t par[]){
  79   // x[0] = phi
  80   // x[1] = r/R
  81   // par[0] = pt
  82   // par[1] = T_fo
  83   // par[2] = s2
  84   // par[3] = rh0_0
  85   // par[4] = rho_a
  86   // par[5] = R_power
  87   // par[6] = mass
  88
  89   Double_t mt = TMath::Sqrt(par[0]*par[0]+par[6]*par[6]);
  90   Double_t rho = (par[3]+par[4]*TMath::Cos(2*x[0]))*TMath::Power(x[1],par[5]);
  91   Double_t alfat = (par[0]/par[1])*TMath::SinH(rho);
  92   Double_t betat = (mt/par[1])*TMath::CosH(rho);
  93   Double_t results;
  94   if(betat < 200) results = TMath::Cos(2*x[0])*TMath::BesselI(2,alfat)*TMath::BesselK(1,betat)*(1+2*par[2]*TMath::Cos(2*x[0])) *x[1];
  95   else results = TMath::Cos(2*x[0])*0.5*TMath::Exp(alfat-betat)/sqrt(alfat*betat)*(1+2*par[2]*TMath::Cos(2*x[0])) *x[1];
  96
  97   return results*par[0]*mt;
  98 }
  99 //------------------------------------------------------------------------------
 100 Int_t AliBlastwaveFit2D::SetParameter(Int_t ipar,Double_t val){
 101   Bool_t status = kTRUE;
 102
 103   if(ipar < 0 || ipar > 4) return status;
 104
 105   if(fFunctionYield){
 106      fFunctionYield->SetParameter(ipar,val);
 107      status = kFALSE;
 108   }
 109   if(fFunctionV2){
 110      fFunctionV2->SetParameter(ipar,val);
 111      status = kFALSE;
 112   }
 113   return status;
 114 }
 115 //------------------------------------------------------------------------------
 116 Int_t AliBlastwaveFit2D::SetNormalization(){
 117   Bool_t status = kTRUE;
 118
 119   TH1 *h = (TH1 *) GetSpectrumObjCopy();
 120   if(fFunctionYield && h){
 121       fFunctionYield->SetParameter(6,1.0);
 122       Double_t intH=0;
 123       Double_t intF=0;
 124       for (Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++) {
 125           Double_t pt = h->GetBinCenter(ibin);
 126           if(pt < GetMaxPt() && pt > GetMinPt() && TMath::Abs(h->GetBinContent(ibin)) > h->GetBinError(ibin)+0.0001){
 127               Double_t width = h->GetBinWidth(ibin);
 128               intH += h->GetBinContent(ibin)*width;
 129               intF += EvalYield(pt)*width;
 130 //            for(Double_t x = -width*0.45;x < width*0.5;x+=width*0.1)
 131 //                intF += EvalYield(pt+x)*width*0.1;
 132           }
 133       }
 134
 135       Double_t norm = 1;
 136       if(intF > 0) norm = intH/intF;
 137       fFunctionYield->SetParameter(6,norm);
 138       status = kFALSE;
 139   }
 140   return status;
 141 }
 142 //------------------------------------------------------------------------------
 143 void AliBlastwaveFit2D::SetMass(Double_t mass) {
 144   fMass=mass;
 145   if(fFunctionV2) fFunctionV2->SetParameter(5,fMass);
 146   if(fFunctionYield) fFunctionYield->SetParameter(5,fMass);
 147 }
 148 //------------------------------------------------------------------------------
 149 void AliBlastwaveFit2D::Initialize(){
 150   char nameFuncYield[100];
 151   char nameFuncV2[100];
 152   snprintf(nameFuncYield,100,"%sFuncYield",GetName());
 153   snprintf(nameFuncV2,100,"%sFuncV2",GetName());
 154
 155   if(fFunctionV2) delete fFunctionV2;
 156   if(fFunctionYield) delete fFunctionYield;
 157   fFunctionV2 = new TF1(nameFuncV2,V2,0,6,6);
 158   fFunctionV2->SetParameter(5,fMass);
 159   fFunctionYield = new TF1(nameFuncYield,Pt,0,6,7);
 160   fFunctionYield->SetParameter(5,fMass);
 161   fFunctionV2->SetNpx(200);
 162   fFunctionYield->SetNpx(200);
 163
 164   for(Int_t i =0;i < 6;i++) fFunctionV2->SetParName(i,fgParName[i]);
 165   for(Int_t i =0;i < 7;i++) fFunctionYield->SetParName(i,fgParName[i]);
 166 }
 167 //------------------------------------------------------------------------------
 168 Double_t AliBlastwaveFit2D::V2(Double_t x[],Double_t par[]){
 169   // par[0] = T_fo
 170   // par[1] = s2
 171   // par[2] = rh0_mean
 172   // par[3] = rho_v2
 173   // par[4] = R_power
 174   // par[5] = mass
 175   Double_t maxRho = par[2]*0.5*(par[4]+2)/(1+2*par[1]*par[3]);
 176
 177   fgFuncIntYield->SetParameter(0,x[0]);
 178   fgFuncIntYield->SetParameter(1,par[0]);
 179   fgFuncIntYield->SetParameter(2,par[1]);
 180   fgFuncIntYield->SetParameter(3,maxRho);
 181   fgFuncIntYield->SetParameter(4,par[3]*maxRho*2);
 182   fgFuncIntYield->SetParameter(5,par[4]);
 183   fgFuncIntYield->SetParameter(6,par[5]);
 184
 185   fgFuncIntV2->SetParameter(0,x[0]);
 186   fgFuncIntV2->SetParameter(1,par[0]);
 187   fgFuncIntV2->SetParameter(2,par[1]);
 188   fgFuncIntV2->SetParameter(3,maxRho);
 189   fgFuncIntV2->SetParameter(4,par[3]*maxRho*2);
 190   fgFuncIntV2->SetParameter(5,par[4]);
 191   fgFuncIntV2->SetParameter(6,par[5]);
 192
 193   Double_t res = fgFuncIntV2->Integral(0.,2*TMath::Pi(),0.,1.);
 194   Double_t den = fgFuncIntYield->Integral(0.,2*TMath::Pi(),0.,1.);
 195
 196   if(den == 0) return 0.0;
 197
 198   return res/den;
 199 }
 200 //------------------------------------------------------------------------------
 201 Double_t AliBlastwaveFit2D::Pt(Double_t x[],Double_t par[]){
 202   // par[0] = T_fo
 203   // par[1] = s2
 204   // par[2] = rh0_mean
 205   // par[3] = rho_v2
 206   // par[4] = R_power
 207   // par[5] = mass
 208
 209   Double_t maxRho = par[2]*0.5*(par[4]+2)/(1+2*par[1]*par[3]);
 210
 211   fgFuncIntYield->SetParameter(0,x[0]);
 212   fgFuncIntYield->SetParameter(1,par[0]);
 213   fgFuncIntYield->SetParameter(2,par[1]);
 214   fgFuncIntYield->SetParameter(3,maxRho);
 215   fgFuncIntYield->SetParameter(4,par[3]*maxRho*2);
 216   fgFuncIntYield->SetParameter(5,par[4]);
 217   fgFuncIntYield->SetParameter(6,par[5]);
 218
 219   fgFuncIntV2->SetParameter(0,x[0]);
 220   fgFuncIntV2->SetParameter(1,par[0]);
 221   fgFuncIntV2->SetParameter(2,par[1]);
 222   fgFuncIntV2->SetParameter(3,maxRho);
 223   fgFuncIntV2->SetParameter(4,par[3]*maxRho*2);
 224   fgFuncIntV2->SetParameter(5,par[4]);
 225   fgFuncIntV2->SetParameter(6,par[5]);
 226
 227   Double_t res = par[6]*1E+5*fgFuncIntYield->Integral(0.,2*TMath::Pi(),0.,1.);
 228
 229   return res;
 230 }
 231 //------------------------------------------------------------------------------
 232 Float_t AliBlastwaveFit2D::GetMeanBeta(){
 233
 234   TF2 fbeta("fbeta","TMath::TanH(([0]+[1]*TMath::Cos(2*x))*TMath::Power(y,[2]))*y*(1+2*[3]*cos(2*x))",0,2*TMath::Pi(),0,1);
 235   fbeta.SetNpx(1000);
 236   fbeta.SetNpy(1000);
 237
 238   fbeta.SetParameter(0,fgFuncIntYield->GetParameter(3));
 239   fbeta.SetParameter(1,fgFuncIntYield->GetParameter(4));
 240   fbeta.SetParameter(2,fgFuncIntYield->GetParameter(5));
 241   fbeta.SetParameter(3,fgFuncIntYield->GetParameter(2));
 242
 243   return fbeta.Integral(0.,2*TMath::Pi(),0.,1.)/TMath::Pi();
 244 }
 245 //------------------------------------------------------------------------------
 246 Float_t AliBlastwaveFit2D::GetMeanBeta(Double_t par[]){
 247   // par[0] = T_fo
 248   // par[1] = s2
 249   // par[2] = rh0_mean
 250   // par[3] = rho_v2
 251   // par[4] = R_power
 252   // par[5] = mass
 253
 254   Double_t maxRho = par[2]*0.5*(par[4]+2)/(1+2*par[1]*par[3]);
 255
 256   TF2 fbeta("fbeta","TMath::TanH(([0]+[1]*TMath::Cos(2*x))*TMath::Power(y,[2]))*y*(1+2*[3]*cos(2*x))",0,2*TMath::Pi(),0,1);
 257   fbeta.SetNpx(1000);
 258   fbeta.SetNpy(1000);
 259
 260   fbeta.SetParameter(0,maxRho);
 261   fbeta.SetParameter(1,2*par[3]*maxRho);
 262   fbeta.SetParameter(2,par[4]);
 263   fbeta.SetParameter(3,par[1]);
 264
 265   return fbeta.Integral(0.,2*TMath::Pi(),0.,1.)/TMath::Pi();
 266 }
 267 //------------------------------------------------------------------------------
 268 void AliBlastwaveFit2D::SwitchOffFlow(TMinuit *m) const{
 269     m->FixParameter(kParS2);
 270     m->FixParameter(kParRho2);
 271
 272     printf("AliBlastwaveFit2D: No flow histos -> Flow parameters put to zero\n");
 273 }