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