[u/mrichter/AliRoot.git] / ITS / AliITSdEdxAnalyzer.cxx

/**************************************************************************
 * 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$ */

///////////////////////////////////////////////////////////////////
//                                                               //
// Implementation of the base class for dEdx analysis            //
// Origin: F.Prino, Torino, prino@to.infn.it                     //
//                                                               //
///////////////////////////////////////////////////////////////////

#include <TString.h>
#include <TFile.h>
#include <TParticlePDG.h>
#include <TDatabasePDG.h>
#include "AliStack.h"
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliStack.h"
#include "AliITSdEdxAnalyzer.h"
#include "AliExternalTrackParam.h"
//#include "AliITSpidESD.h"
#include "AliITSPIDResponse.h"
#include "AliPID.h"

ClassImp(AliITSdEdxAnalyzer)

const Int_t AliITSdEdxAnalyzer::fgkPdgCode[kNParticles]={211,321,2212};
const Int_t AliITSdEdxAnalyzer::fgkLayerCode[kNValuesdEdx]={3,4,5,6,0};

//______________________________________________________________________
AliITSdEdxAnalyzer::AliITSdEdxAnalyzer() :  
  TObject(),
  fNPBins(10),
  fPMin(0.1),
  fPMax(1.1),
  fHistodEdx(0),
  fHistoDeltadEdx(0),
  fHistodEdxVsP(0),
  fThickness(0.03),
  fDensity(2.33),
  fBBmodel(0),
  fMIP(82.),
  fTPCpidCut(0.5)
{
  // default constructor
  BookHistos();
}
//______________________________________________________________________
AliITSdEdxAnalyzer::AliITSdEdxAnalyzer(const Int_t npbins, const Float_t pmin, const Float_t pmax) :  
  TObject(),
  fNPBins(npbins),
  fPMin(pmin),
  fPMax(pmax),
  fHistodEdx(0),
  fHistoDeltadEdx(0),
  fHistodEdxVsP(0),
  fThickness(0.03),
  fDensity(2.33),
  fBBmodel(0),
  fMIP(82.),
  fTPCpidCut(0.5)
{
  // standard constructor
  BookHistos();
}
//______________________________________________________________________
AliITSdEdxAnalyzer::~AliITSdEdxAnalyzer(){
  // destructor
  DeleteHistos();
}
//______________________________________________________________________
void AliITSdEdxAnalyzer::SetMomentumBins(const Int_t npbins, const Float_t pmin, const Float_t pmax){
  // Kill exisiting histos, set new binning, book new histos
  DeleteHistos();
  fNPBins=npbins;
  fPMin=pmin;
  fPMax=pmax;
  BookHistos();
}
//______________________________________________________________________
void AliITSdEdxAnalyzer::DeleteHistos(){
  // deletes the hitograms
  if(fHistodEdx){
    for(Int_t i=0; i<GetNHistos();i++) delete fHistodEdx[i];
    delete [] fHistodEdx;
    fHistodEdx=0;
  }
  if(fHistoDeltadEdx){
    for(Int_t i=0; i<GetNHistos();i++) delete fHistoDeltadEdx[i];
    delete [] fHistoDeltadEdx;
    fHistoDeltadEdx=0;
  }
  if(fHistodEdxVsP){
    for(Int_t i=0; i<GetNHistos2();i++) delete fHistodEdxVsP[i];
    delete [] fHistodEdxVsP;
    fHistodEdxVsP=0;
  }
}
//______________________________________________________________________
void AliITSdEdxAnalyzer::BookHistos(){
  // defines the output histograms 
  fHistodEdx=new TH1F*[GetNHistos()];
  fHistoDeltadEdx=new TH1F*[GetNHistos()];
  for(Int_t iSpecie=0; iSpecie<kNParticles; iSpecie++){
    for(Int_t iPBin=0; iPBin<fNPBins; iPBin++){
      for(Int_t iVal=0; iVal<kNValuesdEdx; iVal++){
	TString hisnam1=Form("hdEdx%dpbin%dpart%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
	TString histit1=Form("dEdx layer %d (keV) pbin%d part%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
	if(iVal==kNValuesdEdx-1) histit1=Form("dEdx trunc. mean (keV) pbin%d part%d",iPBin,fgkPdgCode[iSpecie]);
	fHistodEdx[GetIndex(iSpecie,iPBin,iVal)]=new TH1F(hisnam1.Data(),histit1.Data(),100,0.,500.);
	fHistodEdx[GetIndex(iSpecie,iPBin,iVal)]->GetXaxis()->SetTitle(histit1.Data());
	TString hisnam2=Form("hdeltadEdx%dpbin%dpart%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
	TString histit2=Form("(dEdx_l%d-BetheBloch)/BetheBloch pbin%d part%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
	if(iVal==kNValuesdEdx-1) histit1=Form("(dEdx_truncmean-BetheBloch)/BetheBloch pbin%d part%d",iPBin,fgkPdgCode[iSpecie]);
	fHistoDeltadEdx[GetIndex(iSpecie,iPBin,iVal)]=new TH1F(hisnam2.Data(),histit2.Data(),100,-0.5,0.5);
	fHistoDeltadEdx[GetIndex(iSpecie,iPBin,iVal)]->GetXaxis()->SetTitle(histit2.Data());
      }
    }
  }

  fHistodEdxVsP=new TH2F*[GetNHistos2()];
  for(Int_t iSpecie=0; iSpecie<kNParticles; iSpecie++){
    for(Int_t iVal=0; iVal<kNValuesdEdx; iVal++){
      TString hisnam=Form("hdEdx%dVsPpart%d",fgkLayerCode[iVal],fgkPdgCode[iSpecie]);
      TString histit=Form("dEdx layer %d (keV) vs P part%d",fgkLayerCode[iVal],fgkPdgCode[iSpecie]);
      if(iVal==kNValuesdEdx-1) histit=Form("dEdx trunc. mean (keV) vs P part%d",fgkPdgCode[iSpecie]);      
      fHistodEdxVsP[GetIndex2(iSpecie,iVal)]=new TH2F(hisnam.Data(),histit.Data(),50,fPMin,fPMax,50,0.,500.);
      histit.ReplaceAll(" vs P "," ");
      fHistodEdxVsP[GetIndex2(iSpecie,iVal)]->GetYaxis()->SetTitle(histit.Data());
      fHistodEdxVsP[GetIndex2(iSpecie,iVal)]->GetXaxis()->SetTitle("Momentum (GeV/c)");
    }
  }
}
//______________________________________________________________________
void AliITSdEdxAnalyzer:: WriteHistos(TString filename) const {
  // stores output histogrma to file
  TFile *outfile=new TFile(filename.Data(),"recreate");
  outfile->cd();
  for(Int_t i=0; i<GetNHistos();i++){ 
    fHistodEdx[i]->Write();
    fHistoDeltadEdx[i]->Write();
  }
  for(Int_t i=0; i<GetNHistos2();i++) fHistodEdxVsP[i]->Write();
  outfile->Close();
  delete outfile;
}
//______________________________________________________________________
void AliITSdEdxAnalyzer::ReadEvent(const AliESDEvent* esd, AliStack* stack){
  // Fill histos 
  // if stack is !=0 MC truth is used to define particle specie
  // else TPC pid is used to define particle specie

  if(!esd) AliFatal("Bad ESD event");
  for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
    AliESDtrack* track = esd->GetTrack(iTrack);
    Double_t trmean=track->GetITSsignal();
    Double_t sig[4];
    track->GetITSdEdxSamples(sig);
    Double_t preco=track->P();
    Int_t iPBin=GetMomentumBin(preco);
    if(iPBin==-1) continue;
    Int_t iSpecie=-1;
    Double_t dedxbb=0;
    if(stack){
      Int_t lab=track->GetLabel();
      if(lab<0) continue;
      TParticle* part=(TParticle*)stack->Particle(lab);
      Int_t absPdgCode=TMath::Abs(part->GetPdgCode());
      iSpecie=GetSpecieBin(absPdgCode);
      if(iSpecie==-1) continue;
      dedxbb=BetheBloch(part);
    }else{   
      iSpecie=GetPaticleIdFromTPC(track);
      if(iSpecie==-1) continue;
      TParticlePDG* p=TDatabasePDG::Instance()->GetParticle(fgkPdgCode[iSpecie]);    
      dedxbb=BetheBloch(preco,p->Mass());
    }
    for(Int_t ilay=0;ilay<4;ilay++){    
      if(sig[ilay]>0.){
	fHistodEdx[GetIndex(iSpecie,iPBin,ilay)]->Fill(sig[ilay]);
	fHistoDeltadEdx[GetIndex(iSpecie,iPBin,ilay)]->Fill((sig[ilay]-dedxbb)/dedxbb);
	fHistodEdxVsP[GetIndex2(iSpecie,ilay)]->Fill(preco,sig[ilay]);
      }
    }
    if(trmean>0.){
      fHistodEdx[GetIndex(iSpecie,iPBin,4)]->Fill(trmean);
      fHistoDeltadEdx[GetIndex(iSpecie,iPBin,4)]->Fill((trmean-dedxbb)/dedxbb);
      fHistodEdxVsP[GetIndex2(iSpecie,4)]->Fill(preco,trmean);
    }
  }
}
//______________________________________________________________________
Int_t AliITSdEdxAnalyzer::GetPaticleIdFromTPC(const AliESDtrack* track) const {
  // Returns the particle specie as identified by TPC
  Double_t tpcpid[AliPID::kSPECIES];
  track->GetTPCpid(tpcpid);
  Int_t maxPos=-1;
  Float_t maxProb=0.;
  for(Int_t is=0; is<AliPID::kSPECIES; is++){
    if(tpcpid[is]>maxProb){
      maxProb=tpcpid[is];
      maxPos=is;
    }
  }
  Int_t iSpecie=-1;
  if(maxProb>fTPCpidCut){
    if(maxPos==AliPID::kPion) iSpecie=0;
    else if(maxPos==AliPID::kKaon) iSpecie=1;
    else if(maxPos==AliPID::kProton) iSpecie=2;    
  }
  return iSpecie;
}
//______________________________________________________________________
Double_t AliITSdEdxAnalyzer::BetheBloch(const Float_t p, const Float_t m) const {
  // Computes expected dE/dx value for given particle specie and momentum
  static AliITSPIDResponse pidResponse;
  Double_t dedxbb=0.;
  if(fBBmodel==0){
    Double_t betagamma=p/m;
    Double_t conv=fDensity*1E6*fThickness/116.24*fMIP;
    dedxbb=conv*AliExternalTrackParam::BetheBlochSolid(betagamma);
  }else if(fBBmodel==1){
    dedxbb=pidResponse.Bethe(p,m);
  }
  return dedxbb;
}
//______________________________________________________________________
TGraph* AliITSdEdxAnalyzer::GetBetheBlochGraph(const Int_t pdgcode) const {
 // Fills a TGraph with Bethe-Bloch expe
 TGraph* g=new TGraph(0);
 TParticlePDG* part=TDatabasePDG::Instance()->GetParticle(pdgcode);
 Float_t mass=part->Mass();
 for(Int_t ip=0; ip<100;ip++){
   Float_t p=fPMin+(ip+0.5)*(fPMax-fPMin)/100.;
   g->SetPoint(ip,p,BetheBloch(p,mass));
 }
 g->GetXaxis()->SetTitle("Momentum (GeV/c)");
 g->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
 return g;
}
Commit	Line	Data
75c74fab	1	/**************************************************************************
	2	* Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	///////////////////////////////////////////////////////////////////
	19	// //
	20	// Implementation of the base class for dEdx analysis //
	21	// Origin: F.Prino, Torino, prino@to.infn.it //
	22	// //
	23	///////////////////////////////////////////////////////////////////
	24
	25	#include <TString.h>
	26	#include <TFile.h>
	27	#include <TParticlePDG.h>
	28	#include <TDatabasePDG.h>
b34fed4e	29	#include "AliStack.h"
b34fed4e	30	#include "AliLog.h"
27199d86	31	#include "AliESDEvent.h"
27199d86	32	#include "AliStack.h"
75c74fab	33	#include "AliITSdEdxAnalyzer.h"
75c74fab	34	#include "AliExternalTrackParam.h"
10d100d4	35	//#include "AliITSpidESD.h"
10d100d4	36	#include "AliITSPIDResponse.h"
75c74fab	37	#include "AliPID.h"
	38
	39	ClassImp(AliITSdEdxAnalyzer)
	40
	41	const Int_t AliITSdEdxAnalyzer::fgkPdgCode[kNParticles]={211,321,2212};
	42	const Int_t AliITSdEdxAnalyzer::fgkLayerCode[kNValuesdEdx]={3,4,5,6,0};
	43
	44	//______________________________________________________________________
	45	AliITSdEdxAnalyzer::AliITSdEdxAnalyzer() :
	46	TObject(),
	47	fNPBins(10),
	48	fPMin(0.1),
	49	fPMax(1.1),
	50	fHistodEdx(0),
	51	fHistoDeltadEdx(0),
	52	fHistodEdxVsP(0),
	53	fThickness(0.03),
	54	fDensity(2.33),
	55	fBBmodel(0),
	56	fMIP(82.),
	57	fTPCpidCut(0.5)
	58	{
	59	// default constructor
	60	BookHistos();
	61	}
	62	//______________________________________________________________________
	63	AliITSdEdxAnalyzer::AliITSdEdxAnalyzer(const Int_t npbins, const Float_t pmin, const Float_t pmax) :
	64	TObject(),
	65	fNPBins(npbins),
	66	fPMin(pmin),
	67	fPMax(pmax),
	68	fHistodEdx(0),
	69	fHistoDeltadEdx(0),
	70	fHistodEdxVsP(0),
	71	fThickness(0.03),
	72	fDensity(2.33),
	73	fBBmodel(0),
	74	fMIP(82.),
	75	fTPCpidCut(0.5)
	76	{
	77	// standard constructor
	78	BookHistos();
	79	}
	80	//______________________________________________________________________
	81	AliITSdEdxAnalyzer::~AliITSdEdxAnalyzer(){
	82	// destructor
	83	DeleteHistos();
	84	}
	85	//______________________________________________________________________
	86	void AliITSdEdxAnalyzer::SetMomentumBins(const Int_t npbins, const Float_t pmin, const Float_t pmax){
	87	// Kill exisiting histos, set new binning, book new histos
	88	DeleteHistos();
	89	fNPBins=npbins;
	90	fPMin=pmin;
	91	fPMax=pmax;
	92	BookHistos();
	93	}
	94	//______________________________________________________________________
	95	void AliITSdEdxAnalyzer::DeleteHistos(){
b34fed4e	96	// deletes the hitograms
75c74fab	97	if(fHistodEdx){
	98	for(Int_t i=0; i<GetNHistos();i++) delete fHistodEdx[i];
	99	delete [] fHistodEdx;
	100	fHistodEdx=0;
	101	}
	102	if(fHistoDeltadEdx){
	103	for(Int_t i=0; i<GetNHistos();i++) delete fHistoDeltadEdx[i];
	104	delete [] fHistoDeltadEdx;
	105	fHistoDeltadEdx=0;
	106	}
	107	if(fHistodEdxVsP){
	108	for(Int_t i=0; i<GetNHistos2();i++) delete fHistodEdxVsP[i];
	109	delete [] fHistodEdxVsP;
	110	fHistodEdxVsP=0;
	111	}
	112	}
	113	//______________________________________________________________________
	114	void AliITSdEdxAnalyzer::BookHistos(){
b34fed4e	115	// defines the output histograms
75c74fab	116	fHistodEdx=new TH1F*[GetNHistos()];
	117	fHistoDeltadEdx=new TH1F*[GetNHistos()];
	118	for(Int_t iSpecie=0; iSpecie<kNParticles; iSpecie++){
	119	for(Int_t iPBin=0; iPBin<fNPBins; iPBin++){
	120	for(Int_t iVal=0; iVal<kNValuesdEdx; iVal++){
	121	TString hisnam1=Form("hdEdx%dpbin%dpart%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
	122	TString histit1=Form("dEdx layer %d (keV) pbin%d part%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
	123	if(iVal==kNValuesdEdx-1) histit1=Form("dEdx trunc. mean (keV) pbin%d part%d",iPBin,fgkPdgCode[iSpecie]);
	124	fHistodEdx[GetIndex(iSpecie,iPBin,iVal)]=new TH1F(hisnam1.Data(),histit1.Data(),100,0.,500.);
	125	fHistodEdx[GetIndex(iSpecie,iPBin,iVal)]->GetXaxis()->SetTitle(histit1.Data());
	126	TString hisnam2=Form("hdeltadEdx%dpbin%dpart%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
	127	TString histit2=Form("(dEdx_l%d-BetheBloch)/BetheBloch pbin%d part%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
	128	if(iVal==kNValuesdEdx-1) histit1=Form("(dEdx_truncmean-BetheBloch)/BetheBloch pbin%d part%d",iPBin,fgkPdgCode[iSpecie]);
	129	fHistoDeltadEdx[GetIndex(iSpecie,iPBin,iVal)]=new TH1F(hisnam2.Data(),histit2.Data(),100,-0.5,0.5);
	130	fHistoDeltadEdx[GetIndex(iSpecie,iPBin,iVal)]->GetXaxis()->SetTitle(histit2.Data());
	131	}
	132	}
	133	}
	134
	135	fHistodEdxVsP=new TH2F*[GetNHistos2()];
	136	for(Int_t iSpecie=0; iSpecie<kNParticles; iSpecie++){
	137	for(Int_t iVal=0; iVal<kNValuesdEdx; iVal++){
	138	TString hisnam=Form("hdEdx%dVsPpart%d",fgkLayerCode[iVal],fgkPdgCode[iSpecie]);
	139	TString histit=Form("dEdx layer %d (keV) vs P part%d",fgkLayerCode[iVal],fgkPdgCode[iSpecie]);
	140	if(iVal==kNValuesdEdx-1) histit=Form("dEdx trunc. mean (keV) vs P part%d",fgkPdgCode[iSpecie]);
	141	fHistodEdxVsP[GetIndex2(iSpecie,iVal)]=new TH2F(hisnam.Data(),histit.Data(),50,fPMin,fPMax,50,0.,500.);
	142	histit.ReplaceAll(" vs P "," ");
	143	fHistodEdxVsP[GetIndex2(iSpecie,iVal)]->GetYaxis()->SetTitle(histit.Data());
	144	fHistodEdxVsP[GetIndex2(iSpecie,iVal)]->GetXaxis()->SetTitle("Momentum (GeV/c)");
	145	}
	146	}
	147	}
	148	//______________________________________________________________________
	149	void AliITSdEdxAnalyzer:: WriteHistos(TString filename) const {
b34fed4e	150	// stores output histogrma to file
75c74fab	151	TFile *outfile=new TFile(filename.Data(),"recreate");
	152	outfile->cd();
	153	for(Int_t i=0; i<GetNHistos();i++){
	154	fHistodEdx[i]->Write();
	155	fHistoDeltadEdx[i]->Write();
	156	}
	157	for(Int_t i=0; i<GetNHistos2();i++) fHistodEdxVsP[i]->Write();
	158	outfile->Close();
	159	delete outfile;
	160	}
	161	//______________________________________________________________________
b34fed4e	162	void AliITSdEdxAnalyzer::ReadEvent(const AliESDEvent* esd, AliStack* stack){
75c74fab	163	// Fill histos
	164	// if stack is !=0 MC truth is used to define particle specie
	165	// else TPC pid is used to define particle specie
	166
	167	if(!esd) AliFatal("Bad ESD event");
	168	for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
	169	AliESDtrack* track = esd->GetTrack(iTrack);
	170	Double_t trmean=track->GetITSsignal();
	171	Double_t sig[4];
	172	track->GetITSdEdxSamples(sig);
	173	Double_t preco=track->P();
	174	Int_t iPBin=GetMomentumBin(preco);
	175	if(iPBin==-1) continue;
	176	Int_t iSpecie=-1;
	177	Double_t dedxbb=0;
	178	if(stack){
	179	Int_t lab=track->GetLabel();
	180	if(lab<0) continue;
	181	TParticle* part=(TParticle*)stack->Particle(lab);
	182	Int_t absPdgCode=TMath::Abs(part->GetPdgCode());
	183	iSpecie=GetSpecieBin(absPdgCode);
	184	if(iSpecie==-1) continue;
	185	dedxbb=BetheBloch(part);
	186	}else{
	187	iSpecie=GetPaticleIdFromTPC(track);
	188	if(iSpecie==-1) continue;
	189	TParticlePDG* p=TDatabasePDG::Instance()->GetParticle(fgkPdgCode[iSpecie]);
	190	dedxbb=BetheBloch(preco,p->Mass());
	191	}
	192	for(Int_t ilay=0;ilay<4;ilay++){
	193	if(sig[ilay]>0.){
	194	fHistodEdx[GetIndex(iSpecie,iPBin,ilay)]->Fill(sig[ilay]);
	195	fHistoDeltadEdx[GetIndex(iSpecie,iPBin,ilay)]->Fill((sig[ilay]-dedxbb)/dedxbb);
	196	fHistodEdxVsP[GetIndex2(iSpecie,ilay)]->Fill(preco,sig[ilay]);
	197	}
	198	}
	199	if(trmean>0.){
	200	fHistodEdx[GetIndex(iSpecie,iPBin,4)]->Fill(trmean);
	201	fHistoDeltadEdx[GetIndex(iSpecie,iPBin,4)]->Fill((trmean-dedxbb)/dedxbb);
	202	fHistodEdxVsP[GetIndex2(iSpecie,4)]->Fill(preco,trmean);
	203	}
	204	}
	205	}
	206	//______________________________________________________________________
	207	Int_t AliITSdEdxAnalyzer::GetPaticleIdFromTPC(const AliESDtrack* track) const {
b34fed4e	208	// Returns the particle specie as identified by TPC
75c74fab	209	Double_t tpcpid[AliPID::kSPECIES];
	210	track->GetTPCpid(tpcpid);
	211	Int_t maxPos=-1;
	212	Float_t maxProb=0.;
	213	for(Int_t is=0; is<AliPID::kSPECIES; is++){
	214	if(tpcpid[is]>maxProb){
	215	maxProb=tpcpid[is];
	216	maxPos=is;
	217	}
	218	}
	219	Int_t iSpecie=-1;
	220	if(maxProb>fTPCpidCut){
	221	if(maxPos==AliPID::kPion) iSpecie=0;
	222	else if(maxPos==AliPID::kKaon) iSpecie=1;
	223	else if(maxPos==AliPID::kProton) iSpecie=2;
	224	}
	225	return iSpecie;
	226	}
	227	//______________________________________________________________________
	228	Double_t AliITSdEdxAnalyzer::BetheBloch(const Float_t p, const Float_t m) const {
b34fed4e	229	// Computes expected dE/dx value for given particle specie and momentum
10d100d4	230	static AliITSPIDResponse pidResponse;
75c74fab	231	Double_t dedxbb=0.;
26964b7f	232	if(fBBmodel==0){
75c74fab	233	Double_t betagamma=p/m;
b89be497	234	Double_t conv=fDensity1E6fThickness/116.24*fMIP;
75c74fab	235	dedxbb=conv*AliExternalTrackParam::BetheBlochSolid(betagamma);
26964b7f	236	}else if(fBBmodel==1){
10d100d4	237	dedxbb=pidResponse.Bethe(p,m);
75c74fab	238	}
	239	return dedxbb;
	240	}
	241	//______________________________________________________________________
	242	TGraph* AliITSdEdxAnalyzer::GetBetheBlochGraph(const Int_t pdgcode) const {
b246a937	243	// Fills a TGraph with Bethe-Bloch expe
	244	TGraph* g=new TGraph(0);
	245	TParticlePDG* part=TDatabasePDG::Instance()->GetParticle(pdgcode);
	246	Float_t mass=part->Mass();
	247	for(Int_t ip=0; ip<100;ip++){
	248	Float_t p=fPMin+(ip+0.5)*(fPMax-fPMin)/100.;
	249	g->SetPoint(ip,p,BetheBloch(p,mass));
	250	}
	251	g->GetXaxis()->SetTitle("Momentum (GeV/c)");
	252	g->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
	253	return g;
75c74fab	254	}