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