Added more functions and fixed blast wave
authormfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 14 Oct 2011 13:41:01 +0000 (13:41 +0000)
committermfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 14 Oct 2011 13:41:01 +0000 (13:41 +0000)
PWG2/SPECTRA/Fit/AliBWFunc.cxx
PWG2/SPECTRA/Fit/AliBWFunc.h
PWG2/SPECTRA/Fit/FitParticle.C

index 0063de4841b802970aab74a887e05803976ea17e..00c43409ff7c4ec7335b5e68b35e1c53f5e2a522 100644 (file)
@@ -43,6 +43,7 @@ TF1 * AliBWFunc::GetHistoFunc(TH1 * h, const char * name) {
 
   // Regardless of the variable type, this returns a function made
   // from the histo * a multiplicative normalization.
+  // This uses a bad hack...
 
   fLastFunc = new TF1 (name, StaticHistoFunc, 0.0, 10, 2);
   fLastFunc->SetParameter(0,1);
@@ -52,20 +53,36 @@ TF1 * AliBWFunc::GetHistoFunc(TH1 * h, const char * name) {
   return fLastFunc;
   
 
+
+}
+TF1 * AliBWFunc::GetGraphFunc(TGraph * g, const char * name) {
+
+  // Regardless of the variable type, this returns a function made
+  // from the graph * a multiplicative normalization.
+  // This uses a bad hack...
+
+  fLastFunc = new TF1 (name, StaticHistoFunc, 0.0, 10, 2);
+  fLastFunc->SetParameter(0,1);
+  fLastFunc->FixParameter(1,Double_t(Long64_t(g)));
+  fLastFunc->SetParNames("norm", "pointer to histo");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+  
+
 }
 
 
 TF1 * AliBWFunc::GetBGBW(Double_t mass, Double_t beta, Double_t T,
-             Double_t norm, const char * name){
+                        Double_t n, Double_t norm, const char * name){
 
   // Boltzmann-Gibbs blast wave
 
   switch (fVarType) {
   case kdNdpt:
-    return GetBGBWdNdptTimesPt(mass,beta,T,norm,name);
+    return GetBGBWdNdptTimesPt(mass,beta,T,n,norm,name);
     break;
   case kOneOverPtdNdpt:
-    return GetBGBWdNdpt(mass,beta,T,norm,name);
+    return GetBGBWdNdpt(mass,beta,T,n,norm,name);
     break;
   case kOneOverMtdNdmt:
     AliFatal("Not implemented");
@@ -144,6 +161,52 @@ TF1 * AliBWFunc::GetMTExp(Double_t mass, Double_t T, Double_t norm, const char *
 
 }
 
+
+TF1 * AliBWFunc::GetBoseEinstein(Double_t mass, Double_t T, Double_t norm, const char * name){
+
+  // Bose einstein
+  switch (fVarType) {
+  case kdNdpt:
+    return GetBoseEinsteindNdptTimesPt(mass,T,norm,name);
+    break;
+  case kOneOverPtdNdpt:
+    return GetBoseEinsteindNdpt(mass,T,norm,name);
+    break;
+  case kOneOverMtdNdmt:
+    AliFatal("Not implemented");
+    break;
+  default:
+    AliFatal("Not implemented");
+  }
+  
+  return 0;
+
+
+}
+
+TF1 * AliBWFunc::GetFermiDirac(Double_t mass, Double_t T, Double_t norm, const char * name){
+
+  // Simple exponential in 1/mt*MT
+  switch (fVarType) {
+  case kdNdpt:
+    return GetFermiDiracdNdptTimesPt(mass,T,norm,name);
+    break;
+  case kOneOverPtdNdpt:
+    return GetFermiDiracdNdpt(mass,T,norm,name);
+    break;
+  case kOneOverMtdNdmt:
+    AliFatal("Not implemented");
+    break;
+  default:
+    AliFatal("Not implemented");
+  }
+  
+  return 0;
+
+
+}
+
+
 TF1 * AliBWFunc::GetPTExp(Double_t T, Double_t norm, const char * name){
 
   // Simple exponential in 1/mt*MT
@@ -269,18 +332,26 @@ Double_t AliBWFunc::StaticHistoFunc(const double * x, const double* p){
 
   double norm = p[0];
   
-  TH1 * h     = (TH1*) Long64_t(p[1]);
+  TObject * h     = (TObject*) Long64_t(p[1]);
 
 //    Int_t bin = h->FindBin(x[0]);
 //    double value = h->GetBinContent(bin);
 
-  if (h->FindBin(x[0]) > h->GetNbinsX()) return 0;
 
   // static TH1 * oldptr = 0;
   // static TSpline3 * spl = 0;
   // if (h!=oldptr) {
   // FIXME: recheck static pointers
-  TSpline3 * spl  = new TSpline3(h);
+  TSpline3 * spl  = 0;
+  if(h->InheritsFrom("TH1")) {
+    if ( ((TH1*)h)->FindBin(x[0]) > ((TH1*)h)->GetNbinsX()) return 0;
+    spl= new TSpline3((TH1*)h);
+  }
+  else if(h->InheritsFrom("TGraph")) spl= new TSpline3("fGraph",(TGraph*)h);
+  else {
+    Printf("AliBWFunc::StaticHistoFunc: Unsupported type");
+    return 0;
+  }
     //  }
   double value = spl->Eval(x[0]);
   delete spl;
@@ -329,21 +400,36 @@ Double_t AliBWFunc::StaticUA1Func(const double * x, const double* p) {
 
 Double_t AliBWFunc::IntegrandBG(const double * x, const double* p){
   // integrand for boltzman-gibbs blast wave
+     // x[0] -> r (radius)
+     // p[0] -> mass
+     // p[1] -> pT (transverse momentum)
+     // p[2] -> beta_max (surface velocity)
+     // p[3] -> T (freezout temperature)
+     // p[4] -> n (velocity profile)
+
 
   double x0 = x[0]; 
   
-  double mass = p[0];
-  double pT   = p[1];
-  double beta = p[2];
-  double temp    = p[3];
-  
+  double mass     = p[0];
+  double pT       = p[1];
+  double beta_max = p[2];
+  double temp     = p[3];
+  Double_t n      = p[4];
+
+  // Keep beta within reasonable limits
+  Double_t beta = beta_max * TMath::Power(x0, n);
+  if (beta > 0.9999999999999999) beta = 0.9999999999999999;
+
   double mT      = TMath::Sqrt(mass*mass+pT*pT);
 
-  double rho0   = TMath::ATanH(beta*x0);  
+  double rho0   = TMath::ATanH(beta);  
   double arg00 = pT*TMath::SinH(rho0)/temp;
+  if (arg00 > 700.) arg00 = 700.; // avoid FPE
   double arg01 = mT*TMath::CosH(rho0)/temp;
   double f0 = x0*mT*TMath::BesselI0(arg00)*TMath::BesselK1(arg01);
 
+  //  printf("r=%f, pt=%f, beta_max=%f, temp=%f, n=%f, mt=%f, beta=%f, rho=%f, argI0=%f, argK1=%f\n", x0, pT, beta_max, temp, n, mT, beta, rho0, arg00, arg01);
+
   return f0;
 }
 
@@ -356,17 +442,20 @@ Double_t AliBWFunc::StaticBGdNdPt(const double * x, const double* p) {
   double pT = x[0];;
   
 
-  double mass = p[0];
-  double beta = p[1];
+  double mass    = p[0];
+  double beta    = p[1];
   double temp    = p[2];
+  double n       = p[3];
+  double norm    = p[4];
 
   static TF1 * fIntBG = 0;
   if(!fIntBG)
-    fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 4);
+    fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 5);
 
-  fIntBG->SetParameters(mass, pT, beta, temp);
+  fIntBG->SetParameters(mass, pT, beta, temp,n);
   double result = fIntBG->Integral(0,1);
-  return result*p[3];//*1e30;;
+  //  printf ("[%4.4f], Int :%f\n", pT, result);
+  return result*norm;//*1e30;;
 
 }
 
@@ -377,13 +466,14 @@ Double_t AliBWFunc::StaticBGdNdPtTimesPt(const double * x, const double* p) {
 
 
 TF1 * AliBWFunc::GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t temp,
-                           Double_t norm, const char * name){
+                             Double_t n, Double_t norm, const char * name){
   
   // BGBW 1/pt dNdpt
 
-  fLastFunc = new TF1 (name, StaticBGdNdPt, 0.0, 10, 4);
-  fLastFunc->SetParameters(mass,beta,temp,norm);    
-  fLastFunc->SetParNames("mass", "#beta", "T", "norm");
+  fLastFunc = new TF1 (name, StaticBGdNdPt, 0.0, 10, 5);
+  fLastFunc->SetParameters(mass,beta,temp,n,norm);    
+  fLastFunc->FixParameter(0,mass);
+  fLastFunc->SetParNames("mass", "#beta", "T", "n", "norm");
   fLastFunc->SetLineWidth(fLineWidth);
   return fLastFunc;
   
@@ -484,14 +574,15 @@ TF1 * AliBWFunc::GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t temp,
 
 // Times Pt funcs
 // Boltzmann-Gibbs Blast Wave
-TF1 * AliBWFunc::GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp,
+TF1 * AliBWFunc::GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t temp, Double_t n,
                                     Double_t norm, const char * name){
 
   // BGBW, dNdpt
 
-  fLastFunc = new TF1 (name, StaticBGdNdPtTimesPt, 0.0, 10, 4);
-  fLastFunc->SetParameters(mass,beta,temp,norm);    
-  fLastFunc->SetParNames("mass", "#beta", "temp", "norm");
+  fLastFunc = new TF1 (name, StaticBGdNdPtTimesPt, 0.0, 10, 5);
+  fLastFunc->SetParameters(mass,beta,temp,n,norm);    
+  fLastFunc->FixParameter(0,mass);
+  fLastFunc->SetParNames("mass", "#beta", "temp", "n", "norm");
   fLastFunc->SetLineWidth(fLineWidth);
   return fLastFunc;
 
@@ -531,6 +622,37 @@ TF1 * AliBWFunc::GetMTExpdNdptTimesPt(Double_t mass, Double_t temp, Double_t nor
 
 }
 
+TF1 * AliBWFunc::GetBoseEinsteindNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
+
+  // Bose einstein distribution as a function of dNdpt
+  char formula[500];
+  snprintf(formula,500,"[0]*x*1./(exp(sqrt(x**2+%f**2)/[1])-1)", mass);
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, temp);
+  fLastFunc->SetParLimits(1, 0.01, 10);
+  fLastFunc->SetParNames("norm", "T");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+TF1 * AliBWFunc::GetFermiDiracdNdptTimesPt(Double_t mass, Double_t temp, Double_t norm, const char * name){
+
+  // Bose einstein distribution as a function of dNdpt
+  char formula[500];
+  snprintf(formula,500,"[0]*x*1./(exp(sqrt(x**2+%f**2)/[1])+1)", mass);
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, temp);
+  fLastFunc->SetParLimits(1, 0.01, 10);
+  fLastFunc->SetParNames("norm", "T");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+
 
 TF1 * AliBWFunc::GetPTExpdNdptTimesPt(Double_t temp, Double_t norm, const char * name){
 
@@ -748,6 +870,30 @@ TF1 * AliBWFunc::GetMTExpdNdpt(Double_t mass, Double_t temp, Double_t norm, cons
   return fLastFunc;
 }
 
+TF1 * AliBWFunc::GetBoseEinsteindNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
+  // bose einstein
+  char formula[500];
+  snprintf(formula,500,"[0]*1./(exp(sqrt(x**2+%f**2)/[1])-1)", mass);
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, temp);
+  fLastFunc->SetParLimits(1, 0.01, 10);
+  fLastFunc->SetParNames("norm", "T");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+}
+
+TF1 * AliBWFunc::GetFermiDiracdNdpt(Double_t mass, Double_t temp, Double_t norm, const char * name){
+  // bose einstein
+  char formula[500];
+  snprintf(formula,500,"[0]*1./(exp(sqrt(x**2+%f**2)/[1])+1)", mass);
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, temp);
+  fLastFunc->SetParLimits(1, 0.01, 10);
+  fLastFunc->SetParNames("norm", "T");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+}
+
 
 // // Simple tsallis (a la CMS)
 // TF1 * AliBWFunc::GetTsallisdNdpt(Double_t mass, Double_t temp, Double_t q, Double_t norm, const char * name){
index 27a2f7aadfcd75ac3cb585f1553ceeb6268939c9..d774687fad9aa5ef74321818380ab2b457fde7ac 100644 (file)
@@ -23,7 +23,7 @@
 
 class TF1;
 class TH1;
-
+class TGraph;
 #endif
 
 
@@ -38,9 +38,10 @@ public:
   AliBWFunc();
   ~AliBWFunc();
 
-  // Boltzmann-Gibbs blast wave
-  TF1 * GetBGBW(Double_t mass, Double_t beta, Double_t T,
+  // Boltzmann-Gibbs Blast Wave
+  TF1 * GetBGBW(Double_t mass, Double_t beta, Double_t T, Double_t n,
                Double_t norm, const char * name = "fBGBW");
+
   
   // Boltzmann
   TF1 * GetBoltzmann(Double_t mass, Double_t T, Double_t norm, const char * name ="fBoltzmann");
@@ -50,10 +51,10 @@ public:
                     Double_t norm, Double_t ymax = 0.5, const char * name = "fTsallisBW");
 
   // Simple exponential in 1/mt*dNdmt
-  TF1 * GetMTExp(Double_t mass, Double_t T, Double_t norm, const char * name ="fExp");
+  TF1 * GetMTExp(Double_t mass, Double_t T, Double_t norm, const char * name ="fMtExp");
 
   // Simple exponential in 1/pt*dNdpt
-  TF1 * GetPTExp(Double_t T, Double_t norm, const char * name ="fExp");
+  TF1 * GetPTExp(Double_t T, Double_t norm, const char * name ="fPtExp");
 
   // Tsallis (no BW, a la CMS)
   TF1 * GetTsallis(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name="fTsallis") 
@@ -68,8 +69,18 @@ public:
   // Function derived from a histo
   TF1 * GetHistoFunc(TH1 * h, const char * name = "fHisto");
 
+  // Function derived from a graph
+  TF1 * GetGraphFunc(TGraph * h, const char * name = "fHisto");
+
   // Power law
   TF1 * GetPowerLaw(Double_t pt0, Double_t n, Double_t norm, const char * name="fPowerLaw");
+  
+
+  // Bose-Einstein
+  TF1 * GetBoseEinstein(Double_t mass, Double_t T, Double_t norm, const char * name="fBoseEinstein");
+
+  // Fermi-Dirac
+  TF1 * GetFermiDirac(Double_t mass, Double_t T, Double_t norm, const char * name="fFermiDirac");
 
 
   void SetVarType(VarType_t tp) {fVarType=tp;}
@@ -78,9 +89,10 @@ protected:
 
   // dNdpt here means 1/pt dN/dpt
   
-  // Boltzmann-Gibbs Blast Wave
-  TF1 * GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t T,
-                    Double_t norm, const char * name = "fBGBW");
+
+  // Boltzmann-Gibbs blast wave
+  TF1 * GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t temp,
+                    Double_t n, Double_t norm, const char * name= "fBGBW");
 
   // Tsallis blast wave
   TF1 * GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t T, Double_t q,
@@ -89,6 +101,13 @@ protected:
   // Simple exponential in 1/mt*MT
   TF1 * GetMTExpdNdpt(Double_t mass, Double_t T, Double_t norm, const char * name ="fExp");
 
+  // Bose-Einstein
+  TF1 * GetBoseEinsteindNdpt(Double_t mass, Double_t T, Double_t norm, const char * name="fBoseEinstein");
+
+  // Fermi-Dirac
+  TF1 * GetFermiDiracdNdpt(Double_t mass, Double_t T, Double_t norm, const char * name="fFermiDirac");
+
+
   // Tsallis (no BW, a la CMS)
   TF1 * GetTsallisdNdpt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name="fTsallis");
 
@@ -104,7 +123,7 @@ protected:
   // TimesPt means dNdpt
 
   // Boltzmann-Gibbs Blast Wave
-  TF1 * GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T,
+  TF1 * GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T, Double_t n,
                            Double_t norm, const char * name = "fBGBWTimesPt");
 
   // Tsallis blast wave
@@ -122,6 +141,13 @@ protected:
   // Simple exponential in 1/mt*dNdmT
   TF1 * GetMTExpdNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name ="fMtExpTimesPt");
 
+  // Bose-Einstein
+  TF1 * GetBoseEinsteindNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name="fBoseEinstein");
+
+  // Fermi-Dirac
+  TF1 * GetFermiDiracdNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name="fFermiDirac");
+
+
   // Simple exponential in 1/mp*dNdpT
   TF1 * GetPTExpdNdptTimesPt(Double_t T, Double_t norm, const char * name ="fPtExpTimesPt");
 
index d7ff6ad62121a7e4d9eeb74848b591395eb1ee5d..236854e4f1bb6dafc9690e841f37df874a994b4b 100644 (file)
 #include "TSystem.h"
 #include "THashList.h"
 #include "TMinuit.h"
+#include "TLatex.h"
 
 using namespace std;
 #endif
 
-enum {kFitExpPt, kFitLevi, fFitExpMt, kFitBoltzmann, kFitBlastWave};
+enum {kFitExpPt, kFitLevi, fFitExpMt, kFitBoltzmann, kFitBlastWave, kFitBoseEinstein, kFitFermiDirac};
 
-void FitParticle(const char * file, const char * histo, const char * partName,  const char * listname=0, Float_t min = 0, Float_t max =3, Float_t scaleHisto = -1., Int_t fitFunc = kFitLevi, Int_t vartype = AliBWFunc::kdNdpt) {
+void FitParticle(const char * file, const char * histo, const char * partName,  const char * listname=0, Float_t min = 0, Float_t max =3, Float_t scaleHisto = -1., Int_t fitFunc = kFitLevi, Int_t vartype = AliBWFunc::kdNdpt, const char * fileOut = 0, Bool_t wait = 0) {
 
   // Generic Macro to fit any particle using the PWG2/SPECTRA/Fit macros
   //
   // You have to provide:
-  // - file: file name 
+  // - file: file name. If not set, uses gFile
   // - histo: histogram name name (assumend to be in the main folder of the file)
   // - listname: if different from 0, histo is looked in this tlist
   // - partName: it is used to get the mass of the particle from
@@ -78,28 +79,69 @@ void FitParticle(const char * file, const char * histo, const char * partName,
   Double_t mass = TDatabasePDG::Instance()->GetParticle(partName)->Mass();
   cout << "Fitting ["<<partName<<"], mass: " << mass << endl;
   
+
+  AliLatexTable * table;
+
   TF1* func = 0;
   Int_t normPar = -1;
+  Int_t slopePar=0;
   if (fitFunc == kFitLevi)   {
     func = fm->GetLevi (mass, 0.4, 20,3);
-    func->SetParLimits(1,1,100);
+    func->SetParLimits(1,0.0001,20000);
+    table = new AliLatexTable(9,"c|cccccccc");
+    table->InsertCustomRow("Part & Func & Yield & FD Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
     normPar = 0;
+    slopePar = 2;
   }
   if (fitFunc == kFitExpPt)  {
     func = fm->GetPTExp(0.2, 20);
+    table = new AliLatexTable(9,"c|cccccccc");
+    table->InsertCustomRow("Part & Func & Yield & pT Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
+    normPar = 0;
+    slopePar = 1;
   }
   if (fitFunc == fFitExpMt)  {
     func = fm->GetMTExp(mass,0.2, 20);
+    table = new AliLatexTable(9,"c|cccccccc");
+    table->InsertCustomRow("Part & Func & Yield & mT Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
+    normPar = 0;
+    slopePar = 1;
   }
   if (fitFunc == kFitBoltzmann)  {
     func = fm->GetBoltzmann(mass, 0.2, 20);
+    table = new AliLatexTable(9,"c|cccccccc");
+    table->InsertCustomRow("Part & Func & Yield & Bz Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
+    normPar = 0;
+    slopePar = 1;
+
   }
   if (fitFunc == kFitBlastWave)  {
-    func = fm->GetBGBW(mass,0.6,0.3, 20);
+    func = fm->GetBGBW(mass,0.6,0.3, 1, 1e5);// beta, T, n, norm 
+    table = new AliLatexTable(11,"cc|ccccccccccc");
+    table->InsertCustomRow("Part & Func & Yield & T & beta & n & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
+    func->SetParLimits(1, 0.1, 0.99);
+    func->SetParLimits(2, 0.01, 1);
+    func->SetParLimits(3, 0.01, 2);
+    normPar = 4;
+    slopePar = 2;
+  }
+  if (fitFunc == kFitBoseEinstein)  {
+    func = fm->GetBoseEinstein(mass,0.3,20);
+    table = new AliLatexTable(9,"c|cccccccc");
+    table->InsertCustomRow("Part & Func & Yield & BE Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
+    normPar = 0;
+    slopePar = 1;
+  }
+  if (fitFunc == kFitFermiDirac)  {
+    func = fm->GetFermiDirac(mass,0.3,20);
+    table = new AliLatexTable(9,"c|cccccccc");
+    table->InsertCustomRow("Part & Func & Yield & FD Slope & $\\Chi^2$/NDF & Min X & Fit Min & Fit Max & Frac Below ");
+    normPar = 0;
+    slopePar = 1;
   }
 
-
-  TFile * f = new TFile(file);  
+  
+  TFile * f = file ? new TFile(file) : gFile;  
   TH1 * h = 0;
   if(listname){
     TList * l = (TList*) gDirectory->Get(listname);
@@ -132,45 +174,59 @@ void FitParticle(const char * file, const char * histo, const char * partName,
   func->Draw("same");
   cout << "4" << endl;
 
-
   // Print results nicely
-  AliLatexTable table(9,"c|ccccccc");
-  table.InsertCustomRow("Part & Integral & T Slope & n & $\\Chi^2$/NDF & Min X & Frac Above & \\langle p_{t} \\rangle  & \\langle p_{t}^{2} \\rangle");
   // populate table
-  // Float_t yield  = func->Integral(0,100);
-  // Float_t yieldE = func->IntegralError(0,100);
+
   Double_t yield=0,yieldE=0;
   cout << "Y" << endl;
   AliBWTools::GetYield(h,func,yield,yieldE);  
   cout << "YE" << endl;
+  TLatex * l = new TLatex(2,(h->GetMaximum()+h->GetMinimum())/2,Form("%3.3f #pm %3.3f (%s)", yield,yieldE, func->GetName()));
+  l->Draw();
+
+  if(wait)  c->WaitPrimitive();
+  if(fileOut) {
+    TFile * fout = new TFile(fileOut, "update");    
+    c->SetName(Form("c_%s_%s_fit_%2.2f_%2.2f",histo, func->GetName(), min, max));
+    c->Write();
+    fout->Purge(); // remove old canvases
+    fout->Close();
+  }
+
 //   Float_t yield  = func->Integral(0.45,1.05);
 //   Float_t yieldE = func->IntegralError(0.45,1.05);
-  Double_t tslope   = func->GetParameter(2);
-  Double_t tslopeE = func->GetParError(2);     
-
-  table.SetNextCol(partName);
-  table.SetNextCol(yield,yieldE,-4);
-  table.SetNextCol(tslope,tslopeE,-4);
-  table.SetNextCol(func->GetParameter(1),func->GetParError(1),-4);
-  table.SetNextCol(Form("%2.2f/%d",func->GetChisquare(),func->GetNDF()));
-  cout << "5" << endl;
+  cout << "Slope Par: "<< slopePar << endl;
+
+  Double_t tslope   = func->GetParameter(slopePar);
+  Double_t tslopeE = func->GetParError(slopePar);      
+
+  table->SetNextCol(Form("%s (%s)", partName, h->GetName()));
+  table->SetNextCol(func->GetName());
+  table->SetNextCol(yield,yieldE,-4);
+  table->SetNextCol(tslope,tslopeE,-4);
+  if(fitFunc == kFitBlastWave) {
+    table->SetNextCol(func->GetParameter(1),func->GetParError(1),-4);
+    table->SetNextCol(func->GetParameter(3),func->GetParError(3),-4);
+  }
+  //  table->SetNextCol(func->GetParameter(1),func->GetParError(1),-4);
+  table->SetNextCol(Form("%2.2f/%d",func->GetChisquare(),func->GetNDF()));
   Float_t lowestPoint = TMath::Max(AliBWTools::GetLowestNotEmptyBinEdge(h),min);
-  cout << "6" << endl;
-  //Float_t yieldAbove = 0;
-  Float_t yieldAbove  = func->Integral(lowestPoint,100);
-  table.SetNextCol(lowestPoint,-3);
-  table.SetNextCol(yieldAbove/yield,-3);
-  Float_t mean=0, meane=0;
-  Float_t mean2=0, mean2e=0;
-  cout << "6" << endl;
-  AliBWTools::GetMean      (func, mean,  meane , 0.,100., normPar);
-  AliBWTools::GetMeanSquare(func, mean2, mean2e, 0.,100., normPar);
-  cout << "8" << endl;
-  table.SetNextCol(mean,  meane ,-4);
-  table.SetNextCol(mean2, mean2e,-4);
+  //  Float_t yieldAbove  = func->Integral(lowestPoint,100);
+  Float_t yieldBelow  = func->Integral(0,lowestPoint);
+  table->SetNextCol(lowestPoint,-3);
+  table->SetNextCol(min,-2);
+  table->SetNextCol(max,-2);
+  table->SetNextCol(yieldBelow/yield,-3);
+
+  // Float_t mean=0, meane=0;
+  // Float_t mean2=0, mean2e=0;
+  // AliBWTools::GetMean      (func, mean,  meane , 0.,100., normPar);
+  // AliBWTools::GetMeanSquare(func, mean2, mean2e, 0.,100., normPar);
+  // table->SetNextCol(mean,  meane ,-4);
+  // table->SetNextCol(mean2, mean2e,-4);
   //                    fMean2->IntegralError(0,100)/func->Integral(0,100),-7);
-  table.InsertRow();
-  table.PrintTable("ASCII");
+  table->InsertRow();
+  table->PrintTable("ASCII");