]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEdca.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixing conding violations (Matthieu)
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEdca.cxx
CommitLineData
70da6c5a 1/*************************************************************************
2* Copyright(c) 1998-1999, 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// Class for impact parameter (DCA) of charged particles
17// Study resolution and pull: prepare for beauty study
18//
19// Authors:
20// Hongyan Yang <hongyan@physi.uni-heidelberg.de>
21// Carlo Bombonati <carlo.bombonati@cern.ch>
22//
23
24#include "TMath.h"
25#include "TH1F.h"
26#include "TList.h"
27#include <TParticle.h>
28#include "AliMCParticle.h"
29#include "AliESDtrack.h"
30#include "AliESDEvent.h"
31#include "AliMCEvent.h"
32
33#include "AliHFEdca.h"
34
35ClassImp(AliHFEdca)
36
37//________________________________________________________________________
38const Char_t* AliHFEdca::fgkParticles[12] = {
39 // particles name
40 "electron", "muonMinus","pionMinus", "kaonMinus", "protonMinus",
41 "positron", "muonPlus", "pionPlus", "kaonPlus", "protonPlus",
42 "allNegative", "allPositive"
43};
44//________________________________________________________________________
45const Int_t AliHFEdca::fgkColorPart[12] = {
46 // colors assigned to particles
47 kRed, kBlue, kGreen+2, kYellow+2, kMagenta,
48 kRed+2, kBlue+2, kGreen+4, kYellow+4, kMagenta+2,
49 kBlack, kGray+1
50};
51
52
53
54//________________________________________________________________________
55const Float_t AliHFEdca::fgkPtIntv[44] = {
56 // define pT bins
57 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
58 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8,
59 3.0, 3.3, 3.6, 3.9, 4.2, 4.5, 5.0, 5.5, 6.0, 6.5,
60 7.0, 8.0, 9.0, 10., 11., 12., 13., 14., 15., 16.,
61 17., 18., 19., 20.};
62
63//________________________________________________________________________
64const Char_t *AliHFEdca::fgkDcaVar[2] = {
65 "deltaDcaXY", "deltaDcaZ"};
66
67//________________________________________________________________________
68const Char_t *AliHFEdca::fgkDcaVarTitle[2] ={
69 ";residual #Delta(d_{xy}) [#mum];couts", ";residual #Delta(d_{z}) [#mum];counts"};
70
71
72//________________________________________________________________________
73const Char_t *AliHFEdca::fgkPullDcaVar[2] = {
74 "pullDcaXY", "pullDcaZ"
75};
76
77//________________________________________________________________________
78const Char_t *AliHFEdca::fgkPullDcaVarTitle[2] = {
79 ";(dca_{xy}^{ESD}-dca_{xy}^{MC})/(#sigma_{track}#oplus#sigma_{Vxy});counts",
80 ";(dca_{z}^{ESD}-dca_{z}^{MC})/(#sigma_{track}#oplus#sigma_{Vz}); counts"
81};
82
83//________________________________________________________________________
84AliHFEdca::AliHFEdca():
85 fResidualList(0x0)
86 , fPullList(0x0)
87
88{
89 // default constructor
90
91}
92
93//________________________________________________________________________
94AliHFEdca::AliHFEdca(const AliHFEdca &dca):
95 TObject(dca)
96 , fResidualList(0x0)
97 , fPullList(0x0)
98
99{
100 // copy constructor
101
102}
103//_______________________________________________________________________________________________
104AliHFEdca&
105AliHFEdca::operator=(const AliHFEdca &)
106{
107 //
108 // Assignment operator
109 //
110
111 Printf("Not yet implemented.");
112 return *this;
113}
114
115//________________________________________________________________________
116AliHFEdca::~AliHFEdca()
117{
118 // default destructor
119
120 for(Int_t j=0; j<kNParticles; j++){
121 for(Int_t i=0; i<kNPtBins; i++){
122 if(fHistDcaXYRes[j][i]) delete fHistDcaXYRes[j][i];
123 if(fHistDcaZRes[j][i]) delete fHistDcaZRes[j][i];
124 if(fHistDcaXYPull[j][i]) delete fHistDcaXYPull[j][i];
125 if(fHistDcaXYPull[j][i]) delete fHistDcaZPull[j][i];
126 }
127 }
128
129 if(fResidualList) delete fResidualList;
130 if(fPullList) delete fPullList;
131
132 Printf("analysis done\n");
133}
134
135//________________________________________________________________________
136void AliHFEdca::InitAnalysis(){
137
138 Printf("initialize analysis\n");
139
140}
141
142
143//________________________________________________________________________
144void AliHFEdca::PostAnalysis() const
145{
146 // do fit
147 // moved to dcaPostAnalysis.C
148
149}
150//________________________________________________________________________
151void AliHFEdca::CreateHistogramsResidual(TList *residualList){
152 // define histogram
153 // 1. residual
154
155 // for residuals
156 fHistDcaXYRes[kNParticles][kNPtBins]=0x0;
157 fHistDcaZRes[kNParticles][kNPtBins]=0x0;
158
159 const Int_t nBins = 1000;
160 const Float_t maxXYBin = 1000.;
161 const Float_t maxZBin = 1000.;
162
163
164 for(Int_t k=0; k<kNDcaVar; k++){
165 TString histTitle((const char*)fgkDcaVarTitle[k]);
166
167 for(Int_t j=0; j<kNParticles; j++){
168 for(Int_t i=0; i<kNPtBins; i++){
169
170 TString histName((const char*)fgkParticles[j]);
171
172 histName += Form("_%s_pT-%.1f-%.1f", (const char*)fgkDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
173
174 if(k==0){
175 fHistDcaXYRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxXYBin, maxXYBin);
176 fHistDcaXYRes[j][i]->SetLineColor((const int)fgkColorPart[j]);
177 }
178 if(k==1){
179 fHistDcaZRes[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, -maxZBin, maxZBin);
180 fHistDcaZRes[j][i]->SetLineColor((const int)fgkColorPart[j]);
181 }
182 } // 43 pt bins
183 } //12 nparticles
184 } // 2 dca var
185
186 // TList *fResidualList = 0;
187 residualList->SetOwner();
188 residualList->SetName("residual");
189 for(Int_t iPart=0; iPart<kNParticles; iPart++){
190 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
191 residualList->Add(fHistDcaXYRes[iPart][iPtBin]);
192 residualList->Add(fHistDcaZRes[iPart][iPtBin]);
193 } // loop over pt bins
194 } // loop over particles (pos, neg)
195
196
197
198
199}
200
201
202//________________________________________________________________________
203void AliHFEdca::CreateHistogramsPull(TList *pullList){
204 // define histogram
205 // 2. pull
206
207 const Int_t nBins = 1000;
208 const Float_t maxXYBin = 20.;
209 const Float_t maxZBin = 20.;
210
211
212 // for pull -----------------------------------------------------------------------
213 fHistDcaXYPull[kNParticles][kNPtBins]=0x0;
214 fHistDcaZPull[kNParticles][kNPtBins]=0x0;
215
216
217 for(Int_t k=0; k<kNDcaVar; k++){
218 TString histTitle((const char*)fgkPullDcaVarTitle[k]);
219
220 for(Int_t j=0; j<kNParticles; j++){
221 for(Int_t i=0; i<kNPtBins; i++){
222
223 TString histName((const char*)fgkParticles[j]);
224
225 histName += Form("_%s_pT-%.1f-%.1f", (const char*)fgkPullDcaVar[k], fgkPtIntv[i], fgkPtIntv[i+1]);
226
227 if(k==0){
228 fHistDcaXYPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxXYBin, 1+maxXYBin);
229 fHistDcaXYPull[j][i]->SetLineColor((const int)fgkColorPart[j]);
230 }
231 if(k==1){
232 fHistDcaZPull[j][i] = new TH1F((const char*)histName, (const char*)histTitle, nBins, 1-maxZBin, 1+maxZBin);
233 fHistDcaZPull[j][i]->SetLineColor((const int)fgkColorPart[j]);
234 }
235 } // 43 pt bins
236 } //6 nparticles
237 } // 2 dca var
238
239 // TList *fPullList = 0;
240 pullList->SetOwner();
241 pullList->SetName("pull");
242 for(Int_t iPart=0; iPart<kNParticles; iPart++){
243 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
244 pullList->Add(fHistDcaXYPull[iPart][iPtBin]);
245 pullList->Add(fHistDcaZPull[iPart][iPtBin]);
246 } // loop over pt bins
247 } // loop over particles (pos, neg)
248
249
250
251}
252
253
254//_______________________________________________________________________________________________
255void AliHFEdca::FillHistograms(AliESDEvent * const esdEvent, AliESDtrack * const track, AliMCEvent * const mcEvent)
256{
257// filling historgams track by track
258
259// obtaining reconstructed dca ------------------------------------------------------------------
260 Float_t esdpx = track->Px();
261 Float_t esdpy = track->Py();
262 Float_t esdpt = TMath::Sqrt(esdpx*esdpx+esdpy*esdpy);
263 Float_t b[2]; // dca in cm
264 Float_t bCov[3]; // covariance matrix
265 track->GetImpactParameters(b,bCov);
266
267// obtaining errors of dca ------------------------------------------------------------------
268 const AliESDVertex *primVtx = esdEvent->GetPrimaryVertex();
269 Float_t magneticField = 5; // initialized as 5kG
270 magneticField = esdEvent->GetMagneticField(); // in kG
271
272 Double_t dz[2]; // error of dca in cm
273 Double_t covardz[3];
274 track->PropagateToDCA(primVtx,magneticField, 1000., dz, covardz);
275
276 // calculate mcDca ------------------------------------------------------------------
277
278 AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(TMath::Abs(track->GetLabel())));
279 TParticle *part = mctrack->Particle();
280 Int_t pdg = part->GetPdgCode();
281 Int_t charge = 1;
282 if(pdg==kPDGelectron || pdg==kPDGmuon
283 || pdg==-kPDGpion || pdg==-kPDGkaon || pdg==-kPDGproton) charge = -1;
284
285 Float_t vx = part->Vx(); // in cm
286 Float_t vy = part->Vy(); // in cm
287 Float_t vz = part->Vz(); // in cm
288
289 Float_t vxy = TMath::Sqrt(vx*vx+vy*vy);
290
291 Float_t mcpx = part->Px();
292 Float_t mcpy = part->Py();
293 Float_t mcpt = TMath::Sqrt(mcpx*mcpx+mcpy*mcpy);
294
295 const Float_t conv[2] = {1.783/1.6, 2.99792458};
296 Float_t radiusMc = mcpt/(TMath::Abs(magneticField)/10.)*conv[0]*conv[1]; // pt in GeV/c, magnetic field in Tesla
297
298 Float_t nx = esdpx/mcpt;
299 Float_t ny = esdpy/mcpt;
300
301 Float_t radius;
302 radius = TMath::Abs(radiusMc);
303 Float_t mcDcaXY = (radius - TMath::Sqrt(vxy*vxy/100./100. + radius*radius + 2*radius*charge*(vx*ny-vy*nx)/100.)) ; // in meters
304
305// printf("magnetic Field = %.3f \t", magneticField);
306// printf("pt=esd %.3f/mc %.3f GeV/c, radius=esd %.3f/mc %.3f meter \n", esdpt, mcpt, radiusEsd, radiusMc);
307// printf("mcDcaXY=%.5f micron, esdDcaXY=%.5f micron \t ", mcDcaXY*1.e6, b[0]*1e4);
308// printf("mcDcaZ=%.5f micron, esdDcaZ=%.5f micron\n\n", vz*1e6, b[1]*1e4);
309
310 Double_t mcDca[2] = {mcDcaXY*100, vz}; // in cm
311
312 Double_t residual[2] = {0, 0};
313 Double_t pull[2] = {0, 0};
314 Double_t error[2] ={TMath::Sqrt(covardz[0]), TMath::Sqrt(covardz[2])};
315 for(Int_t i=0; i<2; i++){
316 residual[i] = dz[i] - mcDca[i]; // in centimeters
317 if(error[i]!=0)pull[i] = residual[i]/error[i]; // unitless
318 // printf("error[%d]=%.6f residual[%d]=%.6f pull[%d]=%.6f\n", i, error[i], i, residual[i], i, pull[i]);
319 }
320
321 Int_t fPdgParticle[10] = {
322 kPDGelectron, kPDGmuon, -kPDGpion, -kPDGkaon, -kPDGproton,
323 -kPDGelectron, -kPDGmuon, kPDGpion, kPDGkaon, kPDGproton};
324
325 for(Int_t iPart=0; iPart<kNParticles-2; iPart++){
326
327 // identified ones
328 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
329 if(pdg==fPdgParticle[iPart] && (esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1])) {
330 fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1.0e4); // in microns
331 fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1.0e4); // in microns
332 fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
333 fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
334 }
335 else
336 continue;
337 }
338 }
339
340 // for charged particles
341 for(Int_t iPtBin=0; iPtBin<kNPtBins; iPtBin++){
342 if(esdpt>fgkPtIntv[iPtBin] && esdpt<=fgkPtIntv[iPtBin+1]){
343 Int_t iPart = 10;
344 if(charge>0) iPart = 11;
345 fHistDcaXYRes[iPart][iPtBin]->Fill(residual[0]*1e4);
346 fHistDcaZRes[iPart][iPtBin]->Fill(residual[1]*1e4);
347 fHistDcaXYPull[iPart][iPtBin]->Fill(pull[0]);
348 fHistDcaZPull[iPart][iPtBin]->Fill(pull[1]);
349 }
350 else
351 continue;
352 }
353
354}
355