# SHLIBS - Shared Libraries and objects for linking (Executables only) #
#--------------------------------------------------------------------------------#
-set ( SRCS SPECTRA/AliProtonAnalysisBase.cxx SPECTRA/AliProtonAnalysis.cxx SPECTRA/AliProtonQAAnalysis.cxx SPECTRA/AliAnalysisTaskProtons.cxx SPECTRA/AliAnalysisTaskProtonsQA.cxx SPECTRA/AliAnalysisTaskChargedHadronSpectra.cxx SPECTRA/AliAnalysisTaskCheckV0.cxx SPECTRA/AliAnalysisTaskCheckCascade.cxx SPECTRA/AliAnalysisTaskCheckPerformanceCascade.cxx SPECTRA/AliAnalysisTaskStrange.cxx SPECTRA/AliAnalysisTaskPerformanceStrange.cxx SPECTRA/AliProtonFeedDownAnalysisTask.cxx SPECTRA/AliProtonFeedDownAnalysis.cxx SPECTRA/AliProtonAbsorptionCorrection.cxx SPECTRA/AliProtonSpectraCorrection.cxx SPECTRA/AliProtonCorrectionAnalysisTask.cxx SPECTRA/AliAnalysisCentralCutMC.cxx SPECTRA/AliAnalysisCentralCutESD.cxx SPECTRA/AliAnalysisCentralCutEvtMC.cxx SPECTRA/AliAnalysisCentralCutEvtESD.cxx SPECTRA/AliAnalysisCentralExtrapolate.cxx SPECTRA/AliAnalysisTaskCentral.cxx SPECTRA/AliAnalysisTaskSEITSsaSpectra.cxx SPECTRA/Fit/AliBWTools.cxx SPECTRA/Fit/AliBWFunc.cxx SPECTRA/Fit/AliLatexTable.cxx SPECTRA/AliAnalysisChargedHadronSpectraITSTruncatedMeanTask.cxx)
+set ( SRCS SPECTRA/AliProtonAnalysisBase.cxx SPECTRA/AliProtonAnalysis.cxx SPECTRA/AliProtonQAAnalysis.cxx SPECTRA/AliAnalysisTaskProtons.cxx SPECTRA/AliAnalysisTaskProtonsQA.cxx SPECTRA/AliAnalysisTaskChargedHadronSpectra.cxx SPECTRA/AliAnalysisTaskCheckV0.cxx SPECTRA/AliAnalysisTaskCheckCascade.cxx SPECTRA/AliAnalysisTaskCheckPerformanceCascade.cxx SPECTRA/AliAnalysisTaskStrange.cxx SPECTRA/AliAnalysisTaskPerformanceStrange.cxx SPECTRA/AliProtonFeedDownAnalysisTask.cxx SPECTRA/AliProtonFeedDownAnalysis.cxx SPECTRA/AliProtonAbsorptionCorrection.cxx SPECTRA/AliProtonSpectraCorrection.cxx SPECTRA/AliProtonCorrectionAnalysisTask.cxx SPECTRA/AliAnalysisCentralCutMC.cxx SPECTRA/AliAnalysisCentralCutESD.cxx SPECTRA/AliAnalysisCentralCutEvtMC.cxx SPECTRA/AliAnalysisCentralCutEvtESD.cxx SPECTRA/AliAnalysisCentralExtrapolate.cxx SPECTRA/AliAnalysisTaskCentral.cxx SPECTRA/AliAnalysisTaskSEITSsaSpectra.cxx SPECTRA/AliITSsadEdxFitter.cxx SPECTRA/Fit/AliBWTools.cxx SPECTRA/Fit/AliBWFunc.cxx SPECTRA/Fit/AliLatexTable.cxx SPECTRA/AliAnalysisChargedHadronSpectraITSTruncatedMeanTask.cxx)
string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
#pragma link C++ class AliAnalysisCentralExtrapolate+;
#pragma link C++ class AliAnalysisTaskCentral+;
#pragma link C++ class AliAnalysisTaskSEITSsaSpectra+;
+#pragma link C++ class AliITSsadEdxFitter+;
#pragma link C++ class AliBWTools+;
#pragma link C++ class AliBWFunc+;
#pragma link C++ class AliLatexTable+;
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id: $ */
+
+///////////////////////////////////////////////////////////////////////
+// Class with the fits algorithms to be used in the identified //
+// spectra analysis using the ITS in stand-alone mode //
+// Author: E.Biolcati,biolcati@to.infn.it //
+///////////////////////////////////////////////////////////////////////
+
+#include <Riostream.h>
+#include <TLatex.h>
+#include <TH1F.h>
+#include <TF1.h>
+#include <TH1D.h>
+#include <TLine.h>
+#include <TH2F.h>
+#include <TMath.h>
+#include <TGraph.h>
+#include <TGraphErrors.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+#include <TStyle.h>
+#include <Rtypes.h>
+#include "AliITSsadEdxFitter.h"
+
+
+ClassImp(AliITSsadEdxFitter)
+//______________________________________________________________________
+AliITSsadEdxFitter::AliITSsadEdxFitter():TObject(){
+ // standard constructor
+ for(Int_t i=0; i<5; i++) fFitpar[i] = 0.;
+ for(Int_t i=0; i<5; i++) fFitparErr[i] = 0.;
+ SetRangeStep1();
+ SetRangeStep2();
+ SetRangeStep3();
+ SetRangeFinalStep();
+ SetLimitsOnSigmaPion();
+ SetLimitsOnSigmaKaon();
+ SetLimitsOnSigmaProton();
+};
+
+//________________________________________________________
+Double_t AliITSsadEdxFitter::CalcSigma(Int_t code,Float_t x,Bool_t mc){
+ // calculate the sigma 12-ott-2010
+ Double_t p[2]={0.};
+ if(mc){
+ if(code==211){
+ p[0] = -1.20337e-04;
+ p[1] = 1.13060e-01;
+ }
+ else if(code==321 && x>0.15){
+ p[0] = -2.39631e-03;
+ p[1] = 1.15723e-01;
+ }
+ else if(code==2212 && x>0.3){
+ p[0] = -8.34576e-03;
+ p[1] = 1.34237e-01;
+ }
+ else return -1;
+ }
+ else{
+ if(code==211){
+ p[0] = -6.55200e-05;
+ p[1] = 1.26657e-01;
+ }
+ else if(code==321 && x>0.15){
+ p[0] = -6.22639e-04;
+ p[1] = 1.43289e-01;
+ }
+ else if(code==2212 && x>0.3){
+ p[0] = -2.13243e-03;
+ p[1] = 1.68614e-01;
+ }
+ else return -1;
+ }
+ return p[0]/(x*x)*TMath::Log(x)+p[1];
+}
+
+//_______________________________________________________
+Int_t AliITSsadEdxFitter::CalcMean(Int_t code, Float_t x, Float_t mean0, Float_t &mean1, Float_t &mean2){
+ // calculate the mean 12-ott-2010
+ Double_t p1[4]={0.};
+ Double_t p2[4]={0.};
+ if(code==211){
+ p1[0] = 1.77049;
+ p1[1] = -2.65469;
+ p2[0] = 0.890856;
+ p2[1] = -0.276719;
+ mean1 = mean0 + p1[0]+ p1[1]*x + p1[2]*x*x + p1[3]*x*x*x;
+ mean2 = mean1 + p2[0]+ p2[1]*x + p2[2]*x*x + p2[3]*x*x*x;
+ }
+ else if(code==321){
+ p2[0] = 1.57664;
+ p2[1] = -6.88635;
+ p2[2] = 18.702;
+ p2[3] = -16.3385;
+ mean1 = 0.;
+ mean2 = mean1 + p2[0]+ p2[1]*x + p2[2]*x*x + p2[3]*x*x*x;
+ }
+ else if(code==2212){
+ p1[0] = 4.24861;
+ p1[1] = -19.178;
+ p1[2] = 31.5947;
+ p1[3] = -18.178;
+ mean1 = mean0 + p1[0]+ p1[1]*x + p1[2]*x*x + p1[3]*x*x*x;
+ mean2 = 0.;
+ }
+ else return -1;
+ return 0;
+}
+
+//________________________________________________________
+Bool_t AliITSsadEdxFitter::IsGoodBin(Int_t bin,Int_t code){
+ //method to select the bins used for the analysis only
+ Bool_t retvalue=kTRUE;
+ TLine *l[2];
+ l[0]=new TLine(-2.1,0,2.,100.);
+ l[1]=new TLine(-1.9,120,2.,0.);
+ for(Int_t j=0;j<2;j++){
+ l[j]->SetLineColor(4);
+ l[j]->SetLineWidth(5);
+ if(code==211 && (bin>14 || bin<1)){
+ l[j]->Draw("same");
+ retvalue=kFALSE;
+ }
+ if(code==321 && (bin>12 || bin<5)){
+ l[j]->Draw("same");
+ retvalue=kFALSE;
+ }
+ if(code==2212 && (bin<8 || bin>16)){
+ l[j]->Draw("same");
+ retvalue=kFALSE;
+ }
+ }
+ return retvalue;
+}
+
+//________________________________________________________
+Double_t SingleGausTail(const Double_t *x, const Double_t *par){
+ //single gaussian with exponential tail
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
+ Double_t xx = x[0];
+ Double_t mean1 = par[1];
+ Double_t sigma1 = par[2];
+ Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);
+ Double_t sigmaSquare1 = sigma1 * sigma1;
+ Double_t xdiv1 = mean1 + par[3] * sigma1;
+ Double_t y1=0.0;
+ if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);
+ else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));
+ return y1;
+}
+
+//________________________________________________________
+Double_t DoubleGausTail(const Double_t *x, const Double_t *par){
+ //sum of two gaussians with exponential tail
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
+ Double_t xx = x[0];
+ Double_t mean1 = par[1];
+ Double_t sigma1 = par[2];
+ Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);
+ Double_t sigmaSquare1 = sigma1 * sigma1;
+ Double_t xdiv1 = mean1 + par[3] * sigma1;
+ Double_t y1=0.0;
+ if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);
+ else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));
+
+ Double_t mean2 = par[6];
+ Double_t sigma2 = par[7];
+ Double_t xNormSquare2 = (xx - mean2) * (xx - mean2);
+ Double_t sigmaSquare2 = sigma2 * sigma2;
+ Double_t xdiv2 = mean2 + par[8] * sigma2;
+ Double_t y2=0.0;
+ if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2);
+ else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8]));
+ return y1+y2;
+}
+
+//________________________________________________________
+Double_t FinalGausTail(const Double_t *x, const Double_t *par){
+ //sum of three gaussians with exponential tail
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
+ Double_t xx = x[0];
+ Double_t mean1 = par[1];
+ Double_t sigma1 = par[2];
+ Double_t xNormSquare1 = (xx - mean1) * (xx - mean1);
+ Double_t sigmaSquare1 = sigma1 * sigma1;
+ Double_t xdiv1 = mean1 + par[3] * sigma1;
+ Double_t y1=0.0;
+ if(xx < xdiv1) y1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNormSquare1 / sigmaSquare1);
+ else y1 = TMath::Exp(-par[4]*(xx-xdiv1)) * par[0] / (s2pi*par[2]) * TMath::Exp(-0.5*(par[3]*par[3]));
+
+ Double_t mean2 = par[6];
+ Double_t sigma2 = par[7];
+ Double_t xNormSquare2 = (xx - mean2) * (xx - mean2);
+ Double_t sigmaSquare2 = sigma2 * sigma2;
+ Double_t xdiv2 = mean2 + par[8] * sigma2;
+ Double_t y2=0.0;
+ if(xx < xdiv2) y2 = par[5]/(s2pi*par[7]) * TMath::Exp(-0.5 * xNormSquare2 / sigmaSquare2);
+ else y2 = TMath::Exp(-par[9]*(xx-xdiv2)) * par[5] / (s2pi*par[7]) * TMath::Exp(-0.5*(par[8]*par[8]));
+
+ Double_t mean3 = par[11];
+ Double_t sigma3 = par[12];
+ Double_t xNormSquare3 = (xx - mean3) * (xx - mean3);
+ Double_t sigmaSquare3 = sigma3 * sigma3;
+ Double_t xdiv3 = mean3 + par[13] * sigma3;
+ Double_t y3=0.0;
+ if(xx < xdiv3) y3 = par[10]/(s2pi*par[12]) * TMath::Exp(-0.5 * xNormSquare3 / sigmaSquare3);
+ else y3 = TMath::Exp(-par[14]*(xx-xdiv3)) * par[10] / (s2pi*par[12]) * TMath::Exp(-0.5*(par[13]*par[13]));
+ return y1+y2+y3;
+}
+
+//______________________________________________________________________
+void AliITSsadEdxFitter::CalcResidual(TH1F *h,TF1 *fun,TGraph *gres) const{
+ //code to calculate the difference fit function and histogram point (residual)
+ //to use as goodness test for the fit
+ Int_t ipt=0;
+ Double_t x=0.,yhis=0.,yfun=0.;
+ for(Int_t i=0;i<h->GetNbinsX();i++){
+ x=(h->GetBinLowEdge(i+2)+h->GetBinLowEdge(i+1))/2;
+ yfun=fun->Eval(x);
+ yhis=h->GetBinContent(i+1);
+ if(yhis>0. && yfun>0.) {
+ gres->SetPoint(ipt,x,(yhis-yfun)/yhis);
+ ipt++;
+ }
+ }
+ return;
+}
+
+
+//________________________________________________________
+Double_t SingleGausStep(const Double_t *x, const Double_t *par){
+ //single normalized gaussian
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
+ Double_t xx = x[0];
+ Double_t mean1 = par[1];
+ Double_t sigma1 = par[2];
+ Double_t xNorm1Square = (xx - mean1) * (xx - mean1);
+ Double_t sigma1Square = sigma1 * sigma1;
+ Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp(-0.5 * xNorm1Square / sigma1Square);
+ return step1;
+}
+
+//________________________________________________________
+Double_t DoubleGausStep(const Double_t *x, const Double_t *par){
+ //sum of two normalized gaussians
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
+ Double_t xx = x[0];
+ Double_t mean1 = par[1];
+ Double_t sigma1 = par[2];
+ Double_t xNorm1Square = (xx - mean1) * (xx - mean1);
+ Double_t sigma1Square = sigma1 * sigma1;
+ Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square );
+ Double_t mean2 = par[4];
+ Double_t sigma2 = par[5];
+ Double_t xNorm2Square = (xx - mean2) * (xx - mean2);
+ Double_t sigma2Square = sigma2 * sigma2;
+ Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square );
+ return step1+step2;
+}
+
+//________________________________________________________
+Double_t FinalGausStep(const Double_t *x, const Double_t *par){
+ //sum of three normalized gaussians
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
+ Double_t xx = x[0];
+ Double_t mean1 = par[1];
+ Double_t sigma1 = par[2];
+ Double_t xNorm1Square = (xx - mean1) * (xx - mean1);
+ Double_t sigma1Square = sigma1 * sigma1;
+ Double_t step1 = par[0]/(s2pi*par[2]) * TMath::Exp( - 0.5 * xNorm1Square / sigma1Square );
+ Double_t mean2 = par[4];
+ Double_t sigma2 = par[5];
+ Double_t xNorm2Square = (xx - mean2) * (xx - mean2);
+ Double_t sigma2Square = sigma2 * sigma2;
+ Double_t step2 = par[3]/(s2pi*par[5]) * TMath::Exp( - 0.5 * xNorm2Square / sigma2Square );
+ Double_t mean3 = par[7];
+ Double_t sigma3 = par[8];
+ Double_t xNorm3Square = (xx - mean3) * (xx - mean3);
+ Double_t sigma3Square = sigma3 * sigma3;
+ Double_t step3 = par[6]/(s2pi*par[8]) * TMath::Exp( - 0.5 * xNorm3Square / sigma3Square);
+ return step1+step2+step3;
+}
+
+//______________________________________________________________________
+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{
+ //gaussian with an exponential tail on the right side
+ Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi()));
+ Double_t returnvalue=0.0;
+ 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);
+ if (x<mean+cut*rms) returnvalue=TMath::Exp(-1.0*(x-mean)*(x-mean)/(2.0*rms*rms))*factor/TMath::Abs(rms);
+ else returnvalue=TMath::Exp(slope*(mean+cut*rms-x))*TMath::Exp(-cut*cut*0.5)*factor/TMath::Abs(rms);
+ return c*returnvalue/n;
+}
+
+
+//______________________________________________________________________
+Double_t AliITSsadEdxFitter::GausOnBackground(const Double_t* x, const Double_t *par) const {
+ //gaussian with a background parametrisation
+ //cout<<par[0]<<" "<<par[1]<<" "<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<<" "<<x[0]<< endl;
+ Double_t returnvalue=0.0;
+ Double_t factor=1.0/(TMath::Sqrt(2.0*TMath::Pi()));
+ 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]);
+ else returnvalue+=GausPlusTail(x[0], par[0], par[1],par[2], par[6], 1.2);
+ returnvalue+=par[3]*TMath::Exp((par[5]-x[0])*par[4]);
+ return returnvalue;
+}
+
+//______________________________________________________________________
+void AliITSsadEdxFitter::DrawFitFunction(TF1 *fun) const {
+ //code to draw the final fit function and the single gaussians used
+ //to extract the yields for the 3 species
+ TF1 *fdraw1=new TF1("fdraw1",SingleGausStep,-3.5,3.5,3);
+ TF1 *fdraw2=new TF1("fdraw2",SingleGausStep,-3.5,3.5,3);
+ TF1 *fdraw3=new TF1("fdraw3",SingleGausStep,-3.5,3.5,3);
+ fdraw1->SetParameter(0,fun->GetParameter(0));
+ fdraw1->SetParameter(1,fun->GetParameter(1));
+ fdraw1->SetParameter(2,fun->GetParameter(2));
+ fdraw2->SetParameter(0,fun->GetParameter(3));
+ fdraw2->SetParameter(1,fun->GetParameter(4));
+ fdraw2->SetParameter(2,fun->GetParameter(5));
+ fdraw3->SetParameter(0,fun->GetParameter(6));
+ fdraw3->SetParameter(1,fun->GetParameter(7));
+ fdraw3->SetParameter(2,fun->GetParameter(8));
+
+ fdraw1->SetLineColor(6);//color code
+ fdraw2->SetLineColor(2);
+ fdraw3->SetLineColor(4);
+
+ fdraw1->SetLineStyle(2);//dot lines
+ fdraw2->SetLineStyle(2);
+ fdraw3->SetLineStyle(2);
+
+ fun->SetLineWidth(3);
+ fdraw1->DrawCopy("same");
+ fdraw2->DrawCopy("same");
+ fdraw3->DrawCopy("same");
+ fun->DrawCopy("same");
+
+ TLatex *ltx[3];
+ for(Int_t j=0;j<3;j++){
+ ltx[0]=new TLatex(0.13,0.25,"pions");
+ ltx[1]=new TLatex(0.13,0.20,"kaons");
+ ltx[2]=new TLatex(0.13,0.15,"protons");
+ ltx[0]->SetTextColor(6);
+ ltx[1]->SetTextColor(2);
+ ltx[2]->SetTextColor(4);
+ ltx[j]->SetNDC();
+ ltx[j]->Draw();
+ }
+ return;
+}
+
+//________________________________________________________
+void AliITSsadEdxFitter::DoFit(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, TGraph *gres){
+ // 3-gaussian fit to log(dedx)-log(dedxBB) histogram
+ // pt bin from 0 to 20, code={211,321,2212}
+ // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed
+ // final step: refit all using the parameters and tollerance limits (+-20%)
+
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
+ TF1 *fstep1, *fstep2, *fstep3, *fstepTot;
+ Double_t initialParametersStepTot[9];
+
+ Int_t code=TMath::Abs(signedcode);
+ //************************ drawing and label *******
+ 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};
+ Double_t pt=(xbins[bin+1]+xbins[bin])/2;
+ h->SetTitle(Form("p_{t}=[%1.2f,%1.2f], code=%d",xbins[bin],xbins[bin+1],signedcode));
+ h->GetXaxis()->SetTitle("[ln dE/dx]_{meas} - [ln dE/dx(i)]_{calc}");
+ h->GetYaxis()->SetTitle("counts");
+ h->Draw("e");
+ h->SetFillColor(11);
+ Int_t xmax=-1,ymax=-1,zmax=-1;
+ if(!IsGoodBin(bin,code)) return;
+ h->GetMaximumBin(xmax,ymax,zmax);
+ TString modfit = "M0R+";
+
+ printf("\n$$$$$$$$$$$$$$$$$$$$$$$ BIN %d - hypothesis %d $$$$$$$$$$$$$$$$$$$$$$$$$$$\n",bin,code);
+ Double_t ampl = h->GetMaximum()/(h->GetRMS()*s2pi);
+ Double_t mean = h->GetBinLowEdge(xmax); //expected mean values
+ Float_t expKaonMean=0., expProtonMean=0.;
+ Int_t calcmean=CalcMean(code,pt,mean,expKaonMean,expProtonMean);
+ if(calcmean<0) return;
+ printf("mean -> %f %f %f\n",mean,expKaonMean,expProtonMean);
+ printf("integration range -> (%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]);
+
+ Float_t expPionSig = CalcSigma(211,pt,mc); //expected sigma values
+ Float_t expKaonSig = CalcSigma(321,pt,mc);
+ Float_t expProtonSig = CalcSigma(2212,pt,mc);
+ printf("sigma -> %f %f %f\n",expPionSig,expKaonSig,expProtonSig);
+ printf("sigma range -> (%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]);
+
+ printf("___________________________________________________________________\n");
+ printf("First Step: pions\n");
+ fstep1 = new TF1("step1",SingleGausStep,-3.5,3.5,3);
+ fstep1->SetParameter(0,ampl);//initial
+ fstep1->SetParameter(1,mean);
+ fstep1->SetParLimits(0,0,ampl*1.2);
+ fstep1->SetParLimits(1,mean+fRangeStep1[0],mean+fRangeStep1[1]);
+ fstep1->SetParLimits(2,expPionSig*fLimitsOnSigmaPion[0],expPionSig*fLimitsOnSigmaPion[1]);
+ //fstep1->FixParameter(2,expPionSig);
+ if(expPionSig>0) h->Fit(fstep1,modfit.Data(),"",mean+fRangeStep1[0],mean+fRangeStep1[1]);//first fit
+ else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);
+
+ printf("\n___________________________________________________________________\n");
+ printf("Second Step: kaons\n");
+ fstep2 = new TF1("fstep2",DoubleGausStep,-3.5,3.5,6);
+ fstep2->FixParameter(0,fstep1->GetParameter(0));//fixed
+ fstep2->FixParameter(1,fstep1->GetParameter(1));
+ fstep2->FixParameter(2,fstep1->GetParameter(2));
+ fstep2->SetParameter(3,fstep1->GetParameter(0)/8.);//initial
+ fstep2->SetParameter(4,expKaonMean);
+ fstep2->SetParLimits(3,fstep1->GetParameter(0)/100.,fstep1->GetParameter(0));//limits
+ fstep2->SetParLimits(4,expKaonMean+fRangeStep2[0],expKaonMean+fRangeStep2[1]);
+ fstep2->SetParLimits(5,expKaonSig*fLimitsOnSigmaKaon[0],expKaonSig*fLimitsOnSigmaKaon[1]);
+ //fstep2->FixParameter(5,expKaonSig);
+
+ if(expKaonSig>0) h->Fit(fstep2,modfit.Data(),"",expKaonMean+fRangeStep2[0],expKaonMean+fRangeStep2[1]);
+ else for(Int_t npar=3;npar<6;npar++) fstep2->FixParameter(npar,0.00001);
+
+
+ TLine *l[3];
+ l[0] = new TLine(expKaonMean,0,expKaonMean,10000);
+ l[1] = new TLine(expKaonMean+fRangeStep2[0],0,expKaonMean+fRangeStep2[0],10000);
+ l[2] = new TLine(expKaonMean+fRangeStep2[1],0,expKaonMean+fRangeStep2[1],10000);
+ for(Int_t dp=0;dp<3;dp++) {
+ l[dp]->Draw("same");
+ l[dp]->SetLineColor(4);
+ l[dp]->SetLineWidth(4);
+ }
+
+
+ printf("\n____________________________________________________________________\n");
+ printf("Third Step: protons \n");
+ fstep3= new TF1("fstep3",FinalGausStep,-3.5,3.5,9);
+ fstep3->FixParameter(0,fstep1->GetParameter(0));//fixed
+ fstep3->FixParameter(1,fstep1->GetParameter(1));
+ fstep3->FixParameter(2,fstep1->GetParameter(2));
+ fstep3->FixParameter(3,fstep2->GetParameter(3));
+ fstep3->FixParameter(4,fstep2->GetParameter(4));
+ fstep3->FixParameter(5,fstep2->GetParameter(5));
+ fstep3->SetParameter(6,fstep2->GetParameter(0)/16);//initial
+ fstep3->SetParameter(7,expProtonMean);
+ fstep3->SetParLimits(6,fstep2->GetParameter(0)/300,fstep2->GetParameter(0));//limits
+ fstep3->SetParLimits(7,expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);
+ fstep3->SetParLimits(8,expProtonSig*fLimitsOnSigmaProton[0],expProtonSig*fLimitsOnSigmaProton[1]);
+ //fstep3->FixParameter(8,expProtonSig);
+ if(expProtonSig>0) h->Fit(fstep3,modfit.Data(),"",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);
+ else for(Int_t npar=6;npar<9;npar++) fstep3->FixParameter(npar,0.00001);
+
+ printf("\n_____________________________________________________________________\n");
+ printf("Final Step: refit all \n");
+ fstepTot = new TF1("funztot",FinalGausStep,-3.5,3.5,9);
+ fstepTot->SetLineColor(1);
+ initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian
+ initialParametersStepTot[1] = fstep1->GetParameter(1);
+ initialParametersStepTot[2] = fstep1->GetParameter(2);
+
+ initialParametersStepTot[3] = fstep2->GetParameter(3);//second gaussian
+ initialParametersStepTot[4] = fstep2->GetParameter(4);
+ initialParametersStepTot[5] = fstep2->GetParameter(5);
+
+ initialParametersStepTot[6] = fstep3->GetParameter(6);//third gaussian
+ initialParametersStepTot[7] = fstep3->GetParameter(7);
+ initialParametersStepTot[8] = fstep3->GetParameter(8);
+
+ fstepTot->SetParameters(initialParametersStepTot);//initial parameter
+
+ fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1);//tolerance limit
+ fstepTot->SetParLimits(1,initialParametersStepTot[1]*0.9,initialParametersStepTot[1]*1.1);
+ fstepTot->SetParLimits(2,initialParametersStepTot[2]*0.9,initialParametersStepTot[2]*1.1);
+ //fstepTot->FixParameter(2,initialParametersStepTot[2]);
+ fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1);
+ fstepTot->SetParLimits(4,initialParametersStepTot[4]*0.9,initialParametersStepTot[4]*1.1);
+ fstepTot->SetParLimits(5,initialParametersStepTot[5]*0.9,initialParametersStepTot[5]*1.1);
+ //fstepTot->FixParameter(5,initialParametersStepTot[5]);
+ fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1);
+ fstepTot->SetParLimits(7,initialParametersStepTot[7]*0.9,initialParametersStepTot[7]*1.1);
+ fstepTot->SetParLimits(8,initialParametersStepTot[8]*0.9,initialParametersStepTot[8]*1.1);
+ //fstepTot->FixParameter(8,initialParametersStepTot[8]);
+
+ h->Fit(fstepTot,modfit.Data(),"",fRangeStep4[0],fRangeStep4[1]);//refit all
+
+
+ //************************************* storing parameter to calculate the yields *******
+ Int_t chpa=0;
+ if(code==321) chpa=3;
+ if(code==2212) chpa=6;
+ for(Int_t j=0;j<3;j++) {
+ fFitpar[j] = fstepTot->GetParameter(j+chpa);
+ fFitparErr[j] = fstepTot->GetParError(j+chpa);
+ }
+
+ DrawFitFunction(fstepTot);
+ CalcResidual(h,fstepTot,gres);
+ return;
+}
+
+//________________________________________________________
+void AliITSsadEdxFitter::DoFitProton(TH1F *h, Int_t bin, Int_t signedcode, Bool_t mc, TGraph *gres){
+ // 3-gaussian fit to log(dedx)-log(dedxBB) histogram
+ // pt bin from 0 to 20, code={211,321,2212}
+ // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed
+ // final step: refit all using the parameters and tollerance limits (+-20%)
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
+ TF1 *fstep1, *fstep2, *fstep3, *fstepTot;
+ Double_t initialParametersStepTot[9];
+
+ Int_t code=TMath::Abs(signedcode);
+ //************************ drawing and label *******
+ 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};
+ Double_t pt=(xbins[bin+1]+xbins[bin])/2;
+ h->SetTitle(Form("p_{t}=[%1.2f,%1.2f], code=%d",xbins[bin],xbins[bin+1],signedcode));
+ h->GetXaxis()->SetTitle("[ln dE/dx]_{meas} - [ln dE/dx(i)]_{calc}");
+ h->GetYaxis()->SetTitle("counts");
+ h->Draw("e");
+ h->SetFillColor(11);
+ Int_t xmax=-1,ymax=-1,zmax=-1;
+ h->GetMaximumBin(xmax,ymax,zmax);
+ if(!IsGoodBin(bin,code)) return;
+
+ printf("\n$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ BIN %d - hypothesis %d $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n",bin,code);
+ Double_t ampl = h->GetMaximum()/(h->GetRMS()*s2pi);
+ Double_t mean = h->GetBinCenter(xmax); //expected mean values
+ Float_t expKaonMean=0., expProtonMean=0.;
+ Int_t calcmean=CalcMean(code,pt,mean,expKaonMean,expProtonMean);
+ if(calcmean<0) return;
+ printf("mean -> %f %f %f\n",mean,expKaonMean,expProtonMean);
+ printf("integration range -> (%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]);
+
+ Float_t expPionSig = CalcSigma(211,pt,mc); //expected sigma values
+ Float_t expKaonSig = CalcSigma(321,pt,mc);
+ Float_t expProtonSig = CalcSigma(2212,pt,mc);
+ printf("sigma -> %f %f %f\n",expPionSig,expKaonSig,expProtonSig);
+ printf("sigma range -> (%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]);
+
+ printf("___________________________________________________________________\n");
+ printf("First Step: pions\n\n");
+ fstep1 = new TF1("step1",SingleGausStep,-3.5,3.5,3);
+ fstep1->SetParameter(0,ampl);//initial
+ fstep1->SetParameter(1,mean);
+ fstep1->SetParLimits(0,0,ampl*1.2);
+ fstep1->SetParLimits(1,mean+fRangeStep1[0],mean+fRangeStep1[1]);
+ fstep1->SetParLimits(2,expPionSig*fLimitsOnSigmaPion[0],expPionSig*fLimitsOnSigmaPion[1]);
+ //fstep1->FixParameter(2,expPionSig);
+ if(expPionSig>0) h->Fit(fstep1,"M0R+","",mean+fRangeStep1[0],mean+fRangeStep1[1]);//first fit
+ else for(Int_t npar=0;npar<3;npar++) fstep1->FixParameter(npar,0.00001);
+
+ printf("\n___________________________________________________________________\n");
+ printf("Second Step: protons\n\n");
+ fstep2 = new TF1("step2",SingleGausStep,-3.5,3.5,3);
+ fstep2->SetParameter(0,fstep1->GetParameter(0)/8.);//initial
+ fstep2->SetParameter(1,expProtonMean);
+ fstep2->SetParLimits(0,fstep1->GetParameter(0)/100.,fstep1->GetParameter(0));//limits
+ fstep2->SetParLimits(1,-3.5,3.5);
+ fstep2->SetParLimits(2,expProtonSig*fLimitsOnSigmaProton[0],expProtonSig*fLimitsOnSigmaProton[1]);
+ //fstep2->FixParameter(2,expProtonSig);
+ if(expProtonSig>0) h->Fit(fstep2,"M0R+","",expProtonMean+fRangeStep3[0],expProtonMean+fRangeStep3[1]);
+ else for(Int_t npar=0;npar<3;npar++) fstep2->FixParameter(npar,0.00001);
+
+ printf("\n____________________________________________________________________\n");
+ printf("Third Step: kaons \n\n");
+ fstep3= new TF1("fstep3",FinalGausStep,-3.5,3.5,9);
+ fstep3->FixParameter(0,fstep1->GetParameter(0));//fixed
+ fstep3->FixParameter(1,fstep1->GetParameter(1));
+ fstep3->FixParameter(2,fstep1->GetParameter(2));
+ fstep3->SetParameter(3,fstep1->GetParameter(0)/10);//initial
+ fstep3->SetParameter(4,expKaonMean);
+ fstep3->FixParameter(6,fstep2->GetParameter(0));
+ fstep3->FixParameter(7,fstep2->GetParameter(1));
+ fstep3->FixParameter(8,fstep2->GetParameter(2));
+ fstep3->SetParLimits(3,fstep1->GetParameter(0)/100.,fstep1->GetParameter(0));//limits
+ fstep3->SetParLimits(4,-3.5,3.5);
+ fstep3->SetParLimits(5,expKaonSig*fLimitsOnSigmaKaon[0],expKaonSig*fLimitsOnSigmaKaon[1]);
+ //fstep3->FixParameter(5,expKaonSig);
+
+ TLine *l[3];
+ l[0] = new TLine(expProtonMean,0,expProtonMean,10000);
+ l[1] = new TLine(expProtonMean+fRangeStep3[0],0,expProtonMean+fRangeStep3[0],10000);
+ l[2] = new TLine(expProtonMean+fRangeStep3[1],0,expProtonMean+fRangeStep3[1],10000);
+ for(Int_t dp=0;dp<3;dp++) {
+ l[dp]->Draw("same");
+ l[dp]->SetLineColor(2);
+ l[dp]->SetLineWidth(4);
+ }
+
+ if(expKaonSig>0) h->Fit(fstep3,"M0R+","",expKaonMean+fRangeStep2[0],expKaonMean+fRangeStep2[1]);
+ else for(Int_t npar=3;npar<6;npar++) fstep3->FixParameter(npar,0.00001);
+
+ printf("\n_____________________________________________________________________\n");
+ printf("Final Step: refit all \n\n");
+ fstepTot = new TF1("funztot",FinalGausStep,-3.5,3.5,9);
+ fstepTot->SetLineColor(1);
+ initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian
+ initialParametersStepTot[1] = fstep1->GetParameter(1);
+ initialParametersStepTot[2] = fstep1->GetParameter(2);
+
+ initialParametersStepTot[3] = fstep3->GetParameter(3);//second gaussian
+ initialParametersStepTot[4] = fstep3->GetParameter(4);
+ initialParametersStepTot[5] = fstep3->GetParameter(5);
+
+ initialParametersStepTot[6] = fstep2->GetParameter(0);//third gaussian
+ initialParametersStepTot[7] = fstep2->GetParameter(1);
+ initialParametersStepTot[8] = fstep2->GetParameter(2);
+
+ fstepTot->SetParameters(initialParametersStepTot);//initial parameter
+
+ fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1);//tolerance limit
+ fstepTot->SetParLimits(1,initialParametersStepTot[1]*0.9,initialParametersStepTot[1]*1.1);
+ fstepTot->SetParLimits(2,initialParametersStepTot[2]*0.9,initialParametersStepTot[2]*1.1);
+ //fstepTot->FixParameter(2,initialParametersStepTot[2]);
+ fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1);
+ fstepTot->SetParLimits(4,initialParametersStepTot[4]*0.9,initialParametersStepTot[4]*1.1);
+ fstepTot->SetParLimits(5,initialParametersStepTot[5]*0.9,initialParametersStepTot[5]*1.1);
+ //fstepTot->FixParameter(5,initialParametersStepTot[5]);
+ fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1);
+ fstepTot->SetParLimits(7,initialParametersStepTot[7]*0.9,initialParametersStepTot[7]*1.1);
+ fstepTot->SetParLimits(8,initialParametersStepTot[8]*0.9,initialParametersStepTot[8]*1.1);
+ //fstepTot->FixParameter(8,initialParametersStepTot[8]);
+
+ h->Fit(fstepTot,"M0R+","",fRangeStep4[0],fRangeStep4[1]);//refit all
+
+
+ //************************************* storing parameter to calculate the yields *******
+ Int_t chpa=0;
+ if(code==321) chpa=3;
+ if(code==2212) chpa=6;
+ for(Int_t j=0;j<3;j++) {
+ fFitpar[j] = fstepTot->GetParameter(j+chpa);
+ fFitparErr[j] = fstepTot->GetParError(j+chpa);
+ }
+
+ DrawFitFunction(fstepTot);
+ CalcResidual(h,fstepTot,gres);
+ return;
+}
+
+//________________________________________________________
+void AliITSsadEdxFitter::DoFitTail(TH1F *h, Int_t bin, Int_t code){
+ // 3-gaussian fit to log(dedx)-log(dedxBB) histogram
+ // pt bin from 0 to 20, code={211,321,2212}
+ // first step: all free, second step: pion gaussian fixed, third step: kaon gaussian fixed
+ // final step: refit all using the parameters and tollerance limits (+-20%)
+ // WARNING: exponential tail added in the right of the Gaussian shape
+ Double_t s2pi=TMath::Sqrt(2*TMath::Pi());
+ TF1 *fstep1, *fstep2, *fstep3, *fstepTot;
+ Double_t initialParametersStepTot[15];
+
+
+ //************************ drawing and label *******
+ 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};
+ h->SetTitle(Form("p_{t}=[%1.2f,%1.2f], code=%d",xbins[bin],xbins[bin+1],code));
+ h->GetXaxis()->SetTitle("[ln dE/dx]_{meas} - [ln dE/dx(i)]_{calc}");
+ h->GetYaxis()->SetTitle("counts");
+ h->Draw("");
+ h->SetFillColor(11);
+ Int_t xmax=-1,ymax=-1,zmax=-1;
+ Int_t hmax=h->GetMaximumBin(xmax,ymax,zmax);
+ hmax++;
+ if(!IsGoodBin(bin,code)) return;
+ Double_t ampl = h->GetMaximum()/(h->GetRMS()*s2pi);
+ Double_t mean = h->GetBinLowEdge(xmax);
+
+ printf("\n___________________________________________________________________\n First Step: pions\n\n");
+ fstep1 = new TF1("step1",SingleGausTail,-3.5,3.5,5);
+ fstep1->SetParameter(0,ampl);//initial
+ fstep1->SetParameter(1,mean);
+ fstep1->SetParameter(3,1.2);
+ fstep1->SetParameter(4,10.);
+
+ fstep1->SetParLimits(0,0,ampl*1.2);
+ fstep1->SetParLimits(1,-3.5,3.5);
+ fstep1->SetParLimits(2,0.1,0.25);
+ fstep1->SetParLimits(4,5.,20.);
+ if(bin<8) fstep1->SetParLimits(4,13.,25.);
+
+ h->Fit(fstep1,"M0R+","",mean-0.45,mean+0.45);//first fit
+
+ printf("\n___________________________________________________________________\n Second Step: kaons\n\n");
+ fstep2 = new TF1("fstep2",DoubleGausTail,-3.5,3.5,10);
+ fstep2->FixParameter(0,fstep1->GetParameter(0));//fixed
+ fstep2->FixParameter(1,fstep1->GetParameter(1));
+ fstep2->FixParameter(2,fstep1->GetParameter(2));
+ fstep2->FixParameter(3,fstep1->GetParameter(3));
+ fstep2->FixParameter(4,fstep1->GetParameter(4));
+
+ fstep2->SetParameter(5,fstep1->GetParameter(0)/8);//initial
+ //fstep2->SetParameter(6,CalcP(code,322,bin));
+ fstep2->SetParameter(7,0.145);
+ fstep2->FixParameter(8,1.2);
+ fstep2->SetParameter(9,13.);
+
+ fstep2->SetParLimits(5,fstep1->GetParameter(0)/100,fstep1->GetParameter(0));//limits
+ fstep2->SetParLimits(6,-3.5,3.5);
+ fstep2->SetParLimits(7,0.12,0.2);
+ fstep2->SetParLimits(9,9.,20.);
+ if(bin<9) fstep2->SetParLimits(9,13.,25.);
+
+ //h->Fit(fstep2,"M0R+","",CalcP(code,321,bin)-0.3,CalcP(code,321,bin)+0.3);//second fit
+ if(bin<6 || bin>12) for(Int_t npar=5;npar<10;npar++) fstep2->FixParameter(npar,-0.0000000001);
+
+ printf("\n____________________________________________________________________\n Third Step: protons \n\n");
+ fstep3= new TF1("fstep3",FinalGausTail,-3.5,3.5,15);
+ fstep3->FixParameter(0,fstep1->GetParameter(0));//fixed
+ fstep3->FixParameter(1,fstep1->GetParameter(1));
+ fstep3->FixParameter(2,fstep1->GetParameter(2));
+ fstep3->FixParameter(3,fstep1->GetParameter(3));
+ fstep3->FixParameter(4,fstep1->GetParameter(4));
+ fstep3->FixParameter(5,fstep2->GetParameter(5));
+ fstep3->FixParameter(6,fstep2->GetParameter(6));
+ fstep3->FixParameter(7,fstep2->GetParameter(7));
+ fstep3->FixParameter(8,fstep2->GetParameter(8));
+ fstep3->FixParameter(9,fstep2->GetParameter(9));
+
+ fstep3->SetParameter(10,fstep2->GetParameter(0)/8);//initial
+ //fstep3->SetParameter(11,CalcP(code,2212,bin));
+ fstep3->SetParameter(12,0.145);
+ fstep3->FixParameter(13,1.2);
+ fstep3->SetParameter(14,10.);
+
+ fstep3->SetParLimits(10,fstep2->GetParameter(0)/100,fstep2->GetParameter(0));//limits
+ fstep3->SetParLimits(11,-3.5,3.5);
+ fstep3->SetParLimits(12,0.12,0.2);
+ fstep3->SetParLimits(14,11.,25.);
+
+ //h->Fit(fstep3,"M0R+","",CalcP(code,2212,bin)-0.3,CalcP(code,2212,bin)+0.3);//third fit
+
+ printf("\n_____________________________________________________________________\n Final Step: refit all \n\n");
+ fstepTot = new TF1("funztot",FinalGausTail,-3.5,3.5,15);
+ fstepTot->SetLineColor(1);
+ initialParametersStepTot[0] = fstep1->GetParameter(0);//first gaussian
+ initialParametersStepTot[1] = fstep1->GetParameter(1);
+ initialParametersStepTot[2] = fstep1->GetParameter(2);
+ initialParametersStepTot[3] = fstep1->GetParameter(3);
+ initialParametersStepTot[4] = fstep1->GetParameter(4);
+
+ initialParametersStepTot[5] = fstep2->GetParameter(5);//second gaussian
+ initialParametersStepTot[6] = fstep2->GetParameter(6);
+ initialParametersStepTot[7] = fstep2->GetParameter(7);
+ initialParametersStepTot[8] = fstep2->GetParameter(8);
+ initialParametersStepTot[9] = fstep2->GetParameter(9);
+
+ initialParametersStepTot[10] = fstep3->GetParameter(10);//third gaussian
+ initialParametersStepTot[11] = fstep3->GetParameter(11);
+ initialParametersStepTot[12] = fstep3->GetParameter(12);
+ initialParametersStepTot[13] = fstep3->GetParameter(13);
+ initialParametersStepTot[14] = fstep3->GetParameter(14);
+
+ fstepTot->SetParameters(initialParametersStepTot);//initial parameter
+
+
+ fstepTot->SetParLimits(0,initialParametersStepTot[0]*0.9,initialParametersStepTot[0]*1.1);//tollerance limit
+ fstepTot->SetParLimits(1,initialParametersStepTot[1]*0.9,initialParametersStepTot[1]*1.1);
+ fstepTot->SetParLimits(2,initialParametersStepTot[2]*0.9,initialParametersStepTot[2]*1.1);
+ fstepTot->SetParLimits(3,initialParametersStepTot[3]*0.9,initialParametersStepTot[3]*1.1);
+ fstepTot->SetParLimits(4,initialParametersStepTot[4]*0.9,initialParametersStepTot[4]*1.1);
+ fstepTot->SetParLimits(5,initialParametersStepTot[5]*0.9,initialParametersStepTot[5]*1.1);
+ fstepTot->SetParLimits(6,initialParametersStepTot[6]*0.9,initialParametersStepTot[6]*1.1);
+ fstepTot->SetParLimits(7,initialParametersStepTot[7]*0.9,initialParametersStepTot[7]*1.1);
+ fstepTot->SetParLimits(8,initialParametersStepTot[8]*0.9,initialParametersStepTot[8]*1.1);
+ fstepTot->SetParLimits(9,initialParametersStepTot[9]*0.9,initialParametersStepTot[9]*1.1);
+ fstepTot->SetParLimits(10,initialParametersStepTot[10]*0.9,initialParametersStepTot[10]*1.1);
+ fstepTot->SetParLimits(11,initialParametersStepTot[11]*0.9,initialParametersStepTot[11]*1.1);
+ fstepTot->SetParLimits(12,initialParametersStepTot[12]*0.9,initialParametersStepTot[12]*1.1);
+ fstepTot->SetParLimits(13,initialParametersStepTot[13]*0.9,initialParametersStepTot[13]*1.1);
+ fstepTot->SetParLimits(14,initialParametersStepTot[14]*0.9,initialParametersStepTot[14]*1.1);
+
+ if(bin<9) for(Int_t npar=10;npar<15;npar++) fstepTot->FixParameter(npar,-0.00000000001);
+
+ h->Fit(fstepTot,"M0R+","",-3.5,3.5); //refit all
+
+
+ //************************************* storing parameter to calculate the yields *******
+ Int_t chpa=0;
+ if(code==321) chpa=5;
+ if(code==2212) chpa=10;
+ for(Int_t j=0;j<5;j++) {
+ fFitpar[j] = fstepTot->GetParameter(j+chpa);
+ fFitparErr[j] = fstepTot->GetParError(j+chpa);
+ }
+
+ DrawFitFunction(fstepTot);
+ return;
+}
+
+//________________________________________________________
+void AliITSsadEdxFitter::FillHisto(TH1F *hsps, Int_t bin, Float_t binsize, Int_t code){
+ // fill the spectra histo calculating the yield
+ // first bin has to be 1
+ Double_t yield = 0;
+ Double_t err = 0;
+ Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin);
+
+ if(IsGoodBin(bin-1,code)) {
+ yield = fFitpar[0] / ptbin / binsize;
+ err = fFitparErr[0] / ptbin / binsize;
+ }
+
+ hsps->SetBinContent(bin,yield);
+ hsps->SetBinError(bin,err);
+ return;
+}
+
+//________________________________________________________
+void AliITSsadEdxFitter::FillHistoMC(TH1F *hsps, Int_t bin, Int_t code, TH1F *h){
+ // fill the spectra histo calculating the yield (using the MC truth)
+ // first bin has to be 1
+ Double_t yield = 0;
+ Double_t erryield=0;
+ Double_t ptbin = hsps->GetBinLowEdge(bin+1) - hsps->GetBinLowEdge(bin);
+
+ if(IsGoodBin(bin-1,code)){
+ yield = h->GetEntries() / ptbin;
+ erryield=TMath::Sqrt(h->GetEntries()) / ptbin;
+ }
+ hsps->SetBinContent(bin,yield);
+ hsps->SetBinError(bin,erryield);
+ return;
+}
+
+//________________________________________________________
+void AliITSsadEdxFitter::GetFitPar(Double_t *fitpar, Double_t *fitparerr) const {
+ // getter of the fit parameters and the relative errors
+ for(Int_t i=0;i<5;i++) {
+ fitpar[i] = fFitpar[i];
+ fitparerr[i] = fFitparErr[i];
+ }
+ return;
+}
+
+
+//________________________________________________________
+void AliITSsadEdxFitter::PrintAll() const{
+ printf("Range 1 = %f %f\n",fRangeStep1[0],fRangeStep1[1]);
+ printf("Range 2 = %f %f\n",fRangeStep2[0],fRangeStep2[1]);
+ printf("Range 3 = %f %f\n",fRangeStep3[0],fRangeStep3[1]);
+ printf("Range F = %f %f\n",fRangeStep4[0],fRangeStep4[1]);
+ printf(" Sigma1 = %f %f\n",fLimitsOnSigmaPion[0],fLimitsOnSigmaPion[1]);
+ printf(" Sigma2 = %f %f\n",fLimitsOnSigmaKaon[0],fLimitsOnSigmaKaon[1]);
+ printf(" Sigma3 = %f %f\n",fLimitsOnSigmaProton[0],fLimitsOnSigmaProton[1]);
+}
--- /dev/null
+#ifndef ALIITSSADEDXFITTER_H
+#define ALIITSSADEDXFITTER_H
+/* Copyright(c) 2007-2011, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+/* $Id: $ */
+
+////////////////////////////////////////////////////
+//class to perform different gaussian fits to the //
+//dEdx distribution, using different approach, //
+//for the ITS stand-alone track spectra analysis //
+//E. Biolcati, F. Prino //
+////////////////////////////////////////////////////
+
+#include <TObject.h>
+class TGraph;
+
+class AliITSsadEdxFitter : public TObject {
+
+ public:
+ AliITSsadEdxFitter();
+ virtual ~AliITSsadEdxFitter(){};
+
+ static Double_t CalcSigma(Int_t code,Float_t x, Bool_t mc);
+ static Int_t CalcMean(Int_t code,Float_t x, Float_t mean0, Float_t &mean1, Float_t &mean2);
+
+ void GetFitPar(Double_t *fitpar, Double_t *fitparerr) const;
+ void DoFitTail(TH1F *h, Int_t bin, Int_t code);
+ void DoFit(TH1F *h, Int_t bin, Int_t code, Bool_t mc, TGraph *gres);
+ void DoFitProton(TH1F *h, Int_t bin, Int_t code, Bool_t mc, TGraph *gres);
+ void FillHisto(TH1F *hsps, Int_t bin, Float_t binsize, Int_t code);
+ void FillHistoMC(TH1F *hsps, Int_t bin, Int_t code, TH1F *h);
+ Bool_t IsGoodBin(Int_t bin,Int_t code);
+
+ void SetRangeStep1(Double_t dxlow=-0.2, Double_t dxup=0.3){
+ fRangeStep1[0]=dxlow;
+ fRangeStep1[1]=dxup;
+ }
+ void SetRangeStep2(Double_t dxlow=-0.1, Double_t dxup=0.3){
+ fRangeStep2[0]=dxlow;
+ fRangeStep2[1]=dxup;
+ }
+ void SetRangeStep3(Double_t dxlow=-0.1, Double_t dxup=2.5){
+ fRangeStep3[0]=dxlow;
+ fRangeStep3[1]=dxup;
+ }
+ void SetRangeFinalStep(Double_t dxlow=-3.5, Double_t dxup=3.5){
+ fRangeStep4[0]=dxlow;
+ fRangeStep4[1]=dxup;
+ }
+ void SetLimitsOnSigmaPion(Double_t smin=0.98, Double_t smax=1.02){
+ fLimitsOnSigmaPion[0]=smin;
+ fLimitsOnSigmaPion[1]=smax;
+ }
+ void SetLimitsOnSigmaKaon(Double_t smin=0.98, Double_t smax=1.02){
+ fLimitsOnSigmaKaon[0]=smin;
+ fLimitsOnSigmaKaon[1]=smax;
+ }
+ void SetLimitsOnSigmaProton(Double_t smin=0.98, Double_t smax=1.02){
+ fLimitsOnSigmaProton[0]=smin;
+ fLimitsOnSigmaProton[1]=smax;
+ }
+
+ void PrintAll() const;
+ void CalcResidual(TH1F *h,TF1 *fun,TGraph *gres) const;
+ Double_t GausPlusTail(const Double_t x, const Double_t mean, Double_t rms, Double_t c, Double_t slope, Double_t cut ) const;
+ Double_t GausOnBackground(const Double_t* x, const Double_t *par) const;
+ void DrawFitFunction(TF1 *fun) const;
+
+ private:
+ Double_t fFitpar[5]; // array with fit parameters
+ Double_t fFitparErr[5]; // array with fit parameter errors
+ Double_t fRangeStep1[2]; // Range for Step1 (w.r.t pion peak)
+ Double_t fRangeStep2[2]; // Range for Step2 (w.r.t kaon/proton peak)
+ Double_t fRangeStep3[2]; // Range for Step3 (w.r.t proton/kaon peak)
+ Double_t fRangeStep4[2]; // Range for Last Fit
+ Double_t fLimitsOnSigmaPion[2]; // limits on sigma pions
+ Double_t fLimitsOnSigmaKaon[2]; // limits on sigma pions
+ Double_t fLimitsOnSigmaProton[2]; // limits on sigma protons
+
+ ClassDef(AliITSsadEdxFitter,1);
+};
+
+#endif
+
--- /dev/null
+/////////////////////////////////////////////////////////
+//Macro to read the results of the 3gausfit for spectra
+//E. Biolcati, 13-feb-2010
+/////////////////////////////////////////////////////////
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include <TLatex.h>
+#include <TH1F.h>
+#include <TH2D.h>
+#include <TF1.h>
+#include <TMath.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+#include <TGraphErrors.h>
+#include <TGraphAsymmErrors.h>
+#include <TTree.h>
+#include <TCanvas.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include <TStyle.h>
+#include <TPad.h>
+#endif
+
+void Labella(Int_t i);
+void Legenda(TH1 *h2, TH1 *h3, TH1 *h4);
+
+//_______________________________________________________
+void MakeCorrectedITSsaSpectraMultiBin(Int_t multibin=0){
+
+ TString dirNameSIM, dirNameDATA;
+ TFile *fiSIM = new TFile(Form("../gridmultiplicitybins/LHC10d1_1.5sigma_7DCA_negmag/Spectra_MC_negmag_MultBin%d/SpectraReco.root",multibin));
+ TFile *fiDATA = new TFile(Form("../gridmultiplicitybins/data_1.5sigma_7DCA_negmag/Spectra_data_negmag_MultBin%d/SpectraReco.root",multibin));
+ TFile *fout = new TFile(Form("ITSsaSpectraCorr_MultiBin%d.root",multibin),"recreate");
+
+ //dca correction
+ TFile *fPanosPosP= new TFile(Form("RootFilesPanosCorr/PanosCorr_1.5sigma_7DCA_PosP_expostrange_MultBin%d.root",multibin));
+ TFile *fPanosNegP= new TFile(Form("RootFilesPanosCorr/PanosCorr_1.5sigma_7DCA_NegP_expostrange_MultBin%d.root",multibin));
+ TH1F *hPanosPosP= (TH1F*)fPanosPosP->Get("fHistSecTOTCorrDATAMC");
+ TH1F *hPanosNegP= (TH1F*)fPanosNegP->Get("fHistSecTOTCorrDATAMC");
+ for(Int_t pbin=0;pbin<=hPanosPosP->GetNbinsX();pbin++)hPanosPosP->SetBinError(pbin,0);
+ for(Int_t pbin=0;pbin<=hPanosNegP->GetNbinsX();pbin++)hPanosNegP->SetBinError(pbin,0);
+
+ //nevts
+ TH1F *hstat = (TH1F*)fiDATA->Get("fHistNEvents");
+ Double_t NEvts = hstat->GetBinContent(hstat->FindBin(1.));
+ cout<<"Event number used for the normalization "<<NEvts<<endl;
+
+ //canvas
+ gStyle->SetOptStat(0);
+ TCanvas *cs=new TCanvas("cs","cs",1200,700);
+ cs->Divide(2,1);
+ TCanvas *cmix=new TCanvas("cmix","cmix",800,900);
+ cmix->Divide(1,3,0,0);
+
+ TH1F *hPieff[2];
+ TH1F *hKeff[2];
+ TH1F *hPeff[2];
+ TH1F *hPiDATA[2];
+ TH1F *hKDATA[2];
+ TH1F *hPDATA[2];
+ TH1F *hPiDATANorm[2];
+ TH1F *hKDATANorm[2];
+ TH1F *hPDATANorm[2];
+ TH1F *hKsuPi[2];
+ TH1F *hPsuPi[2];
+ TH1F *hPsuK[2];
+
+ hPieff[0] = (TH1F*)fiSIM->Get("hCorrFacNeg0");
+ hKeff[0] = (TH1F*)fiSIM->Get("hCorrFacNeg1");
+ hPeff[0] = (TH1F*)fiSIM->Get("hCorrFacNeg2");
+ hPieff[1] = (TH1F*)fiSIM->Get("hCorrFacPos0");
+ hKeff[1] = (TH1F*)fiSIM->Get("hCorrFacPos1");
+ hPeff[1] = (TH1F*)fiSIM->Get("hCorrFacPos2");
+ TCanvas *ceff=new TCanvas("ceff","ceff",800,500);
+ ceff->Divide(2,1);
+ ceff->cd(1);
+ gPad->SetGridy();
+ hPieff[0]->Draw();
+ hKeff[0]->Draw("same");
+ hPeff[0]->Draw("same");
+ ceff->cd(2);
+ gPad->SetGridy();
+ hPieff[1]->Draw();
+ hKeff[1]->Draw("same");
+ hPeff[1]->Draw("same");
+
+ hPiDATA[0] = (TH1F*)fiDATA->Get("hSpectraNeg0");
+ hKDATA[0] = (TH1F*)fiDATA->Get("hSpectraNeg1");
+ hPDATA[0] = (TH1F*)fiDATA->Get("hSpectraNeg2");
+ hPiDATA[1] = (TH1F*)fiDATA->Get("hSpectraPos0");
+ hKDATA[1] = (TH1F*)fiDATA->Get("hSpectraPos1");
+ hPDATA[1] = (TH1F*)fiDATA->Get("hSpectraPos2");
+
+ hPiDATANorm[0] = new TH1F(*hPiDATA[0]);
+ hPiDATANorm[0]->SetName("hSpectraPiPlusN");
+ hPiDATANorm[1] = new TH1F(*hPiDATA[1]);
+ hPiDATANorm[1]->SetName("hSpectraPiMinusN");
+ hKDATANorm[0] = new TH1F(*hKDATA[0]);
+ hKDATANorm[0]->SetName("hSpectraKPlusN");
+ hKDATANorm[1] = new TH1F(*hKDATA[1]);
+ hKDATANorm[1]->SetName("hSpectraKMinusN");
+ hPDATANorm[0] = new TH1F(*hPDATA[0]);
+ hPDATANorm[0]->SetName("hSpectraPPlusN");
+ hPDATANorm[1] = new TH1F(*hPDATA[1]);
+ hPDATANorm[1]->SetName("hSpectraPMinusN");
+
+ for(Int_t i=0;i<2;i++){ //0==> negative, 1==>positive
+ //line colors
+ hPiDATA[i] -> SetLineColor(2);
+ hKDATA[i] -> SetLineColor(3);
+ hPDATA[i] -> SetLineColor(4);
+ hPiDATA[i] -> SetMarkerStyle(27);
+ hKDATA[i] -> SetMarkerStyle(23);
+ hPDATA[i] -> SetMarkerStyle(24);
+ hPiDATA[i] -> SetMarkerColor(2);
+ hKDATA[i] -> SetMarkerColor(3);
+ hPDATA[i] -> SetMarkerColor(4);
+ hPiDATA[i] -> SetMarkerStyle(23);
+ hKDATA[i] -> SetMarkerStyle(23);
+ hPDATA[i] -> SetMarkerStyle(23);
+
+ hPiDATANorm[i] -> SetLineColor(2);
+ hKDATANorm[i] -> SetLineColor(3);
+ hPDATANorm[i] -> SetLineColor(4);
+ hPiDATANorm[i] -> SetMarkerStyle(27);
+ hKDATANorm[i] -> SetMarkerStyle(23);
+ hPDATANorm[i] -> SetMarkerStyle(24);
+ hPiDATANorm[i] -> SetMarkerColor(2);
+ hKDATANorm[i] -> SetMarkerColor(3);
+ hPDATANorm[i] -> SetMarkerColor(4);
+ hPiDATANorm[i] -> SetMarkerStyle(23);
+ hKDATANorm[i] -> SetMarkerStyle(23);
+ hPDATANorm[i] -> SetMarkerStyle(23);
+
+ //division for efficiency
+ hPiDATA[i] -> Divide(hPieff[i]);
+ hKDATA[i] -> Divide(hKeff[i]);
+ hPDATA[i] -> Divide(hPeff[i]);
+ hPiDATANorm[i] -> Divide(hPieff[i]);
+ hKDATANorm[i] -> Divide(hKeff[i]);
+ hPDATANorm[i] -> Divide(hPeff[i]);
+
+ //normalization number of events
+ hPiDATANorm[i] -> Scale(1./NEvts);
+ hKDATANorm[i] -> Scale(1./NEvts);
+ hPDATANorm[i] -> Scale(1./NEvts);
+
+ //correction factor based on fit to DCA distr for P and Pibar
+ hPDATA[0]->Multiply(hPanosPosP);
+ hPDATA[1]->Multiply(hPanosNegP);
+
+ //drawing
+ cs->cd(i+1);
+ gPad->SetLogy();
+ gPad->SetGridy();
+ gPad->SetGridx();
+ Float_t minim=0.01;
+ Float_t maxim=10.;
+ hPiDATANorm[i]->Draw("e");
+ hPiDATANorm[i]->GetYaxis()->SetRangeUser(minim,maxim);
+ hPiDATANorm[i]->GetYaxis()->SetTitle("#frac{1}{N}#frac{d^{2}N}{dp_{t}dy}");
+ hPiDATANorm[i]->SetTitle("DATA from ITSsa - corrected");
+ hKDATANorm[i]->Draw("esames");
+ hPDATANorm[i]->Draw("esames");
+ Labella(i);
+ Legenda(hPiDATA[i],hKDATA[i],hPDATA[i]);
+ cs->Update();
+ }
+
+ //fluka correction for pi
+ TFile *fGeanFlukaPi= new TFile(Form("RootFilesGeantFlukaCorrection/correctionForCrossSection.211.root"));
+ TH1F *hGeantFlukaPiPos=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionParticles");
+ TH1F *hGeantFlukaPiNeg=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionAntiParticles");
+ for(Int_t binPi=0;binPi<=hPiDATA[0]->GetNbinsX();binPi++){
+ Float_t FlukaCorrPiPos=hGeantFlukaPiPos->GetBinContent(hGeantFlukaPiPos->FindBin(hPiDATA[0]->GetBinCenter(binPi)));
+ Float_t FlukaCorrPiNeg=hGeantFlukaPiNeg->GetBinContent(hGeantFlukaPiNeg->FindBin(hPiDATA[1]->GetBinCenter(binPi)));
+ //cout<<"PiPos "<<FlukaCorrPiPos<<" "<<hPiDATA[0]->GetBinCenter(binPi)<<" " <<binPi <<endl;
+ //cout<<"PiNeg "<<FlukaCorrPiNeg<<" "<<hPiDATA[0]->GetBinCenter(binPi)<<" " <<binPi <<endl;
+ hPiDATA[0]->SetBinContent(binPi,hPiDATA[0]->GetBinContent(binPi)*FlukaCorrPiPos);
+ hPiDATA[1]->SetBinContent(binPi,hPiDATA[1]->GetBinContent(binPi)*FlukaCorrPiNeg);
+ }
+ //fluka correction for pi
+ TFile *fGeanFlukaK= new TFile(Form("RootFilesGeantFlukaCorrection/correctionForCrossSection.321.root"));
+ TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
+ TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
+ for(Int_t binK=0;binK<=hKDATA[0]->GetNbinsX();binK++){
+ Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(hKDATA[0]->GetBinCenter(binK)));
+ Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(hKDATA[1]->GetBinCenter(binK)));
+ //cout<<"KPos :"<<FlukaCorrKPos<<" "<<hKDATA[0]->GetBinCenter(binK)<<" " <<binK <<endl;
+ //cout<<"KNeg :"<<FlukaCorrKNeg<<" "<<hKDATA[0]->GetBinCenter(binK)<<" " <<binK <<endl;
+ hKDATA[0]->SetBinContent(binK,hKDATA[0]->GetBinContent(binK)*FlukaCorrKPos);
+ hKDATA[1]->SetBinContent(binK,hKDATA[1]->GetBinContent(binK)*FlukaCorrKNeg);
+ }
+ //fluka correction for P
+ //ITS specific file for protons/antiprotons
+ Int_t kPos=0;
+ Int_t kNeg=1;
+ TFile* fITS = new TFile ("RootFilesGeantFlukaCorrection/correctionForCrossSectionITS_20100719.root");
+ TH2D * hCorrFlukaITS[2];
+ hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
+ hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
+ for(Int_t icharge = 0; icharge < 2; icharge++){
+ Int_t nbins = hPDATA[0]->GetNbinsX();
+ Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
+ for(Int_t ibin = 0; ibin < nbins; ibin++){
+ Float_t pt = hPDATA[0]->GetBinCenter(ibin);
+ Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
+ Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
+ if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
+ if (pt > maxPtCorrection) pt = maxPtCorrection;
+ Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
+ if(icharge==0){
+ if (correction != 0) {// If the bin is empty this is a 0
+ hPDATA[0]->SetBinContent(ibin,hPDATA[0]->GetBinContent(ibin)*correction);
+ hPDATA[0]->SetBinError(ibin,hPDATA[0]->GetBinError (ibin)*correction);
+ }else if (hPDATA[0]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
+ cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << endl;
+ cout << " Bin content: " << hPDATA[0]->GetBinContent(ibin) << endl;
+ }
+ }
+ if(icharge==1){
+ if (correction != 0) {// If the bin is empty this is a 0
+ hPDATA[1]->SetBinContent(ibin,hPDATA[1]->GetBinContent(ibin)*correction);
+ hPDATA[1]->SetBinError(ibin,hPDATA[1]->GetBinError (ibin)*correction);
+ }else if (hPDATA[1]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
+ cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, ITS, " << endl;
+ cout << " Bin content: " << hPDATA[1]->GetBinContent(ibin) << endl;
+ }
+ }
+ }
+ }
+
+
+ //mixed particle ratios
+ for(Int_t i=0;i<2;i++){
+ hKsuPi[i] = (TH1F*)hKDATA[i]->Clone("KsuPi");
+ hPsuPi[i] = (TH1F*)hPDATA[i]->Clone("PsuPi");
+ hPsuK[i] = (TH1F*)hPDATA[i]->Clone("PsuK");
+ hKsuPi[i] -> Divide(hPiDATA[i]);
+ hPsuPi[i] -> Divide(hPiDATA[i]);
+ hPsuK[i] -> Divide(hKDATA[i]);
+ }
+
+ //positive/negative ratios
+ TH1F *hPiratio = (TH1F*)hPiDATA[1]->Clone("PionsRatio");
+ TH1F *hKratio = (TH1F*)hKDATA[1]->Clone("KaonsRatio");
+ TH1F * hPratio = (TH1F*)hPDATA[1]->Clone("ProtonsRatio");
+ hPiratio -> Divide(hPiDATA[0]);
+ hKratio -> Divide(hKDATA[0]);
+ hPratio -> Divide(hPDATA[0]);
+
+
+ //drawing positive/negative ratios
+ TCanvas *cratio=new TCanvas("cratio","",980,600);
+ cratio->Divide(1,3,0,0);
+ cratio->SetBottomMargin(0.08);
+ cratio->cd(1);
+ gPad->SetGridy();
+ hPiratio->GetYaxis()->SetTitle("");
+ hPiratio->GetYaxis()->SetRangeUser(0.7,1.3);
+ hPiratio->SetTitle("");
+ hPiratio->Draw("mp");
+ TLatex *ll1=new TLatex(0.7,0.7,"#pi^{+}/#pi^{-}");
+ ll1->SetNDC();
+ ll1->SetTextSize(0.14);
+ ll1->Draw();
+ cratio->cd(2);
+ gPad->SetGridy();
+ hKratio->SetTitle("");
+ hKratio->GetYaxis()->SetRangeUser(0.7,1.3);
+ hKratio->Draw("mp");
+ TLatex *ll2=new TLatex(0.7,0.7,"K^{+}/K^{-}");
+ ll2->SetNDC();
+ ll2->SetTextSize(0.14);
+ ll2->Draw();
+ cratio->cd(3);
+ gPad->SetGridy();
+ hPratio->SetTitle("");
+ hPratio->GetYaxis()->SetRangeUser(0.7,1.3);
+ hPratio->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ hPratio->Draw("mp");
+ TLatex *ll3=new TLatex(0.7,0.7,"p/#bar{p}");
+ ll3->SetNDC();
+ ll3->SetTextSize(0.144);
+ ll3->Draw();
+
+ //drawing mixed particle ratios
+ gStyle->SetOptTitle(0);
+ cmix->cd(1);
+ hKsuPi[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ hKsuPi[0]->Draw();
+ hKsuPi[1]->Draw("same");
+ hKsuPi[0]->SetMinimum(0);
+ hKsuPi[0]->SetMarkerStyle(23);
+ hKsuPi[1]->SetMarkerStyle(24);
+ TLegend *legm1=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
+ legm1->AddEntry(hKsuPi[0],"K^{-}/#pi^{-}","p");
+ legm1->AddEntry(hKsuPi[1],"K^{+}/#pi^{+}","p");
+ legm1->SetFillColor(0);
+ legm1->SetBorderSize(0);
+ legm1->Draw();
+ hKsuPi[0]->SetMarkerColor(2);
+ hKsuPi[1]->SetMarkerColor(4);
+ hKsuPi[0]->SetLineColor(2);
+ hKsuPi[1]->SetLineColor(4);
+
+ cmix->cd(2);
+ hPsuPi[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ hPsuPi[0]->Draw();
+ hPsuPi[1]->Draw("same");
+ hPsuPi[0]->SetMinimum(0);
+ hPsuPi[0]->SetMarkerStyle(23);
+ hPsuPi[1]->SetMarkerStyle(24);
+ hPsuPi[0]->SetMarkerColor(2);
+ hPsuPi[1]->SetMarkerColor(4);
+ hPsuPi[0]->SetLineColor(2);
+ hPsuPi[1]->SetLineColor(4);
+ TLegend *legm2=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
+ legm2->AddEntry(hPsuPi[0],"#bar{p}/#pi^{-}","p");
+ legm2->AddEntry(hPsuPi[1],"p/#pi^{+}","p");
+ legm2->SetFillColor(0);
+ legm2->SetBorderSize(0);
+ legm2->Draw();
+
+ cmix->cd(3);
+ hPsuK[0]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ hPsuK[0]->Draw();
+ hPsuK[1]->Draw("same");
+ hPsuK[0]->SetMinimum(0);
+ hPsuK[0]->SetMarkerStyle(23);
+ hPsuK[1]->SetMarkerStyle(24);
+ hPsuK[0]->SetMarkerColor(2);
+ hPsuK[1]->SetMarkerColor(4);
+ hPsuK[0]->SetLineColor(2);
+ hPsuK[1]->SetLineColor(4);
+ TLegend *legm3=new TLegend(0.2,0.6,0.39,0.89,NULL,"brNDC");
+ legm3->AddEntry(hPsuPi[0],"#bar{p}/K^{-}","p");
+ legm3->AddEntry(hPsuPi[1],"p/K^{+}","p");
+ legm3->SetFillColor(0);
+ legm3->SetBorderSize(0);
+ legm3->Draw();
+
+ //save histograms in the root files
+ fout->cd();
+ hPiDATA[0]->Write();
+ hKDATA[0] ->Write();
+ hPDATA[0] ->Write();
+ hPiDATA[1]->Write();
+ hKDATA[1] ->Write();
+ hPDATA[1] ->Write();
+ hPiDATANorm[0]->Write();
+ hKDATANorm[0] ->Write();
+ hPDATANorm[0] ->Write();
+ hPiDATANorm[1]->Write();
+ hKDATANorm[1] ->Write();
+ hPDATANorm[1] ->Write();
+ //fout->Close();
+ return;
+
+}//end of the main
+
+
+//_______________________________________________________
+void Labella(Int_t i){
+ Char_t txt[50];
+ if(i==0) sprintf(txt,"negative particles");
+ else sprintf(txt,"positive particles");
+ TLatex *ltx=new TLatex(0.4,0.3,txt);
+ ltx->SetNDC();
+ ltx->SetTextColor(6);
+ ltx->SetTextFont(22);
+ ltx->Draw();
+ return;
+}
+
+//_______________________________________________________
+void Legenda(TH1 *h2, TH1 *h3, TH1 *h4){
+ TLegend *leg=new TLegend(0.51,0.11,0.84,0.25,NULL,"brNDC");
+ leg->SetFillColor(0);
+ leg->SetBorderSize(0);
+ TLegendEntry *entry2=leg->AddEntry(h2,"pions","p");
+ entry2->SetTextColor(2);
+ TLegendEntry *entry3=leg->AddEntry(h3,"kaons","p");
+ entry3->SetTextColor(3);
+ TLegendEntry *entry4=leg->AddEntry(h4,"protons","p");
+ entry4->SetTextColor(4);
+ leg->Draw("same");
+}
+
+//EOF
+
--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include <TLatex.h>
+#include <TImage.h>
+#include <TSystem.h>
+#include <TPaveText.h>
+#include <TH1F.h>
+#include <TF1.h>
+#include <TH1D.h>
+#include <TH2F.h>
+#include <TMath.h>
+#include <TNtuple.h>
+#include <TGraphErrors.h>
+#include <TList.h>
+#include <TLegend.h>
+#include <TLegendEntry.h>
+#include <TCanvas.h>
+#include <TFile.h>
+#include <TStyle.h>
+#include <Rtypes.h>
+#endif
+
+#include "AliITSspectra.h"
+
+Double_t BetheBloch(Double_t *mom, Double_t *mass);
+void Logo();
+//
+
+//______________________________________________________________________
+void MakeRawITSsaSpectraMultiBin(Bool_t optMC=kTRUE, Int_t multibin=0){
+
+ //AliITSspectraObject
+ AliITSspectra *ITSsa=new AliITSspectra();
+
+ TString filename, dirName, ps0, ps1, ps2;
+ Char_t listname[50];
+ if(optMC) dirName=(Form("../gridmultiplicitybins/LHC10d1_1.5sigma_7DCA_negmag/Spectra_MC_negmag_MultBin%d",multibin));
+ else dirName=(Form("../gridmultiplicitybins/data_1.5sigma_7DCA_negmag/Spectra_data_negmag_MultBin%d",multibin));
+ switch(multibin){
+ case 0:
+ sprintf(listname,"clistITSsaMult0to9999");
+ break;
+ case 1:
+ sprintf(listname,"clistITSsaMult0to5");
+ break;
+ case 2:
+ sprintf(listname,"clistITSsaMult6to9");
+ break;
+ case 3:
+ sprintf(listname,"clistITSsaMult10to14");
+ break;
+ case 4:
+ sprintf(listname,"clistITSsaMult15to22");
+ break;
+ case 5:
+ sprintf(listname,"clistITSsaMult23to9999");
+ break;
+ }
+ if(optMC){
+ filename=Form("%s/AnalysisResults.root",dirName.Data());
+ ps0 = Form("%s/outSIM.ps[",dirName.Data());
+ ps1 = Form("%s/outSIM.ps",dirName.Data());
+ ps2 = Form("%s/outSIM.ps]",dirName.Data());
+ }
+ else{
+ filename=Form("%s/AnalysisResults.root",dirName.Data());
+ ps0 = Form("%s/outDATA.ps[",dirName.Data());
+ ps1 = Form("%s/outDATA.ps",dirName.Data());
+ ps2 = Form("%s/outDATA.ps]",dirName.Data());
+ }
+ TString openfilename="./";
+ openfilename+=filename;
+ cout<<openfilename<<endl;
+ TFile *fi=new TFile(openfilename.Data());
+ if(!fi){
+ cout<<"TFile loading failed"<<endl;
+ return;
+ }
+ TDirectoryFile *di=(TDirectoryFile*) fi->Get("PWG2SpectraITSsa");
+ if(!di){
+ cout<<"TDirectory loading failed!"<<endl;
+ return;
+ }
+ TList *li=(TList*)di->Get(listname);
+ if(!li){
+ cout<<"TList loading failed"<<endl;
+ return;
+ }
+ cout<<"File loaded"<<endl;
+
+ TCanvas *cdummy=new TCanvas("dummy","IntegralMethod",1000,800);
+ cdummy->Print(ps0.Data());
+
+ //binning
+ const Int_t nbins = 22;
+ Double_t xbins[nbins+1]={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};
+
+ //histograms
+ TH1F *fHistMCPosPi[nbins];
+ TH1F *fHistMCPosK[nbins];
+ TH1F *fHistMCPosP[nbins];
+ TH1F *fHistMCNegPi[nbins];
+ TH1F *fHistMCNegK[nbins];
+ TH1F *fHistMCNegP[nbins];
+ TH1F *fHistPosPi[nbins];
+ TH1F *fHistPosK[nbins];
+ TH1F *fHistPosP[nbins];
+ TH1F *fHistNegPi[nbins];
+ TH1F *fHistNegK[nbins];
+ TH1F *fHistNegP[nbins];
+ for(Int_t m=0;m<nbins;m++){
+ fHistMCNegPi[m] = (TH1F*)li->FindObject(Form("fHistMCNegPi%d",m));
+ fHistMCNegK[m] = (TH1F*)li->FindObject(Form("fHistMCNegK%d",m));
+ fHistMCNegP[m] = (TH1F*)li->FindObject(Form("fHistMCNegP%d",m));
+ fHistMCPosPi[m] = (TH1F*)li->FindObject(Form("fHistMCPosPi%d",m));
+ fHistMCPosK[m] = (TH1F*)li->FindObject(Form("fHistMCPosK%d",m));
+ fHistMCPosP[m] = (TH1F*)li->FindObject(Form("fHistMCPosP%d",m));
+ fHistPosPi[m] = (TH1F*)li->FindObject(Form("fHistPosPi%d",m));
+ fHistPosK[m] = (TH1F*)li->FindObject(Form("fHistPosK%d",m));
+ fHistPosP[m] = (TH1F*)li->FindObject(Form("fHistPosP%d",m));
+ fHistNegPi[m] = (TH1F*)li->FindObject(Form("fHistNegPi%d",m));
+ fHistNegK[m] = (TH1F*)li->FindObject(Form("fHistNegK%d",m));
+ fHistNegP[m] = (TH1F*)li->FindObject(Form("fHistNegP%d",m));
+ }
+
+ TH1F *fHistNEvents = (TH1F*)li->FindObject("fHistNEvents");
+ TH2F *fHistDEDX = (TH2F*)li->FindObject("fHistDEDX");
+ TH2F *fHistDEDXdouble = (TH2F*)li->FindObject("fHistDEDXdouble");
+ TH1F *fHistCharge[4];
+ for(Int_t j=0;j<4;j++) fHistCharge[j] = (TH1F*)li->FindObject((Form("fHistChargeLay%d",j)));
+
+ TH1F *hEffPos[3];
+ TH1F *hEffNeg[3];
+ TH1F *hCorrFacPos[3];
+ TH1F *hCorrFacNeg[3];
+ TH1F *hEffMCPIDPos[3];
+ TH1F *hEffMCPIDNeg[3];
+ TH1F *hCorrFacMCPIDNeg[3];
+ TH1F *hCorrFacMCPIDPos[3];
+ TH1F *hSpectraPrimPosMC[3];
+ TH1F *hSpectraPrimNegMC[3];
+ TH1F *hSpectraPrimPosMCBefEvSel[3];
+ TH1F *hSpectraPrimNegMCBefEvSel[3];
+ TH1F *hSpectraPos[3];
+ TH1F *hSpectraNeg[3];
+ TH1F *hSpectraMCPIDPos[3];
+ TH1F *hSpectraMCPIDNeg[3];
+ TH1F* hMeanPos[3];
+ TH1F* hMeanNeg[3];
+ TH1F* hSigmaPos[3];
+ TH1F* hSigmaNeg[3];
+ TGraph *gres[6][22];
+
+ for(Int_t i=0; i<3; i++){
+ hSpectraPrimPosMC[i]=(TH1F*)li->FindObject(Form("fHistPrimMCpos%d",i));
+ hSpectraPrimNegMC[i]=(TH1F*)li->FindObject(Form("fHistPrimMCneg%d",i));
+ hSpectraPrimPosMCBefEvSel[i]=(TH1F*)li->FindObject(Form("fHistPrimMCposBefEvSel%d",i));
+ hSpectraPrimNegMCBefEvSel[i]=(TH1F*)li->FindObject(Form("fHistPrimMCnegBefEvSel%d",i));
+ hSpectraMCPIDPos[i]=new TH1F(Form("hSpectraMCPIDPos%d",i),Form("hSpectraMCPIDPos%d",i),nbins,xbins);
+ hSpectraMCPIDNeg[i]=new TH1F(Form("hSpectraMCPIDNeg%d",i),Form("hSpectraMCPIDNeg%d",i),nbins,xbins);
+ hSpectraPos[i]=new TH1F(Form("hSpectraPos%d",i),Form("hSpectraPos%d",i),nbins,xbins);
+ hSpectraNeg[i]=new TH1F(Form("hSpectraNeg%d",i),Form("hSpectraNeg%d",i),nbins,xbins);
+ hMeanPos[i]=new TH1F(Form("hMeanPos%d",i),Form("hMeanPos%d",i),nbins,xbins);
+ hMeanNeg[i]=new TH1F(Form("hMeanNeg%d",i),Form("hMeanNeg%d",i),nbins,xbins);
+ hSigmaPos[i]=new TH1F(Form("hSigmaPos%d",i),Form("hSigmaPos%d",i),nbins,xbins);
+ hSigmaNeg[i]=new TH1F(Form("hSigmaNeg%d",i),Form("hSigmaNeg%d",i),nbins,xbins);
+ }
+
+ //division for DeltaPt (for MC spectra)
+ for(Int_t ipt=0;ipt<3;ipt++){
+ for(Int_t bin=1; bin <= hSpectraPrimPosMC[ipt]->GetNbinsX(); bin++){
+ Float_t binSize=hSpectraPrimPosMC[ipt]->GetBinLowEdge(bin+1) - hSpectraPrimPosMC[ipt]->GetBinLowEdge(bin);
+ hSpectraPrimPosMC[ipt]->SetBinContent(bin, hSpectraPrimPosMC[ipt]->GetBinContent(bin) / binSize);
+ hSpectraPrimPosMC[ipt]->SetBinError(bin, 0);
+ hSpectraPrimNegMC[ipt]->SetBinContent(bin, hSpectraPrimNegMC[ipt]->GetBinContent(bin) / binSize);
+ hSpectraPrimNegMC[ipt]->SetBinError(bin, 0);
+ if(hSpectraPrimPosMCBefEvSel[ipt]){
+ hSpectraPrimPosMCBefEvSel[ipt]->SetBinContent(bin, hSpectraPrimPosMCBefEvSel[ipt]->GetBinContent(bin) / binSize);
+ hSpectraPrimPosMCBefEvSel[ipt]->SetBinError(bin, 0);
+ }
+ if(hSpectraPrimNegMCBefEvSel[ipt]){
+ hSpectraPrimNegMCBefEvSel[ipt]->SetBinContent(bin, hSpectraPrimNegMCBefEvSel[ipt]->GetBinContent(bin) / binSize);
+ hSpectraPrimNegMCBefEvSel[ipt]->SetBinError(bin, 0);
+ }
+ }
+ }
+ cout<<"All plots loaded"<<endl;
+
+ //open output file to store the dedx distribution and the spectra
+ TString savename=filename.Data();
+ savename.ReplaceAll("AnalysisResults","SpectraReco");
+ TFile *fout=new TFile(savename.Data(),"recreate");
+ fout->cd();
+
+ gStyle->SetOptStat(0);
+ gStyle->SetOptFit(111);
+ gStyle->SetPalette(1);
+
+ //propaganda plot
+ Float_t pdgmass[5]={0.13957,0.493677,0.938272,1.875612762,0.00596}; //mass for pi, K, P, d (Gev/c^2)
+ TF1 *funPos[5];
+ TF1 *funNeg[5];
+ for(Int_t m=0;m<5;m++){
+ funPos[m] = new TF1(Form("funPos%d",m),BetheBloch,0.02,5,1);
+ funPos[m]->SetParameter(0,pdgmass[m]);
+ funPos[m]->SetLineWidth(2);
+ funNeg[m] = new TF1(Form("funNeg%d",m),BetheBloch,-5,0.02,1);
+ funNeg[m]->SetParameter(0,-pdgmass[m]);
+ funNeg[m]->SetLineWidth(2);
+ }
+
+ TCanvas *cEvents=new TCanvas("cEvents","cEvents",500,400);
+ cEvents->cd();
+ fHistNEvents->Draw("text");
+ fHistNEvents->Write();
+
+ TCanvas *cDEDX=new TCanvas("cDEDX","DEDX",1000,800);
+ cDEDX->cd();
+ gPad->SetLogx();
+ gPad->SetLogz();
+ fHistDEDX->GetXaxis()->SetRangeUser(0.08,5);
+ fHistDEDX->GetYaxis()->SetRangeUser(0.,700);
+ fHistDEDX->GetXaxis()->SetTitle("momentum [GeV/c]");
+ fHistDEDX->GetYaxis()->SetTitle("dE [keV/300#mum]");
+ fHistDEDX->Draw("colz");
+ fHistDEDX->Write();
+ for(Int_t m=0;m<3;m++) funPos[m]->Draw("same");
+ Logo();
+
+ TCanvas *cDEDXdouble=new TCanvas("cDEDXdouble","DEDXdouble",1000,800);
+ cDEDXdouble->cd();
+ gPad->SetLogz();
+ fHistDEDXdouble->GetXaxis()->SetRangeUser(-3,3);
+ fHistDEDXdouble->GetXaxis()->SetTitle("momentum * sign [GeV/c]");
+ fHistDEDXdouble->GetYaxis()->SetTitle("dE [keV/300#mum]");
+ fHistDEDXdouble->Draw("colz");
+ fHistDEDXdouble->Write();
+ for(Int_t m=0;m<3;m++) {
+ funPos[m]->Draw("same");
+ funNeg[m]->Draw("same");
+ }
+ Logo();
+
+ //calibration check histo
+ TCanvas *cs=new TCanvas("cs","cs",1000,800);
+ cs->Divide(2,2);
+ for(Int_t j=0;j<4;j++){
+ cs->cd(j+1);
+ fHistCharge[j]->GetXaxis()->SetTitle("Charge [keV]");
+ if(j==0) fHistCharge[j]->SetTitle("Drift inner layer");
+ if(j==1) fHistCharge[j]->SetTitle("Drift outer layer");
+ if(j==2) fHistCharge[j]->SetTitle("Strip inner layer");
+ if(j==3) fHistCharge[j]->SetTitle("Strip inner layer");
+ fHistCharge[j]->Draw();
+ fHistCharge[j]->SetFillColor(7);
+ fHistCharge[j]->Fit("gaus","QR","",70,100);
+ fHistCharge[j]->Write();
+ }
+
+ //canvas MC spectra
+ if(optMC){
+ TCanvas *cspectraMC=new TCanvas("cspectraMC","cspectraMC",1000,800);
+ cspectraMC->Divide(2,1);
+ cspectraMC->cd(1);
+ gPad->SetLogy();
+ gPad->SetGridy();
+ for(Int_t i=0;i<3;i++){
+ if(i==0){
+ hSpectraPrimPosMC[i]->Draw("");
+ hSpectraPrimPosMC[i]->SetTitle("ITS EvMC truth positive");
+ hSpectraPrimPosMC[i]->SetMinimum(100);
+ hSpectraPrimPosMC[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
+ hSpectraPrimPosMC[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ }
+ hSpectraPrimPosMC[i]->Draw("same");
+ hSpectraPrimPosMC[i]->SetLineColor(i+2);
+ hSpectraPrimPosMC[i]->SetMarkerColor(i+2);
+ hSpectraPrimPosMC[i]->SetMarkerStyle(21+i);
+ hSpectraPrimPosMC[i]->Write();
+ }
+ cspectraMC->cd(2);
+ gPad->SetLogy();
+ gPad->SetGridy();
+ for(Int_t i=0;i<3;i++){
+ if(i==0){
+ hSpectraPrimNegMC[i]->Draw("");
+ hSpectraPrimNegMC[i]->SetTitle("ITS MC truth negative");
+ hSpectraPrimNegMC[i]->SetMinimum(100);
+ hSpectraPrimNegMC[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
+ hSpectraPrimNegMC[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ }
+ hSpectraPrimNegMC[i]->Draw("same");
+ hSpectraPrimNegMC[i]->SetLineColor(i+2);
+ hSpectraPrimNegMC[i]->SetMarkerColor(i+2);
+ hSpectraPrimNegMC[i]->SetMarkerStyle(21+i);
+ hSpectraPrimNegMC[i]->Write();
+ }
+ cspectraMC->Print(ps0.Data());
+ }
+
+ //dedx distribution to be fitted
+ TCanvas *cgausPipos=new TCanvas("cgausPipos","PIONS pos",1000,800);
+ cgausPipos->Divide(5,4,0.001,0.001);
+ TCanvas *cgausKpos=new TCanvas("cgausKpos","KAONS pos",1000,800);
+ cgausKpos->Divide(5,4,0.001,0.001);
+ TCanvas *cgausPpos=new TCanvas("cgausPpos","PROTONS pos",1000,800);
+ cgausPpos->Divide(5,4,0.001,0.001);
+ TCanvas *cgausPineg=new TCanvas("cgausPineg","PIONS neg",1000,800);
+ cgausPineg->Divide(5,4,0.001,0.001);
+ TCanvas *cgausKneg=new TCanvas("cgausKneg","KAONS neg",1000,800);
+ cgausKneg->Divide(5,4,0.001,0.001);
+ TCanvas *cgausPneg=new TCanvas("cgausPneg","PROTONS neg",1000,800);
+ cgausPneg->Divide(5,4,0.001,0.001);
+
+ //binsize for one histo only because are all equal
+ Float_t binsize = fHistPosPi[0]->GetBinWidth(1);
+ Double_t fpar[5],efpar[5];
+ for(Int_t i=0; i<nbins-2; i++){
+ //------------------- positive particles
+ //pions
+ cgausPipos->cd(i+1);
+ gPad->SetLogy();
+ fHistPosPi[i]->Write();
+ gres[0][i]=new TGraph();
+ ITSsa->DoFit(fHistPosPi[i],i,211,optMC,gres[0][i]);
+ ITSsa->FillHisto(hSpectraPos[0],i+1,binsize,211);
+ if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDPos[0],i+1,211,fHistMCPosPi[i]);
+ if(ITSsa->IsGoodBin(i,211)){
+ ITSsa->GetFitPar(fpar,efpar);
+ hMeanPos[0]->SetBinContent(i+1,fpar[1]);
+ hMeanPos[0]->SetBinError(i+1,efpar[1]);
+ hSigmaPos[0]->SetBinContent(i+1,fpar[2]);
+ hSigmaPos[0]->SetBinError(i+1,efpar[2]);
+ }
+ cgausPipos->Update();
+ //kaons
+ cgausKpos->cd(i+1);
+ gPad->SetLogy();
+ fHistPosK[i]->Write();
+ gres[1][i]=new TGraph();
+ ITSsa->DoFit(fHistPosK[i],i,321,optMC,gres[1][i]);
+ ITSsa->FillHisto(hSpectraPos[1],i+1,binsize,321);
+ if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDPos[1],i+1,321,fHistMCPosK[i]);
+ if(ITSsa->IsGoodBin(i,321)){
+ ITSsa->GetFitPar(fpar,efpar);
+ hMeanPos[1]->SetBinContent(i+1,fpar[1]);
+ hMeanPos[1]->SetBinError(i+1,efpar[1]);
+ hSigmaPos[1]->SetBinContent(i+1,fpar[2]);
+ hSigmaPos[1]->SetBinError(i+1,efpar[2]);
+ }
+ cgausKpos->Update();
+ //protons
+ cgausPpos->cd(i+1);
+ gPad->SetLogy();
+ fHistPosP[i]->Write();
+ fHistPosP[i]->SetFillColor(16);
+ gres[2][i]=new TGraph();
+ ITSsa->DoFitProton(fHistPosP[i],i,2212,optMC,gres[2][i]);
+ ITSsa->FillHisto(hSpectraPos[2],i+1,binsize,2212);
+ if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDPos[2],i+1,2212,fHistMCPosP[i]);
+ if(ITSsa->IsGoodBin(i,2212)){
+ ITSsa->GetFitPar(fpar,efpar);
+ hMeanPos[2]->SetBinContent(i+1,fpar[1]);
+ hMeanPos[2]->SetBinError(i+1,efpar[1]);
+ hSigmaPos[2]->SetBinContent(i+1,fpar[2]);
+ hSigmaPos[2]->SetBinError(i+1,efpar[2]);
+ }
+ cgausPpos->Update();
+
+ //------------------- negative particles
+ //pions
+ cgausPineg->cd(i+1);
+ gPad->SetLogy();
+ fHistNegPi[i]->Write();
+ gres[3][i]=new TGraph();
+ ITSsa->DoFit(fHistNegPi[i],i,211,optMC,gres[3][i]);
+ ITSsa->FillHisto(hSpectraNeg[0],i+1,binsize,211);
+ if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDNeg[0],i+1,211,fHistMCNegPi[i]);
+ if(ITSsa->IsGoodBin(i,211)){
+ ITSsa->GetFitPar(fpar,efpar);
+ hMeanNeg[0]->SetBinContent(i+1,fpar[1]);
+ hMeanNeg[0]->SetBinError(i+1,efpar[1]);
+ hSigmaNeg[0]->SetBinContent(i+1,fpar[2]);
+ hSigmaNeg[0]->SetBinError(i+1,efpar[2]);
+ }
+ cgausPineg->Update();
+ //kaons
+ cgausKneg->cd(i+1);
+ gPad->SetLogy();
+ fHistNegK[i]->Write();
+ gres[4][i]=new TGraph();
+ ITSsa->DoFit(fHistNegK[i],i,321,optMC,gres[4][i]);
+ ITSsa->FillHisto(hSpectraNeg[1],i+1,binsize,321);
+ if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDNeg[1],i+1,321,fHistMCNegK[i]);
+ if(ITSsa->IsGoodBin(i,321)){
+ ITSsa->GetFitPar(fpar,efpar);
+ hMeanNeg[1]->SetBinContent(i+1,fpar[1]);
+ hMeanNeg[1]->SetBinError(i+1,efpar[1]);
+ hSigmaNeg[1]->SetBinContent(i+1,fpar[2]);
+ hSigmaNeg[1]->SetBinError(i+1,efpar[2]);
+ }
+ cgausKneg->Update();
+ //protons
+ cgausPneg->cd(i+1);
+ gPad->SetLogy();
+ fHistNegP[i]->Write();
+ gres[5][i]=new TGraph();
+ ITSsa->DoFitProton(fHistNegP[i],i,2212,optMC,gres[5][i]);
+ ITSsa->FillHisto(hSpectraNeg[2],i+1,binsize,2212);
+ if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDNeg[2],i+1,2212,fHistMCNegP[i]);
+ if(ITSsa->IsGoodBin(i,2212)){
+ ITSsa->GetFitPar(fpar,efpar);
+ hMeanNeg[2]->SetBinContent(i+1,fpar[1]);
+ hMeanNeg[2]->SetBinError(i+1,efpar[1]);
+ hSigmaNeg[2]->SetBinContent(i+1,fpar[2]);
+ hSigmaNeg[2]->SetBinError(i+1,efpar[2]);
+ }
+ cgausPneg->Update();
+ }
+
+ //save histograms in the ps file
+ cgausPipos->Print(ps1.Data());
+ cgausKpos->Print(ps1.Data());
+ cgausPpos->Print(ps1.Data());
+ cgausPineg->Print(ps1.Data());
+ cgausKneg->Print(ps1.Data());
+ cgausPneg->Print(ps1.Data());
+
+ //spectra REC
+ TCanvas *cspecREC=new TCanvas("cspecREC","SPECTRA rec",1000,800);
+ cspecREC->Divide(2,1);
+ cspecREC->cd(1);
+ gPad->SetLogy();
+ gPad->SetGridy();
+ for(Int_t i=0;i<3;i++){
+ if(i==0){
+ hSpectraPos[i]->Draw("text");
+ hSpectraPos[i]->SetTitle("ITSsa RAW positive");
+ hSpectraPos[i]->SetMinimum(10);
+ hSpectraPos[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
+ hSpectraPos[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ }
+ hSpectraPos[i]->Draw("esame");
+ hSpectraPos[i]->SetLineColor(i+2);
+ hSpectraPos[i]->SetMarkerColor(i+2);
+ hSpectraPos[i]->SetMarkerStyle(21+i);
+ hSpectraPos[i]->Write();
+ }
+ TLegend *leg=new TLegend(0.51,0.11,0.84,0.35,NULL,"brNDC");
+ leg->SetFillColor(0);
+ leg->SetBorderSize(0);
+ TLegendEntry *entry0=leg->AddEntry(hSpectraPos[0],"pions","p");
+ entry0->SetTextColor(2);
+ TLegendEntry *entry2=leg->AddEntry(hSpectraPos[1],"kaons","p");
+ entry2->SetTextColor(3);
+ TLegendEntry *entry4=leg->AddEntry(hSpectraPos[2],"protons","p");
+ entry4->SetTextColor(4);
+ leg->Draw("same");
+
+ cspecREC->cd(2);
+ gPad->SetLogy();
+ gPad->SetGridy();
+ for(Int_t i=0;i<3;i++){
+ if(i==0){
+ hSpectraNeg[i]->Draw("e");
+ hSpectraNeg[i]->SetTitle("ITSsa RAW negative");
+ hSpectraNeg[i]->SetMinimum(10);
+ hSpectraNeg[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
+ hSpectraNeg[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ }
+ hSpectraNeg[i]->Draw("esame");
+ hSpectraNeg[i]->SetLineColor(i+2);
+ hSpectraNeg[i]->SetMarkerColor(i+2);
+ hSpectraNeg[i]->SetMarkerStyle(21+i);
+ hSpectraNeg[i]->Write();
+ }
+
+
+ //FitParameters
+ TCanvas *cFitPar=new TCanvas("cFitPar","FitParameters",1000,800);
+ cFitPar->Divide(2,3);
+ TLatex** tplus=new TLatex*[3];
+ TLatex** tminus=new TLatex*[3];
+ tplus[0]=new TLatex(0.7,0.75,"#pi^{+}");
+ tplus[1]=new TLatex(0.7,0.75,"K^{+}");
+ tplus[2]=new TLatex(0.7,0.75,"p");
+ tminus[0]=new TLatex(0.7,0.65,"#pi^{-}");
+ tminus[1]=new TLatex(0.7,0.65,"K^{-}");
+ tminus[2]=new TLatex(0.7,0.65,"#bar{p}");
+
+ for(Int_t ipart=0; ipart<3; ipart++){
+ tplus[ipart]->SetNDC();
+ tplus[ipart]->SetTextColor(1);
+ tplus[ipart]->SetTextSize(0.05);
+ tminus[ipart]->SetNDC();
+ tminus[ipart]->SetTextColor(4);
+ tminus[ipart]->SetTextSize(0.05);
+ cFitPar->cd(1+2*ipart);
+ gPad->SetGridy();
+ hMeanPos[ipart]->SetMarkerStyle(20);
+ hMeanPos[ipart]->Draw("e");
+ hMeanPos[ipart]->GetYaxis()->SetRangeUser(-1,1);
+ hMeanNeg[ipart]->SetMarkerStyle(24);
+ hMeanNeg[ipart]->SetMarkerColor(4);
+ hMeanNeg[ipart]->SetLineColor(4);
+ hMeanNeg[ipart]->Draw("esame");
+ tplus[ipart]->Draw();
+ tminus[ipart]->Draw();
+ cFitPar->cd(2+2*ipart);
+ gPad->SetGridy();
+ hSigmaPos[ipart]->SetMarkerStyle(20);
+ hSigmaPos[ipart]->Draw("e");
+ hSigmaPos[ipart]->GetYaxis()->SetRangeUser(0.1,0.3);
+ hSigmaNeg[ipart]->SetMarkerStyle(24);
+ hSigmaNeg[ipart]->SetMarkerColor(4);
+ hSigmaNeg[ipart]->SetLineColor(4);
+ hSigmaNeg[ipart]->Draw("esame");
+ hSigmaPos[ipart]->Write();
+ hSigmaNeg[ipart]->Write();
+ tplus[ipart]->Draw();
+ tminus[ipart]->Draw();
+ }
+
+ //plus/minus ratio plots
+ TH1F* hRatioPions=new TH1F(*hSpectraPos[0]);
+ hRatioPions->Divide(hSpectraNeg[0]);
+ TH1F* hRatioKaons=new TH1F(*hSpectraPos[1]);
+ hRatioKaons->Divide(hSpectraNeg[1]);
+ TH1F* hRatioProtons=new TH1F(*hSpectraPos[2]);
+ hRatioProtons->Divide(hSpectraNeg[2]);
+
+ TCanvas *cratios=new TCanvas("cratios","Ratios +/-",1000,800);
+ cratios->Divide(1,3);
+ cratios->cd(1);
+ hRatioPions->SetMinimum(0.7);
+ hRatioPions->SetMaximum(1.3);
+ hRatioPions->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ hRatioPions->GetYaxis()->SetTitle("#pi^{+}/#pi^{-}");
+ hRatioPions->Draw("e");
+ cratios->cd(2);
+ hRatioKaons->SetMinimum(0.7);
+ hRatioKaons->SetMaximum(1.3);
+ hRatioKaons->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ hRatioKaons->GetYaxis()->SetTitle("K^{+}/K^{-}");
+ hRatioKaons->Draw("e");
+ cratios->cd(3);
+ hRatioProtons->SetMinimum(0.7);
+ hRatioProtons->SetMaximum(1.3);
+ hRatioProtons->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ hRatioProtons->GetYaxis()->SetTitle("p/#bar{p}");
+ hRatioProtons->Draw("e");
+
+ //efficiency and correction factor histograms
+ if(optMC){
+ for(Int_t i=0;i<3;i++){
+ hEffPos[i] = (TH1F*)hSpectraPos[i]->Clone(Form("hEffPos%d",i));
+ hEffPos[i]->SetTitle(Form("hEffPos%d",i));
+ hEffPos[i]->Divide(hEffPos[i], hSpectraPrimPosMC[i], 1.0, 1.0, "B");//binomial errors
+ hEffPos[i]->SetLineColor(i+2);
+ hEffPos[i]->SetMarkerColor(i+2);
+ hEffPos[i]->SetMarkerStyle(21+i);
+ hEffPos[i]->Write();
+
+ if(hSpectraPrimPosMCBefEvSel[i]){
+ hCorrFacPos[i] = (TH1F*)hSpectraPos[i]->Clone(Form("hCorrFacPos%d",i));
+ hCorrFacPos[i]->SetTitle(Form("hCorrFacPos%d",i));
+ hCorrFacPos[i]->Divide(hCorrFacPos[i], hSpectraPrimPosMCBefEvSel[i], 1.0, 1.0, "B");//binomial errors
+ hCorrFacPos[i]->SetLineColor(i+2);
+ hCorrFacPos[i]->SetMarkerColor(i+2);
+ hCorrFacPos[i]->SetMarkerStyle(25+i);
+ hCorrFacPos[i]->Write();
+ }
+ hEffMCPIDPos[i] = (TH1F*)hSpectraMCPIDPos[i]->Clone(Form("hEffMCPIDPos%d",i));
+ hEffMCPIDPos[i]->SetTitle(Form("hEffMCPIDPos%d",i));
+ hEffMCPIDPos[i]->Divide(hEffMCPIDPos[i], hSpectraPrimPosMC[i], 1.0, 1.0, "B");//binomial errors
+ hEffMCPIDPos[i]->SetLineColor(i+2);
+ hEffMCPIDPos[i]->SetMarkerColor(i+2);
+ hEffMCPIDPos[i]->SetMarkerStyle(25+i);
+ hEffMCPIDPos[i]->Write();
+
+ if(hSpectraPrimPosMCBefEvSel[i]){
+ hCorrFacMCPIDPos[i] = (TH1F*)hSpectraMCPIDPos[i]->Clone(Form("hCorrFacMCPIDPos%d",i));
+ hCorrFacMCPIDPos[i]->SetTitle(Form("hCorrFacMCPIDPos%d",i));
+ hCorrFacMCPIDPos[i]->Divide(hCorrFacMCPIDPos[i], hSpectraPrimPosMCBefEvSel[i], 1.0, 1.0, "B");//binomial errors
+ hCorrFacMCPIDPos[i]->SetLineColor(i+2);
+ hCorrFacMCPIDPos[i]->SetMarkerColor(i+2);
+ hCorrFacMCPIDPos[i]->SetMarkerStyle(25+i);
+ hCorrFacMCPIDPos[i]->Write();
+ }
+
+ hEffNeg[i] = (TH1F*)hSpectraNeg[i]->Clone(Form("hEffNeg%d",i));
+ hEffNeg[i]->SetTitle(Form("hEffNeg%d",i));
+ hEffNeg[i]->Divide(hEffNeg[i], hSpectraPrimNegMC[i], 1.0, 1.0, "B");//binomial errors
+ hEffNeg[i]->SetLineColor(i+2);
+ hEffNeg[i]->SetMarkerColor(i+2);
+ hEffNeg[i]->SetMarkerStyle(21+i);
+ hEffNeg[i]->Write();
+
+ if(hSpectraPrimNegMCBefEvSel[i]){
+ hCorrFacNeg[i] = (TH1F*)hSpectraNeg[i]->Clone(Form("hCorrFacNeg%d",i));
+ hCorrFacNeg[i]->SetTitle(Form("hCorrFacNeg%d",i));
+ hCorrFacNeg[i]->Divide(hCorrFacNeg[i], hSpectraPrimNegMCBefEvSel[i], 1.0, 1.0, "B");//binomial errors
+ hCorrFacNeg[i]->SetLineColor(i+2);
+ hCorrFacNeg[i]->SetMarkerColor(i+2);
+ hCorrFacNeg[i]->SetMarkerStyle(25+i);
+ hCorrFacNeg[i]->Write();
+ }
+
+ hEffMCPIDNeg[i] = (TH1F*)hSpectraMCPIDNeg[i]->Clone(Form("hEffMCPIDNeg%d",i));
+ hEffMCPIDNeg[i]->SetTitle(Form("hEffMCPIDNeg%d",i));
+ hEffMCPIDNeg[i]->Divide(hEffMCPIDNeg[i], hSpectraPrimNegMC[i], 1.0, 1.0, "B");//binomial errors
+ hEffMCPIDNeg[i]->SetLineColor(i+2);
+ hEffMCPIDNeg[i]->SetMarkerColor(i+2);
+ hEffMCPIDNeg[i]->SetMarkerStyle(25+i);
+ hEffMCPIDNeg[i]->Write();
+
+ if(hSpectraPrimPosMCBefEvSel[i]){
+ hCorrFacMCPIDNeg[i] = (TH1F*)hSpectraMCPIDPos[i]->Clone(Form("hCorrFacMCPIDNeg%d",i));
+ hCorrFacMCPIDNeg[i]->SetTitle(Form("hCorrFacMCPIDNeg%d",i));
+ hCorrFacMCPIDNeg[i]->Divide(hCorrFacMCPIDNeg[i], hSpectraPrimNegMCBefEvSel[i], 1.0, 1.0, "B");//binomial errors
+ hCorrFacMCPIDNeg[i]->SetLineColor(i+2);
+ hCorrFacMCPIDNeg[i]->SetMarkerColor(i+2);
+ hCorrFacMCPIDNeg[i]->SetMarkerStyle(25+i);
+ hCorrFacMCPIDNeg[i]->Write();
+ }
+
+ hSpectraMCPIDPos[i]->SetLineColor(i+2);
+ hSpectraMCPIDPos[i]->SetMarkerColor(i+2);
+ hSpectraMCPIDPos[i]->SetMarkerStyle(25+i);
+ hSpectraMCPIDPos[i]->Write();
+
+ hSpectraMCPIDNeg[i]->SetLineColor(i+2);
+ hSpectraMCPIDNeg[i]->SetMarkerColor(i+2);
+ hSpectraMCPIDNeg[i]->SetMarkerStyle(25+i);
+ hSpectraMCPIDNeg[i]->Write();
+ }
+
+ //spectra REC Ideal PID
+ TCanvas *cspecREC2=new TCanvas("cspecIdPid","SPECTRA rec Ideal PID",1000,800);
+ cspecREC2->Divide(2,1);
+ cspecREC2->cd(1);
+ gPad->SetLogy();
+ gPad->SetGridy();
+ for(Int_t i=0;i<3;i++){
+ if(i==0){
+ hSpectraMCPIDPos[i]->Draw("text");
+ hSpectraMCPIDPos[i]->SetTitle("ITSsa RAW positive - Ideal PID");
+ hSpectraMCPIDPos[i]->SetMinimum(10);
+ hSpectraMCPIDPos[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
+ hSpectraMCPIDPos[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ }
+ hSpectraMCPIDPos[i]->Draw("esame");
+ }
+ leg->Draw("same");
+ cspecREC2->cd(2);
+ gPad->SetLogy();
+ gPad->SetGridy();
+ for(Int_t i=0;i<3;i++){
+ if(i==0){
+ hSpectraMCPIDNeg[i]->Draw("e");
+ hSpectraMCPIDNeg[i]->SetTitle("ITSsa RAW negative -Ideal PID");
+ hSpectraMCPIDNeg[i]->SetMinimum(10);
+ hSpectraMCPIDNeg[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
+ hSpectraMCPIDNeg[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ }
+ hSpectraMCPIDNeg[i]->Draw("esame");
+ }
+
+ // Efficiency
+ TCanvas *ceff=new TCanvas("ceff","ceff",1000,800);
+ ceff->Divide(2,1);
+ ceff->cd(1);
+ gPad->SetGridy();
+ for(Int_t i=0;i<3;i++){
+ if(i==0){
+ hEffPos[i]->SetTitle("ITSsa efficiency (positive)");
+ hEffPos[i]->Draw("e");
+ hEffPos[i]->GetYaxis()->SetRangeUser(0.1,0.9);
+ hEffPos[i]->GetYaxis()->SetTitle("#epsilon");
+ hEffPos[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ }
+ hEffPos[i]->Draw("esame");
+ hEffMCPIDPos[i]->Draw("esame");
+ }
+ TLegend *legeff=new TLegend(0.51,0.11,0.89,0.5,NULL,"brNDC");
+ legeff->SetBorderSize(0);
+ legeff->SetFillColor(0);
+ TLegendEntry *entryeff0=legeff->AddEntry(hEffPos[0],"pions from fit","p");
+ entryeff0->SetTextColor(2);
+ TLegendEntry *entryeff2=legeff->AddEntry(hEffPos[1],"kaons from fit","p");
+ entryeff2->SetTextColor(3);
+ TLegendEntry *entryeff4=legeff->AddEntry(hEffPos[2],"protons from fit","p");
+ entryeff4->SetTextColor(4);
+ TLegendEntry *entryeff1=legeff->AddEntry(hEffMCPIDPos[0],"pions ideal PID","p");
+ entryeff1->SetTextColor(2);
+ TLegendEntry *entryeff3=legeff->AddEntry(hEffMCPIDPos[1],"kaons ideal PID","p");
+ entryeff3->SetTextColor(3);
+ TLegendEntry *entryeff5=legeff->AddEntry(hEffMCPIDPos[2],"protons ideal PID","p");
+ entryeff5->SetTextColor(4);
+ legeff->Draw("same");
+
+ ceff->cd(2);
+ gPad->SetGridy();
+ for(Int_t i=0;i<3;i++){
+ if(i==0){
+ hEffNeg[i]->SetTitle("ITSsa efficiency (negative)");
+ hEffNeg[i]->Draw("e");
+ hEffNeg[i]->GetYaxis()->SetRangeUser(0.1,0.9);
+ hEffNeg[i]->GetYaxis()->SetTitle("#epsilon");
+ hEffNeg[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ }
+ hEffNeg[i]->Draw("esame");
+ hEffMCPIDNeg[i]->Draw("esame");
+ }
+ ceff->Update();
+
+ // Correction factor
+ TCanvas *ccf=new TCanvas("ccf","Corr Factor",1000,800);
+ ccf->Divide(2,1);
+ ccf->cd(1);
+ gPad->SetGridy();
+ for(Int_t i=0;i<3;i++){
+ if(i==0){
+ hEffPos[i]->SetTitle("ITSsa efficiency (positive)");
+ hEffPos[i]->Draw("e");
+ hEffPos[i]->GetYaxis()->SetRangeUser(0.1,0.9);
+ hEffPos[i]->GetYaxis()->SetTitle("#epsilon");
+ hEffPos[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ }
+ hEffPos[i]->Draw("esame");
+ if(hSpectraPrimPosMCBefEvSel[i]) hCorrFacPos[i]->Draw("esame");
+ }
+ TLegend *legcf=new TLegend(0.51,0.11,0.89,0.5,NULL,"brNDC");
+ legcf->SetBorderSize(0);
+ legcf->SetFillColor(0);
+ TLegendEntry *entrycf0=legcf->AddEntry(hEffPos[0],"pions - eff.","p");
+ entrycf0->SetTextColor(2);
+ TLegendEntry *entrycf2=legcf->AddEntry(hEffPos[1],"kaons - eff.","p");
+ entrycf2->SetTextColor(3);
+ TLegendEntry *entrycf4=legcf->AddEntry(hEffPos[2],"protons - eff.","p");
+ entrycf4->SetTextColor(4);
+ TLegendEntry *entrycf1=legcf->AddEntry(hCorrFacPos[0],"pions - Corr. Fac.","p");
+ entrycf1->SetTextColor(2);
+ TLegendEntry *entrycf3=legcf->AddEntry(hCorrFacPos[1],"kaons - Corr. Fac.","p");
+ entrycf3->SetTextColor(3);
+ TLegendEntry *entrycf5=legcf->AddEntry(hCorrFacPos[2],"protons - Corr. Fac.","p");
+ entrycf5->SetTextColor(4);
+ legcf->Draw("same");
+
+ ccf->cd(2);
+ gPad->SetGridy();
+ for(Int_t i=0;i<3;i++){
+ if(i==0){
+ hEffNeg[i]->SetTitle("ITSsa efficiency (negative)");
+ hEffNeg[i]->Draw("e");
+ hEffNeg[i]->GetYaxis()->SetRangeUser(0.1,0.9);
+ hEffNeg[i]->GetYaxis()->SetTitle("#epsilon");
+ hEffNeg[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
+ }
+ hEffNeg[i]->Draw("esame");
+ if(hSpectraPrimNegMCBefEvSel[i]) hCorrFacNeg[i]->Draw("esame");
+ }
+ ccf->Update();
+
+ cspecREC2->Print(ps1.Data());
+ ceff->Print(ps1.Data());
+ ccf->Print(ps1.Data());
+ }
+
+ //save canvas in the ps file
+ cratios->Print(ps1.Data());
+ cspecREC->Print(ps1.Data());
+ cFitPar->Print(ps1.Data());
+ cs->Print(ps1.Data());
+ cDEDX->Print(ps1.Data());
+ cDEDXdouble->Print(ps1.Data());
+ cdummy->Print(ps2.Data());
+ // fout->Close();
+ return;
+
+}//end of the main
+
+
+//______________________________________________________________________
+Double_t BetheBloch(Double_t *mom, Double_t *mass) {
+ // BB PHobos parametrization fine tuned by G. Ortona
+ // on data and MC in June 2010
+ Double_t bg=mom[0]/mass[0];
+ Double_t par[5];
+ Bool_t MC=kFALSE;
+ if(MC){
+ par[0]=1.39126e+02;//tuned on LHC10d (sim pass4)
+ par[1]=2.33589e+01;
+ par[2]=6.05219e-02;
+ par[3]=2.04336e-01;
+ par[4]=-4.99854e-04;
+ }
+ else{
+ par[0]=5.33458e+04;//tuned on data pass6
+ par[1]=1.65303e+01;
+ par[2]=2.60065e-03;
+ par[3]=3.59533e-04;
+ par[4]=7.51168e-05;
+ }
+ Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
+ Double_t gamma=bg/beta;
+ Double_t eff=1.0;
+ if(bg<par[2]) eff=(bg-par[3])*(bg-par[3])+par[4];
+ else eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
+ return (par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
+}
+
+//______________________________________________________________________
+void Logo(){
+ //method to plot the ALICE logo in the propaganda plot
+ TPaveText *tpave=new TPaveText(0.3,0.7,0.7,0.89,"brNDC");
+ tpave->SetBorderSize(0);
+ tpave->SetFillStyle(0);
+ tpave->SetFillColor(0);
+ TText *txt1=tpave->AddText("ALICE performance");
+ TText *txt2=tpave->AddText("ITS stand-alone tracks");
+ TText *txt3=tpave->AddText("pp @ #sqrt{s} = 7 TeV");
+ txt1->SetTextFont(62);
+ txt2->SetTextFont(62);
+ txt3->SetTextFont(62);
+ tpave->Draw();
+ TPad *pad1=new TPad("pad1", "pad1", 0.75, 0.7, 0.89, 0.89);
+ pad1->Draw();
+ pad1->cd();
+ TImage *img1 = TImage::Open("ALICElogo.gif");
+ img1->FillRectangle(0);
+ img1->Draw();
+}
+
+//EOF
+
SPECTRA/AliAnalysisCentralExtrapolate.cxx \
SPECTRA/AliAnalysisTaskCentral.cxx \
SPECTRA/AliAnalysisTaskSEITSsaSpectra.cxx \
+ SPECTRA/AliITSsadEdxFitter.cxx \
SPECTRA/Fit/AliBWTools.cxx \
SPECTRA/Fit/AliBWFunc.cxx \
SPECTRA/Fit/AliLatexTable.cxx \
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff