[u/mrichter/AliRoot.git] / ITS / ParametrizeITSBetheBloch.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TCanvas.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TGraphAsymmErrors.h>
#include <TH1.h>
#include <TF1.h>
#include <TH2.h>
#include <TFile.h>
#include <TList.h>
#include <TString.h>
#include <TDirectoryFile.h>
#include <TLatex.h>
#include <TStyle.h>
#include <TMath.h>
#include <TDatabasePDG.h>
#include "AliITSPIDResponse.h"
#endif

/*  $Id$    */


// Macro to extract the parameterizations of dE/dx in ITS vs. p
// using Phobos BB parameterization


Double_t BetheBlochPhobos(const Double_t *x, const Double_t *par);
Double_t BetheBlochPhobosBG(const Double_t *x, const Double_t *par);
Double_t BetheBlochAleph(const Double_t *x, const Double_t *par);


void ParametrizeITSBetheBloch(Bool_t optFromITSsa=kFALSE){

  TFile *fil=new TFile("AnalysisResults.root");
  TH2F *hsigvsp=0x0;
  TString dfName="PWG2SpectraITSsa";
  TString listName="clistITSsaMult-1to-1";
  TString hisName="fHistDEDX";

  if(optFromITSsa){
    dfName="ITSsaTracks";
    listName="clistITSsaTracks";
    hisName="hdedxvsPt4clsITSpureSA";
  }
  TDirectoryFile* df=(TDirectoryFile*)fil->Get(dfName.Data());
  if(!df){
    printf("ERROR: directory file %s not found\n",dfName.Data());
    return;
  }
  TList* l=(TList*)df->Get(listName.Data());
  if(!l){
    printf("ERROR: TList %s not found\n",listName.Data());
    return;
  }
  hsigvsp=(TH2F*)l->FindObject(hisName.Data());
  if(!hsigvsp){ 
    printf("ERROR: historgam %s not found\n",hisName.Data());
    return;
  }
  hsigvsp->GetXaxis()->SetTitle("momentum (GeV/c)");
  hsigvsp->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");

  Double_t minpPionFit=0.1;
  Double_t maxpPionFit=0.8;
  Double_t minbgGlobFit=0.3;
  Double_t maxbgGlobFit=2.5;
  Double_t minbgPions=2.;  
  Double_t maxbgPions=10.;  
  Double_t minpKaons=0.2;
  Double_t maxpKaons=0.4;
  Double_t maxbgKaons=0.7;
  Double_t minpProtons=0.35;
  Double_t maxpProtons=0.57;
  Double_t maxbgProtons=0.7;

  Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
  Double_t massKaon=TDatabasePDG::Instance()->GetParticle(321)->Mass();
  Double_t massProton=TDatabasePDG::Instance()->GetParticle(2212)->Mass();
  Double_t massDeuteron=1.8756;

  Int_t iPad=1;
  Int_t iCanvas=0;
  Int_t nbinsX=hsigvsp->GetNbinsX();
  printf("Total Number of bins =%d\n",nbinsX);
  Int_t group=2;
 
  Int_t totPads=(Int_t)((Float_t)nbinsX/(Float_t)group+0.5);
  Int_t totCanvas=totPads/16;
  TCanvas** c0=new TCanvas*[totCanvas];
  c0[iCanvas]=new TCanvas(Form("c0_%d",iCanvas),Form("c0_%d",iCanvas));
  c0[iCanvas]->Divide(4,4);
  TCanvas** c0k=new TCanvas*[totCanvas];
  c0k[iCanvas]=new TCanvas(Form("c0k_%d",iCanvas),Form("c0k_%d",iCanvas));
  c0k[iCanvas]->Divide(4,4);
  TCanvas** c0p=new TCanvas*[totCanvas];
  c0p[iCanvas]=new TCanvas(Form("c0p_%d",iCanvas),Form("c0p_%d",iCanvas));
  c0p[iCanvas]->Divide(4,4);

  TGraphErrors* gdedxpi=new TGraphErrors(0);
  TGraphErrors* gdedxk=new TGraphErrors(0);
  TGraphErrors* gdedxp=new TGraphErrors(0);
  TGraphErrors* gdedxbg=new TGraphErrors(0);
  gdedxbg->SetTitle("");

  Int_t nPpi=0;
  Int_t nPk=0;
  Int_t nPp=0;
  Int_t nPbg=0;
  for(Int_t iBin=1; iBin<=nbinsX; iBin+=group){
    if(iPad>16){
      iCanvas++;
      c0[iCanvas]=new TCanvas(Form("c0_%d",iCanvas),Form("c0_%d",iCanvas));
      c0[iCanvas]->Divide(4,4);
      c0k[iCanvas]=new TCanvas(Form("c0k_%d",iCanvas),Form("c0k_%d",iCanvas));
      c0k[iCanvas]->Divide(4,4);
      c0p[iCanvas]=new TCanvas(Form("c0p_%d",iCanvas),Form("c0p_%d",iCanvas));
      c0p[iCanvas]->Divide(4,4);
      iPad=1;
    }

    Float_t sum=0;
    Float_t sumwei=0;
    Double_t pmed=0.;
    for(Int_t ibx=iBin; ibx<iBin+group;ibx++){
      for(Int_t iby=1; iby<=hsigvsp->GetNbinsY();iby++){
	sum+=hsigvsp->GetXaxis()->GetBinCenter(ibx)*hsigvsp->GetBinContent(ibx,iby);
	sumwei+=hsigvsp->GetBinContent(ibx,iby);
      }
    }
    if(sumwei>0.) pmed=sum/sumwei;
    if(pmed<0.1) continue;
    Int_t iPmed=(Int_t)(pmed*1000);
    TH1D* hProj=hsigvsp->TH2F::ProjectionY(Form("proj%d",iPmed),iBin,iBin+group-1);
    TH1D* hKaon=(TH1D*)hProj->Clone(Form("projk%d",iPmed));
    c0[iCanvas]->cd(iPad);
    hProj->SetLineColor(2);
    hProj->Draw();    

    Float_t peakpos=hProj->GetBinCenter(hProj->GetMaximumBin());
    Float_t peaksig=0.1*peakpos;
    hProj->Fit("gaus","LL","",peakpos-3*peaksig,peakpos+3.*peaksig);
    TF1* fgaus=(TF1*)hProj->GetListOfFunctions()->FindObject("gaus");
    gdedxpi->SetPoint(nPpi,sum/sumwei,fgaus->GetParameter(1));
    gdedxpi->SetPointError(nPpi,0.,fgaus->GetParError(1));
    nPpi++;
    if(pmed/massPion>minbgPions && pmed/massPion<maxbgPions){
      gdedxbg->SetPoint(nPbg,pmed/massPion,fgaus->GetParameter(1));
      gdedxbg->SetPointError(nPbg,0.,fgaus->GetParError(1));
      nPbg++;
    }
    for(Int_t iBinK=1; iBinK<=hKaon->GetNbinsX(); iBinK++){
      Double_t cont=hKaon->GetBinContent(iBinK)-fgaus->Eval(hKaon->GetBinCenter(iBinK));
      if(cont<0) cont=0.;
      if(hKaon->GetBinCenter(iBinK)<fgaus->GetParameter(1)) cont=0;
      hKaon->SetBinContent(iBinK,cont);
    }
    c0k[iCanvas]->cd(iPad);
    hKaon->SetLineColor(4);
    hKaon->Draw();
    TH1D* hProton=(TH1D*)hKaon->Clone(Form("projp%d",iPmed));
    TF1* fgausk=0x0;
    Double_t mink=250-600*(pmed-0.2);
    Double_t maxk=mink+mink;
    hKaon->Fit("gaus","LL","",mink,maxk);
    fgausk=(TF1*)hKaon->GetListOfFunctions()->FindObject("gaus");
    if(pmed>minpKaons && pmed<maxpKaons){
      gdedxk->SetPoint(nPk,pmed,fgausk->GetParameter(1));
      gdedxk->SetPointError(nPk,0.,fgausk->GetParError(1));
      nPk++;
      if(pmed/massKaon<maxbgKaons){
	gdedxbg->SetPoint(nPbg,pmed/massKaon,fgausk->GetParameter(1));
	gdedxbg->SetPointError(nPbg,0.,fgausk->GetParError(1));
	nPbg++;
      }    
    }
    if(fgausk){
      for(Int_t iBinP=1; iBinP<=hProton->GetNbinsX(); iBinP++){
	Double_t cont=hKaon->GetBinContent(iBinP)-fgausk->Eval(hKaon->GetBinCenter(iBinP));
	if(cont<0) cont=0.;
	if(hProton->GetBinCenter(iBinP)<fgausk->GetParameter(1)) cont=0;
	hProton->SetBinContent(iBinP,cont);
      }
    }
    c0p[iCanvas]->cd(iPad);
    hProton->SetLineColor(kGreen+1);
    hProton->Draw();
    if(pmed>minpProtons && pmed<maxpProtons){
      Double_t minP=300-600*(pmed-0.35);
      Double_t maxP=minP+minP;
      hProton->Fit("gaus","LL","",minP,maxP);
      TF1* fgausp=(TF1*)hProton->GetListOfFunctions()->FindObject("gaus");
      gdedxp->SetPoint(nPp,pmed,fgausp->GetParameter(1));
      gdedxp->SetPointError(nPp,0.,fgausp->GetParError(1));
      nPp++;
      if(pmed/massProton<maxbgProtons){
	gdedxbg->SetPoint(nPbg,pmed/massProton,fgausp->GetParameter(1));
	gdedxbg->SetPointError(nPbg,0.,fgausp->GetParError(1));
	nPbg++;
      }    
    }
    ++iPad;
  }

  gStyle->SetPalette(1);
  hsigvsp->GetXaxis()->SetRangeUser(0.07,1.8);
  hsigvsp->SetStats(0);
  TCanvas* c1=new TCanvas("c1","Pion Fit");
  c1->SetLogz();
  hsigvsp->Draw("col");
  gdedxpi->Draw("Psame");
  gdedxk->Draw("Psame");
  gdedxp->Draw("Psame");

  TF1 *fPhobos=new TF1("fPhobos",BetheBlochPhobos,minpPionFit,maxpPionFit,5);
  fPhobos->SetParLimits(0,0.,1.E9);
  fPhobos->SetParLimits(1,0.,100.);
  fPhobos->SetParLimits(2,0.,1.);
  fPhobos->SetParLimits(3,0.,1.);
  fPhobos->SetParLimits(4,0.,1.);
  gdedxpi->Fit("fPhobos","R");

  TCanvas* c1g=new TCanvas("c1g","Global Fit");
  TF1 *fPhobosbg=new TF1("fPhobosbg",BetheBlochPhobosBG,minbgGlobFit,maxbgGlobFit,5);
  fPhobosbg->SetParLimits(0,0.,1.E9);
  fPhobosbg->SetParLimits(1,0.,100.);
  fPhobosbg->SetParLimits(2,0.,1.);
  fPhobosbg->SetParLimits(3,0.,1.);
  fPhobosbg->SetParLimits(4,0.,1.);
  gdedxbg->Draw("AP");
  gdedxbg->GetXaxis()->SetTitle("#beta#gamma");
  gdedxbg->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
  gdedxbg->Fit("fPhobosbg","R");
  c1g->Update();
  

  AliITSPIDResponse* itsresGlobFit=new AliITSPIDResponse(kFALSE);
  AliITSPIDResponse* itsres=new AliITSPIDResponse(kFALSE);
  Double_t parameters[5],parametersGlob[5];

  printf("-------- Pion Fit Parameters -----------\n");
  for(Int_t iPar=0; iPar<5; iPar++){
    parameters[iPar]=fPhobos->GetParameter(iPar);
    printf("parameters[%d]=%g;\n",iPar,parameters[iPar]);
  }
  printf("-------- Global Fit Parameters -----------\n");
  for(Int_t iPar=0; iPar<5; iPar++){
    parametersGlob[iPar]=fPhobosbg->GetParameter(iPar);
    printf("parametersGlob[%d]=%g;\n",iPar,parametersGlob[iPar]);
  }

  itsresGlobFit->SetBetheBlochParamsITSsa(parametersGlob);
  itsres->SetBetheBlochParamsITSsa(parameters);
  
  TGraph* gPion=new TGraph(0);
  TGraph* gKaon=new TGraph(0);
  TGraph* gProton=new TGraph(0);
  TGraph* gPionGlobFit=new TGraph(0);
  TGraph* gKaonGlobFit=new TGraph(0);
  TGraph* gProtonGlobFit=new TGraph(0);

  TGraphAsymmErrors* gBBpions=new TGraphAsymmErrors(0);
  TGraphAsymmErrors* gBBkaons=new TGraphAsymmErrors(0);
  TGraphAsymmErrors* gBBprotons=new TGraphAsymmErrors(0);
  TGraphAsymmErrors* gBBpionsGlobFit=new TGraphAsymmErrors(0);
  TGraphAsymmErrors* gBBkaonsGlobFit=new TGraphAsymmErrors(0);
  TGraphAsymmErrors* gBBprotonsGlobFit=new TGraphAsymmErrors(0);

  Double_t numberofsigmapions=2.;
  Float_t resodedx[4];
  resodedx[0]=0.0001; //default value
  resodedx[1]=0.0001; //default value
  resodedx[2]=0.116;
  resodedx[3]=0.103;
  

  for(Int_t ip=0; ip<400; ip++){
    Double_t mom=0.1+0.02*ip;
    Double_t dedxpions=itsres->Bethe(mom,massPion,kTRUE);
    Double_t dedxkaons=itsres->Bethe(mom,massKaon,kTRUE);
    Double_t dedxprotons=itsres->Bethe(mom,massProton,kTRUE);
    Double_t dedxdeutons=itsres->Bethe(mom,massDeuteron,kTRUE);
    Double_t dedxpionsGlobFit=itsresGlobFit->Bethe(mom,massPion,kTRUE);
    Double_t dedxkaonsGlobFit=itsresGlobFit->Bethe(mom,massKaon,kTRUE);
    Double_t dedxprotonsGlobFit=itsresGlobFit->Bethe(mom,massProton,kTRUE);
    Double_t dedxdeutonsGlobFit=itsresGlobFit->Bethe(mom,massDeuteron,kTRUE);

    gPion->SetPoint(ip,mom,dedxpions);
    gKaon->SetPoint(ip,mom,dedxkaons);
    gProton->SetPoint(ip,mom,dedxprotons);
    gPionGlobFit->SetPoint(ip,mom,dedxpionsGlobFit);
    gKaonGlobFit->SetPoint(ip,mom,dedxkaonsGlobFit);
    gProtonGlobFit->SetPoint(ip,mom,dedxprotonsGlobFit);

    gBBpions->SetPoint(ip,mom,dedxpions);
    gBBpions->SetPointError(ip,0,0,(dedxpions*resodedx[3]*numberofsigmapions),TMath::Abs((dedxpions-dedxkaons)/2));
    gBBkaons->SetPoint(ip,mom,dedxkaons);
    gBBkaons->SetPointError(ip,0,0,TMath::Abs((dedxpions-dedxkaons)/2),TMath::Abs((dedxprotons-dedxkaons)/2));
    gBBprotons->SetPoint(ip,mom,dedxprotons);
    gBBprotons->SetPointError(ip,0,0,TMath::Abs((dedxprotons-dedxkaons)/2),TMath::Abs((dedxprotons-dedxdeutons)/2));

    gBBpionsGlobFit->SetPoint(ip,mom,dedxpionsGlobFit);
    gBBpionsGlobFit->SetPointError(ip,0,0,(dedxpionsGlobFit*resodedx[3]*numberofsigmapions),TMath::Abs((dedxpions-dedxkaonsGlobFit)/2));
    gBBkaonsGlobFit->SetPoint(ip,mom,dedxkaonsGlobFit);
    gBBkaonsGlobFit->SetPointError(ip,0,0,TMath::Abs((dedxpionsGlobFit-dedxkaonsGlobFit)/2),TMath::Abs((dedxprotonsGlobFit-dedxkaonsGlobFit)/2));
    gBBprotonsGlobFit->SetPoint(ip,mom,dedxprotonsGlobFit);
    gBBprotonsGlobFit->SetPointError(ip,0,0,TMath::Abs((dedxprotonsGlobFit-dedxkaonsGlobFit)/2),TMath::Abs((dedxprotonsGlobFit-dedxdeutonsGlobFit)/2));
  }


  TCanvas* c2=new TCanvas("c2","Parameterizations");
  c2->SetLogz();
  hsigvsp->SetMinimum(20);
  hsigvsp->Draw("col");
  gPion->Draw("LSAME");
  gKaon->Draw("LSAME");
  gProton->Draw("LSAME");
  gPionGlobFit->SetLineColor(kRed+1);
  gKaonGlobFit->SetLineColor(kRed+1);
  gProtonGlobFit->SetLineColor(kRed+1);
  gPionGlobFit->Draw("LSAME");
  gKaonGlobFit->Draw("LSAME");
  gProtonGlobFit->Draw("LSAME");
  TLatex* t1=new TLatex(0.7,0.85,"Pion Fit");
  t1->SetNDC();
  t1->Draw();
  TLatex* tG=new TLatex(0.7,0.75,"Global Fit");
  tG->SetNDC();
  tG->SetTextColor(kRed+1);
  tG->Draw();


  gBBpions->SetLineColor(2);
  gBBpions->SetLineWidth(2);
  gBBpions->SetFillColor(2);
  gBBpions->SetFillStyle(3001);
  gBBkaons->SetLineColor(3);
  gBBkaons->SetLineWidth(2);
  gBBkaons->SetFillColor(3);
  gBBkaons->SetFillStyle(3001);
  gBBprotons->SetLineColor(4);
  gBBprotons->SetLineWidth(2);
  gBBprotons->SetFillColor(4);
  gBBprotons->SetFillStyle(3001);

  TCanvas* c3=new TCanvas("c3","Asymm Bands, PionFit");
  c3->SetLogz();
  c3->SetLogx();
  hsigvsp->SetMinimum(20);
  hsigvsp->Draw("col");
  gBBpions->Draw("L4same");
  gBBkaons->Draw("L4same");
  gBBprotons->Draw("L4same");

  gBBpionsGlobFit->SetLineColor(2);
  gBBpionsGlobFit->SetLineWidth(2);
  gBBpionsGlobFit->SetFillColor(2);
  gBBpionsGlobFit->SetFillStyle(3001);
  gBBkaonsGlobFit->SetLineColor(3);
  gBBkaonsGlobFit->SetLineWidth(2);
  gBBkaonsGlobFit->SetFillColor(3);
  gBBkaonsGlobFit->SetFillStyle(3001);
  gBBprotonsGlobFit->SetLineColor(4);
  gBBprotonsGlobFit->SetLineWidth(2);
  gBBprotonsGlobFit->SetFillColor(4);
  gBBprotonsGlobFit->SetFillStyle(3001);

  TCanvas* c3g=new TCanvas("c3g","Asymm Bands, GlobFit");
  c3g->SetLogz();
  c3g->SetLogx();
  hsigvsp->SetMinimum(20);
  hsigvsp->Draw("col");
  gBBpionsGlobFit->Draw("L4same");
  gBBkaonsGlobFit->Draw("L4same");
  gBBprotonsGlobFit->Draw("L4same");

}


