Code to combine and fit ID spectra (pi/K/p).
authormfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 8 Jul 2010 16:08:14 +0000 (16:08 +0000)
committermfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 8 Jul 2010 16:08:14 +0000 (16:08 +0000)
Warning: introduces a dependece from libMinuit on LibPWGspectra

12 files changed:
PWG2/PWG2spectraLinkDef.h
PWG2/SPECTRA/Fit/ALICEWorkInProgress.C [new file with mode: 0644]
PWG2/SPECTRA/Fit/AliBWFunc.cxx [new file with mode: 0644]
PWG2/SPECTRA/Fit/AliBWFunc.h [new file with mode: 0644]
PWG2/SPECTRA/Fit/AliBWTools.cxx [new file with mode: 0644]
PWG2/SPECTRA/Fit/AliBWTools.h [new file with mode: 0644]
PWG2/SPECTRA/Fit/AliLatexTable.cxx [new file with mode: 0644]
PWG2/SPECTRA/Fit/AliLatexTable.h [new file with mode: 0644]
PWG2/SPECTRA/Fit/CombineSpectra.C [new file with mode: 0644]
PWG2/SPECTRA/Fit/GetE735Ratios.C [new file with mode: 0644]
PWG2/SPECTRA/Fit/StarPPSpectra.C [new file with mode: 0644]
PWG2/libPWG2spectra.pkg

index 78bef61336f324a2947b062ce5f84c9cd936db50..3567be8d727c44dd5c0f66f1703f5b487dce4c4f 100644 (file)
@@ -26,4 +26,7 @@
 #pragma link C++ class AliAnalysisCentralCutEvtESD+;
 #pragma link C++ class AliAnalysisCentralExtrapolate+;
 #pragma link C++ class AliAnalysisTaskCentral+;
+#pragma link C++ class AliBWTools+;
+#pragma link C++ class AliBWFunc+;
+#pragma link C++ class AliLatexTable+;
 #endif
diff --git a/PWG2/SPECTRA/Fit/ALICEWorkInProgress.C b/PWG2/SPECTRA/Fit/ALICEWorkInProgress.C
new file mode 100644 (file)
index 0000000..fa8b880
--- /dev/null
@@ -0,0 +1,32 @@
+void ALICEWorkInProgress(TCanvas *c,TString today="11/05/2010", TString label = "ALICE performance"){
+
+  TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.72,0.75,0.82,0.89);
+  myPadLogo->SetFillColor(0); 
+  myPadLogo->SetBorderMode(0);
+  myPadLogo->SetBorderSize(2);
+  myPadLogo->SetFrameBorderMode(0);
+  myPadLogo->SetLeftMargin(0.0);
+  myPadLogo->SetTopMargin(0.0);
+  myPadLogo->SetBottomMargin(0.0);
+  myPadLogo->SetRightMargin(0.0);
+  myPadLogo->SetFillStyle(0);
+  myPadLogo->Draw();
+  myPadLogo->cd();
+  TASImage *myAliceLogo = new TASImage("~/WORK/ALICE/ANALYSIS/macros/alice_logo.png");
+  myAliceLogo->Draw();
+  c->cd();
+  TPaveText* t1=new TPaveText(0.65,0.7,0.89,0.75,"NDC");
+  t1->SetFillStyle(0);
+  t1->SetBorderSize(0);
+  t1->AddText(0.,0.,label);
+  t1->SetTextColor(kRed);
+  t1->SetTextFont(42);
+  t1->Draw();
+  TPaveText* t2=new TPaveText(0.65,0.65,0.89,0.7,"NDC");
+  t2->SetFillStyle(0);
+  t2->SetBorderSize(0);
+  t2->SetTextColor(kRed);
+  t2->SetTextFont(52);
+  t2->AddText(0.,0.,today.Data());
+  t2->Draw();
+}
diff --git a/PWG2/SPECTRA/Fit/AliBWFunc.cxx b/PWG2/SPECTRA/Fit/AliBWFunc.cxx
new file mode 100644 (file)
index 0000000..8892cbf
--- /dev/null
@@ -0,0 +1,749 @@
+//----------------------------------------------------------------------
+//                              AliBWFunc
+//
+// This class implements several function useful to fit pt spectra,
+// including but not limited to blast wave models.
+//
+// It can return the same functional for as a function of different
+// variables: dNdpt vs pt, 1/pt dNdpt vs pt, 1/mt dNdmt vs mt. 
+//
+// Before getting the function you need, you have to chose the
+// variable you want to use calling AliBWFunc::SetVarType with one of
+// the elements of the VarType_t enum.
+//
+// Warning: not all variables are iplemented for a 
+//
+// Author: M. Floris, CERN
+//----------------------------------------------------------------------
+
+#include "AliBWFunc.h"
+#include "TMath.h"
+#include "TF1.h"
+#include "TF3.h"
+#include "TH1.h"
+#include "TSpline.h"
+
+ClassImp(AliBWFunc)
+
+AliBWFunc::AliBWFunc () : fLastFunc(0),  fLineWidth(1), fVarType(kdNdpt) {
+
+  // ctor
+  fLineWidth = 1;
+}
+AliBWFunc::~AliBWFunc(){
+  
+  // dtor
+  if (fLastFunc) delete fLastFunc;
+
+}
+
+
+TF1 * AliBWFunc::GetHistoFunc(TH1 * h, const char * name) {
+
+  // Regardless of the variable type, this returns a function made
+  // from the histo * a multiplicative normalization.
+
+  fLastFunc = new TF1 (name, StaticHistoFunc, 0.0, 10, 2);
+  fLastFunc->SetParameter(0,1);
+  fLastFunc->FixParameter(1,Double_t(Long64_t(h)));
+  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){
+
+  // Boltzmann-Gibbs blast wave
+
+  switch (fVarType) {
+  case kdNdpt:
+    return GetBGBWdNdptTimesPt(mass,beta,T,norm,name);
+    break;
+  case kOneOverPtdNdpt:
+    return GetBGBWdNdpt(mass,beta,T,norm,name);
+    break;
+  case kOneOverMtdNdmt:
+    AliFatal("Not implemented");
+    break;
+  default:
+    AliFatal("Not implemented");
+  }
+  
+  return 0;
+
+}
+  
+
+TF1 * AliBWFunc::GetBoltzmann(Double_t mass, Double_t T, Double_t norm, const char * name){
+  // Boltzmann
+  switch (fVarType) {
+  case kdNdpt:
+    return GetBoltzmanndNdptTimesPt(mass, T, norm, name);
+  case kOneOverPtdNdpt:
+    AliFatal("Not implemented");
+    break;
+  case kOneOverMtdNdmt:
+    AliFatal("Not implemented");
+    break;
+  default:
+    AliFatal("Not implemented");
+  }
+  
+  return 0;
+
+}
+
+
+TF1 * AliBWFunc::GetTsallisBW(Double_t mass, Double_t beta, Double_t T, Double_t q,
+                             Double_t norm, Double_t ymax, const char * name){
+
+  // Tsallis blast wave
+  switch (fVarType) {
+  case kdNdpt:
+    return GetTsallisBWdNdptTimesPt(mass,beta,T,q,norm,ymax,name);
+    break;
+  case kOneOverPtdNdpt:
+    return GetTsallisBWdNdpt(mass,beta,T,q,norm,ymax,name);
+    break;
+  case kOneOverMtdNdmt:
+    AliFatal("Not implemented");
+    break;
+  default:
+    AliFatal("Not implemented");
+  }
+  
+  return 0;
+  
+}
+
+
+TF1 * AliBWFunc::GetMTExp(Double_t mass, Double_t T, Double_t norm, const char * name){
+
+  // Simple exponential in 1/mt*MT
+  switch (fVarType) {
+  case kdNdpt:
+    return GetMTExpdNdptTimesPt(mass,T,norm,name);
+    break;
+  case kOneOverPtdNdpt:
+    return GetMTExpdNdpt(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
+  switch (fVarType) {
+  case kdNdpt:
+    return GetPTExpdNdptTimesPt(T,norm,name);
+    break;
+  case kOneOverPtdNdpt:
+    AliFatal("Not implemented");
+    break;
+  case kOneOverMtdNdmt:
+    AliFatal("Not implemented");
+    break;
+  default:
+    AliFatal("Not implemented");
+  }
+  
+  return 0;
+
+
+}
+
+
+TF1 * AliBWFunc::GetLevi(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
+  // Levi function (aka Tsallis)
+  switch (fVarType) {
+  case kdNdpt:
+    return GetLevidNdptTimesPt(mass,T,n,norm,name);
+    break;
+  case kOneOverPtdNdpt:
+    return GetLevidNdpt(mass,T,n,norm,name);
+    break;
+  case kOneOverMtdNdmt:
+    return GetLevidNdmt(mass,T,n,norm,name);
+    break;
+  default:
+    AliFatal("Not implemented");
+  }
+  
+  return 0;
+
+
+}
+
+TF1 * AliBWFunc::GetPowerLaw(Double_t pt0, Double_t n, Double_t norm, const char * name){
+  // power law Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
+  // This is sometimes also called Hagedorn or modified Hagedorn
+
+  switch (fVarType) {
+  case kdNdpt:
+    return GetPowerLawdNdptTimesPt(pt0,n,norm,name);
+    break;
+  case kOneOverPtdNdpt:
+    return GetPowerLawdNdpt(pt0,n,norm,name);
+    break;
+  case kOneOverMtdNdmt:
+    AliFatal("Not Implemented");
+    //    return GetUA1dNdmt(mass,T,n,norm,name);
+    break;
+  default:
+    AliFatal("Not implemented");
+  }
+  
+  return 0;
+
+
+}
+
+TF1 * AliBWFunc::GetUA1(Double_t mass, Double_t p0star, Double_t pt0, Double_t n, Double_t T, Double_t norm, const char * name) {
+  // UA1 parametrization Nuclear Physics B, Vol. 335, No. 2. (7 May 1990), pp. 261-287.
+
+  switch (fVarType) {
+  case kdNdpt:
+
+    fLastFunc = new TF1 (name, StaticUA1Func, 0.0, 10, 6);
+    fLastFunc->FixParameter(0,mass);
+    fLastFunc->SetParameter(1,p0star);
+    fLastFunc->SetParameter(2,pt0);
+    fLastFunc->SetParameter(3,n);
+    fLastFunc->SetParameter(4,T);
+    fLastFunc->SetParameter(5,norm);
+    fLastFunc->SetParLimits(1,0.01,1);
+    fLastFunc->SetParLimits(2,0.01,100);
+    fLastFunc->SetParLimits(3,0.01,100);
+    fLastFunc->SetParLimits(4,0.01,100);
+    fLastFunc->SetParNames("mass","p0star","pt0","n","T","norm");
+    fLastFunc->SetNpx(5000);
+    fLastFunc->SetLineWidth(fLineWidth);
+    return fLastFunc;
+
+    break;
+  case kOneOverPtdNdpt:
+    AliFatal("Not Implemented");
+    break;
+  case kOneOverMtdNdmt:
+    AliFatal("Not Implemented");
+    //    return GetUA1dNdmt(mass,T,n,norm,name);
+    break;
+  default:
+    AliFatal("Not implemented");
+  }
+
+  return 0;
+}
+
+
+
+
+// ________________________________________________________________________
+
+// Backend (private functions and support functions for numerical integration)
+
+Double_t AliBWFunc::StaticHistoFunc(double * x, double* p){
+
+  // provides a function interpolating a histo with a spline; 
+  // using double to store a pointer... This is a bad hack. To be replaced
+
+  double norm = p[0];
+  
+  TH1 * h     = (TH1*) 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;
+  TSpline3 * spl = 0;
+  if (h!=oldptr) {
+    spl  = new TSpline3(h);
+  }
+  double value = spl->Eval(x[0]);
+  delete spl;
+
+  return value*norm;
+  
+}
+
+Double_t AliBWFunc::StaticUA1Func(double * x, double* p) {
+  
+
+  // "mass","p0star","pt0","n","T","norm"
+  Double_t mass   = p[0];
+  Double_t p0star = p[1];
+  Double_t pt0    = p[2];
+  Double_t n      = p[3];
+  Double_t T      = p[4];
+  Double_t norm   = p[5];
+  
+  Double_t xx = x[0];
+
+  static AliBWFunc * self = new AliBWFunc;
+  static TF1 * fPLaw   = self->GetPowerLawdNdptTimesPt(pt0, n, norm, "fLocalPLawUA1");
+  static TF1 * fPMTExp = self->GetMTExpdNdptTimesPt   (mass, T, norm, "fLocalMTexpUA1");
+
+  fPLaw->SetParameters(norm,pt0,n);
+  fPMTExp->SetParameters(1,T);
+  
+
+  Double_t normMT =fPMTExp->Eval(p0star) >0 ? fPLaw->Eval(p0star) / fPMTExp->Eval(p0star) *  fPMTExp->GetParameter(0) : 1;
+  fPMTExp->SetParameter(0,normMT);
+  
+  
+  if (TMath::Abs(fPMTExp->Eval(p0star) - fPLaw->Eval(p0star)) > 0.0001 ) {
+    Printf("AliBWFunc::StaticUA1Func - Wrong norm") ; 
+    Printf(" p0* %f  NMT: %f  N: %f  PL: %f  MT: %f", p0star, normMT, norm, fPLaw->Eval(p0star), fPMTExp->Eval(p0star));
+  }
+
+  if (xx > p0star)  return fPLaw->Eval(xx);
+  return fPMTExp->Eval(xx);    
+  
+  
+
+}
+
+
+Double_t AliBWFunc::IntegrandBG(double * x, double* p){
+  // integrand for boltzman-gibbs blast wave
+
+  double x0 = x[0]; 
+  
+  double mass = p[0];
+  double pT   = p[1];
+  double beta = p[2];
+  double T    = p[3];
+  
+  double mT      = TMath::Sqrt(mass*mass+pT*pT);
+
+  double rho0   = TMath::ATanH(beta*x0);  
+  double arg0_0 = pT*TMath::SinH(rho0)/T;
+  double arg0_1 = mT*TMath::CosH(rho0)/T;
+  double f0 = x0*mT*TMath::BesselI0(arg0_0)*TMath::BesselK1(arg0_1);
+
+  return f0;
+}
+
+
+
+Double_t AliBWFunc::StaticBGdNdPt(double * x, double* p) {
+
+  // implementation of BGBW (1/pt dNdpt)
+
+  double pT = x[0];;
+  
+
+  double mass = p[0];
+  double beta = p[1];
+  double T    = p[2];
+
+  static TF1 * fIntBG = 0;
+  if(!fIntBG)
+    fIntBG = new TF1 ("fIntBG", IntegrandBG, 0, 1, 4);
+
+  fIntBG->SetParameters(mass, pT, beta, T);
+  double result = fIntBG->Integral(0,1);
+  return result*p[3];//*1e30;;
+
+}
+
+Double_t AliBWFunc::StaticBGdNdPtTimesPt(double * x, double* p) {
+  // BGBW dNdpt implementation
+  return x[0]*StaticBGdNdPt(x,p);
+}
+
+
+TF1 * AliBWFunc::GetBGBWdNdpt(Double_t mass, Double_t beta, Double_t T,
+                           Double_t norm, const char * name){
+  
+  // BGBW 1/pt dNdpt
+
+  fLastFunc = new TF1 (name, StaticBGdNdPt, 0.0, 10, 4);
+  fLastFunc->SetParameters(mass,beta,T,norm);    
+  fLastFunc->SetParNames("mass", "#beta", "T", "norm");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+  
+}
+
+
+//_____________________________________________________________________
+// Tsallis
+
+Double_t AliBWFunc::IntegrandTsallis(double * x, double* p){
+
+  // integrand for numerical integration (tsallis)
+
+  Double_t r   = x[0]; 
+  Double_t phi = x[1];
+  Double_t y   = x[2];
+
+  Double_t mass = p[0];
+  Double_t pt   = p[1];
+  Double_t beta = p[2];
+  Double_t T    = p[3];
+  Double_t q    = p[4];
+  
+  Double_t mt      = TMath::Sqrt(mass*mass+pt*pt);
+
+  Double_t rho    = TMath::ATanH(beta*r); // TODO: implement different velocity profiles  
+
+  Double_t res = mt*
+    r*TMath::CosH(y) *TMath::Power( (
+                                    1+(q-1)/T * (
+                                                 mt*TMath::CosH(y)*TMath::CosH(rho) -
+                                                 pt*TMath::SinH(rho)*TMath::Cos(phi)
+                                                 )
+                                    ),
+                                      -1/(q-1)
+                                   );                  
+
+
+  return res;
+}
+
+
+
+Double_t AliBWFunc::StaticTsallisdNdPt(double * x, double* p) {
+
+  // tsallis BW implementation 1/pt dNdpt
+
+  double pT = x[0];;
+  
+
+  double mass = p[0];
+  double beta = p[1];
+  double T    = p[2];
+  double q    = p[3];
+
+  Double_t ymax = p[5];
+
+
+  static TF3 * fInt = 0;
+  if(!fInt){
+    fInt = new TF3 ("fIntTsa", IntegrandTsallis, 0, 1, -TMath::Pi(), TMath::Pi(), -ymax, ymax, 5);
+//     fInt->SetNpx(10000);
+//     fInt->SetNpy(10000);
+//     fInt->SetNpz(10000);
+  }
+  
+  fInt->SetParameters(mass, pT, beta, T, q);
+  double result = fInt->Integral(0,1, -TMath::Pi(), TMath::Pi(), -ymax, ymax);
+  //  double result = fInt->Integral(0,1, -2, 2, -ymax, ymax);
+  
+  return result*p[4];//*1e30;;
+
+}
+
+Double_t AliBWFunc::StaticTsallisdNdPtTimesPt(double * x, double* p) {
+
+  // tsallis BW , implementatio of dNdpt
+  return x[0]*StaticTsallisdNdPt(x,p);
+
+}
+
+TF1 * AliBWFunc::GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t T, Double_t q,
+                                  Double_t norm, Double_t ymax,const char * name){
+  
+
+  // tsallis BW, 1/pt dNdpt
+
+  fLastFunc = new TF1 (name, StaticTsallisdNdPt, 0.0, 10, 6);
+  fLastFunc->SetParameters(mass,beta,T,q,norm,ymax);
+  fLastFunc->SetParLimits(1,0.0,0.99);
+  fLastFunc->SetParLimits(2,0.01,0.99);
+  fLastFunc->SetParLimits(3,1.0001,1.9);
+  fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+  
+}
+
+// Times Pt funcs
+// Boltzmann-Gibbs Blast Wave
+TF1 * AliBWFunc::GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T,
+                                    Double_t norm, const char * name){
+
+  // BGBW, dNdpt
+
+  fLastFunc = new TF1 (name, StaticBGdNdPtTimesPt, 0.0, 10, 4);
+  fLastFunc->SetParameters(mass,beta,T,norm);    
+  fLastFunc->SetParNames("mass", "#beta", "T", "norm");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+
+
+TF1 * AliBWFunc::GetTsallisBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T, Double_t q,
+                                         Double_t norm, Double_t ymax, const char * name){
+
+// Tsallis blast wave, dNdpt
+
+  fLastFunc = new TF1 (name, StaticTsallisdNdPtTimesPt, 0.0, 10, 6);
+  fLastFunc->SetParameters(mass,beta,T,q,norm,ymax);    
+  fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+}
+
+
+
+TF1 * AliBWFunc::GetMTExpdNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name){
+
+  // Simple exponential in 1/mt*MT, as a function of dNdpt
+  char formula[500];
+  sprintf(formula,"[0]*x*exp(-sqrt(x**2+%f**2)/[1])", mass);
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, T);
+  fLastFunc->SetParLimits(1, 0.01, 10);
+  fLastFunc->SetParNames("norm", "T");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+
+TF1 * AliBWFunc::GetPTExpdNdptTimesPt(Double_t T, Double_t norm, const char * name){
+
+  // Simple exponential in 1/pt*dNdpT, as a function of dNdpt
+  char formula[500];
+  sprintf(formula,"[0]*x*exp(-x/[1])");
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, T);
+  fLastFunc->SetParLimits(1, 0.01, 10);
+  fLastFunc->SetParNames("norm", "T");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+
+TF1 * AliBWFunc::GetBoltzmanndNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name){
+  // Boltzmann (exp in 1/mt*dNdmT times mt) as a function of dNdpt
+ char formula[500];
+ sprintf(formula,"[0]*x*sqrt(x**2+%f**2)*exp(-sqrt(x**2+%f**2)/[1])", mass,mass);
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, T);
+  fLastFunc->SetParLimits(1, 0.01, 10);
+  fLastFunc->SetParNames("norm", "T");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+
+// Tsallis (no BW, a la CMS)
+// TF1 * AliBWFunc::GetTsallisdNdptTimesPt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name){
+
+//   char formula[500];
+//   //  sprintf(formula,"[0]*x*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2)-%f)),(-1/([2]-1)))", mass, mass); //CMS
+//   sprintf(formula,"[0]*x*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass);  // STAR
+//   //sprintf(formula,"[0]*x*sqrt(x**2+%f**2)*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass,mass);  // STAR * mt
+//   fLastFunc=new TF1(name,formula,0,10);
+//   fLastFunc->SetParameters(norm, T, q);
+//   fLastFunc->SetParLimits(1, 0.001, 10);
+//   fLastFunc->SetParNames("norm", "T", "q");
+//   fLastFunc->SetLineWidth(fLineWidth);
+//   return fLastFunc;
+
+
+// }
+
+
+TF1 * AliBWFunc::GetLevidNdptTimesPt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
+
+  // Levi function, dNdpt
+  char formula[500];
+
+  sprintf(formula,"( x*[0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2])  )^(-[1])");
+  //  sprintf(formula,"( x*[0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (sqrt([3]*[3]+x*x))/([1]*[2])  )^(-[1])");
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, n, T,mass);
+  fLastFunc->SetParLimits(2, 0.01, 10);
+  fLastFunc->SetParNames("norm (dN/dy)", "n", "T", "mass");
+  fLastFunc->FixParameter(3,mass);
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+TF1 * AliBWFunc::GetPowerLawdNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name){
+
+  // PowerLaw function, dNdpt
+  char formula[500];
+
+  sprintf(formula,"x*[0]*( 1 + x/[1] )^(-[2])");
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, pt0, n);
+  fLastFunc->SetParLimits(1, 0.01, 10);
+  //fLastFunc->SetParLimits(2, 0.01, 50);
+  fLastFunc->SetParNames("norm", "pt0", "n");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+TF1 * AliBWFunc::GetPowerLawdNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name){
+
+  // PowerLaw function, 1/pt dNdpt
+  char formula[500];
+
+  sprintf(formula," [0]*( 1 + x/[1] )^(-[2])");
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, pt0, n);
+  //  fLastFunc->SetParLimits(2, 0.01, 10);
+  fLastFunc->SetParNames("norm", "pt0", "n");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+
+TF1 * AliBWFunc::GetLevidNdpt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
+
+  // Levi function, dNdpt
+  char formula[500];
+
+  sprintf(formula,"( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (sqrt([3]*[3]+x*x) -[3])/([1]*[2])  )^(-[1])");
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, n, T,mass);
+  fLastFunc->SetParLimits(2, 0.01, 10);
+  fLastFunc->SetParNames("norm (dN/dy)", "n", "T", "mass");
+  fLastFunc->FixParameter(3,mass);
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+TF1 * AliBWFunc::GetLevidNdmt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name){
+
+  // Levi function, dNdmt
+  char formula[500];
+
+  //  sprintf(formula,"( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + (x -[3])/([1]*[2])  )^(-[1])");
+  sprintf(formula,"( [0]*([1]-1)*([1]-2)  )/( [1]*[2]*( [1]*[2]+[3]*([1]-2) )  ) * ( 1 + x/([1]*[2])  )^(-[1])");
+  //  sprintf(formula,"[0] * ( 1 + x/([1]*[2])  )^(-[1])");
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, n, T,mass);
+  fLastFunc->SetParLimits(2, 0.01, 10);
+  fLastFunc->SetParNames("norm", "n", "T", "mass");
+  fLastFunc->FixParameter(3,mass);
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+
+
+}
+
+
+
+// Test Function
+Double_t AliBWFunc::IntegrandTest(double * x, double* p){
+
+  // test function
+
+  Double_t y = x[0];
+
+  Double_t mass = p[0];
+  Double_t pt   = p[1];
+  Double_t T    = p[2];
+
+  Double_t mt      = TMath::Sqrt(mass*mass+pt*pt);    
+  
+  return mt*TMath::CosH(y)*TMath::Exp(-mt*TMath::CosH(y)/T);
+
+}
+
+Double_t AliBWFunc::StaticTest(double * x, double* p) {
+
+  // test function
+
+  double pT = x[0];;
+  
+
+  double mass = p[0];
+  double T    = p[1];
+  Double_t ymax = p[3];
+
+
+  static TF3 * fIntTest = 0;
+  if(!fIntTest){
+    fIntTest = new TF3 ("fIntTest", IntegrandTest, 0, 1, -TMath::Pi(), TMath::Pi(), -ymax, ymax, 5);
+    //    fInt->SetNpx(10000);
+  }
+  
+  fIntTest->SetParameters(mass, pT, T);
+  double result = fIntTest->Integral(-ymax, ymax);
+  
+  return result*p[2];//*1e30;;
+
+}
+
+TF1 * AliBWFunc::GetTestFunc(Double_t mass, Double_t T, Double_t norm, Double_t ymax, const char * name){
+  
+  // test function
+  
+  fLastFunc = new TF1 (name, StaticTest, 0.0, 10, 4);
+  fLastFunc->SetParameters(mass,T,norm,ymax);    
+  fLastFunc->SetParNames("mass", "#beta", "T", "q", "norm", "ymax");
+  fLastFunc->SetLineWidth(fLineWidth);
+  return fLastFunc;
+  
+}
+
+
+//___________________________________________________________
+
+
+TF1 * AliBWFunc::GetMTExpdNdpt(Double_t mass, Double_t T, Double_t norm, const char * name){
+  // Simple exp in 1/mt dNdmt, as a function of dNdpt
+  // mt scaling
+  char formula[500];
+  sprintf(formula,"[0]*exp(-sqrt(x**2+%f**2)/[1])", mass);
+  fLastFunc=new TF1(name,formula,0,10);
+  fLastFunc->SetParameters(norm, T);
+  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 T, Double_t q, Double_t norm, const char * name){
+  
+//   char formula[500];
+//   sprintf(formula,"[0]*sqrt(x**2+%f**2)*pow((1+(([2]-1)/[1])*(sqrt(x**2+%f**2))),(-1/([2]-1)))", mass,mass); 
+//   fLastFunc=new TF1(name,formula,0,10);
+//   fLastFunc->SetParameters(norm, T, q);
+//   fLastFunc->SetParLimits(1, 0.01, 10);
+//   fLastFunc->SetParNames("norm", "T", "q");
+//   fLastFunc->SetLineWidth(fLineWidth);
+//   return fLastFunc;
+// }
diff --git a/PWG2/SPECTRA/Fit/AliBWFunc.h b/PWG2/SPECTRA/Fit/AliBWFunc.h
new file mode 100644 (file)
index 0000000..cd6b26c
--- /dev/null
@@ -0,0 +1,184 @@
+//----------------------------------------------------------------------
+//                              AliBWFunc
+//
+// This class implements several function useful to fit pt spectra,
+// including but not limited to blast wave models.
+//
+// It can return the same functional for as a function of different
+// variables: dNdpt vs pt, 1/pt dNdpt vs pt, 1/mt dNdmt vs mt. 
+//
+// Before getting the function you need, you have to chose the
+// variable you want to use calling AliBWFunc::SetVarType with one of
+// the elements of the VarType_t enum.
+//
+// Author: M. Floris, CERN
+//----------------------------------------------------------------------
+
+#ifndef ALIBWFUNC_H
+#define ALIBWFUNC_H
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "TObject.h"
+#include "AliLog.h"
+
+class TF1;
+class TH1;
+
+#endif
+
+
+
+class AliBWFunc : public TObject {
+
+
+public:
+  // define the variables used for the function
+  typedef enum {kdNdpt,kOneOverPtdNdpt,kOneOverMtdNdmt} VarType_t;
+
+  AliBWFunc();
+  ~AliBWFunc();
+
+  // Boltzmann-Gibbs blast wave
+  TF1 * GetBGBW(Double_t mass, Double_t beta, Double_t T,
+               Double_t norm, const char * name = "fBGBW");
+  
+  // Boltzmann
+  TF1 * GetBoltzmann(Double_t mass, Double_t T, Double_t norm, const char * name ="fBoltzmann");
+
+  // Tsallis blast wave
+  TF1 * GetTsallisBW(Double_t mass, Double_t beta, Double_t T, Double_t q,
+                    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");
+
+  // Simple exponential in 1/pt*dNdpt
+  TF1 * GetPTExp(Double_t T, Double_t norm, const char * name ="fExp");
+
+  // Tsallis (no BW, a la CMS)
+  TF1 * GetTsallis(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name="fTsallis") 
+  {return GetLevi (mass,T,1/(q-1),norm,name);}
+  
+  // Levi function (aka Tsallis)
+  TF1 * GetLevi(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name="fLevi");
+
+  // UA1 function 
+  TF1 * GetUA1(Double_t mass, Double_t p0star, Double_t pt0, Double_t n, Double_t T, Double_t norm, const char * name="fUA1");
+
+  // Function derived from a histo
+  TF1 * GetHistoFunc(TH1 * h, const char * name = "fHisto");
+
+  // Power law
+  TF1 * GetPowerLaw(Double_t pt0, Double_t n, Double_t norm, const char * name="fPowerLaw");
+
+
+  void SetVarType(VarType_t tp) {fVarType=tp;}
+
+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");
+
+  // Tsallis blast wave
+  TF1 * GetTsallisBWdNdpt(Double_t mass, Double_t beta, Double_t T, Double_t q,
+                         Double_t norm, Double_t ymax = 0.5, const char * name = "fTsallisBW");
+
+  // Simple exponential in 1/mt*MT
+  TF1 * GetMTExpdNdpt(Double_t mass, Double_t T, Double_t norm, const char * name ="fExp");
+
+  // Tsallis (no BW, a la CMS)
+  TF1 * GetTsallisdNdpt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name="fTsallis");
+
+  // Levi function
+  TF1 * GetLevidNdpt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name="fLevi");
+
+  // Power Law function
+  TF1 * GetPowerLawdNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name="fLevi");
+
+  // UA1 function
+  TF1 * GetUA1dNdpt(Double_t pt0, Double_t n, Double_t norm, const char * name="fLevi");
+
+  // TimesPt means dNdpt
+
+  // Boltzmann-Gibbs Blast Wave
+  TF1 * GetBGBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T,
+                           Double_t norm, const char * name = "fBGBWTimesPt");
+
+  // Tsallis blast wave
+  TF1 * GetTsallisBWdNdptTimesPt(Double_t mass, Double_t beta, Double_t T, Double_t q,
+                                Double_t norm, Double_t ymax = 0.5, const char * name = "fTsallisBWTimesPt");
+  // Levi function
+  TF1 * GetLevidNdptTimesPt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name="fLevi");
+
+  // UA1 function
+  TF1 * GetUA1dNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name="fLevi");
+
+  // PowerLaw function
+  TF1 * GetPowerLawdNdptTimesPt(Double_t pt0, Double_t n, Double_t norm, const char * name="fLevi");
+
+  // Simple exponential in 1/mt*dNdmT
+  TF1 * GetMTExpdNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name ="fMtExpTimesPt");
+
+  // Simple exponential in 1/mp*dNdpT
+  TF1 * GetPTExpdNdptTimesPt(Double_t T, Double_t norm, const char * name ="fPtExpTimesPt");
+
+  // Boltzmann (exp in 1/mt*dNdmT times mt)
+  TF1 * GetBoltzmanndNdptTimesPt(Double_t mass, Double_t T, Double_t norm, const char * name ="fBoltzmannTimesPt");
+
+  // Tsallis (no BW, a la CMS)
+  TF1 * GetTsallisdNdptTimesPt(Double_t mass, Double_t T, Double_t q, Double_t norm, const char * name="fTsallisTimesPt");
+
+  // 1/mt dNdmt
+  
+  // Levi
+  TF1 * GetLevidNdmt(Double_t mass, Double_t T, Double_t n, Double_t norm, const char * name = "fLeviMt");
+
+
+  // gereral setters
+  void SetLineWidth(Width_t width) { fLineWidth = width;}
+
+  TF1 * GetTestFunc(Double_t mass, Double_t T, Double_t norm, Double_t ymax, const char * name ="fTest") ;
+
+  // static functions for TF1
+  // Boltzmann-Gibbs Blast Wave
+  static Double_t StaticBGdNdPt(double * x, double* p);
+  static Double_t StaticBGdNdPtTimesPt(double * x, double* p);
+  // Tsallis blast wave
+  static Double_t StaticTsallisdNdPt(double * x, double* p);
+  static Double_t StaticTsallisdNdPtTimesPt(double * x, double* p);
+  // Helper funcs for numeric integration
+  static Double_t IntegrandBG(double * x, double* p);
+  static Double_t IntegrandTsallis(double * x, double* p);
+
+  // Test func
+  static Double_t StaticTest(double * x, double* p);
+  static Double_t IntegrandTest(double * x, double* p);
+
+  // histo func
+  static Double_t StaticHistoFunc(double * x, double* p);
+
+  // UA1 parametrization
+  static Double_t StaticUA1Func(double * x, double* p);
+  
+
+private:
+
+
+  TF1 * fLastFunc;
+  Width_t fLineWidth;
+  VarType_t fVarType;
+  
+  AliBWFunc(const AliBWFunc&);            // not implemented
+  AliBWFunc& operator=(const AliBWFunc&); // not implemented
+
+
+  ClassDef(AliBWFunc, 1)
+
+
+};
+
+#endif
diff --git a/PWG2/SPECTRA/Fit/AliBWTools.cxx b/PWG2/SPECTRA/Fit/AliBWTools.cxx
new file mode 100644 (file)
index 0000000..b77e744
--- /dev/null
@@ -0,0 +1,780 @@
+// ----------------------------------------------------------------------
+//                     AliBWTools
+// 
+// This class provides some tools which can be useful in the analsis
+// of spectra, to fit or transform histograms. See the comments of the
+// individual methods for details
+//
+// Author: M. Floris (CERN)
+// ----------------------------------------------------------------------
+
+#include "AliBWTools.h"
+#include "TH1D.h"
+#include "TF1.h"
+#include "TH1.h"
+#include "TMath.h"
+#include "TGraphErrors.h"
+#include "TVirtualFitter.h"
+#include "TMinuit.h"
+#include "AliLog.h"
+#include <iostream>
+
+using namespace std;
+
+TF1 * AliBWTools::fFuncForNormalized = 0;
+
+ClassImp(AliBWTools)
+
+AliBWTools::AliBWTools() {
+  // ctor
+}
+
+AliBWTools::~AliBWTools(){
+  // dtor
+}
+
+TH1 * AliBWTools::GetdNdmtFromdNdpt(TH1 * hpt, Double_t mass) {
+  // convert the x axis from pt to mt. Assumes you have 1/pt dNdpt in the histo you start with
+
+  Int_t nbins = hpt->GetNbinsX();
+  Float_t * xbins = new Float_t[nbins+1];
+  for(Int_t ibins = 0; ibins <= nbins; ibins++){
+    xbins[ibins] = TMath::Sqrt(hpt->GetBinLowEdge(ibins+1)*hpt->GetBinLowEdge(ibins+1) +
+                              mass *mass) - mass;
+//     xbins[ibins] = TMath::Sqrt(hpt->GetBinLowEdge(ibins+1)*hpt->GetBinLowEdge(ibins+1) +
+//                            mass *mass);
+    //    cout << ibins << " "<< xbins[ibins]  << endl;
+
+  }
+
+  TH1D * hmt = new TH1D(TString(hpt->GetName())+"_mt",
+                     TString(hpt->GetName())+"_mt",
+                     nbins, xbins);
+  for(Int_t ibins = 1; ibins <= nbins; ibins++){
+    hmt->SetBinContent(ibins, hpt->GetBinContent(ibins));
+    hmt->SetBinError(ibins,   hpt->GetBinError(ibins));
+
+  }
+
+  hmt->SetXTitle("m_{t} - m_{0} (GeV/c^{2})");
+  hmt->SetYTitle("1/m_{t} dN/dm_{t} (a.u.)");
+  
+  return hmt;
+
+}
+
+TH1 * AliBWTools::GetdNdPtFromOneOverPt(TH1 * h1Pt) {
+
+  // convert an histo from 1/pt dNdpt to dNdpt, using the pt at the center of the bin
+
+
+  TH1 * hPt = (TH1 *) h1Pt->Clone((TString(h1Pt->GetName()) + "_inv").Data());
+  hPt->Reset();
+
+  Int_t nbinx = hPt->GetNbinsX();
+
+  for(Int_t ibinx = 1; ibinx <= nbinx; ibinx++){
+
+    Double_t cont = h1Pt->GetBinContent(ibinx);
+    Double_t err  = h1Pt->GetBinError(ibinx);
+    
+    Double_t pt   = h1Pt->GetBinCenter(ibinx);
+    
+    if(pt != 0) {
+      cont *= pt;
+      err  *= pt;
+    } else {
+      cont = 0;
+      err  = 0;
+    }
+
+    hPt->SetBinContent(ibinx, cont);
+    hPt->SetBinError(ibinx, err);
+    
+  }
+
+  hPt->SetXTitle("p_{t} (GeV)");
+  hPt->SetYTitle("dN/dp_{t} (GeV^{-2})");
+
+  return hPt;    
+
+}
+
+
+
+
+TH1 * AliBWTools::GetOneOverPtdNdPt(TH1 * hPt) {
+
+  // convert an histo from dNdpt to 1/pt dNdpt, using the pt at the center of the bin
+
+  TH1 * h1Pt = (TH1 *) hPt->Clone((TString(hPt->GetName()) + "_inv").Data());
+  h1Pt->Reset();
+
+  Int_t nbinx = h1Pt->GetNbinsX();
+
+  for(Int_t ibinx = 1; ibinx <= nbinx; ibinx++){
+
+    Double_t cont = hPt->GetBinContent(ibinx);
+    Double_t err  = hPt->GetBinError(ibinx);
+    
+    Double_t pt   = hPt->GetBinCenter(ibinx);
+    
+    if(pt != 0) {
+      cont /= pt;
+      err  /= pt;
+    } else {
+      cont = 0;
+      err  = 0;
+    }
+
+    h1Pt->SetBinContent(ibinx, cont);
+    h1Pt->SetBinError(ibinx, err);
+    
+  }
+
+  h1Pt->SetXTitle("p_{t} (GeV)");
+  h1Pt->SetYTitle("1/p_{t} dN/dp_{t} (GeV^{-2})");
+
+  return h1Pt;    
+
+}
+
+
+TGraphErrors * AliBWTools::GetGraphFromHisto(TH1F * h, Bool_t binWidth) {
+  // Convert a histo to a graph
+  // if binWidth is true ex is set to the bin width of the histos, otherwise it is set to zero
+  Int_t nbin = h->GetNbinsX();
+
+  TGraphErrors * g = new TGraphErrors();
+  Int_t ipoint = 0;
+  for(Int_t ibin = 1; ibin <= nbin; ibin++){
+    Double_t xerr = binWidth ? h->GetBinWidth(ibin)/2 : 0;
+    if (h->GetBinContent(ibin)!=0) {
+      g->SetPoint     (ipoint,   h->GetBinCenter(ibin), h->GetBinContent(ibin));
+      g->SetPointError(ipoint,   xerr,  h->GetBinError(ibin));
+      ipoint++;
+    }
+  }
+  
+  g->SetMarkerStyle(h->GetMarkerStyle());
+  g->SetMarkerColor(h->GetMarkerColor());
+  g->SetLineColor(h->GetLineColor());
+
+  g->SetTitle(h->GetTitle());
+  g->SetName(TString("g_")+h->GetName());
+
+  return g;
+
+}
+
+TH1F * AliBWTools::GetHistoFromGraph(TGraphErrors * g, TH1F* hTemplate) {
+
+  // convert a graph to histo with the binning of hTemplate.
+  // Warning: the template should be chosen
+  // properly: if you have two graph points in the same histo bin this
+  // won't work!
+
+  TH1F * h = (TH1F*) hTemplate->Clone(TString("h_")+g->GetName());
+  h->Reset();
+  Int_t npoint = g->GetN();
+  //  g->Print();
+  for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
+    Float_t x  = g->GetX() [ipoint];
+    Float_t y  = g->GetY() [ipoint];
+    Float_t ey = g->GetEY()[ipoint];
+    Int_t bin = h->FindBin(x);
+    //    cout << "bin: "<< bin << " -> " << x << ", "<< y <<", " << ey << endl;
+    
+    h->SetBinContent(bin,y);
+    h->SetBinError  (bin,ey);
+  }
+  h->SetMarkerStyle(g->GetMarkerStyle());
+  h->SetMarkerColor(g->GetMarkerColor());
+  h->SetLineColor  (g->GetLineColor());
+
+  return h;
+}
+
+TGraphErrors * AliBWTools::ConcatenateGraphs(TGraphErrors * g1,TGraphErrors * g2){
+
+  // concatenates two graphs
+
+  Int_t npoint1=g1->GetN();
+  Int_t npoint2=g2->GetN();
+
+  TGraphErrors * gClone = (TGraphErrors*) g1->Clone();
+
+//   for(Int_t ipoint = 0; ipoint < npoint1; ipoint++){
+//     gClone->SetPointError(ipoint,0,g1->GetEY()[ipoint]);
+
+//   }
+  for(Int_t ipoint = 0; ipoint < npoint2; ipoint++){
+    gClone->SetPoint(ipoint+npoint1,g2->GetX()[ipoint],g2->GetY()[ipoint]);
+    gClone->SetPointError(ipoint+npoint1,g2->GetEX()[ipoint],g2->GetEY()[ipoint]);
+    //    gClone->SetPointError(ipoint+npoint1,0,g2->GetEY()[ipoint]);
+  }
+
+  gClone->GetHistogram()->GetXaxis()->SetTimeDisplay(1);
+  gClone->SetTitle(TString(gClone->GetTitle())+" + "+g2->GetTitle());
+  gClone->SetName(TString(gClone->GetName())+"_"+g2->GetName());
+
+  return gClone;
+}
+
+
+TH1F * AliBWTools::CombineHistos(TH1 * h1, TH1 * h2, TH1* htemplate, Float_t renorm1){
+
+  // Combine two histos. This assumes the histos have the same binning
+  // in the overlapping region. It computes the arithmetic mean in the
+  // overlapping region and assigns as an error the relative error
+  // h2. TO BE IMPROVED
+
+  TH1F * hcomb = (TH1F*) htemplate->Clone(TString(h1->GetName())+"_"+h2->GetName());
+
+  TH1F * h1local = (TH1F*) h1->Clone();
+  h1local->Scale(renorm1);
+  
+  Int_t nBinComb = hcomb->GetNbinsX();
+  for(Int_t ibin = 1; ibin <= nBinComb; ibin++){
+    Int_t ibin1 = h1local->FindBin(hcomb->GetBinCenter(ibin));
+    Int_t ibin2 = h2->FindBin(hcomb->GetBinCenter(ibin));
+    
+      if (h1local->GetBinContent(ibin1) == 0 && h2->GetBinContent(ibin2) == 0) {
+       // None has data: go to next bin
+       hcomb->SetBinContent(ibin,0);
+       hcomb->SetBinError  (ibin,0);   
+      } else if(h1local->GetBinContent(ibin1) != 0 && h2->GetBinContent(ibin2) == 0) {
+       // take data from h1local:
+       hcomb->SetBinContent(ibin,h1local->GetBinContent(ibin1));
+       hcomb->SetBinError  (ibin,h1local->GetBinError(ibin1));
+      } else if(h1local->GetBinContent(ibin1) == 0 && h2->GetBinContent(ibin2) != 0) {
+       // take data from h2:
+       hcomb->SetBinContent(ibin,h2->GetBinContent(ibin2));
+       hcomb->SetBinError  (ibin,h2->GetBinError(ibin2));
+      }  else {
+       hcomb->SetBinContent(ibin,(h1local->GetBinContent(ibin1) +h2->GetBinContent(ibin2))/2);
+       //      hcomb->SetBinError  (ibin,h1local->GetBinError(ibin1)/h1local->GetBinContent(ibin1)*hcomb->GetBinContent(ibin));
+       hcomb->SetBinError  (ibin,h2->GetBinError(ibin2)/h2->GetBinContent(ibin2)*hcomb->GetBinContent(ibin));
+      }
+
+
+  }
+  
+
+  hcomb->SetMarkerStyle(h1local->GetMarkerStyle());
+  hcomb->SetMarkerColor(h1local->GetMarkerColor());
+  hcomb->SetLineColor  (h1local->GetLineColor());
+
+  hcomb->SetXTitle(h1local->GetXaxis()->GetTitle());
+  hcomb->SetYTitle(h1local->GetYaxis()->GetTitle());
+  delete h1local;
+  return hcomb;  
+
+}
+
+
+void AliBWTools::GetFromHistoGraphDifferentX(TH1F * h, TF1 * f, TGraphErrors ** gBarycentre, TGraphErrors ** gXlw) {
+
+  // Computes the baycentre in each bin and the xlw as defined in NIMA
+  // 355 - 541 using f. Returs 2 graphs with the same y content of h
+  // but with a different x (barycentre and xlw)
+
+  Int_t nbin = h->GetNbinsX();
+  
+  (*gBarycentre) = new TGraphErrors();
+  (*gXlw)        = new TGraphErrors();
+
+  Int_t ipoint = 0;
+  for(Int_t ibin = 1; ibin <= nbin; ibin++){
+    Float_t min = h->GetBinLowEdge(ibin);
+    Float_t max = h->GetBinLowEdge(ibin+1);
+    Double_t xerr = 0;
+    Double_t xbar = f->Mean(min,max);
+    // compute x_LW
+    Double_t temp = 1./(max-min) * f->Integral(min,max);
+    Double_t epsilon   = 0.000000001;
+    Double_t increment = 0.0000000001;
+    Double_t x_lw = min;
+
+    while ((f->Eval(x_lw)- temp) > epsilon) {
+      x_lw += increment;
+      if(x_lw > max) {
+       Printf("Cannot find x_lw");
+       break;
+      }
+    }
+      
+    if (h->GetBinContent(ibin)!=0) {
+      (*gBarycentre)->SetPoint     (ipoint,   xbar, h->GetBinContent(ibin));
+      (*gBarycentre)->SetPointError(ipoint,   xerr, h->GetBinError(ibin));
+      (*gXlw)       ->SetPoint     (ipoint,   x_lw, h->GetBinContent(ibin));
+      (*gXlw)       ->SetPointError(ipoint,   xerr, h->GetBinError(ibin));
+      ipoint++;
+    }
+  }
+  
+  (*gBarycentre)->SetMarkerStyle(h->GetMarkerStyle());
+  (*gBarycentre)->SetMarkerColor(h->GetMarkerColor());
+  (*gBarycentre)->SetLineColor(h->GetLineColor());
+
+  (*gBarycentre)->SetTitle(h->GetTitle());
+  (*gBarycentre)->SetName(TString("g_")+h->GetName());
+
+  (*gXlw)->SetMarkerStyle(h->GetMarkerStyle());
+  (*gXlw)->SetMarkerColor(h->GetMarkerColor());
+  (*gXlw)->SetLineColor(h->GetLineColor());
+  (*gXlw)->SetTitle(h->GetTitle());
+  (*gXlw)->SetName(TString("g_")+h->GetName());
+
+
+}
+
+
+Float_t AliBWTools::GetMean(TH1F * h, Float_t min, Float_t max) {
+
+  // Get the mean of histo in a range; root is not reliable in sub ranges with variable binning.
+
+  Int_t minBin = h->FindBin(min);
+  Int_t maxBin = h->FindBin(max);
+
+  Float_t mean = 0 ;
+  Float_t integral = 0;
+  for(Int_t ibin = minBin; ibin < maxBin; ibin++){
+    mean     = mean + (h->GetBinCenter(ibin) *  h->GetBinWidth(ibin)* h->GetBinContent(ibin));
+    integral = integral + h->GetBinContent(ibin) * h->GetBinWidth(ibin);
+  }
+  
+  return mean/integral;
+
+
+}
+
+void AliBWTools::GetMean(TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) {
+
+  // Get the mean of function in a range; 
+
+  return GetMoment("fMean", "x*", func, mean, error, min,max);
+
+}
+
+void AliBWTools::GetMeanSquare(TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) {
+
+  // Get the mean2 of function in a range; 
+
+  return GetMoment("fMean2", "x*x*", func, mean, error, min,max);
+
+
+}
+
+void AliBWTools::GetMoment(TString name, TString var, TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) {
+
+  // returns the integral of a function derived from func by prefixing some operation.
+  // Used as a base method for mean and mean2
+  Printf("AliBWTools::GetMoment: Error on <pt> is not correct!!! It is overestimated, fix required");
+  Int_t npar = func->GetNpar();
+
+  TF1 * f = new TF1(name, var+func->GetName());        // FIXME
+//   fFuncForNormalized = func;// TMP: computing mean pt
+//   TF1 * f = new TF1(name,GetNormalizedFunc,0,10,npar);// FIXME
+//   for(Int_t ipar = 0; ipar < npar; ipar++){ // FIXME
+//     f->SetParameter(ipar,func->GetParameter(ipar)); // FIXME
+//   }
+  
+  
+  // The parameter of the function used to compute the mean should be
+  // the same as the parent function: fixed if needed and they should
+  // also have the same errors.
+
+  //  cout << "npar :" << npar << endl;
+  
+  for(Int_t ipar = 0; ipar < npar; ipar++){
+    Double_t parmin, parmax;
+    Double_t value = func->GetParameter(ipar);
+    func->GetParLimits(ipar, parmin, parmax);
+    if ( parmin == parmax )   {
+      if ( parmin != 0 || (parmin == 1 && value == 0) ) {
+       f->FixParameter(ipar,func->GetParameter(ipar));
+       //      cout << "Fixing " << ipar << "("<<value<<","<<parmin<<","<<parmax<<")"<<endl;
+      }       
+      else {
+       f->SetParError (ipar,func->GetParError(ipar) );
+       //      cout << "Setting Err" << ipar << "("<<func->GetParError(ipar)<<")"<<endl;      
+      }
+    }
+    else {
+      f->SetParError (ipar,func->GetParError(ipar) );
+      //      cout << "Setting Err" << ipar << "("<<func->GetParError(ipar)<<")"<<endl;      
+    }
+  }
+  //  f->Print();
+//   mean  = f->Integral     (min,max)/func->Integral(min,max);
+//   error = f->IntegralError(min,max)/func->Integral(min,max);
+  mean  = f->Integral     (min,max);
+  error = f->IntegralError(min,max);
+//   cout << "Mean " << mean <<"+-"<< error<< endl;
+//   cout << "Integral Error "  << error << endl;
+  
+}
+
+Double_t AliBWTools::GetNormalizedFunc(double * x, double* p){
+
+  // Static function used to provide normalized pointer to a function
+
+  Int_t npar = fFuncForNormalized->GetNpar();
+  for(Int_t ipar = 0; ipar < npar; ipar++){ // FIXME
+    fFuncForNormalized->SetParameter(ipar,p[ipar]); // FIXME
+  }
+
+  return x[0]*fFuncForNormalized->Eval(x[0])/fFuncForNormalized->Integral(0,100);
+  
+}
+
+
+Bool_t AliBWTools::Fit (TH1 * h1, TF1* func, Float_t min, Float_t max) { 
+  
+  // Fits h1 with func, preforms several checks on the quality of the
+  // fit and tries to improve it. If the fit is not good enough, it
+  // returs false.
+
+  Double_t amin; Double_t edm; Double_t errdef; Int_t nvpar; Int_t nparx;
+  TVirtualFitter *fitter;
+  cout << "--- Fitting : " << h1->GetName() << " ["<< h1->GetTitle() <<"] ---"<< endl;
+
+  h1-> Fit(func,"IME0","",min,max);      
+  Int_t fitResult = h1-> Fit(func,"IME0","",min,max);      
+//   h1-> Fit(func,"0","",min,max);      
+//   Int_t fitResult = h1-> Fit(func,"0IE","",min,max);      
+  
+
+// From TH1:
+// The fitStatus is 0 if the fit is OK (i.e no error occurred).  The
+// value of the fit status code is negative in case of an error not
+// connected with the minimization procedure, for example when a wrong
+// function is used.  Otherwise the return value is the one returned
+// from the minimization procedure.  When TMinuit (default case) or
+// Minuit2 are used as minimizer the status returned is : fitStatus =
+// migradResult + 10*minosResult + 100*hesseResult +
+// 1000*improveResult.  TMinuit will return 0 (for migrad, minos,
+// hesse or improve) in case of success and 4 in case of error (see
+// the documentation of TMinuit::mnexcm). So for example, for an error
+// only in Minos but not in Migrad a fitStatus of 40 will be returned.
+// Minuit2 will return also 0 in case of success and different values
+// in migrad minos or hesse depending on the error. See in this case
+// the documentation of Minuit2Minimizer::Minimize for the
+// migradResult, Minuit2Minimizer::GetMinosError for the minosResult
+// and Minuit2Minimizer::Hesse for the hesseResult.  If other
+// minimizers are used see their specific documentation for the status
+// code returned.  For example in the case of Fumili, for the status
+// returned see TFumili::Minimize.
+
+  if( gMinuit->fLimset ) {
+    Printf("ERROR: AliBWTools: Parameters at limits");
+    return kFALSE;
+  } 
+
+
+  ///
+  fitter = TVirtualFitter::GetFitter();   
+  Int_t  fitStat = fitter->GetStats(amin, edm, errdef, nvpar, nparx);  
+
+  if( ( (fitStat < 3  && gMinuit->fCstatu != "UNCHANGED ")|| (edm > 1e6) || (fitResult !=0 && fitResult < 4000) ) && 
+      TString(gMinuit->fCstatu) != "SUCCESSFUL"  &&
+      TString(gMinuit->fCstatu) != "CONVERGED "  ) {
+    if(fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ") {
+      Printf("WARNING: AliBWTools: Cannot properly compute errors");
+      if (fitStat == 0) Printf(" not calculated at all");
+      if (fitStat == 1) Printf(" approximation only, not accurate");
+      if (fitStat == 2) Printf(" full matrix, but forced positive-definite");
+    }
+
+    if (edm > 1e6) {
+
+      Printf("WARNING: AliBWTools: Huge EDM  (%f): Fit probably failed!", edm);
+    }
+    if (fitResult != 0) {
+      Printf("WARNING: AliBWTools: Fit Result (%d)", fitResult);
+    }
+      
+    Printf ("AliBWTools: Trying Again with Strategy = 2");
+
+    gMinuit->Command("SET STRATEGY 2"); // more effort
+    fitResult = 0;
+    fitResult = h1-> Fit(func,"0","",min,max);      
+    fitResult = h1-> Fit(func,"IME0","",min,max);      
+    fitResult = h1-> Fit(func,"IME0","",min,max);      
+      
+    fitter = TVirtualFitter::GetFitter();   
+  
+    fitStat = fitter->GetStats(amin, edm, errdef, nvpar, nparx);  
+
+    if(fitStat < 3 && gMinuit->fCstatu != "UNCHANGED ") {
+      Printf("ERROR: AliBWTools: Cannot properly compute errors");
+      if (fitStat == 0) Printf(" not calculated at all");
+      if (fitStat == 1) Printf(" approximation only, not accurate");
+      if (fitStat == 2) Printf(" full matrix, but forced positive-definite");
+      cout << "[" <<gMinuit->fCstatu<<"]" << endl;
+      return kFALSE;
+    }
+
+    if (edm > 1e6) {
+      Printf("ERROR: AliBWTools: Huge EDM  (%f): Fit probably failed!", edm);
+
+      return kFALSE;
+    }
+    if (fitResult != 0) {
+      Printf("ERROR: AliBWTools: Fit Result (%d)", fitResult);
+
+      return kFALSE;
+    }
+      
+    gMinuit->Command("SET STRATEGY 1"); // back to normal value
+
+  }
+
+  cout << "---- FIT OK ----" << endl;
+  
+  return kTRUE;
+  
+}
+
+Int_t AliBWTools::GetLowestNotEmptyBin(TH1*h) {
+
+  // Return the index of the lowest non empty bin in the histo h
+
+  Int_t nbin = h->GetNbinsX();
+  for(Int_t ibin = 1; ibin <= nbin; ibin++){
+    if(h->GetBinContent(ibin)>0) return ibin;
+  }
+  
+  return -1;
+
+}
+
+Int_t AliBWTools::GetHighestNotEmptyBin(TH1*h) {
+
+  // Return the index of the highest non empty bin in the histo h
+
+  Int_t nbin = h->GetNbinsX();
+  for(Int_t ibin = nbin; ibin > 0; ibin--){
+    if(h->GetBinContent(ibin)>0) return ibin;
+  }
+  
+  return -1;
+
+}
+
+void AliBWTools::GetResiduals(TGraphErrors * gdata, TF1 * func, TH1F ** hres, TGraphErrors ** gres) {
+
+  // Returns a graph of residuals vs point and the res/err distribution
+
+  Int_t npoint = gdata->GetN();
+
+  (*gres) =new TGraphErrors();
+  (*hres) = new TH1F(TString("hres_")+gdata->GetName()+"-"+func->GetName(),
+                  TString("hres_")+gdata->GetName()+"-"+func->GetName(),
+                  20,-5,5);
+
+
+  for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
+    Float_t x   = gdata->GetX()[ipoint];
+    Float_t res = (gdata->GetY()[ipoint] - func->Eval(x))/func->Eval(x);
+    Float_t err = gdata->GetEY()[ipoint]/func->Eval(x);
+    (*hres)->Fill(res/err);
+    (*gres)->SetPoint(ipoint, x, res/err);
+    //    (*gres)->SetPointError(ipoint, 0, err);
+    
+  }
+  
+  (*gres)->SetMarkerStyle(gdata->GetMarkerStyle());
+  (*gres)->SetMarkerColor(gdata->GetMarkerColor());
+  (*gres)->SetLineColor  (gdata->GetLineColor());
+  (*gres)->GetHistogram()->GetYaxis()->SetTitle("(data-function)/function");
+  (*hres)->SetMarkerStyle(gdata->GetMarkerStyle());
+  (*hres)->SetMarkerColor(gdata->GetMarkerColor());
+  (*hres)->SetLineColor  (gdata->GetLineColor());
+
+
+
+}
+
+void AliBWTools::GetResiduals(TH1F* hdata, TF1 * func, TH1F ** hres, TH1F ** hresVsBin) {
+
+  // Returns an histo of residuals bin by bin and the res/err distribution
+
+  if (!func) {
+    Printf("AliBWTools::GetResiduals: No function provided");
+    return;
+  }
+  if (!hdata) {
+    Printf("AliBWTools::GetResiduals: No data provided");
+    return;
+  }
+
+  (*hresVsBin) = (TH1F*) hdata->Clone(TString("hres_")+hdata->GetName());
+  (*hresVsBin)->Reset();
+  (*hres) = new TH1F(TString("hres_")+hdata->GetName()+"-"+func->GetName(),
+                    TString("hres_")+hdata->GetName()+"-"+func->GetName(),
+                    20,-5,5);
+
+  Int_t nbin = hdata->GetNbinsX();
+  for(Int_t ibin = 1; ibin <= nbin; ibin++){
+    if(hdata->GetBinContent(ibin)==0) continue;
+    Float_t res = (hdata->GetBinContent(ibin) - func->Eval(hdata->GetBinCenter(ibin)) ) / 
+      func->Eval(hdata->GetBinCenter(ibin));
+    Float_t err = hdata->GetBinError  (ibin) /  func->Eval(hdata->GetBinCenter(ibin));
+    (*hresVsBin)->SetBinContent(ibin,res);
+    (*hresVsBin)->SetBinError  (ibin,err);
+    (*hres)->Fill(res/err);
+    
+  }
+  
+  (*hresVsBin)->SetMarkerStyle(hdata->GetMarkerStyle());
+  (*hresVsBin)->SetMarkerColor(hdata->GetMarkerColor());
+  (*hresVsBin)->SetLineColor  (hdata->GetLineColor()  );
+  (*hresVsBin)->GetYaxis()->SetTitle("(data-function)/function");
+  (*hres)->SetMarkerStyle(hdata->GetMarkerStyle());
+  (*hres)->SetMarkerColor(hdata->GetMarkerColor());
+  (*hres)->SetLineColor  (hdata->GetLineColor()  );
+
+}
+
+void AliBWTools::GetYield(TH1* h, TF1 * f, Double_t &yield, Double_t &yieldError, Float_t min, Float_t max,
+                         Double_t *partialYields, Double_t *partialYieldsErrors){
+
+  // Returns the yield extracted from the data in the histo where
+  // there are points and from the fit to extrapolate, in the given
+  // range.
+
+  // Partial yields are also returned if the corresponding pointers are non null
+
+  Int_t bin1 = h->FindBin(min);
+  Int_t bin2 = h->FindBin(max);
+  Float_t bin1Edge = GetLowestNotEmptyBinEdge (h);
+  Float_t bin2Edge = GetHighestNotEmptyBinEdge(h);
+
+  Double_t integralFromHistoError ;
+  Double_t integralFromHisto = h->IntegralAndError(bin1,bin2,integralFromHistoError,"width");
+  
+  Double_t integralBelow      = min < bin1Edge ? f->Integral(min,bin1Edge) : 0;
+  Double_t integralBelowError = min < bin1Edge ? f->IntegralError(min,bin1Edge) : 0;
+  Double_t integralAbove      = max > bin2Edge ? f->Integral(bin2Edge,max) : 0;
+  Double_t integralAboveError = max > bin2Edge ? f->Integral(bin2Edge,max) : 0;
+
+  cout << "GetYield INFO" << endl;
+  cout << " " << bin1Edge << " " << bin2Edge << endl;  
+  cout << " " << integralFromHisto      << " " << integralBelow      << " " << integralAbove      << endl;
+  cout << " " << integralFromHistoError << " " << integralBelowError << " " << integralAboveError << endl;
+  
+  if(partialYields) {
+    partialYields[0] = integralFromHisto;
+    partialYields[1] = integralBelow;
+    partialYields[2] = integralAbove;
+  }
+  if(partialYieldsErrors) {
+    partialYieldsErrors[0] = integralFromHistoError;
+    partialYieldsErrors[1] = integralBelowError;
+    partialYieldsErrors[2] = integralAboveError;
+  }
+  yield      = integralFromHisto+integralBelow+integralAbove;
+  yieldError = TMath::Sqrt(integralFromHistoError*integralFromHistoError+
+                          integralBelowError*integralBelowError+
+                          integralAboveError*integralAboveError);
+
+}
+
+TGraphErrors * AliBWTools::DivideGraphByFunc(TGraphErrors * g, TF1 * f, Bool_t invert){ 
+
+  // Divides g/f. If invert == true => f/g
+
+  TGraphErrors * gRatio = new TGraphErrors();
+  Int_t npoint = g->GetN();
+  for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
+    Double_t x = g->GetX()[ipoint];
+    Double_t ratio  = invert ? f->Eval(x)/g->GetY()[ipoint] :g->GetY()[ipoint]/f->Eval(x);
+    gRatio->SetPoint     (ipoint, x, ratio);
+    gRatio->SetPointError(ipoint, 0, g->GetEY()[ipoint]/f->Eval(x));
+    //    cout << x << " " << g->GetY()[ipoint] << " " << f->Eval(x) << endl;
+    
+  }
+  gRatio->SetMarkerStyle(20);
+  //gRatio->Print();
+  return gRatio;
+
+}
+
+TGraphErrors * AliBWTools::DivideGraphByHisto(TGraphErrors * g, TH1 * h, Bool_t invert){ 
+
+  // Divides g/h. If invert == true => h/g
+
+
+  TGraphErrors * gRatio = new TGraphErrors();
+  Int_t npoint = g->GetN();
+  for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
+    Double_t xj  = g->GetX()[ipoint];
+    Double_t yj  = g->GetY()[ipoint];
+    Double_t yje = g->GetEY()[ipoint];
+
+    Int_t binData = h->FindBin(xj);
+    Double_t yd   = h->GetBinContent(binData);
+    Double_t yde  = h->GetBinError(binData);
+    Double_t xd   = h->GetBinCenter(binData);
+    
+    //    cout << TMath::Abs((xd-xj)/xd) << endl;
+    
+
+     
+    if (yd == 0) continue;
+    //    if (binData == 28 ) cout << TMath::Abs(xd-xj)/TMath::Abs(xd) << endl;
+    
+    if (TMath::Abs(xd-xj)/TMath::Abs(xd) > 0.01){
+      Printf( "WARNING: bin center (%f)  and x graph (%f) are more than 1 %% away, skipping",xd,xj );
+      continue;
+      
+    }
+
+    Double_t ratio = invert ? yd/yj : yj/yd;
+
+    gRatio->SetPoint(ipoint, xj, ratio);
+    gRatio->SetPointError(ipoint, 0, TMath::Sqrt(yde*yde/yd/yd + yje*yje/yj/yj)*ratio);
+    //    gRatio->SetPointError(ipoint, 0, yje/yj * ratio);
+  }
+
+  return gRatio;
+
+
+}
+
+TH1F * AliBWTools::DivideHistoByFunc(TH1F * h, TF1 * f, Bool_t invert){ 
+
+  // Divides h/f. If invert == true => f/g
+  // Performs the integral of f on the bin range to perform the ratio
+  // Returns a histo with the same binnig as h
+
+  // Prepare histo for ratio
+  TH1F * hRatio = (TH1F*) h->Clone(TString("hRatio_")+h->GetName()+"_"+f->GetName());
+  hRatio->Reset();
+  // Set y title
+  if(!invert) hRatio->SetYTitle(TString(h->GetName())+"/"+f->GetName());
+  else        hRatio->SetYTitle(TString(f->GetName())+"/"+h->GetName());
+
+  // Loop over all bins
+  Int_t nbin = hRatio->GetNbinsX();
+
+  for(Int_t ibin = 1; ibin <= nbin; ibin++){
+    Double_t yhisto = h->GetBinContent(ibin);
+    Double_t yerror = h->GetBinError(ibin);
+    Double_t xmin   = h->GetBinLowEdge(ibin);
+    Double_t xmax   = h->GetBinLowEdge(ibin+1);
+    Double_t yfunc  = f->Integral(xmin,xmax)/(xmax-xmin);
+    Double_t ratio = invert ? yfunc/yhisto : yhisto/yfunc ;
+    Double_t error = yerror/yfunc  ;
+    hRatio->SetBinContent(ibin,ratio);
+    hRatio->SetBinError  (ibin,error);
+  }
+
+  return hRatio;
+
+}
diff --git a/PWG2/SPECTRA/Fit/AliBWTools.h b/PWG2/SPECTRA/Fit/AliBWTools.h
new file mode 100644 (file)
index 0000000..2f61603
--- /dev/null
@@ -0,0 +1,79 @@
+// ----------------------------------------------------------------------
+//                     AliBWTools
+// 
+// This class provides some tools which can be useful in the analsis
+// of spectra, to fit or transform histograms. See the comments of the
+// individual methods for details
+//
+// Author: M. Floris (CERN)
+// ----------------------------------------------------------------------
+
+#ifndef ALIBWTOOLS_H
+#define ALIBWTOOLS_H
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "TObject.h"
+#include "TH1.h"
+
+class TF1;
+class TH1D;
+class TH1F;
+class TGraphErrors;
+#endif
+
+
+
+class AliBWTools : public TObject {
+
+public:
+
+  AliBWTools();
+  ~AliBWTools();
+
+  static TH1 * GetOneOverPtdNdPt(TH1 * hPt) ;
+  static TH1 * GetdNdmtFromdNdpt(TH1 * hpt, Double_t mass);
+  static TH1 * GetdNdPtFromOneOverPt(TH1 * h1Pt) ;
+
+  static TGraphErrors * ConcatenateGraphs(TGraphErrors * g1,TGraphErrors * g2); 
+  static TH1F *         GetHistoFromGraph(TGraphErrors * g, TH1F* hTemplate) ;  
+  static TGraphErrors * GetGraphFromHisto(TH1F * h, Bool_t binWidth = kTRUE) ; 
+  static TH1F *         CombineHistos(TH1 * h1, TH1 * h2, TH1* htemplate, Float_t renorm1=1.);
+  static void GetFromHistoGraphDifferentX(TH1F * h, TF1 * f, TGraphErrors ** gBarycentre, TGraphErrors ** gXlw); 
+  static Float_t GetMean(TH1F * h, Float_t min, Float_t max) ; 
+
+  static void GetMean(TF1 * func, Float_t &mean, Float_t &error, Float_t min=0, Float_t max=100);
+  static void GetMeanSquare(TF1 * func, Float_t &mean, Float_t &error, Float_t min=0, Float_t max=100) ;
+
+  static Bool_t Fit (TH1 * h, TF1* f, Float_t min, Float_t max) ;
+
+  static Int_t GetLowestNotEmptyBin(TH1*h);
+  static Int_t GetHighestNotEmptyBin(TH1*h);
+  static Float_t GetLowestNotEmptyBinEdge(TH1*h) { return h->GetBinLowEdge(GetLowestNotEmptyBin(h));}
+  static Float_t GetHighestNotEmptyBinEdge(TH1*h) { return h->GetBinLowEdge(GetHighestNotEmptyBin(h)+1);}
+
+  static void GetResiduals(TGraphErrors * gdata, TF1 * func, TH1F ** hres, TGraphErrors ** gres) ;
+  static void GetResiduals(TH1F* hdata, TF1 * func, TH1F ** hres, TH1F ** hresVsBin) ;
+
+  static void GetYield(TH1* h, TF1 * f, Double_t &yield, Double_t &yieldError, Float_t min = 0, 
+                      Float_t max = 100,  Double_t *partialYields=0, Double_t *partialYieldsErrors=0);
+
+  static TGraphErrors * DivideGraphByFunc (TGraphErrors * g, TF1 * f, Bool_t invert = kFALSE);
+  static TGraphErrors * DivideGraphByHisto(TGraphErrors * g, TH1 * h, Bool_t invert = kFALSE);
+  static TH1F         * DivideHistoByFunc (TH1F * h, TF1 * f, Bool_t invert = kFALSE);
+
+private:
+
+  static void GetMoment(TString name, TString var, TF1 * func, Float_t &mean, Float_t &error, Float_t min, Float_t max) ;
+  static Double_t GetNormalizedFunc(double * x, double* p);
+
+  static TF1 * fFuncForNormalized; // Function used in GetNormalizedFunc
+
+  AliBWTools(const AliBWTools&);            // not implemented
+  AliBWTools& operator=(const AliBWTools&); // not implemented
+
+  ClassDef(AliBWTools,1);
+
+};
+
+#endif
diff --git a/PWG2/SPECTRA/Fit/AliLatexTable.cxx b/PWG2/SPECTRA/Fit/AliLatexTable.cxx
new file mode 100644 (file)
index 0000000..771abfd
--- /dev/null
@@ -0,0 +1,463 @@
+// ------------------------------------------------------------------
+//
+//                           AliLatexTable
+//
+// This class is used to produce tables using latex sintax. It can
+// than print the table as Latex source code (to be pasted on a paper)
+// or ASCII.
+//
+// The basic idea is that you add columns one after the other with the
+// SetNextCol method and than you insert this row. The SetNextCol
+// method comes in different flavours to alow the insertion of
+// different type of values (numbers with or without errors,
+// strings...).
+//
+// TODO:
+// 1. Make the class drawable
+// 2. Implement vertical lines in ascii print
+// 3. Print output in HTML format
+//
+// Author: Michele Floris, CERN
+// ------------------------------------------------------------------
+
+#include "AliLatexTable.h"
+#include "TString.h"
+#include "TRegexp.h"
+#include "TObjString.h"
+#include "TObjArray.h"
+#include <stdarg.h>
+#include "snprintf.h"
+#include "Varargs.h"
+#include "TPRegexp.h"
+#include "TMath.h"
+#include <iostream>
+
+using namespace std;
+
+ClassImp(AliLatexTable);
+
+AliLatexTable::AliLatexTable(Int_t ncol, TString format) : fNcol(0), fFormat(""), fRows(0), fCols(0), fNcolReady(0){
+  // constructor, specify number of cols 
+  fNcol = ncol;
+  fFormat = format;
+  fNcolReady = 0;
+  fRows = new TObjArray;
+  fCols = new TObjArray;
+
+  fRows->SetOwner();
+  fCols->SetOwner();
+}
+
+
+AliLatexTable::~AliLatexTable() {
+
+  // dtor
+  if (fRows) delete fRows;
+  if (fCols) delete fCols;
+
+}
+
+void AliLatexTable::SetNextCol(Int_t val){
+  // Set next column in current row - integer
+  char col[200];
+  sprintf(col, " %d ", val);
+  SetNextCol(col);
+
+} 
+
+void AliLatexTable::SetNextCol(Int_t val, Int_t err){
+
+  // Set next column in current row - int +- int
+  char col[200];
+  sprintf(col, " $%d \\pm %d$ ", val, err);
+  SetNextCol(col);
+
+} 
+
+void AliLatexTable::SetNextCol(Double_t val, Int_t scientificNotation){
+
+  // Set next column in current row - double, specify resolution
+
+  char col[200];
+  if(scientificNotation >= 0) {
+    char format[100];
+
+    //    cout << format << endl;    
+    Double_t mantissa, exp;
+    GetMantissaAndExpBase10(val, mantissa, exp);
+    
+    if (exp == 0) {
+      sprintf(format, " $%%%d.%df $", scientificNotation, scientificNotation);
+      sprintf(col, format, mantissa);      
+    }
+    else  {
+      sprintf(format, " $%%%d.%df \\cdot 10^{%%0.0f} $", scientificNotation, scientificNotation);
+      sprintf(col, format, mantissa, exp);
+    }
+
+
+    
+
+    SetNextCol(col);
+
+  } else if (scientificNotation < -1) {
+    char format[100];
+    sprintf(format, " $%%%d.%df $", -scientificNotation,-scientificNotation);    
+    sprintf(col, format , val);
+    SetNextCol(col);
+  } else {
+    sprintf(col, " %f ", val);
+    SetNextCol(col);
+  }
+
+} 
+void AliLatexTable::SetNextCol(Double_t val, Double_t err, Int_t scientificNotation){
+
+  // Set next column in current row - double +- double, specify resolution
+
+  // if scientific notation is != -1 is used to determine number of
+  // digits in error 
+
+  //if it is > 0 exp notation is used 
+
+  //if it is < -1 number is truncated to the number of digits after
+  //point dictated by scientificNotation,
+
+  char col[200];
+  if(scientificNotation >=0 ) {
+
+    Double_t mantissa, exp;
+    GetMantissaAndExpBase10(val, mantissa, exp);
+
+    Double_t mantissa_err, exp_err;
+    GetMantissaAndExpBase10(err, mantissa_err, exp_err);
+
+    Int_t nSigDigits =TMath::Nint(exp - exp_err);
+    if(scientificNotation != 0) nSigDigits =  nSigDigits + scientificNotation - 1;
+
+
+    char format[100];
+    if (exp == 0) {
+      sprintf(format, " $%%%d.%df \\pm %%%d.%df $", nSigDigits, nSigDigits, nSigDigits, nSigDigits);
+      sprintf(col, format, mantissa, mantissa_err/TMath::Power(10,exp - exp_err));
+    }
+    else  {
+      sprintf(format, " $%%%d.%df \\pm %%%d.%df \\cdot 10^{%%0.0f}$", nSigDigits, nSigDigits, nSigDigits, nSigDigits);
+      sprintf(col, format, mantissa, mantissa_err/TMath::Power(10,exp - exp_err), exp);
+    }
+
+
+ //cout << format << endl;
+    
+    SetNextCol(col);
+
+  } else if (scientificNotation < -1) {
+    char format[100];
+    sprintf(format, " $%%%d.%df \\pm %%%d.%df $", -scientificNotation,-scientificNotation,-scientificNotation,-scientificNotation);    
+    sprintf(col, format , val, err);
+    SetNextCol(col);
+  }
+  else {
+    sprintf(col, " $%f \\pm %f$ ", val, err);
+    SetNextCol(col);
+  }
+} 
+
+void AliLatexTable::SetNextCol(Double_t val, Double_t err, Double_t errSyst, Int_t scientificNotation){
+
+  // Set next column in current row - double +- double +- double, specify resolution
+
+  // if scientific notation is != -1 is used to determine number of
+  // digits in error 
+
+  //if it is > 0 exp notation is used 
+
+  //if it is < -1 number is truncated to the number of digits after
+  //point dictated by scientificNotation,
+
+  char col[200];
+  if(scientificNotation >=0 ) {
+
+    Double_t mantissa, exp;
+    GetMantissaAndExpBase10(val, mantissa, exp);
+
+    Double_t mantissa_err, exp_err;
+    GetMantissaAndExpBase10(err, mantissa_err, exp_err);
+
+    Double_t mantissa_errsyst, exp_errsyst;
+    GetMantissaAndExpBase10(errSyst, mantissa_errsyst, exp_errsyst);
+
+    Int_t nSigDigits =TMath::Nint(exp - exp_err);
+    if(scientificNotation != 0) nSigDigits =  nSigDigits + scientificNotation - 1;
+
+
+    char format[100];
+    if (exp == 0) {
+      sprintf(format, " $%%%d.%df \\pm %%%d.%df \\pm %%%d.%df $", 
+             nSigDigits, nSigDigits, nSigDigits, nSigDigits, nSigDigits, nSigDigits);
+      sprintf(col, format, mantissa, 
+             mantissa_err/TMath::Power(10,exp - exp_err), 
+             mantissa_errsyst/TMath::Power(10,exp - exp_errsyst));
+    }
+    else  {
+      sprintf(format, " $%%%d.%df \\pm %%%d.%df  \\pm %%%d.%df \\cdot 10^{%%0.0f}$", 
+             nSigDigits, nSigDigits, nSigDigits, nSigDigits, nSigDigits, nSigDigits);
+      sprintf(col, format, mantissa, 
+             mantissa_err/TMath::Power(10,exp - exp_err), 
+             mantissa_errsyst/TMath::Power(10,exp - exp_errsyst),
+             exp);
+    }
+
+
+ //cout << format << endl;
+    
+    SetNextCol(col);
+
+  } else if (scientificNotation < -1) {
+    char format[100];
+    sprintf(format, " $%%%d.%df \\pm %%%d.%df \\pm %%%d.%df $", -scientificNotation,-scientificNotation,-scientificNotation,-scientificNotation, -scientificNotation,-scientificNotation);    
+    sprintf(col, format , val, err);
+    SetNextCol(col);
+  } 
+  else {
+    sprintf(col, " $%f \\pm %f  \\pm %f$ ", val, err, errSyst);
+    SetNextCol(col);
+  }
+} 
+
+// void AliLatexTable::SetNextColPrintf(const char *va_(fmt), ...){
+
+
+//   const Int_t buf_size = 1024;  // hope it's long enough
+//   char colBuffer[buf_size];
+
+//   va_list ap;
+//   va_start(ap, va_(fmt));
+//   Int_t n = vsnprintf(colBuffer, buf_size, fmt, ap);
+
+//   if (n == -1 || n >= buf_size) {
+//     Error("SetNextCol", "vsnprintf error!");
+//   }
+
+//   va_end(ap);
+     
+// //   for(int iobj=0; iobj<num; iobj++){   //Loop until all objects are added
+// //     TH1 * obj = (TH1 *) va_arg(arguments, void *);
+// //     returnValue = returnValue && AddObject(obj,(char*)fOptMany.Data(), EColor(iobj+1));
+// //   }
+//   SetNextCol(colBuffer);
+// }
+
+
+void AliLatexTable::SetNextCol(TString val){
+
+  // Set next column in current row - tstring
+
+  fCols->Add(new TObjString(val));
+  fNcolReady++;
+
+} 
+
+void AliLatexTable::InsertCustomRow(TString row){
+  // insert a full row from string. Make sure you have all the columns
+  fRows->Add(new TObjString(row));
+}
+
+
+
+void AliLatexTable::InsertRow(){
+
+  // Insert the row, based on all the individual columns
+
+  if ( fNcolReady != fNcol) {
+    Warning("InsertRow", "Wrong number of cols: %d (!= %d)", fNcolReady, fNcol); 
+  }
+
+  TString row = "";
+  for(Int_t icol = 0; icol < fNcol; icol++){
+    row = row + ((TObjString*) fCols->At(icol))->String();
+    if(icol != (fNcol-1)) row = row + " & ";
+  }
+  row+="\\\\";
+
+  fRows->Add(new TObjString(row));
+
+  fNcolReady = 0;
+  fCols->Clear();
+
+}
+
+void AliLatexTable::InsertHline(){
+
+  // insert an horizontal line
+  fRows->Add(new TObjString("\\hline"));
+  
+
+}
+
+Int_t * AliLatexTable::GetColWidths() {
+
+  // Compute the width of columns, for nice ascii printing
+
+  Int_t * col_widths= new Int_t[fNcol];
+  Int_t nrow = fRows->GetEntriesFast();
+
+  for(Int_t icol = 0; icol < fNcol; icol++){
+    col_widths[icol] = 0;
+  }
+  
+
+  for(Int_t irow = 0; irow < nrow; irow++){
+    TString row(((TObjString*) fRows->At(irow))->String());
+    if(row.Contains("\\hline"))continue;
+    
+    StripLatex(row, "ASCII");
+
+    TObjArray * cols = row.Tokenize("&");
+    if(cols->GetEntries() != fNcol) {
+      Error("GetColWidths", "Wrong number of cols in row %s: %d - %d", row.Data(), cols->GetEntries(), fNcol);      
+    }
+    for(Int_t icol = 0; icol < fNcol; icol++){
+      Int_t w = ((TObjString *) cols->At(icol))->String().Length();
+      if (w>col_widths[icol]) col_widths[icol] = w;
+    }
+    
+  }
+
+  return col_widths;
+}
+
+void AliLatexTable::PrintTable(Option_t * opt){
+
+  // Print the table on screen. You can specify the format with opt.
+  // Currently supported:
+  // "" -> LaTeX
+  // "ASCII" -> plain ASCII
+  // "HTML"  -> HTML, to be improved
+
+  if(TString(opt) == "ASCII" || TString(opt)=="HTML") {
+    
+    Int_t nrow = fRows->GetEntriesFast();
+
+    Int_t * col_widths = GetColWidths();
+
+    Int_t total_lenght = 0;
+    for(Int_t icol = 0; icol < fNcol; icol++) total_lenght = total_lenght + col_widths[icol] + 2 ;
+    
+
+    for(Int_t irow = 0; irow < nrow; irow++){
+      TString row = ((TObjString*) fRows->At(irow))->String();
+      if (row.Contains("\\hline") ){
+       for(Int_t il = 0; il < total_lenght; il++) printf("-");
+       printf("\n");     
+       continue;
+      }
+      StripLatex(row, opt);
+      TObjArray * cols = row.Tokenize("&");
+      for(Int_t icol = 0; icol < fNcol; icol++){
+       const char * colstr = ((TObjString *) cols->At(icol))->String().Data();
+       char format [200];
+       sprintf(format, "%%%ds", col_widths[icol] + 2); 
+       printf(format, colstr);
+
+      }
+      printf ("\n");
+    }
+    
+
+    return;
+  }
+  
+
+  cout << "\\begin{tabular}{"<<fFormat<<"}"<<endl;
+
+  Int_t nrow = fRows->GetEntriesFast();
+
+  for(Int_t irow = 0; irow < nrow; irow++){
+    cout << ((TObjString*) fRows->At(irow))->String() << endl;
+  }
+  
+  cout << "\\end{tabular}" << endl;
+
+
+}
+
+// void AliLatexTable::ParseExponent(TString &expo){
+// //    TString parseExponent = col;
+//   TRegexp re = "e[+-][0-9][0-9]";
+//   //  cout << col << endl;
+  
+//   if(expo.Contains(re)){
+//     Int_t index = expo.Index(re);
+//     TString replacement = "\\cdot 10^{";
+//     replacement = replacement + expo(index+1, 3) +"}";
+//     //    cout << replacement <<" --- "<< endl;
+//     replacement.ReplaceAll("+","");
+// //     cout << "B: " << expo << endl;      
+// //     cout << "RE " << expo(re) << endl;
+    
+//     expo.ReplaceAll(expo(re), replacement);
+//     //    cout << "A: " << expo << endl;      
+//     // recursion to parse all exponents
+//     if (expo.Contains(re)) ParseExponent(expo);
+//   }
+//   else Warning("", "Error parsing exponent");
+// }
+
+void AliLatexTable::GetMantissaAndExpBase10(Double_t num, Double_t &man, Double_t &exp) {
+
+  // Helper used to get mantissa and exponent
+
+  exp = TMath::Floor(TMath::Log10(TMath::Abs(num)));  
+  man = num / TMath::Power(10, exp);
+
+//   cout << "" << endl;
+//   cout << num << " = " << man << " x10^{"<<exp<<"} "   << endl;
+  
+
+}
+
+void AliLatexTable::StripLatex(TString &text, TString format) {
+
+  // Strip latex away for ascii and html printing. Replaces latex
+  // command with corresponding text/tags
+
+  text.ReplaceAll("\\cdot", "x");
+  text.ReplaceAll("\\pm", "+-");
+  text.ReplaceAll("$", "");
+  if (format == "ASCII") {
+    text.ReplaceAll("\\right>", ">");
+    text.ReplaceAll("\\left<", "<");
+    text.ReplaceAll("\\rangle", ">");
+    text.ReplaceAll("\\langle", "<");
+  } else if (format == "HTML") {
+    text.ReplaceAll("\\right>", "&rang;");
+    text.ReplaceAll("\\left<",  "&lang;");
+    text.ReplaceAll("\\rangle", "&rang;");
+    text.ReplaceAll("\\langle", "&lang;");
+  } 
+  if(text.Contains("multicolumn")) {
+    //    cout << "col " << text.Data() << endl;
+    // in case of multicol span, replace first column with content and
+    // add n empty cols
+    TObjArray * array = TPRegexp("multicolumn\{(.*)\\}\{.*\\}\{(.*)\\}").MatchS(text); // Added double \\ because gcc 4 triggers hard error on unknown escape sequence. Hope it still works...
+    const TString content = ((TObjString *)array->At(2))->GetString();
+    Int_t nspan   = ((TObjString *)array->At(1))->GetString().Atoi();
+    text = content;
+//     cout << "ns " << nspan <<  ((TObjString *)array->At(1))->GetString() << endl;
+//     cout << "t " << text << endl;
+    
+    for(Int_t ispan = 1; ispan < nspan; ispan++){
+      text+=" & ";
+    }
+    //    cout << "t " << text << endl;
+    
+  
+  }
+
+  text.ReplaceAll("\\mathrm", "");
+
+  text.ReplaceAll("\\", "");
+  text.Strip(TString::EStripType(1), ' ');
+
+}
diff --git a/PWG2/SPECTRA/Fit/AliLatexTable.h b/PWG2/SPECTRA/Fit/AliLatexTable.h
new file mode 100644 (file)
index 0000000..4eca83f
--- /dev/null
@@ -0,0 +1,78 @@
+// This class is used to print a table which can be pasted on a latex
+// document
+
+#ifndef ALILATEXTABLE_H
+#define ALILATEXTABLE_H
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "TObject.h"
+#include "TString.h"
+class TObjArray;
+
+#endif
+
+
+using namespace std;
+
+class AliLatexTable : public TObject {
+
+public:
+  AliLatexTable() : fNcol(1), fFormat("c"), fRows(0), fCols(0), fNcolReady(0) {;}
+  AliLatexTable(Int_t ncol, TString format);
+  ~AliLatexTable();
+
+  // first you set the value of each column
+  void SetNextCol(Int_t val); 
+  void SetNextCol(Int_t val, Int_t err); 
+
+  void SetNextCol(Double_t val, Int_t scientificNotation = -1); // if different from -1 gives significant digits
+  void SetNextCol(Double_t val, Double_t err, Int_t scientificNotation = -1); 
+  void SetNextCol(Double_t val, Double_t err, Double_t errSyst, Int_t scientificNotation = -1); 
+
+  void SetNextCol(TString val); 
+//   // allows to use printf syntax
+//   void SetNextColPrintf(const char *va_(fmt), ...);
+
+  // Then you add the row (it's up to user make sure all columns are
+  // there)
+  void InsertRow();
+
+
+  // insert a row without building it up with the methods. May be
+  // usefull for header or multcol rows
+  void InsertCustomRow(TString row);
+
+
+  void InsertHline();
+
+  void PrintTable(Option_t * opt = "");
+  //  void ParseExponent(TString &expo);
+  void GetMantissaAndExpBase10(Double_t num, Double_t &man, Double_t &exp) ;
+
+  // used to print columns with correct width
+  Int_t * GetColWidths();
+  void StripLatex(TString &row, TString format) ;
+  
+
+private:
+
+  Int_t fNcol;     // number of columns
+  TString fFormat; // latex format (es "c|ccc")
+  
+  TObjArray * fRows; // rows 
+  TObjArray * fCols; // columns
+  
+  Int_t fNcolReady; // number of cols ready to be insert
+
+
+  AliLatexTable(const AliLatexTable&);
+  AliLatexTable& operator=(const AliLatexTable&);
+
+  ClassDef(AliLatexTable, 1)
+
+
+};
+
+#endif
+
diff --git a/PWG2/SPECTRA/Fit/CombineSpectra.C b/PWG2/SPECTRA/Fit/CombineSpectra.C
new file mode 100644 (file)
index 0000000..cb4e76f
--- /dev/null
@@ -0,0 +1,1216 @@
+#include <iostream>
+#include "TH1F.h"
+#include "TGraphErrors.h"
+#include "AliBWFunc.h"
+#include "AliBWTools.h"
+#include "TF1.h"
+#include "TFile.h"
+#include "TDatabasePDG.h"
+#include "TROOT.h"
+#include "TCanvas.h"
+#include "TFolder.h"
+#include "TStyle.h"
+#include "AliLatexTable.h"
+#include "TLegend.h"
+#include "TVirtualFitter.h"
+#include "TMath.h"
+#include "TH2F.h"
+#include "TSystem.h"
+#include "TLine.h"
+#include "TLatex.h"
+#include "TMath.h"
+
+#include "TASImage.h"
+#include "TPaveText.h"
+#include "ALICEWorkInProgress.C"
+#include "StarPPSpectra.C"
+#include "GetE735Ratios.C"
+
+using namespace std;
+
+
+// A bunch of useful enums and constants
+enum {kPion=0,kKaon,kProton,kNPart};
+enum {kTPC=0,kTOF,kITS,kITSTPC,kK0,kKinks,kCombTOFTPC,kNHist};// "k0" listed here as a kind of PID method...
+const Int_t kNDet = kITS+2;
+enum {kPos=0,kNeg,kNCharge};
+enum {kPhojet=0,kPyTuneAtlasCSC, kPyTuneCMS6D6T, kPyTunePerugia0, kNTunes} ;
+enum {kFitLevi=0, kFitUA1, kFitPowerLaw,
+      kFitPhojet, kFitAtlasCSC, kFitCMS6D6T, kFitPerugia0,
+      kNFit};
+
+
+// flags, labels and names
+const char * partFlag[] = {"Pion", "Kaon", "Proton"};
+const char * detFlag[]  = {"TPC", "TOF", "ITS", "ITS Global", "K0", "Kinks", "Combined TOF + TPC"};
+const char * chargeFlag[]  = {"Pos", "Neg"};
+const char * chargeLabel[]  = {"Positive", "Negative"};
+const char * partLabel[kNPart][kNCharge] = {{"#pi^{+}", "#pi^{-}"}, 
+                                           //                                      {"K^{+} (#times 2)", "K^{-} (#times 2)"}, 
+                                           {"K^{+}", "K^{-}"}, 
+                                           {"p" ,  "#bar{p}"}};
+const char * mcTuneName[] = {"Phojet", 
+                            "Pythia - CSC 306", 
+                            "Pythia - D6T 109", 
+                            "Pythia - Perugia0 - 320", };
+const char * funcName[] = { "Levi", "UA1", "PowerLaw", "Phojet", "AtlasCSC", "CMS6D6T", "Perugia0"};
+
+// Style
+//const Int_t marker[] = {25,24,28,20,21}; // set for highlithining marek
+const Int_t marker[] = {20,24,25,28,21}; // standard set
+const Int_t color [] = {1,2,4};
+
+const Int_t mcLineColor[] = {kGreen,kRed,kBlue,kBlack};
+const Int_t mcLineStyle[] = {1,2,3,4};
+
+
+// Template needed to combine different detectors
+const Float_t templBins[] = {0.05,0.1,0.15,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,2,2.2,2.4,2.6};
+Int_t nbinsTempl=31;
+
+TH1F * htemplate = new TH1F("htemplate", "htemplate",nbinsTempl, templBins );
+
+//  Globals
+TH1F * hSpectra[kNHist][kNPart][kNCharge];
+TH1F * hSpectraMC[kNTunes][kNPart][kNCharge];
+Double_t mass[kNPart];
+
+// Functions:
+// Loading
+void LoadSpectra() ;
+void LoadMC() ;
+
+// Additional tasks (may be uncommented below)
+void DrawStar(Int_t icharge);
+void GetITSResiduals();
+void DrawWithModels() ;
+void DrawAllAndKaons();
+void DrawWithJacek();
+void DrawRatioToStar();
+void DrawRatios();
+
+// External stuff
+//void ALICEWorkInProgress(TCanvas *c,TString today, TString label);
+
+// Used to tag plots
+TDatime dt;
+TString today = "";
+
+
+
+// Switches
+Bool_t convertToMT = 0;
+Bool_t doPrint = 0;
+Int_t  fitFuncID = kFitLevi;
+Bool_t scaleKaons =  kFALSE;
+Bool_t correctSecondaries  = 1;
+Bool_t correctGeantFlukaXS = 1;
+
+void CombineSpectra() {
+
+  // This macro is used to combine the 900 GeV spectra from 2009
+
+  // Load stuff
+  gSystem->Load("libTree.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libMinuit.so");
+  gSystem->Load("libSTEERBase.so");
+  gSystem->Load("libESD.so");
+  gSystem->Load("libAOD.so");
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libCORRFW.so");
+  gSystem->Load("libPWG2spectra.so");
+
+  // Set date
+  today = today +  long(dt.GetDay()) +"/" + long(dt.GetMonth()) +"/"+ long(dt.GetYear());
+
+
+  // Set Masses
+  mass[kPion]   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
+  mass[kKaon]   = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
+  mass[kProton] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
+
+  // Load histos
+  LoadSpectra();
+  LoadMC();
+
+  // Additional tasks
+  //DrawStar(icharge);
+  //  GetITSResiduals();
+  //DrawAllAndKaons();  
+  //  DrawWithModels() ;
+  //DrawWithJacek();
+  //DrawRatioToStar();
+  //DrawRatios();
+  //return;
+
+
+  // Draw combined & Fit
+  AliBWFunc * fm = new AliBWFunc;
+  fm->SetVarType(AliBWFunc::kdNdpt);
+  if (convertToMT) fm->SetVarType(AliBWFunc::kOneOverMtdNdmt);
+
+  // table to print results
+  AliLatexTable table(10,"c|cccccccc");
+  if (fitFuncID == kFitLevi) {
+    table.InsertCustomRow("Part & Yield & Yield (FIT) &  T Slope & n & $\\Chi^2$/NDF & Min X & Frac Above & \\langle p_{t} \\rangle  & \\langle p_{t}^{2} \\rangle");
+  }  else if (fitFuncID == kFitPowerLaw) {
+    table.InsertCustomRow("Part & Yield & Norm &  n & pt0 & $\\Chi^2$/NDF & Min X & Frac Above & \\langle p_{t} \\rangle  & \\langle p_{t}^{2} \\rangle");    
+  } else {
+    table.InsertCustomRow("Part & Yield & Par0 & Par2  & Par1 & $\\Chi^2$/NDF & Min X & Frac Above & \\langle p_{t} \\rangle  & \\langle p_{t}^{2} \\rangle");
+
+  }
+
+  AliLatexTable tempTable(4,"c|ccc");
+  tempTable.InsertCustomRow("Part & Yield & Yield Below & Frac Above\\\\");
+
+
+  //  Fit all  
+  for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+    
+    TCanvas * c2 = new TCanvas(TString("cCombined")+chargeFlag[icharge]+"_"+funcName[fitFuncID], TString("cCombined")+chargeFlag[icharge]);
+    TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,2.9, 100, 0.001,5);
+    hempty->SetXTitle("p_{t} (GeV/c)");
+    hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
+    hempty->Draw();
+    //    DrawStar(icharge);
+    c2->SetLogy();
+    TLegend * l = new TLegend(0.516779, 0.729021, 0.89094 ,0.916084, chargeLabel[icharge]);
+    l->SetFillColor(kWhite);
+
+    for(Int_t ipart = 0; ipart < kNPart; ipart++){
+      Float_t fitmin = 0;
+      Float_t fitmax = 3;
+
+      // Get functions
+      TF1 * func = 0;
+      if(fitFuncID == kFitLevi)          {
+       if (ipart == kPion)
+         func = fm->GetLevi(mass[ipart], 0.12, 7, 1.5);
+       if (ipart == kKaon)
+         func = fm->GetLevi(mass[ipart], 0.17, 7, 0.17);
+       if (ipart == kProton)
+         func = fm->GetLevi(mass[ipart], 0.15, 8, 0.08);
+      }
+      else if(fitFuncID == kFitUA1)      func = fm->GetUA1(mass[ipart],0.2,1.25,1.25,0.2,1.5);
+      else if(fitFuncID == kFitPowerLaw) {
+       if (ipart == kPion)
+         func = fm->GetPowerLaw(1.0, 8.6, 7);
+       if (ipart == kKaon)
+         func = fm->GetPowerLaw(3.0, 12, 2.6);
+       if (ipart == kProton)
+         func = fm->GetPowerLaw(24, 72, 0.8);
+      }
+      else if(fitFuncID == kFitPhojet)   func = fm->GetHistoFunc(hSpectraMC[kPhojet]        [ipart][icharge]);
+      else if(fitFuncID == kFitAtlasCSC) func = fm->GetHistoFunc(hSpectraMC[kPyTuneAtlasCSC][ipart][icharge]);
+      else if(fitFuncID == kFitPerugia0) func = fm->GetHistoFunc(hSpectraMC[kPyTunePerugia0][ipart][icharge]);
+      else if(fitFuncID == kFitCMS6D6T)  func = fm->GetHistoFunc(hSpectraMC[kPyTuneCMS6D6T] [ipart][icharge]);
+      else {
+       cout << "Unknown fit ID " << fitFuncID << endl;
+       return; 
+      }
+
+      if(fitFuncID >= kFitPhojet){
+       fitmin = 0.0;
+       fitmax = 1.0;
+      }
+
+      if(!AliBWTools::Fit(hSpectra[kCombTOFTPC][ipart][icharge],func,fitmin,fitmax)) {
+       cout << " FIT ERROR " << endl;
+       return;      
+      }
+      hSpectra[kCombTOFTPC][ipart][icharge]->Draw("same");    
+      hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Draw("same");
+      ((TF1*)hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0))->SetRange(0,4);
+      l->AddEntry(hSpectra[kCombTOFTPC][ipart][icharge], 
+                 scaleKaons && ipart == kKaon ? 
+                 (TString(partLabel[ipart][icharge])+" #times 2").Data() 
+                 : partLabel[ipart][icharge]);
+//       TF1 * fClone = (TF1*) hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Clone();
+//       hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->Add(fClone);
+//       fClone->SetLineStyle(kDashed);
+//       fClone->SetRange(0,100);
+//       fClone->Draw("same");
+      
+      // populate table
+      //     Float_t yield  = func->Integral(0.45,1.05);
+      //     Float_t yieldE = func->IntegralError(0.45,1.05);
+      
+      Float_t yield  = func->Integral(0.,100);
+      //Float_t yieldE = func->IntegralError(0.,100);
+
+      Double_t yieldTools, yieldETools;
+      Double_t partialYields[3],partialYieldsErrors[3]; 
+      AliBWTools::GetYield(hSpectra[kCombTOFTPC][ipart][icharge], func, yieldTools, yieldETools, 
+                          0, 100, partialYields,partialYieldsErrors);
+      Double_t tslope   = func->GetParameter(2);
+      Double_t tslopeE  = func->GetParError(2);        
+      
+      table.SetNextCol(partLabel[ipart][icharge]);
+      //table.SetNextCol(yield,yieldE,-4);
+      table.SetNextCol(yieldTools, yieldETools,-4);
+      table.SetNextCol(func->GetParameter(0));
+      table.SetNextCol(tslope,tslopeE,-4);
+      table.SetNextCol(func->GetParameter(1)); 
+      table.SetNextCol(Form("%2.2f/%d",func->GetChisquare(),func->GetNDF())); 
+      Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kCombTOFTPC][ipart][icharge]);
+      //Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kITS][ipart][icharge]);
+      Float_t yieldAbove  = func->Integral(lowestPoint,100);
+      table.SetNextCol(lowestPoint,-2);
+      table.SetNextCol(yieldAbove/yield,-2);
+      Float_t mean, meane;
+      Float_t mean2, mean2e;
+      AliBWTools::GetMean      (func, mean,  meane );
+      AliBWTools::GetMeanSquare(func, mean2, mean2e);
+      table.SetNextCol(mean,  meane ,-4);
+      table.SetNextCol(mean2, mean2e,-4);
+      
+      //                        fMean2->IntegralError(0,100)/func->Integral(0,100),-7);
+      table.InsertRow();
+      
+
+      /// TEMP TABLE
+      tempTable.SetNextCol(partLabel[ipart][icharge]);
+      tempTable.SetNextCol(yieldTools, yieldETools, -4);
+      tempTable.SetNextCol(partialYields[1], partialYieldsErrors[1], -4);
+      tempTable.SetNextCol(yieldAbove/yield,-2);
+      tempTable.InsertRow();
+    }
+    l->Draw();
+    if (doPrint) {
+      c2->Update();
+      gSystem->ProcessEvents();
+      c2->Print(TString(c2->GetName()) + ".eps");
+      ALICEWorkInProgress(c2,"","#splitline{ALICE Preliminary}{Statistical Error Only}");
+      c2->Print(TString(c2->GetName()) + ".png");
+    }
+    
+  }
+
+  
+  table.PrintTable("ASCII");
+  
+  cout << "" << endl;
+  tempTable.PrintTable("ASCII");
+
+}
+
+void LoadSpectra() {
+
+  TFile * f=0;
+
+  // Load
+
+
+  // TOF
+  //  f = new TFile("./Files/spectra-pos-y.root");
+  f = new TFile("./Files/spectra-pos-y_20100615.root");
+  hSpectra[kTOF][kPion]  [kPos]= (TH1F*) f->Get("hpi");
+  hSpectra[kTOF][kProton][kPos]= (TH1F*) f->Get("hpr");
+  hSpectra[kTOF][kKaon]  [kPos]= (TH1F*) f->Get("hka");
+  hSpectra[kTOF][kPion]  [kPos]->SetName("hpiPos");
+  hSpectra[kTOF][kProton][kPos]->SetName("hprPos");
+  hSpectra[kTOF][kKaon]  [kPos]->SetName("hkaPos");
+  //f = new TFile("./Files/spectra-neg-y.root");
+  f = new TFile("./Files/spectra-neg-y_20100615.root");
+  hSpectra[kTOF][kPion]  [kNeg]= (TH1F*) f->Get("hpi");
+  hSpectra[kTOF][kProton][kNeg]= (TH1F*) f->Get("hpr");
+  hSpectra[kTOF][kKaon]  [kNeg]= (TH1F*) f->Get("hka");
+  hSpectra[kTOF][kPion]  [kNeg]->SetName("hpiNeg");
+  hSpectra[kTOF][kProton][kNeg]->SetName("hprNeg");
+  hSpectra[kTOF][kKaon]  [kNeg]->SetName("hkaNeg");
+
+
+  // Clean UP TPC spectra, removing unwanted points
+  cout << "Cleaning Up TOF spectra" << endl;
+  Int_t nbin =  hSpectra[kTOF][kKaon][kPos]->GetNbinsX();
+  for(Int_t ibin = 1; ibin <= nbin; ibin++){
+    Float_t pt =  hSpectra[kTOF][kKaon][kPos]->GetBinCenter(ibin);
+    for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+      if (pt > 2.35) {
+       hSpectra[kTOF][kKaon][icharge]->SetBinContent(ibin,0);
+       hSpectra[kTOF][kKaon][icharge]->SetBinError  (ibin,0);  
+       hSpectra[kTOF][kProton][icharge]->SetBinContent(ibin,0);
+       hSpectra[kTOF][kProton][icharge]->SetBinError  (ibin,0);        
+      }            
+    }
+  }
+//   cout << "Scaling TOF to TPC" << endl;  
+//   // Scale protons to TPC
+//   hSpectra[kTOF][kProton][kPos]->Scale(1./1.05);
+//   // Scale pbar so that pbar/p is compatible with Panos
+//   hSpectra[kTOF][kProton][kNeg]->Scale(1./1.1);
+  
+  
+  
+
+
+  // ITS SA (Emanuele)
+  //  f = new TFile("./Files/ITSsaSPECTRA_3clusters20100619.root");
+  f = new TFile("./Files/ITSsaSPECTRAperMICHELE_20100703.root");
+  hSpectra[kITS][kPion]  [kPos]= (TH1F*) f->Get("hSpectraPos0");
+  hSpectra[kITS][kKaon]  [kPos]= (TH1F*) f->Get("hSpectraPos1");
+  hSpectra[kITS][kProton][kPos]= (TH1F*) f->Get("hSpectraPos2");
+  hSpectra[kITS][kPion]  [kNeg]= (TH1F*) f->Get("hSpectraNeg0");
+  hSpectra[kITS][kKaon]  [kNeg]= (TH1F*) f->Get("hSpectraNeg1");
+  hSpectra[kITS][kProton][kNeg]= (TH1F*) f->Get("hSpectraNeg2");
+
+  for(Int_t ipart = 0; ipart < kNPart; ipart++){
+    for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+      Int_t nbin = hSpectra[kITS][ipart][icharge]->GetNbinsX();
+      //      hSpectra[kITS][ipart][icharge]->Scale(276004.);// Emanule divided his spectra...
+      for(Int_t ibin = 1; ibin <= nbin; ibin++){
+       if(hSpectra[kITS][ipart][icharge]->GetBinContent(ibin) < 0 ) {
+         hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
+         hSpectra[kITS][ipart][icharge]->SetBinError  (ibin,0);
+       }
+//     if ((ipart == kKaon && ibin >= 12) || (ipart == kProton && ibin >= 20)) {
+//       hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
+//       hSpectra[kITS][ipart][icharge]->SetBinError  (ibin,0);
+//     }
+      }
+      
+    }
+  }
+
+
+  // ITS + TPC (Marek)
+  f = TFile::Open("./Files/SpectraCorrectedITSBeforeProtons_20100618.root");
+  TList * list = (TList*) gDirectory->Get("output");
+  hSpectra[kITSTPC][kPion]  [kPos]= (TH1F*) list->FindObject("Pions");
+  hSpectra[kITSTPC][kKaon]  [kPos]= (TH1F*) list->FindObject("Kaons");
+  hSpectra[kITSTPC][kProton][kPos]= (TH1F*) list->FindObject("Protons");
+  hSpectra[kITSTPC][kPion]  [kNeg]= (TH1F*) list->FindObject("AntiPions");
+  hSpectra[kITSTPC][kKaon]  [kNeg]= (TH1F*) list->FindObject("AntiKaons");
+  hSpectra[kITSTPC][kProton][kNeg]= (TH1F*) list->FindObject("AntiProtons");
+
+  // TPC
+  //  htemplate =  hSpectra[kITS][kProton][kNeg]; //FIXME
+  f = new TFile("./Files/protonSpectra_20100615.root");
+  hSpectra[kTPC][kProton][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonPosClassic"),htemplate);
+  hSpectra[kTPC][kProton][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonNegClassic"),htemplate);
+  f = new TFile("./Files/pionSpectra_20100615.root");
+  hSpectra[kTPC][kPion][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionPosClassic"),htemplate);
+  hSpectra[kTPC][kPion][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionNegClassic"),htemplate);
+  f = new TFile("./Files/kaonSpectra_20100615.root");
+  hSpectra[kTPC][kKaon][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonPosClassic"),htemplate);
+  hSpectra[kTPC][kKaon][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonNegClassic"),htemplate);
+
+  // Clean UP TPC spectra, removing unwanted points
+  cout << "Cleaning Up TPC spectra" << endl;
+  nbin =  hSpectra[kTPC][kKaon][kPos]->GetNbinsX();
+  for(Int_t ibin = 0; ibin < nbin; ibin++){
+    Float_t pt =  hSpectra[kTPC][kKaon][kPos]->GetBinCenter(ibin);
+    if (pt > 0.45) {
+      for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+       hSpectra[kTPC][kKaon][icharge]->SetBinContent(ibin,0);
+       hSpectra[kTPC][kKaon][icharge]->SetBinError  (ibin,0);  
+      }      
+    }
+    if (pt < 0.45) {
+      for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+       hSpectra[kTPC][kProton][icharge]->SetBinContent(ibin,0);
+       hSpectra[kTPC][kProton][icharge]->SetBinError  (ibin,0);        
+      }      
+    }
+  }
+  
+
+
+  
+  // K0s
+  f = new TFile ("./Files/PtSpectraCorrectedK0sOff-dNdy-Jun02.root");
+  //  hSpectra[kK0][kKaon][kPos] = (TH1F*) AliBWTools::GetdNdPtFromOneOverPt((TH1*) gDirectory->Get("hSpectraOff")); 
+  hSpectra[kK0][kKaon][kPos] = (TH1F*) gDirectory->Get("hSpectraOff"); 
+  //  hSpectra[kK0][kKaon][kPos]->Scale(2*TMath::Pi());
+  //  hSpectra[kK0][kKaon][kPos]->Scale(1./272463);
+  hSpectra[kK0][kKaon][kNeg] = hSpectra[kK0][kKaon][kPos];
+
+  // Kinks: TO BE FIXED WITH POSITIVES AND NEGATIVES
+  //  f = new TFile ("./Files/PtAllKaonKinkRap6Apr24.root");
+  f = new TFile ("./Files/PtKaonKinkJune13AllPN_20100615.root");
+  hSpectra[kKinks][kKaon][kPos] = (TH1F*)gDirectory->Get("fptallKPA");
+  hSpectra[kKinks][kKaon][kNeg] = (TH1F*)gDirectory->Get("fptallKNA");
+
+
+  // Apply correction factrs
+  // Secondaries for protons
+  f = new TFile ("./Files/corrFactorProtons_20100615.root");
+  TH1F * hCorrSecondaries = (TH1F*) gDirectory->Get("corrFactorProtons");
+  if(correctSecondaries) {
+    cout << "CORRECTING SECONDARIES" << endl;
+    
+    for(Int_t idet = 0; idet <= kTOF; idet++){ // TPC and TOF only
+      for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+       Int_t ipart = kProton;
+       TH1* h = hSpectra[idet][ipart][icharge]; // lighter notation below
+       if (h){
+         Int_t nbins = h->GetNbinsX();
+         for(Int_t ibin = 0; ibin < nbins; ibin++){
+           //      Float_t pt = convertToMT ? TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) : h->GetBinCenter(ibin);
+           Float_t pt = h->GetBinCenter(ibin);
+           if (icharge == kNeg) pt = -pt;
+           Int_t binCorrection = hCorrSecondaries->FindBin(pt);
+           Float_t correction    = hCorrSecondaries->GetBinContent(binCorrection);
+           //      cout << pt << " " << correction << endl;
+           
+           if (correction != 0) {// If the bin is empty this is a  0
+             h->SetBinContent(ibin,h->GetBinContent(ibin)/correction);
+             h->SetBinError  (ibin,h->GetBinError  (ibin)/correction);
+           } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
+             cout << "Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
+             cout << " Bin content: " << h->GetBinContent(ibin)  << endl;
+             
+           }
+         }     
+       }
+      }
+    }
+  }
+  // geant/fluka absorption
+  if(correctGeantFlukaXS) {
+    cout << "CORRECTING GEANT3/FLUKA" << endl;
+    
+    f = new TFile ("./Files/correctionForCrossSection_20100615.root");
+    TH2D * hCorrFluka[kNCharge];
+    hCorrFluka[kPos] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionProtons");
+    hCorrFluka[kNeg] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionAntiProtons");
+    for(Int_t idet = 0; idet < kNDet; idet++){
+      if (idet == kITS) continue; // skip ITS sa
+      for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+       Int_t ipart = kProton;
+       TH1 * h = hSpectra[idet][ipart][icharge]; // lighter notation below
+       if (h){
+         Int_t nbins = h->GetNbinsX();
+         for(Int_t ibin = 0; ibin < nbins; ibin++){
+//         Float_t pt = convertToMT ? 
+//           TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) :
+//           h->GetBinCenter(ibin);
+           Float_t pt = h->GetBinCenter(ibin);
+           Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
+           Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(h->GetNbinsY()+1);
+           if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
+           if (pt > maxPtCorrection) pt = maxPtCorrection;
+           Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
+           if (idet == kTOF && icharge == kNeg) {
+
+             // Apply parametrized correction computed by francesco
+             // Fitted panos correction and using momentum at the outer radius of the TPC
+
+             Float_t ptav = pt; // Just to use the same name francesco uses...
+
+             // from pT constrained at P.V. (ptav) to pT TPC outer (ptTPCout)        
+             Float_t ptTPCout=ptav*(1-6.81059e-01*TMath::Exp(-ptav*4.20094));
+             
+             // traking correction (fit to Panos)
+             Float_t antiprotonEC = 1 - 0.129758 *TMath::Exp(-ptav*0.679612);
+             
+             // TOF matching efficiency correction (derived from Panos one scaled for M.B.(TOF)/M.B.(TPC)).
+             Float_t antiprotonEC2 = TMath::Power(1 -  0.129758*TMath::Exp(-ptTPCout*0.679612),0.07162/0.03471);
+             correction = antiprotonEC * antiprotonEC2;
+           }
+           //      cout << icharge<< " " << h->GetBinCenter(ibin) << " " << pt << " " << correction << endl;
+           if (correction != 0) {// If the bin is empty this is a  0
+             h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
+             h->SetBinError  (ibin,h->GetBinError  (ibin)*correction);
+//           if (idet == kTOF) {
+//             cout << "CORRECTING TOF TWICE" << endl;
+//             h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
+//             h->SetBinError  (ibin,h->GetBinError  (ibin)*correction);               
+//           }
+           } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
+             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
+             cout << " Bin content: " << h->GetBinContent(ibin)  << endl;
+           }
+           
+         }
+         
+       }
+      }
+    }    
+  }
+
+
+  // Set style and scale
+  for(Int_t idet = 0; idet < kNDet; idet++){
+    for(Int_t ipart = 0; ipart < kNPart; ipart++){
+      for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+       if (hSpectra[idet][ipart][icharge]){
+         hSpectra[idet][ipart][icharge]->SetStats(0); // disable stats
+         hSpectra[idet][ipart][icharge]->SetMarkerColor (color[ipart] );
+         hSpectra[idet][ipart][icharge]->SetLineColor   (color[ipart] );
+         hSpectra[idet][ipart][icharge]->SetMarkerStyle (marker[ipart]);
+         hSpectra[idet][ipart][icharge]->SetXTitle("p_{t} (GeV/c)");
+         hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} dN/dp_{t} (GeV/c)^{-1}");
+         if (convertToMT) {
+           TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[idet][ipart][icharge]);
+           hSpectra[idet][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
+           hSpectra[idet][ipart][icharge]->SetXTitle("m_{t} (GeV)");
+           hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
+         }
+//       if (idet == kTOF || idet == kTPC) {
+//           hSpectra[idet][ipart][icharge]->Scale(1./236763);       
+//       } 
+//       if (idet == kITS){
+//         hSpectra[idet][ipart][icharge]->Scale(202945./236763);                    
+//       }
+         if (scaleKaons && ipart == kKaon) {
+           hSpectra[idet][ipart][icharge]->Scale(2.);              
+         }
+       } else {
+         printf("Cannot load %s,%s,%s\n",detFlag[idet],partFlag[ipart],chargeFlag[icharge]);
+       }
+      }
+    }
+  }
+
+  
+
+
+  // Combine detectors
+  for(Int_t ipart = 0; ipart < kNPart; ipart++){
+    for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+      TH1F * htemplLocal = htemplate; // If we are converting to 1/mt dNdmt we need to convert the template as well...
+      
+      if (convertToMT) {
+       TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(htemplate);
+       htemplLocal = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
+
+      }
+      hSpectra[kCombTOFTPC][ipart][icharge] = AliBWTools::CombineHistos(hSpectra[kTOF][ipart][icharge],
+                                                                       hSpectra[kTPC][ipart][icharge],
+                                                                       htemplLocal,1.);;
+//       if (convertToMT) {
+//     TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[kCombTOFTPC][ipart][icharge]);
+//     hSpectra[kCombTOFTPC][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
+//     hSpectra[kCombTOFTPC][ipart][icharge]->SetXTitle("m_{t} (GeV)");
+//     hSpectra[kCombTOFTPC][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
+//       }
+    }
+  }
+
+
+  // Scale for the number of inelastic collisions and correct for
+  // efficiency losses due to physics selection:
+
+  Double_t effPhysSel[kNPart];
+  effPhysSel[kPion]   = 1.012;
+  effPhysSel[kKaon]   = 1.013;
+  effPhysSel[kProton] = 1.014;
+
+
+  for(Int_t idet = 0; idet < kNHist; idet++){
+    for(Int_t ipart = 0; ipart < kNPart; ipart++){
+      for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+       if(hSpectra[idet][ipart][icharge]) {
+         //      cout << "Scaling!" << endl;
+         hSpectra[idet][ipart][icharge]->Scale(1.*effPhysSel[ipart]/278366.15); // Scale PhysSel tutti? // FIXME
+       }
+      }
+    }
+  }
+
+
+}
+
+void LoadMC() {
+
+  TFile * f = 0;
+  const char * evClass= "INEL";
+  const char * files[] = {"./Files/dndeta_Phojet_10M_900GeV.root", 
+                         "./Files/dndeta_AtlasCSC306_10M_900GeV.root", 
+                         "./Files/dndeta_CMSD6T109_10M_900GeV.root", 
+                         "./Files/dndeta_Perugia0320_10M_900GeV.root", };
+  
+  // Phojet
+  for(Int_t itune = 0; itune < kNTunes; itune++){
+    f = new TFile(files[itune]);
+    for(Int_t ipart = 0; ipart < kNPart; ipart++){
+      for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+       hSpectraMC[itune][ipart][icharge] = (TH1F*) f->Get(Form("fHistPtID_%s_%s%s",evClass,partFlag[ipart],icharge==0 ? "Pos" : "Neg"));
+      }
+    }
+  }
+  
+
+  // Set style
+  for(Int_t itune = 0; itune < kNTunes; itune++){
+    for(Int_t ipart = 0; ipart < kNPart; ipart++){
+      for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+       if (hSpectraMC[itune][ipart][icharge]){
+         hSpectraMC[itune][ipart][icharge]->SetName(TString(hSpectraMC[itune][ipart][icharge]->GetName())+mcTuneName[itune]);
+         hSpectraMC[itune][ipart][icharge]->SetMarkerColor (mcLineColor[itune] );
+         hSpectraMC[itune][ipart][icharge]->SetLineColor   (mcLineColor[itune] );
+         hSpectraMC[itune][ipart][icharge]->SetLineStyle   (mcLineStyle[itune] );
+         hSpectraMC[itune][ipart][icharge]->SetMarkerStyle (1);
+         if (convertToMT) {
+           TH1F * htmp = (TH1F*)AliBWTools::GetOneOverPtdNdPt(hSpectraMC[itune][ipart][icharge]);
+           hSpectraMC[itune][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
+           hSpectraMC[itune][ipart][icharge]->SetXTitle("m_{t} (GeV)");
+           hSpectraMC[itune][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
+         }
+
+       } else {
+         printf("Cannot load MC # %d,%s,%s\n",itune,partFlag[ipart],chargeFlag[icharge]);
+       }
+      }
+    }
+  }
+  
+}
+
+void DrawStar(Int_t icharge) {
+  //  cout << icharge << endl;
+  
+  //  gROOT->LoadMacro("StarPPSpectra.C");
+  TGraphErrors ** gStar = StarPPSpectra();
+  
+  for(Int_t istar = 0; istar < 6; istar++){
+    gStar[istar]->SetMarkerStyle(kOpenStar);
+    if      (icharge==kPos &&  (istar%2) ) gStar[istar]->Draw("P");
+    else if (icharge==kNeg && !(istar%2) ) gStar[istar]->Draw("P");
+    else cout << "Skipping Star " << istar << endl;    
+  }
+}
+
+void GetITSResiduals() {
+
+  for(Int_t ipart = 0; ipart < kNPart; ipart++){
+    for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+      //      cout << "1 " << ipart << " " << icharge << endl;
+      
+//       gSystem->ProcessEvents();
+//       gSystem->Sleep(1000);
+      TF1* f = (TF1*)   hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0);
+      TH1F * hres1, *hres2;
+      AliBWTools::GetResiduals(hSpectra[kITS][ipart][icharge], f, &hres1, &hres2);
+
+      TCanvas * c1 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge],TString(partFlag[ipart])+"_"+chargeFlag[icharge]);
+      c1->SetLogy();
+      hSpectra[kCombTOFTPC][ipart][icharge]->Draw();
+      hSpectra[kITS][ipart][icharge]->SetMarkerStyle(24);
+      hSpectra[kCombTOFTPC][ipart][icharge]->SetMarkerStyle(20);
+      hSpectra[kITS][ipart][icharge]->Draw("same");
+      hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Draw("same");
+      TLegend* l = new TLegend(    0.182886,    0.192308,    0.505034,0.384615, TString(partLabel[ipart][icharge])+" "+chargeFlag[icharge]);
+      l->AddEntry(hSpectra[kCombTOFTPC][ipart][icharge], "TOF+TPC");
+      l->AddEntry(hSpectra[kITS][ipart][icharge],        "ITS");
+      l->AddEntry(f,        "Fit to TOF+TPC");
+      l->Draw();
+
+      TCanvas * c2 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res",
+                                TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res");  
+      c2->SetGridy();
+      hres2->SetMinimum(-1);
+      hres2->SetMaximum(1);
+      hres2->Draw();
+      hres2->GetYaxis()->SetTitleOffset(1.2);
+      Float_t x = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kCombTOFTPC][ipart][icharge]);
+      TLine * line = new TLine(x,-1,x,1);
+      line->SetLineStyle(kDashed);
+      line->Draw("same");
+      
+      if (doPrint) {
+       c1->Update();
+       c2->Update();
+       gSystem->ProcessEvents();
+       c1->Print(TString(c1->GetName()) + ".png");
+       c2->Print(TString(c2->GetName()) + ".png");
+      }
+    }
+  }
+}
+
+void DrawWithModels() {
+
+  for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+    // Draw with models
+    for(Int_t ipart = 0; ipart < kNPart; ipart++){
+      // Multipad canvas
+      TCanvas * c1 = new TCanvas(TString("cSpectra")+partFlag[ipart]+chargeFlag[icharge],TString("cSpectra")+partFlag[ipart]+chargeFlag[icharge] );
+      TPad *p1 = new TPad(TString("p1")+partFlag[ipart]+chargeFlag[icharge], "p1", 0.0, 0.3, 1.0,  0.95, 0, 0, 0);
+      p1->SetBottomMargin(0);
+      p1->Draw();
+      
+      TPad *p2 = new TPad(TString("p2")+partFlag[ipart]+chargeFlag[icharge], "p2", 0.0, 0.05, 1.0,  0.3, 0, 0, 0);
+      p2->SetTopMargin(0);
+      p2->SetBottomMargin(0.4);
+      p2->Draw();
+
+      Float_t scaleFonts = (0.95-0.3)/(0.3-0.05);
+
+      // Draw spectra
+      p1->cd();
+      p1->SetLogy();
+      TH2F * hempty = new TH2F(TString("hempty")+long(ipart),"hempty",100,0.,4, 100, 0.0015,5);
+      hempty->SetXTitle("p_{t} (GeV/c)");
+      hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
+      hempty->Draw();
+      c1->SetLogy();
+
+
+      TLegend * l =new TLegend(       0.543624,  0.431818,  0.892617,0.629371);
+      l->SetFillColor(kWhite);
+      hSpectra[kCombTOFTPC][ipart][icharge]->Draw("same");
+      l->AddEntry(hSpectra[kTOF][ipart][icharge],TString ("Data: ")+partLabel[ipart][icharge]);
+      for(Int_t itune = 0; itune < kNTunes; itune++){      
+       l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
+       hSpectraMC[itune][ipart][icharge]->SetLineWidth(2);    
+       hSpectraMC[itune][ipart][icharge]->Draw("same,chist");    
+      }
+      l->Draw("same");
+
+      // Draw ratio
+      p2->cd();
+      TH2F * hemptyr = new TH2F(TString("hemptyratio")+long(ipart),"hempty",100,0.,4, 100, 0.01,2.99);
+      hemptyr->SetXTitle("p_{t} (GeV/c)");
+      hemptyr->SetYTitle("Data/MC");
+      hemptyr->GetXaxis()->SetLabelSize(0.04*scaleFonts);
+      hemptyr->GetYaxis()->SetLabelSize(0.04*scaleFonts);
+      hemptyr->GetYaxis()->SetTitleSize(0.05*scaleFonts);
+      hemptyr->GetYaxis()->SetTitleOffset(1.4/scaleFonts);
+      hemptyr->GetXaxis()->SetTitleSize(0.05*scaleFonts);
+      hemptyr->SetTickLength(0.03*scaleFonts, "X");
+      hemptyr->SetTickLength(0.02*scaleFonts, "Y");
+      //      hemptyr->GetXaxis()->SetTitleOffset(1.4/scaleFonts);
+      hemptyr->GetYaxis()->SetNdivisions(5);
+      hemptyr->Draw("");
+
+      AliBWFunc fm;
+      for(Int_t itune = 0; itune < kNTunes; itune++){      
+       TF1* f = fm.GetHistoFunc(hSpectraMC[itune][ipart][icharge], TString("f")+mcTuneName[itune]);
+
+       //      l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
+       TH1F* hRatio = AliBWTools::DivideHistoByFunc(hSpectra[kCombTOFTPC][ipart][icharge],f);
+       hRatio->SetLineStyle(hSpectraMC[itune][ipart][icharge]->GetLineStyle());
+       hRatio->SetLineColor(hSpectraMC[itune][ipart][icharge]->GetLineColor());
+       hRatio->SetLineWidth(hSpectraMC[itune][ipart][icharge]->GetLineWidth());
+       hRatio->Draw("lhist,same");
+      }
+
+
+      // print
+      if(doPrint) {
+       c1->Update();
+       gSystem->ProcessEvents();
+       c1->Print(TString(c1->GetName())+".eps");
+       ALICEWorkInProgress(c1,"","#splitline{ALICE Preliminary}{Statistical Error Only}");
+       c1->Print(TString(c1->GetName())+".png");
+       
+      }
+    }
+  }
+
+
+
+}
+
+void DrawAllAndKaons() {
+
+
+  //  gROOT->LoadMacro("ALICEWorkInProgress.C");
+
+  gStyle->SetOptFit(0);
+
+  TH1F * hKaonsAllTPCTOF = (TH1F*) hSpectra[kCombTOFTPC][kKaon][kPos]->Clone();
+  hKaonsAllTPCTOF->Add(hSpectra[kCombTOFTPC][kKaon][kNeg]);
+  
+  TH1F * hK0Scaled    = (TH1F*) hSpectra[kK0][kKaon][kPos]->Clone();
+  hK0Scaled->Add(hSpectra[kK0][kKaon][kPos]);
+
+  hSpectra[kKinks][kKaon][kPos]->SetMarkerStyle(25);
+  hSpectra[kKinks][kKaon][kPos]->SetStats(0);
+  TH1F * hKinksAll = (TH1F*) hSpectra[kKinks][kKaon][kPos]->Clone();
+  hKinksAll->Add(hSpectra[kKinks][kKaon][kNeg]);
+  
+  TCanvas * c1 = new TCanvas("cKaons","cKaons");
+  c1->SetLogy();
+  TH2F * hempty = new TH2F("hempty_allkaons","hempty",100,0.,3, 100, 1e-4,10);
+  hempty->SetXTitle("p_{t} (GeV/c)");
+  hempty->SetYTitle("dN / dp_{t} (A.U.)");
+  hempty->Draw();
+  hKaonsAllTPCTOF->Draw("same");
+  hK0Scaled->Draw("same");
+  hKinksAll->Draw("same");
+
+  TLegend * leg = new TLegend(0.2013423,0.2255245,0.5503356,0.4335664,NULL,"brNDC");
+//    leg->SetBorderSize(0);
+//    leg->SetLineColor(1);
+//    leg->SetLineStyle(1);
+//    leg->SetLineWidth(1);
+//    leg->SetFillColor(19);
+//    leg->SetFillStyle(1001);
+   TLegendEntry *entry=leg->AddEntry("hkaPos_h_tpcKaonPos","K^{+} + K^{-}, TPC+TOF ","lpf");
+   entry=leg->AddEntry("h1PtSpectraOff_inv","K^{0} #times 2","lpf");
+   entry=leg->AddEntry("fptallK","K^{+} + K ^{-}, Kinks","lpf");
+   leg->Draw();
+
+   ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Performance}{Not fully corrected}");
+   TLatex * tex = new TLatex(0.2120805,0.01288336,"Statistical error only");
+   tex->SetTextColor(2);
+   tex->SetTextFont(42);
+   tex->SetTextSize(0.03496503);
+   tex->Draw();
+
+   c1->Update();
+   if(doPrint) c1->Print(TString(c1->GetName())+".png");
+
+  // Draw all "stable" hadrons
+  for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+    TCanvas * c1 = new TCanvas(TString("cAll_")+chargeFlag[icharge],TString("cAll_")+chargeFlag[icharge]);
+    c1->SetLogy();
+    TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,4, 100, 1e-4,10);
+    hempty->SetXTitle("p_{t} (GeV/c)");
+    hempty->SetYTitle("dN / dp_{t} (A.U.)");
+    hempty->Draw();
+    leg = new TLegend(  0.645973,  0.325175,  0.892617,0.636364, NULL,"brNDC");
+    for(Int_t ipart = 0; ipart < kNPart; ipart++) {
+      for(Int_t idet = 0; idet <= kITSTPC; idet++){
+       //      if (idet == kITS) continue;
+       //      if (idet == kITSTPC) hSpectra[idet][ipart][icharge]->SetMarkerColor(kGreen);
+       hSpectra[idet][ipart][icharge]->SetMarkerStyle(marker[idet]);
+       hSpectra[idet][ipart][icharge]->Draw("same");
+       leg->AddEntry(hSpectra[idet][ipart][icharge],TString(partLabel[ipart][icharge])+" (" + detFlag[idet]  + ")","lpf");
+      }
+      //      leg->AddLine();
+    }    
+    leg->Draw();
+    ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Performance}{Not Fully Corrected}");
+    c1->Update();
+    if(doPrint) c1->Print(TString(c1->GetName())+".png");
+  }
+
+  //  Draw ratios (tmp)
+
+//   new TCanvas;
+//   hSpectra[kTPC][kKaon][kNeg]->Divide(hSpectra[kTPC][kKaon][kPos]);
+//   hSpectra[kTPC][kKaon][kNeg]->Draw();
+//   new TCanvas;
+//   hSpectra[kITS][kKaon][kNeg]->Divide(hSpectra[kITS][kKaon][kPos]);
+//   hSpectra[kITS][kKaon][kNeg]->Draw("");
+  //    hSpectra[kTOF][kProton][kPos]->Draw("same");
+  
+//   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+//     TCanvas * c1 = new TCanvas(TString("cAllRatio_")+chargeFlag[icharge],TString("cAllRatio_")+chargeFlag[icharge]);
+//     c1->SetGridy();
+//     TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.2,1, 100, 0.0,2.0);
+//     hempty->SetXTitle("p_{t} (GeV/c)");
+//     hempty->SetYTitle("ITSsa / TPC");
+//     hempty->Draw();
+    
+
+//     for(Int_t ipart = 0; ipart < kNPart; ipart++) {
+//       hSpectra[kITS][ipart][icharge]->Divide(hSpectra[kTPC][ipart][icharge]);
+//       hSpectra[kITS][ipart][icharge]->Draw("same");
+//     }
+//     if(doPrint) c1->Print(TString(c1->GetName())+".png");
+//   }
+    
+// end of tmp
+
+}
+
+void DrawWithJacek() {
+
+  //1. Convert spectra to dNdeta and sum
+  TH1F * hsum = (TH1F*) htemplate->Clone();
+  hsum->Reset();
+  Int_t idet= kCombTOFTPC;
+  for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+    for(Int_t ipart = 0; ipart < kNPart; ipart++){
+      TH1 * h = hSpectra[idet][ipart][icharge];
+      Int_t nbin = h->GetNbinsX();
+      for(Int_t ibin = 1; ibin <= nbin; ibin++){
+       Double_t pt = h->GetBinCenter(ibin);
+       Double_t mt = TMath::Sqrt(pt*pt + mass[ipart]*mass[ipart]);
+       Double_t jacobian = pt/mt;
+       h->SetBinContent(ibin,h->GetBinContent(ibin)*jacobian);
+       h->SetBinError  (ibin,h->GetBinError  (ibin)*jacobian);
+       Int_t ibinSum = hsum->FindBin(pt);
+       Double_t epsilon = 0.0001;
+       if ( h->GetBinContent(ibin) > 0 && 
+            (TMath::Abs(h->GetBinLowEdge(ibin)   - hsum->GetBinLowEdge(ibinSum)) > epsilon || 
+             TMath::Abs(h->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
+            ) {
+         cout << "DISCREPANCY IN BIN RANGES" << endl;
+         cout << pt << " " << ibinSum << " " << ibin  << "; " << h->GetBinContent(ibin) << endl
+              << h->GetBinLowEdge(ibin) << "-"  << h->GetBinLowEdge(ibin+1) << endl
+              << hsum->GetBinLowEdge(ibinSum) << "-"  << hsum->GetBinLowEdge(ibinSum+1) << endl;
+         cout << "" << endl;       
+       }
+       
+
+       hsum->SetBinContent(ibinSum,hsum->GetBinContent(ibinSum)+h->GetBinContent(ibin)); // EROOR FIXME
+       hsum->SetBinError  (ibinSum,0);
+      }
+      //       hsum->Add(h);
+    }      
+  }    
+
+
+  // Load Jacek and Draw both:  
+  new TFile ("./Files/dNdPt_Data_Points_ALICE_900GeV.root");
+  TGraphErrors * gJacek = (TGraphErrors*) gDirectory->Get("inel");
+  gJacek->Draw("AP");
+  hsum->Draw("same");
+
+  TGraphErrors * gRatio = AliBWTools::DivideGraphByHisto(gJacek,hsum);
+  
+  new TCanvas();
+  gRatio->Draw("AP");
+
+  
+
+}
+
+
+void DrawRatioToStar() {
+
+  // Star data
+  //  gROOT->LoadMacro("StarPPSpectra.C");
+  TGraphErrors ** gStar = StarPPSpectra();
+
+  // ALICE, INEL -> NSD
+  Double_t scaleYield = 3.58/3.02; // from paper 2
+  for(Int_t ipart = 0; ipart < kNPart; ipart++){
+    for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+      hSpectra[kCombTOFTPC][ipart][icharge]->Scale(scaleYield);
+    }    
+  }
+  
+    
+  TCanvas * c1 = new TCanvas("cRatioToStarNeg","cRatioToStarNeg");
+  TH2F * hempty = new TH2F(TString("hemptyNeg"),"hemptyNeg",100,0.,1.5, 100, 0.001,1.8);
+  hempty->SetXTitle("p_{t} (GeV/c)");
+  hempty->SetYTitle("ALICE/STAR (NSD)");
+  
+  TLegend *leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Negative","brNDC");
+  leg->SetBorderSize(0);
+  leg->SetLineColor(1);
+  leg->SetLineStyle(1);
+  leg->SetLineWidth(1);
+  leg->SetFillColor(0);
+  leg->SetFillStyle(1001);
+
+
+  hempty->Draw();
+  TGraphErrors * g ;
+  g = AliBWTools::DivideGraphByHisto(gStar[0],hSpectra[kCombTOFTPC][kPion][kNeg],1);
+  g->SetMarkerStyle(kFullCircle);
+  g->SetMarkerColor(kBlack);
+  g->Draw("p");
+  leg->AddEntry(g,"#pi^{-}","lp");
+  g = AliBWTools::DivideGraphByHisto(gStar[2],hSpectra[kCombTOFTPC][kKaon][kNeg],1);
+  g->SetMarkerStyle(kOpenCircle);
+  g->SetMarkerColor(kRed);
+  g->Draw("p");
+  leg->AddEntry(g,"K^{-}","lp");
+  g = AliBWTools::DivideGraphByHisto(gStar[4],hSpectra[kCombTOFTPC][kProton][kNeg],1);
+  g->SetMarkerStyle(kOpenSquare);
+  g->SetMarkerColor(kBlue);
+  g->Draw("p");
+  leg->AddEntry(g,"#bar{p}","lp");  
+  leg->Draw();
+
+
+
+
+  TCanvas * c2 = new TCanvas("cRatioToStarPos","cRatioToStarPos");
+  hempty->Draw(); 
+  leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Positive","brNDC");
+  leg->SetBorderSize(0);
+  leg->SetLineColor(1);
+  leg->SetLineStyle(1);
+  leg->SetLineWidth(1);
+  leg->SetFillColor(0);
+  leg->SetFillStyle(1001);
+ //  TGraphErrors * g ;
+  g = AliBWTools::DivideGraphByHisto(gStar[1],hSpectra[kCombTOFTPC][kPion][kPos],1);
+  g->SetMarkerStyle(kFullCircle);
+  g->SetMarkerColor(kBlack);
+  g->Draw("p");
+  leg->AddEntry(g,"#pi^{+}","lp");
+  g = AliBWTools::DivideGraphByHisto(gStar[3],hSpectra[kCombTOFTPC][kKaon][kPos],1);
+  g->SetMarkerStyle(kOpenCircle);
+  g->SetMarkerColor(kRed);
+  g->Draw("p");
+  leg->AddEntry(g,"K^{+}","lp");
+  g = AliBWTools::DivideGraphByHisto(gStar[5],hSpectra[kCombTOFTPC][kProton][kPos],1);
+  g->SetMarkerStyle(kOpenSquare);
+  g->SetMarkerColor(kBlue);
+  g->Draw("p");
+  leg->AddEntry(g,"p","lp");
+
+
+  leg->Draw();
+
+  c1->Update();
+  c2->Update();
+  gSystem->ProcessEvents();
+  c1->Print(TString(c1->GetName()) + ".eps");
+  c2->Print(TString(c2->GetName()) + ".eps");
+  ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
+  ALICEWorkInProgress(c2,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
+  c1->Update();
+  c2->Update();
+  c1->Print(TString(c1->GetName()) + ".png");
+  c2->Print(TString(c2->GetName()) + ".png");
+
+
+
+
+}
+
+
+
+void DrawRatios() {
+
+  // Draws ratios of combined spectra
+
+  // Compute ratios
+  TH1F * hPosNegRatio[kNPart];
+  
+  for(Int_t ipart = 0; ipart < kNPart; ipart++){
+    hPosNegRatio[ipart] = (TH1F*) hSpectra[kCombTOFTPC][ipart][kPos]->Clone();
+    hPosNegRatio[ipart]->Divide(hSpectra[kCombTOFTPC][ipart][kNeg]);
+    hPosNegRatio[ipart]->SetYTitle(TString(partLabel[ipart][kPos])+"/"+partLabel[ipart][kNeg]);    
+    hPosNegRatio[ipart]->SetMinimum(0.5);
+    hPosNegRatio[ipart]->SetMaximum(1.5);
+  }
+  
+  TH1F * hKPiRatio = (TH1F*) hSpectra[kCombTOFTPC][kKaon][kPos]->Clone();
+  hKPiRatio->Add(hSpectra[kCombTOFTPC][kKaon][kNeg]);
+  TH1F * htmp = (TH1F*) hSpectra[kCombTOFTPC][kPion][kPos]->Clone();
+  htmp->Add(hSpectra[kCombTOFTPC][kPion][kNeg]);
+  hKPiRatio->Divide(htmp);
+  hKPiRatio->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
+
+  TH1F * hKPiRatioMC[kNTunes];
+  for(Int_t itune = 0; itune < kNTunes; itune++){
+    hKPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kKaon][kPos]->Clone();
+    hKPiRatioMC[itune]->Add(hSpectraMC[itune][kKaon][kNeg]);
+    TH1F * htmp = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
+    htmp->Add(hSpectraMC[itune][kPion][kNeg]);
+    hKPiRatioMC[itune]->Divide(htmp);
+    hKPiRatioMC[itune]->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");    
+  }
+  
+
+
+  TH1F * hPPiRatio = (TH1F*) hSpectra[kCombTOFTPC][kProton][kPos]->Clone();
+  hPPiRatio->Add(hSpectra[kCombTOFTPC][kProton][kNeg]);
+  htmp = (TH1F*) hSpectra[kCombTOFTPC][kPion][kPos]->Clone();
+  htmp->Add(hSpectra[kCombTOFTPC][kPion][kNeg]);
+  hPPiRatio->Divide(htmp);
+  hPPiRatio->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
+
+  TH1F * hPPiRatioMC[kNTunes];
+  for(Int_t itune = 0; itune < kNTunes; itune++){
+    hPPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kProton][kPos]->Clone();
+    hPPiRatioMC[itune]->Add(hSpectraMC[itune][kProton][kNeg]);
+    TH1F * htmp = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
+    htmp->Add(hSpectraMC[itune][kPion][kNeg]);
+    hPPiRatioMC[itune]->Divide(htmp);
+    hPPiRatioMC[itune]->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");    
+  }
+  // Draw
+//   TH2F * hempty = new TH2F(TString("hempty"),"hempty",100,0.,1.5, 100, 0.001,1.8);
+//   hempty->SetXTitle("p_{t} (GeV/c)");
+  //  hempty->SetYTitle("");
+
+  // tmp: overlay levi fits
+  AliBWFunc * fm2 = new AliBWFunc;
+  fm2->SetVarType(AliBWFunc::kdNdpt);
+  TF1 * fLevi[kNPart]  [kNCharge];
+  fLevi[kPion]  [kPos] = fm2->GetLevi (mass[0], 0.1243, 7.614785, 1.524167, "fLeviPiPlus");
+  fLevi[kKaon]  [kPos] = fm2->GetLevi (mass[1], 0.1625, 5.869318, 0.186361, "fLeviKPlus");
+  fLevi[kProton][kPos] = fm2->GetLevi (mass[2], 0.1773, 6.918065, 0.086389, "fLeviPPlus");
+  fLevi[kPion]  [kNeg] = fm2->GetLevi (mass[0], 0.1267, 7.979582, 1.515908, "fLeviPiNeg");
+  fLevi[kKaon]  [kNeg] = fm2->GetLevi (mass[1], 0.1721, 6.927956, 0.191140, "fLeviKNeg");
+  fLevi[kProton][kNeg] = fm2->GetLevi (mass[2], 0.1782, 8.160362, 0.082091, "fLeviPNeg");
+  for(Int_t ipart = 0; ipart < kNPart; ipart++){
+    for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+      fLevi[ipart][icharge]->SetRange(0,4);
+    }    
+  }
+  
+
+  for(Int_t ipart = 0; ipart < kNPart; ipart++){
+    TCanvas * c1 = new TCanvas(TString("cRatio_")+partFlag[ipart], TString("cRatio_")+partFlag[ipart]);
+    hPosNegRatio[ipart]->Draw();
+    TF1 * fRatio = new TF1 (TString("fRatio")+partFlag[ipart], TString(fLevi[ipart][kPos]->GetName())+"/"+fLevi[ipart][kNeg]->GetName());
+    //    fRatio->Draw("same");
+    fRatio->SetRange(0,5);
+    if (doPrint) {
+      c1->Update();
+      gSystem->ProcessEvents();
+      c1->Print(TString(c1->GetName()) + ".eps");
+    }
+    
+  }
+
+
+  TCanvas * c2 = new TCanvas(TString("cRatio_KPi"),TString("cRatio_KPi"));  
+  hKPiRatio->Draw();
+  TLegend * lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
+  lMC->SetFillColor(kWhite);
+
+  //  gROOT->LoadMacro("GetE735Ratios.C");
+  GetE735Ratios(0,0)->Draw("EX0,same");
+  GetE735Ratios(0,1)->Draw("EX0,same");
+  GetE735Ratios(0,2)->Draw("EX0,same");
+  GetE735Ratios(0,3)->Draw("EX0,same");
+  hKPiRatio->SetMarkerStyle(20);
+  hKPiRatio->Draw("same");
+  
+  for(Int_t itune = 0; itune < kNTunes; itune++){
+    lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
+    hKPiRatioMC[itune]->SetLineWidth(2);    
+    hKPiRatioMC[itune]->Draw("same,chist");            
+  }
+  
+  lMC->Draw();
+
+
+  TLegend * l = new TLegend(  0.1879,  0.68,  0.54,0.92);
+  l->SetFillColor(kWhite);
+  l->AddEntry(hKPiRatio, "ALICE, #sqrt{s} = 900 GeV","lpf");
+  l->AddEntry(GetE735Ratios(0,0), "E735, #sqrt{s} = 300 GeV","lpf");
+  l->AddEntry(GetE735Ratios(0,1), "E735, #sqrt{s} = 540 GeV","lpf");
+  l->AddEntry(GetE735Ratios(0,2), "E735, #sqrt{s} = 1000 GeV","lpf");
+  l->AddEntry(GetE735Ratios(0,3), "E735, #sqrt{s} = 1800 GeV","lpf");
+  l->Draw();
+
+
+  TCanvas * c3 = new TCanvas(TString("cRatio_PPi"),TString("cRatio_PPi"));  
+  hPPiRatio->Draw();
+  hPPiRatio->SetMaximum(0.39);
+  lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
+  lMC->SetFillColor(kWhite);
+  for(Int_t itune = 0; itune < kNTunes; itune++){
+    lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
+    hKPiRatioMC[itune]->SetLineWidth(2);    
+    hKPiRatioMC[itune]->Draw("same,chist");            
+  }
+  
+  lMC->Draw();
+
+  if (doPrint) {
+    c2->Update();
+    gSystem->ProcessEvents();
+    c2->Print(TString(c2->GetName()) + ".eps");
+    c3->Update();
+    gSystem->ProcessEvents();
+    c3->Print(TString(c3->GetName()) + ".eps");
+  }
+
+
+}
+
diff --git a/PWG2/SPECTRA/Fit/GetE735Ratios.C b/PWG2/SPECTRA/Fit/GetE735Ratios.C
new file mode 100644 (file)
index 0000000..e363f5f
--- /dev/null
@@ -0,0 +1,112 @@
+TH1F* GetE735Ratios(Int_t ratio, Int_t energy) {
+
+  Int_t nbins = 10;
+
+  const Float_t binsKpi [] = {0.25,0.30, 0.35, 0.40,0.50,0.60,0.70,0.95,1.15,1.35,1.55};
+  const Float_t binsPpi [] = {0.30, 0.35, 0.40,0.50,0.60,0.70,0.95,1.15,1.35,1.55,1.75};
+  enum {kEn300,kEn540,kEn1000,kEn1800,kNenergy};
+
+  
+  TH1F * hKpi[kNenergy] ;
+  TH1F * hPpi[kNenergy] ;
+  
+  const Int_t markers[] = {24,25,27,28};
+
+  for(Int_t iener = 0; iener < kNenergy; iener++){
+    hKpi[iener] = new TH1F(TString("hKpi")+long(iener),TString("hKpi")+long(iener),nbins,binsKpi);
+    hPpi[iener] = new TH1F(TString("hPpi")+long(iener),TString("hPpi")+long(iener),nbins,binsPpi);
+    hKpi[iener]->SetMarkerStyle(markers[iener]);
+    hPpi[iener]->SetMarkerStyle(markers[iener]);    
+  }
+  
+
+
+  // 300
+  hKpi[kEn300]  ->SetBinContent(0 , 8.1  /100);
+  hKpi[kEn300]  ->SetBinContent(1 , 7.5  /100);
+  hKpi[kEn300]  ->SetBinContent(2 , 9.5  /100);
+  hKpi[kEn300]  ->SetBinContent(3 , 16.1 /100);
+  hKpi[kEn300]  ->SetBinContent(4 , 16.4 /100);
+  hKpi[kEn300]  ->SetBinContent(5 , 17.9 /100);
+  hKpi[kEn300]  ->SetBinContent(6 , 19.4 /100);
+  hKpi[kEn300]  ->SetBinContent(7 , 23.8 /100);
+  hKpi[kEn300]  ->SetBinContent(8 , 25.7 /100);
+  hKpi[kEn300]  ->SetBinContent(9 , 0    /100);
+  hKpi[kEn300]  ->SetBinError  (0 , 1.7/100);
+  hKpi[kEn300]  ->SetBinError  (1 , 1.4/100);
+  hKpi[kEn300]  ->SetBinError  (2 , 1.6/100);
+  hKpi[kEn300]  ->SetBinError  (3 , 1.5/100);
+  hKpi[kEn300]  ->SetBinError  (4 , 2.5/100);
+  hKpi[kEn300]  ->SetBinError  (5 , 3.5/100);
+  hKpi[kEn300]  ->SetBinError  (6 , 2.4/100);
+  hKpi[kEn300]  ->SetBinError  (7 , 4.9/100);
+  hKpi[kEn300]  ->SetBinError  (8 , 7.5/100);
+  hKpi[kEn300]  ->SetBinError  (9 , 0  /100);
+
+  hKpi[kEn540]  ->SetBinContent(0 , 3.3 /100);
+  hKpi[kEn540]  ->SetBinContent(1 , 12.7/100);
+  hKpi[kEn540]  ->SetBinContent(2 , 10.6/100);
+  hKpi[kEn540]  ->SetBinContent(3 , 13.1/100);
+  hKpi[kEn540]  ->SetBinContent(4 , 16.9/100);
+  hKpi[kEn540]  ->SetBinContent(5 , 17.3/100);
+  hKpi[kEn540]  ->SetBinContent(6 , 20.2/100);
+  hKpi[kEn540]  ->SetBinContent(7 , 21.2/100);
+  hKpi[kEn540]  ->SetBinContent(8 , 25.0/100);
+  hKpi[kEn540]  ->SetBinContent(9 , 31.2/100);
+  hKpi[kEn540]  ->SetBinError  (0 , 0.9/100);
+  hKpi[kEn540]  ->SetBinError  (1 , 1.1/100);
+  hKpi[kEn540]  ->SetBinError  (2 , 0.9/100);
+  hKpi[kEn540]  ->SetBinError  (3 , 0.7/100);
+  hKpi[kEn540]  ->SetBinError  (4 , 1.6/100);
+  hKpi[kEn540]  ->SetBinError  (5 , 1.7/100);
+  hKpi[kEn540]  ->SetBinError  (6 , 1.4/100);
+  hKpi[kEn540]  ->SetBinError  (7 , 2.1/100);
+  hKpi[kEn540]  ->SetBinError  (8 , 3.2/100);
+  hKpi[kEn540]  ->SetBinError  (9 , 4.9/100);
+
+  hKpi[kEn1000] ->SetBinContent(0 , 9.2 /100);
+  hKpi[kEn1000] ->SetBinContent(1 , 10.4/100);
+  hKpi[kEn1000] ->SetBinContent(2 , 9.3 /100);
+  hKpi[kEn1000] ->SetBinContent(3 , 11.8/100);
+  hKpi[kEn1000] ->SetBinContent(4 , 15.3/100);
+  hKpi[kEn1000] ->SetBinContent(5 , 19.8/100);
+  hKpi[kEn1000] ->SetBinContent(6 , 19.8/100);
+  hKpi[kEn1000] ->SetBinContent(7 , 25.9/100);
+  hKpi[kEn1000] ->SetBinContent(8 , 27.8/100);
+  hKpi[kEn1000] ->SetBinContent(9 , 32.8/100);
+  hKpi[kEn1000] ->SetBinError  (0 , 0.7/100);
+  hKpi[kEn1000] ->SetBinError  (1 , 0.7/100);
+  hKpi[kEn1000] ->SetBinError  (2 , 0.7/100);
+  hKpi[kEn1000] ->SetBinError  (3 , 0.7/100);
+  hKpi[kEn1000] ->SetBinError  (4 , 1.2/100);
+  hKpi[kEn1000] ->SetBinError  (5 , 1.5/100);
+  hKpi[kEn1000] ->SetBinError  (6 , 1.2/100);
+  hKpi[kEn1000] ->SetBinError  (7 , 1.9/100);
+  hKpi[kEn1000] ->SetBinError  (8 , 2.6/100);
+  hKpi[kEn1000] ->SetBinError  (9 , 3.9/100);
+
+  hKpi[kEn1800] ->SetBinContent(0 , 7.8 /100);
+  hKpi[kEn1800] ->SetBinContent(1 , 9.1 /100);
+  hKpi[kEn1800] ->SetBinContent(2 , 9.7 /100);
+  hKpi[kEn1800] ->SetBinContent(3 , 13.9/100);
+  hKpi[kEn1800] ->SetBinContent(4 , 16.2/100);
+  hKpi[kEn1800] ->SetBinContent(5 , 18.1/100);
+  hKpi[kEn1800] ->SetBinContent(6 , 21.6/100);
+  hKpi[kEn1800] ->SetBinContent(7 , 24.6/100);
+  hKpi[kEn1800] ->SetBinContent(8 , 25.8/100);
+  hKpi[kEn1800] ->SetBinContent(9 , 25.2/100);
+  hKpi[kEn1800] ->SetBinError  (0 , 0.3/100);
+  hKpi[kEn1800] ->SetBinError  (1 , 0.3/100);
+  hKpi[kEn1800] ->SetBinError  (2 , 0.3/100);
+  hKpi[kEn1800] ->SetBinError  (3 , 0.3/100);
+  hKpi[kEn1800] ->SetBinError  (4 , 0.6/100);
+  hKpi[kEn1800] ->SetBinError  (5 , 0.7/100);
+  hKpi[kEn1800] ->SetBinError  (6 , 0.7/100);
+  hKpi[kEn1800] ->SetBinError  (7 , 0.8/100);
+  hKpi[kEn1800] ->SetBinError  (8 , 0.9/100);
+  hKpi[kEn1800] ->SetBinError  (9 , 1.0/100);
+
+  if (ratio == 0) return hKpi[energy];
+  else            return hPpi[energy];
+
+}
diff --git a/PWG2/SPECTRA/Fit/StarPPSpectra.C b/PWG2/SPECTRA/Fit/StarPPSpectra.C
new file mode 100644 (file)
index 0000000..5d4e726
--- /dev/null
@@ -0,0 +1,330 @@
+TGraphErrors ** StarPPSpectra() {
+
+  TGraphErrors ** gStar = new TGraphErrors*[6];
+  const Int_t kNMaxPoints = 19;
+  Double_t pt[kNMaxPoints];
+  Double_t errpt[kNMaxPoints] = {0};
+  Double_t pim [kNMaxPoints];
+  Double_t pime[kNMaxPoints];
+  Double_t pip [kNMaxPoints];
+  Double_t pipe[kNMaxPoints];
+  Double_t km [kNMaxPoints];
+  Double_t kme[kNMaxPoints];
+  Double_t kp [kNMaxPoints];
+  Double_t kpe[kNMaxPoints];
+  Double_t pm [kNMaxPoints];
+  Double_t pme[kNMaxPoints];
+  Double_t pp [kNMaxPoints];
+  Double_t ppe[kNMaxPoints];
+
+  Int_t ipoint = 0;
+  pt[ipoint++] = 0.225;
+  pt[ipoint++] = 0.275;
+  pt[ipoint++] = 0.325;
+  pt[ipoint++] = 0.375;
+  pt[ipoint++] = 0.425;
+  pt[ipoint++] = 0.475;
+  pt[ipoint++] = 0.525;
+  pt[ipoint++] = 0.575;
+  pt[ipoint++] = 0.625;
+  pt[ipoint++] = 0.675;
+  pt[ipoint++] = 0.725;
+  pt[ipoint++] = 0.775;
+  pt[ipoint++] = 0.825;
+  pt[ipoint++] = 0.875;
+  pt[ipoint++] = 0.925;
+  pt[ipoint++] = 0.975;
+  pt[ipoint++] = 1.025;
+  pt[ipoint++] = 1.075;
+  pt[ipoint++] = 1.125;
+
+
+  // pi-
+  ipoint = 0;
+  pim[ipoint++]  = 2.02;
+  pim[ipoint++]  = 1.52;
+  pim[ipoint++]  = 1.13;
+  pim[ipoint++]  = 8.44*1e-1;
+  pim[ipoint++]  = 6.35*1e-1;
+  pim[ipoint++]  = 4.69*1e-1;
+  pim[ipoint++]  = 3.54*1e-1;
+  pim[ipoint++]  = 2.67*1e-1;
+  pim[ipoint++]  = 2.02*1e-1;
+  pim[ipoint++]  = 1.53*1e-1;
+  pim[ipoint++]  = 1.16*1e-1;
+  pim[ipoint++]  = 8.65*1e-2;
+
+  ipoint = 0;
+  pime[ipoint++] = 0.06;
+  pime[ipoint++] = 0.03;
+  pime[ipoint++] = 0.02;
+  pime[ipoint++] = 0.09*1e-1;
+  pime[ipoint++] = 0.07*1e-1;
+  pime[ipoint++] = 0.05*1e-1;
+  pime[ipoint++] = 0.04*1e-1;
+  pime[ipoint++] = 0.03*1e-1;
+  pime[ipoint++] = 0.04*1e-1;
+  pime[ipoint++] = 0.03*1e-1;
+  pime[ipoint++] = 0.03*1e-1;
+  pime[ipoint++] = 0.30*1e-2;
+
+  gStar[0] = new TGraphErrors(ipoint,pt,pim,errpt,pime);
+
+  // pi+
+  ipoint = 0;
+  pip[ipoint++]  = 2.07;
+  pip[ipoint++]  = 1.54;
+  pip[ipoint++]  = 1.14;
+  pip[ipoint++]  = 8.57*1e-1;
+  pip[ipoint++]  = 6.38*1e-1;
+  pip[ipoint++]  = 4.76*1e-1;
+  pip[ipoint++]  = 3.59*1e-1;
+  pip[ipoint++]  = 2.73*1e-1;
+  pip[ipoint++]  = 2.07*1e-1;
+  pip[ipoint++]  = 1.55*1e-1;
+  pip[ipoint++]  = 1.16*1e-1;
+  pip[ipoint++]  = 8.95*1e-2;
+
+  ipoint = 0;
+  pipe[ipoint++] = 0.06;
+  pipe[ipoint++] = 0.03;
+  pipe[ipoint++] = 0.02;
+  pipe[ipoint++] = 0.09*1e-1;
+  pipe[ipoint++] = 0.07*1e-1;
+  pipe[ipoint++] = 0.05*1e-1;
+  pipe[ipoint++] = 0.04*1e-1;
+  pipe[ipoint++] = 0.03*1e-1;
+  pipe[ipoint++] = 0.04*1e-1;
+  pipe[ipoint++] = 0.03*1e-1;
+  pipe[ipoint++] = 0.03*1e-1;
+  pipe[ipoint++] = 0.31*1e-2;
+
+  gStar[1] = new TGraphErrors(ipoint,pt,pip,errpt,pipe);
+
+  // k-
+  ipoint = 0;
+  km[ipoint++]  = 1.43*1e-1;
+  km[ipoint++]  = 1.26*1e-1;
+  km[ipoint++]  = 1.08*1e-1;
+  km[ipoint++]  = 8.77*1e-2;
+  km[ipoint++]  = 7.34*1e-2;
+  km[ipoint++]  = 6.17*1e-2;
+  km[ipoint++]  = 4.87*1e-2;
+  km[ipoint++]  = 4.11*1e-2;
+  km[ipoint++]  = 3.70*1e-2;
+  km[ipoint++]  = 2.86*1e-2;
+  km[ipoint++]  = 2.42*1e-2;
+
+  ipoint = 0;
+  kme[ipoint++] = 0.11*1e-1;
+  kme[ipoint++] = 0.05*1e-1;
+  kme[ipoint++] = 0.02*1e-1;
+  kme[ipoint++] = 0.24*1e-2;
+  kme[ipoint++] = 0.26*1e-2;
+  kme[ipoint++] = 0.63*1e-2;
+  kme[ipoint++] = 0.50*1e-2;
+  kme[ipoint++] = 0.42*1e-2;
+  kme[ipoint++] = 0.28*1e-2;
+  kme[ipoint++] = 0.22*1e-2;
+  kme[ipoint++] = 0.27*1e-2;
+
+  gStar[2] = new TGraphErrors(ipoint,pt,km,errpt,kme);
+
+
+  // k+
+  ipoint = 0;
+  kp[ipoint++]  = 1.52*1e-1;
+  kp[ipoint++]  = 1.30*1e-1;
+  kp[ipoint++]  = 1.08*1e-1;
+  kp[ipoint++]  = 9.16*1e-2;
+  kp[ipoint++]  = 7.47*1e-2;
+  kp[ipoint++]  = 6.26*1e-2;
+  kp[ipoint++]  = 5.26*1e-2;
+  kp[ipoint++]  = 4.41*1e-2;
+  kp[ipoint++]  = 3.81*1e-2;
+  kp[ipoint++]  = 3.06*1e-2;
+  kp[ipoint++]  = 2.49*1e-2;
+
+  ipoint = 0;
+  kpe[ipoint++] = 0.11*1e-1;
+  kpe[ipoint++] = 0.05*1e-1;
+  kpe[ipoint++] = 0.02*1e-1;
+  kpe[ipoint++] = 0.25*1e-2;
+  kpe[ipoint++] = 0.27*1e-2;
+  kpe[ipoint++] = 0.64*1e-2;
+  kpe[ipoint++] = 0.54*1e-2;
+  kpe[ipoint++] = 0.45*1e-2;
+  kpe[ipoint++] = 0.28*1e-2;
+  kpe[ipoint++] = 0.24*1e-2;
+  kpe[ipoint++] = 0.28*1e-2;
+
+  gStar[3] = new TGraphErrors(ipoint,pt,kp,errpt,kpe);
+
+
+  // pbar
+  ipoint = 0;
+  pm[ipoint++]  = 0;
+  pm[ipoint++]  = 0;
+  pm[ipoint++]  = 0;
+  pm[ipoint++]  = 5.54*1e-2;
+  pm[ipoint++]  = 4.81*1e-2;
+  pm[ipoint++]  = 4.25*1e-2;
+  pm[ipoint++]  = 3.77*1e-2;
+  pm[ipoint++]  = 3.38*1e-2;
+  pm[ipoint++]  = 2.78*1e-2;
+  pm[ipoint++]  = 2.49*1e-2;
+  pm[ipoint++]  = 2.03*1e-2;
+  pm[ipoint++]  = 1.75*1e-2;
+  pm[ipoint++]  = 1.52*1e-2;
+  pm[ipoint++]  = 1.27*1e-2;
+  pm[ipoint++]  = 1.05*1e-2;
+  pm[ipoint++]  = 8.95*1e-3;
+  pm[ipoint++]  = 7.36*1e-3;
+  pm[ipoint++]  = 6.59*1e-3;
+  pm[ipoint++]  = 5.21*1e-3;
+
+  ipoint = 0;
+  pme[ipoint++] = 0;
+  pme[ipoint++] = 0;
+  pme[ipoint++] = 0;
+  pme[ipoint++] = 0.13*1e-2;
+  pme[ipoint++] = 0.11*1e-2;
+  pme[ipoint++] = 0.10*1e-2;
+  pme[ipoint++] = 0.09*1e-2;
+  pme[ipoint++] = 0.08*1e-2;
+  pme[ipoint++] = 0.07*1e-2;
+  pme[ipoint++] = 0.06*1e-2;
+  pme[ipoint++] = 0.06*1e-2;
+  pme[ipoint++] = 0.06*1e-2;
+  pme[ipoint++] = 0.05*1e-2;
+  pme[ipoint++] = 0.04*1e-2;
+  pme[ipoint++] = 0.04*1e-2;
+  pme[ipoint++] = 0.34*1e-3;
+  pme[ipoint++] = 0.34*1e-3;
+  pme[ipoint++] = 0.32*1e-3;
+  pme[ipoint++] = 0.29*1e-3;
+
+  gStar[4] =  new TGraphErrors(ipoint,pt,pm,errpt,pme);
+
+
+  // p
+  ipoint = 0;
+  pp[ipoint++]  = 0;
+  pp[ipoint++]  = 0;
+  pp[ipoint++]  = 0;
+  pp[ipoint++]  = 0;
+  pp[ipoint++]  = 0;
+  pp[ipoint++]  = 5.07*1e-2;
+  pp[ipoint++]  = 4.69*1e-2;
+  pp[ipoint++]  = 4.08*1e-2;
+  pp[ipoint++]  = 3.42*1e-2;
+  pp[ipoint++]  = 2.87*1e-2;
+  pp[ipoint++]  = 2.41*1e-2;
+  pp[ipoint++]  = 2.13*1e-2;
+  pp[ipoint++]  = 1.82*1e-2;
+  pp[ipoint++]  = 1.54*1e-2;
+  pp[ipoint++]  = 1.31*1e-2;
+  pp[ipoint++]  = 1.11*1e-2;
+  pp[ipoint++]  = 9.78*1e-3;
+  pp[ipoint++]  = 8.56*1e-3;
+  pp[ipoint++]  = 7.38*1e-3;
+
+  ipoint = 0;
+  ppe[ipoint++] = 0;
+  ppe[ipoint++] = 0;
+  ppe[ipoint++] = 0;
+  ppe[ipoint++] = 0;
+  ppe[ipoint++] = 0;
+  ppe[ipoint++] = 0.25*1e-2;
+  ppe[ipoint++] = 0.19*1e-2;
+  ppe[ipoint++] = 0.14*1e-2;
+  ppe[ipoint++] = 0.10*1e-2;
+  ppe[ipoint++] = 0.07*1e-2;
+  ppe[ipoint++] = 0.07*1e-2;
+  ppe[ipoint++] = 0.06*1e-2;
+  ppe[ipoint++] = 0.05*1e-2;
+  ppe[ipoint++] = 0.05*1e-2;
+  ppe[ipoint++] = 0.04*1e-2;
+  ppe[ipoint++] = 0.04*1e-2;
+  ppe[ipoint++] = 0.40*1e-3;
+  ppe[ipoint++] = 0.37*1e-3;
+  ppe[ipoint++] = 0.38*1e-3;
+
+
+
+  gStar[5] = new TGraphErrors(ipoint,pt,pp,errpt,ppe);
+
+  gStar[0]->SetMarkerStyle(kOpenTriangleUp);
+  gStar[1]->SetMarkerStyle(kFullTriangleUp);
+  gStar[2]->SetMarkerStyle(kOpenCircle);
+  gStar[3]->SetMarkerStyle(kFullCircle);
+  gStar[4]->SetMarkerStyle(kOpenSquare);
+  gStar[5]->SetMarkerStyle(kFullSquare);
+
+  gStar[0]->SetMarkerColor(kBlack);
+  gStar[1]->SetMarkerColor(kBlack);
+  gStar[2]->SetMarkerColor(kRed);
+  gStar[3]->SetMarkerColor(kRed);
+  gStar[4]->SetMarkerColor(kBlue);
+  gStar[5]->SetMarkerColor(kBlue);
+
+  gStar[0]->SetLineColor (kBlack);
+  gStar[1]->SetLineColor (kBlack);
+  gStar[2]->SetLineColor (kRed);
+  gStar[3]->SetLineColor (kRed);
+  gStar[4]->SetLineColor (kBlue);
+  gStar[5]->SetLineColor (kBlue);
+
+  // Scale all for 2*pi and pt in order to go to dNdpt
+  for(Int_t istar = 0; istar < 6; istar++){
+    Int_t npoint = gStar[istar]->GetN();
+    for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
+      Double_t ptbin = gStar[istar]->GetX()[ipoint];
+      gStar[istar]->SetPoint(ipoint,
+                            gStar[istar]->GetX()[ipoint],
+                            gStar[istar]->GetY()[ipoint]*2*TMath::Pi()*ptbin);
+      gStar[istar]->SetPointError(ipoint,
+                                 gStar[istar]->GetEX()[ipoint],
+                                 gStar[istar]->GetEY()[ipoint]*2*TMath::Pi()*ptbin);
+    }
+    
+  }
+  
+
+
+//   gStar[0]->Draw("AP");
+//   gStar[1]->Draw("P");
+//   gStar[2]->Draw("P");
+//   gStar[3]->Draw("P");
+//   gStar[4]->Draw("P");
+//   gStar[5]->Draw("P");
+
+
+  return gStar;
+
+
+  
+  // table from PRC79  79 (2009) 34909
+//         pi-                   pi+               k-                 K+               pbar              p
+//  0.225  2.02+-0.06            2.07+-0.06        (1.43+-0.11)*10-1 (1.52+-0.11)*10-1
+//  0.275  1.52+-0.03            1.54+-0.03        (1.26+-0.05)*10-1 (1.30+-0.05)*10-1
+//  0.325  1.13+-0.02            1.14+-0.02        (1.08+-0.02)*10-1 (1.08+-0.02)*10-1
+//  0.375  (8.44+-0.09)*10-1     (8.57+-0.09)*10-1 (8.77+-0.24)*10-2 (9.16+-0.25)*10-2 (5.54+-0.13)*10-2
+//  0.425  (6.35+-0.07)*10-1     (6.38+-0.07)*10-1 (7.34+-0.26)*10-2 (7.47+-0.27)*10-2 (4.81+-0.11)*10-2
+//  0.475  (4.69+-0.05)*10-1     (4.76+-0.05)*10-1 (6.17+-0.63)*10-2 (6.26+-0.64)*10-2 (4.25+-0.10)*10-2 (5.07+-0.25)*10-2
+//  0.525  (3.54+-0.04)*10-1     (3.59+-0.04)*10-1 (4.87+-0.50)*10-2 (5.26+-0.54)*10-2 (3.77+-0.09)*10-2 (4.69+-0.19)*10-2
+//  0.575  (2.67+-0.03)*10-1     (2.73+-0.03)*10-1 (4.11+-0.42)*10-2 (4.41+-0.45)*10-2 (3.38+-0.08)*10-2 (4.08+-0.14)*10-2
+//  0.625  (2.02+-0.04)*10-1     (2.07+-0.04)*10-1 (3.70+-0.28)*10-2 (3.81+-0.28)*10-2 (2.78+-0.07)*10-2 (3.42+-0.10)*10-2
+//  0.675  (1.53+-0.03)*10-1     (1.55+-0.03)*10-1 (2.86+-0.22)*10-2 (3.06+-0.24)*10-2 (2.49+-0.06)*10-2 (2.87+-0.07)*10-2
+//  0.725  (1.16+-0.03)*10-1     (1.16+-0.03)*10-1 (2.42+-0.27)*10-2 (2.49+-0.28)*10-2 (2.03+-0.06)*10-2 (2.41+-0.07)*10-2
+//  0.775  (8.65+-0.30)*10-2     (8.95+-0.31)*10-2                                     (1.75+-0.06)*10-2 (2.13+-0.06)*10-2
+//  0.825                                                                              (1.52+-0.05)*10-2 (1.82+-0.05)*10-2
+//  0.875                                                                              (1.27+-0.04)*10-2 (1.54+-0.05)*10-2
+//  0.925                                                                              (1.05+-0.04)*10-2 (1.31+-0.04)*10-2
+//  0.975                                                                              (8.95+-0.34)*10-3 (1.11+-0.04)*10-2
+//  1.025                                                                              (7.36+-0.34)*10-3 (9.78+-0.40)*10-3
+//  1.075                                                                              (6.59+-0.32)*10-3 (8.56+-0.37)*10-3
+//  1.125                                                                              (5.21+-0.29)*10-3 (7.38+-0.38)*10-3
+
+
+  }
index 6c34388573fa7c05721c8d4cc3c1cf85896a0e56..d9b65facccd217b36ca6ea252c102b9080c4dbfe 100644 (file)
@@ -21,7 +21,10 @@ SRCS= SPECTRA/AliProtonAnalysisBase.cxx \
       SPECTRA/AliAnalysisCentralCutEvtMC.cxx \
       SPECTRA/AliAnalysisCentralCutEvtESD.cxx \
       SPECTRA/AliAnalysisCentralExtrapolate.cxx \
-      SPECTRA/AliAnalysisTaskCentral.cxx
+      SPECTRA/AliAnalysisTaskCentral.cxx \
+      SPECTRA/Fit/AliBWTools.cxx \
+      SPECTRA/Fit/AliBWFunc.cxx \
+      SPECTRA/Fit/AliLatexTable.cxx 
 
 HDRS= $(SRCS:.cxx=.h)