PWGCF/FLOW/blastwave/AliBlastwaveFitSpectra.cxx

   1 #include<stdio.h>
   2
   3 #include"AliBlastwaveFitSpectra.h"
   4
   5 #include"TMath.h"
   6
   7 ClassImp(AliBlastwaveFitSpectra);
   8
   9
  10 TF1 *AliBlastwaveFitSpectra::fgFuncIntYield = NULL;
  11 const char *AliBlastwaveFitSpectra::fgParName[5] = {"T_{FO}","mean #rho","#gamma","mass","norm"};
  12 Float_t AliBlastwaveFitSpectra::fgStartValues[3] = {0.1,1.2,1.1};
  13 const Float_t AliBlastwaveFitSpectra::fgStepValues[3] = {0.001,0.001,0.001};
  14 const Float_t AliBlastwaveFitSpectra::fgMinValues[3] = {0.0001,0.01,0.5};
  15 const Float_t AliBlastwaveFitSpectra::fgMaxValues[3] = {1.0,10.0,5};
  16
  17 AliBlastwaveFitSpectra::AliBlastwaveFitSpectra(const char *name,Double_t mass) :
  18   AliBlastwaveFit(name,mass)
  19 {
  20   if(!fgFuncIntYield){
  21      fgFuncIntYield = new TF1("fFuncInt2Yield",AliBlastwaveFitSpectra::FunctionIntYield,0,1,5);
  22      fgFuncIntYield->SetNpx(100);
  23   }
  24
  25   Initialize();
  26 }
  27 //------------------------------------------------------------------------------
  28 AliBlastwaveFitSpectra::AliBlastwaveFitSpectra() :
  29   AliBlastwaveFit("BlastwaveFitSpectra",0.0)
  30 {
  31   if(!fgFuncIntYield){
  32      fgFuncIntYield = new TF1("fFuncInt2Yield",AliBlastwaveFitSpectra::FunctionIntYield,0,1,5);
  33      fgFuncIntYield->SetNpx(100);
  34   }
  35
  36   Initialize();
  37 }
  38 //------------------------------------------------------------------------------
  39 AliBlastwaveFitSpectra::~AliBlastwaveFitSpectra(){
  40 }
  41 //------------------------------------------------------------------------------
  42 Double_t AliBlastwaveFitSpectra::FunctionIntYield(Double_t x[],Double_t par[]){
  43   // x[0] = r/R
  44   // par[0] = pt
  45   // par[1] = T
  46   // par[2] = rho_av
  47   // par[3] = R_power
  48   // par[4] = mass
  49
  50   Double_t maxrho = par[2]*0.5*(par[3]+2);
  51
  52   Double_t mt = TMath::Sqrt(par[0]*par[0]+par[4]*par[4]);
  53   Double_t rho = maxrho*TMath::Power(x[0],par[3]);
  54   Double_t alfat = (par[0]/par[1])*TMath::SinH(rho);
  55   Double_t betat = (mt/par[1])*TMath::CosH(rho);
  56   Double_t results;
  57   if(betat < 200) results = TMath::BesselI(0,alfat)*TMath::BesselK(1,betat) *x[0];
  58   else results = 0.5*TMath::Exp(alfat-betat)/sqrt(alfat*betat) *x[0];
  59
  60   return results*par[0]*mt;
  61 }
  62 //------------------------------------------------------------------------------
  63 Int_t AliBlastwaveFitSpectra::SetParameter(Int_t ipar,Double_t val){
  64   Bool_t status = kTRUE;
  65
  66   if(ipar < 0 || ipar > 4) return status;
  67
  68   if(fFunctionYield){
  69      fFunctionYield->SetParameter(ipar,val);
  70      status = kFALSE;
  71   }
  72   return status;
  73 }
  74 //------------------------------------------------------------------------------
  75 Int_t AliBlastwaveFitSpectra::SetNormalization(){
  76   Bool_t status = kTRUE;
  77
  78   TH1 *h = (TH1 *) GetSpectrumObjCopy();
  79   if(fFunctionYield && h){
  80       fFunctionYield->SetParameter(4,1.0);
  81       Double_t intH=0;
  82       Double_t intF=0;
  83       for (Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++) {
  84           Double_t pt = h->GetBinCenter(ibin);
  85           if(pt < GetMaxPt() && pt > GetMinPt() && TMath::Abs(h->GetBinContent(ibin)) > h->GetBinError(ibin)+0.0001){
  86               Double_t width = h->GetBinWidth(ibin);
  87               intH += h->GetBinContent(ibin)*width;
  88               for(Double_t x = -width*0.45;x < width*0.5;x+=width*0.1)
  89                   intF += EvalYield(pt+x)*width*0.1;
  90           }
  91       }
  92
  93       Double_t norm = 1;
  94       if(intF > 0) norm = intH/intF;
  95       fFunctionYield->SetParameter(4,norm);
  96       status = kFALSE;
  97   }
  98
  99   return status;
 100 }
 101 //------------------------------------------------------------------------------
 102 void AliBlastwaveFitSpectra::SetMass(Double_t mass) {
 103   fMass=mass;
 104   if(fFunctionYield) fFunctionYield->SetParameter(3,fMass);
 105 }
 106 //------------------------------------------------------------------------------
 107 void AliBlastwaveFitSpectra::Initialize(){
 108   char nameFuncYield[100];
 109   snprintf(nameFuncYield,100,"%sFuncYield",GetName());
 110
 111   if(fFunctionYield) delete fFunctionYield;
 112   fFunctionYield = new TF1(nameFuncYield,Pt,0,6,7);
 113   fFunctionYield->SetParameter(3,fMass);
 114   fFunctionYield->SetNpx(200);
 115
 116   for(Int_t i =0;i < 5;i++) fFunctionYield->SetParName(i,fgParName[i]);
 117 }
 118 //------------------------------------------------------------------------------
 119 Double_t AliBlastwaveFitSpectra::Pt(Double_t x[],Double_t par[]){
 120   // par[0] = T_fo
 121   // par[1] = rho_av
 122   // par[2] = R_power
 123   // par[3] = mass
 124
 125   fgFuncIntYield->SetParameter(0,x[0]);
 126   fgFuncIntYield->SetParameter(1,par[0]);
 127   fgFuncIntYield->SetParameter(2,par[1]);
 128   fgFuncIntYield->SetParameter(3,par[2]);
 129   fgFuncIntYield->SetParameter(4,par[3]);
 130
 131
 132   Double_t res = par[4]*1E+5*fgFuncIntYield->Integral(0.,1.);
 133
 134   return res;
 135 }
 136 //------------------------------------------------------------------------------
 137 Float_t AliBlastwaveFitSpectra::GetMeanBeta(){
 138     TF1 fbeta("fbeta","TMath::TanH([0]*TMath::Power(x,[1]))*x",0,1);
 139     fbeta.SetNpx(1000);
 140
 141     Double_t maxrho = fgFuncIntYield->GetParameter(2)*0.5*(fgFuncIntYield->GetParameter(3)+2);
 142
 143     fbeta.SetParameter(0,maxrho);
 144     fbeta.SetParameter(1,fgFuncIntYield->GetParameter(3));
 145
 146     return fbeta.Integral(0.,1.)*2;
 147 }
 148 //------------------------------------------------------------------------------
 149 Float_t AliBlastwaveFitSpectra::GetMeanBeta(Double_t par[]){
 150   // par[0] = T_fo
 151   // par[1] = rho_av
 152   // par[2] = R_power
 153
 154     TF1 fbeta("fbeta","TMath::TanH([0]*TMath::Power(x,[1]))*x",0,1);
 155     fbeta.SetNpx(1000);
 156
 157     Double_t maxrho = par[1]*0.5*(par[2]+2);
 158
 159     fbeta.SetParameter(0,maxrho);
 160     fbeta.SetParameter(1,par[2]);
 161
 162     return fbeta.Integral(0.,1.)*2;
 163 }