[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);

Int_t AliBlastwaveFitter::fgNparReal = 0;
Int_t AliBlastwaveFitter::fgNDGF = 0;
Int_t AliBlastwaveFitter::fgNfunctionCurrent=0;
AliBlastwaveFit *AliBlastwaveFitter::fgFuncC[50];
Int_t AliBlastwaveFitter::fgSpectraFlagC[50];
Int_t AliBlastwaveFitter::fgV2FlagC[50];
Double_t AliBlastwaveFitter::fgChi2 = 0;

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