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

#include<stdio.h>

#include"AliBlastwaveFitter.h"

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

ClassImp(AliBlastwaveFitter);

AliBlastwaveFitter::AliBlastwaveFitter(const char *name) :
    TNamed(name,name),
    fNfunction(0),
    fMinuit(new TMinuit()),
    fMinos(kFALSE)
{
}
//------------------------------------------------------------------------------
AliBlastwaveFitter::AliBlastwaveFitter() :
    TNamed("BlastwaveFitter","BlastwaveFitter"),
    fNfunction(0),
    fMinuit(new TMinuit()),
    fMinos(kFALSE)
{  
}
//------------------------------------------------------------------------------
AliBlastwaveFitter::~AliBlastwaveFitter(){
}
//------------------------------------------------------------------------------
Int_t AliBlastwaveFitter::CheckAvailability(){
    Int_t something = 2;
    Int_t consistency = 0;
    Int_t nparameters = -1;
    for(Int_t i=0;i < fNfunction;i++){
	if(!fFunc[i]){
	    fSpectraFlag[i] = 0;
	    fV2Flag[i] = 0;
	    continue;
	}
	if(fFunc[i]->GetSpectraFit() && fFunc[i]->GetSpectrumObj()){
	    TString classe(fFunc[i]->GetSpectrumObj()->ClassName());
	    if(classe.Contains("TH1")) fSpectraFlag[i] = 1;
	    else if(classe.Contains("TGraphE")) fSpectraFlag[i] = 2;
	    else fSpectraFlag[i] = 0;
	}
	else fSpectraFlag[i] = 0;
	if(fFunc[i]->GetV2Fit() && fFunc[i]->GetV2Obj()){
	    TString classe(fFunc[i]->GetV2Obj()->ClassName());
	    if(classe.Contains("TH1")) fV2Flag[i] = 1;
	    else if(classe.Contains("TGraphErr")) fV2Flag[i] = 2;
	    else if(classe.Contains("TGraphAsymmErr")) fV2Flag[i] = 3;
	    else fV2Flag[i] = 0;
	}
	else fV2Flag[i] = 0;

	if(fSpectraFlag[i] || fV2Flag[i]) something = 0;
	else continue;

	Int_t npar = fFunc[i]->GetNpar();
	if(nparameters == -1) nparameters = npar;
	else if(npar != nparameters) consistency = 1;
    }

    return (consistency + something);
}
//------------------------------------------------------------------------------
Int_t AliBlastwaveFitter::PrepareToFit(){
    Int_t check = CheckAvailability();

    Int_t npar=0;

    if(check%2){
	printf("Some problems with number of parameters are found\n");
	return check;
    }
    if((check/2)%2){
	printf("Nothing to fit\n");
	return check;
    }
    printf("Summary of fit functions and histos\n");

    Bool_t kFlowAvailable=kFALSE;

    for(Int_t i=0;i < fNfunction;i++){
	if(!fFunc[i]) continue;
	printf("------------------------------------------------------------------\n");
	printf("%2i) SpectraFlag = %i -- V2Flag = %i for \"%s\" (class:%s)  with mass = %6.3f GeV/c^2\n",i,fSpectraFlag[i],fV2Flag[i],fFunc[i]->GetName(),fFunc[i]->ClassName(),fFunc[i]->GetMass());
	if(fSpectraFlag[i]){
	    printf("   data spectra:\"%s\" (class:%s)\n",fFunc[i]->GetSpectrumObj()->GetName(),fFunc[i]->GetSpectrumObj()->ClassName());
	    npar = fFunc[i]->GetNpar();
	}
	if(fV2Flag[i]){
	    printf("   data v2     :\"%s\" (class:%s)\n",fFunc[i]->GetV2Obj()->GetName(),fFunc[i]->GetV2Obj()->ClassName());
	    npar = fFunc[i]->GetNpar();
	    kFlowAvailable=kTRUE;	
	}
    }
    printf("------------------------------------------------------------------\n");

    fgNparReal = npar;

    fgNDGF = 0;
    fgNfunctionCurrent = fNfunction;
    for(Int_t i=0;i < fNfunction;i++){
	fgFuncC[i] = fFunc[i];
	fgSpectraFlagC[i] = fSpectraFlag[i];
	fgV2FlagC[i] = fV2Flag[i];
    }
    // prepare minuit
    if(fgNparReal >= 25){
	fMinuit->BuildArrays(fgNparReal);
    }
    fMinuit->SetFCN(AliBlastwaveFitter::FCN);
    Double_t arglist[3];
    Int_t ierflg = 0;
    arglist[0] = 1;
    fMinuit->mnexcm("SET ERR", arglist, 1, ierflg);

    // par limits of the medium
    for(Int_t i=0;i < npar;i++){
      fMinuit->mnparm(i/*#par*/, fFunc[0]->GetParName(i)/*name*/, fFunc[0]->GetParStart(i)/*startvalue*/, fFunc[0]->GetParStep(i)/*step*/, fFunc[0]->GetParMin(i)/*min*/,fFunc[0]->GetParMax(i)/*max*/, ierflg);
    }
    
    if(!kFlowAvailable){ // if there is no flow histos
	for(Int_t i=0;i < fNfunction;i++){
	    if(fSpectraFlag[i]){
		fFunc[i]->SwitchOffFlow(fMinuit);
		i = fNfunction;
	    }
	}
    }

    // set strategy
    arglist[0] = 1; // !!!
    fMinuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
  
    return 0;
}
//------------------------------------------------------------------------------
Int_t AliBlastwaveFitter::Fit(){
    Int_t npar = fFunc[0]->GetNpar();
    /* start MIGRAD minimization */
    Double_t arglist[20];
    Int_t ierflg = 0;
    arglist[0] = 500000;
    arglist[1] = 1.;
    fMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);

    /* set strategy */
    arglist[0] = 2;
    fMinuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
    
    /* start MIGRAD minimization */
    arglist[0] = 500000;
    arglist[1] = 1.;
    fMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);

    /* start IMPROVE minimization */
    arglist[0] = 500000;
    fMinuit->mnexcm("IMPROVE", arglist, 1, ierflg);
    
    /* start MINOS */
    if(fMinos){
      arglist[0] = 500000;
      for(Int_t i=0;i < npar;i++){
  	arglist[i+1] = i;
      }
      fMinuit->mnexcm("MINOS", arglist, npar+1, ierflg);
    }

    /* print results */
    Double_t amin,edm,errdef;
    Int_t nvpar,nparx,icstat;
    fMinuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
    fMinuit->mnprin(4, amin);

    for(Int_t i=0;i < fgNfunctionCurrent;i++){
	fgFuncC[i]->Terminate();
    }

    return 0;
}
//------------------------------------------------------------------------------
void AliBlastwaveFitter::FCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){
    Double_t chi = 0., val, vale, pull,pt;
    Int_t iNorm = 0,nparFunc;

    if(gin[0]==1){}

//    printf("%i) %f %f %f %f %f %f %f\n",npar,par[0],par[1],par[2],par[3],par[4],par[5],par[6]);


    fgNDGF=0;  
  
    for(Int_t i=0;i < fgNfunctionCurrent;i++){
	if(!fgFuncC[i]) continue;
	
	// Set parameter in the current function
  	nparFunc = fgFuncC[i]->GetNpar();
	for(Int_t j=0;j<nparFunc;j++)
	    fgFuncC[i]->SetParameter(j,par[j]);

	// Set normalization in the current function
	if(fgSpectraFlagC[i]){// && nparFunc+iNorm < npar){
	    fgFuncC[i]->SetNormalization();//par[nparFunc+iNorm]);
	    iNorm++;
	}

	if(fgSpectraFlagC[i]){
	    // spectra comparison
	    TH1 *hForFit = (TH1 *) fgFuncC[i]->GetSpectrumObj();
	    TGraphErrors *gForFit = (TGraphErrors *) fgFuncC[i]->GetSpectrumObj();

	    if(fgSpectraFlagC[i]==1){
		for (Int_t ibin = 1; ibin <= hForFit->GetNbinsX(); ibin++) {
		    pt = hForFit->GetBinCenter(ibin);
		    if(pt < fgFuncC[i]->GetMaxPt() && pt > fgFuncC[i]->GetMinPt()){
			val = hForFit->GetBinContent(ibin);
			vale = hForFit->GetBinError(ibin);
			if(vale>0 && TMath::Abs(val) > vale+0.00001){
			  pull = (val - fgFuncC[i]->EvalYield(pt))/vale;
			  chi += pull * pull;
			  fgNDGF++;
			}
		    }
		} 
	    }
	    else if(fgSpectraFlagC[i]==2){
		for (Int_t ibin = 0; ibin < gForFit->GetN(); ibin++) {
		    pt = gForFit->GetX()[ibin];
		    if(pt < fgFuncC[i]->GetMaxPt() && pt > fgFuncC[i]->GetMinPt()){
			val = gForFit->GetY()[ibin];
			vale = gForFit->GetEY()[ibin];
			if(vale>0){
			  pull = (val - fgFuncC[i]->EvalYield(pt))/vale;
			  chi += pull * pull;
			  fgNDGF++;
			}
		    }		    
		}		
	    }
	}
	if(fgV2FlagC[i]){
	    // v2 comparison
	    TH1 *hForFit = (TH1 *) fgFuncC[i]->GetV2Obj();
	    TGraphErrors *gForFit = (TGraphErrors *) fgFuncC[i]->GetV2Obj();
	    TGraphAsymmErrors *gForFitA = (TGraphAsymmErrors *) fgFuncC[i]->GetV2Obj();

	    if(fgV2FlagC[i]==1){
		for (Int_t ibin = 1; ibin <= hForFit->GetNbinsX(); ibin++) {
		    pt = hForFit->GetBinCenter(ibin);
		    if(pt < fgFuncC[i]->GetMaxPt() && pt > fgFuncC[i]->GetMinPt()){
			val = hForFit->GetBinContent(ibin);
			vale = hForFit->GetBinError(ibin);
			if(vale>0 && TMath::Abs(val) > vale+0.00001){
			  pull = (val - fgFuncC[i]->EvalV2(pt))/vale;
			  chi += pull * pull;
			  fgNDGF++;
			}
		    }
		} 
	    }
	    else if(fgV2FlagC[i]==2){
		for (Int_t ibin = 0; ibin < gForFit->GetN(); ibin++) {
		    pt = gForFit->GetX()[ibin];
		    if(pt < fgFuncC[i]->GetMaxPt() && pt > fgFuncC[i]->GetMinPt()){
			val = gForFit->GetY()[ibin];
			vale = gForFit->GetEY()[ibin];
			if(vale>0){
			  pull = (val - fgFuncC[i]->EvalV2(pt))/vale;
			  chi += pull * pull;
			  fgNDGF++;
			}
		    }		    
		}		
	    }
	    else if(fgV2FlagC[i]==3){
		for (Int_t ibin = 0; ibin < gForFitA->GetN(); ibin++) {
		    pt = gForFit->GetX()[ibin];
		    if(pt < fgFuncC[i]->GetMaxPt() && pt > fgFuncC[i]->GetMinPt()){
			val = gForFit->GetY()[ibin];
			if(val - fgFuncC[i]->EvalV2(pt) > 0) vale = gForFit->GetEYlow()[ibin];
			else vale = gForFit->GetEYhigh()[ibin];
			if(vale>0){
			  pull = (val - fgFuncC[i]->EvalV2(pt))/vale;
			  chi += pull * pull;
			  fgNDGF++;
			}
		    }		    
		}		
	    }
	}
    }

    fgNDGF -= npar + iNorm;

    if(fgNDGF > 0) fgChi2 = chi/fgNDGF;

    f = chi;
    iflag = 0;
}
//------------------------------------------------------------------------------
TGraph*  AliBlastwaveFitter::DoContour(Int_t np,Int_t ip1,Int_t ip2,Float_t nsigma){

  TGraph *gCorr=NULL;
  
  Double_t par1,par2,err;
  fMinuit->GetParameter(ip1, par1,err);
  fMinuit->GetParameter(ip2, par2,err);

  Int_t trial = 0;
  Double_t scale = 1;
  while((! gCorr || gCorr->GetN() < np/4+5) && trial < 10){
    fMinuit->SetErrorDef(nsigma*nsigma*scale); //3 sigma contour
    gCorr = (TGraph *) fMinuit->Contour(np,ip1,ip2);
    scale *= 0.5;
    trial++;
  }
  scale *= 2;
  scale = TMath::Sqrt(scale);

  if(scale < 0.7 && gCorr){
    Double_t *x = gCorr->GetX();
    Double_t *y = gCorr->GetY();
    for(Int_t i=0;i < gCorr->GetN();i++){
      gCorr->SetPoint(i,par1+(x[i]-par1)/scale,par2+(y[i]-par2)/scale);
    }
  }

  if(trial > 1) printf("Contour realized with %3.1f sigma instead of %3.1f and then rescaled\n",nsigma*TMath::Sqrt(scale),nsigma);

  if(gCorr) gCorr->SetMarkerStyle(20);
  
  return gCorr;

}
 //------------------------------------------------------------------------------
