prepare PID LQ ref maker for production
[u/mrichter/AliRoot.git] / PWG1 / TRD / info / AliTRDv0Info.cxx
CommitLineData
1ee39b3a 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/* $Id: AliTRDv0Info.cxx 27496 2008-07-22 08:35:45Z cblume $ */
17
18////////////////////////////////////////////////////////////////////////////
19// //
20// Reconstruction QA //
21// //
22// Gathers all information necessary for reference data selection about //
23// the track and (in case) its corresponding V0. //
24// Carries out the selection of electrons (from gamma conversions), //
25// pions (from K0s decays) and protons (from Lambda and Anti-Lambda //
26// decays) by cuts specific for the respective decay and particle //
27// species. //
28// (M.Heide, 2009/10/06) //
29// //
30// Authors: //
31// Alex Bercuci <A.Bercuci@gsi.de> //
32// Alex Wilk <wilka@uni-muenster.de> //
33// Markus Heide <mheide@uni-muenster.de> //
34// //
35////////////////////////////////////////////////////////////////////////////
36
37#include "TMath.h"
38
39#include "AliESDtrack.h"
40#include "AliESDv0.h"
41#include "AliESDInputHandler.h"
42#include "AliAnalysisManager.h"
43#include "AliLog.h"
44
45#include "AliTRDv0Info.h"
46#include "AliTRDtrackInfo.h"
47#include "AliTRDtrackInfo.h"
48
49ClassImp(AliTRDv0Info)
50
51//_________________________________________________
52AliTRDv0Info::AliTRDv0Info()
53 : TObject()
54 ,fESD(0x0)
55 ,fHasV0(0)
56 ,fQuality(0)
57 ,fMomentum(0)
58 ,fDCA(10)
59 ,fPointingAngle(10)
60 ,fOpenAngle(10)
61 ,fPsiPair(99)
62 ,fMagField(0)
63 ,fRadius(0)
64 ,fTrackID(0)
65 ,fV0Momentum(0)
66 ,fTrackP(0x0)
67 ,fTrackN(0x0)
68 ,fTrack(0x0)
69 ,fNindex(0)
70 ,fPindex(0)
71{
72 //
73 // Default constructor
74 //
75
76 memset(fPplus, 0, 2*kNlayer*sizeof(Float_t));
77 memset(fPminus, 0, 2*kNlayer*sizeof(Float_t));
78 memset(fDetPID, 0, 2*kNDaughters*kNDetectors*AliPID::kSPECIES*sizeof(Float_t));
79 memset(fInvMass, 0, kNMomBins*kNDecays*sizeof(Double_t));
80
81 /////////////////////////////////////////////////////////////////////////////
82 //Set Cut values: First specify decay in brackets, then the actual cut value!
83 /////////////////////////////////////////////////////////////////////////////
84
85 //Upper limit for distance of closest approach of two daughter tracks :
86 fUpDCA[kGamma] = 0.25;
87 fUpDCA[kK0s] = 0.25;
88 fUpDCA[kLambda] = 0.25;
89 fUpDCA[kAntiLambda] = 0.25;
90
91 //Upper limit for pointing angle (= angle between between vector from primary to secondary vertex and reconstructed momentum of V0 mother particle) :
92 fUpPointingAngle[kGamma] = 0.03;
93 fUpPointingAngle[kK0s] = 0.03;
94 fUpPointingAngle[kLambda] = 0.03;
95 fUpPointingAngle[kAntiLambda] = 0.03;
96
97 //Upper limit for invariant mass of V0 mother :
98 fUpInvMass[kGamma][0] = 0.04;// second pair of brackets is for momentum bin: 0: below mother momentm of 2.5 GeV
99 fUpInvMass[kGamma][1] = 0.07;//1: above 2.5 GeV
100 fUpInvMass[kK0s][0] = fUpInvMass[kK0s][1] = 0.51;
101 fUpInvMass[kLambda][0] = fUpInvMass[kLambda][1] = 1.22;
102 fUpInvMass[kAntiLambda][0] = fUpInvMass[kAntiLambda][1] = 1.22;
103
104 //Lower limit for invariant mass of V0 mother :
105 fDownInvMass[kGamma] = -1.;
106 fDownInvMass[kK0s] = 0.49;
107 fDownInvMass[kLambda] = 1.;
108 fDownInvMass[kAntiLambda] = 1.;
109
110 //Lower limit for distance from secondary vertex to primary vertex in x-y plane :
111 fDownRadius[kGamma] = 5.7;
112 fDownRadius[kK0s] = 0.;
113 fDownRadius[kLambda] = 10.;
114 fDownRadius[kAntiLambda] = 10.;
115
116 //Upper limit for distance from secondary vertex to primary vertex in x-y plane :
117 fUpRadius[kGamma] = 1000.;
118 fUpRadius[kK0s] = 1000.;
119 fUpRadius[kLambda] = 1000.;
120 fUpRadius[kAntiLambda] = 1000.;
121
122 //Upper limit for opening angle between two daughter tracks (characteristically near zero for conversions) :
123 fUpOpenAngle[kGamma] = 0.1;
124 fUpOpenAngle[kK0s] = 3.15;
125 fUpOpenAngle[kLambda] = 3.15;
126 fUpOpenAngle[kAntiLambda] = 3.15;
127
128 //Upper limit for angle between daughter momentum plane and plane perpendicular to magnetic field (characteristically around zero for conversions) :
129 fUpPsiPair[kGamma] = 0.1;
130 fUpPsiPair[kK0s] = 1.6;
131 fUpPsiPair[kLambda] = 1.6;
132 fUpPsiPair[kAntiLambda] = 1.6;
133
134 //Lower limit for likelihood value of TPC PID :
135 fDownTPCPIDneg[AliPID::kElectron] = 0.21;
136 fDownTPCPIDpos[AliPID::kElectron] = 0.21;
137
138 fDownTPCPIDneg[AliPID::kMuon] = 0.21;
139 fDownTPCPIDpos[AliPID::kMuon] = 0.21;
140
141 fDownTPCPIDneg[AliPID::kPion] = 0.21;
142 fDownTPCPIDpos[AliPID::kPion] = 0.21;
143
144 fDownTPCPIDneg[AliPID::kKaon] = 0.21;
145 fDownTPCPIDpos[AliPID::kKaon] = 0.21;
146
147 fDownTPCPIDneg[AliPID::kProton] = 0.21;
148 fDownTPCPIDpos[AliPID::kProton] = 0.21;
149 //////////////////////////////////////////////////////////////////////////////////
150
151}
152
153//_________________________________________________
154void AliTRDv0Info::GetESDv0Info(AliESDv0 *esdv0)
155{//Gets values of ESDv0 and daughter track properties
156 //See header file for description of variables
157
158 Int_t part1 = -1;
159 Int_t part2 = -1;
160
161 fQuality = Quality(esdv0);//Attributes an Int_t to the V0 due to quality cuts (= 1 if V0 is accepted, other integers depending on cut which excludes the vertex)
162
163 fRadius = Radius(esdv0);//distance from secondary vertex to primary vertex in x-y plane
164
165 fDCA = esdv0->GetDcaV0Daughters();//distance of closest approach of two daughter tracks
166
167 fPointingAngle = TMath::ACos(esdv0->GetV0CosineOfPointingAngle());// pointing angle (= angle between between vector from primary to secondary vertex and reconstructed momentum of V0 mother particle)
168
169 fOpenAngle = OpenAngle(esdv0);//Opening angle between two daughter tracks
170
171 fPsiPair = PsiPair(esdv0);//Angle between daughter momentum plane and plane perpendicular to magnetic field
172
173 fV0Momentum = V0Momentum(esdv0);//Reconstructed momentum of the mother particle
174
175 for(Int_t idecay = 0; idecay < kNDecays; idecay++)//4 decay types : conversions, K0s, Lambda, Anti-Lambda
176 //five particle types: electrons, muons, pions, kaons, protons (muons and kaons not involved)
177 {
178 if(idecay == kLambda)//protons and pions from Lambda
179 {
180 part1 = AliPID::kProton;
181 part2 = AliPID::kPion;
182 }
183 else if(idecay == kAntiLambda)//antiprotons and pions from Anti-Lambda
184 {
185 part1 = AliPID::kPion;
186 part2 = AliPID::kProton;
187 }
188 else if(idecay == kK0s)//pions from K0s
189 part1 = part2 = AliPID::kPion;
190 else if(idecay == kGamma)//electrons from conversions
191 part1 = part2 = AliPID::kElectron;
192
193 fInvMass[idecay] = InvMass(part1, part2, esdv0);//Calculate invariant mass for all of our four supposed decays
194 }
195 GetDetectorPID();//Gets all likelihood values from TPC, TOF and ITS PID for the fDetPID[kNDaughters][kNDetectors][AliPID::kSPECIES] array
196
197
198}
199//_________________________________________________
200Float_t AliTRDv0Info::V0Momentum(AliESDv0 *esdv0) const
201{
202 //
203 // Reconstructed momentum of V0 mother particle
204 //
205
206 Double_t mn[3] = {0,0,0};
207 Double_t mp[3] = {0,0,0};
208
209
210 esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
211 esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
212
213
214 return TMath::Sqrt((mn[0]+mp[0])*(mn[0]+mp[0]) + (mn[1]+mp[1])*(mn[1]+mp[1])+(mn[2]+mp[2])*(mn[2]+mp[2]));
215}
216
217//_________________________________________________
218Double_t AliTRDv0Info::InvMass(Int_t part1, Int_t part2, AliESDv0 *esdv0) const
219{
220 //
221 // Invariant mass of reconstructed V0 mother
222 //
223
224 const Double_t kpmass[5] = {AliPID::ParticleMass(AliPID::kElectron),AliPID::ParticleMass(AliPID::kMuon),AliPID::ParticleMass(AliPID::kPion),AliPID::ParticleMass(AliPID::kKaon),AliPID::ParticleMass(AliPID::kProton)};
225 //Masses of electrons, muons, pions, kaons and protons, as implemented in ROOT
226
227
228 Double_t mn[3] = {0,0,0};
229 Double_t mp[3] = {0,0,0};
230
231 esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
232 esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
233
234 Double_t mass1 = kpmass[part1];//sets supposed rest masses for both daughters
235 Double_t mass2 = kpmass[part2];
236
237 //Calculate daughters' energies :
238 Double_t e1 = TMath::Sqrt(mass1*mass1+
239 mp[0]*mp[0]+
240 mp[1]*mp[1]+
241 mp[2]*mp[2]);
242 Double_t e2 = TMath::Sqrt(mass2*mass2+
243 mn[0]*mn[0]+
244 mn[1]*mn[1]+
245 mn[2]*mn[2]);
246
247 //Sum of daughter momenta :
248 Double_t momsum =
249 (mn[0]+mp[0])*(mn[0]+mp[0])+
250 (mn[1]+mp[1])*(mn[1]+mp[1])+
251 (mn[2]+mp[2])*(mn[2]+mp[2]);
252
253 //invariant mass :
254 Double_t mInv = TMath::Sqrt((e1+e2)*(e1+e2)-momsum);
255
256 return mInv;
257
258}
259//_________________________________________________
260Float_t AliTRDv0Info::OpenAngle(AliESDv0 *esdv0)
261{//Opening angle between two daughter tracks
262 Double_t mn[3] = {0,0,0};
263 Double_t mp[3] = {0,0,0};
264
265
266 esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
267 esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
268
269
270 fOpenAngle = TMath::ACos((mp[0]*mn[0] + mp[1]*mn[1] + mp[2]*mn[2])/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1] + mp[2]*mp[2])*TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1] + mn[2]*mn[2])));
271
272 return fOpenAngle;
273}
274
275//_________________________________________________
276Float_t AliTRDv0Info::PsiPair(AliESDv0 *esdv0)
277{//Angle between daughter momentum plane and plane perpendicular to magnetic field
278 Double_t x, y, z;
279 esdv0->GetXYZ(x,y,z);//Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
280
281 Double_t mn[3] = {0,0,0};
282 Double_t mp[3] = {0,0,0};
283
284
285 esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
286 esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
287
288
289 Double_t deltat = 1.;
290 deltat = TMath::ATan(mp[2]/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1])+1.e-13)) - TMath::ATan(mn[2]/(TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
291
292 Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
293
294 Double_t momPosProp[3];
295 Double_t momNegProp[3];
296
297 AliExternalTrackParam nt(*fTrackN), pt(*fTrackP);
298
299 fPsiPair = 4.;
300
301 if(nt.PropagateTo(radiussum,fMagField) == 0)//propagate tracks to the outside
302 fPsiPair = -5.;
303 if(pt.PropagateTo(radiussum,fMagField) == 0)
304 fPsiPair = -5.;
305 pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
306 nt.GetPxPyPz(momNegProp);
307
308 Double_t pEle =
309 TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
310 Double_t pPos =
311 TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
312
313 Double_t scalarproduct =
314 momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
315
316 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
317
318 fPsiPair = TMath::Abs(TMath::ASin(deltat/chipair));
319
320 return fPsiPair;
321
322}
323
324//_________________________________________________
325void AliTRDv0Info::V0fromTrack(AliTRDtrackInfo * const track, Int_t ivertex)
326{//Checks if track is a secondary vertex daughter (due to V0 finder)
327
328 fMagField = fESD->GetMagneticField();
329
330 fTrackID = track->GetTrackId();//index of the track
331
332 fTrack = fESD->GetTrack(fTrackID);//sets track information
333
334 fHasV0 = 0;
335
336 //comparing index of track with indices of pos./neg. V0 daughter :
337 AliESDv0 * esdv0 = fESD->GetV0(ivertex);
338 if((esdv0->GetIndex(0) == fTrackID)||(esdv0->GetIndex(1) == fTrackID))
339 {
340 fHasV0 = 1;//track belongs to vertex found by V0 finder!
341 fNindex = esdv0->GetIndex(0);
342 fPindex = esdv0->GetIndex(1);
343 fTrackN = fESD->GetTrack(esdv0->GetIndex(0));//providing information about the other of the two daughter tracks
344 fTrackP = fESD->GetTrack(esdv0->GetIndex(1));
345 GetESDv0Info(esdv0);//gets all the relevant information about our V0
346 }
347}
348//_________________________________________________
349void AliTRDv0Info::GetDetectorPID()
350{//PID likelihoods from TPC, TOF, and ITS, for all particle species
351
352 fTrackN->GetTPCpid(fDetPID[kNeg][kTPC]);
353 fTrackP->GetTPCpid(fDetPID[kPos][kTPC]);
354 fTrackN->GetTOFpid(fDetPID[kNeg][kTOF]);
355 fTrackP->GetTOFpid(fDetPID[kPos][kTOF]);
356 fTrackN->GetITSpid(fDetPID[kNeg][kITS]);
357 fTrackP->GetITSpid(fDetPID[kPos][kITS]);
358
359}
360
361//_________________________________________________
362Float_t AliTRDv0Info::Radius(AliESDv0 *esdv0)
363{//distance from secondary vertex to primary vertex in x-y plane
364 Double_t x, y, z;
365 esdv0->GetXYZ(x,y,z); //Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
366 fRadius = TMath::Sqrt(x*x + y*y);
367 return fRadius;
368
369}
370
371//_________________________________________________
372Int_t AliTRDv0Info::Quality(AliESDv0 *const esdv0)
373{
374 //
375 // Checking track and V0 quality status in order to exclude vertices based on poor information
376 //
377
378 Float_t nClsN;
379 nClsN = fTrackN->GetTPCNcls();//number of found clusters in TPC for negative track
380 Float_t nClsFN;
381 nClsFN = fTrackN->GetTPCNclsF();//number of findable clusters in TPC for negative track
382 Float_t nClsP;
383 nClsP = fTrackP->GetTPCNcls();//number of found clusters in TPC for positive track
384 Float_t nClsFP;
385 nClsFP = fTrackP->GetTPCNclsF();//number of findable clusters in TPC for positive track
386
387 fQuality = 0;
388
389
390 Float_t clsRatioN;
391 Float_t clsRatioP;
392
393 if((nClsFN == 0) || (nClsFP == 0))
394 return 2;
395
396 clsRatioN = nClsN/nClsFN; //ratios of found to findable clusters in TPC
397 clsRatioP = nClsP/nClsFP;
398
399 if (!(esdv0->GetOnFlyStatus()))//accept only vertices from online V0 finder
400 return 3;
401 if (!((fTrackP->GetStatus() &
402 AliESDtrack::kTPCrefit)))//accept only vertices in which both tracks have TPC refit
403 return 4;
404 if (!((fTrackN->GetStatus() &
405 AliESDtrack::kTPCrefit)))
406 return 5;
407 if (fTrackP->GetKinkIndex(0)>0 ||
408 fTrackN->GetKinkIndex(0)>0 )//exclude tracks with kinks
409 return 6;
410 if((clsRatioN < 0.6)||(clsRatioP < 0.6))//exclude tracks with low ratio of found to findable TPC clusters
411 return 7;
412 fQuality = 1;
413 return fQuality;
414}
415//_________________________________________________
416Bool_t AliTRDv0Info::GetV0PID(Int_t ipart, AliTRDtrackInfo *track)
417{//decides if track is accepted for one of the reference data samples
418
419 Int_t iDecay = -1;
420
421 //decide which decay has to be considered for which particle sample (Anti-Lambda will be treated separately)
422 if(ipart == AliPID::kElectron)
423 iDecay = kGamma;
424 else if(ipart == AliPID::kPion)
425 iDecay = kK0s;
426 else if(ipart == AliPID::kProton)
427 iDecay = kLambda;
428
429 Int_t iPSlot;//Mother momentum slots above/below 2.5 GeV
430
431
432 Bool_t pid = 0;//return value for V0 PID decision
433
434 if(!(track))
435 {
436 AliError("AliTRDv0Info::GetV0PID(Int_t ipart, AliTRDtrackInfo *track) : No track info found.\n");
437 return 0;
438 }
439
440 AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
441 if(!esdH)
442 {
443 AliError("AliTRDv0Info::GetV0PID(Int_t ipart, AliTRDtrackInfo *track) : ERROR - ESD input handler not found");
444 return 0;
445 }
446
447
448 fESD = esdH->GetEvent();
449
450 for(Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++)
451 {
452
453 if(pid == 0)
454 {
455 V0fromTrack(track, ivertex);//Get the V0 corresponding to the track (if there is a V0)
456
457 if(fV0Momentum > 2.5)//divide into slots according to reconstructed momentum of the mother particle
458 {iPSlot = 1;}
459 else
460 {iPSlot = 0;}
461 //Accept track for a sample only if...
462
463 if(!(fHasV0))//... there is a V0 found for it
464 continue;
465 if(!(fQuality == 1))//... it fulfills our quality criteria
466 continue;
467 if(!(fDCA < fUpDCA[iDecay]))//... distance of closest approach between daughters is reasonably small
468 continue;
469 if(!(fPointingAngle < fUpPointingAngle[iDecay]))//... pointing angle between momentum of mother particle and vector from prim. to sec. vertex is small
470 continue;
471 if(!(fRadius > fDownRadius[iDecay]))//... x-y plane distance of decay point to prim. vertex is bigger than a certain minimum value (for conversions)
472 continue;
473 if(!(fOpenAngle < fUpOpenAngle[iDecay]))//... opening angle is close enough to zero (for conversions)
474 continue;
475 if(!(TMath::Abs(fPsiPair) < fUpPsiPair[iDecay]))//... Psi-pair angle is close enough to zero(for conversions)
476 continue;
477
478 //specific cut criteria :
479 if(ipart == AliPID::kProton)
480 {//for proton sample: separate treatment of Lamba and Anti-Lambda decays:
481 //for Anti-Lambda:
482 //TPC PID likelihoods high enough for pi+ and anti-proton ; invariant mass calculated postulating these two particle species...
483 if((fDetPID[kNeg][kTPC][AliPID::kProton] > fDownTPCPIDneg[AliPID::kProton]) && (fDetPID[kPos][kTPC][AliPID::kPion] > fDownTPCPIDpos[AliPID::kPion]))
484 {
485 if(fNindex == fTrackID)
486 {
487 if((fInvMass[kAntiLambda] < fUpInvMass[kAntiLambda][iPSlot]) && (fInvMass[kAntiLambda] > fDownInvMass[kAntiLambda]))
488 {
489 pid = 1;
490 }
491 }
492 }
493 //for Lambda:
494 //TPC PID likelihoods high enough for pi- and proton ; invariant mass calculated accordingly
495 if((fDetPID[kNeg][kTPC][AliPID::kPion] > fDownTPCPIDneg[AliPID::kPion]) && (fDetPID[kPos][kTPC][AliPID::kProton] > fDownTPCPIDpos[AliPID::kProton]))
496 {
497 if(fPindex == fTrackID)
498 {
499 if((fInvMass[kLambda] < fUpInvMass[kLambda][iPSlot]) && (fInvMass[kLambda] > fDownInvMass[kLambda]))
500 {
501 pid = 1;
502 }
503 }
504 }
505 }
506 //for photon and K0s decays: equal TPC PID likelihood criteria for both daughters ; invariant mass calculated postulating two electrons/two pions
507 else
508 if((fDetPID[kNeg][kTPC][ipart] > fDownTPCPIDneg[ipart]) && (fDetPID[kPos][kTPC][ipart] > fDownTPCPIDpos[ipart]))
509 {
510 if((fInvMass[iDecay] < fUpInvMass[iDecay][iPSlot]) && (fInvMass[iDecay] > fDownInvMass[iDecay]))
511 {
512 pid = 1;
513 }
514 }
515 }
516 }
517
518 return pid;
519
520}
521//_________________________________________________
522void AliTRDv0Info::Print(Option_t */*opt*/) const
523{
524
525}