PWGLF/SPECTRA/PiKaPr/ITSsa/AliITSsadEdxFitter.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////
  19 // Class with the fits algorithms to be used in the identified       //
  20 // spectra analysis using the ITS in stand-alone mode                //
  21 // Author: E.Biolcati, biolcati@to.infn.it                           //
  22 //         F.Prino, prino@to.infn.it                                 //
  23 ///////////////////////////////////////////////////////////////////////
  24
  25 #include <Riostream.h>
  26 #include <TLatex.h>
  27 #include <TH1F.h>
  28 #include <TF1.h>
  29 #include <TH1D.h>
  30 #include <TLine.h>
  31 #include <TH2F.h>
  32 #include <TMath.h>
  33 #include <TGraph.h>
  34 #include <TGraphErrors.h>
  35 #include <TLegend.h>
  36 #include <TLegendEntry.h>
  37 #include <TStyle.h>
  38 #include <Rtypes.h>
  39 #include "AliITSsadEdxFitter.h"
  40
  41
  42 ClassImp(AliITSsadEdxFitter)
  43 //______________________________________________________________________
  44 AliITSsadEdxFitter::AliITSsadEdxFitter():TObject(),
  45   fExpPionMean(0.),
  46   fExpKaonMean(0.),
  47   fExpProtonMean(0.),
  48   fExpPionMeanRange(0.),
  49   fExpKaonMeanRange(0.),
  50   fExpProtonMeanRange(0.),
  51   fExpPionSigma(0.),
  52   fExpKaonSigma(0.),
  53   fExpProtonSigma(0.),
  54   fExpPionAmpl(0.),
  55   fIsMC(kFALSE),
  56   fOptInit(0),
  57   fFitOption("M0R+"),
  58   fITSpid(0)
  59 {
  60   // standard constructor
  61   for(Int_t i=0; i<9; i++)  fFitPars[i] = 0.;
  62   for(Int_t i=0; i<9; i++)  fFitParErrors[i] = 0.;
  63   for(Int_t i=0; i<3; i++)  fFitpar[i] = 0.;
  64   for(Int_t i=0; i<3; i++)  fFitparErr[i] = 0.;
  65   SetRangeStep1();
  66   SetRangeStep2();
  67   SetRangeStep3();
  68   SetRangeFinalStep();
  69   SetLimitsOnSigmaPion();
  70   SetLimitsOnSigmaKaon();
  71   SetLimitsOnSigmaProton();
  72   SetBinsUsedPion();
  73   SetBinsUsedKaon();
  74   SetBinsUsedProton();
  75   fITSpid=new AliITSPIDResponse(kFALSE);
  76 };
  77 //______________________________________________________________________
  78 AliITSsadEdxFitter::AliITSsadEdxFitter(Bool_t isMC):TObject(),
  79   fExpPionMean(0.),
  80   fExpKaonMean(0.),
  81   fExpProtonMean(0.),
  82   fExpPionMeanRange(0.),
  83   fExpKaonMeanRange(0.),
  84   fExpProtonMeanRange(0.),
  85   fExpPionSigma(0.),
  86   fExpKaonSigma(0.),
  87   fExpProtonSigma(0.),
  88   fExpPionAmpl(0.),
  89   fIsMC(isMC),
  90   fOptInit(0),
  91   fFitOption("M0R+"),
  92   fITSpid(0)
  93 {
  94   // standard constructor
  95   for(Int_t i=0; i<9; i++)  fFitPars[i] = 0.;
  96   for(Int_t i=0; i<9; i++)  fFitParErrors[i] = 0.;
  97   for(Int_t i=0; i<3; i++)  fFitpar[i] = 0.;
  98   for(Int_t i=0; i<3; i++)  fFitparErr[i] = 0.;
  99   SetRangeStep1();
 100   SetRangeStep2();
 101   SetRangeStep3();
 102   SetRangeFinalStep();
 103   SetLimitsOnSigmaPion();
 104   SetLimitsOnSigmaKaon();
 105   SetLimitsOnSigmaProton();
 106   SetBinsUsedPion();
 107   SetBinsUsedKaon();
 108   SetBinsUsedProton();
 109   fITSpid=new AliITSPIDResponse(isMC);
 110 };
 111
 112 //________________________________________________________
 113 Double_t AliITSsadEdxFitter::CalcSigma(Int_t code,Float_t x) const {
 114   // calculate the sigma 12-ott-2010
 115   Double_t p[2]={0.};
 116   Double_t xMinKaon=0.15; //minimum pt value to consider the kaon peak
 117   Double_t xMinProton=0.3;//minimum pt value to consider the proton peak
 118   if(fIsMC){
 119     if(code==211){
 120       p[0] = -1.20337e-04;
 121       p[1] = 1.13060e-01;
 122     }
 123     else if(code==321 && x>xMinKaon){
 124       p[0] = -2.39631e-03;
 125       p[1] = 1.15723e-01;
 126     }
 127     else if(code==2212 && x>xMinProton){
 128       p[0] = -8.34576e-03;
 129       p[1] = 1.34237e-01;
 130     }
 131     else return -1;
 132   }
 133   else{
 134     if(code==211){
 135       p[0] = -6.55200e-05;
 136       p[1] = 1.26657e-01;
 137     }
 138     else if(code==321 && x>xMinKaon){
 139       p[0] = -6.22639e-04;
 140       p[1] = 1.43289e-01;
 141     }
 142     else if(code==2212 && x>xMinProton){
 143       p[0] = -2.13243e-03;
 144       p[1] = 1.68614e-01;
 145     }
 146     else return -1;
 147   }
 148   return p[0]/(x*x)*TMath::Log(x)+p[1];
 149 }
 150
 151 //_______________________________________________________
 152 Int_t AliITSsadEdxFitter::InitializeMeanPosParam(Int_t code, Float_t x){
 153   // Set initial values for peak positions from ad-hoc parameterizations (Melinda, 5 april 2010)
 154
 155   Double_t p1[5]={0.};
 156   Double_t p2[5]={0.};
 157   Double_t p3[5]={0.};
 158   Double_t ep1[5]={0.};
 159   Double_t ep2[5]={0.};
 160   Double_t ep3[5]={0.};
 161   if(!fIsMC){ // data parameters
 162     if(code==211){
 163
 164       p1[0]=-4.21001e-04; p1[1]=-3.41133e-02; p1[2]=0.; p1[3]=0.; p1[4]=0.;
 165       ep1[0]= 3.62105e-06; ep1[1]= 1.52077e-04; ep1[2]= 0.; ep1[3]= 0.; ep1[4]= 0.;
 166
 167       p2[0] = 8.04038e+00; p2[1] = 1.51139e-01; p2[2] = 8.89984e-02; p2[3]= 0.; p2[4]= 0.;
 168       ep2[0]= 6.47330e-02; ep2[1]= 3.09501e-03; ep2[2]= 6.03642e-04; ep2[3]= 0.; ep2[4]= 0.;
 169
 170       p3[0] = 7.41939e+00; p3[1] =-1.31031e+00; p3[2] = 9.44093e-01; p3[3]= 0.; p3[4]= 0.;
 171       ep3[0]= 1.29862e-01; ep3[1]= 1.37370e-02; ep3[2]= 4.74995e-03; ep3[3]= 0.; ep3[4]= 0.;
 172
 173       fExpPionMean=p1[0]/(x*x)*TMath::Log(x)+p1[1];
 174       fExpPionMeanRange=TMath::Sqrt(TMath::Power(1/(x*x)*TMath::Log(x),2.)*ep1[0]*ep1[0]+ep1[1]*ep1[1] );
 175       fExpKaonMean=p2[0]*TMath::Landau(x,p2[1],p2[2],kFALSE);
 176       fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]);
 177       fExpProtonMean=-9999.;
 178       fExpProtonMeanRange=0.;
 179       if(x>0.2){
 180         fExpProtonMean=p3[0]*TMath::Exp(-0.5*TMath::Power((x-p3[1])/p3[2],2.)) ;
 181         fExpProtonMeanRange=TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]);
 182       }
 183
 184     }else if(code==321){
 185
 186       p1[0] = 8.35645e-01; p1[1] = 3.61638e+00; p1[2] =-2.51068e-01; p1[3]= 0.; p1[4]= 0.;
 187       ep1[0]= 6.63101e-04; ep1[1]= 2.33140e-02; ep1[2]= 1.93276e-03; ep1[3]= 0.; ep1[4]= 0.;
 188
 189       p2[0] = -4.75573e-02; p2[1] = -7.05764e-03; p2[2] =  7.96719e+00; p2[3]=  -7.48333e-02; p2[4]= 0.;
 190       ep2[0]= 1.70759e-02; ep2[1]= 2.95064e-04; ep2[2]= 2.70493e+00; ep2[3]= 1.14466e-02; ep2[4]= 0.;
 191
 192       p3[0] = 4.20188e+00; p3[1] = 3.95855e-01; p3[2] = 2.01159e-01; p3[3]=  0.; p3[4]=  0.;
 193       ep3[0]= 2.18608e-02; ep3[1]= 4.15293e-03; ep3[2]= 8.91846e-03; ep3[3]= 0.; ep3[4]= 0.;
 194
 195       fExpPionMean= p1[0]*TMath::Log(x)*TMath::Exp(-x*x*p1[1])+p1[2];
 196       fExpPionMeanRange=TMath::Sqrt(ep1[0]*ep1[0]+(x*x)*ep1[1]*ep1[1]+(x*x*x*x)*ep1[2]*ep1[2]);
 197       fExpKaonMean = p2[0]*TMath::Log(x)+p2[1]*TMath::Exp(-x*x*p2[2])/(x*x)+p2[3];
 198       fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]);
 199       fExpProtonMean = p3[0]*TMath::Landau(x, p3[1], p3[2],kFALSE);
 200       fExpProtonMeanRange=TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]);
 201
 202     }else if(code==2212){
 203
 204       p1[0] =1.48277e+00; p1[1] =1.14050e+00; p1[2] =1.01751e+02; p1[3]= -2.52768e-01; p1[4]= 0.;
 205       ep1[0]= 1.85491e-03; ep1[1]= 4.23200e-02; ep1[2]= 2.66855e+00; ep1[3]= 1.82240e-03; ep1[4]= 0.;
 206
 207       p2[0]= 3.80541e-01 ; p2[1]=-3.79158e-02; p2[2]=-4.81653e-01; p2[3]= 0.; p2[4]= 0.;
 208       ep2[0]= 1.43643e-01; ep2[1]= 2.47156e-02; ep2[2]= 1.06836e-01; ep2[3]= 0.; ep2[4]= 0.;
 209
 210       p3[0] = 1.81168e-01; p3[1] = -6.69364e-01; p3[2] =  2.10685e+01; p3[3]=   5.14868e-02; p3[4]=   0.;
 211       ep3[0]= 4.61718e-02; ep3[1]= 8.21949e-02; ep3[2]= 4.46134e+00; ep3[3]= 2.69751e-02; ep3[4]= 0.;
 212
 213       fExpPionMean= p1[0]*TMath::Log(x)+p1[1]*TMath::Exp(-x*x*p1[2])+p1[3];
 214       fExpPionMeanRange=TMath::Sqrt(ep1[0]*ep1[0]+(x*x)*ep1[1]*ep1[1]+(x*x*x*x)*ep1[2]*ep1[2]);
 215       fExpKaonMean = p2[0]*TMath::Log(x)+p2[1]/(x)+p2[2];
 216       fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]+(x*x*x*x*x*x)*ep2[3]*ep2[3]+(x*x*x*x*x*x*x*x)*ep2[4]*ep2[4]);
 217       fExpProtonMean= p3[0]*TMath::Log(x)+p3[1]*TMath::Exp(-x*x*p3[2])+p3[3];
 218       fExpProtonMeanRange=TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]+(x*x*x*x*x*x)*ep3[3]*ep3[3]+(x*x*x*x*x*x*x*x)*ep3[4]*ep3[4]);
 219     }else{
 220       printf("ERROR: Wrong particle code %d\n",code);
 221       return -1;
 222     }
 223   }else{ // MC parameters
 224     if(code==211){
 225
 226       p1[0]= -3.28744e-04; p1[1]= -1.78123e-02; p1[2]= 0.; p1[3]= 0.; p1[4]= 0.;
 227       ep1[0]= 2.03102e-06; ep1[1]= 7.65859e-05; ep1[2]= 0.; ep1[3]= 0.; ep1[4]= 0.;
 228
 229       p2[0] = 8.84991e+00 ; p2[1] = 1.32096e-01 ; p2[2] = 9.42394e-02 ; p2[3]=0.; p2[4]=0.;
 230       ep2[0]= 4.78123e-02; ep2[1]= 1.99037e-03; ep2[2]= 2.85136e-04; ep2[3]= 0.; ep2[4]= 0.;
 231
 232       p3[0] =  6.25025e+00; p3[1] = -1.08315e+00; p3[2] =  8.86474e-01; p3[3]= 0.; p3[4]= 0.;
 233       ep3[0]= 1.42379e+00; ep3[1]= 2.34807e-01; ep3[2]= 6.78015e-02; ep3[3]= 0.; ep3[4]= 0.;
 234
 235       fExpPionMean= p1[0]/(x*x)*TMath::Log(x)+p1[1];
 236       fExpPionMeanRange=TMath::Sqrt(TMath::Power(1/(x*x)*TMath::Log(x),2.)*ep1[0]*ep1[0]+ep1[1]*ep1[1] );
 237       fExpKaonMean = p2[0]*TMath::Landau(x,p2[1],p2[2],kFALSE);
 238       fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]);
 239       fExpProtonMean=-9999.;
 240       fExpProtonMeanRange=0.;
 241       if(x>0.2){
 242         fExpProtonMean = p3[0]*TMath::Exp(-0.5*TMath::Power((x-p3[1])/p3[2],2.));
 243         fExpProtonMeanRange = TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]);
 244       }
 245
 246     }else if(code==321){
 247
 248       p1[0] = 9.19679e-01; p1[1] = 5.07077e-03; p1[2] = 5.48262e-05; p1[3]=  8.65674e-02; p1[4]= 0.;
 249       ep1[0]= 4.03325e-04; ep1[1]= 8.85421e-06; ep1[2]= 1.77437e+00; ep1[3]= 9.00492e-03; ep1[4]= 0.;
 250
 251       p2[0] = -2.05127e-01; p2[1] = -2.05303e-03; p2[2] =  3.10090e-06; p2[3]=  -2.19936e-01; p2[4]= 0.;
 252       ep2[0]= 3.41607e-03; ep2[1]= 1.03661e-04; ep2[2]= 9.50776e-01; ep2[3]= 3.41148e-03; ep2[4]= 0.;
 253
 254       p3[0] = 4.68022e+00; p3[1] = 3.02178e-01; p3[2] = 2.47185e-01; p3[3]= 0.; p3[4]= 0.;
 255       ep3[0]= 1.63926e-02; ep3[1]= 6.42072e-03; ep3[2]= 7.24079e-03; ep3[3]= 0.; ep3[4]= 0.;
 256
 257       fExpPionMean=p1[0]*TMath::Log(x)+p1[1]*TMath::Exp(-x*x*p1[2])/(x*x)+p1[3];
 258       fExpPionMeanRange=TMath::Sqrt(ep1[0]*ep1[0]+(x*x)*ep1[1]*ep1[1]+(x*x*x*x)*ep1[2]*ep1[2]);
 259       fExpKaonMean=p2[0]*TMath::Log(x)+p2[1]*TMath::Exp(-x*x*p2[2])/(x*x+p2[3]);
 260       fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]);
 261       fExpProtonMean= p3[0]*TMath::Landau(x, p3[1], p3[2]);
 262       fExpProtonMeanRange=TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]);
 263
 264     }else if(code==2212){
 265
 266       p1[0] = 1.05826e+00; p1[1]=8.59088e-01; p1[2]=9.39446e+01; p1[3]= -4.86736e-01; p1[4]=  0.;
 267       ep1[0]= 1.07813e-03; ep1[1]= 1.15536e-02; ep1[2]= 9.74843e-01; ep1[3]= 1.24205e-03; ep1[4]= 0.;
 268
 269       p2[0]=4.91075e-01; p2[1]=6.27855e-01; p2[2] = 1.09965e+01; p2[3] =-4.24656e-01; p2[4] = 0.;
 270       ep2[0]=3.00310e-02; ep2[1]= 4.42650e-02; ep2[2]= 3.40324e-01; ep2[3]= 2.46448e-02; ep2[4]= 0. ;
 271
 272       p3[0] = 6.75789e-01; p3[1] = 1.23129e+00; p3[2] = 3.46294e+00; p3[3]= -2.48336e-02; p3[4]= 0.;
 273       ep3[0]= 2.90342e-02; ep3[1]= 4.34524e-02; ep3[2]= 7.18743e-02; ep3[3]= 1.08694e-02; ep3[4]= 0.;
 274
 275       fExpPionMean= p1[0]*TMath::Log(x)+p1[1]*TMath::Exp(-x*x*p1[2])+p1[3];
 276       fExpPionMeanRange=TMath::Sqrt(ep1[0]*ep1[0]+(x*x)*ep1[1]*ep1[1]+(x*x*x*x)*ep1[2]*ep1[2]);
 277       fExpKaonMean = p2[0]*TMath::Log(x)+p2[1]*TMath::Exp(-x*x*p2[2])+p2[3];
 278       fExpKaonMeanRange=TMath::Sqrt(ep2[0]*ep2[0]+(x*x)*ep2[1]*ep2[1]+(x*x*x*x)*ep2[2]*ep2[2]+(x*x*x*x*x*x)*ep2[3]*ep2[3]+(x*x*x*x*x*x*x*x)*ep2[4]*ep2[4]);
 279       fExpProtonMean=  p3[0]*TMath::Log(x)+p3[1]*TMath::Exp(-x*x*p3[2])+p3[3];
 280       fExpProtonMeanRange=TMath::Sqrt(ep3[0]*ep3[0]+(x*x)*ep3[1]*ep3[1]+(x*x*x*x)*ep3[2]*ep3[2]+(x*x*x*x*x*x)*ep3[3]*ep3[3]+(x*x*x*x*x*x*x*x)*ep3[4]*ep3[4]);
 281     }else{
 282       printf("ERROR: Wrong particle code %d\n",code);
 283       return -1;
 284     }
 285   }
 286   return 0;
 287 }
 288
 289 //_______________________________________________________
 290 Int_t AliITSsadEdxFitter::InitializeMeanPosBB(Int_t code, Float_t pt){
 291   // computes peak positions from parametrized BB
 292
 293   Double_t massp=AliPID::ParticleMass(AliPID::kProton);
 294   Double_t massk=AliPID::ParticleMass(AliPID::kKaon);
 295   Double_t masspi=AliPID::ParticleMass(AliPID::kPion);
 296   Bool_t isSA=kTRUE;
 297   Float_t sinthmed=0.8878;
 298   Float_t p=pt/sinthmed;
 299   Double_t bethep=fITSpid->Bethe(p,massp,isSA);
 300   Double_t bethek=fITSpid->Bethe(p,massk,isSA);
 301   Double_t bethepi=fITSpid->Bethe(p,masspi,isSA);
 302   Double_t betheref=bethepi;
 303   if(TMath::Abs(code)==321) betheref=bethek;
 304   else if(TMath::Abs(code)==2212) betheref=bethep;
 305   fExpPionMean=TMath::Log(bethepi)-TMath::Log(betheref);
 306   fExpKaonMean=TMath::Log(bethek)-TMath::Log(betheref);
 307   fExpProtonMean=TMath::Log(bethep)-TMath::Log(betheref);
 308   return 0;
 309 }
 310 //________________________________________________________
 311 Bool_t AliITSsadEdxFitter::IsGoodBin(Int_t bin,Int_t code){
 312   //method to select the bins used for the analysis only
 313   Bool_t retvalue=kTRUE;
 314   TLine *l[2]; //for the cross
 315   l[0]=new TLine(-2.1,0,2.,100.);
 316   l[1]=new TLine(-1.9,120,2.,0.);
 317   for(Int_t j=0;j<2;j++){
 318     l[j]->SetLineColor(4);
 319     l[j]->SetLineWidth(5);
 320   }
 321
 322   if(code==211 && (bin<fBinsUsedPion[0] || bin>fBinsUsedPion[1])){
 323     for(Int_t j=0;j<2;j++) l[j]->Draw("same");
 324     retvalue=kFALSE;
 325   }
 326   if(code==321 && (bin<fBinsUsedKaon[0] || bin>fBinsUsedKaon[1])){
 327     for(Int_t j=0;j<2;j++) l[j]->Draw("same");
 328     retvalue=kFALSE;
 329   }
 330   if(code==2212 && (bin<fBinsUsedProton[0] || bin>fBinsUsedProton[1])){
 331     for(Int_t j=0;j<2;j++) l[j]->Draw("same");
 332     retvalue=kFALSE;
 333   }
 334   return retvalue;
 335 }
 336
 337 //________________________________________________________
 338 Double_t SingleGausTail(const Double_t *x, const Double_t *par){
 339   //single gaussian with exponential tail
 340   Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
 341   Double_t xx = x[0];
 342   Double_t mean1 = par[1];
 343   Double_t sigma1 = par[2];
 344   Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);
 345   Double_t sigmaSquare1 = sigma1 * sigma1;
 346   Double_t xdiv1 = mean1 + par[3] * sigma1;
 347   Double_t y1=0.0;
 348   if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);
 349   else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));
 350   return y1;
 351 }
 352
 353 //________________________________________________________
 354 Double_t DoubleGausTail(const Double_t *x, const Double_t *par){
 355   //sum of two gaussians with exponential tail
 356   Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
 357   Double_t xx = x[0];
 358   Double_t mean1 = par[1];
 359   Double_t sigma1 = par[2];
 360   Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);
 361   Double_t sigmaSquare1 = sigma1 * sigma1;
 362   Double_t xdiv1 = mean1 + par[3] * sigma1;
 363   Double_t y1=0.0;
 364   if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);
 365   else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));
 366
 367   Double_t mean2 = par[6];
 368   Double_t sigma2 = par[7];
 369   Double_t xNormSquare2 = (xx - mean2) * (xx - mean2);
 370   Double_t sigmaSquare2 = sigma2 * sigma2;
 371   Double_t xdiv2 = mean2 + par[8] * sigma2;
 372   Double_t y2=0.0;
 373   if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2);
 374   else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8]));
 375   return y1+y2;
 376 }
 377
 378 //________________________________________________________
 379 Double_t FinalGausTail(const Double_t *x, const Double_t *par){
 380   //sum of three gaussians with exponential tail
 381   Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
 382   Double_t xx = x[0];
 383   Double_t mean1 = par[1];
 384   Double_t sigma1 = par[2];
 385   Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);
 386   Double_t sigmaSquare1 = sigma1 * sigma1;
 387   Double_t xdiv1 = mean1 + par[3] * sigma1;
 388   Double_t y1=0.0;
 389   if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);
 390   else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));
 391
 392   Double_t mean2 = par[6];
 393   Double_t sigma2 = par[7];
 394   Double_t xNormSquare2 = (xx - mean2) * (xx - mean2);
 395   Double_t sigmaSquare2 = sigma2 * sigma2;
 396   Double_t xdiv2 = mean2 + par[8] * sigma2;
 397   Double_t y2=0.0;
 398   if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2);
 399   else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8]));
 400
 401   Double_t mean3 = par[11];
 402   Double_t sigma3 = par[12];
 403   Double_t xNormSquare3 = (xx - mean3) * (xx - mean3);
 404   Double_t sigmaSquare3 = sigma3 * sigma3;
 405   Double_t xdiv3 = mean3 + par[13] * sigma3;
 406   Double_t y3=0.0;
 407   if(xx < xdiv3) y3 = par[10]/(s2pi*par[12]) * TMath::Exp(-0.5 * xNormSquare3 / sigmaSquare3);
 408   else y3 = TMath::Exp(-par[14]*(xx-xdiv3)) * par[10] / (s2pi*par[12]) * TMath::Exp(-0.5*(par[13]*par[13]));
 409   return y1+y2+y3;
 410 }
 411
 412 //______________________________________________________________________
 413 void AliITSsadEdxFitter::CalcResidual(TH1F *h,TF1 *fun,TGraph *gres) const{
 414   //code to calculate the difference fit function and histogram point (residual)
 415   //to use as goodness test for the fit
 416   Int_t ipt=0;
 417   Double_t x=0.,yhis=0.,yfun=0.;
 418   for(Int_t i=0;i<h->GetNbinsX();i++){
 419     x=(h->GetBinLowEdge(i+2)+h->GetBinLowEdge(i+1))/2;
 420     yfun=fun->Eval(x);
 421     yhis=h->GetBinContent(i+1);
 422     if(yhis>0. && yfun>0.) {
 423       gres->SetPoint(ipt,x,(yhis-yfun)/yhis);
 424       ipt++;
 425     }
 426   }
 427   return;
 428 }
 429
 430
 431 //________________________________________________________
 432 Double_t SingleGausStep(const Double_t *x, const Double_t *par){
 433   //single normalized gaussian
 434   Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
 435   Double_t xx = x[0];
 436   Double_t mean1 = par[1];
 437   Double_t sigma1 = par[2];
 438   Double_t xNorm1Square = (xx - mean1) * (xx - mean1);
 439   Double_t sigma1Square = sigma1 * sigma1;
 440   Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNorm1Square / sigma1Square);
 441   return step1;
 442 }
 443
 444 //________________________________________________________
 445 Double_t DoubleGausStep(const Double_t *x, const Double_t *par){
 446   //sum of two normalized gaussians
 447   Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
 448   Double_t xx = x[0];
 449   Double_t mean1 = par[1];
 450   Double_t sigma1 = par[2];
 451   Double_t xNorm1Square = (xx - mean1) * (xx - mean1);
 452   Double_t sigma1Square = sigma1 * sigma1;
 453   Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square );
 454   Double_t mean2 = par[4];
 455   Double_t sigma2 = par[5];
 456   Double_t xNorm2Square = (xx - mean2) * (xx - mean2);
 457   Double_t sigma2Square = sigma2 * sigma2;
 458   Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square );
 459   return step1+step2;
 460 }
 461
 462 //________________________________________________________
 463 Double_t FinalGausStep(const Double_t *x, const Double_t *par){
 464   //sum of three normalized gaussians
 465   Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
 466   Double_t xx = x[0];
 467   Double_t mean1 = par[1];
 468   Double_t sigma1 = par[2];
 469   Double_t xNorm1Square = (xx - mean1) * (xx - mean1);
 470   Double_t sigma1Square = sigma1 * sigma1;
 471   Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square );
 472   Double_t mean2 = par[4];
 473   Double_t sigma2 = par[5];
 474   Double_t xNorm2Square = (xx - mean2) * (xx - mean2);
 475   Double_t sigma2Square = sigma2 * sigma2;
 476   Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square );
 477   Double_t mean3 = par[7];
 478   Double_t sigma3 = par[8];
 479   Double_t xNorm3Square = (xx - mean3) * (xx - mean3);
 480   Double_t sigma3Square = sigma3 * sigma3;
 481   Double_t step3 = par[6]/(s2pi*par[8]) * TMath::Exp( - 0.5 * xNorm3Square / sigma3Square);
 482   return step1+step2+step3;
 483 }
 484
 485 //______________________________________________________________________
 486 Double_t AliITSsadEdxFitter::GausPlusTail(const Double_t x, const Double_t mean, Double_t rms, Double_t c, Double_t slope, Double_t cut ) const{
 487   //gaussian with an exponential tail on the right side
 488   Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi()));
 489   Double_t returnvalue=0.0;
 490   Double_t n=0.5*(1.0+TMath::Erf(cut/TMath::Sqrt(2.0)))+TMath::Exp(-cut*cut*0.5)*factor/(TMath::Abs(rms)*slope);
 491   if (x<mean+cut*rms) returnvalue=TMath::Exp(-1.0*(x-mean)*(x-mean)/(2.0*rms*rms))*factor/TMath::Abs(rms);
 492   else returnvalue=TMath::Exp(slope*(mean+cut*rms-x))*TMath::Exp(-cut*cut*0.5)*factor/TMath::Abs(rms);
 493   return c*returnvalue/n;
 494 }
 495
 496
 497 //______________________________________________________________________
 498 Double_t AliITSsadEdxFitter::GausOnBackground(const Double_t* x, const Double_t *par) const {
 499   //gaussian with a background parametrisation
 500   //cout<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<<" "<<x[0]<< endl;
 501   Double_t returnvalue=0.0;
 502   Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi()));
 503   if(par[6]<0.0) returnvalue+=TMath::Exp(-1.0*(x[0]-par[0])*(x[0]-par[0])/(2.0*par[1]*par[1]))*par[2]* factor/TMath::Abs(par[1]);
 504   else returnvalue+=GausPlusTail(x[0], par[0], par[1],par[2], par[6], 1.2);
 505   returnvalue+=par[3]*TMath::Exp((par[5]-x[0])*par[4]);
 506   return returnvalue;
 507 }
 508
 509 //______________________________________________________________________
 510 void AliITSsadEdxFitter::DrawFitFunction(TF1 *fun) const {
 511   //code to draw the final fit function and the single gaussians used
 512   //to extract the yields for the 3 species
 513   TF1 *fdraw1=new TF1("fdraw1",SingleGausStep,-3.5,3.5,3);
 514   TF1 *fdraw2=new TF1("fdraw2",SingleGausStep,-3.5,3.5,3);
 515   TF1 *fdraw3=new TF1("fdraw3",SingleGausStep,-3.5,3.5,3);
 516   fdraw1->SetParameter(0,fun->GetParameter(0));
 517   fdraw1->SetParameter(1,fun->GetParameter(1));
 518   fdraw1->SetParameter(2,fun->GetParameter(2));
 519   fdraw2->SetParameter(0,fun->GetParameter(3));
 520   fdraw2->SetParameter(1,fun->GetParameter(4));
 521   fdraw2->SetParameter(2,fun->GetParameter(5));
 522   fdraw3->SetParameter(0,fun->GetParameter(6));
 523   fdraw3->SetParameter(1,fun->GetParameter(7));
 524   fdraw3->SetParameter(2,fun->GetParameter(8));
 525
 526   fdraw1->SetLineColor(6);//color code
 527   fdraw2->SetLineColor(2);
 528   fdraw3->SetLineColor(4);
 529
 530   fdraw1->SetLineStyle(2);//dot lines
 531   fdraw2->SetLineStyle(2);
 532   fdraw3->SetLineStyle(2);
 533
 534   fun->SetLineWidth(3);
 535   fdraw1->DrawCopy("same");
 536   fdraw2->DrawCopy("same");
 537   fdraw3->DrawCopy("same");
 538   fun->DrawCopy("same");
 539
 540   TLatex *ltx[3];
 541   for(Int_t j=0;j<3;j++){
 542     ltx[0]=new TLatex(0.13,0.25,"pions");
 543     ltx[1]=new TLatex(0.13,0.20,"kaons");
 544     ltx[2]=new TLatex(0.13,0.15,"protons");
 545     ltx[0]->SetTextColor(6);
 546     ltx[1]->SetTextColor(2);
 547     ltx[2]->SetTextColor(4);
 548     ltx[j]->SetNDC();
 549     ltx[j]->Draw();
 550   }
 551   return;
 552 }
 553
 554 //______________________________________________________________________
 555 void AliITSsadEdxFitter::Initialize(TH1F* h, Int_t code,Int_t bin){
 556   //code to get the expected values to use for fitting
 557
 558   Double_t xbins[23]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};
 559   Double_t pt=(xbins[bin+1]+xbins[bin])/2;
 560
 561   //draw and label
 562   h->SetTitle(Form("p_{t}=[%1.2f,%1.2f], code=%d",xbins[bin],xbins[bin+1],code));
 563   h->GetXaxis()->SetTitle("[ln dE/dx]_{meas} - [ln dE/dx(i)]_{calc}");
 564   h->GetYaxis()->SetTitle("counts");
 565   h->Draw("e");
 566   h->SetFillColor(11);
 567
 568   //expected values
 569   fExpPionAmpl = h->GetMaximum()/(h->GetRMS()*TMath::Sqrt(2*TMath::Pi()));
 570   if(fOptInit==0) InitializeMeanPosBB(code, pt);
 571   else if(fOptInit==1) InitializeMeanPosParam(code, pt);
 572   fExpPionSigma = CalcSigma(211,pt); //expected sigma values
 573   fExpKaonSigma = CalcSigma(321,pt);
 574   fExpProtonSigma = CalcSigma(2212,pt);
 575
 576 }
 577
 578
 579 //________________________________________________________
 580 void AliITSsadEdxFitter::DoFit(TH1F *h, Int_t bin, Int_t signedcode, TGraph *gres){
 581   // 3-gaussian fit to log(dedx)-log(dedxBB) histogram
 582   // pt bin from 0 to 20, code={211,321,2212}
 583   // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed
 584   // final step: refit all using the parameters and tollerance limits (+-20%)
 585
 586   TF1 *fstep1, *fstep2, *fstep3, *fstepTot;
 587   Int_t code=TMath::Abs(signedcode);
 588   Initialize(h,code,bin);
 589   PrintInitialValues();
 590   if(!IsGoodBin(bin,code)) return;
 591
 592   printf("___________________________________________________________________ First Step: pions\n");
 593   fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);
 594   fstep1->SetParameter(0,fExpPionAmpl);       //initial ampl pion
 595   fstep1->SetParameter(1,fExpPionMean);       //initial mean pion
 596   fstep1->SetParameter(2,fExpPionSigma); //initial sigma pion
 597   fstep1->SetParLimits(0,0.,fExpPionAmpl*1.2);    //limits ampl pion
 598   fstep1->SetParLimits(1,fRangeStep4[0],fRangeStep4[1]);                                     //limits mean pion (dummy)
 599   fstep1->SetParLimits(2,fExpPionSigma*fLimitsOnSigmaPion[0],fExpPionSigma*fLimitsOnSigmaPion[1]); //limits sigma pion
 600
 601   if(fExpPionSigma>0) h->Fit(fstep1,fFitOption.Data(),"",fExpPionMean+fRangeStep1[0],fExpPionMean+fRangeStep1[1]);//first fit
 602   else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);
 603
 604   printf("___________________________________________________________________ Second Step: kaons\n");
 605   fstep2 = new TF1("fstep2",DoubleGausStep,fRangeStep4[0],fRangeStep4[1],6);
 606   fstep2->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion
 607   fstep2->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion
 608   fstep2->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion
 609   fstep2->SetParameter(3,fstep1->GetParameter(0)/8.); //initial ampl kaon
 610   fstep2->SetParameter(4,fExpKaonMean);                //initial mean kaon
 611   fstep2->SetParameter(3,fExpKaonSigma);                 //initial sigma kaon
 612   fstep2->SetParLimits(3,0.,fstep1->GetParameter(0));                                         //limits ampl kaon
 613   fstep2->SetParLimits(4,fstep1->GetParameter(1),fRangeStep4[1]);                             //limits mean kaon
 614   fstep2->SetParLimits(5,fExpKaonSigma*fLimitsOnSigmaKaon[0],fExpKaonSigma*fLimitsOnSigmaKaon[1]);  //limits sigma kaon
 615
 616   if(fExpKaonSigma>0) h->Fit(fstep2,fFitOption.Data(),"",fExpKaonMean+fRangeStep2[0],fExpKaonMean+fRangeStep2[1]);//second fit
 617   else for(Int_t npar=3;npar<6;npar++) fstep2->FixParameter(npar,0.00001);
 618
 619
 620   printf("___________________________________________________________________ Third Step: protons\n");
 621   fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);
 622   fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion
 623   fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion
 624   fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion
 625   fstep3->FixParameter(3,fstep2->GetParameter(3)); //fixed ampl kaon
 626   fstep3->FixParameter(4,fstep2->GetParameter(4)); //fixed mean kaon
 627   fstep3->FixParameter(5,fstep2->GetParameter(5)); //fidex sigma kaon
 628   fstep3->SetParameter(6,fstep2->GetParameter(0)/16.); //initial ampl proton
 629   fstep3->SetParameter(7,fExpProtonMean);               //initial mean proton
 630   fstep3->SetParameter(8,fExpProtonSigma);                //initial sigma proton
 631   fstep3->SetParLimits(6,0.,fstep2->GetParameter(0));                                                //limits ampl proton
 632   fstep3->SetParLimits(7,fstep2->GetParameter(4),fRangeStep4[1]);                                    //limits mean proton
 633   fstep3->SetParLimits(8,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton
 634
 635   if(fExpProtonSigma>0) h->Fit(fstep3,fFitOption.Data(),"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//third fit
 636   else for(Int_t npar=6;npar<9;npar++) fstep3->FixParameter(npar,0.00001);
 637
 638   printf("___________________________________________________________________ Final Step: refit all\n");
 639   fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);
 640   fstepTot->SetLineColor(1);
 641
 642   Double_t initialParametersStepTot[9];
 643   initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian
 644   initialParametersStepTot[1] = fstep1->GetParameter(1);
 645   initialParametersStepTot[2] = fstep1->GetParameter(2);
 646
 647   initialParametersStepTot[3] = fstep2->GetParameter(3);//second gaussian
 648   initialParametersStepTot[4] = fstep2->GetParameter(4);
 649   initialParametersStepTot[5] = fstep2->GetParameter(5);
 650
 651   initialParametersStepTot[6] = fstep3->GetParameter(6);//third gaussian
 652   initialParametersStepTot[7] = fstep3->GetParameter(7);
 653   initialParametersStepTot[8] = fstep3->GetParameter(8);
 654
 655   fstepTot->SetParameters(initialParametersStepTot); //initial parameters
 656   fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl pion
 657   fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion
 658   fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion
 659   fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon
 660   fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon
 661   fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon
 662   fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl proton
 663   fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton
 664   fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton
 665
 666   h->Fit(fstepTot,fFitOption.Data(),"",fRangeStep4[0],fRangeStep4[1]);//refit all
 667
 668   for(Int_t j=0;j<9;j++) {
 669     fFitPars[j] = fstepTot->GetParameter(j);
 670     fFitParErrors[j] = fstepTot->GetParError(j);
 671   }
 672
 673   //************************************* storing parameter to calculate the yields *******
 674
 675   Int_t chpa=0;
 676   if(code==321) chpa=3;
 677   if(code==2212) chpa=6;
 678   for(Int_t j=0;j<3;j++) {
 679     fFitpar[j] = fstepTot->GetParameter(j+chpa);
 680     fFitparErr[j] = fstepTot->GetParError(j+chpa);
 681   }
 682
 683   DrawFitFunction(fstepTot);
 684   CalcResidual(h,fstepTot,gres);
 685   return;
 686 }
 687
 688 //________________________________________________________
 689 void AliITSsadEdxFitter::DoFitProton(TH1F *h, Int_t bin, Int_t signedcode, TGraph *gres){
 690   // 3-gaussian fit to log(dedx)-log(dedxBB) histogram
 691   // pt bin from 0 to 20, code={211,321,2212}
 692   // first step: pion peak, second step: proton peak, third step: kaon peak
 693   // final step: refit all using the parameters
 694   TF1 *fstep1, *fstep2, *fstep3, *fstepTot;
 695   Int_t code=TMath::Abs(signedcode);
 696   Initialize(h,code,bin);
 697   PrintInitialValues();
 698   if(!IsGoodBin(bin,code)) return;
 699
 700   printf("___________________________________________________________________ First Step: pion\n");
 701   fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);
 702   fstep1->SetParameter(0,fExpPionAmpl);       //initial ampl pion
 703   fstep1->SetParameter(1,fExpPionMean);       //initial mean pion
 704   fstep1->SetParameter(2,fExpPionSigma); //initial sigma pion
 705   fstep1->SetParLimits(0,0,fExpPionAmpl*1.2);                                                          //limits ampl pion
 706   fstep1->SetParLimits(1,fRangeStep4[0],fRangeStep4[1]);                                       //limits mean pion (dummy)
 707   fstep1->SetParLimits(2,fExpPionSigma*fLimitsOnSigmaPion[0],fExpPionSigma*fLimitsOnSigmaPion[1]);   //limits sigma pion
 708
 709   if(fExpPionSigma>0)  h->Fit(fstep1,fFitOption.Data(),"",fExpPionMean+fRangeStep1[0],fExpPionMean+fRangeStep1[1]);//first fit
 710   else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);
 711
 712   printf("___________________________________________________________________ Second Step: proton\n");
 713   fstep2 = new TF1("step2",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);
 714   fstep2->SetParameter(0,fstep1->GetParameter(0)/16.);//initial ampl proton
 715   fstep2->SetParameter(1,fExpProtonMean);              //initial mean proton
 716   fstep2->SetParameter(2,fExpProtonSigma);               //initial sigma proton
 717   fstep2->SetParLimits(0,0.,fstep1->GetParameter(0));                                                //limits ampl proton
 718   fstep2->SetParLimits(1,fstep1->GetParameter(1),fRangeStep4[1]);                                    //limits mean proton
 719   fstep2->SetParLimits(2,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton
 720
 721   if(fExpProtonSigma>0) h->Fit(fstep2,fFitOption.Data(),"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//second fit
 722   else for(Int_t npar=0;npar<3;npar++) fstep2->FixParameter(npar,0.00001);
 723
 724   printf("___________________________________________________________________ Third Step: kaon\n");
 725   fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);
 726   fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion
 727   fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion
 728   fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion
 729   fstep3->FixParameter(6,fstep2->GetParameter(0)); //fixed ampl proton
 730   fstep3->FixParameter(7,fstep2->GetParameter(1)); //fixed mean proton
 731   fstep3->FixParameter(8,fstep2->GetParameter(2)); //fixed sigma proton
 732   fstep3->SetParameter(3,fstep1->GetParameter(0)/8.); //initial ampl kaon
 733   fstep3->SetParameter(4,fExpKaonMean);                //initial mean kaon
 734   fstep3->SetParameter(5,fExpKaonSigma);                 //initial sigma kaon
 735   fstep3->SetParLimits(3,fstep2->GetParameter(0),fstep1->GetParameter(0));                   //limits ampl kaon
 736   fstep3->SetParLimits(4,fstep1->GetParameter(1),fstep2->GetParameter(1));                   //limits mean kaon
 737   fstep3->SetParLimits(5,fExpKaonSigma*fLimitsOnSigmaKaon[0],fExpKaonSigma*fLimitsOnSigmaKaon[1]); //limits sigma kaon
 738   if(fExpKaonSigma>0) h->Fit(fstep3,fFitOption.Data(),"",fExpKaonMean+fRangeStep2[0],fExpKaonMean+fRangeStep2[1]);//third fit
 739   else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);
 740
 741   printf("___________________________________________________________________ Final Step: refit all\n");
 742   fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);
 743   fstepTot->SetLineColor(1);
 744
 745   Double_t initialParametersStepTot[9];
 746   initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian
 747   initialParametersStepTot[1] = fstep1->GetParameter(1);
 748   initialParametersStepTot[2] = fstep1->GetParameter(2);
 749
 750   initialParametersStepTot[3] = fstep3->GetParameter(3);//second gaussian
 751   initialParametersStepTot[4] = fstep3->GetParameter(4);
 752   initialParametersStepTot[5] = fstep3->GetParameter(5);
 753
 754   initialParametersStepTot[6] = fstep2->GetParameter(0);//third gaussian
 755   initialParametersStepTot[7] = fstep2->GetParameter(1);
 756   initialParametersStepTot[8] = fstep2->GetParameter(2);
 757
 758   fstepTot->SetParameters(initialParametersStepTot); //initial parameters
 759   fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl pion
 760   fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion
 761   fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion
 762   fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon
 763   fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon
 764   fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon
 765   fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl proton
 766   fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton
 767   fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton
 768
 769   h->Fit(fstepTot,fFitOption.Data(),"",fRangeStep4[0],fRangeStep4[1]);//refit all
 770
 771   for(Int_t j=0;j<9;j++) {
 772     fFitPars[j] = fstepTot->GetParameter(j);
 773     fFitParErrors[j] = fstepTot->GetParError(j);
 774   }
 775   //************************************* storing parameter to calculate the yields *******
 776   Int_t chpa=0;
 777   if(code==321) chpa=3;
 778   if(code==2212) chpa=6;
 779   for(Int_t j=0;j<3;j++) {
 780     fFitpar[j] = fstepTot->GetParameter(j+chpa);
 781     fFitparErr[j] = fstepTot->GetParError(j+chpa);
 782   }
 783
 784   DrawFitFunction(fstepTot);
 785   CalcResidual(h,fstepTot,gres);
 786   return;
 787 }
 788
 789 //________________________________________________________
 790 void AliITSsadEdxFitter::DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode, TGraph *gres){
 791   // 3-gaussian fit to log(dedx)-log(dedxBB) histogram
 792   // pt bin from 0 to 20, code={211,321,2212}
 793   // first step: proton peak, second step: pion peak, third step: kaon peak
 794   // final step: refit all using the parameters
 795   TF1 *fstep1, *fstep2, *fstep3, *fstepTot;
 796   Int_t code=TMath::Abs(signedcode);
 797   Initialize(h,code,bin);
 798   PrintInitialValues();
 799   if(!IsGoodBin(bin,code)) return;
 800
 801   printf("___________________________________________________________________ First Step: proton\n");
 802   fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);
 803   fstep1->SetParameter(0,fExpPionAmpl/16.);       //initial ampl proton`
 804   fstep1->SetParameter(1,fExpProtonMean);  //initial mean proton
 805   fstep1->SetParameter(2,fExpProtonSigma);   //initial sigma proton
 806   fstep1->SetParLimits(0,0,fExpPionAmpl);                                                                    //limits ampl proton
 807   fstep1->SetParLimits(1,fExpPionMean,fRangeStep4[1]);                                                       //limits mean proton (dummy)
 808   fstep1->SetParLimits(2,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton
 809
 810   if(fExpProtonSigma>0)  h->Fit(fstep1,fFitOption.Data(),"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//first fit
 811   else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);
 812
 813   printf("___________________________________________________________________ Second Step: pion\n");
 814   fstep2 = new TF1("step2",DoubleGausStep,fRangeStep4[0],fRangeStep4[1],6);
 815   fstep2->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl proton
 816   fstep2->FixParameter(1,fstep1->GetParameter(1)); //fixed mean proton
 817   fstep2->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma proton
 818   fstep2->SetParameter(3,fExpPionAmpl);             //initial ampl pion
 819   fstep2->SetParameter(4,fExpPionMean);             //initial mean pion
 820   fstep2->SetParameter(5,fExpPionSigma);       //initial sigma pion
 821   fstep2->SetParLimits(3,0.,1.2*fExpPionAmpl);                                                               //limits ampl pion
 822   fstep2->SetParLimits(4,fRangeStep4[0],fstep1->GetParameter(1));                                //limits mean pion
 823   fstep2->SetParLimits(5,fExpPionSigma*fLimitsOnSigmaPion[0],fExpPionSigma*fLimitsOnSigmaPion[1]);     //limits sigma pion
 824
 825   if(fExpPionSigma>0) h->Fit(fstep2,fFitOption.Data(),"",fExpPionMean+fRangeStep1[0],fExpPionMean+fRangeStep1[1]);//second fit
 826   else for(Int_t npar=0;npar<3;npar++) fstep2->FixParameter(npar,0.00001);
 827
 828   printf("___________________________________________________________________ Third Step: kaon\n");
 829   fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);
 830   fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl proton
 831   fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean proton
 832   fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma proton
 833   fstep3->FixParameter(3,fstep2->GetParameter(3)); //fixed ampl pion
 834   fstep3->FixParameter(4,fstep2->GetParameter(4)); //fixed mean pion
 835   fstep3->FixParameter(5,fstep2->GetParameter(5)); //fixed sigma pion
 836   fstep3->SetParameter(6,fstep2->GetParameter(0)/8.); //initial ampl kaon
 837   fstep3->SetParameter(7,fExpKaonMean);                //initial mean kaon
 838   fstep3->SetParameter(8,fExpKaonSigma);                 //initial sigma kaon
 839   fstep3->SetParLimits(6,fstep1->GetParameter(0),fstep2->GetParameter(3));                   //limits ampl kaon
 840   fstep3->SetParLimits(7,fstep2->GetParameter(4),fstep1->GetParameter(1));                   //limits mean kaon
 841   fstep3->SetParLimits(8,fExpKaonSigma*fLimitsOnSigmaKaon[0],fExpKaonSigma*fLimitsOnSigmaKaon[1]); //limits sigma kaon
 842   if(fExpKaonSigma>0) h->Fit(fstep3,fFitOption.Data(),"",fExpKaonMean+fRangeStep2[0],fExpKaonMean+fRangeStep2[1]);//third fit
 843   else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);
 844
 845   printf("___________________________________________________________________ Final Step: refit all\n");
 846   fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);
 847   fstepTot->SetLineColor(1);
 848
 849   Double_t initialParametersStepTot[9];
 850   initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian
 851   initialParametersStepTot[1] = fstep1->GetParameter(1);
 852   initialParametersStepTot[2] = fstep1->GetParameter(2);
 853
 854   initialParametersStepTot[3] = fstep3->GetParameter(6);//second gaussian
 855   initialParametersStepTot[4] = fstep3->GetParameter(7);
 856   initialParametersStepTot[5] = fstep3->GetParameter(8);
 857
 858   initialParametersStepTot[6] = fstep2->GetParameter(3);//third gaussian
 859   initialParametersStepTot[7] = fstep2->GetParameter(4);
 860   initialParametersStepTot[8] = fstep2->GetParameter(5);
 861
 862   fstepTot->SetParameters(initialParametersStepTot); //initial parameters
 863   fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl proton
 864   fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion
 865   fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion
 866   fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon
 867   fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon
 868   fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon
 869   fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl pion
 870   fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton
 871   fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton
 872
 873   h->Fit(fstepTot,fFitOption.Data(),"",fRangeStep4[0],fRangeStep4[1]);//refit all
 874
 875   for(Int_t j=0;j<9;j++) {
 876     fFitPars[j] = fstepTot->GetParameter(j);
 877     fFitParErrors[j] = fstepTot->GetParError(j);
 878   }
 879   //************************************* storing parameter to calculate the yields *******
 880   Int_t chpa=0;
 881   if(code==321) chpa=3;
 882   if(code==211) chpa=6;
 883   for(Int_t j=0;j<3;j++) {
 884     fFitpar[j] = fstepTot->GetParameter(j+chpa);
 885     fFitparErr[j] = fstepTot->GetParError(j+chpa);
 886   }
 887
 888   DrawFitFunction(fstepTot);
 889   CalcResidual(h,fstepTot,gres);
 890   return;
 891 }
 892
 893
 894 //________________________________________________________
 895 void AliITSsadEdxFitter::DoFitOnePeak(TH1F *h, Int_t bin, Int_t signedcode){
 896   // single-gaussian fit to log(dedx)-log(dedxBB) histogram
 897   TF1 *fstep1;
 898   Int_t code=TMath::Abs(signedcode);
 899   Initialize(h,code,bin);
 900   PrintInitialValues();
 901   if(!IsGoodBin(bin,code)) return;
 902
 903   printf("___________________________________________________________________ Single Step\n");
 904   fstep1 = new TF1("step2",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);
 905   fstep1->SetParameter(0,fExpPionAmpl/16.);                   //initial ampl
 906   fstep1->SetParameter(1,fExpProtonMean);              //initial mean
 907   fstep1->SetParameter(2,fExpProtonSigma);               //initial sigma
 908   fstep1->SetParLimits(0,0.,fExpPionAmpl);                                                                   //limits ampl proton
 909   fstep1->SetParLimits(1,fExpPionMean,fRangeStep4[1]);                                                       //limits mean proton
 910   //fstep1->SetParLimits(2,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton
 911
 912   if(fExpProtonSigma>0) h->Fit(fstep1,fFitOption.Data(),"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//fit
 913   else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);
 914
 915   fstep1->SetLineColor(1);
 916   fstep1->Draw("same");
 917
 918   fFitpar[0] = fstep1->GetParameter(0);
 919   fFitparErr[0] = fstep1->GetParError(0);
 920   return;
 921 }
 922
 923 //________________________________________________________
 924 void AliITSsadEdxFitter::DoFitTail(TH1F *h, Int_t bin, Int_t signedcode){
 925   // 3-gaussian fit to log(dedx)-log(dedxBB) histogram
 926   // pt bin from 0 to 20, code={211,321,2212}
 927   // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed
 928   // final step: refit all using the parameters and tollerance limits (+-20%)
 929   // WARNING: exponential tail added in the right of the Gaussian shape
 930   Int_t code=TMath::Abs(signedcode);
 931   if(!IsGoodBin(bin,code)) return;
 932
 933   TF1 *fstep1, *fstep2, *fstep3, *fstepTot;
 934   Initialize(h,code,bin);
 935   PrintInitialValues();
 936
 937   printf("\n___________________________________________________________________\n First Step: pions\n\n");
 938   fstep1 = new TF1("step1",SingleGausTail,-3.5,3.5,5);
 939   fstep1->SetParameter(0,fExpPionAmpl);//initial
 940   fstep1->SetParameter(1,fExpPionMean);
 941   fstep1->SetParameter(3,1.2);
 942   fstep1->SetParameter(4,10.);
 943
 944   fstep1->SetParLimits(0,0,fExpPionAmpl*1.2);
 945   fstep1->SetParLimits(1,-3.5,3.5);
 946   fstep1->SetParLimits(2,0.1,0.25);
 947   fstep1->SetParLimits(4,5.,20.);
 948   if(bin<8) fstep1->SetParLimits(4,13.,25.);
 949
 950   h->Fit(fstep1,fFitOption.Data(),"",fExpPionMean-0.45,fExpPionMean+0.45);//first fit
 951
 952   printf("\n___________________________________________________________________\n Second Step: kaons\n\n");
 953   fstep2 = new TF1("fstep2",DoubleGausTail,-3.5,3.5,10);
 954   fstep2->FixParameter(0,fstep1->GetParameter(0));//fixed
 955   fstep2->FixParameter(1,fstep1->GetParameter(1));
 956   fstep2->FixParameter(2,fstep1->GetParameter(2));
 957   fstep2->FixParameter(3,fstep1->GetParameter(3));
 958   fstep2->FixParameter(4,fstep1->GetParameter(4));
 959
 960   fstep2->SetParameter(5,fstep1->GetParameter(0)/8);//initial
 961   //fstep2->SetParameter(6,CalcP(code,322,bin));
 962   fstep2->SetParameter(7,0.145);
 963   fstep2->FixParameter(8,1.2);
 964   fstep2->SetParameter(9,13.);
 965
 966   fstep2->SetParLimits(5,fstep1->GetParameter(0)/100,fstep1->GetParameter(0));//limits
 967   fstep2->SetParLimits(6,-3.5,3.5);
 968   fstep2->SetParLimits(7,0.12,0.2);
 969   fstep2->SetParLimits(9,9.,20.);
 970   if(bin<9) fstep2->SetParLimits(9,13.,25.);
 971
 972   if(bin<6 || bin>12) for(Int_t npar=5;npar<10;npar++) fstep2->FixParameter(npar,-0.0000000001);
 973
 974   printf("\n____________________________________________________________________\n Third Step: protons \n\n");
 975   fstep3= new TF1("fstep3",FinalGausTail,-3.5,3.5,15);
 976   fstep3->FixParameter(0,fstep1->GetParameter(0));//fixed
 977   fstep3->FixParameter(1,fstep1->GetParameter(1));
 978   fstep3->FixParameter(2,fstep1->GetParameter(2));
 979   fstep3->FixParameter(3,fstep1->GetParameter(3));
 980   fstep3->FixParameter(4,fstep1->GetParameter(4));
 981   fstep3->FixParameter(5,fstep2->GetParameter(5));
 982   fstep3->FixParameter(6,fstep2->GetParameter(6));
 983   fstep3->FixParameter(7,fstep2->GetParameter(7));
 984   fstep3->FixParameter(8,fstep2->GetParameter(8));
 985   fstep3->FixParameter(9,fstep2->GetParameter(9));
 986
 987   fstep3->SetParameter(10,fstep2->GetParameter(0)/8);//initial
 988   //fstep3->SetParameter(11,CalcP(code,2212,bin));
 989   fstep3->SetParameter(12,0.145);
 990   fstep3->FixParameter(13,1.2);
 991   fstep3->SetParameter(14,10.);
 992
 993   fstep3->SetParLimits(10,fstep2->GetParameter(0)/100,fstep2->GetParameter(0));//limits
 994   fstep3->SetParLimits(11,-3.5,3.5);
 995   fstep3->SetParLimits(12,0.12,0.2);
 996   fstep3->SetParLimits(14,11.,25.);
 997
 998   printf("\n_____________________________________________________________________\n Final Step: refit all \n\n");
 999   fstepTot = new TF1("funztot",FinalGausTail,-3.5,3.5,15);
1000   fstepTot->SetLineColor(1);
1001
1002   Double_t initialParametersStepTot[15];
1003   initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian
1004   initialParametersStepTot[1] = fstep1->GetParameter(1);
1005   initialParametersStepTot[2] = fstep1->GetParameter(2);
1006   initialParametersStepTot[3] = fstep1->GetParameter(3);
1007   initialParametersStepTot[4] = fstep1->GetParameter(4);
1008
1009   initialParametersStepTot[5] = fstep2->GetParameter(5);//second gaussian
1010   initialParametersStepTot[6] = fstep2->GetParameter(6);
1011   initialParametersStepTot[7] = fstep2->GetParameter(7);
1012   initialParametersStepTot[8] = fstep2->GetParameter(8);
1013   initialParametersStepTot[9] = fstep2->GetParameter(9);
1014
1015   initialParametersStepTot[10] = fstep3->GetParameter(10);//third gaussian
1016   initialParametersStepTot[11] = fstep3->GetParameter(11);
1017   initialParametersStepTot[12] = fstep3->GetParameter(12);
1018   initialParametersStepTot[13] = fstep3->GetParameter(13);
1019   initialParametersStepTot[14] = fstep3->GetParameter(14);
1020
1021   fstepTot->SetParameters(initialParametersStepTot);//initial parameter
1022
1023
1024   fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1);//tollerance limit
1025   fstepTot->SetParLimits(1,initialParametersStepTot[1]*0.9,initialParametersStepTot[1]*1.1);
1026   fstepTot->SetParLimits(2,initialParametersStepTot[2]*0.9,initialParametersStepTot[2]*1.1);
1027   fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1);
1028   fstepTot->SetParLimits(4,initialParametersStepTot[4]*0.9,initialParametersStepTot[4]*1.1);
1029   fstepTot->SetParLimits(5,initialParametersStepTot[5]*0.9,initialParametersStepTot[5]*1.1);
1030   fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1);
1031   fstepTot->SetParLimits(7,initialParametersStepTot[7]*0.9,initialParametersStepTot[7]*1.1);
1032   fstepTot->SetParLimits(8,initialParametersStepTot[8]*0.9,initialParametersStepTot[8]*1.1);
1033   fstepTot->SetParLimits(9,initialParametersStepTot[9]*0.9,initialParametersStepTot[9]*1.1);
1034   fstepTot->SetParLimits(10,initialParametersStepTot[10]*0.9,initialParametersStepTot[10]*1.1);
1035   fstepTot->SetParLimits(11,initialParametersStepTot[11]*0.9,initialParametersStepTot[11]*1.1);
1036   fstepTot->SetParLimits(12,initialParametersStepTot[12]*0.9,initialParametersStepTot[12]*1.1);
1037   fstepTot->SetParLimits(13,initialParametersStepTot[13]*0.9,initialParametersStepTot[13]*1.1);
1038   fstepTot->SetParLimits(14,initialParametersStepTot[14]*0.9,initialParametersStepTot[14]*1.1);
1039
1040   if(bin<9) for(Int_t npar=10;npar<15;npar++) fstepTot->FixParameter(npar,-0.00000000001);
1041   h->Fit(fstepTot,fFitOption.Data(),"",-3.5,3.5); //refit all
1042
1043   for(Int_t j=0;j<9;j++) {
1044     fFitPars[j] = fstepTot->GetParameter(j);
1045     fFitParErrors[j] = fstepTot->GetParError(j);
1046   }
1047
1048   //************************************* storing parameter to calculate the yields *******
1049   Int_t chpa=0;
1050   if(code==321) chpa=5;
1051   if(code==2212) chpa=10;
1052   for(Int_t j=0;j<3;j++) {
1053     fFitpar[j] = fstepTot->GetParameter(j+chpa);
1054     fFitparErr[j] = fstepTot->GetParError(j+chpa);
1055   }
1056
1057   DrawFitFunction(fstepTot);
1058   return;
1059 }
1060
1061 //________________________________________________________
1062 void AliITSsadEdxFitter::FillHisto(TH1F *hsps, Int_t bin, Float_t binsize, Int_t code){
1063   // fill the spectra histo calculating the yield
1064   // first bin has to be 1
1065   Double_t yield = 0;
1066   Double_t err = 0;
1067   Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin);
1068
1069   if(IsGoodBin(bin-1,code)) {
1070     yield = fFitpar[0] / ptbin / binsize;
1071     err = fFitparErr[0] / ptbin / binsize;
1072   }
1073
1074   hsps->SetBinContent(bin,yield);
1075   hsps->SetBinError(bin,err);
1076   return;
1077 }
1078
1079 //________________________________________________________
1080 void AliITSsadEdxFitter::FillHistoMC(TH1F *hsps, Int_t bin, Int_t code, TH1F *h){
1081   // fill the spectra histo calculating the yield (using the MC truth)
1082   // first bin has to be 1
1083   Double_t yield = 0;
1084   Double_t erryield=0;
1085   Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin);
1086
1087   if(IsGoodBin(bin-1,code)){
1088     yield = h->GetEntries() / ptbin;
1089     erryield=TMath::Sqrt(h->GetEntries()) / ptbin;
1090   }
1091   hsps->SetBinContent(bin,yield);
1092   hsps->SetBinError(bin,erryield);
1093   return;
1094 }
1095
1096 //________________________________________________________
1097 void AliITSsadEdxFitter::GetFitPar(Double_t *fitpar, Double_t *fitparerr) const {
1098   // getter of the fit parameters and the relative errors
1099   for(Int_t i=0;i<3;i++) {
1100     fitpar[i] = fFitpar[i];
1101     fitparerr[i] = fFitparErr[i];
1102   }
1103   return;
1104 }
1105
1106
1107 //________________________________________________________
1108 void AliITSsadEdxFitter::PrintAll() const{
1109   // print parameters
1110
1111   printf("Range 1 = %f %f\n",fRangeStep1[0],fRangeStep1[1]);
1112   printf("Range 2 = %f %f\n",fRangeStep2[0],fRangeStep2[1]);
1113   printf("Range 3 = %f %f\n",fRangeStep3[0],fRangeStep3[1]);
1114   printf("Range F = %f %f\n",fRangeStep4[0],fRangeStep4[1]);
1115   printf(" Sigma1 = %f %f\n",fLimitsOnSigmaPion[0],fLimitsOnSigmaPion[1]);
1116   printf(" Sigma2 = %f %f\n",fLimitsOnSigmaKaon[0],fLimitsOnSigmaKaon[1]);
1117   printf(" Sigma3 = %f %f\n",fLimitsOnSigmaProton[0],fLimitsOnSigmaProton[1]);
1118 }
1119 //______________________________________________________________________
1120 void AliITSsadEdxFitter::PrintInitialValues() const{
1121   // print initial values
1122
1123   printf("mean values        -> %f %f %f\n",fExpPionMean,fExpKaonMean,fExpProtonMean);
1124   printf("integration ranges -> (%1.2f,%1.2f) (%1.2f,%1.2f) (%1.2f,%1.2f)\n",
1125          fRangeStep1[0],fRangeStep1[1],fRangeStep2[0],fRangeStep2[1],fRangeStep3[0],fRangeStep3[1]);
1126   printf("sigma values -> %f %f %f\n",fExpPionSigma,fExpKaonSigma,fExpProtonSigma);
1127   printf("sigma ranges -> (%1.2f,%1.2f) (%1.2f,%1.2f) (%1.2f,%1.2f)\n",
1128          fLimitsOnSigmaPion[0],fLimitsOnSigmaPion[1],
1129          fLimitsOnSigmaKaon[0],fLimitsOnSigmaKaon[1],
1130          fLimitsOnSigmaProton[0],fLimitsOnSigmaProton[1]);
1131 }