TGraph* AliBlastwaveFitter::DoContourBetaT(Int_t np,Int_t iBoostOrBeta,Int_t iT,Float_t nsigma){


  TGraph *gCorr=NULL;
  
  Double_t par1,par2,err;
  fMinuit->GetParameter(iBoostOrBeta, par1,err);
  fMinuit->GetParameter(iT, par2,err);

  const Int_t maxnpar = 20;
  Double_t par[maxnpar];

  Int_t npar = fMinuit->GetNumPars();
  for(Int_t i=0;i < TMath::Min(npar,maxnpar);i++){
    fMinuit->GetParameter(i, par[i],err);
  }

  Int_t trial = 0;
  Double_t scale = 1;
  while((! gCorr || gCorr->GetN() < np/4+5) && trial < 10){
    fMinuit->SetErrorDef(nsigma*nsigma*scale); //3 sigma contour
    gCorr = (TGraph *) fMinuit->Contour(np,iBoostOrBeta,iT);
    scale *= 0.5;
    trial++;
  }
  scale *= 2;
  scale = TMath::Sqrt(scale);

  if(gCorr){
    Double_t *x = gCorr->GetX();
    Double_t *y = gCorr->GetY();
    for(Int_t i=0;i < gCorr->GetN();i++){
      Float_t meanboost = par1+(x[i]-par1)/scale;
      par[iBoostOrBeta] = meanboost;
      par[iT] = y[i];
      Float_t meanbeta = fgFuncC[0]->GetMeanBeta(par);
      gCorr->SetPoint(i,meanbeta,par2+(y[i]-par2)/scale);
    }
  }

  if(trial > 1) printf("Contour realized with %3.1f sigma instead of %3.1f and then rescaled\n",nsigma*TMath::Sqrt(scale),nsigma);

  if(gCorr) gCorr->SetMarkerStyle(20);
  
  return gCorr;

}
 //------------------------------------------------------------------------------
