\r
\r
ClassImp(AliITSsadEdxFitter)\r
- //______________________________________________________________________\r
- AliITSsadEdxFitter::AliITSsadEdxFitter():TObject(){\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
- };\r
+//______________________________________________________________________\r
+AliITSsadEdxFitter::AliITSsadEdxFitter():TObject(){\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
+};\r
\r
//________________________________________________________\r
Double_t AliITSsadEdxFitter::CalcSigma(Int_t code,Float_t x,Bool_t mc){\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(code==211){\r
- p[0] = -1.20337e-04;\r
- p[1] = 1.13060e-01;\r
- } \r
- else if(code==321 && x>xMinKaon){\r
- p[0] = -2.39631e-03;\r
- p[1] = 1.15723e-01;\r
- } \r
- else if(code==2212 && x>xMinProton){\r
- p[0] = -8.34576e-03;\r
- p[1] = 1.34237e-01;\r
- } \r
- else return -1;\r
- }\r
- else{\r
- if(code==211){\r
- p[0] = -6.55200e-05;\r
- p[1] = 1.26657e-01;\r
- } \r
- else if(code==321 && x>xMinKaon){\r
- p[0] = -6.22639e-04;\r
- p[1] = 1.43289e-01;\r
- } \r
- else if(code==2212 && x>xMinProton){\r
- p[0] = -2.13243e-03;\r
- p[1] = 1.68614e-01;\r
- } \r
- else return -1;\r
- }\r
- return p[0]/(x*x)*TMath::Log(x)+p[1];\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(code==211){\r
+ p[0] = -1.20337e-04;\r
+ p[1] = 1.13060e-01;\r
+ } \r
+ else if(code==321 && x>xMinKaon){\r
+ p[0] = -2.39631e-03;\r
+ p[1] = 1.15723e-01;\r
+ } \r
+ else if(code==2212 && x>xMinProton){\r
+ p[0] = -8.34576e-03;\r
+ p[1] = 1.34237e-01;\r
+ } \r
+ else return -1;\r
+ }\r
+ else{\r
+ if(code==211){\r
+ p[0] = -6.55200e-05;\r
+ p[1] = 1.26657e-01;\r
+ } \r
+ else if(code==321 && x>xMinKaon){\r
+ p[0] = -6.22639e-04;\r
+ p[1] = 1.43289e-01;\r
+ } \r
+ else if(code==2212 && x>xMinProton){\r
+ p[0] = -2.13243e-03;\r
+ p[1] = 1.68614e-01;\r
+ } \r
+ else return -1;\r
+ }\r
+ return p[0]/(x*x)*TMath::Log(x)+p[1];\r
}\r
\r
//_______________________________________________________\r
Int_t AliITSsadEdxFitter::CalcMean(Int_t code, Float_t x, Float_t mean0, Float_t &mean1, Float_t &mean2){\r
- // calculate the mean 12-ott-2010 \r
- Double_t p1[4]={0.};\r
- Double_t p2[4]={0.};\r
- if(code==211){\r
- p1[0] = 1.77049;\r
- p1[1] = -2.65469;\r
- p2[0] = 0.890856;\r
- p2[1] = -0.276719;\r
- mean1 = mean0 + p1[0]+ p1[1]*x + p1[2]*x*x + p1[3]*x*x*x;\r
- mean2 = mean1 + p2[0]+ p2[1]*x + p2[2]*x*x + p2[3]*x*x*x;\r
- }\r
- else if(code==321){\r
- p2[0] = 1.57664;\r
- p2[1] = -6.88635;\r
- p2[2] = 18.702;\r
- p2[3] = -16.3385;\r
- mean1 = 0.;\r
- mean2 = mean1 + p2[0]+ p2[1]*x + p2[2]*x*x + p2[3]*x*x*x;\r
- }\r
- else if(code==2212){\r
- p1[0] = 4.24861; \r
- p1[1] = -19.178;\r
- p1[2] = 31.5947;\r
- p1[3] = -18.178;\r
- mean1 = mean0 + p1[0]+ p1[1]*x + p1[2]*x*x + p1[3]*x*x*x;\r
- mean2 = 0.; \r
- }\r
- else return -1;\r
- return 0;\r
+ // calculate the mean 12-ott-2010 \r
+ Double_t p1[4]={0.};\r
+ Double_t p2[4]={0.};\r
+ if(code==211){\r
+ p1[0] = 1.77049;\r
+ p1[1] = -2.65469;\r
+ p2[0] = 0.890856;\r
+ p2[1] = -0.276719;\r
+ mean1 = mean0 + p1[0]+ p1[1]*x + p1[2]*x*x + p1[3]*x*x*x;\r
+ mean2 = mean1 + p2[0]+ p2[1]*x + p2[2]*x*x + p2[3]*x*x*x;\r
+ }\r
+ else if(code==321){\r
+ p2[0] = 1.57664;\r
+ p2[1] = -6.88635;\r
+ p2[2] = 18.702;\r
+ p2[3] = -16.3385;\r
+ mean1 = 0.;\r
+ mean2 = mean1 + p2[0]+ p2[1]*x + p2[2]*x*x + p2[3]*x*x*x;\r
+ }\r
+ else if(code==2212){\r
+ p1[0] = 4.24861; \r
+ p1[1] = -19.178;\r
+ p1[2] = 31.5947;\r
+ p1[3] = -18.178;\r
+ mean1 = mean0 + p1[0]+ p1[1]*x + p1[2]*x*x + p1[3]*x*x*x;\r
+ mean2 = 0.; \r
+ }\r
+ else return -1;\r
+ return 0;\r
}\r
\r
//________________________________________________________\r
Bool_t AliITSsadEdxFitter::IsGoodBin(Int_t bin,Int_t code){\r
- //method to select the bins used for the analysis only\r
- Bool_t retvalue=kTRUE;\r
- TLine *l[2]; //for the cross \r
- l[0]=new TLine(-2.1,0,2.,100.);\r
- l[1]=new TLine(-1.9,120,2.,0.);\r
- for(Int_t j=0;j<2;j++){\r
- l[j]->SetLineColor(4);\r
- l[j]->SetLineWidth(5);\r
- }\r
-\r
- if(code==211 && (bin<fBinsUsedPion[0] || bin>fBinsUsedPion[1])){\r
- for(Int_t j=0;j<2;j++) l[j]->Draw("same"); \r
- retvalue=kFALSE;\r
- }\r
- if(code==321 && (bin<fBinsUsedKaon[0] || bin>fBinsUsedKaon[1])){\r
- for(Int_t j=0;j<2;j++) l[j]->Draw("same"); \r
- retvalue=kFALSE;\r
- }\r
- if(code==2212 && (bin<fBinsUsedProton[0] || bin>fBinsUsedProton[1])){\r
- for(Int_t j=0;j<2;j++) l[j]->Draw("same"); \r
- retvalue=kFALSE;\r
- }\r
- return retvalue;\r
+ //method to select the bins used for the analysis only\r
+ Bool_t retvalue=kTRUE;\r
+ TLine *l[2]; //for the cross \r
+ l[0]=new TLine(-2.1,0,2.,100.);\r
+ l[1]=new TLine(-1.9,120,2.,0.);\r
+ for(Int_t j=0;j<2;j++){\r
+ l[j]->SetLineColor(4);\r
+ l[j]->SetLineWidth(5);\r
+ }\r
+\r
+ if(code==211 && (bin<fBinsUsedPion[0] || bin>fBinsUsedPion[1])){\r
+ for(Int_t j=0;j<2;j++) l[j]->Draw("same"); \r
+ retvalue=kFALSE;\r
+ }\r
+ if(code==321 && (bin<fBinsUsedKaon[0] || bin>fBinsUsedKaon[1])){\r
+ for(Int_t j=0;j<2;j++) l[j]->Draw("same"); \r
+ retvalue=kFALSE;\r
+ }\r
+ if(code==2212 && (bin<fBinsUsedProton[0] || bin>fBinsUsedProton[1])){\r
+ for(Int_t j=0;j<2;j++) l[j]->Draw("same"); \r
+ retvalue=kFALSE;\r
+ }\r
+ return retvalue;\r
}\r
\r
//________________________________________________________\r
Double_t SingleGausTail(const Double_t *x, const Double_t *par){\r
- //single gaussian with exponential tail\r
- Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
- Double_t xx = x[0];\r
- Double_t mean1 = par[1];\r
- Double_t sigma1 = par[2];\r
- Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);\r
- Double_t sigmaSquare1 = sigma1 * sigma1;\r
- Double_t xdiv1 = mean1 + par[3] * sigma1;\r
- Double_t y1=0.0;\r
- if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);\r
- else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));\r
- return y1;\r
+ //single gaussian with exponential tail\r
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
+ Double_t xx = x[0];\r
+ Double_t mean1 = par[1];\r
+ Double_t sigma1 = par[2];\r
+ Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);\r
+ Double_t sigmaSquare1 = sigma1 * sigma1;\r
+ Double_t xdiv1 = mean1 + par[3] * sigma1;\r
+ Double_t y1=0.0;\r
+ if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);\r
+ else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));\r
+ return y1;\r
}\r
\r
//________________________________________________________\r
Double_t DoubleGausTail(const Double_t *x, const Double_t *par){\r
- //sum of two gaussians with exponential tail\r
- Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
- Double_t xx = x[0];\r
- Double_t mean1 = par[1];\r
- Double_t sigma1 = par[2];\r
- Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);\r
- Double_t sigmaSquare1 = sigma1 * sigma1;\r
- Double_t xdiv1 = mean1 + par[3] * sigma1;\r
- Double_t y1=0.0;\r
- if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);\r
- else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));\r
-\r
- Double_t mean2 = par[6];\r
- Double_t sigma2 = par[7];\r
- Double_t xNormSquare2 = (xx - mean2) * (xx - mean2);\r
- Double_t sigmaSquare2 = sigma2 * sigma2;\r
- Double_t xdiv2 = mean2 + par[8] * sigma2;\r
- Double_t y2=0.0;\r
- if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2);\r
- else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8]));\r
- return y1+y2;\r
+ //sum of two gaussians with exponential tail\r
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
+ Double_t xx = x[0];\r
+ Double_t mean1 = par[1];\r
+ Double_t sigma1 = par[2];\r
+ Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);\r
+ Double_t sigmaSquare1 = sigma1 * sigma1;\r
+ Double_t xdiv1 = mean1 + par[3] * sigma1;\r
+ Double_t y1=0.0;\r
+ if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);\r
+ else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));\r
+\r
+ Double_t mean2 = par[6];\r
+ Double_t sigma2 = par[7];\r
+ Double_t xNormSquare2 = (xx - mean2) * (xx - mean2);\r
+ Double_t sigmaSquare2 = sigma2 * sigma2;\r
+ Double_t xdiv2 = mean2 + par[8] * sigma2;\r
+ Double_t y2=0.0;\r
+ if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2);\r
+ else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8]));\r
+ return y1+y2;\r
}\r
\r
//________________________________________________________\r
Double_t FinalGausTail(const Double_t *x, const Double_t *par){\r
- //sum of three gaussians with exponential tail\r
- Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
- Double_t xx = x[0];\r
- Double_t mean1 = par[1];\r
- Double_t sigma1 = par[2];\r
- Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);\r
- Double_t sigmaSquare1 = sigma1 * sigma1;\r
- Double_t xdiv1 = mean1 + par[3] * sigma1;\r
- Double_t y1=0.0;\r
- if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);\r
- else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));\r
-\r
- Double_t mean2 = par[6];\r
- Double_t sigma2 = par[7];\r
- Double_t xNormSquare2 = (xx - mean2) * (xx - mean2);\r
- Double_t sigmaSquare2 = sigma2 * sigma2;\r
- Double_t xdiv2 = mean2 + par[8] * sigma2;\r
- Double_t y2=0.0;\r
- if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2);\r
- else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8]));\r
-\r
- Double_t mean3 = par[11];\r
- Double_t sigma3 = par[12];\r
- Double_t xNormSquare3 = (xx - mean3) * (xx - mean3);\r
- Double_t sigmaSquare3 = sigma3 * sigma3;\r
- Double_t xdiv3 = mean3 + par[13] * sigma3;\r
- Double_t y3=0.0;\r
- if(xx < xdiv3) y3 = par[10]/(s2pi*par[12]) * TMath::Exp(-0.5 * xNormSquare3 / sigmaSquare3);\r
- else y3 = TMath::Exp(-par[14]*(xx-xdiv3)) * par[10] / (s2pi*par[12]) * TMath::Exp(-0.5*(par[13]*par[13]));\r
- return y1+y2+y3;\r
+ //sum of three gaussians with exponential tail\r
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
+ Double_t xx = x[0];\r
+ Double_t mean1 = par[1];\r
+ Double_t sigma1 = par[2];\r
+ Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);\r
+ Double_t sigmaSquare1 = sigma1 * sigma1;\r
+ Double_t xdiv1 = mean1 + par[3] * sigma1;\r
+ Double_t y1=0.0;\r
+ if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);\r
+ else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));\r
+\r
+ Double_t mean2 = par[6];\r
+ Double_t sigma2 = par[7];\r
+ Double_t xNormSquare2 = (xx - mean2) * (xx - mean2);\r
+ Double_t sigmaSquare2 = sigma2 * sigma2;\r
+ Double_t xdiv2 = mean2 + par[8] * sigma2;\r
+ Double_t y2=0.0;\r
+ if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2);\r
+ else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8]));\r
+\r
+ Double_t mean3 = par[11];\r
+ Double_t sigma3 = par[12];\r
+ Double_t xNormSquare3 = (xx - mean3) * (xx - mean3);\r
+ Double_t sigmaSquare3 = sigma3 * sigma3;\r
+ Double_t xdiv3 = mean3 + par[13] * sigma3;\r
+ Double_t y3=0.0;\r
+ if(xx < xdiv3) y3 = par[10]/(s2pi*par[12]) * TMath::Exp(-0.5 * xNormSquare3 / sigmaSquare3);\r
+ else y3 = TMath::Exp(-par[14]*(xx-xdiv3)) * par[10] / (s2pi*par[12]) * TMath::Exp(-0.5*(par[13]*par[13]));\r
+ return y1+y2+y3;\r
}\r
\r
//______________________________________________________________________\r
void AliITSsadEdxFitter::CalcResidual(TH1F *h,TF1 *fun,TGraph *gres) const{\r
- //code to calculate the difference fit function and histogram point (residual)\r
- //to use as goodness test for the fit\r
- Int_t ipt=0;\r
- Double_t x=0.,yhis=0.,yfun=0.;\r
- for(Int_t i=0;i<h->GetNbinsX();i++){\r
- x=(h->GetBinLowEdge(i+2)+h->GetBinLowEdge(i+1))/2;\r
- yfun=fun->Eval(x);\r
- yhis=h->GetBinContent(i+1);\r
- if(yhis>0. && yfun>0.) {\r
- gres->SetPoint(ipt,x,(yhis-yfun)/yhis);\r
- ipt++;\r
- }\r
- }\r
- return;\r
+ //code to calculate the difference fit function and histogram point (residual)\r
+ //to use as goodness test for the fit\r
+ Int_t ipt=0;\r
+ Double_t x=0.,yhis=0.,yfun=0.;\r
+ for(Int_t i=0;i<h->GetNbinsX();i++){\r
+ x=(h->GetBinLowEdge(i+2)+h->GetBinLowEdge(i+1))/2;\r
+ yfun=fun->Eval(x);\r
+ yhis=h->GetBinContent(i+1);\r
+ if(yhis>0. && yfun>0.) {\r
+ gres->SetPoint(ipt,x,(yhis-yfun)/yhis);\r
+ ipt++;\r
+ }\r
+ }\r
+ return;\r
}\r
\r
\r
//________________________________________________________\r
Double_t SingleGausStep(const Double_t *x, const Double_t *par){\r
- //single normalized gaussian\r
- Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
- Double_t xx = x[0];\r
- Double_t mean1 = par[1];\r
- Double_t sigma1 = par[2];\r
- Double_t xNorm1Square = (xx - mean1) * (xx - mean1);\r
- Double_t sigma1Square = sigma1 * sigma1;\r
- Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNorm1Square / sigma1Square);\r
- return step1;\r
+ //single normalized gaussian\r
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
+ Double_t xx = x[0];\r
+ Double_t mean1 = par[1];\r
+ Double_t sigma1 = par[2];\r
+ Double_t xNorm1Square = (xx - mean1) * (xx - mean1);\r
+ Double_t sigma1Square = sigma1 * sigma1;\r
+ Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNorm1Square / sigma1Square);\r
+ return step1;\r
}\r
\r
//________________________________________________________\r
Double_t DoubleGausStep(const Double_t *x, const Double_t *par){\r
- //sum of two normalized gaussians\r
- Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
- Double_t xx = x[0];\r
- Double_t mean1 = par[1];\r
- Double_t sigma1 = par[2];\r
- Double_t xNorm1Square = (xx - mean1) * (xx - mean1);\r
- Double_t sigma1Square = sigma1 * sigma1;\r
- Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square );\r
- Double_t mean2 = par[4];\r
- Double_t sigma2 = par[5];\r
- Double_t xNorm2Square = (xx - mean2) * (xx - mean2);\r
- Double_t sigma2Square = sigma2 * sigma2;\r
- Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square );\r
- return step1+step2;\r
+ //sum of two normalized gaussians\r
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
+ Double_t xx = x[0];\r
+ Double_t mean1 = par[1];\r
+ Double_t sigma1 = par[2];\r
+ Double_t xNorm1Square = (xx - mean1) * (xx - mean1);\r
+ Double_t sigma1Square = sigma1 * sigma1;\r
+ Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square );\r
+ Double_t mean2 = par[4];\r
+ Double_t sigma2 = par[5];\r
+ Double_t xNorm2Square = (xx - mean2) * (xx - mean2);\r
+ Double_t sigma2Square = sigma2 * sigma2;\r
+ Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square );\r
+ return step1+step2;\r
}\r
\r
//________________________________________________________\r
Double_t FinalGausStep(const Double_t *x, const Double_t *par){\r
- //sum of three normalized gaussians\r
- Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
- Double_t xx = x[0];\r
- Double_t mean1 = par[1];\r
- Double_t sigma1 = par[2];\r
- Double_t xNorm1Square = (xx - mean1) * (xx - mean1);\r
- Double_t sigma1Square = sigma1 * sigma1;\r
- Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square );\r
- Double_t mean2 = par[4];\r
- Double_t sigma2 = par[5];\r
- Double_t xNorm2Square = (xx - mean2) * (xx - mean2);\r
- Double_t sigma2Square = sigma2 * sigma2;\r
- Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square );\r
- Double_t mean3 = par[7];\r
- Double_t sigma3 = par[8];\r
- Double_t xNorm3Square = (xx - mean3) * (xx - mean3);\r
- Double_t sigma3Square = sigma3 * sigma3;\r
- Double_t step3 = par[6]/(s2pi*par[8]) * TMath::Exp( - 0.5 * xNorm3Square / sigma3Square);\r
- return step1+step2+step3;\r
+ //sum of three normalized gaussians\r
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());\r
+ Double_t xx = x[0];\r
+ Double_t mean1 = par[1];\r
+ Double_t sigma1 = par[2];\r
+ Double_t xNorm1Square = (xx - mean1) * (xx - mean1);\r
+ Double_t sigma1Square = sigma1 * sigma1;\r
+ Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square );\r
+ Double_t mean2 = par[4];\r
+ Double_t sigma2 = par[5];\r
+ Double_t xNorm2Square = (xx - mean2) * (xx - mean2);\r
+ Double_t sigma2Square = sigma2 * sigma2;\r
+ Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square );\r
+ Double_t mean3 = par[7];\r
+ Double_t sigma3 = par[8];\r
+ Double_t xNorm3Square = (xx - mean3) * (xx - mean3);\r
+ Double_t sigma3Square = sigma3 * sigma3;\r
+ Double_t step3 = par[6]/(s2pi*par[8]) * TMath::Exp( - 0.5 * xNorm3Square / sigma3Square);\r
+ return step1+step2+step3;\r
}\r
\r
//______________________________________________________________________\r
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{\r
- //gaussian with an exponential tail on the right side\r
- Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi()));\r
- Double_t returnvalue=0.0;\r
- 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);\r
- if (x<mean+cut*rms) returnvalue=TMath::Exp(-1.0*(x-mean)*(x-mean)/(2.0*rms*rms))*factor/TMath::Abs(rms);\r
- else returnvalue=TMath::Exp(slope*(mean+cut*rms-x))*TMath::Exp(-cut*cut*0.5)*factor/TMath::Abs(rms);\r
- return c*returnvalue/n;\r
+ //gaussian with an exponential tail on the right side\r
+ Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi()));\r
+ Double_t returnvalue=0.0;\r
+ 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);\r
+ if (x<mean+cut*rms) returnvalue=TMath::Exp(-1.0*(x-mean)*(x-mean)/(2.0*rms*rms))*factor/TMath::Abs(rms);\r
+ else returnvalue=TMath::Exp(slope*(mean+cut*rms-x))*TMath::Exp(-cut*cut*0.5)*factor/TMath::Abs(rms);\r
+ return c*returnvalue/n;\r
}\r
\r
\r
//______________________________________________________________________\r
Double_t AliITSsadEdxFitter::GausOnBackground(const Double_t* x, const Double_t *par) const {\r
- //gaussian with a background parametrisation \r
-
//cout<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<<" "<<x[0]<< endl;
\r- Double_t returnvalue=0.0;\r
- Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi()));\r
- 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]);\r
- else returnvalue+=GausPlusTail(x[0], par[0], par[1],par[2], par[6], 1.2);\r
- returnvalue+=par[3]*TMath::Exp((par[5]-x[0])*par[4]);\r
- return returnvalue;\r
+ //gaussian with a background parametrisation \r
+ //cout<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<<" "<<x[0]<< endl;\r
+ Double_t returnvalue=0.0;\r
+ Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi()));\r
+ 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]);\r
+ else returnvalue+=GausPlusTail(x[0], par[0], par[1],par[2], par[6], 1.2);\r
+ returnvalue+=par[3]*TMath::Exp((par[5]-x[0])*par[4]);\r
+ return returnvalue;\r
}\r
\r
//______________________________________________________________________\r
void AliITSsadEdxFitter::DrawFitFunction(TF1 *fun) const {\r
- //code to draw the final fit function and the single gaussians used\r
- //to extract the yields for the 3 species\r
- TF1 *fdraw1=new TF1("fdraw1",SingleGausStep,-3.5,3.5,3);\r
- TF1 *fdraw2=new TF1("fdraw2",SingleGausStep,-3.5,3.5,3);\r
- TF1 *fdraw3=new TF1("fdraw3",SingleGausStep,-3.5,3.5,3);\r
- fdraw1->SetParameter(0,fun->GetParameter(0));\r
- fdraw1->SetParameter(1,fun->GetParameter(1));\r
- fdraw1->SetParameter(2,fun->GetParameter(2));\r
- fdraw2->SetParameter(0,fun->GetParameter(3));\r
- fdraw2->SetParameter(1,fun->GetParameter(4));\r
- fdraw2->SetParameter(2,fun->GetParameter(5));\r
- fdraw3->SetParameter(0,fun->GetParameter(6));\r
- fdraw3->SetParameter(1,fun->GetParameter(7));\r
- fdraw3->SetParameter(2,fun->GetParameter(8));\r
-\r
- fdraw1->SetLineColor(6);//color code\r
- fdraw2->SetLineColor(2);\r
- fdraw3->SetLineColor(4); \r
-\r
- fdraw1->SetLineStyle(2);//dot lines\r
- fdraw2->SetLineStyle(2);\r
- fdraw3->SetLineStyle(2);\r
-\r
- fun->SetLineWidth(3);\r
- fdraw1->DrawCopy("same");\r
- fdraw2->DrawCopy("same");\r
- fdraw3->DrawCopy("same");\r
- fun->DrawCopy("same");\r
-\r
- TLatex *ltx[3];\r
- for(Int_t j=0;j<3;j++){\r
- ltx[0]=new TLatex(0.13,0.25,"pions");\r
- ltx[1]=new TLatex(0.13,0.20,"kaons");\r
- ltx[2]=new TLatex(0.13,0.15,"protons");\r
- ltx[0]->SetTextColor(6);\r
- ltx[1]->SetTextColor(2);\r
- ltx[2]->SetTextColor(4);\r
- ltx[j]->SetNDC();\r
- ltx[j]->Draw();\r
- }\r
- return;\r
+ //code to draw the final fit function and the single gaussians used\r
+ //to extract the yields for the 3 species\r
+ TF1 *fdraw1=new TF1("fdraw1",SingleGausStep,-3.5,3.5,3);\r
+ TF1 *fdraw2=new TF1("fdraw2",SingleGausStep,-3.5,3.5,3);\r
+ TF1 *fdraw3=new TF1("fdraw3",SingleGausStep,-3.5,3.5,3);\r
+ fdraw1->SetParameter(0,fun->GetParameter(0));\r
+ fdraw1->SetParameter(1,fun->GetParameter(1));\r
+ fdraw1->SetParameter(2,fun->GetParameter(2));\r
+ fdraw2->SetParameter(0,fun->GetParameter(3));\r
+ fdraw2->SetParameter(1,fun->GetParameter(4));\r
+ fdraw2->SetParameter(2,fun->GetParameter(5));\r
+ fdraw3->SetParameter(0,fun->GetParameter(6));\r
+ fdraw3->SetParameter(1,fun->GetParameter(7));\r
+ fdraw3->SetParameter(2,fun->GetParameter(8));\r
+\r
+ fdraw1->SetLineColor(6);//color code\r
+ fdraw2->SetLineColor(2);\r
+ fdraw3->SetLineColor(4); \r
+\r
+ fdraw1->SetLineStyle(2);//dot lines\r
+ fdraw2->SetLineStyle(2);\r
+ fdraw3->SetLineStyle(2);\r
+\r
+ fun->SetLineWidth(3);\r
+ fdraw1->DrawCopy("same");\r
+ fdraw2->DrawCopy("same");\r
+ fdraw3->DrawCopy("same");\r
+ fun->DrawCopy("same");\r
+\r
+ TLatex *ltx[3];\r
+ for(Int_t j=0;j<3;j++){\r
+ ltx[0]=new TLatex(0.13,0.25,"pions");\r
+ ltx[1]=new TLatex(0.13,0.20,"kaons");\r
+ ltx[2]=new TLatex(0.13,0.15,"protons");\r
+ ltx[0]->SetTextColor(6);\r
+ ltx[1]->SetTextColor(2);\r
+ ltx[2]->SetTextColor(4);\r
+ ltx[j]->SetNDC();\r
+ ltx[j]->Draw();\r
+ }\r
+ return;\r
}\r
\r
//______________________________________________________________________\r
void AliITSsadEdxFitter::GetInitialParam(TH1F* h,Bool_t mc,Int_t code,Int_t bin, Float_t &pt, Float_t &l, Float_t &mean1, Float_t &mean2, Float_t &mean3, Float_t &sigma1, Float_t &sigma2, Float_t &sigma3){\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
-\r
- //draw and label\r
- h->SetTitle(Form("p_{t}=[%1.2f,%1.2f], code=%d",xbins[bin],xbins[bin+1],code));\r
- h->GetXaxis()->SetTitle("[ln dE/dx]_{meas} - [ln dE/dx(i)]_{calc}");\r
- h->GetYaxis()->SetTitle("counts");\r
- h->Draw("e");\r
- h->SetFillColor(11);\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
+\r
+ //draw and label\r
+ h->SetTitle(Form("p_{t}=[%1.2f,%1.2f], code=%d",xbins[bin],xbins[bin+1],code));\r
+ h->GetXaxis()->SetTitle("[ln dE/dx]_{meas} - [ln dE/dx(i)]_{calc}");\r
+ h->GetYaxis()->SetTitle("counts");\r
+ h->Draw("e");\r
+ 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
- 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,sigma1,sigma1);\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
+ //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
+ 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,sigma1,sigma1);\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
- // 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
- 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
- Int_t code=TMath::Abs(signedcode);\r
- GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\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->SetParLimits(1,fRangeStep4[0],fRangeStep4[1]); //limits mean pion (dummy)\r
- fstep1->SetParLimits(2,expPionSig*fLimitsOnSigmaPion[0],expPionSig*fLimitsOnSigmaPion[1]); //limits sigma pion\r
-\r
- if(expPionSig>0) h->Fit(fstep1,modfit.Data(),"",mean+fRangeStep1[0],mean+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
- fstep2 = new TF1("fstep2",DoubleGausStep,fRangeStep4[0],fRangeStep4[1],6);\r
- fstep2->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion\r
- 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->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
-\r
- if(expKaonSig>0) h->Fit(fstep2,modfit.Data(),"",expKaonMean+fRangeStep2[0],expKaonMean+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
- fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion\r
- fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion\r
- fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion\r
- fstep3->FixParameter(3,fstep2->GetParameter(3)); //fixed ampl kaon\r
- 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->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
-\r
- if(expProtonSig>0) h->Fit(fstep3,modfit.Data(),"",expProtonMean+fRangeStep3[0],expProtonMean+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
- fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
- fstepTot->SetLineColor(1);\r
-\r
- Double_t initialParametersStepTot[9];\r
- initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
- initialParametersStepTot[1] = fstep1->GetParameter(1);\r
- initialParametersStepTot[2] = fstep1->GetParameter(2);\r
-\r
- initialParametersStepTot[3] = fstep2->GetParameter(3);//second gaussian\r
- initialParametersStepTot[4] = fstep2->GetParameter(4);\r
- initialParametersStepTot[5] = fstep2->GetParameter(5);\r
-\r
- initialParametersStepTot[6] = fstep3->GetParameter(6);//third gaussian\r
- initialParametersStepTot[7] = fstep3->GetParameter(7);\r
- initialParametersStepTot[8] = fstep3->GetParameter(8); \r
-\r
- fstepTot->SetParameters(initialParametersStepTot); //initial parameters\r
- fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl pion\r
- fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion\r
- fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion\r
- fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon\r
- fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon\r
- fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon\r
- fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl proton\r
- fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton\r
- fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton\r
-\r
- h->Fit(fstepTot,modfit.Data(),"",fRangeStep4[0],fRangeStep4[1]);//refit all\r
-\r
- //************************************* storing parameter to calculate the yields *******\r
- Int_t chpa=0;\r
- if(code==321) chpa=3;\r
- if(code==2212) chpa=6;\r
- for(Int_t j=0;j<3;j++) {\r
- fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
- fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
- }\r
-\r
- DrawFitFunction(fstepTot);\r
- CalcResidual(h,fstepTot,gres);\r
- return;\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
+ 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
+ Int_t code=TMath::Abs(signedcode);\r
+ GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\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->SetParLimits(1,fRangeStep4[0],fRangeStep4[1]); //limits mean pion (dummy)\r
+ fstep1->SetParLimits(2,expPionSig*fLimitsOnSigmaPion[0],expPionSig*fLimitsOnSigmaPion[1]); //limits sigma pion\r
+\r
+ if(expPionSig>0) h->Fit(fstep1,modfit.Data(),"",mean+fRangeStep1[0],mean+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
+ fstep2 = new TF1("fstep2",DoubleGausStep,fRangeStep4[0],fRangeStep4[1],6);\r
+ fstep2->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion\r
+ 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->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
+\r
+ if(expKaonSig>0) h->Fit(fstep2,modfit.Data(),"",expKaonMean+fRangeStep2[0],expKaonMean+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
+ fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion\r
+ fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion\r
+ fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion\r
+ fstep3->FixParameter(3,fstep2->GetParameter(3)); //fixed ampl kaon\r
+ 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->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
+\r
+ if(expProtonSig>0) h->Fit(fstep3,modfit.Data(),"",expProtonMean+fRangeStep3[0],expProtonMean+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
+ fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
+ fstepTot->SetLineColor(1);\r
+\r
+ Double_t initialParametersStepTot[9];\r
+ initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
+ initialParametersStepTot[1] = fstep1->GetParameter(1);\r
+ initialParametersStepTot[2] = fstep1->GetParameter(2);\r
+\r
+ initialParametersStepTot[3] = fstep2->GetParameter(3);//second gaussian\r
+ initialParametersStepTot[4] = fstep2->GetParameter(4);\r
+ initialParametersStepTot[5] = fstep2->GetParameter(5);\r
+\r
+ initialParametersStepTot[6] = fstep3->GetParameter(6);//third gaussian\r
+ initialParametersStepTot[7] = fstep3->GetParameter(7);\r
+ initialParametersStepTot[8] = fstep3->GetParameter(8); \r
+\r
+ fstepTot->SetParameters(initialParametersStepTot); //initial parameters\r
+ fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl pion\r
+ fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion\r
+ fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion\r
+ fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon\r
+ fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon\r
+ fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon\r
+ fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl proton\r
+ fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton\r
+ fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton\r
+\r
+ h->Fit(fstepTot,modfit.Data(),"",fRangeStep4[0],fRangeStep4[1]);//refit all\r
+\r
+ //************************************* storing parameter to calculate the yields *******\r
+ Int_t chpa=0;\r
+ if(code==321) chpa=3;\r
+ if(code==2212) chpa=6;\r
+ for(Int_t j=0;j<3;j++) {\r
+ fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
+ fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
+ }\r
+\r
+ DrawFitFunction(fstepTot);\r
+ CalcResidual(h,fstepTot,gres);\r
+ return;\r
}\r
\r
//________________________________________________________\r
void AliITSsadEdxFitter::DoFitProton(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, 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
- Int_t code=TMath::Abs(signedcode);\r
- GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\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->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
-\r
- if(expPionSig>0) h->Fit(fstep1,modfit,"",mean+fRangeStep1[0],mean+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->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
-\r
- if(expProtonSig>0) h->Fit(fstep2,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+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
- fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
- fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion\r
- fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion\r
- fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion\r
- fstep3->FixParameter(6,fstep2->GetParameter(0)); //fixed ampl proton\r
- 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->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
- else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);\r
-\r
- printf("___________________________________________________________________ Final Step: refit all\n");\r
- fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
- fstepTot->SetLineColor(1);\r
-\r
- Double_t initialParametersStepTot[9];\r
- initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
- initialParametersStepTot[1] = fstep1->GetParameter(1);\r
- initialParametersStepTot[2] = fstep1->GetParameter(2);\r
-\r
- initialParametersStepTot[3] = fstep3->GetParameter(3);//second gaussian\r
- initialParametersStepTot[4] = fstep3->GetParameter(4);\r
- initialParametersStepTot[5] = fstep3->GetParameter(5);\r
-\r
- initialParametersStepTot[6] = fstep2->GetParameter(0);//third gaussian\r
- initialParametersStepTot[7] = fstep2->GetParameter(1);\r
- initialParametersStepTot[8] = fstep2->GetParameter(2); \r
-\r
- fstepTot->SetParameters(initialParametersStepTot); //initial parameters\r
- fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl pion\r
- fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion\r
- fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion\r
- fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon\r
- fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon\r
- fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon\r
- fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl proton\r
- fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton\r
- fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton\r
-\r
- h->Fit(fstepTot,modfit,"",fRangeStep4[0],fRangeStep4[1]);//refit all\r
-\r
- //************************************* storing parameter to calculate the yields *******\r
- Int_t chpa=0;\r
- if(code==321) chpa=3;\r
- if(code==2212) chpa=6;\r
- for(Int_t j=0;j<3;j++) {\r
- fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
- fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
- }\r
-\r
- DrawFitFunction(fstepTot);\r
- CalcResidual(h,fstepTot,gres);\r
- return;\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
+ Int_t code=TMath::Abs(signedcode);\r
+ GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\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->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
+\r
+ if(expPionSig>0) h->Fit(fstep1,modfit,"",mean+fRangeStep1[0],mean+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->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
+\r
+ if(expProtonSig>0) h->Fit(fstep2,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+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
+ fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
+ fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl pion\r
+ fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean pion\r
+ fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma pion\r
+ fstep3->FixParameter(6,fstep2->GetParameter(0)); //fixed ampl proton\r
+ 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->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
+ else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);\r
+\r
+ printf("___________________________________________________________________ Final Step: refit all\n");\r
+ fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
+ fstepTot->SetLineColor(1);\r
+\r
+ Double_t initialParametersStepTot[9];\r
+ initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
+ initialParametersStepTot[1] = fstep1->GetParameter(1);\r
+ initialParametersStepTot[2] = fstep1->GetParameter(2);\r
+\r
+ initialParametersStepTot[3] = fstep3->GetParameter(3);//second gaussian\r
+ initialParametersStepTot[4] = fstep3->GetParameter(4);\r
+ initialParametersStepTot[5] = fstep3->GetParameter(5);\r
+\r
+ initialParametersStepTot[6] = fstep2->GetParameter(0);//third gaussian\r
+ initialParametersStepTot[7] = fstep2->GetParameter(1);\r
+ initialParametersStepTot[8] = fstep2->GetParameter(2); \r
+\r
+ fstepTot->SetParameters(initialParametersStepTot); //initial parameters\r
+ fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl pion\r
+ fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion\r
+ fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion\r
+ fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon\r
+ fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon\r
+ fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon\r
+ fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl proton\r
+ fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton\r
+ fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton\r
+\r
+ h->Fit(fstepTot,modfit,"",fRangeStep4[0],fRangeStep4[1]);//refit all\r
+\r
+ //************************************* storing parameter to calculate the yields *******\r
+ Int_t chpa=0;\r
+ if(code==321) chpa=3;\r
+ if(code==2212) chpa=6;\r
+ for(Int_t j=0;j<3;j++) {\r
+ fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
+ fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
+ }\r
+\r
+ DrawFitFunction(fstepTot);\r
+ CalcResidual(h,fstepTot,gres);\r
+ return;\r
}\r
\r
//________________________________________________________\r
void AliITSsadEdxFitter::DoFitProtonFirst(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, 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
- Int_t code=TMath::Abs(signedcode);\r
- GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\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->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
-\r
- if(expProtonSig>0) h->Fit(fstep1,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+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
- fstep2 = new TF1("step2",DoubleGausStep,fRangeStep4[0],fRangeStep4[1],6);\r
- fstep2->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl proton \r
- 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->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
-\r
- if(expPionSig>0) h->Fit(fstep2,modfit,"",mean+fRangeStep1[0],mean+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
- fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
- fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl proton\r
- fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean proton\r
- fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma proton\r
- fstep3->FixParameter(3,fstep2->GetParameter(3)); //fixed ampl pion\r
- 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->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
- else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);\r
-\r
- printf("___________________________________________________________________ Final Step: refit all\n");\r
- fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
- fstepTot->SetLineColor(1);\r
-\r
- Double_t initialParametersStepTot[9];\r
- initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
- initialParametersStepTot[1] = fstep1->GetParameter(1);\r
- initialParametersStepTot[2] = fstep1->GetParameter(2);\r
-\r
- initialParametersStepTot[3] = fstep3->GetParameter(6);//second gaussian\r
- initialParametersStepTot[4] = fstep3->GetParameter(7);\r
- initialParametersStepTot[5] = fstep3->GetParameter(8);\r
-\r
- initialParametersStepTot[6] = fstep2->GetParameter(3);//third gaussian\r
- initialParametersStepTot[7] = fstep2->GetParameter(4);\r
- initialParametersStepTot[8] = fstep2->GetParameter(5); \r
-\r
- fstepTot->SetParameters(initialParametersStepTot); //initial parameters\r
- fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl proton\r
- fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion\r
- fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion\r
- fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon\r
- fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon\r
- fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon\r
- fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl pion\r
- fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton\r
- fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton\r
-\r
- h->Fit(fstepTot,modfit,"",fRangeStep4[0],fRangeStep4[1]);//refit all\r
-\r
- //************************************* storing parameter to calculate the yields *******\r
- Int_t chpa=0;\r
- if(code==321) chpa=3;\r
- if(code==211) chpa=6;\r
- for(Int_t j=0;j<3;j++) {\r
- fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
- fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
- }\r
-\r
- DrawFitFunction(fstepTot);\r
- CalcResidual(h,fstepTot,gres);\r
- return;\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
+ Int_t code=TMath::Abs(signedcode);\r
+ GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\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->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
+\r
+ if(expProtonSig>0) h->Fit(fstep1,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+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
+ fstep2 = new TF1("step2",DoubleGausStep,fRangeStep4[0],fRangeStep4[1],6);\r
+ fstep2->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl proton \r
+ 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->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
+\r
+ if(expPionSig>0) h->Fit(fstep2,modfit,"",mean+fRangeStep1[0],mean+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
+ fstep3= new TF1("fstep3",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
+ fstep3->FixParameter(0,fstep1->GetParameter(0)); //fixed ampl proton\r
+ fstep3->FixParameter(1,fstep1->GetParameter(1)); //fixed mean proton\r
+ fstep3->FixParameter(2,fstep1->GetParameter(2)); //fixed sigma proton\r
+ fstep3->FixParameter(3,fstep2->GetParameter(3)); //fixed ampl pion\r
+ 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->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
+ else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);\r
+\r
+ printf("___________________________________________________________________ Final Step: refit all\n");\r
+ fstepTot = new TF1("funztot",FinalGausStep,fRangeStep4[0],fRangeStep4[1],9);\r
+ fstepTot->SetLineColor(1);\r
+\r
+ Double_t initialParametersStepTot[9];\r
+ initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
+ initialParametersStepTot[1] = fstep1->GetParameter(1);\r
+ initialParametersStepTot[2] = fstep1->GetParameter(2);\r
+\r
+ initialParametersStepTot[3] = fstep3->GetParameter(6);//second gaussian\r
+ initialParametersStepTot[4] = fstep3->GetParameter(7);\r
+ initialParametersStepTot[5] = fstep3->GetParameter(8);\r
+\r
+ initialParametersStepTot[6] = fstep2->GetParameter(3);//third gaussian\r
+ initialParametersStepTot[7] = fstep2->GetParameter(4);\r
+ initialParametersStepTot[8] = fstep2->GetParameter(5); \r
+\r
+ fstepTot->SetParameters(initialParametersStepTot); //initial parameters\r
+ fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1); //tolerance limits ampl proton\r
+ fstepTot->FixParameter(1,initialParametersStepTot[1]); //fixed mean pion\r
+ fstepTot->FixParameter(2,initialParametersStepTot[2]); //fixed sigma pion\r
+ fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1); //tolerance limits ampl kaon\r
+ fstepTot->FixParameter(4,initialParametersStepTot[4]); //fixed mean kaon\r
+ fstepTot->FixParameter(5,initialParametersStepTot[5]); //fixed sigma kaon\r
+ fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1); //tolerance limits ampl pion\r
+ fstepTot->FixParameter(7,initialParametersStepTot[7]); //fixed mean proton\r
+ fstepTot->FixParameter(8,initialParametersStepTot[8]); //fixed sigma proton\r
+\r
+ h->Fit(fstepTot,modfit,"",fRangeStep4[0],fRangeStep4[1]);//refit all\r
+\r
+ //************************************* storing parameter to calculate the yields *******\r
+ Int_t chpa=0;\r
+ if(code==321) chpa=3;\r
+ if(code==211) chpa=6;\r
+ for(Int_t j=0;j<3;j++) {\r
+ fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
+ fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
+ }\r
+\r
+ DrawFitFunction(fstepTot);\r
+ CalcResidual(h,fstepTot,gres);\r
+ return;\r
}\r
\r
\r
//________________________________________________________\r
void AliITSsadEdxFitter::DoFitOnePeak(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc){\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
- Int_t code=TMath::Abs(signedcode);\r
- GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\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->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
-\r
- if(expProtonSig>0) h->Fit(fstep1,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);//fit\r
- else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);\r
-\r
- fstep1->SetLineColor(1);\r
- fstep1->Draw("same");\r
-\r
- fFitpar[0] = fstep1->GetParameter(0);\r
- fFitparErr[0] = fstep1->GetParError(0);\r
- return;\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
+ Int_t code=TMath::Abs(signedcode);\r
+ GetInitialParam(h,mc,code,bin,pt,ampl,mean,expKaonMean,expProtonMean,expPionSig,expKaonSig,expProtonSig);\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->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
+\r
+ if(expProtonSig>0) h->Fit(fstep1,modfit,"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);//fit\r
+ else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);\r
+\r
+ fstep1->SetLineColor(1);\r
+ fstep1->Draw("same");\r
+\r
+ fFitpar[0] = fstep1->GetParameter(0);\r
+ fFitparErr[0] = fstep1->GetParError(0);\r
+ return;\r
}\r
\r
//________________________________________________________\r
void AliITSsadEdxFitter::DoFitTail(TH1F *h, Int_t bin, Int_t signedcode){\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
- // 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
-\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(3,1.2);\r
- fstep1->SetParameter(4,10.);\r
-\r
- fstep1->SetParLimits(0,0,ampl*1.2);\r
- fstep1->SetParLimits(1,-3.5,3.5);\r
- fstep1->SetParLimits(2,0.1,0.25);\r
- 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
-\r
- printf("\n___________________________________________________________________\n Second Step: kaons\n\n"); \r
- fstep2 = new TF1("fstep2",DoubleGausTail,-3.5,3.5,10);\r
- fstep2->FixParameter(0,fstep1->GetParameter(0));//fixed\r
- fstep2->FixParameter(1,fstep1->GetParameter(1));\r
- fstep2->FixParameter(2,fstep1->GetParameter(2));\r
- fstep2->FixParameter(3,fstep1->GetParameter(3));\r
- fstep2->FixParameter(4,fstep1->GetParameter(4));\r
-\r
- fstep2->SetParameter(5,fstep1->GetParameter(0)/8);//initial\r
- //fstep2->SetParameter(6,CalcP(code,322,bin));\r
- fstep2->SetParameter(7,0.145);\r
- fstep2->FixParameter(8,1.2);\r
- fstep2->SetParameter(9,13.);\r
-\r
- fstep2->SetParLimits(5,fstep1->GetParameter(0)/100,fstep1->GetParameter(0));//limits\r
- fstep2->SetParLimits(6,-3.5,3.5);\r
- fstep2->SetParLimits(7,0.12,0.2);\r
- fstep2->SetParLimits(9,9.,20.);\r
- if(bin<9) fstep2->SetParLimits(9,13.,25.);\r
-\r
- //h->Fit(fstep2,"M0R+","",CalcP(code,321,bin)-0.3,CalcP(code,321,bin)+0.3);//second fit\r
- if(bin<6 || bin>12) for(Int_t npar=5;npar<10;npar++) fstep2->FixParameter(npar,-0.0000000001);\r
-\r
- printf("\n____________________________________________________________________\n Third Step: protons \n\n");\r
- fstep3= new TF1("fstep3",FinalGausTail,-3.5,3.5,15);\r
- fstep3->FixParameter(0,fstep1->GetParameter(0));//fixed\r
- fstep3->FixParameter(1,fstep1->GetParameter(1));\r
- fstep3->FixParameter(2,fstep1->GetParameter(2));\r
- fstep3->FixParameter(3,fstep1->GetParameter(3));\r
- fstep3->FixParameter(4,fstep1->GetParameter(4));\r
- fstep3->FixParameter(5,fstep2->GetParameter(5));\r
- fstep3->FixParameter(6,fstep2->GetParameter(6));\r
- fstep3->FixParameter(7,fstep2->GetParameter(7));\r
- fstep3->FixParameter(8,fstep2->GetParameter(8));\r
- fstep3->FixParameter(9,fstep2->GetParameter(9));\r
-\r
- fstep3->SetParameter(10,fstep2->GetParameter(0)/8);//initial\r
- //fstep3->SetParameter(11,CalcP(code,2212,bin));\r
- fstep3->SetParameter(12,0.145);\r
- fstep3->FixParameter(13,1.2);\r
- fstep3->SetParameter(14,10.);\r
-\r
- fstep3->SetParLimits(10,fstep2->GetParameter(0)/100,fstep2->GetParameter(0));//limits\r
- fstep3->SetParLimits(11,-3.5,3.5);\r
- fstep3->SetParLimits(12,0.12,0.2);\r
- fstep3->SetParLimits(14,11.,25.);\r
-\r
- printf("\n_____________________________________________________________________\n Final Step: refit all \n\n"); \r
- fstepTot = new TF1("funztot",FinalGausTail,-3.5,3.5,15);\r
- fstepTot->SetLineColor(1);\r
-\r
- Double_t initialParametersStepTot[15];\r
- initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
- initialParametersStepTot[1] = fstep1->GetParameter(1);\r
- initialParametersStepTot[2] = fstep1->GetParameter(2);\r
- initialParametersStepTot[3] = fstep1->GetParameter(3);\r
- initialParametersStepTot[4] = fstep1->GetParameter(4);\r
-\r
- initialParametersStepTot[5] = fstep2->GetParameter(5);//second gaussian\r
- initialParametersStepTot[6] = fstep2->GetParameter(6);\r
- initialParametersStepTot[7] = fstep2->GetParameter(7);\r
- initialParametersStepTot[8] = fstep2->GetParameter(8);\r
- initialParametersStepTot[9] = fstep2->GetParameter(9);\r
-\r
- initialParametersStepTot[10] = fstep3->GetParameter(10);//third gaussian\r
- initialParametersStepTot[11] = fstep3->GetParameter(11);\r
- initialParametersStepTot[12] = fstep3->GetParameter(12); \r
- initialParametersStepTot[13] = fstep3->GetParameter(13); \r
- initialParametersStepTot[14] = fstep3->GetParameter(14); \r
-\r
- fstepTot->SetParameters(initialParametersStepTot);//initial parameter\r
-\r
-\r
- fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1);//tollerance limit\r
- fstepTot->SetParLimits(1,initialParametersStepTot[1]*0.9,initialParametersStepTot[1]*1.1);\r
- fstepTot->SetParLimits(2,initialParametersStepTot[2]*0.9,initialParametersStepTot[2]*1.1);\r
- fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1);\r
- fstepTot->SetParLimits(4,initialParametersStepTot[4]*0.9,initialParametersStepTot[4]*1.1);\r
- fstepTot->SetParLimits(5,initialParametersStepTot[5]*0.9,initialParametersStepTot[5]*1.1);\r
- fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1);\r
- fstepTot->SetParLimits(7,initialParametersStepTot[7]*0.9,initialParametersStepTot[7]*1.1);\r
- fstepTot->SetParLimits(8,initialParametersStepTot[8]*0.9,initialParametersStepTot[8]*1.1);\r
- fstepTot->SetParLimits(9,initialParametersStepTot[9]*0.9,initialParametersStepTot[9]*1.1);\r
- fstepTot->SetParLimits(10,initialParametersStepTot[10]*0.9,initialParametersStepTot[10]*1.1);\r
- fstepTot->SetParLimits(11,initialParametersStepTot[11]*0.9,initialParametersStepTot[11]*1.1);\r
- fstepTot->SetParLimits(12,initialParametersStepTot[12]*0.9,initialParametersStepTot[12]*1.1);\r
- fstepTot->SetParLimits(13,initialParametersStepTot[13]*0.9,initialParametersStepTot[13]*1.1);\r
- fstepTot->SetParLimits(14,initialParametersStepTot[14]*0.9,initialParametersStepTot[14]*1.1);\r
-\r
- if(bin<9) for(Int_t npar=10;npar<15;npar++) fstepTot->FixParameter(npar,-0.00000000001);\r
- h->Fit(fstepTot,modfit,"",-3.5,3.5); //refit all\r
-\r
-\r
- //************************************* storing parameter to calculate the yields *******\r
- Int_t chpa=0;\r
- if(code==321) chpa=5;\r
- if(code==2212) chpa=10;\r
- for(Int_t j=0;j<5;j++) {\r
- fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
- fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
- }\r
-\r
- DrawFitFunction(fstepTot);\r
- return;\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
+ // 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
+\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(3,1.2);\r
+ fstep1->SetParameter(4,10.);\r
+\r
+ fstep1->SetParLimits(0,0,ampl*1.2);\r
+ fstep1->SetParLimits(1,-3.5,3.5);\r
+ fstep1->SetParLimits(2,0.1,0.25);\r
+ 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
+\r
+ printf("\n___________________________________________________________________\n Second Step: kaons\n\n"); \r
+ fstep2 = new TF1("fstep2",DoubleGausTail,-3.5,3.5,10);\r
+ fstep2->FixParameter(0,fstep1->GetParameter(0));//fixed\r
+ fstep2->FixParameter(1,fstep1->GetParameter(1));\r
+ fstep2->FixParameter(2,fstep1->GetParameter(2));\r
+ fstep2->FixParameter(3,fstep1->GetParameter(3));\r
+ fstep2->FixParameter(4,fstep1->GetParameter(4));\r
+\r
+ fstep2->SetParameter(5,fstep1->GetParameter(0)/8);//initial\r
+ //fstep2->SetParameter(6,CalcP(code,322,bin));\r
+ fstep2->SetParameter(7,0.145);\r
+ fstep2->FixParameter(8,1.2);\r
+ fstep2->SetParameter(9,13.);\r
+\r
+ fstep2->SetParLimits(5,fstep1->GetParameter(0)/100,fstep1->GetParameter(0));//limits\r
+ fstep2->SetParLimits(6,-3.5,3.5);\r
+ fstep2->SetParLimits(7,0.12,0.2);\r
+ fstep2->SetParLimits(9,9.,20.);\r
+ if(bin<9) fstep2->SetParLimits(9,13.,25.);\r
+\r
+ //h->Fit(fstep2,"M0R+","",CalcP(code,321,bin)-0.3,CalcP(code,321,bin)+0.3);//second fit\r
+ if(bin<6 || bin>12) for(Int_t npar=5;npar<10;npar++) fstep2->FixParameter(npar,-0.0000000001);\r
+\r
+ printf("\n____________________________________________________________________\n Third Step: protons \n\n");\r
+ fstep3= new TF1("fstep3",FinalGausTail,-3.5,3.5,15);\r
+ fstep3->FixParameter(0,fstep1->GetParameter(0));//fixed\r
+ fstep3->FixParameter(1,fstep1->GetParameter(1));\r
+ fstep3->FixParameter(2,fstep1->GetParameter(2));\r
+ fstep3->FixParameter(3,fstep1->GetParameter(3));\r
+ fstep3->FixParameter(4,fstep1->GetParameter(4));\r
+ fstep3->FixParameter(5,fstep2->GetParameter(5));\r
+ fstep3->FixParameter(6,fstep2->GetParameter(6));\r
+ fstep3->FixParameter(7,fstep2->GetParameter(7));\r
+ fstep3->FixParameter(8,fstep2->GetParameter(8));\r
+ fstep3->FixParameter(9,fstep2->GetParameter(9));\r
+\r
+ fstep3->SetParameter(10,fstep2->GetParameter(0)/8);//initial\r
+ //fstep3->SetParameter(11,CalcP(code,2212,bin));\r
+ fstep3->SetParameter(12,0.145);\r
+ fstep3->FixParameter(13,1.2);\r
+ fstep3->SetParameter(14,10.);\r
+\r
+ fstep3->SetParLimits(10,fstep2->GetParameter(0)/100,fstep2->GetParameter(0));//limits\r
+ fstep3->SetParLimits(11,-3.5,3.5);\r
+ fstep3->SetParLimits(12,0.12,0.2);\r
+ fstep3->SetParLimits(14,11.,25.);\r
+\r
+ printf("\n_____________________________________________________________________\n Final Step: refit all \n\n"); \r
+ fstepTot = new TF1("funztot",FinalGausTail,-3.5,3.5,15);\r
+ fstepTot->SetLineColor(1);\r
+\r
+ Double_t initialParametersStepTot[15];\r
+ initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian\r
+ initialParametersStepTot[1] = fstep1->GetParameter(1);\r
+ initialParametersStepTot[2] = fstep1->GetParameter(2);\r
+ initialParametersStepTot[3] = fstep1->GetParameter(3);\r
+ initialParametersStepTot[4] = fstep1->GetParameter(4);\r
+\r
+ initialParametersStepTot[5] = fstep2->GetParameter(5);//second gaussian\r
+ initialParametersStepTot[6] = fstep2->GetParameter(6);\r
+ initialParametersStepTot[7] = fstep2->GetParameter(7);\r
+ initialParametersStepTot[8] = fstep2->GetParameter(8);\r
+ initialParametersStepTot[9] = fstep2->GetParameter(9);\r
+\r
+ initialParametersStepTot[10] = fstep3->GetParameter(10);//third gaussian\r
+ initialParametersStepTot[11] = fstep3->GetParameter(11);\r
+ initialParametersStepTot[12] = fstep3->GetParameter(12); \r
+ initialParametersStepTot[13] = fstep3->GetParameter(13); \r
+ initialParametersStepTot[14] = fstep3->GetParameter(14); \r
+\r
+ fstepTot->SetParameters(initialParametersStepTot);//initial parameter\r
+\r
+\r
+ fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1);//tollerance limit\r
+ fstepTot->SetParLimits(1,initialParametersStepTot[1]*0.9,initialParametersStepTot[1]*1.1);\r
+ fstepTot->SetParLimits(2,initialParametersStepTot[2]*0.9,initialParametersStepTot[2]*1.1);\r
+ fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1);\r
+ fstepTot->SetParLimits(4,initialParametersStepTot[4]*0.9,initialParametersStepTot[4]*1.1);\r
+ fstepTot->SetParLimits(5,initialParametersStepTot[5]*0.9,initialParametersStepTot[5]*1.1);\r
+ fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1);\r
+ fstepTot->SetParLimits(7,initialParametersStepTot[7]*0.9,initialParametersStepTot[7]*1.1);\r
+ fstepTot->SetParLimits(8,initialParametersStepTot[8]*0.9,initialParametersStepTot[8]*1.1);\r
+ fstepTot->SetParLimits(9,initialParametersStepTot[9]*0.9,initialParametersStepTot[9]*1.1);\r
+ fstepTot->SetParLimits(10,initialParametersStepTot[10]*0.9,initialParametersStepTot[10]*1.1);\r
+ fstepTot->SetParLimits(11,initialParametersStepTot[11]*0.9,initialParametersStepTot[11]*1.1);\r
+ fstepTot->SetParLimits(12,initialParametersStepTot[12]*0.9,initialParametersStepTot[12]*1.1);\r
+ fstepTot->SetParLimits(13,initialParametersStepTot[13]*0.9,initialParametersStepTot[13]*1.1);\r
+ fstepTot->SetParLimits(14,initialParametersStepTot[14]*0.9,initialParametersStepTot[14]*1.1);\r
+\r
+ if(bin<9) for(Int_t npar=10;npar<15;npar++) fstepTot->FixParameter(npar,-0.00000000001);\r
+ h->Fit(fstepTot,modfit,"",-3.5,3.5); //refit all\r
+\r
+\r
+ //************************************* storing parameter to calculate the yields *******\r
+ Int_t chpa=0;\r
+ if(code==321) chpa=5;\r
+ if(code==2212) chpa=10;\r
+ for(Int_t j=0;j<3;j++) {\r
+ fFitpar[j] = fstepTot->GetParameter(j+chpa);\r
+ fFitparErr[j] = fstepTot->GetParError(j+chpa);\r
+ }\r
+\r
+ DrawFitFunction(fstepTot);\r
+ return;\r
}\r
\r
//________________________________________________________\r
void AliITSsadEdxFitter::FillHisto(TH1F *hsps, Int_t bin, Float_t binsize, Int_t code){\r
- // fill the spectra histo calculating the yield\r
- // first bin has to be 1\r
- Double_t yield = 0;\r
- Double_t err = 0;\r
- Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin); \r
-\r
- if(IsGoodBin(bin-1,code)) {\r
- yield = fFitpar[0] / ptbin / binsize;\r
- err = fFitparErr[0] / ptbin / binsize;\r
- }\r
-\r
- hsps->SetBinContent(bin,yield);\r
- hsps->SetBinError(bin,err);\r
- return;\r
+ // fill the spectra histo calculating the yield\r
+ // first bin has to be 1\r
+ Double_t yield = 0;\r
+ Double_t err = 0;\r
+ Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin); \r
+\r
+ if(IsGoodBin(bin-1,code)) {\r
+ yield = fFitpar[0] / ptbin / binsize;\r
+ err = fFitparErr[0] / ptbin / binsize;\r
+ }\r
+\r
+ hsps->SetBinContent(bin,yield);\r
+ hsps->SetBinError(bin,err);\r
+ return;\r
}\r
\r
//________________________________________________________\r
void AliITSsadEdxFitter::FillHistoMC(TH1F *hsps, Int_t bin, Int_t code, TH1F *h){\r
- // fill the spectra histo calculating the yield (using the MC truth)\r
- // first bin has to be 1\r
- Double_t yield = 0;\r
- Double_t erryield=0;\r
- Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin);\r
-\r
- if(IsGoodBin(bin-1,code)){\r
- yield = h->GetEntries() / ptbin;\r
- erryield=TMath::Sqrt(h->GetEntries()) / ptbin;\r
- }\r
- hsps->SetBinContent(bin,yield);\r
- hsps->SetBinError(bin,erryield);\r
- return;\r
+ // fill the spectra histo calculating the yield (using the MC truth)\r
+ // first bin has to be 1\r
+ Double_t yield = 0;\r
+ Double_t erryield=0;\r
+ Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin);\r
+\r
+ if(IsGoodBin(bin-1,code)){\r
+ yield = h->GetEntries() / ptbin;\r
+ erryield=TMath::Sqrt(h->GetEntries()) / ptbin;\r
+ }\r
+ hsps->SetBinContent(bin,yield);\r
+ hsps->SetBinError(bin,erryield);\r
+ return;\r
}\r
\r
//________________________________________________________\r
void AliITSsadEdxFitter::GetFitPar(Double_t *fitpar, Double_t *fitparerr) const {\r
- // getter of the fit parameters and the relative errors\r
- for(Int_t i=0;i<3;i++) {\r
- fitpar[i] = fFitpar[i];\r
- fitparerr[i] = fFitparErr[i];\r
- }\r
- return;\r
+ // getter of the fit parameters and the relative errors\r
+ for(Int_t i=0;i<3;i++) {\r
+ fitpar[i] = fFitpar[i];\r
+ fitparerr[i] = fFitparErr[i];\r
+ }\r
+ return;\r
}\r
\r
\r
//________________________________________________________\r
void AliITSsadEdxFitter::PrintAll() const{\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
- printf("Range F = %f %f\n",fRangeStep4[0],fRangeStep4[1]);\r
- printf(" Sigma1 = %f %f\n",fLimitsOnSigmaPion[0],fLimitsOnSigmaPion[1]);\r
- printf(" Sigma2 = %f %f\n",fLimitsOnSigmaKaon[0],fLimitsOnSigmaKaon[1]);\r
- printf(" Sigma3 = %f %f\n",fLimitsOnSigmaProton[0],fLimitsOnSigmaProton[1]);\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
+ printf("Range F = %f %f\n",fRangeStep4[0],fRangeStep4[1]);\r
+ printf(" Sigma1 = %f %f\n",fLimitsOnSigmaPion[0],fLimitsOnSigmaPion[1]);\r
+ printf(" Sigma2 = %f %f\n",fLimitsOnSigmaKaon[0],fLimitsOnSigmaKaon[1]);\r
+ printf(" Sigma3 = %f %f\n",fLimitsOnSigmaProton[0],fLimitsOnSigmaProton[1]);\r
}\r
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff