New class for dE/dx analysis (comparison with Bethe-Bloch, check of response function...
[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 "AliITSdEdxAnalyzer.h"
30 #include "AliExternalTrackParam.h"
31 #include "AliITSpidESD.h"
32 #include "AliPID.h"
33
34 ClassImp(AliITSdEdxAnalyzer)
35
36 const Int_t AliITSdEdxAnalyzer::fgkPdgCode[kNParticles]={211,321,2212};
37 const Int_t AliITSdEdxAnalyzer::fgkLayerCode[kNValuesdEdx]={3,4,5,6,0};
38
39 //______________________________________________________________________
40 AliITSdEdxAnalyzer::AliITSdEdxAnalyzer() :  
41   TObject(),
42   fNPBins(10),
43   fPMin(0.1),
44   fPMax(1.1),
45   fHistodEdx(0),
46   fHistoDeltadEdx(0),
47   fHistodEdxVsP(0),
48   fThickness(0.03),
49   fDensity(2.33),
50   fBBmodel(0),
51   fMIP(82.),
52   fTPCpidCut(0.5)
53 {
54   // default constructor
55   BookHistos();
56 }
57 //______________________________________________________________________
58 AliITSdEdxAnalyzer::AliITSdEdxAnalyzer(const Int_t npbins, const Float_t pmin, const Float_t pmax) :  
59   TObject(),
60   fNPBins(npbins),
61   fPMin(pmin),
62   fPMax(pmax),
63   fHistodEdx(0),
64   fHistoDeltadEdx(0),
65   fHistodEdxVsP(0),
66   fThickness(0.03),
67   fDensity(2.33),
68   fBBmodel(0),
69   fMIP(82.),
70   fTPCpidCut(0.5)
71 {
72   // standard constructor
73   BookHistos();
74 }
75 //______________________________________________________________________
76 AliITSdEdxAnalyzer::~AliITSdEdxAnalyzer(){
77   // destructor
78   DeleteHistos();
79 }
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
83   DeleteHistos();
84   fNPBins=npbins;
85   fPMin=pmin;
86   fPMax=pmax;
87   BookHistos();
88 }
89 //______________________________________________________________________
90 void AliITSdEdxAnalyzer::DeleteHistos(){
91   //
92   if(fHistodEdx){
93     for(Int_t i=0; i<GetNHistos();i++) delete fHistodEdx[i];
94     delete [] fHistodEdx;
95     fHistodEdx=0;
96   }
97   if(fHistoDeltadEdx){
98     for(Int_t i=0; i<GetNHistos();i++) delete fHistoDeltadEdx[i];
99     delete [] fHistoDeltadEdx;
100     fHistoDeltadEdx=0;
101   }
102   if(fHistodEdxVsP){
103     for(Int_t i=0; i<GetNHistos2();i++) delete fHistodEdxVsP[i];
104     delete [] fHistodEdxVsP;
105     fHistodEdxVsP=0;
106   }
107 }
108 //______________________________________________________________________
109 void AliITSdEdxAnalyzer::BookHistos(){
110   //
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());
126       }
127     }
128   }
129
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)");
140     }
141   }
142 }
143 //______________________________________________________________________
144 void AliITSdEdxAnalyzer:: WriteHistos(TString filename) const {
145   TFile *outfile=new TFile(filename.Data(),"recreate");
146   outfile->cd();
147   for(Int_t i=0; i<GetNHistos();i++){ 
148     fHistodEdx[i]->Write();
149     fHistoDeltadEdx[i]->Write();
150   }
151   for(Int_t i=0; i<GetNHistos2();i++) fHistodEdxVsP[i]->Write();
152   outfile->Close();
153   delete outfile;
154 }
155 //______________________________________________________________________
156 void AliITSdEdxAnalyzer::ReadEvent(AliESDEvent* esd, AliStack* stack){
157   // Fill histos 
158   // if stack is !=0 MC truth is used to define particle specie
159   // else TPC pid is used to define particle specie
160
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();
165     Double_t sig[4];
166     track->GetITSdEdxSamples(sig);
167     Double_t preco=track->P();
168     Int_t iPBin=GetMomentumBin(preco);
169     if(iPBin==-1) continue;
170     Int_t iSpecie=-1;
171     Double_t dedxbb=0;
172     if(stack){
173       Int_t lab=track->GetLabel();
174       if(lab<0) continue;
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);
180     }else{   
181       iSpecie=GetPaticleIdFromTPC(track);
182       if(iSpecie==-1) continue;
183       TParticlePDG* p=TDatabasePDG::Instance()->GetParticle(fgkPdgCode[iSpecie]);    
184       dedxbb=BetheBloch(preco,p->Mass());
185     }
186     for(Int_t ilay=0;ilay<4;ilay++){    
187       if(sig[ilay]>0.){
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]);
191       }
192     }
193     if(trmean>0.){
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);
197     }
198   }
199 }
200 //______________________________________________________________________
201 Int_t AliITSdEdxAnalyzer::GetPaticleIdFromTPC(const AliESDtrack* track) const {
202   //
203   Double_t tpcpid[AliPID::kSPECIES];
204   track->GetTPCpid(tpcpid);
205   Int_t maxPos=-1;
206   Float_t maxProb=0.;
207   for(Int_t is=0; is<AliPID::kSPECIES; is++){
208     if(tpcpid[is]>maxProb){
209       maxProb=tpcpid[is];
210       maxPos=is;
211     }
212   }
213   Int_t iSpecie=-1;
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;    
218   }
219   return iSpecie;
220 }
221 //______________________________________________________________________
222 Double_t AliITSdEdxAnalyzer::BetheBloch(const Float_t p, const Float_t m) const {
223   //
224   Double_t dedxbb=0.;
225   switch(fBBmodel){
226   case 0:
227     Double_t betagamma=p/m;
228     Double_t conv=fDensity*1E6*fThickness/116.31*fMIP;
229     dedxbb=conv*AliExternalTrackParam::BetheBlochSolid(betagamma);
230     break;
231   case 1:
232     dedxbb=fMIP*AliITSpidESD::Bethe(p,m);
233     break;
234   }
235   return dedxbb;
236 }
237 //______________________________________________________________________
238 TGraph* AliITSdEdxAnalyzer::GetBetheBlochGraph(const Int_t pdgcode) const {
239   // Fills a TGraph with Bethe-Bloch expe
240   TGraph* g=new TGraph(0);
241   TParticlePDG* p=TDatabasePDG::Instance()->GetParticle(pdgcode);
242   Float_t mass=p->Mass();
243   for(Int_t ip=0; ip<100;ip++){
244     Float_t p=fPMin+(ip+0.5)*(fPMax-fPMin)/100.;
245     g->SetPoint(ip,p,BetheBloch(p,mass));
246   }
247   g->GetXaxis()->SetTitle("Momentum (GeV/c)");
248   g->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
249   return g;
250 }