[u/mrichter/AliRoot.git] / PWGCF / FLOW / blastwave / AliBlastwaveFit2D.cxx

#include<stdio.h>

#include"AliBlastwaveFit2D.h"

#include"TMath.h"
#include"TH1.h"
#include"TGraphErrors.h"

ClassImp(AliBlastwaveFit2D);


TF2 *AliBlastwaveFit2D::fgFuncIntYield = NULL;
TF2 *AliBlastwaveFit2D::fgFuncIntV2 = NULL;
const char *AliBlastwaveFit2D::fgParName[7] = {"T_{FO}","s_{2}","mean #rho_{0}","#rho_{2}","#gamma","mass","norm"};
Float_t AliBlastwaveFit2D::fgStartValues[5] = {0.1,0.057,1.2,0.025,1.1};
const Float_t AliBlastwaveFit2D::fgStepValues[5] = {0.001,0.001,0.001,0.001,0.001};
const Float_t AliBlastwaveFit2D::fgMinValues[5] = {0.0001,-0.45,0.01,-0.45,0.1};
const Float_t AliBlastwaveFit2D::fgMaxValues[5] = {1.0,0.45,10.0,0.45,5};

AliBlastwaveFit2D::AliBlastwaveFit2D(const char *name,Double_t mass) :
  AliBlastwaveFit(name,mass)
{
  if(!fgFuncIntYield){
     fgFuncIntYield = new TF2("fFuncInt2Yield",AliBlastwaveFit2D::FunctionIntYield,0,2*TMath::Pi(),0,1,7);
     fgFuncIntYield->SetNpx(100);
     fgFuncIntYield->SetNpy(100);
  }
  if(!fgFuncIntV2){
     fgFuncIntV2 = new TF2("fFuncInt2V2",AliBlastwaveFit2D::FunctionIntV2,0,2*TMath::Pi(),0,1,7);
     fgFuncIntV2->SetNpx(100);
     fgFuncIntV2->SetNpy(100);
  }

  Initialize();
}
//------------------------------------------------------------------------------
AliBlastwaveFit2D::AliBlastwaveFit2D() :
  AliBlastwaveFit()
{  
  if(!fgFuncIntYield){
     fgFuncIntYield = new TF2("fFuncInt2Yield",AliBlastwaveFit2D::FunctionIntYield,0,2*TMath::Pi(),0,1,7);
     fgFuncIntYield->SetNpx(100);
     fgFuncIntYield->SetNpy(100);
  }
  if(!fgFuncIntV2){
     fgFuncIntV2 = new TF2("fFuncInt2V2",AliBlastwaveFit2D::FunctionIntV2,0,2*TMath::Pi(),0,1,7);
     fgFuncIntV2->SetNpx(100);
     fgFuncIntV2->SetNpy(100);
  }
  Initialize();
}
//------------------------------------------------------------------------------
AliBlastwaveFit2D::~AliBlastwaveFit2D(){
}
//------------------------------------------------------------------------------
Double_t AliBlastwaveFit2D::FunctionIntYield(Double_t x[],Double_t par[]){
  // x[0] = phi
  // x[1] = r/R
  // par[0] = pt
  // par[1] = T
  // par[2] = s2
  // par[3] = rh0_0
  // par[4] = rho_a
  // par[5] = R_power
  // par[6] = mass 

  Double_t mt = TMath::Sqrt(par[0]*par[0]+par[6]*par[6]);
  Double_t rho = (par[3]+par[4]*TMath::Cos(2*x[0]))*TMath::Power(x[1],par[5]);
  Double_t alfat = (par[0]/par[1])*TMath::SinH(rho);
  Double_t betat = (mt/par[1])*TMath::CosH(rho);
  Double_t results;
  if(betat < 200) results = TMath::BesselI(0,alfat)*TMath::BesselK(1,betat)*(1+2*par[2]*TMath::Cos(2*x[0])) *x[1];
  else results = 0.5*TMath::Exp(alfat-betat)/sqrt(alfat*betat)*(1+2*par[2]*TMath::Cos(2*x[0])) *x[1];

  return results*par[0]*mt;
}
//------------------------------------------------------------------------------
Double_t AliBlastwaveFit2D::FunctionIntV2(Double_t x[],Double_t par[]){
  // x[0] = phi
  // x[1] = r/R
  // par[0] = pt
  // par[1] = T_fo
  // par[2] = s2
  // par[3] = rh0_0
  // par[4] = rho_a
  // par[5] = R_power
  // par[6] = mass

  Double_t mt = TMath::Sqrt(par[0]*par[0]+par[6]*par[6]);
  Double_t rho = (par[3]+par[4]*TMath::Cos(2*x[0]))*TMath::Power(x[1],par[5]);
  Double_t alfat = (par[0]/par[1])*TMath::SinH(rho);
  Double_t betat = (mt/par[1])*TMath::CosH(rho);
  Double_t results;
  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];
  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];

  return results*par[0]*mt;
}
//------------------------------------------------------------------------------
Int_t AliBlastwaveFit2D::SetParameter(Int_t ipar,Double_t val){
  Bool_t status = kTRUE;

  if(ipar < 0 || ipar > 4) return status;

  if(fFunctionYield){
     fFunctionYield->SetParameter(ipar,val);
     status = kFALSE;
  }
  if(fFunctionV2){
     fFunctionV2->SetParameter(ipar,val);
     status = kFALSE;
  }
  return status;
}
//------------------------------------------------------------------------------
Int_t AliBlastwaveFit2D::SetNormalization(){
  Bool_t status = kTRUE;

  TH1 *h = (TH1 *) GetSpectrumObjCopy();
  if(fFunctionYield && h){
      fFunctionYield->SetParameter(6,1.0);
      Double_t intH=0;
      Double_t intF=0;
      for (Int_t ibin = 1; ibin <= h->GetNbinsX(); ibin++) {
	  Double_t pt = h->GetBinCenter(ibin);
	  if(pt < GetMaxPt() && pt > GetMinPt() && TMath::Abs(h->GetBinContent(ibin)) > h->GetBinError(ibin)+0.0001){
	      Double_t width = h->GetBinWidth(ibin);
	      intH += h->GetBinContent(ibin)*width;
	      intF += EvalYield(pt)*width;
//  	      for(Double_t x = -width*0.45;x < width*0.5;x+=width*0.1)
//  		  intF += EvalYield(pt+x)*width*0.1;
	  }
      }

      Double_t norm = 1;
      if(intF > 0) norm = intH/intF;
      fFunctionYield->SetParameter(6,norm);
      status = kFALSE;
  }
  return status;
}
//------------------------------------------------------------------------------
void AliBlastwaveFit2D::SetMass(Double_t mass) {
  fMass=mass;
  if(fFunctionV2) fFunctionV2->SetParameter(5,fMass);
  if(fFunctionYield) fFunctionYield->SetParameter(5,fMass);   
}
//------------------------------------------------------------------------------
void AliBlastwaveFit2D::Initialize(){
  char nameFuncYield[100];
  char nameFuncV2[100];
  snprintf(nameFuncYield,100,"%sFuncYield",GetName());
  snprintf(nameFuncV2,100,"%sFuncV2",GetName());

  if(fFunctionV2) delete fFunctionV2;
  if(fFunctionYield) delete fFunctionYield;
  fFunctionV2 = new TF1(nameFuncV2,V2,0,6,6);
  fFunctionV2->SetParameter(5,fMass);  
  fFunctionYield = new TF1(nameFuncYield,Pt,0,6,7);
  fFunctionYield->SetParameter(5,fMass);
  fFunctionV2->SetNpx(200);
  fFunctionYield->SetNpx(200);

  for(Int_t i =0;i < 6;i++) fFunctionV2->SetParName(i,fgParName[i]);
  for(Int_t i =0;i < 7;i++) fFunctionYield->SetParName(i,fgParName[i]);
}
//------------------------------------------------------------------------------
Double_t AliBlastwaveFit2D::V2(Double_t x[],Double_t par[]){
  // par[0] = T_fo
  // par[1] = s2
  // par[2] = rh0_mean
  // par[3] = rho_v2
  // par[4] = R_power
  // par[5] = mass
  Double_t maxRho = par[2]*0.5*(par[4]+2)/(1+2*par[1]*par[3]);

  fgFuncIntYield->SetParameter(0,x[0]);
  fgFuncIntYield->SetParameter(1,par[0]);
  fgFuncIntYield->SetParameter(2,par[1]);
  fgFuncIntYield->SetParameter(3,maxRho);
  fgFuncIntYield->SetParameter(4,par[3]*maxRho*2);
  fgFuncIntYield->SetParameter(5,par[4]);
  fgFuncIntYield->SetParameter(6,par[5]);

  fgFuncIntV2->SetParameter(0,x[0]);
  fgFuncIntV2->SetParameter(1,par[0]);
  fgFuncIntV2->SetParameter(2,par[1]);
  fgFuncIntV2->SetParameter(3,maxRho);
  fgFuncIntV2->SetParameter(4,par[3]*maxRho*2);
  fgFuncIntV2->SetParameter(5,par[4]);
  fgFuncIntV2->SetParameter(6,par[5]);

  Double_t res = fgFuncIntV2->Integral(0.,2*TMath::Pi(),0.,1.);
  Double_t den = fgFuncIntYield->Integral(0.,2*TMath::Pi(),0.,1.);
  
  if(den == 0) return 0.0;

  return res/den;
}
//------------------------------------------------------------------------------
Double_t AliBlastwaveFit2D::Pt(Double_t x[],Double_t par[]){
  // par[0] = T_fo
  // par[1] = s2
  // par[2] = rh0_mean
  // par[3] = rho_v2
  // par[4] = R_power
  // par[5] = mass

  Double_t maxRho = par[2]*0.5*(par[4]+2)/(1+2*par[1]*par[3]);

  fgFuncIntYield->SetParameter(0,x[0]);
  fgFuncIntYield->SetParameter(1,par[0]);
  fgFuncIntYield->SetParameter(2,par[1]);
  fgFuncIntYield->SetParameter(3,maxRho);
  fgFuncIntYield->SetParameter(4,par[3]*maxRho*2);
  fgFuncIntYield->SetParameter(5,par[4]);
  fgFuncIntYield->SetParameter(6,par[5]);

  fgFuncIntV2->SetParameter(0,x[0]);
  fgFuncIntV2->SetParameter(1,par[0]);
  fgFuncIntV2->SetParameter(2,par[1]);
  fgFuncIntV2->SetParameter(3,maxRho);
  fgFuncIntV2->SetParameter(4,par[3]*maxRho*2);
  fgFuncIntV2->SetParameter(5,par[4]);
  fgFuncIntV2->SetParameter(6,par[5]);

  Double_t res = par[6]*1E+5*fgFuncIntYield->Integral(0.,2*TMath::Pi(),0.,1.);

  return res;
}
//------------------------------------------------------------------------------
Float_t AliBlastwaveFit2D::GetMeanBeta(){

  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);
  fbeta.SetNpx(1000);
  fbeta.SetNpy(1000);

  fbeta.SetParameter(0,fgFuncIntYield->GetParameter(3));
  fbeta.SetParameter(1,fgFuncIntYield->GetParameter(4));
  fbeta.SetParameter(2,fgFuncIntYield->GetParameter(5));
  fbeta.SetParameter(3,fgFuncIntYield->GetParameter(2));

  return fbeta.Integral(0.,2*TMath::Pi(),0.,1.)/TMath::Pi();
}
//------------------------------------------------------------------------------
Float_t AliBlastwaveFit2D::GetMeanBeta(Double_t par[]){
  // par[0] = T_fo
  // par[1] = s2
  // par[2] = rh0_mean
  // par[3] = rho_v2
  // par[4] = R_power
  // par[5] = mass

  Double_t maxRho = par[2]*0.5*(par[4]+2)/(1+2*par[1]*par[3]);

  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);
  fbeta.SetNpx(1000);
  fbeta.SetNpy(1000);

  fbeta.SetParameter(0,maxRho);
  fbeta.SetParameter(1,2*par[3]*maxRho);
  fbeta.SetParameter(2,par[4]);
  fbeta.SetParameter(3,par[1]);

  return fbeta.Integral(0.,2*TMath::Pi(),0.,1.)/TMath::Pi();
}
//------------------------------------------------------------------------------
void AliBlastwaveFit2D::SwitchOffFlow(TMinuit *m) const{
    m->FixParameter(kParS2);
    m->FixParameter(kParRho2);

    printf("AliBlastwaveFit2D: No flow histos -> Flow parameters put to zero\n");
}
Commit	Line	Data
3416bb60	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;
b4912369	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};
3416bb60	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(2x[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+2par[2]TMath::Cos(2x[0])) x[1];
	73	else results = 0.5TMath::Exp(alfat-betat)/sqrt(alfatbetat)(1+2par[2]TMath::Cos(2x[0])) *x[1];
	74
	75	return resultspar[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(2x[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(2x[0])TMath::BesselI(2,alfat)TMath::BesselK(1,betat)(1+2par[2]TMath::Cos(2x[0])) x[1];
95	else results = TMath::Cos(2x[0])0.5TMath::Exp(alfat-betat)/sqrt(alfatbetat)(1+2par[2]TMath::Cos(2x[0])) *x[1];
96
97	return resultspar[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 = -width0.45;x < width0.5;x+=width*0.1)
131	// intF += EvalYield(pt+x)width0.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+2par[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]maxRho2);
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]maxRho2);
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+2par[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]maxRho2);
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]maxRho2);
224	fgFuncIntV2->SetParameter(5,par[4]);
225	fgFuncIntV2->SetParameter(6,par[5]);
226
227	Double_t res = par[6]1E+5fgFuncIntYield->Integral(0.,2*TMath::Pi(),0.,1.);
228
229	return res;
230	}
231	//------------------------------------------------------------------------------
e931d3ab	232	Float_t AliBlastwaveFit2D::GetMeanBeta(){
3416bb60	233
	234	TF2 fbeta("fbeta","TMath::TanH(([0]+[1]TMath::Cos(2x))TMath::Power(y,[2]))y(1+2[3]cos(2x))",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	//------------------------------------------------------------------------------
e931d3ab	246	Float_t AliBlastwaveFit2D::GetMeanBeta(Double_t par[]){
3416bb60	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+2par[1]par[3]);
	255
	256	TF2 fbeta("fbeta","TMath::TanH(([0]+[1]TMath::Cos(2x))TMath::Power(y,[2]))y(1+2[3]cos(2x))",0,2*TMath::Pi(),0,1);
	257	fbeta.SetNpx(1000);
	258	fbeta.SetNpy(1000);
	259
	260	fbeta.SetParameter(0,maxRho);
	261	fbeta.SetParameter(1,2par[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	}