fix in calling of gaussian spread function
[u/mrichter/AliRoot.git] / ITS / AliITSdEdxAnalyzer.cxx
CommitLineData
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"
30#include "AliLog.h"
27199d86 31#include "AliESDEvent.h"
32#include "AliStack.h"
75c74fab 33#include "AliITSdEdxAnalyzer.h"
34#include "AliExternalTrackParam.h"
10d100d4 35//#include "AliITSpidESD.h"
36#include "AliITSPIDResponse.h"
75c74fab 37#include "AliPID.h"
38
39ClassImp(AliITSdEdxAnalyzer)
40
41const Int_t AliITSdEdxAnalyzer::fgkPdgCode[kNParticles]={211,321,2212};
42const Int_t AliITSdEdxAnalyzer::fgkLayerCode[kNValuesdEdx]={3,4,5,6,0};
43
44//______________________________________________________________________
45AliITSdEdxAnalyzer::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//______________________________________________________________________
63AliITSdEdxAnalyzer::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//______________________________________________________________________
81AliITSdEdxAnalyzer::~AliITSdEdxAnalyzer(){
82 // destructor
83 DeleteHistos();
84}
85//______________________________________________________________________
86void 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//______________________________________________________________________
95void 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//______________________________________________________________________
114void 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//______________________________________________________________________
149void 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 162void 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//______________________________________________________________________
207Int_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//______________________________________________________________________
228Double_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=fDensity*1E6*fThickness/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//______________________________________________________________________
242TGraph* 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}