Coding conventions
[u/mrichter/AliRoot.git] / ITS / AliITSdEdxAnalyzer.cxx
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>
29 #include "AliStack.h"
30 #include "AliLog.h"
31 #include "AliITSdEdxAnalyzer.h"
32 #include "AliExternalTrackParam.h"
33 //#include "AliITSpidESD.h"
34 #include "AliITSPIDResponse.h"
35 #include "AliPID.h"
36
37 ClassImp(AliITSdEdxAnalyzer)
38
39 const Int_t AliITSdEdxAnalyzer::fgkPdgCode[kNParticles]={211,321,2212};
40 const Int_t AliITSdEdxAnalyzer::fgkLayerCode[kNValuesdEdx]={3,4,5,6,0};
41
42 //______________________________________________________________________
43 AliITSdEdxAnalyzer::AliITSdEdxAnalyzer() :  
44   TObject(),
45   fNPBins(10),
46   fPMin(0.1),
47   fPMax(1.1),
48   fHistodEdx(0),
49   fHistoDeltadEdx(0),
50   fHistodEdxVsP(0),
51   fThickness(0.03),
52   fDensity(2.33),
53   fBBmodel(0),
54   fMIP(82.),
55   fTPCpidCut(0.5)
56 {
57   // default constructor
58   BookHistos();
59 }
60 //______________________________________________________________________
61 AliITSdEdxAnalyzer::AliITSdEdxAnalyzer(const Int_t npbins, const Float_t pmin, const Float_t pmax) :  
62   TObject(),
63   fNPBins(npbins),
64   fPMin(pmin),
65   fPMax(pmax),
66   fHistodEdx(0),
67   fHistoDeltadEdx(0),
68   fHistodEdxVsP(0),
69   fThickness(0.03),
70   fDensity(2.33),
71   fBBmodel(0),
72   fMIP(82.),
73   fTPCpidCut(0.5)
74 {
75   // standard constructor
76   BookHistos();
77 }
78 //______________________________________________________________________
79 AliITSdEdxAnalyzer::~AliITSdEdxAnalyzer(){
80   // destructor
81   DeleteHistos();
82 }
83 //______________________________________________________________________
84 void AliITSdEdxAnalyzer::SetMomentumBins(const Int_t npbins, const Float_t pmin, const Float_t pmax){
85   // Kill exisiting histos, set new binning, book new histos
86   DeleteHistos();
87   fNPBins=npbins;
88   fPMin=pmin;
89   fPMax=pmax;
90   BookHistos();
91 }
92 //______________________________________________________________________
93 void AliITSdEdxAnalyzer::DeleteHistos(){
94   // deletes the hitograms
95   if(fHistodEdx){
96     for(Int_t i=0; i<GetNHistos();i++) delete fHistodEdx[i];
97     delete [] fHistodEdx;
98     fHistodEdx=0;
99   }
100   if(fHistoDeltadEdx){
101     for(Int_t i=0; i<GetNHistos();i++) delete fHistoDeltadEdx[i];
102     delete [] fHistoDeltadEdx;
103     fHistoDeltadEdx=0;
104   }
105   if(fHistodEdxVsP){
106     for(Int_t i=0; i<GetNHistos2();i++) delete fHistodEdxVsP[i];
107     delete [] fHistodEdxVsP;
108     fHistodEdxVsP=0;
109   }
110 }
111 //______________________________________________________________________
112 void AliITSdEdxAnalyzer::BookHistos(){
113   // defines the output histograms 
114   fHistodEdx=new TH1F*[GetNHistos()];
115   fHistoDeltadEdx=new TH1F*[GetNHistos()];
116   for(Int_t iSpecie=0; iSpecie<kNParticles; iSpecie++){
117     for(Int_t iPBin=0; iPBin<fNPBins; iPBin++){
118       for(Int_t iVal=0; iVal<kNValuesdEdx; iVal++){
119         TString hisnam1=Form("hdEdx%dpbin%dpart%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
120         TString histit1=Form("dEdx layer %d (keV) pbin%d part%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
121         if(iVal==kNValuesdEdx-1) histit1=Form("dEdx trunc. mean (keV) pbin%d part%d",iPBin,fgkPdgCode[iSpecie]);
122         fHistodEdx[GetIndex(iSpecie,iPBin,iVal)]=new TH1F(hisnam1.Data(),histit1.Data(),100,0.,500.);
123         fHistodEdx[GetIndex(iSpecie,iPBin,iVal)]->GetXaxis()->SetTitle(histit1.Data());
124         TString hisnam2=Form("hdeltadEdx%dpbin%dpart%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
125         TString histit2=Form("(dEdx_l%d-BetheBloch)/BetheBloch pbin%d part%d",fgkLayerCode[iVal],iPBin,fgkPdgCode[iSpecie]);
126         if(iVal==kNValuesdEdx-1) histit1=Form("(dEdx_truncmean-BetheBloch)/BetheBloch pbin%d part%d",iPBin,fgkPdgCode[iSpecie]);
127         fHistoDeltadEdx[GetIndex(iSpecie,iPBin,iVal)]=new TH1F(hisnam2.Data(),histit2.Data(),100,-0.5,0.5);
128         fHistoDeltadEdx[GetIndex(iSpecie,iPBin,iVal)]->GetXaxis()->SetTitle(histit2.Data());
129       }
130     }
131   }
132
133   fHistodEdxVsP=new TH2F*[GetNHistos2()];
134   for(Int_t iSpecie=0; iSpecie<kNParticles; iSpecie++){
135     for(Int_t iVal=0; iVal<kNValuesdEdx; iVal++){
136       TString hisnam=Form("hdEdx%dVsPpart%d",fgkLayerCode[iVal],fgkPdgCode[iSpecie]);
137       TString histit=Form("dEdx layer %d (keV) vs P part%d",fgkLayerCode[iVal],fgkPdgCode[iSpecie]);
138       if(iVal==kNValuesdEdx-1) histit=Form("dEdx trunc. mean (keV) vs P part%d",fgkPdgCode[iSpecie]);      
139       fHistodEdxVsP[GetIndex2(iSpecie,iVal)]=new TH2F(hisnam.Data(),histit.Data(),50,fPMin,fPMax,50,0.,500.);
140       histit.ReplaceAll(" vs P "," ");
141       fHistodEdxVsP[GetIndex2(iSpecie,iVal)]->GetYaxis()->SetTitle(histit.Data());
142       fHistodEdxVsP[GetIndex2(iSpecie,iVal)]->GetXaxis()->SetTitle("Momentum (GeV/c)");
143     }
144   }
145 }
146 //______________________________________________________________________
147 void AliITSdEdxAnalyzer:: WriteHistos(TString filename) const {
148   // stores output histogrma to file
149   TFile *outfile=new TFile(filename.Data(),"recreate");
150   outfile->cd();
151   for(Int_t i=0; i<GetNHistos();i++){ 
152     fHistodEdx[i]->Write();
153     fHistoDeltadEdx[i]->Write();
154   }
155   for(Int_t i=0; i<GetNHistos2();i++) fHistodEdxVsP[i]->Write();
156   outfile->Close();
157   delete outfile;
158 }
159 //______________________________________________________________________
160 void AliITSdEdxAnalyzer::ReadEvent(const AliESDEvent* esd, AliStack* stack){
161   // Fill histos 
162   // if stack is !=0 MC truth is used to define particle specie
163   // else TPC pid is used to define particle specie
164
165   if(!esd) AliFatal("Bad ESD event");
166   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
167     AliESDtrack* track = esd->GetTrack(iTrack);
168     Double_t trmean=track->GetITSsignal();
169     Double_t sig[4];
170     track->GetITSdEdxSamples(sig);
171     Double_t preco=track->P();
172     Int_t iPBin=GetMomentumBin(preco);
173     if(iPBin==-1) continue;
174     Int_t iSpecie=-1;
175     Double_t dedxbb=0;
176     if(stack){
177       Int_t lab=track->GetLabel();
178       if(lab<0) continue;
179       TParticle* part=(TParticle*)stack->Particle(lab);
180       Int_t absPdgCode=TMath::Abs(part->GetPdgCode());
181       iSpecie=GetSpecieBin(absPdgCode);
182       if(iSpecie==-1) continue;
183       dedxbb=BetheBloch(part);
184     }else{   
185       iSpecie=GetPaticleIdFromTPC(track);
186       if(iSpecie==-1) continue;
187       TParticlePDG* p=TDatabasePDG::Instance()->GetParticle(fgkPdgCode[iSpecie]);    
188       dedxbb=BetheBloch(preco,p->Mass());
189     }
190     for(Int_t ilay=0;ilay<4;ilay++){    
191       if(sig[ilay]>0.){
192         fHistodEdx[GetIndex(iSpecie,iPBin,ilay)]->Fill(sig[ilay]);
193         fHistoDeltadEdx[GetIndex(iSpecie,iPBin,ilay)]->Fill((sig[ilay]-dedxbb)/dedxbb);
194         fHistodEdxVsP[GetIndex2(iSpecie,ilay)]->Fill(preco,sig[ilay]);
195       }
196     }
197     if(trmean>0.){
198       fHistodEdx[GetIndex(iSpecie,iPBin,4)]->Fill(trmean);
199       fHistoDeltadEdx[GetIndex(iSpecie,iPBin,4)]->Fill((trmean-dedxbb)/dedxbb);
200       fHistodEdxVsP[GetIndex2(iSpecie,4)]->Fill(preco,trmean);
201     }
202   }
203 }
204 //______________________________________________________________________
205 Int_t AliITSdEdxAnalyzer::GetPaticleIdFromTPC(const AliESDtrack* track) const {
206   // Returns the particle specie as identified by TPC
207   Double_t tpcpid[AliPID::kSPECIES];
208   track->GetTPCpid(tpcpid);
209   Int_t maxPos=-1;
210   Float_t maxProb=0.;
211   for(Int_t is=0; is<AliPID::kSPECIES; is++){
212     if(tpcpid[is]>maxProb){
213       maxProb=tpcpid[is];
214       maxPos=is;
215     }
216   }
217   Int_t iSpecie=-1;
218   if(maxProb>fTPCpidCut){
219     if(maxPos==AliPID::kPion) iSpecie=0;
220     else if(maxPos==AliPID::kKaon) iSpecie=1;
221     else if(maxPos==AliPID::kProton) iSpecie=2;    
222   }
223   return iSpecie;
224 }
225 //______________________________________________________________________
226 Double_t AliITSdEdxAnalyzer::BetheBloch(const Float_t p, const Float_t m) const {
227   // Computes expected dE/dx value for given particle specie and momentum
228   static AliITSPIDResponse pidResponse;
229   Double_t dedxbb=0.;
230   if(fBBmodel==0){
231     Double_t betagamma=p/m;
232     Double_t conv=fDensity*1E6*fThickness/116.24*fMIP;
233     dedxbb=conv*AliExternalTrackParam::BetheBlochSolid(betagamma);
234   }else if(fBBmodel==1){
235     dedxbb=pidResponse.Bethe(p,m);
236   }
237   return dedxbb;
238 }
239 //______________________________________________________________________
240 TGraph* AliITSdEdxAnalyzer::GetBetheBlochGraph(const Int_t pdgcode) const {
241  // Fills a TGraph with Bethe-Bloch expe
242  TGraph* g=new TGraph(0);
243  TParticlePDG* part=TDatabasePDG::Instance()->GetParticle(pdgcode);
244  Float_t mass=part->Mass();
245  for(Int_t ip=0; ip<100;ip++){
246    Float_t p=fPMin+(ip+0.5)*(fPMax-fPMin)/100.;
247    g->SetPoint(ip,p,BetheBloch(p,mass));
248  }
249  g->GetXaxis()->SetTitle("Momentum (GeV/c)");
250  g->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
251  return g;
252 }