1 /**************************************************************************
2 * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ///////////////////////////////////////////////////////////////////
20 // Implementation of the base class for dEdx analysis //
21 // Origin: F.Prino, Torino, prino@to.infn.it //
23 ///////////////////////////////////////////////////////////////////
27 #include <TParticlePDG.h>
28 #include <TDatabasePDG.h>
29 #include "AliITSdEdxAnalyzer.h"
30 #include "AliExternalTrackParam.h"
31 #include "AliITSpidESD.h"
34 ClassImp(AliITSdEdxAnalyzer)
36 const Int_t AliITSdEdxAnalyzer::fgkPdgCode[kNParticles]={211,321,2212};
37 const Int_t AliITSdEdxAnalyzer::fgkLayerCode[kNValuesdEdx]={3,4,5,6,0};
39 //______________________________________________________________________
40 AliITSdEdxAnalyzer::AliITSdEdxAnalyzer() :
54 // default constructor
57 //______________________________________________________________________
58 AliITSdEdxAnalyzer::AliITSdEdxAnalyzer(const Int_t npbins, const Float_t pmin, const Float_t pmax) :
72 // standard constructor
75 //______________________________________________________________________
76 AliITSdEdxAnalyzer::~AliITSdEdxAnalyzer(){
80 //______________________________________________________________________
81 void AliITSdEdxAnalyzer::SetMomentumBins(const Int_t npbins, const Float_t pmin, const Float_t pmax){
82 // Kill exisiting histos, set new binning, book new histos
89 //______________________________________________________________________
90 void AliITSdEdxAnalyzer::DeleteHistos(){
93 for(Int_t i=0; i<GetNHistos();i++) delete fHistodEdx[i];
98 for(Int_t i=0; i<GetNHistos();i++) delete fHistoDeltadEdx[i];
99 delete [] fHistoDeltadEdx;
103 for(Int_t i=0; i<GetNHistos2();i++) delete fHistodEdxVsP[i];
104 delete [] fHistodEdxVsP;
108 //______________________________________________________________________
109 void AliITSdEdxAnalyzer::BookHistos(){
111 fHistodEdx=new TH1F*[GetNHistos()];
112 fHistoDeltadEdx=new TH1F*[GetNHistos()];
113 for(Int_t iSpecie=0; iSpecie<kNParticles; iSpecie++){
114 for(Int_t iPBin=0; iPBin<fNPBins; iPBin++){
115 for(Int_t iVal=0; iVal<kNValuesdEdx; iVal++){
116 TString hisnam1=Form("hdEdx%dpbin%dpart%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
117 TString histit1=Form("dEdx layer %d (keV) pbin%d part%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
118 if(iVal==kNValuesdEdx-1) histit1=Form("dEdx trunc. mean (keV) pbin%d part%d",iPBin,fgkPdgCode[iSpecie]);
119 fHistodEdx[GetIndex(iSpecie,iPBin,iVal)]=new TH1F(hisnam1.Data(),histit1.Data(),100,0.,500.);
120 fHistodEdx[GetIndex(iSpecie,iPBin,iVal)]->GetXaxis()->SetTitle(histit1.Data());
121 TString hisnam2=Form("hdeltadEdx%dpbin%dpart%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
122 TString histit2=Form("(dEdx_l%d-BetheBloch)/BetheBloch pbin%d part%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
123 if(iVal==kNValuesdEdx-1) histit1=Form("(dEdx_truncmean-BetheBloch)/BetheBloch pbin%d part%d",iPBin,fgkPdgCode[iSpecie]);
124 fHistoDeltadEdx[GetIndex(iSpecie,iPBin,iVal)]=new TH1F(hisnam2.Data(),histit2.Data(),100,-0.5,0.5);
125 fHistoDeltadEdx[GetIndex(iSpecie,iPBin,iVal)]->GetXaxis()->SetTitle(histit2.Data());
130 fHistodEdxVsP=new TH2F*[GetNHistos2()];
131 for(Int_t iSpecie=0; iSpecie<kNParticles; iSpecie++){
132 for(Int_t iVal=0; iVal<kNValuesdEdx; iVal++){
133 TString hisnam=Form("hdEdx%dVsPpart%d",fgkLayerCode[iVal],fgkPdgCode[iSpecie]);
134 TString histit=Form("dEdx layer %d (keV) vs P part%d",fgkLayerCode[iVal],fgkPdgCode[iSpecie]);
135 if(iVal==kNValuesdEdx-1) histit=Form("dEdx trunc. mean (keV) vs P part%d",fgkPdgCode[iSpecie]);
136 fHistodEdxVsP[GetIndex2(iSpecie,iVal)]=new TH2F(hisnam.Data(),histit.Data(),50,fPMin,fPMax,50,0.,500.);
137 histit.ReplaceAll(" vs P "," ");
138 fHistodEdxVsP[GetIndex2(iSpecie,iVal)]->GetYaxis()->SetTitle(histit.Data());
139 fHistodEdxVsP[GetIndex2(iSpecie,iVal)]->GetXaxis()->SetTitle("Momentum (GeV/c)");
143 //______________________________________________________________________
144 void AliITSdEdxAnalyzer:: WriteHistos(TString filename) const {
145 TFile *outfile=new TFile(filename.Data(),"recreate");
147 for(Int_t i=0; i<GetNHistos();i++){
148 fHistodEdx[i]->Write();
149 fHistoDeltadEdx[i]->Write();
151 for(Int_t i=0; i<GetNHistos2();i++) fHistodEdxVsP[i]->Write();
155 //______________________________________________________________________
156 void AliITSdEdxAnalyzer::ReadEvent(AliESDEvent* esd, AliStack* stack){
158 // if stack is !=0 MC truth is used to define particle specie
159 // else TPC pid is used to define particle specie
161 if(!esd) AliFatal("Bad ESD event");
162 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
163 AliESDtrack* track = esd->GetTrack(iTrack);
164 Double_t trmean=track->GetITSsignal();
166 track->GetITSdEdxSamples(sig);
167 Double_t preco=track->P();
168 Int_t iPBin=GetMomentumBin(preco);
169 if(iPBin==-1) continue;
173 Int_t lab=track->GetLabel();
175 TParticle* part=(TParticle*)stack->Particle(lab);
176 Int_t absPdgCode=TMath::Abs(part->GetPdgCode());
177 iSpecie=GetSpecieBin(absPdgCode);
178 if(iSpecie==-1) continue;
179 dedxbb=BetheBloch(part);
181 iSpecie=GetPaticleIdFromTPC(track);
182 if(iSpecie==-1) continue;
183 TParticlePDG* p=TDatabasePDG::Instance()->GetParticle(fgkPdgCode[iSpecie]);
184 dedxbb=BetheBloch(preco,p->Mass());
186 for(Int_t ilay=0;ilay<4;ilay++){
188 fHistodEdx[GetIndex(iSpecie,iPBin,ilay)]->Fill(sig[ilay]);
189 fHistoDeltadEdx[GetIndex(iSpecie,iPBin,ilay)]->Fill((sig[ilay]-dedxbb)/dedxbb);
190 fHistodEdxVsP[GetIndex2(iSpecie,ilay)]->Fill(preco,sig[ilay]);
194 fHistodEdx[GetIndex(iSpecie,iPBin,4)]->Fill(trmean);
195 fHistoDeltadEdx[GetIndex(iSpecie,iPBin,4)]->Fill((trmean-dedxbb)/dedxbb);
196 fHistodEdxVsP[GetIndex2(iSpecie,4)]->Fill(preco,trmean);
200 //______________________________________________________________________
201 Int_t AliITSdEdxAnalyzer::GetPaticleIdFromTPC(const AliESDtrack* track) const {
203 Double_t tpcpid[AliPID::kSPECIES];
204 track->GetTPCpid(tpcpid);
207 for(Int_t is=0; is<AliPID::kSPECIES; is++){
208 if(tpcpid[is]>maxProb){
214 if(maxProb>fTPCpidCut){
215 if(maxPos==AliPID::kPion) iSpecie=0;
216 else if(maxPos==AliPID::kKaon) iSpecie=1;
217 else if(maxPos==AliPID::kProton) iSpecie=2;
221 //______________________________________________________________________
222 Double_t AliITSdEdxAnalyzer::BetheBloch(const Float_t p, const Float_t m) const {
226 Double_t betagamma=p/m;
227 Double_t conv=fDensity*1E6*fThickness/116.31*fMIP;
228 dedxbb=conv*AliExternalTrackParam::BetheBlochSolid(betagamma);
229 }else if(fBBmodel==1){
230 dedxbb=fMIP*AliITSpidESD::Bethe(p,m);
234 //______________________________________________________________________
235 TGraph* AliITSdEdxAnalyzer::GetBetheBlochGraph(const Int_t pdgcode) const {
236 // Fills a TGraph with Bethe-Bloch expe
237 TGraph* g=new TGraph(0);
238 TParticlePDG* p=TDatabasePDG::Instance()->GetParticle(pdgcode);
239 Float_t mass=p->Mass();
240 for(Int_t ip=0; ip<100;ip++){
241 Float_t p=fPMin+(ip+0.5)*(fPMax-fPMin)/100.;
242 g->SetPoint(ip,p,BetheBloch(p,mass));
244 g->GetXaxis()->SetTitle("Momentum (GeV/c)");
245 g->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");