TGraph* AliBlastwaveFitter::ConvertContourFromBoostToBeta(TGraph *g,Int_t iBoostOrBeta,Int_t iT){

  if(!g) return NULL;

  TGraph *gCorr=new TGraph(g->GetN());
  
  const Int_t maxnpar = 20;
  Double_t par[maxnpar],err;

  Int_t npar = fMinuit->GetNumPars();
  for(Int_t i=0;i < TMath::Min(npar,maxnpar);i++){
    fMinuit->GetParameter(i, par[i],err);
  }

  if(gCorr){
    Double_t *x = g->GetX();
    Double_t *y = g->GetY();
    for(Int_t i=0;i < gCorr->GetN();i++){
      par[iBoostOrBeta] = x[i];
      par[iT] = y[i];
      Float_t meanbeta = fgFuncC[0]->GetMeanBeta(par);
      gCorr->SetPoint(i,meanbeta,y[i]);
    }
  }

  if(gCorr) gCorr->SetMarkerStyle(20);
  
  return gCorr;

}
Commit	Line	Data
3416bb60	1	#include<stdio.h>
	2
	3	#include"AliBlastwaveFitter.h"
	4
	5	#include"TMath.h"
	6	#include"TH1.h"
	7	#include"TGraphErrors.h"
	8	#include"TGraphAsymmErrors.h"
	9
	10	ClassImp(AliBlastwaveFitter);
	11
	12	AliBlastwaveFitter::AliBlastwaveFitter(const char *name) :
	13	TNamed(name,name),
	14	fNfunction(0),
	15	fMinuit(new TMinuit()),
	16	fMinos(kFALSE)
	17	{
	18	}
	19	//------------------------------------------------------------------------------
	20	AliBlastwaveFitter::AliBlastwaveFitter() :
	21	TNamed("BlastwaveFitter","BlastwaveFitter"),
	22	fNfunction(0),
	23	fMinuit(new TMinuit()),
	24	fMinos(kFALSE)
	25	{
	26	}
	27	//------------------------------------------------------------------------------
	28	AliBlastwaveFitter::~AliBlastwaveFitter(){
	29	}
	30	//------------------------------------------------------------------------------
	31	Int_t AliBlastwaveFitter::CheckAvailability(){
	32	Int_t something = 2;
	33	Int_t consistency = 0;
	34	Int_t nparameters = -1;
	35	for(Int_t i=0;i < fNfunction;i++){
	36	if(!fFunc[i]){
	37	fSpectraFlag[i] = 0;
	38	fV2Flag[i] = 0;
	39	continue;
	40	}
	41	if(fFunc[i]->GetSpectraFit() && fFunc[i]->GetSpectrumObj()){
	42	TString classe(fFunc[i]->GetSpectrumObj()->ClassName());
	43	if(classe.Contains("TH1")) fSpectraFlag[i] = 1;
	44	else if(classe.Contains("TGraphE")) fSpectraFlag[i] = 2;
	45	else fSpectraFlag[i] = 0;
	46	}
	47	else fSpectraFlag[i] = 0;
	48	if(fFunc[i]->GetV2Fit() && fFunc[i]->GetV2Obj()){
	49	TString classe(fFunc[i]->GetV2Obj()->ClassName());
	50	if(classe.Contains("TH1")) fV2Flag[i] = 1;
	51	else if(classe.Contains("TGraphErr")) fV2Flag[i] = 2;
	52	else if(classe.Contains("TGraphAsymmErr")) fV2Flag[i] = 3;
	53	else fV2Flag[i] = 0;
	54	}
	55	else fV2Flag[i] = 0;
	56
	57	if(fSpectraFlag[i] \|\| fV2Flag[i]) something = 0;
	58	else continue;
	59
	60	Int_t npar = fFunc[i]->GetNpar();
	61	if(nparameters == -1) nparameters = npar;
	62	else if(npar != nparameters) consistency = 1;
	63	}
	64
65	return (consistency + something);
66	}
67	//------------------------------------------------------------------------------
68	Int_t AliBlastwaveFitter::PrepareToFit(){
69	Int_t check = CheckAvailability();
70
71	Int_t npar=0;
72
73	if(check%2){
74	printf("Some problems with number of parameters are found\n");
75	return check;
76	}
77	if((check/2)%2){
78	printf("Nothing to fit\n");
79	return check;
80	}
81	printf("Summary of fit functions and histos\n");
82
83	Bool_t kFlowAvailable=kFALSE;
84
85	for(Int_t i=0;i < fNfunction;i++){
86	if(!fFunc[i]) continue;
87	printf("------------------------------------------------------------------\n");
88	printf("%2i) SpectraFlag = %i -- V2Flag = %i for \"%s\" (class:%s) with mass = %6.3f GeV/c^2\n",i,fSpectraFlag[i],fV2Flag[i],fFunc[i]->GetName(),fFunc[i]->ClassName(),fFunc[i]->GetMass());
89	if(fSpectraFlag[i]){
90	printf(" data spectra:\"%s\" (class:%s)\n",fFunc[i]->GetSpectrumObj()->GetName(),fFunc[i]->GetSpectrumObj()->ClassName());
91	npar = fFunc[i]->GetNpar();
92	}
93	if(fV2Flag[i]){
94	printf(" data v2 :\"%s\" (class:%s)\n",fFunc[i]->GetV2Obj()->GetName(),fFunc[i]->GetV2Obj()->ClassName());
95	npar = fFunc[i]->GetNpar();
96	kFlowAvailable=kTRUE;
97	}
98	}
99	printf("------------------------------------------------------------------\n");
100
101	fgNparReal = npar;
102
103	fgNDGF = 0;
104	fgNfunctionCurrent = fNfunction;
105	for(Int_t i=0;i < fNfunction;i++){
106	fgFuncC[i] = fFunc[i];
107	fgSpectraFlagC[i] = fSpectraFlag[i];
108	fgV2FlagC[i] = fV2Flag[i];
109	}
110	// prepare minuit
111	if(fgNparReal >= 25){
112	fMinuit->BuildArrays(fgNparReal);
113	}
114	fMinuit->SetFCN(AliBlastwaveFitter::FCN);
115	Double_t arglist[3];
116	Int_t ierflg = 0;
117	arglist[0] = 1;
118	fMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
119
120	// par limits of the medium
121	for(Int_t i=0;i < npar;i++){
122	fMinuit->mnparm(i/#par/, fFunc[0]->GetParName(i)/name/, fFunc[0]->GetParStart(i)/startvalue/, fFunc[0]->GetParStep(i)/step/, fFunc[0]->GetParMin(i)/min/,fFunc[0]->GetParMax(i)/max/, ierflg);
123	}
124
125	if(!kFlowAvailable){ // if there is no flow histos
126	for(Int_t i=0;i < fNfunction;i++){
127	if(fSpectraFlag[i]){
128	fFunc[i]->SwitchOffFlow(fMinuit);
129	i = fNfunction;
130	}
131	}
132	}
133
134	// set strategy
135	arglist[0] = 1; // !!!
136	fMinuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
137
138	return 0;
139	}
140	//------------------------------------------------------------------------------
141	Int_t AliBlastwaveFitter::Fit(){
142	Int_t npar = fFunc[0]->GetNpar();
143	/* start MIGRAD minimization */
144	Double_t arglist[20];
145	Int_t ierflg = 0;
146	arglist[0] = 500000;
147	arglist[1] = 1.;
148	fMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);
149
150	/* set strategy */
151	arglist[0] = 2;
152	fMinuit->mnexcm("SET STRATEGY", arglist, 1, ierflg);
153
154	/* start MIGRAD minimization */
155	arglist[0] = 500000;
156	arglist[1] = 1.;
157	fMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);
158
159	/* start IMPROVE minimization */
160	arglist[0] = 500000;
161	fMinuit->mnexcm("IMPROVE", arglist, 1, ierflg);
162
163	/* start MINOS */
164	if(fMinos){
165	arglist[0] = 500000;
166	for(Int_t i=0;i < npar;i++){
167	arglist[i+1] = i;
168	}
169	fMinuit->mnexcm("MINOS", arglist, npar+1, ierflg);
170	}
171
172	/* print results */
173	Double_t amin,edm,errdef;
174	Int_t nvpar,nparx,icstat;
175	fMinuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
176	fMinuit->mnprin(4, amin);
177
178	for(Int_t i=0;i < fgNfunctionCurrent;i++){
179	fgFuncC[i]->Terminate();
180	}
181
182	return 0;
183	}
184	//------------------------------------------------------------------------------
185	void AliBlastwaveFitter::FCN(Int_t &npar, Double_t gin, Double_t &f, Double_t par, Int_t iflag){
186	Double_t chi = 0., val, vale, pull,pt;
187	Int_t iNorm = 0,nparFunc;
188
189	if(gin[0]==1){}
190
191	// printf("%i) %f %f %f %f %f %f %f\n",npar,par[0],par[1],par[2],par[3],par[4],par[5],par[6]);
192
193
194	fgNDGF=0;
195
196	for(Int_t i=0;i < fgNfunctionCurrent;i++){
197	if(!fgFuncC[i]) continue;
198
199	// Set parameter in the current function
200	nparFunc = fgFuncC[i]->GetNpar();
201	for(Int_t j=0;j<nparFunc;j++)
202	fgFuncC[i]->SetParameter(j,par[j]);
203
204	// Set normalization in the current function
205	if(fgSpectraFlagC[i]){// && nparFunc+iNorm < npar){
206	fgFuncC[i]->SetNormalization();//par[nparFunc+iNorm]);
207	iNorm++;
208	}
209
210	if(fgSpectraFlagC[i]){
211	// spectra comparison
212	TH1 hForFit = (TH1 ) fgFuncC[i]->GetSpectrumObj();
213	TGraphErrors gForFit = (TGraphErrors ) fgFuncC[i]->GetSpectrumObj();
214
215	if(fgSpectraFlagC[i]==1){
216	for (Int_t ibin = 1; ibin <= hForFit->GetNbinsX(); ibin++) {
217	pt = hForFit->GetBinCenter(ibin);
218	if(pt < fgFuncC[i]->GetMaxPt() && pt > fgFuncC[i]->GetMinPt()){
219	val = hForFit->GetBinContent(ibin);
220	vale = hForFit->GetBinError(ibin);
221	if(vale>0 && TMath::Abs(val) > vale+0.00001){
222	pull = (val - fgFuncC[i]->EvalYield(pt))/vale;
223	chi += pull * pull;
224	fgNDGF++;
225	}
226	}
227	}
228	}
229	else if(fgSpectraFlagC[i]==2){
230	for (Int_t ibin = 0; ibin < gForFit->GetN(); ibin++) {
231	pt = gForFit->GetX()[ibin];
232	if(pt < fgFuncC[i]->GetMaxPt() && pt > fgFuncC[i]->GetMinPt()){
233	val = gForFit->GetY()[ibin];
234	vale = gForFit->GetEY()[ibin];
235	if(vale>0){
236	pull = (val - fgFuncC[i]->EvalYield(pt))/vale;
237	chi += pull * pull;
238	fgNDGF++;
239	}
240	}
241	}
242	}
243	}
244	if(fgV2FlagC[i]){
245	// v2 comparison
246	TH1 hForFit = (TH1 ) fgFuncC[i]->GetV2Obj();
247	TGraphErrors gForFit = (TGraphErrors ) fgFuncC[i]->GetV2Obj();
248	TGraphAsymmErrors gForFitA = (TGraphAsymmErrors ) fgFuncC[i]->GetV2Obj();
249
250	if(fgV2FlagC[i]==1){
251	for (Int_t ibin = 1; ibin <= hForFit->GetNbinsX(); ibin++) {
252	pt = hForFit->GetBinCenter(ibin);
253	if(pt < fgFuncC[i]->GetMaxPt() && pt > fgFuncC[i]->GetMinPt()){
254	val = hForFit->GetBinContent(ibin);
255	vale = hForFit->GetBinError(ibin);
256	if(vale>0 && TMath::Abs(val) > vale+0.00001){
257	pull = (val - fgFuncC[i]->EvalV2(pt))/vale;
258	chi += pull * pull;
259	fgNDGF++;
260	}
261	}
262	}
263	}
264	else if(fgV2FlagC[i]==2){
265	for (Int_t ibin = 0; ibin < gForFit->GetN(); ibin++) {
266	pt = gForFit->GetX()[ibin];
267	if(pt < fgFuncC[i]->GetMaxPt() && pt > fgFuncC[i]->GetMinPt()){
268	val = gForFit->GetY()[ibin];
269	vale = gForFit->GetEY()[ibin];
270	if(vale>0){
271	pull = (val - fgFuncC[i]->EvalV2(pt))/vale;
272	chi += pull * pull;
273	fgNDGF++;
274	}
275	}
276	}
277	}
278	else if(fgV2FlagC[i]==3){
279	for (Int_t ibin = 0; ibin < gForFitA->GetN(); ibin++) {
280	pt = gForFit->GetX()[ibin];
281	if(pt < fgFuncC[i]->GetMaxPt() && pt > fgFuncC[i]->GetMinPt()){
282	val = gForFit->GetY()[ibin];
283	if(val - fgFuncC[i]->EvalV2(pt) > 0) vale = gForFit->GetEYlow()[ibin];
284	else vale = gForFit->GetEYhigh()[ibin];
285	if(vale>0){
286	pull = (val - fgFuncC[i]->EvalV2(pt))/vale;
287	chi += pull * pull;
288	fgNDGF++;
289	}
290	}
291	}
292	}
293	}
294	}
295
296	fgNDGF -= npar + iNorm;
297
298	if(fgNDGF > 0) fgChi2 = chi/fgNDGF;
299
300	f = chi;
301	iflag = 0;
302	}
303	//------------------------------------------------------------------------------
304	TGraph* AliBlastwaveFitter::DoContour(Int_t np,Int_t ip1,Int_t ip2,Float_t nsigma){
305
306	TGraph *gCorr=NULL;
307
308	Double_t par1,par2,err;
309	fMinuit->GetParameter(ip1, par1,err);
310	fMinuit->GetParameter(ip2, par2,err);
311
312	Int_t trial = 0;
313	Double_t scale = 1;
314	while((! gCorr \|\| gCorr->GetN() < np/4+5) && trial < 10){
315	fMinuit->SetErrorDef(nsigmansigmascale); //3 sigma contour
316	gCorr = (TGraph *) fMinuit->Contour(np,ip1,ip2);
317	scale *= 0.5;
318	trial++;
319	}
320	scale *= 2;
321	scale = TMath::Sqrt(scale);
322
323	if(scale < 0.7 && gCorr){
324	Double_t *x = gCorr->GetX();
325	Double_t *y = gCorr->GetY();
326	for(Int_t i=0;i < gCorr->GetN();i++){
327	gCorr->SetPoint(i,par1+(x[i]-par1)/scale,par2+(y[i]-par2)/scale);
328	}
329	}
330
331	if(trial > 1) printf("Contour realized with %3.1f sigma instead of %3.1f and then rescaled\n",nsigma*TMath::Sqrt(scale),nsigma);
332
333	if(gCorr) gCorr->SetMarkerStyle(20);
334
335	return gCorr;
336
337	}
338	//------------------------------------------------------------------------------
339	TGraph* AliBlastwaveFitter::DoContourBetaT(Int_t np,Int_t iBoostOrBeta,Int_t iT,Float_t nsigma){
340
341
342	TGraph *gCorr=NULL;
343
344	Double_t par1,par2,err;
345	fMinuit->GetParameter(iBoostOrBeta, par1,err);
346	fMinuit->GetParameter(iT, par2,err);
347
348	const Int_t maxnpar = 20;
349	Double_t par[maxnpar];
350
351	Int_t npar = fMinuit->GetNumPars();
352	for(Int_t i=0;i < TMath::Min(npar,maxnpar);i++){
353	fMinuit->GetParameter(i, par[i],err);
354	}
355
356	Int_t trial = 0;
357	Double_t scale = 1;
358	while((! gCorr \|\| gCorr->GetN() < np/4+5) && trial < 10){
359	fMinuit->SetErrorDef(nsigmansigmascale); //3 sigma contour
360	gCorr = (TGraph *) fMinuit->Contour(np,iBoostOrBeta,iT);
361	scale *= 0.5;
362	trial++;
363	}
364	scale *= 2;
365	scale = TMath::Sqrt(scale);
366
367	if(gCorr){
368	Double_t *x = gCorr->GetX();
369	Double_t *y = gCorr->GetY();
370	for(Int_t i=0;i < gCorr->GetN();i++){
371	Float_t meanboost = par1+(x[i]-par1)/scale;
372	par[iBoostOrBeta] = meanboost;
373	par[iT] = y[i];
374	Float_t meanbeta = fgFuncC[0]->GetMeanBeta(par);
375	gCorr->SetPoint(i,meanbeta,par2+(y[i]-par2)/scale);
376	}
377	}
378
379	if(trial > 1) printf("Contour realized with %3.1f sigma instead of %3.1f and then rescaled\n",nsigma*TMath::Sqrt(scale),nsigma);
380
381	if(gCorr) gCorr->SetMarkerStyle(20);
382
383	return gCorr;
384
385	}
386	//------------------------------------------------------------------------------
387	TGraph* AliBlastwaveFitter::ConvertContourFromBoostToBeta(TGraph *g,Int_t iBoostOrBeta,Int_t iT){
388
389	if(!g) return NULL;
390
391	TGraph *gCorr=new TGraph(g->GetN());
392
393	const Int_t maxnpar = 20;
394	Double_t par[maxnpar],err;
395
396	Int_t npar = fMinuit->GetNumPars();
397	for(Int_t i=0;i < TMath::Min(npar,maxnpar);i++){
398	fMinuit->GetParameter(i, par[i],err);
399	}
400
401	if(gCorr){
402	Double_t *x = g->GetX();
403	Double_t *y = g->GetY();
404	for(Int_t i=0;i < gCorr->GetN();i++){
405	par[iBoostOrBeta] = x[i];
406	par[iT] = y[i];
407	Float_t meanbeta = fgFuncC[0]->GetMeanBeta(par);
408	gCorr->SetPoint(i,meanbeta,y[i]);
409	}
410	}
411
412	if(gCorr) gCorr->SetMarkerStyle(20);
413
414	return gCorr;
415
416	}