Double_t BetheBlochPhobos(const Double_t *x, const Double_t *par){

  //
  // returns AliExternalTrackParam::BetheBloch normalized to 
  // fgMIP at the minimum
  //

  Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
  Double_t bg=x[0]/massPion;
  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];

  Double_t bb=0.;
  if(gamma>=0. && beta>0.){
    bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
  }
  return bb;
}

Double_t BetheBlochPhobosBG(const Double_t *x, const Double_t *par){

  //
  // returns AliExternalTrackParam::BetheBloch normalized to 
  // fgMIP at the minimum
  //

  Double_t bg=x[0];
  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];

  Double_t bb=0.;
  if(gamma>=0. && beta>0.){
    bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
  }
  return bb;
}

Double_t BetheBlochAleph(const Double_t *x, const Double_t *par){

  //
  // This is the empirical ALEPH parameterization of the Bethe-Bloch formula.
  // It is normalized to 1 at the minimum.
  //
  // bg - beta*gamma
  // 
  // The default values for the kp* parameters are for ALICE TPC.
  // The returned value is in MIP units
  //

  Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
  Double_t bg=x[0]/massPion;
  Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);

  Double_t aa = TMath::Power(beta,par[3]);
  Double_t bb = TMath::Power(1./bg,par[4]);

  bb=TMath::Log(par[2]+bb);
  
  return (par[1]-aa-bb)*par[0]/aa;
}
Commit	Line	Data
fba679b6	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include <TCanvas.h>
	3	#include <TGraph.h>
	4	#include <TGraphErrors.h>
	5	#include <TGraphAsymmErrors.h>
	6	#include <TH1.h>
	7	#include <TF1.h>
	8	#include <TH2.h>
	9	#include <TFile.h>
	10	#include <TList.h>
	11	#include <TString.h>
	12	#include <TDirectoryFile.h>
	13	#include <TLatex.h>
	14	#include <TStyle.h>
	15	#include <TMath.h>
	16	#include <TDatabasePDG.h>
	17	#include "AliITSPIDResponse.h"
	18	#endif
	19
	20	/* $Id$ */
	21
	22
	23	// Macro to extract the parameterizations of dE/dx in ITS vs. p
	24	// using Phobos BB parameterization
	25
	26
	27	Double_t BetheBlochPhobos(const Double_t x, const Double_t par);
	28	Double_t BetheBlochPhobosBG(const Double_t x, const Double_t par);
	29	Double_t BetheBlochAleph(const Double_t x, const Double_t par);
	30
	31
	32	void ParametrizeITSBetheBloch(Bool_t optFromITSsa=kFALSE){
	33
	34	TFile *fil=new TFile("AnalysisResults.root");
	35	TH2F *hsigvsp=0x0;
	36	TString dfName="PWG2SpectraITSsa";
	37	TString listName="clistITSsaMult-1to-1";
	38	TString hisName="fHistDEDX";
	39
	40	if(optFromITSsa){
	41	dfName="ITSsaTracks";
	42	listName="clistITSsaTracks";
	43	hisName="hdedxvsPt4clsITSpureSA";
	44	}
	45	TDirectoryFile* df=(TDirectoryFile*)fil->Get(dfName.Data());
	46	if(!df){
	47	printf("ERROR: directory file %s not found\n",dfName.Data());
	48	return;
	49	}
	50	TList* l=(TList*)df->Get(listName.Data());
	51	if(!l){
	52	printf("ERROR: TList %s not found\n",listName.Data());
	53	return;
	54	}
	55	hsigvsp=(TH2F*)l->FindObject(hisName.Data());
	56	if(!hsigvsp){
	57	printf("ERROR: historgam %s not found\n",hisName.Data());
	58	return;
	59	}
	60	hsigvsp->GetXaxis()->SetTitle("momentum (GeV/c)");
	61	hsigvsp->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
	62
	63	Double_t minpPionFit=0.1;
	64	Double_t maxpPionFit=0.8;
65	Double_t minbgGlobFit=0.3;
66	Double_t maxbgGlobFit=2.5;
67	Double_t minbgPions=2.;
68	Double_t maxbgPions=10.;
69	Double_t minpKaons=0.2;
70	Double_t maxpKaons=0.4;
71	Double_t maxbgKaons=0.7;
72	Double_t minpProtons=0.35;
73	Double_t maxpProtons=0.57;
74	Double_t maxbgProtons=0.7;
75
76	Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
77	Double_t massKaon=TDatabasePDG::Instance()->GetParticle(321)->Mass();
78	Double_t massProton=TDatabasePDG::Instance()->GetParticle(2212)->Mass();
79	Double_t massDeuteron=1.8756;
80
81	Int_t iPad=1;
82	Int_t iCanvas=0;
83	Int_t nbinsX=hsigvsp->GetNbinsX();
84	printf("Total Number of bins =%d\n",nbinsX);
85	Int_t group=2;
86
87	Int_t totPads=(Int_t)((Float_t)nbinsX/(Float_t)group+0.5);
88	Int_t totCanvas=totPads/16;
89	TCanvas** c0=new TCanvas*[totCanvas];
90	c0[iCanvas]=new TCanvas(Form("c0_%d",iCanvas),Form("c0_%d",iCanvas));
91	c0[iCanvas]->Divide(4,4);
92	TCanvas** c0k=new TCanvas*[totCanvas];
93	c0k[iCanvas]=new TCanvas(Form("c0k_%d",iCanvas),Form("c0k_%d",iCanvas));
94	c0k[iCanvas]->Divide(4,4);
95	TCanvas** c0p=new TCanvas*[totCanvas];
96	c0p[iCanvas]=new TCanvas(Form("c0p_%d",iCanvas),Form("c0p_%d",iCanvas));
97	c0p[iCanvas]->Divide(4,4);
98
99	TGraphErrors* gdedxpi=new TGraphErrors(0);
100	TGraphErrors* gdedxk=new TGraphErrors(0);
101	TGraphErrors* gdedxp=new TGraphErrors(0);
102	TGraphErrors* gdedxbg=new TGraphErrors(0);
103	gdedxbg->SetTitle("");
104
105	Int_t nPpi=0;
106	Int_t nPk=0;
107	Int_t nPp=0;
108	Int_t nPbg=0;
109	for(Int_t iBin=1; iBin<=nbinsX; iBin+=group){
110	if(iPad>16){
111	iCanvas++;
112	c0[iCanvas]=new TCanvas(Form("c0_%d",iCanvas),Form("c0_%d",iCanvas));
113	c0[iCanvas]->Divide(4,4);
114	c0k[iCanvas]=new TCanvas(Form("c0k_%d",iCanvas),Form("c0k_%d",iCanvas));
115	c0k[iCanvas]->Divide(4,4);
116	c0p[iCanvas]=new TCanvas(Form("c0p_%d",iCanvas),Form("c0p_%d",iCanvas));
117	c0p[iCanvas]->Divide(4,4);
118	iPad=1;
119	}
120
121	Float_t sum=0;
122	Float_t sumwei=0;
123	Double_t pmed=0.;
124	for(Int_t ibx=iBin; ibx<iBin+group;ibx++){
125	for(Int_t iby=1; iby<=hsigvsp->GetNbinsY();iby++){
126	sum+=hsigvsp->GetXaxis()->GetBinCenter(ibx)*hsigvsp->GetBinContent(ibx,iby);
127	sumwei+=hsigvsp->GetBinContent(ibx,iby);
128	}
129	}
130	if(sumwei>0.) pmed=sum/sumwei;
131	if(pmed<0.1) continue;
132	Int_t iPmed=(Int_t)(pmed*1000);
133	TH1D* hProj=hsigvsp->TH2F::ProjectionY(Form("proj%d",iPmed),iBin,iBin+group-1);
134	TH1D* hKaon=(TH1D*)hProj->Clone(Form("projk%d",iPmed));
135	c0[iCanvas]->cd(iPad);
136	hProj->SetLineColor(2);
137	hProj->Draw();
138
139	Float_t peakpos=hProj->GetBinCenter(hProj->GetMaximumBin());
140	Float_t peaksig=0.1*peakpos;
141	hProj->Fit("gaus","LL","",peakpos-3peaksig,peakpos+3.peaksig);
142	TF1* fgaus=(TF1*)hProj->GetListOfFunctions()->FindObject("gaus");
143	gdedxpi->SetPoint(nPpi,sum/sumwei,fgaus->GetParameter(1));
144	gdedxpi->SetPointError(nPpi,0.,fgaus->GetParError(1));
145	nPpi++;
146	if(pmed/massPion>minbgPions && pmed/massPion<maxbgPions){
147	gdedxbg->SetPoint(nPbg,pmed/massPion,fgaus->GetParameter(1));
148	gdedxbg->SetPointError(nPbg,0.,fgaus->GetParError(1));
149	nPbg++;
150	}
151	for(Int_t iBinK=1; iBinK<=hKaon->GetNbinsX(); iBinK++){
152	Double_t cont=hKaon->GetBinContent(iBinK)-fgaus->Eval(hKaon->GetBinCenter(iBinK));
153	if(cont<0) cont=0.;
154	if(hKaon->GetBinCenter(iBinK)<fgaus->GetParameter(1)) cont=0;
155	hKaon->SetBinContent(iBinK,cont);
156	}
157	c0k[iCanvas]->cd(iPad);
158	hKaon->SetLineColor(4);
159	hKaon->Draw();
160	TH1D* hProton=(TH1D*)hKaon->Clone(Form("projp%d",iPmed));
161	TF1* fgausk=0x0;
162	Double_t mink=250-600*(pmed-0.2);
163	Double_t maxk=mink+mink;
164	hKaon->Fit("gaus","LL","",mink,maxk);
165	fgausk=(TF1*)hKaon->GetListOfFunctions()->FindObject("gaus");
166	if(pmed>minpKaons && pmed<maxpKaons){
167	gdedxk->SetPoint(nPk,pmed,fgausk->GetParameter(1));
168	gdedxk->SetPointError(nPk,0.,fgausk->GetParError(1));
169	nPk++;
170	if(pmed/massKaon<maxbgKaons){
171	gdedxbg->SetPoint(nPbg,pmed/massKaon,fgausk->GetParameter(1));
172	gdedxbg->SetPointError(nPbg,0.,fgausk->GetParError(1));
173	nPbg++;
174	}
175	}
176	if(fgausk){
177	for(Int_t iBinP=1; iBinP<=hProton->GetNbinsX(); iBinP++){
178	Double_t cont=hKaon->GetBinContent(iBinP)-fgausk->Eval(hKaon->GetBinCenter(iBinP));
179	if(cont<0) cont=0.;
180	if(hProton->GetBinCenter(iBinP)<fgausk->GetParameter(1)) cont=0;
181	hProton->SetBinContent(iBinP,cont);
182	}
183	}
184	c0p[iCanvas]->cd(iPad);
185	hProton->SetLineColor(kGreen+1);
186	hProton->Draw();
187	if(pmed>minpProtons && pmed<maxpProtons){
188	Double_t minP=300-600*(pmed-0.35);
189	Double_t maxP=minP+minP;
190	hProton->Fit("gaus","LL","",minP,maxP);
191	TF1* fgausp=(TF1*)hProton->GetListOfFunctions()->FindObject("gaus");
192	gdedxp->SetPoint(nPp,pmed,fgausp->GetParameter(1));
193	gdedxp->SetPointError(nPp,0.,fgausp->GetParError(1));
194	nPp++;
195	if(pmed/massProton<maxbgProtons){
196	gdedxbg->SetPoint(nPbg,pmed/massProton,fgausp->GetParameter(1));
197	gdedxbg->SetPointError(nPbg,0.,fgausp->GetParError(1));
198	nPbg++;
199	}
200	}
201	++iPad;
202	}
203
204	gStyle->SetPalette(1);
205	hsigvsp->GetXaxis()->SetRangeUser(0.07,1.8);
206	hsigvsp->SetStats(0);
207	TCanvas* c1=new TCanvas("c1","Pion Fit");
208	c1->SetLogz();
209	hsigvsp->Draw("col");
210	gdedxpi->Draw("Psame");
211	gdedxk->Draw("Psame");
212	gdedxp->Draw("Psame");
213
214	TF1 *fPhobos=new TF1("fPhobos",BetheBlochPhobos,minpPionFit,maxpPionFit,5);
215	fPhobos->SetParLimits(0,0.,1.E9);
216	fPhobos->SetParLimits(1,0.,100.);
217	fPhobos->SetParLimits(2,0.,1.);
218	fPhobos->SetParLimits(3,0.,1.);
219	fPhobos->SetParLimits(4,0.,1.);
220	gdedxpi->Fit("fPhobos","R");
221
222	TCanvas* c1g=new TCanvas("c1g","Global Fit");
223	TF1 *fPhobosbg=new TF1("fPhobosbg",BetheBlochPhobosBG,minbgGlobFit,maxbgGlobFit,5);
224	fPhobosbg->SetParLimits(0,0.,1.E9);
225	fPhobosbg->SetParLimits(1,0.,100.);
226	fPhobosbg->SetParLimits(2,0.,1.);
227	fPhobosbg->SetParLimits(3,0.,1.);
228	fPhobosbg->SetParLimits(4,0.,1.);
229	gdedxbg->Draw("AP");
230	gdedxbg->GetXaxis()->SetTitle("#beta#gamma");
231	gdedxbg->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
232	gdedxbg->Fit("fPhobosbg","R");
233	c1g->Update();
234
235
236	AliITSPIDResponse* itsresGlobFit=new AliITSPIDResponse(kFALSE);
237	AliITSPIDResponse* itsres=new AliITSPIDResponse(kFALSE);
238	Double_t parameters[5],parametersGlob[5];
239
240	printf("-------- Pion Fit Parameters -----------\n");
241	for(Int_t iPar=0; iPar<5; iPar++){
242	parameters[iPar]=fPhobos->GetParameter(iPar);
243	printf("parameters[%d]=%g;\n",iPar,parameters[iPar]);
244	}
245	printf("-------- Global Fit Parameters -----------\n");
246	for(Int_t iPar=0; iPar<5; iPar++){
247	parametersGlob[iPar]=fPhobosbg->GetParameter(iPar);
248	printf("parametersGlob[%d]=%g;\n",iPar,parametersGlob[iPar]);
249	}
250
251	itsresGlobFit->SetBetheBlochParamsITSsa(parametersGlob);
252	itsres->SetBetheBlochParamsITSsa(parameters);
253
254	TGraph* gPion=new TGraph(0);
255	TGraph* gKaon=new TGraph(0);
256	TGraph* gProton=new TGraph(0);
257	TGraph* gPionGlobFit=new TGraph(0);
258	TGraph* gKaonGlobFit=new TGraph(0);
259	TGraph* gProtonGlobFit=new TGraph(0);
260
261	TGraphAsymmErrors* gBBpions=new TGraphAsymmErrors(0);
262	TGraphAsymmErrors* gBBkaons=new TGraphAsymmErrors(0);
263	TGraphAsymmErrors* gBBprotons=new TGraphAsymmErrors(0);
264	TGraphAsymmErrors* gBBpionsGlobFit=new TGraphAsymmErrors(0);
265	TGraphAsymmErrors* gBBkaonsGlobFit=new TGraphAsymmErrors(0);
266	TGraphAsymmErrors* gBBprotonsGlobFit=new TGraphAsymmErrors(0);
267
268	Double_t numberofsigmapions=2.;
269	Float_t resodedx[4];
270	resodedx[0]=0.0001; //default value
271	resodedx[1]=0.0001; //default value
272	resodedx[2]=0.116;
273	resodedx[3]=0.103;
274
275
276	for(Int_t ip=0; ip<400; ip++){
277	Double_t mom=0.1+0.02*ip;
278	Double_t dedxpions=itsres->Bethe(mom,massPion,kTRUE);
279	Double_t dedxkaons=itsres->Bethe(mom,massKaon,kTRUE);
280	Double_t dedxprotons=itsres->Bethe(mom,massProton,kTRUE);
281	Double_t dedxdeutons=itsres->Bethe(mom,massDeuteron,kTRUE);
282	Double_t dedxpionsGlobFit=itsresGlobFit->Bethe(mom,massPion,kTRUE);
283	Double_t dedxkaonsGlobFit=itsresGlobFit->Bethe(mom,massKaon,kTRUE);
284	Double_t dedxprotonsGlobFit=itsresGlobFit->Bethe(mom,massProton,kTRUE);
285	Double_t dedxdeutonsGlobFit=itsresGlobFit->Bethe(mom,massDeuteron,kTRUE);
286
287	gPion->SetPoint(ip,mom,dedxpions);
288	gKaon->SetPoint(ip,mom,dedxkaons);
289	gProton->SetPoint(ip,mom,dedxprotons);
290	gPionGlobFit->SetPoint(ip,mom,dedxpionsGlobFit);
291	gKaonGlobFit->SetPoint(ip,mom,dedxkaonsGlobFit);
292	gProtonGlobFit->SetPoint(ip,mom,dedxprotonsGlobFit);
293
294	gBBpions->SetPoint(ip,mom,dedxpions);
295	gBBpions->SetPointError(ip,0,0,(dedxpionsresodedx[3]numberofsigmapions),TMath::Abs((dedxpions-dedxkaons)/2));
296	gBBkaons->SetPoint(ip,mom,dedxkaons);
297	gBBkaons->SetPointError(ip,0,0,TMath::Abs((dedxpions-dedxkaons)/2),TMath::Abs((dedxprotons-dedxkaons)/2));
298	gBBprotons->SetPoint(ip,mom,dedxprotons);
299	gBBprotons->SetPointError(ip,0,0,TMath::Abs((dedxprotons-dedxkaons)/2),TMath::Abs((dedxprotons-dedxdeutons)/2));
300
301	gBBpionsGlobFit->SetPoint(ip,mom,dedxpionsGlobFit);
302	gBBpionsGlobFit->SetPointError(ip,0,0,(dedxpionsGlobFitresodedx[3]numberofsigmapions),TMath::Abs((dedxpions-dedxkaonsGlobFit)/2));
303	gBBkaonsGlobFit->SetPoint(ip,mom,dedxkaonsGlobFit);
304	gBBkaonsGlobFit->SetPointError(ip,0,0,TMath::Abs((dedxpionsGlobFit-dedxkaonsGlobFit)/2),TMath::Abs((dedxprotonsGlobFit-dedxkaonsGlobFit)/2));
305	gBBprotonsGlobFit->SetPoint(ip,mom,dedxprotonsGlobFit);
306	gBBprotonsGlobFit->SetPointError(ip,0,0,TMath::Abs((dedxprotonsGlobFit-dedxkaonsGlobFit)/2),TMath::Abs((dedxprotonsGlobFit-dedxdeutonsGlobFit)/2));
307	}
308
309
310	TCanvas* c2=new TCanvas("c2","Parameterizations");
311	c2->SetLogz();
312	hsigvsp->SetMinimum(20);
313	hsigvsp->Draw("col");
314	gPion->Draw("LSAME");
315	gKaon->Draw("LSAME");
316	gProton->Draw("LSAME");
317	gPionGlobFit->SetLineColor(kRed+1);
318	gKaonGlobFit->SetLineColor(kRed+1);
319	gProtonGlobFit->SetLineColor(kRed+1);
320	gPionGlobFit->Draw("LSAME");
321	gKaonGlobFit->Draw("LSAME");
322	gProtonGlobFit->Draw("LSAME");
323	TLatex* t1=new TLatex(0.7,0.85,"Pion Fit");
324	t1->SetNDC();
325	t1->Draw();
326	TLatex* tG=new TLatex(0.7,0.75,"Global Fit");
327	tG->SetNDC();
328	tG->SetTextColor(kRed+1);
329	tG->Draw();
330
331
332	gBBpions->SetLineColor(2);
333	gBBpions->SetLineWidth(2);
334	gBBpions->SetFillColor(2);
335	gBBpions->SetFillStyle(3001);
336	gBBkaons->SetLineColor(3);
337	gBBkaons->SetLineWidth(2);
338	gBBkaons->SetFillColor(3);
339	gBBkaons->SetFillStyle(3001);
340	gBBprotons->SetLineColor(4);
341	gBBprotons->SetLineWidth(2);
342	gBBprotons->SetFillColor(4);
343	gBBprotons->SetFillStyle(3001);
344
345	TCanvas* c3=new TCanvas("c3","Asymm Bands, PionFit");
346	c3->SetLogz();
347	c3->SetLogx();
348	hsigvsp->SetMinimum(20);
349	hsigvsp->Draw("col");
350	gBBpions->Draw("L4same");
351	gBBkaons->Draw("L4same");
352	gBBprotons->Draw("L4same");
353
354	gBBpionsGlobFit->SetLineColor(2);
355	gBBpionsGlobFit->SetLineWidth(2);
356	gBBpionsGlobFit->SetFillColor(2);
357	gBBpionsGlobFit->SetFillStyle(3001);
358	gBBkaonsGlobFit->SetLineColor(3);
359	gBBkaonsGlobFit->SetLineWidth(2);
360	gBBkaonsGlobFit->SetFillColor(3);
361	gBBkaonsGlobFit->SetFillStyle(3001);
362	gBBprotonsGlobFit->SetLineColor(4);
363	gBBprotonsGlobFit->SetLineWidth(2);
364	gBBprotonsGlobFit->SetFillColor(4);
365	gBBprotonsGlobFit->SetFillStyle(3001);
366
367	TCanvas* c3g=new TCanvas("c3g","Asymm Bands, GlobFit");
368	c3g->SetLogz();
369	c3g->SetLogx();
370	hsigvsp->SetMinimum(20);
371	hsigvsp->Draw("col");
372	gBBpionsGlobFit->Draw("L4same");
373	gBBkaonsGlobFit->Draw("L4same");
374	gBBprotonsGlobFit->Draw("L4same");
375
376	}
377
378
379
380
381
382
383	Double_t BetheBlochPhobos(const Double_t x, const Double_t par){
384
385	//
386	// returns AliExternalTrackParam::BetheBloch normalized to
387	// fgMIP at the minimum
388	//
389
390	Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
391	Double_t bg=x[0]/massPion;
392	Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
393	Double_t gamma=bg/beta;
394	Double_t eff=1.0;
395	if(bg<par[2])
396	eff=(bg-par[3])*(bg-par[3])+par[4];
397	else
398	eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
399
400	Double_t bb=0.;
401	if(gamma>=0. && beta>0.){
402	bb=(par[1]+2.0TMath::Log(gamma)-betabeta)(par[0]/(betabeta))*eff;
403	}
404	return bb;
405	}
406
407	Double_t BetheBlochPhobosBG(const Double_t x, const Double_t par){
408
409	//
410	// returns AliExternalTrackParam::BetheBloch normalized to
411	// fgMIP at the minimum
412	//
413
414	Double_t bg=x[0];
415	Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
416	Double_t gamma=bg/beta;
417	Double_t eff=1.0;
418	if(bg<par[2])
419	eff=(bg-par[3])*(bg-par[3])+par[4];
420	else
421	eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
422
423	Double_t bb=0.;
424	if(gamma>=0. && beta>0.){
425	bb=(par[1]+2.0TMath::Log(gamma)-betabeta)(par[0]/(betabeta))*eff;
426	}
427	return bb;
428	}
429
430	Double_t BetheBlochAleph(const Double_t x, const Double_t par){
431
432	//
433	// This is the empirical ALEPH parameterization of the Bethe-Bloch formula.
434	// It is normalized to 1 at the minimum.
435	//
436	// bg - beta*gamma
437	//
438	// The default values for the kp* parameters are for ALICE TPC.
439	// The returned value is in MIP units
440	//
441
442	Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
443	Double_t bg=x[0]/massPion;
444	Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
445
446	Double_t aa = TMath::Power(beta,par[3]);
447	Double_t bb = TMath::Power(1./bg,par[4]);
448
449	bb=TMath::Log(par[2]+bb);
450
451	return (par[1]-aa-bb)*par[0]/aa;
452	}