Updated parameterizations of peak positions (R. Russo, F. Prino)

diff --git a/PWG2/SPECTRA/AliITSsadEdxFitter.cxx b/PWG2/SPECTRA/AliITSsadEdxFitter.cxx
index cd71681f36ab3cac184eca8a80859008531826cf..cc673695fc686525af14a3670e6a7e7600b3245d 100644 (file)
--- a/PWG2/SPECTRA/AliITSsadEdxFitter.cxx
+++ b/PWG2/SPECTRA/AliITSsadEdxFitter.cxx
@@ -41,7 +41,16 @@
 \r
 ClassImp(AliITSsadEdxFitter)\r
 //______________________________________________________________________\r
-AliITSsadEdxFitter::AliITSsadEdxFitter():TObject(){\r
+AliITSsadEdxFitter::AliITSsadEdxFitter():TObject(),\r
+  fExpPionMean(0.),\r
+  fExpKaonMean(0.),\r
+  fExpProtonMean(0.),\r
+  fExpPionSigma(0.),\r
+  fExpKaonSigma(0.),\r
+  fExpProtonSigma(0.),\r
+  fIsMC(kFALSE),\r
+  fITSpid(0)\r
+{\r
   // standard constructor\r
   for(Int_t i=0; i<3; i++)  fFitpar[i] = 0.;\r
   for(Int_t i=0; i<3; i++)  fFitparErr[i] = 0.;\r
@@ -55,15 +64,42 @@ AliITSsadEdxFitter::AliITSsadEdxFitter():TObject(){
   SetBinsUsedPion();\r
   SetBinsUsedKaon();\r
   SetBinsUsedProton();\r
+  fITSpid=new AliITSPIDResponse(kFALSE);\r
+};\r
+//______________________________________________________________________\r
+AliITSsadEdxFitter::AliITSsadEdxFitter(Bool_t isMC):TObject(),\r
+  fExpPionMean(0.),\r
+  fExpKaonMean(0.),\r
+  fExpProtonMean(0.),\r
+  fExpPionSigma(0.),\r
+  fExpKaonSigma(0.),\r
+  fExpProtonSigma(0.),\r
+  fIsMC(isMC),\r
+  fITSpid(0)\r
+{\r
+  // standard constructor\r
+  for(Int_t i=0; i<3; i++)  fFitpar[i] = 0.;\r
+  for(Int_t i=0; i<3; i++)  fFitparErr[i] = 0.;\r
+  SetRangeStep1();\r
+  SetRangeStep2();\r
+  SetRangeStep3();\r
+  SetRangeFinalStep();\r
+  SetLimitsOnSigmaPion();\r
+  SetLimitsOnSigmaKaon();\r
+  SetLimitsOnSigmaProton();\r
+  SetBinsUsedPion();\r
+  SetBinsUsedKaon();\r
+  SetBinsUsedProton();\r
+  fITSpid=new AliITSPIDResponse(isMC);\r
 };\r
 \r
 //________________________________________________________\r
-Double_t AliITSsadEdxFitter::CalcSigma(Int_t code,Float_t x,Bool_t mc){\r
+Double_t AliITSsadEdxFitter::CalcSigma(Int_t code,Float_t x) const {\r
   // calculate the sigma 12-ott-2010  \r
   Double_t p[2]={0.};\r
   Double_t xMinKaon=0.15; //minimum pt value to consider the kaon peak\r
   Double_t xMinProton=0.3;//minimum pt value to consider the proton peak\r
-  if(mc){\r
+  if(fIsMC){\r
     if(code==211){\r
       p[0] = -1.20337e-04;\r
       p[1] = 1.13060e-01;\r
@@ -97,7 +133,7 @@ Double_t AliITSsadEdxFitter::CalcSigma(Int_t code,Float_t x,Bool_t mc){
 }\r
 \r
 //_______________________________________________________\r
-Int_t AliITSsadEdxFitter::CalcMean(Int_t code, Float_t x, Float_t mean0, Float_t &mean1, Float_t &mean2){\r
+Int_t AliITSsadEdxFitter::CalcMean(Int_t code, Float_t x, Float_t mean0, Float_t &mean1, Float_t &mean2) const{\r
   // calculate the mean 12-ott-2010  \r
   Double_t p1[4]={0.};\r
   Double_t p2[4]={0.};\r
@@ -374,7 +410,7 @@ void AliITSsadEdxFitter::DrawFitFunction(TF1 *fun) const {
 }\r
 \r
 //______________________________________________________________________\r
-void AliITSsadEdxFitter::GetInitialParam(TH1F* h,Bool_t mc,Int_t code,Int_t bin, Float_t &pt, Float_t &ampl, Float_t &mean1, Float_t &mean2, Float_t &mean3, Float_t &sigma1, Float_t &sigma2, Float_t &sigma3){\r
+void AliITSsadEdxFitter::GetInitialParam(TH1F* h, Int_t code,Int_t bin, Float_t &pt, Float_t &ampl){\r
   //code to get the expected values to use for fitting\r
   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};\r
   pt=(xbins[bin+1]+xbins[bin])/2;\r
@@ -387,47 +423,59 @@ void AliITSsadEdxFitter::GetInitialParam(TH1F* h,Bool_t mc,Int_t code,Int_t bin,
   h->SetFillColor(11);\r
        \r
   //expected values\r
-  Int_t xmax=-1,ymax=-1,zmax=-1;\r
-  h->GetMaximumBin(xmax,ymax,zmax);\r
-  printf("\n---------------------------------- BIN %d - hypothesis %d ----------------------------------\n",bin,code);\r
   Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
   ampl = h->GetMaximum()/(h->GetRMS()*s2pi);\r
-  mean1 = h->GetBinLowEdge(xmax); //expected mean values\r
-  Int_t calcmean=CalcMean(code,pt,mean1,mean2,mean3);\r
-  if(calcmean<0) cout<<"Error during mean calculation"<<endl;\r
-  printf("mean values        -> %f %f %f\n",mean1,mean2,mean3);\r
+\r
+  Double_t massp=AliPID::ParticleMass(AliPID::kProton);\r
+  Double_t massk=AliPID::ParticleMass(AliPID::kKaon);\r
+  Double_t masspi=AliPID::ParticleMass(AliPID::kPion);\r
+  Bool_t isSA=kTRUE;\r
+  Float_t sinthmed=0.8878;\r
+  Float_t p=pt/sinthmed;\r
+  Double_t bethep=fITSpid->Bethe(p,massp,isSA);\r
+  Double_t bethek=fITSpid->Bethe(p,massk,isSA);\r
+  Double_t bethepi=fITSpid->Bethe(p,masspi,isSA);\r
+  Double_t betheref=bethepi;\r
+  if(TMath::Abs(code)==321) betheref=bethek;\r
+  else if(TMath::Abs(code)==2212) betheref=bethep;\r
+  fExpPionMean=TMath::Log(bethepi)-TMath::Log(betheref);\r
+  fExpKaonMean=TMath::Log(bethek)-TMath::Log(betheref);\r
+  fExpProtonMean=TMath::Log(bethep)-TMath::Log(betheref);\r
+\r
+  printf("mean values        -> %f %f %f\n",fExpPionMean,fExpKaonMean,fExpProtonMean);\r
   printf("integration ranges -> (%1.2f,%1.2f) (%1.2f,%1.2f) (%1.2f,%1.2f)\n",fRangeStep1[0],fRangeStep1[1],fRangeStep2[0],fRangeStep2[1],fRangeStep3[0],fRangeStep3[1]);\r
-  sigma1 = CalcSigma(211,pt,mc); //expected sigma values\r
-  sigma2 = CalcSigma(321,pt,mc);\r
-  sigma3 = CalcSigma(2212,pt,mc);\r
-  printf("sigma values -> %f %f %f\n",sigma1,sigma2,sigma3);\r
+  fExpPionSigma = CalcSigma(211,pt); //expected sigma values\r
+  fExpKaonSigma = CalcSigma(321,pt);\r
+  fExpProtonSigma = CalcSigma(2212,pt);\r
+  printf("sigma values -> %f %f %f\n",fExpPionSigma,fExpKaonSigma,fExpProtonSigma);\r
   printf("sigma ranges -> (%1.2f,%1.2f) (%1.2f,%1.2f) (%1.2f,%1.2f)\n",fLimitsOnSigmaPion[0],fLimitsOnSigmaPion[1],fLimitsOnSigmaKaon[0],fLimitsOnSigmaKaon[1],fLimitsOnSigmaProton[0],fLimitsOnSigmaProton[1]);\r
   return;\r
 }\r
 \r
 //________________________________________________________\r
-void AliITSsadEdxFitter::DoFit(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, TGraph *gres){\r
+void AliITSsadEdxFitter::DoFit(TH1F *h, Int_t bin, Int_t signedcode, TGraph *gres){\r
   // 3-gaussian fit to log(dedx)-log(dedxBB) histogram\r
   // pt bin from 0 to 20, code={211,321,2212} \r
   // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed\r
   // final step: refit all using the parameters and tollerance limits (+-20%)\r
+\r
   TF1 *fstep1, *fstep2, *fstep3, *fstepTot;\r
   TString modfit = "M0R+";\r
-  Float_t pt=0., ampl=0., mean=0., expKaonMean=0., expProtonMean=0., expPionSig=0., expKaonSig=0., expProtonSig=0.;\r
+  Float_t pt=0., ampl=0.;\r
   Int_t code=TMath::Abs(signedcode);\r
-  GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\r
+  GetInitialParam(h,code,bin,pt,ampl);\r
   if(!IsGoodBin(bin,code)) return;\r
 \r
   printf("___________________________________________________________________ First Step: pions\n");\r
   fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);\r
   fstep1->SetParameter(0,ampl);       //initial ampl pion\r
-  fstep1->SetParameter(1,mean);       //initial mean pion\r
-  fstep1->SetParameter(2,expPionSig); //initial sigma pion\r
-  fstep1->SetParLimits(0,0.,ampl*1.2);                                                       //limits ampl pion\r
+  fstep1->SetParameter(1,fExpPionMean);       //initial mean pion\r
+  fstep1->SetParameter(2,fExpPionSigma); //initial sigma pion\r
+  fstep1->SetParLimits(0,0.,ampl*1.2);    //limits ampl pion\r
   fstep1->SetParLimits(1,fRangeStep4[0],fRangeStep4[1]);                                     //limits mean pion (dummy)\r
-  fstep1->SetParLimits(2,expPionSig*fLimitsOnSigmaPion[0],expPionSig*fLimitsOnSigmaPion[1]); //limits sigma pion\r
+  fstep1->SetParLimits(2,fExpPionSigma*fLimitsOnSigmaPion[0],fExpPionSigma*fLimitsOnSigmaPion[1]); //limits sigma pion\r
 \r
-  if(expPionSig>0) h->Fit(fstep1,modfit.Data(),"",mean+fRangeStep1[0],mean+fRangeStep1[1]);//first fit\r
+  if(fExpPionSigma>0) h->Fit(fstep1,modfit.Data(),"",fExpPionMean+fRangeStep1[0],fExpPionMean+fRangeStep1[1]);//first fit\r
   else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);\r
 \r
   printf("___________________________________________________________________ Second Step: kaons\n");\r
@@ -436,24 +484,15 @@ void AliITSsadEdxFitter::DoFit(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc,
   fstep2->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion\r
   fstep2->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion\r
   fstep2->SetParameter(3,fstep1->GetParameter(0)/8.); //initial ampl kaon\r
-  fstep2->SetParameter(4,expKaonMean);                //initial mean kaon\r
-  fstep2->SetParameter(3,expKaonSig);                 //initial sigma kaon\r
+  fstep2->SetParameter(4,fExpKaonMean);                //initial mean kaon\r
+  fstep2->SetParameter(3,fExpKaonSigma);                 //initial sigma kaon\r
   fstep2->SetParLimits(3,0.,fstep1->GetParameter(0));                                         //limits ampl kaon\r
   fstep2->SetParLimits(4,fstep1->GetParameter(1),fRangeStep4[1]);                             //limits mean kaon \r
-  fstep2->SetParLimits(5,expKaonSig*fLimitsOnSigmaKaon[0],expKaonSig*fLimitsOnSigmaKaon[1]);  //limits sigma kaon\r
+  fstep2->SetParLimits(5,fExpKaonSigma*fLimitsOnSigmaKaon[0],fExpKaonSigma*fLimitsOnSigmaKaon[1]);  //limits sigma kaon\r
 \r
-  if(expKaonSig>0) h->Fit(fstep2,modfit.Data(),"",expKaonMean+fRangeStep2[0],expKaonMean+fRangeStep2[1]);//second fit\r
+  if(fExpKaonSigma>0) h->Fit(fstep2,modfit.Data(),"",fExpKaonMean+fRangeStep2[0],fExpKaonMean+fRangeStep2[1]);//second fit\r
   else for(Int_t npar=3;npar<6;npar++) fstep2->FixParameter(npar,0.00001);\r
 \r
-  /*TLine *l[3];\r
-    l[0] = new TLine(expKaonMean,0,expKaonMean,10000);\r
-    l[1] = new TLine(expKaonMean+fRangeStep2[0],0,expKaonMean+fRangeStep2[0],10000);\r
-    l[2] = new TLine(expKaonMean+fRangeStep2[1],0,expKaonMean+fRangeStep2[1],10000);\r
-    for(Int_t dp=0;dp<3;dp++) {\r
-    l[dp]->Draw("same");\r
-    l[dp]->SetLineColor(4);\r
-    l[dp]->SetLineWidth(3);\r
-    }*/\r
 \r
   printf("___________________________________________________________________ Third Step: protons\n");\r
   fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
@@ -464,13 +503,13 @@ void AliITSsadEdxFitter::DoFit(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc,
   fstep3->FixParameter(4,fstep2->GetParameter(4)); //fixed mean kaon\r
   fstep3->FixParameter(5,fstep2->GetParameter(5)); //fidex sigma kaon\r
   fstep3->SetParameter(6,fstep2->GetParameter(0)/16.); //initial ampl proton\r
-  fstep3->SetParameter(7,expProtonMean);               //initial mean proton\r
-  fstep3->SetParameter(8,expProtonSig);                //initial sigma proton\r
+  fstep3->SetParameter(7,fExpProtonMean);               //initial mean proton\r
+  fstep3->SetParameter(8,fExpProtonSigma);                //initial sigma proton\r
   fstep3->SetParLimits(6,0.,fstep2->GetParameter(0));                                                //limits ampl proton\r
   fstep3->SetParLimits(7,fstep2->GetParameter(4),fRangeStep4[1]);                                    //limits mean proton\r
-  fstep3->SetParLimits(8,expProtonSig*fLimitsOnSigmaProton[0],expProtonSig*fLimitsOnSigmaProton[1]); //limits sigma proton\r
+  fstep3->SetParLimits(8,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton\r
 \r
-  if(expProtonSig>0) h->Fit(fstep3,modfit.Data(),"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);//third fit\r
+  if(fExpProtonSigma>0) h->Fit(fstep3,modfit.Data(),"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//third fit\r
   else for(Int_t npar=6;npar<9;npar++) fstep3->FixParameter(npar,0.00001);\r
 \r
   printf("___________________________________________________________________ Final Step: refit all\n");\r
@@ -518,40 +557,40 @@ void AliITSsadEdxFitter::DoFit(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc,
 }\r
 \r
 //________________________________________________________\r
-void AliITSsadEdxFitter::DoFitProton(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, TGraph *gres){\r
+void AliITSsadEdxFitter::DoFitProton(TH1F *h, Int_t bin, Int_t signedcode, TGraph *gres){\r
   // 3-gaussian fit to log(dedx)-log(dedxBB) histogram\r
   // pt bin from 0 to 20, code={211,321,2212} \r
   // first step: pion peak, second step: proton peak, third step: kaon peak\r
   // final step: refit all using the parameters\r
   TF1 *fstep1, *fstep2, *fstep3, *fstepTot;\r
   TString modfit = "M0R+";\r
-  Float_t pt=0., ampl=0., mean=0., expKaonMean=0., expProtonMean=0., expPionSig=0., expKaonSig=0., expProtonSig=0.;\r
+  Float_t pt=0., ampl=0.;\r
   Int_t code=TMath::Abs(signedcode);\r
-  GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\r
+  GetInitialParam(h,code,bin,pt,ampl);\r
   if(!IsGoodBin(bin,code)) return;\r
 \r
   printf("___________________________________________________________________ First Step: pion\n");\r
   fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);\r
   fstep1->SetParameter(0,ampl);       //initial ampl pion\r
-  fstep1->SetParameter(1,mean);       //initial mean pion\r
-  fstep1->SetParameter(2,expPionSig); //initial sigma pion\r
+  fstep1->SetParameter(1,fExpPionMean);       //initial mean pion\r
+  fstep1->SetParameter(2,fExpPionSigma); //initial sigma pion\r
   fstep1->SetParLimits(0,0,ampl*1.2);                                                          //limits ampl pion\r
   fstep1->SetParLimits(1,fRangeStep4[0],fRangeStep4[1]);                                       //limits mean pion (dummy)\r
-  fstep1->SetParLimits(2,expPionSig*fLimitsOnSigmaPion[0],expPionSig*fLimitsOnSigmaPion[1]);   //limits sigma pion\r
+  fstep1->SetParLimits(2,fExpPionSigma*fLimitsOnSigmaPion[0],fExpPionSigma*fLimitsOnSigmaPion[1]);   //limits sigma pion\r
 \r
-  if(expPionSig>0)  h->Fit(fstep1,modfit,"",mean+fRangeStep1[0],mean+fRangeStep1[1]);//first fit\r
+  if(fExpPionSigma>0)  h->Fit(fstep1,modfit,"",fExpPionMean+fRangeStep1[0],fExpPionMean+fRangeStep1[1]);//first fit\r
   else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);\r
 \r
   printf("___________________________________________________________________ Second Step: proton\n");\r
   fstep2 = new TF1("step2",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);\r
   fstep2->SetParameter(0,fstep1->GetParameter(0)/16.);//initial ampl proton\r
-  fstep2->SetParameter(1,expProtonMean);              //initial mean proton\r
-  fstep2->SetParameter(2,expProtonSig);               //initial sigma proton\r
+  fstep2->SetParameter(1,fExpProtonMean);              //initial mean proton\r
+  fstep2->SetParameter(2,fExpProtonSigma);               //initial sigma proton\r
   fstep2->SetParLimits(0,0.,fstep1->GetParameter(0));                                                //limits ampl proton\r
   fstep2->SetParLimits(1,fstep1->GetParameter(1),fRangeStep4[1]);                                    //limits mean proton\r
-  fstep2->SetParLimits(2,expProtonSig*fLimitsOnSigmaProton[0],expProtonSig*fLimitsOnSigmaProton[1]); //limits sigma proton\r
+  fstep2->SetParLimits(2,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton\r
 \r
-  if(expProtonSig>0) h->Fit(fstep2,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);//second fit\r
+  if(fExpProtonSigma>0) h->Fit(fstep2,modfit,"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//second fit\r
   else for(Int_t npar=0;npar<3;npar++) fstep2->FixParameter(npar,0.00001);\r
 \r
   printf("___________________________________________________________________ Third Step: kaon\n");\r
@@ -563,21 +602,12 @@ void AliITSsadEdxFitter::DoFitProton(TH1F *h, Int_t bin, Int_t signedcode, Bool_
   fstep3->FixParameter(7,fstep2->GetParameter(1)); //fixed mean proton\r
   fstep3->FixParameter(8,fstep2->GetParameter(2)); //fixed sigma proton\r
   fstep3->SetParameter(3,fstep1->GetParameter(0)/8.); //initial ampl kaon\r
-  fstep3->SetParameter(4,expKaonMean);                //initial mean kaon\r
-  fstep3->SetParameter(5,expKaonSig);                 //initial sigma kaon\r
+  fstep3->SetParameter(4,fExpKaonMean);                //initial mean kaon\r
+  fstep3->SetParameter(5,fExpKaonSigma);                 //initial sigma kaon\r
   fstep3->SetParLimits(3,fstep2->GetParameter(0),fstep1->GetParameter(0));                   //limits ampl kaon\r
   fstep3->SetParLimits(4,fstep1->GetParameter(1),fstep2->GetParameter(1));                   //limits mean kaon\r
-  fstep3->SetParLimits(5,expKaonSig*fLimitsOnSigmaKaon[0],expKaonSig*fLimitsOnSigmaKaon[1]); //limits sigma kaon\r
-  /*TLine *l[3];\r
-    l[0] = new TLine(expProtonMean,0,expProtonMean,10000);\r
-    l[1] = new TLine(expProtonMean+fRangeStep3[0],0,expProtonMean+fRangeStep3[0],10000);\r
-    l[2] = new TLine(expProtonMean+fRangeStep3[1],0,expProtonMean+fRangeStep3[1],10000);\r
-    for(Int_t dp=0;dp<3;dp++) {\r
-    l[dp]->Draw("same");\r
-    l[dp]->SetLineColor(2);\r
-    l[dp]->SetLineWidth(4);\r
-    }*/\r
-  if(expKaonSig>0) h->Fit(fstep3,modfit,"",expKaonMean+fRangeStep2[0],expKaonMean+fRangeStep2[1]);//third fit\r
+  fstep3->SetParLimits(5,fExpKaonSigma*fLimitsOnSigmaKaon[0],fExpKaonSigma*fLimitsOnSigmaKaon[1]); //limits sigma kaon\r
+  if(fExpKaonSigma>0) h->Fit(fstep3,modfit,"",fExpKaonMean+fRangeStep2[0],fExpKaonMean+fRangeStep2[1]);//third fit\r
   else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);\r
 \r
   printf("___________________________________________________________________ Final Step: refit all\n");\r
@@ -625,28 +655,28 @@ void AliITSsadEdxFitter::DoFitProton(TH1F *h, Int_t bin, Int_t signedcode, Bool_
 }\r
 \r
 //________________________________________________________\r
-void AliITSsadEdxFitter::DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, TGraph *gres){\r
+void AliITSsadEdxFitter::DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode, TGraph *gres){\r
   // 3-gaussian fit to log(dedx)-log(dedxBB) histogram\r
   // pt bin from 0 to 20, code={211,321,2212} \r
   // first step: proton peak, second step: pion peak, third step: kaon peak\r
   // final step: refit all using the parameters\r
   TF1 *fstep1, *fstep2, *fstep3, *fstepTot;\r
   TString modfit = "M0R+";\r
-  Float_t pt=0., ampl=0., mean=0., expKaonMean=0., expProtonMean=0., expPionSig=0., expKaonSig=0., expProtonSig=0.;\r
+  Float_t pt=0., ampl=0.;\r
   Int_t code=TMath::Abs(signedcode);\r
-  GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\r
+  GetInitialParam(h,code,bin,pt,ampl);\r
   if(!IsGoodBin(bin,code)) return;\r
 \r
   printf("___________________________________________________________________ First Step: proton\n");\r
   fstep1 = new TF1("step1",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);\r
   fstep1->SetParameter(0,ampl/16.);       //initial ampl proton`\r
-  fstep1->SetParameter(1,expProtonMean);  //initial mean proton\r
-  fstep1->SetParameter(2,expProtonSig);   //initial sigma proton\r
+  fstep1->SetParameter(1,fExpProtonMean);  //initial mean proton\r
+  fstep1->SetParameter(2,fExpProtonSigma);   //initial sigma proton\r
   fstep1->SetParLimits(0,0,ampl);                                                                    //limits ampl proton\r
-  fstep1->SetParLimits(1,mean,fRangeStep4[1]);                                                       //limits mean proton (dummy)\r
-  fstep1->SetParLimits(2,expProtonSig*fLimitsOnSigmaProton[0],expProtonSig*fLimitsOnSigmaProton[1]); //limits sigma proton\r
+  fstep1->SetParLimits(1,fExpPionMean,fRangeStep4[1]);                                                       //limits mean proton (dummy)\r
+  fstep1->SetParLimits(2,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton\r
 \r
-  if(expProtonSig>0)  h->Fit(fstep1,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);//first fit\r
+  if(fExpProtonSigma>0)  h->Fit(fstep1,modfit,"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//first fit\r
   else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);\r
 \r
   printf("___________________________________________________________________ Second Step: pion\n");\r
@@ -655,13 +685,13 @@ void AliITSsadEdxFitter::DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode,
   fstep2->FixParameter(1,fstep1->GetParameter(1)); //fixed mean proton\r
   fstep2->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma proton\r
   fstep2->SetParameter(3,ampl);             //initial ampl pion\r
-  fstep2->SetParameter(4,mean);             //initial mean pion\r
-  fstep2->SetParameter(5,expPionSig);       //initial sigma pion\r
+  fstep2->SetParameter(4,fExpPionMean);             //initial mean pion\r
+  fstep2->SetParameter(5,fExpPionSigma);       //initial sigma pion\r
   fstep2->SetParLimits(3,0.,ampl);                                                               //limits ampl pion\r
   fstep2->SetParLimits(4,fRangeStep4[0],fstep1->GetParameter(1));                                //limits mean pion\r
-  fstep2->SetParLimits(5,expPionSig*fLimitsOnSigmaPion[0],expPionSig*fLimitsOnSigmaPion[1]);     //limits sigma pion\r
+  fstep2->SetParLimits(5,fExpPionSigma*fLimitsOnSigmaPion[0],fExpPionSigma*fLimitsOnSigmaPion[1]);     //limits sigma pion\r
 \r
-  if(expPionSig>0) h->Fit(fstep2,modfit,"",mean+fRangeStep1[0],mean+fRangeStep1[1]);//second fit\r
+  if(fExpPionSigma>0) h->Fit(fstep2,modfit,"",fExpPionMean+fRangeStep1[0],fExpPionMean+fRangeStep1[1]);//second fit\r
   else for(Int_t npar=0;npar<3;npar++) fstep2->FixParameter(npar,0.00001);\r
 \r
   printf("___________________________________________________________________ Third Step: kaon\n");\r
@@ -673,21 +703,12 @@ void AliITSsadEdxFitter::DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode,
   fstep3->FixParameter(4,fstep2->GetParameter(4)); //fixed mean pion\r
   fstep3->FixParameter(5,fstep2->GetParameter(5)); //fixed sigma pion\r
   fstep3->SetParameter(6,fstep2->GetParameter(0)/8.); //initial ampl kaon\r
-  fstep3->SetParameter(7,expKaonMean);                //initial mean kaon\r
-  fstep3->SetParameter(8,expKaonSig);                 //initial sigma kaon\r
+  fstep3->SetParameter(7,fExpKaonMean);                //initial mean kaon\r
+  fstep3->SetParameter(8,fExpKaonSigma);                 //initial sigma kaon\r
   fstep3->SetParLimits(6,fstep1->GetParameter(0),fstep2->GetParameter(3));                   //limits ampl kaon\r
   fstep3->SetParLimits(7,fstep2->GetParameter(4),fstep1->GetParameter(1));                   //limits mean kaon\r
-  fstep3->SetParLimits(8,expKaonSig*fLimitsOnSigmaKaon[0],expKaonSig*fLimitsOnSigmaKaon[1]); //limits sigma kaon\r
-  /*TLine *l[3];\r
-    l[0] = new TLine(expProtonMean,0,expProtonMean,10000);\r
-    l[1] = new TLine(expProtonMean+fRangeStep3[0],0,expProtonMean+fRangeStep3[0],10000);\r
-    l[2] = new TLine(expProtonMean+fRangeStep3[1],0,expProtonMean+fRangeStep3[1],10000);\r
-    for(Int_t dp=0;dp<3;dp++) {\r
-    l[dp]->Draw("same");\r
-    l[dp]->SetLineColor(2);\r
-    l[dp]->SetLineWidth(4);\r
-    }*/\r
-  if(expKaonSig>0) h->Fit(fstep3,modfit,"",expKaonMean+fRangeStep2[0],expKaonMean+fRangeStep2[1]);//third fit\r
+  fstep3->SetParLimits(8,fExpKaonSigma*fLimitsOnSigmaKaon[0],fExpKaonSigma*fLimitsOnSigmaKaon[1]); //limits sigma kaon\r
+  if(fExpKaonSigma>0) h->Fit(fstep3,modfit,"",fExpKaonMean+fRangeStep2[0],fExpKaonMean+fRangeStep2[1]);//third fit\r
   else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);\r
 \r
   printf("___________________________________________________________________ Final Step: refit all\n");\r
@@ -736,25 +757,25 @@ void AliITSsadEdxFitter::DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode,
 \r
 \r
 //________________________________________________________\r
-void AliITSsadEdxFitter::DoFitOnePeak(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc){\r
+void AliITSsadEdxFitter::DoFitOnePeak(TH1F *h, Int_t bin, Int_t signedcode){\r
   // single-gaussian fit to log(dedx)-log(dedxBB) histogram\r
   TF1 *fstep1;\r
   TString modfit = "M0R+";\r
-  Float_t pt=0., ampl=0., mean=0., expKaonMean=0., expProtonMean=0., expPionSig=0., expKaonSig=0., expProtonSig=0.;\r
+  Float_t pt=0., ampl=0.;\r
   Int_t code=TMath::Abs(signedcode);\r
-  GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\r
+  GetInitialParam(h,code,bin,pt,ampl);\r
   if(!IsGoodBin(bin,code)) return;\r
 \r
   printf("___________________________________________________________________ Single Step\n");\r
   fstep1 = new TF1("step2",SingleGausStep,fRangeStep4[0],fRangeStep4[1],3);\r
   fstep1->SetParameter(0,ampl/16.);                   //initial ampl \r
-  fstep1->SetParameter(1,expProtonMean);              //initial mean \r
-  fstep1->SetParameter(2,expProtonSig);               //initial sigma \r
+  fstep1->SetParameter(1,fExpProtonMean);              //initial mean \r
+  fstep1->SetParameter(2,fExpProtonSigma);               //initial sigma \r
   fstep1->SetParLimits(0,0.,ampl);                                                                   //limits ampl proton\r
-  fstep1->SetParLimits(1,mean,fRangeStep4[1]);                                                       //limits mean proton\r
-  //fstep1->SetParLimits(2,expProtonSig*fLimitsOnSigmaProton[0],expProtonSig*fLimitsOnSigmaProton[1]); //limits sigma proton\r
+  fstep1->SetParLimits(1,fExpPionMean,fRangeStep4[1]);                                                       //limits mean proton\r
+  //fstep1->SetParLimits(2,fExpProtonSigma*fLimitsOnSigmaProton[0],fExpProtonSigma*fLimitsOnSigmaProton[1]); //limits sigma proton\r
 \r
-  if(expProtonSig>0) h->Fit(fstep1,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);//fit\r
+  if(fExpProtonSigma>0) h->Fit(fstep1,modfit,"",fExpProtonMean+fRangeStep3[0],fExpProtonMean+fRangeStep3[1]);//fit\r
   else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);\r
 \r
   fstep1->SetLineColor(1);\r
@@ -772,19 +793,18 @@ void AliITSsadEdxFitter::DoFitTail(TH1F *h, Int_t bin, Int_t signedcode){
   // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed\r
   // final step: refit all using the parameters and tollerance limits (+-20%)\r
   // WARNING: exponential tail added in the right of the Gaussian shape\r
-  Bool_t mc=kFALSE;\r
   Int_t code=TMath::Abs(signedcode);\r
   if(!IsGoodBin(bin,code)) return;\r
 \r
   TF1 *fstep1, *fstep2, *fstep3, *fstepTot;\r
   TString modfit = "M0R+";\r
-  Float_t pt=0., ampl=0., mean=0., expKaonMean=0., expProtonMean=0., expPionSig=0., expKaonSig=0., expProtonSig=0.;\r
-  GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\r
+  Float_t pt=0., ampl=0.;\r
+  GetInitialParam(h,code,bin,pt,ampl);\r
 \r
   printf("\n___________________________________________________________________\n First Step: pions\n\n");\r
   fstep1 = new TF1("step1",SingleGausTail,-3.5,3.5,5);\r
   fstep1->SetParameter(0,ampl);//initial \r
-  fstep1->SetParameter(1,mean);\r
+  fstep1->SetParameter(1,fExpPionMean);\r
   fstep1->SetParameter(3,1.2);\r
   fstep1->SetParameter(4,10.);\r
 \r
@@ -794,7 +814,7 @@ void AliITSsadEdxFitter::DoFitTail(TH1F *h, Int_t bin, Int_t signedcode){
   fstep1->SetParLimits(4,5.,20.);\r
   if(bin<8) fstep1->SetParLimits(4,13.,25.);\r
 \r
-  h->Fit(fstep1,modfit,"",mean-0.45,mean+0.45);//first fit\r
+  h->Fit(fstep1,modfit,"",fExpPionMean-0.45,fExpPionMean+0.45);//first fit\r
 \r
   printf("\n___________________________________________________________________\n Second Step: kaons\n\n"); \r
   fstep2 = new TF1("fstep2",DoubleGausTail,-3.5,3.5,10);\r
@@ -950,7 +970,8 @@ void AliITSsadEdxFitter::GetFitPar(Double_t *fitpar, Double_t *fitparerr) const
 \r
 //________________________________________________________\r
 void AliITSsadEdxFitter::PrintAll() const{\r
-  //\r
+  // print parameters\r
+\r
   printf("Range 1 = %f %f\n",fRangeStep1[0],fRangeStep1[1]);\r
   printf("Range 2 = %f %f\n",fRangeStep2[0],fRangeStep2[1]);\r
   printf("Range 3 = %f %f\n",fRangeStep3[0],fRangeStep3[1]);\r

diff --git a/PWG2/SPECTRA/AliITSsadEdxFitter.h b/PWG2/SPECTRA/AliITSsadEdxFitter.h
index a353d65024e63a37bf2a175fb098f9518e1afa26..65f2b0eedd8a8963a26729f960b5ffe4597e2912 100644 (file)
--- a/PWG2/SPECTRA/AliITSsadEdxFitter.h
+++ b/PWG2/SPECTRA/AliITSsadEdxFitter.h
@@ -13,28 +13,44 @@
 ///////////////////////////////////////////////////////////////////////\r
 \r
 #include <TObject.h>\r
+#include "AliITSPIDResponse.h"\r
+\r
 class TGraph;\r
 \r
 class AliITSsadEdxFitter  : public TObject {\r
 \r
  public:\r
-  AliITSsadEdxFitter();  \r
-  virtual ~AliITSsadEdxFitter(){};\r
+  AliITSsadEdxFitter();\r
+  AliITSsadEdxFitter(Bool_t isMC);\r
+  virtual ~AliITSsadEdxFitter(){\r
+    delete fITSpid;\r
+  };\r
 \r
-  static Double_t CalcSigma(Int_t code,Float_t x, Bool_t mc);\r
-  static Int_t CalcMean(Int_t code,Float_t x, Float_t mean0, Float_t &mean1, Float_t &mean2);\r
+  Double_t CalcSigma(Int_t code,Float_t x) const;\r
+  Int_t CalcMean(Int_t code,Float_t x, Float_t mean0, Float_t &mean1, Float_t &mean2) const;\r
 \r
   void GetFitPar(Double_t *fitpar, Double_t *fitparerr) const;\r
   void DoFitTail(TH1F *h, Int_t bin, Int_t code);\r
-  void DoFit(TH1F *h, Int_t bin, Int_t code, Bool_t mc, TGraph *gres);\r
-  void DoFitProton(TH1F *h, Int_t bin, Int_t code, Bool_t mc, TGraph *gres);\r
-       void DoFitOnePeak(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc);\r
-       void DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, TGraph *gres);\r
-       void GetInitialParam(TH1F* h,Bool_t mc,Int_t code,Int_t bin, Float_t &pt, Float_t &ampl, Float_t &mean1, Float_t &mean2, Float_t &mean3, Float_t &sigma1, Float_t &sigma2, Float_t &sigma3);\r
+  void DoFit(TH1F *h, Int_t bin, Int_t code, TGraph *gres);\r
+  void DoFitProton(TH1F *h, Int_t bin, Int_t code, TGraph *gres);\r
+  void DoFitOnePeak(TH1F *h, Int_t bin, Int_t signedcode);\r
+  void DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode, TGraph *gres);\r
+  void GetInitialParam(TH1F* h,Int_t code,Int_t bin, Float_t &pt, Float_t &ampl); \r
   void FillHisto(TH1F *hsps, Int_t bin, Float_t binsize, Int_t code);\r
   void FillHistoMC(TH1F *hsps, Int_t bin, Int_t code, TH1F *h);\r
   Bool_t IsGoodBin(Int_t bin,Int_t code);\r
 \r
+  void SetMCConfig(){\r
+    fIsMC=kTRUE;\r
+    if(fITSpid) delete fITSpid;\r
+    fITSpid=new AliITSPIDResponse(kTRUE);        \r
+  }\r
+  void SetDataConfig(){\r
+    fIsMC=kFALSE;\r
+    if(fITSpid) delete fITSpid;\r
+    fITSpid=new AliITSPIDResponse(kFALSE);       \r
+  }\r
+\r
   void SetRangeStep1(Double_t dxlow=-0.2, Double_t dxup=0.3){\r
     fRangeStep1[0]=dxlow;\r
     fRangeStep1[1]=dxup;\r
@@ -63,18 +79,18 @@ class AliITSsadEdxFitter  : public TObject {
     fLimitsOnSigmaProton[0]=smin;\r
     fLimitsOnSigmaProton[1]=smax;\r
   }\r
-       void SetBinsUsedPion(Int_t bmin=2, Int_t bmax=14){\r
-               fBinsUsedPion[0]=bmin;\r
-               fBinsUsedPion[1]=bmax;\r
-       }\r
-       void SetBinsUsedKaon(Int_t bmin=7, Int_t bmax=12){\r
-               fBinsUsedKaon[0]=bmin;\r
-               fBinsUsedKaon[1]=bmax;\r
-       }\r
-       void SetBinsUsedProton(Int_t bmin=8, Int_t bmax=17){\r
-               fBinsUsedProton[0]=bmin;\r
-               fBinsUsedProton[1]=bmax;\r
-       }\r
+  void SetBinsUsedPion(Int_t bmin=2, Int_t bmax=14){\r
+    fBinsUsedPion[0]=bmin;\r
+    fBinsUsedPion[1]=bmax;\r
+  }\r
+  void SetBinsUsedKaon(Int_t bmin=7, Int_t bmax=12){\r
+    fBinsUsedKaon[0]=bmin;\r
+    fBinsUsedKaon[1]=bmax;\r
+  }\r
+  void SetBinsUsedProton(Int_t bmin=8, Int_t bmax=17){\r
+    fBinsUsedProton[0]=bmin;\r
+    fBinsUsedProton[1]=bmax;\r
+  }\r
 \r
   void PrintAll() const;\r
   void CalcResidual(TH1F *h,TF1 *fun,TGraph *gres) const;\r
@@ -83,6 +99,9 @@ class AliITSsadEdxFitter  : public TObject {
   void DrawFitFunction(TF1 *fun) const;\r
 \r
  private:\r
+  AliITSsadEdxFitter(const AliITSsadEdxFitter& s);\r
+  AliITSsadEdxFitter& operator=(const AliITSsadEdxFitter& s);\r
+\r
   Double_t fFitpar[3];     // array with fit parameters\r
   Double_t fFitparErr[3];  // array with fit parameter errors \r
   Double_t fRangeStep1[2]; // Range for Step1 (w.r.t pion peak)\r
@@ -92,11 +111,22 @@ class AliITSsadEdxFitter  : public TObject {
   Double_t fLimitsOnSigmaPion[2]; // limits on sigma pions\r
   Double_t fLimitsOnSigmaKaon[2]; // limits on sigma kaons\r
   Double_t fLimitsOnSigmaProton[2]; // limits on sigma protons\r
+  \r
+  Double_t fExpPionMean;    // expected mean for pion peak\r
+  Double_t fExpKaonMean;    // expected mean for kaon peak\r
+  Double_t fExpProtonMean;  // expected mean for proton peak\r
+  Double_t fExpPionSigma;   // expected sigma for pion peak\r
+  Double_t fExpKaonSigma;   // expected sigma for kaon peak\r
+  Double_t fExpProtonSigma; // expected sigma for proton peak\r
+\r
   Int_t fBinsUsedPion[2];   // limits on bins used pions\r
   Int_t fBinsUsedKaon[2];   // limits on bins used kaons\r
   Int_t fBinsUsedProton[2]; // limits on bins used protons\r
 \r
-  ClassDef(AliITSsadEdxFitter,2);\r
+  Bool_t fIsMC;                 // flag MC/data\r
+  AliITSPIDResponse* fITSpid;   // ITS pid object\r
+\r
+  ClassDef(AliITSsadEdxFitter,3);\r
 };\r
 \r
 #endif\r