1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /* $Id: AliTRDv0Info.cxx 27496 2008-07-22 08:35:45Z cblume $ */
18 ////////////////////////////////////////////////////////////////////////////
20 // Reconstruction QA //
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 //
28 // (M.Heide, 2009/10/06) //
31 // Alex Bercuci <A.Bercuci@gsi.de> //
32 // Alex Wilk <wilka@uni-muenster.de> //
33 // Markus Heide <mheide@uni-muenster.de> //
35 ////////////////////////////////////////////////////////////////////////////
38 #include "TDatabasePDG.h"
40 #include "AliESDtrack.h"
44 #include "AliKFParticle.h"
45 #include "AliKFVertex.h"
47 #include "AliTRDv0Info.h"
48 #include "AliTRDtrackInfo.h"
49 #include "AliTRDtrackInfo.h"
51 ClassImp(AliTRDv0Info)
53 //_________________________________________________
54 AliTRDv0Info::AliTRDv0Info()
72 // Default constructor
75 memset(fDetPID, 0, 2*kNDaughters*kNDetectors*AliPID::kSPECIES*sizeof(Float_t));
76 memset(fComPID, 0, 2*kNDaughters*AliPID::kSPECIES*sizeof(Float_t));
77 memset(fInvMass, 0, kNDecays*sizeof(Double_t));
78 memset(fArmenteros, 0, kNDecays*sizeof(Bool_t));
79 memset(fTPCdEdx, 0, kNDaughters*sizeof(Float_t));
80 memset(fChi2ndf, 0, kNDecays*sizeof(Double_t));
82 /////////////////////////////////////////////////////////////////////////////
83 //Set Cut values: First specify decay in brackets, then the actual cut value!
84 /////////////////////////////////////////////////////////////////////////////
86 //Upper limit for distance of closest approach of two daughter tracks :
87 fUpDCA[kGamma] = 1000.;
89 fUpDCA[kLambda] = 0.2;
90 fUpDCA[kAntiLambda] = 0.2;
92 //Upper limit for pointing angle (= angle between between vector from primary to secondary vertex and reconstructed momentum of V0 mother particle) :
93 fUpPointingAngle[kGamma] = 0.03;
94 fUpPointingAngle[kK0s] = 0.03;
95 fUpPointingAngle[kLambda] = 0.04;
96 fUpPointingAngle[kAntiLambda] = 0.04;
98 //Upper limit for invariant mass of V0 mother :
99 fUpInvMass[kGamma][0] = 0.05;// second pair of brackets is for momentum bin: 0: below mother momentm of 2.5 GeV
100 fUpInvMass[kGamma][1] = 0.07;//1: above 2.5 GeV
101 fUpInvMass[kK0s][0] = fUpInvMass[kK0s][1] = 0.50265;
102 fUpInvMass[kLambda][0] = fUpInvMass[kLambda][1] = 1.1207;
103 fUpInvMass[kAntiLambda][0] = fUpInvMass[kAntiLambda][1] = 1.1207;
105 //Lower limit for invariant mass of V0 mother :
106 fDownInvMass[kGamma] = -1.;
107 fDownInvMass[kK0s] = 0.49265;
108 fDownInvMass[kLambda] = 1.107;
109 fDownInvMass[kAntiLambda] = 1.107;
111 //Upper limit for KF Chi2/NDF value;
112 fUpChi2ndf[kGamma] = 10000.;//7.;
113 fUpChi2ndf[kK0s] = 10000.;//5.;
114 fUpChi2ndf[kLambda] = 10000.;//5.;
115 fUpChi2ndf[kAntiLambda] = 10000.;//5.;
117 //Lower limit for distance from secondary vertex to primary vertex in x-y plane :
118 fDownRadius[kGamma] = 6.;
119 fDownRadius[kK0s] = 0.;
120 fDownRadius[kLambda] = 0.;
121 fDownRadius[kAntiLambda] = 0.;
123 //Upper limit for distance from secondary vertex to primary vertex in x-y plane :
124 fUpRadius[kGamma] = 1000.;
125 fUpRadius[kK0s] = 20.;
126 fUpRadius[kLambda] = 1000.;
127 fUpRadius[kAntiLambda] = 1000.;
129 //Upper limit for opening angle between two daughter tracks (characteristically near zero for conversions) :
130 fUpOpenAngle[kGamma] = 0.1;
131 fUpOpenAngle[kK0s] = 3.15;
132 fUpOpenAngle[kLambda] = 3.15;
133 fUpOpenAngle[kAntiLambda] = 3.15;
135 //Upper limit for angle between daughter momentum plane and plane perpendicular to magnetic field (characteristically around zero for conversions) :
136 fUpPsiPair[kGamma] = 0.05;
137 fUpPsiPair[kK0s] = 1.6;
138 fUpPsiPair[kLambda] = 1.6;
139 fUpPsiPair[kAntiLambda] = 1.6;
141 //Lower limit for likelihood value of TPC PID :
142 fDownTPCPIDneg[AliPID::kElectron] = 0.;
143 fDownTPCPIDpos[AliPID::kElectron] = 0.;
145 fDownTPCPIDneg[AliPID::kMuon] = 0.;
146 fDownTPCPIDpos[AliPID::kMuon] = 0.;
148 fDownTPCPIDneg[AliPID::kPion] = 0.;
149 fDownTPCPIDpos[AliPID::kPion] = 0.;
151 fDownTPCPIDneg[AliPID::kKaon] = 0.;
152 fDownTPCPIDpos[AliPID::kKaon] = 0.;
154 fDownTPCPIDneg[AliPID::kProton] = 0.;
155 fDownTPCPIDpos[AliPID::kProton] = 0.;
157 //Lower limit for likelihood value of combined PID :
158 fDownComPIDneg[AliPID::kElectron] = 0.;
159 fDownComPIDpos[AliPID::kElectron] = 0.;
161 fDownComPIDneg[AliPID::kMuon] = 0.;
162 fDownComPIDpos[AliPID::kMuon] = 0.;
164 fDownComPIDneg[AliPID::kPion] = 0.;
165 fDownComPIDpos[AliPID::kPion] = 0.;
167 fDownComPIDneg[AliPID::kKaon] = 0.;
168 fDownComPIDpos[AliPID::kKaon] = 0.;
170 fDownComPIDneg[AliPID::kProton] = 0.;
171 fDownComPIDpos[AliPID::kProton] = 0.;
173 //Lower limit for likelihood value of combined PID for daughter track which doesn't enter reference data (here: pion daughters from Lambda decays:
174 fDownComPIDnegPart[AliPID::kElectron] = 0.;
175 fDownComPIDposPart[AliPID::kElectron] = 0.;
177 fDownComPIDnegPart[AliPID::kMuon] = 0.;
178 fDownComPIDposPart[AliPID::kMuon] = 0.;
180 fDownComPIDnegPart[AliPID::kPion] = 0.;
181 fDownComPIDposPart[AliPID::kPion] = 0.;
183 fDownComPIDnegPart[AliPID::kKaon] = 0.;
184 fDownComPIDposPart[AliPID::kKaon] = 0.;
186 fDownComPIDnegPart[AliPID::kProton] = 0.;
187 fDownComPIDposPart[AliPID::kProton] = 0.;
189 //Parameters for data with well-calibrated PID (after usage of tender):
190 /* //Lower limit for likelihood value of TPC PID :
191 fDownTPCPIDneg[AliPID::kElectron] = 0.21;
192 fDownTPCPIDpos[AliPID::kElectron] = 0.21;
194 fDownTPCPIDneg[AliPID::kMuon] = 0.21;
195 fDownTPCPIDpos[AliPID::kMuon] = 0.21;
197 fDownTPCPIDneg[AliPID::kPion] = 0.21;
198 fDownTPCPIDpos[AliPID::kPion] = 0.21;
200 fDownTPCPIDneg[AliPID::kKaon] = 0.21;
201 fDownTPCPIDpos[AliPID::kKaon] = 0.21;
203 fDownTPCPIDneg[AliPID::kProton] = 0.21;
204 fDownTPCPIDpos[AliPID::kProton] = 0.21;
206 //Lower limit for likelihood value of combined PID :
207 fDownComPIDneg[AliPID::kElectron] = 0.21;
208 fDownComPIDpos[AliPID::kElectron] = 0.21;
210 fDownComPIDneg[AliPID::kMuon] = 0.21;
211 fDownComPIDpos[AliPID::kMuon] = 0.21;
213 fDownComPIDneg[AliPID::kPion] = 0.9;
214 fDownComPIDpos[AliPID::kPion] = 0.9;
216 fDownComPIDneg[AliPID::kKaon] = 0.21;
217 fDownComPIDpos[AliPID::kKaon] = 0.21;
219 fDownComPIDneg[AliPID::kProton] = 0.9;
220 fDownComPIDpos[AliPID::kProton] = 0.9;
222 //Lower limit for likelihood value of combined PID for daughter track which doesn't enter reference data (here: pion daughters from Lambda decays:
223 fDownComPIDnegPart[AliPID::kElectron] = 0.05;
224 fDownComPIDposPart[AliPID::kElectron] = 0.05;
226 fDownComPIDnegPart[AliPID::kMuon] = 0.05;
227 fDownComPIDposPart[AliPID::kMuon] = 0.05;
229 fDownComPIDnegPart[AliPID::kPion] = 0.05;
230 fDownComPIDposPart[AliPID::kPion] = 0.05;
232 fDownComPIDnegPart[AliPID::kKaon] = 0.05;
233 fDownComPIDposPart[AliPID::kKaon] = 0.05;
235 fDownComPIDnegPart[AliPID::kProton] = 0.05;
236 fDownComPIDposPart[AliPID::kProton] = 0.05;*/
239 //_________________________________________________
240 AliTRDv0Info::AliTRDv0Info(const AliTRDv0Info &ref)
242 ,fQuality(ref.fQuality)
244 ,fPointingAngle(ref.fPointingAngle)
245 ,fOpenAngle(ref.fOpenAngle)
246 ,fPsiPair(ref.fPsiPair)
247 ,fMagField(ref.fMagField)
248 ,fRadius(ref.fRadius)
249 ,fV0Momentum(ref.fV0Momentum)
250 ,fNindex(ref.fNindex)
251 ,fPindex(ref.fPindex)
252 ,fInputEvent(ref.fInputEvent)
253 ,fPrimaryVertex(ref.fPrimaryVertex)
254 ,fTrackP(ref.fTrackP)
255 ,fTrackN(ref.fTrackN)
261 memcpy(fDetPID, ref.fDetPID, 2*kNDaughters*kNDetectors*AliPID::kSPECIES*sizeof(Float_t));
262 memcpy(fComPID, ref.fComPID, 2*kNDaughters*AliPID::kSPECIES*sizeof(Float_t));
263 memcpy(fInvMass, ref.fInvMass, kNDecays*sizeof(Double_t));
264 memcpy(fArmenteros, ref.fArmenteros, kNDecays*sizeof(Bool_t));
265 memcpy(fChi2ndf, ref.fChi2ndf, kNDecays*sizeof(Double_t));
266 memcpy(fTPCdEdx, ref.fTPCdEdx, kNDaughters*sizeof(Float_t));
268 //Upper limit for distance of closest approach of two daughter tracks :
269 memcpy(fUpDCA, ref.fUpDCA, kNDecays*sizeof(Float_t));
270 memcpy(fUpPointingAngle, ref.fUpPointingAngle, kNDecays*sizeof(Float_t));
271 memcpy(fUpOpenAngle, ref.fUpOpenAngle, kNDecays*sizeof(Float_t));
272 memcpy(fDownOpenAngle, ref.fDownOpenAngle, kNDecays*sizeof(Float_t));
273 memcpy(fUpPsiPair, ref.fUpPsiPair, kNDecays*sizeof(Float_t));
274 memcpy(fDownPsiPair, ref.fDownPsiPair, kNDecays*sizeof(Float_t));
275 memcpy(fUpInvMass, ref.fUpInvMass, kNDecays*kNMomBins*sizeof(Double_t));
276 memcpy(fDownInvMass, ref.fDownInvMass, kNDecays*sizeof(Double_t));
277 memcpy(fUpChi2ndf, ref.fUpChi2ndf, kNDecays*sizeof(Double_t));
278 memcpy(fUpRadius, ref.fUpRadius, kNDecays*sizeof(Float_t));
279 memcpy(fDownRadius, ref.fDownRadius, kNDecays*sizeof(Float_t));
280 memcpy(fDownTPCPIDneg, ref.fDownTPCPIDneg, AliPID::kSPECIES*sizeof(Float_t));
281 memcpy(fDownTPCPIDpos, ref.fDownTPCPIDpos, AliPID::kSPECIES*sizeof(Float_t));
282 memcpy(fDownComPIDneg, ref.fDownComPIDneg, AliPID::kSPECIES*sizeof(Float_t));
283 memcpy(fDownComPIDpos, ref.fDownComPIDpos, AliPID::kSPECIES*sizeof(Float_t));
284 memcpy(fDownComPIDnegPart, ref.fDownComPIDnegPart, AliPID::kSPECIES*sizeof(Float_t));
285 memcpy(fDownComPIDposPart, ref.fDownComPIDposPart, AliPID::kSPECIES*sizeof(Float_t));
288 //_________________________________________________
289 void AliTRDv0Info::SetV0Info(const AliESDv0 *esdv0)
291 //Gets values of ESDv0 and daughter track properties
292 //See header file for description of variables
294 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)
296 fRadius = Radius(esdv0);//distance from secondary vertex to primary vertex in x-y plane
298 fDCA = esdv0->GetDcaV0Daughters();//distance of closest approach of two daughter tracks
300 fPointingAngle = TMath::ACos(esdv0->GetV0CosineOfPointingAngle());// pointing angle (= angle between between vector from primary to secondary vertex and reconstructed momentum of V0 mother particle)
302 fOpenAngle = OpenAngle(esdv0);//Opening angle between two daughter tracks
304 fPsiPair = PsiPair(esdv0);//Angle between daughter momentum plane and plane perpendicular to magnetic field
306 fV0Momentum = V0Momentum(esdv0);//Reconstructed momentum of the mother particle
308 //4 decay types : conversions, K0s, Lambda, Anti-Lambda
309 //five particle types: electrons, muons, pions, kaons, protons (muons and kaons not involved)
310 for(Int_t idecay(0), part1(-1), part2(-1); idecay < kNDecays; idecay++){
312 fArmenteros[idecay]=Armenteros(esdv0, idecay);//Attribute the Armenteros yes/no decision for every decay type
313 if(idecay == kLambda){ //protons and pions from Lambda
314 part1 = AliPID::kProton;
315 part2 = AliPID::kPion;
316 } else if(idecay == kAntiLambda) { //antiprotons and pions from Anti-Lambda
317 part1 = AliPID::kPion;
318 part2 = AliPID::kProton;
319 } else if(idecay == kK0s) {//pions from K0s
320 part1 = part2 = AliPID::kPion;
321 } else if(idecay == kGamma) {//electrons from conversions
322 part1 = part2 = AliPID::kElectron;
324 fInvMass[idecay] = InvMass(part1, part2, esdv0);//Calculate invariant mass for all of our four supposed decays
326 // Comment out until bug fix is provided
327 // A.Bercuci 14. July 2010
328 //fChi2ndf[idecay] = KFChi2ndf(part1, part2,idecay);
331 //Gets all likelihood values from TPC, TOF and ITS PID for the fDetPID[kNDaughters][kNDetectors][AliPID::kSPECIES] array
333 //Bayesian combination of likelihoods from TPC and TOF
335 //TPC dE/dx values for both tracks
339 //_________________________________________________
340 Float_t AliTRDv0Info::V0Momentum(const AliESDv0 *esdv0) const
343 // Reconstructed momentum of V0 mother particle
346 Double_t mn[3] = {0,0,0};
347 Double_t mp[3] = {0,0,0};
350 esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
351 esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
354 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]));
357 //_________________________________________________
358 Double_t AliTRDv0Info::InvMass(Int_t part1, Int_t part2, const AliESDv0 *esdv0) const
361 // Invariant mass of reconstructed V0 mother
364 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)};
365 //Masses of electrons, muons, pions, kaons and protons, as implemented in ROOT
368 Double_t mn[3] = {0,0,0};
369 Double_t mp[3] = {0,0,0};
371 esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
372 esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
374 Double_t mass1 = kpmass[part1];//sets supposed rest masses for both daughters: positive
375 Double_t mass2 = kpmass[part2];//negative
377 //Calculate daughters' energies :
378 Double_t e1 = TMath::Sqrt(mass1*mass1+
382 Double_t e2 = TMath::Sqrt(mass2*mass2+
387 //Sum of daughter momenta :
389 (mn[0]+mp[0])*(mn[0]+mp[0])+
390 (mn[1]+mp[1])*(mn[1]+mp[1])+
391 (mn[2]+mp[2])*(mn[2]+mp[2]);
394 Double_t mInv = TMath::Sqrt((e1+e2)*(e1+e2)-momsum);
399 //_________________________________________________
400 Float_t AliTRDv0Info::OpenAngle(const AliESDv0 *esdv0)
402 //Opening angle between two daughter tracks
403 Double_t mn[3] = {0,0,0};
404 Double_t mp[3] = {0,0,0};
407 esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
408 esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
411 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])));
416 //_________________________________________________
417 Float_t AliTRDv0Info::PsiPair(const AliESDv0 *esdv0)
419 //Angle between daughter momentum plane and plane perpendicular to magnetic field
421 esdv0->GetXYZ(x,y,z);//Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
423 Double_t mn[3] = {0,0,0};
424 Double_t mp[3] = {0,0,0};
427 esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
428 esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
431 Double_t deltat = 1.;
432 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
434 Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
436 Double_t momPosProp[3];
437 Double_t momNegProp[3];
439 AliExternalTrackParam nt(*fTrackN), pt(*fTrackP);
443 if(nt.PropagateTo(radiussum,fMagField) == 0)//propagate tracks to the outside
445 if(pt.PropagateTo(radiussum,fMagField) == 0)
447 pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
448 nt.GetPxPyPz(momNegProp);
451 TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
453 TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
455 Double_t scalarproduct =
456 momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
458 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
460 fPsiPair = TMath::Abs(TMath::ASin(deltat/chipair));
465 //_________________________________________________
466 Double_t AliTRDv0Info::KFChi2ndf(Int_t part1, Int_t part2,Int_t decay){
467 //Calculates Kalman filter Chi2/NDF
468 Int_t mothers[4]={22,310,3122,3122};
470 const Double_t partMass=TDatabasePDG::Instance()->GetParticle(mothers[decay])->Mass();
471 const Double_t massWidth[4] = {0.001, 0., 0., 0.};
473 AliKFParticle *kfMother = CreateMotherParticle(fTrackP, fTrackN, part1, part2);
480 // production vertex is set in the 'CreateMotherParticle' function
481 kfMother->SetMassConstraint(partMass, massWidth[decay]);
483 Double_t chi2ndf = (kfMother->GetChi2()/kfMother->GetNDF());
485 if(kfMother)delete kfMother;
488 //________________________________________________________________
489 AliKFParticle *AliTRDv0Info::CreateMotherParticle(const AliESDtrack *pdaughter, const AliESDtrack *ndaughter, Int_t pspec, Int_t nspec){
491 // Creates a mother particle
493 AliKFParticle pkfdaughter(*pdaughter, pspec);
494 AliKFParticle nkfdaughter(*ndaughter, nspec);
497 // Create the mother particle
498 AliKFParticle *m = new AliKFParticle(pkfdaughter, nkfdaughter);
500 AliKFVertex improvedVertex = *fPrimaryVertex;
501 improvedVertex += *m;
502 m->SetProductionVertex(improvedVertex);
507 //_________________________________________________
508 Int_t AliTRDv0Info::HasTrack(const AliTRDtrackInfo * const track) const
510 //Checks if track is a secondary vertex daughter (due to V0 finder)
513 if(!fTrackP->GetID()) return 0;
514 if(!fTrackN->GetID()) return 0;
516 Int_t trackID(track->GetTrackId());//index of the track
517 return HasTrack(trackID);
520 //_________________________________________________
521 Int_t AliTRDv0Info::HasTrack(Int_t trackID) const
523 //comparing index of track with indices of pos./neg. V0 daughter :
524 if(fNindex==trackID) return -1;
525 else if(fPindex==trackID) return 1;
529 //_________________________________________________
530 void AliTRDv0Info::GetDetectorPID()
532 //PID likelihoods from TPC, TOF, and ITS, for all particle species
534 fTrackN->GetTPCpid(fDetPID[kNeg][kTPC]);
535 fTrackP->GetTPCpid(fDetPID[kPos][kTPC]);
536 fTrackN->GetTOFpid(fDetPID[kNeg][kTOF]);
537 fTrackP->GetTOFpid(fDetPID[kPos][kTOF]);
538 fTrackN->GetITSpid(fDetPID[kNeg][kITS]);
539 fTrackP->GetITSpid(fDetPID[kPos][kITS]);
541 Long_t statusN = fTrackN->GetStatus();
542 Long_t statusP = fTrackP->GetStatus();
544 if(!(statusN & AliESDtrack::kTPCpid)){
545 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
546 fDetPID[kNeg][kTPC][iPart] = 0.2;
549 if(!(statusN & AliESDtrack::kTOFpid)){
550 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
551 fDetPID[kNeg][kTOF][iPart] = 0.2;
555 if(!(statusN & AliESDtrack::kITSpid)){
556 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
557 fDetPID[kNeg][kITS][iPart] = 0.2;
560 if(!(statusP & AliESDtrack::kTPCpid)){
561 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
562 fDetPID[kPos][kTPC][iPart] = 0.2;
565 if(!(statusP & AliESDtrack::kTOFpid)){
566 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
567 fDetPID[kPos][kTOF][iPart] = 0.2;
571 if(!(statusP & AliESDtrack::kITSpid)){
572 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++){
573 fDetPID[kPos][kITS][iPart] = 0.2;
578 //____________________________________________________________________________________
579 void AliTRDv0Info::CombinePID()
581 //combined bayesian PID from TPC and TOF
582 Double_t partrat[AliPID::kSPECIES] = {0.208, 0.010, 0.662, 0.019, 0.101};
584 for(Int_t iSign = 0; iSign < kNDaughters; iSign++)
586 for(Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++)
588 fComPID[iSign][iPart] = (partrat[iPart]*fDetPID[iSign][kTPC][iPart]*fDetPID[iSign][kTOF][iPart])/((partrat[0]*fDetPID[iSign][kTPC][0]*fDetPID[iSign][kTOF][0])+(partrat[1]*fDetPID[iSign][kTPC][1]*fDetPID[iSign][kTOF][1])+(partrat[2]*fDetPID[iSign][kTPC][2]*fDetPID[iSign][kTOF][2])+(partrat[3]*fDetPID[iSign][kTPC][3]*fDetPID[iSign][kTOF][3])+(partrat[4]*fDetPID[iSign][kTPC][4]*fDetPID[iSign][kTOF][4]));
593 //_________________________________________________
594 Bool_t AliTRDv0Info::GetTPCdEdx()
596 //gets the TPC dE/dx for both daughter tracks
597 if(!fTrackP->GetID()) return 0;
598 if(!fTrackN->GetID()) return 0;
600 fTPCdEdx[kNeg] = fTrackN->GetTPCsignal();
601 fTPCdEdx[kPos] = fTrackP->GetTPCsignal();
605 //_________________________________________________
606 Bool_t AliTRDv0Info::TPCdEdxCuts(Int_t part, const AliTRDtrackInfo * const track)
608 //applies cuts on TPC dE/dx according to particle species; cutting lines are drawn shifted to the Bethe-Bloch paremeterization
609 if(!fTrackP->GetID()) return 0;
610 if(!fTrackN->GetID()) return 0;
613 Double_t alephParameters[5];
616 alephParameters[0] = 0.0283086;
617 alephParameters[1] = 2.63394e+01;
618 alephParameters[2] = 5.04114e-11;
619 alephParameters[3] = 2.12543e+00;
620 alephParameters[4] = 4.88663e+00;
623 Double_t deposit = 0;
625 if(HasTrack(track) == 1){
627 deposit = fTPCdEdx[kPos];
629 else if(HasTrack(track) == -1){
631 deposit = fTPCdEdx[kNeg];
634 printf("No track found");
639 Float_t upLimits[5]={85,1000,50*AliExternalTrackParam::BetheBlochAleph(x/0.13957, alephParameters[0], alephParameters[1], alephParameters[2], alephParameters[3], alephParameters[4])+6,1000,50*AliExternalTrackParam::BetheBlochAleph(x/0.93827, alephParameters[0], alephParameters[1], alephParameters[2], alephParameters[3], alephParameters[4])+10};
640 Float_t downLimits[5]={62,40,50*AliExternalTrackParam::BetheBlochAleph(x/0.13957, alephParameters[0], alephParameters[1], alephParameters[2], alephParameters[3], alephParameters[4])-6,40,50*AliExternalTrackParam::BetheBlochAleph(x/0.93827, alephParameters[0], alephParameters[1], alephParameters[2], alephParameters[3], alephParameters[4])-11};
654 if(deposit < downLimits[part])
656 if(deposit > upLimits[part])
663 //_________________________________________________
664 Float_t AliTRDv0Info::Radius(const AliESDv0 *esdv0)
666 //distance from secondary vertex to primary vertex in x-y plane
668 esdv0->GetXYZ(x,y,z); //Reconstructed coordinates of V0
669 fRadius = TMath::Sqrt(x*x + y*y);
674 //_________________________________________________
675 Int_t AliTRDv0Info::Quality(const AliESDv0 *const esdv0)
678 // Checking track and V0 quality status in order to exclude vertices based on poor information
682 nClsN = fTrackN->GetTPCNcls();//number of found clusters in TPC for negative track
684 nClsFN = fTrackN->GetTPCNclsF();//number of findable clusters in TPC for negative track
686 nClsP = fTrackP->GetTPCNcls();//number of found clusters in TPC for positive track
688 nClsFP = fTrackP->GetTPCNclsF();//number of findable clusters in TPC for positive track
693 if (!(esdv0->GetOnFlyStatus()))//accept only vertices from online V0 finder
699 if((nClsFN < 80) || (nClsFP < 80)) return -2;//reject all V0s where at least one track has less than 80 TPC clusters
701 // Chi2 per TPC cluster
702 Int_t nTPCclustersP = fTrackP->GetTPCclusters(0);
703 Int_t nTPCclustersN = fTrackN->GetTPCclusters(0);
704 Float_t chi2perTPCclusterP = fTrackP->GetTPCchi2()/Float_t(nTPCclustersP);
705 Float_t chi2perTPCclusterN = fTrackN->GetTPCchi2()/Float_t(nTPCclustersN);
707 if((chi2perTPCclusterN > 3.5)||(chi2perTPCclusterP > 3.5)) return -3;//reject all V0s where at least one track has a chi2 above 3.5
709 clsRatioN = nClsN/nClsFN; //ratios of found to findable clusters in TPC
710 clsRatioP = nClsP/nClsFP;
712 if((clsRatioN < 0.6)||(clsRatioP < 0.6))//exclude tracks with low ratio of found to findable TPC clusters
715 if (!((fTrackP->GetStatus() &
716 AliESDtrack::kTPCrefit)))//accept only vertices in which both tracks have TPC refit
718 if (!((fTrackN->GetStatus() &
719 AliESDtrack::kTPCrefit)))
721 if (fTrackP->GetKinkIndex(0)>0 ||
722 fTrackN->GetKinkIndex(0)>0 )//exclude tracks with kinks
730 //________________________________________________________________
731 Bool_t AliTRDv0Info::V0SignCheck(){
733 // Check if v0 daughters really carry opposite charges
736 Int_t qP = fTrackP->Charge();
737 Int_t qN = fTrackN->Charge();
739 if((qP*qN) != -1) return kFALSE;
743 //___________________________________________________________________
744 Bool_t AliTRDv0Info::Armenteros(const AliESDv0 *esdv0, Int_t decay){
746 // computes the Armenteros variables for given V0
748 Double_t mn[3] = {0,0,0};
749 Double_t mp[3] = {0,0,0};
750 Double_t mm[3] = {0,0,0};
753 esdv0->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
754 esdv0->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
757 esdv0->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
758 esdv0->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
760 esdv0->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
762 TVector3 vecN(mn[0],mn[1],mn[2]);
763 TVector3 vecP(mp[0],mp[1],mp[2]);
764 TVector3 vecM(mm[0],mm[1],mm[2]);
766 Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
767 Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
769 Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
770 ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
771 Double_t qt = vecP.Mag()*sin(thetaP);
777 Double_t lCutAP[2];//Lambda/Anti-Lambda cuts
780 const Double_t cutAlpha[2] = {0.35, 0.45}; // [0.35, 0.45]
781 const Double_t cutQT = 0.015;
782 if(TMath::Abs(ap[0]) > cutAlpha[0] && TMath::Abs(ap[0]) < cutAlpha[1]) return kFALSE;
784 if(ap[1] > cutQT) return kFALSE;
788 const Double_t cutQT = 0.1075;
789 const Double_t cutAP = 0.22 * TMath::Sqrt( TMath::Abs( (1-ap[0]*ap[0]/(0.92*0.92)) ) );
790 if(ap[1] < cutQT) return kFALSE;
791 if(ap[1] > cutAP) return kFALSE;
794 const Double_t cutQT = 0.03;
795 const Double_t cutAlpha = 0.7; // VERY strong - should supress the overlap with K0
796 lCutAP[0] = 1.0 - (ap[0]-0.7 * ap[0]-0.7)*1.1 - 0.87;
797 if(TMath::Abs(ap[0]) > cutAlpha) return kFALSE;
798 if(ap[1] < cutQT) return kFALSE;
799 if(ap[1] > lCutAP[0]) return kFALSE;
803 const Double_t cutQT = 0.03;
804 const Double_t cutAlpha = 0.7; // VERY strong - should supress the overlap with K0
805 lCutAP[1] = 1.0 - (ap[0]+0.7 * ap[0]+0.7)*1.1 - 0.87;
806 if(TMath::Abs(ap[0]) > cutAlpha) return kFALSE;
807 if(ap[1] < cutQT) return kFALSE;
808 if(ap[1] > lCutAP[1]) return kFALSE;
812 //_________________________________________________
813 Int_t AliTRDv0Info::GetPID(Int_t ipart, AliTRDtrackInfo *track)
815 // Decides if track is accepted for one of the reference data samples
818 AliError("No track info");
821 if(!HasTrack(track)){
822 AliDebug(2, "Track not attached to v0.");
826 //translate ipart to decay (Anti-Lambda will be treated separately)
829 case AliPID::kElectron: iDecay = kGamma; break;
830 case AliPID::kPion: iDecay = kK0s; break;
831 case AliPID::kProton: iDecay = kLambda; break;
833 AliDebug(1, Form("Hypothesis \"ipart=%d\" not handled", ipart));
837 //... it fulfills our quality criteria
838 if(!(fQuality == 1)) return -4;
839 //... distance of closest approach between daughters is reasonably small
840 if((fDCA > fUpDCA[iDecay])) return -5;
841 //... pointing angle between momentum of mother particle and vector from prim. to sec. vertex is small
842 if((fPointingAngle > fUpPointingAngle[iDecay])) return -6;
843 //... x-y plane distance of decay point to prim. vertex is bigger than a certain minimum value (for conversions)
844 if((fRadius < fDownRadius[iDecay])) return -7;
845 //...or smaller than a maximum value (for K0s)
846 if((fRadius > fUpRadius[iDecay])) return -8;
847 //... opening angle is close enough to zero (for conversions)
848 if((fOpenAngle > fUpOpenAngle[iDecay])) return -9;
849 //... Psi-pair angle is close enough to zero(for conversions)
850 if((TMath::Abs(fPsiPair) > fUpPsiPair[iDecay])) return -10;
854 //Mother momentum slots above/below 2.5 GeV
855 Int_t iPSlot(fV0Momentum > 2.5);
856 Int_t trackID(track->GetTrackId());
858 //specific cut criteria :
859 if(ipart == AliPID::kProton) {
860 if((fInvMass[kK0s] < fUpInvMass[kK0s][iPSlot]) && (fInvMass[kK0s] > fDownInvMass[kK0s])) return -11;//explicit exclusion of K0s decays
862 if(fOpenAngle < (0.3 - 0.2*fV0Momentum))return -9;
866 //for proton sample: separate treatment of Lamba and Anti-Lambda decays:
868 //Combined PID likelihoods high enough for pi+ and anti-proton ; invariant mass calculated postulating these two particle species...
869 //if((fComPID[kNeg][AliPID::kProton] > fDownComPIDneg[AliPID::kProton]) && (fComPID[kPos][AliPID::kPion] > fDownComPIDposPart[AliPID::kPion])) {
870 //if((fDetPID[kNeg][kTPC][AliPID::kProton] > fDownTPCPIDneg[AliPID::kProton]) && (fDetPID[kPos][kTPC][AliPID::kPion] > fDownTPCPIDpos[AliPID::kPion])){
871 if((TPCdEdxCuts(ipart, track))){//momentary solution: direct cut on TPC dE/dx
872 if(fNindex == trackID) {//we're only interested in the anti-proton
873 if(fArmenteros[kAntiLambda]){//Armenteros condition has to be fulfilled
874 if(fChi2ndf[kAntiLambda] < fUpChi2ndf[kAntiLambda]){//Kalman filter Chi2/NDF not allowed to be too large
875 if((fInvMass[kAntiLambda] < fUpInvMass[kAntiLambda][iPSlot]) && (fInvMass[kAntiLambda] > fDownInvMass[kAntiLambda])){
877 } else cutCode = -15;
886 //TPC PID likelihoods high enough for pi- and proton ; invariant mass calculated accordingly
887 //if((fComPID[kNeg][AliPID::kPion] > fDownComPIDnegPart[AliPID::kPion]) && (fComPID[kPos][AliPID::kProton] > fDownComPIDpos[AliPID::kProton])) {
888 //if((fDetPID[kNeg][kTPC][AliPID::kPion] > fDownTPCPIDneg[AliPID::kPion]) && (fDetPID[kPos][kTPC][AliPID::kProton] > fDownTPCPIDpos[AliPID::kProton])){
889 if((TPCdEdxCuts(ipart, track))){//momentary solution: direct TPC dE/dx cuts
890 if(fPindex == trackID) {
891 if(fArmenteros[kLambda]){
892 if(fChi2ndf[kLambda] < fUpChi2ndf[kLambda]){
893 if((fInvMass[kLambda] < fUpInvMass[kLambda][iPSlot]) && (fInvMass[kLambda] > fDownInvMass[kLambda])){
895 } else cutCode = -15;
906 //for K0s decays: equal TPC PID likelihood criteria for both daughters ; invariant mass calculated postulating two pions
907 if(ipart == AliPID::kPion) {
909 if(fOpenAngle < (1.0/(fV0Momentum + 0.3) - 0.1))
912 //explicit exclusion of Lambda decays
913 if((fInvMass[kLambda] < fUpInvMass[kLambda][iPSlot]) && (fInvMass[kLambda] > fDownInvMass[kLambda])) return -11;
914 //explicit exclusion of Anti-Lambda decays
915 if((fInvMass[kAntiLambda] < fUpInvMass[kAntiLambda][iPSlot]) && (fInvMass[kAntiLambda] > fDownInvMass[kAntiLambda])) return -11;
917 //if((fDetPID[kNeg][kTPC][ipart] < fDownTPCPIDneg[ipart]) || (fDetPID[kPos][kTPC][ipart] < fDownTPCPIDpos[ipart])) return -12;
918 if(!(TPCdEdxCuts(ipart, track))){//momentary solution: direct TPC dE/dx cuts
924 //for photon conversions: equal combined PID likelihood criteria for both daughters ; invariant mass calculated postulating two electrons
925 //No Lambda/K0s exclusion is provided, since these contributions hardly ever interfere with gamma invariant mass!
926 //Float_t momentum(track->GetESDinfo()->GetOuterParam()->P());
927 if(ipart == AliPID::kElectron) {
928 //if(momentum > 1.75) {//since combined PID performs a little worse in simulations than TPC standalone for higher momenta, ONLY TPC PID is used here
929 //if((fDetPID[kNeg][kTPC][ipart] < fDownTPCPIDneg[ipart]) || (fDetPID[kPos][kTPC][ipart] < fDownTPCPIDpos[ipart])) return -12;
930 //} else {//for low momenta, combined PID from TOF and TPC is used to get rid of proton contamination
931 //if((fComPID[kNeg][ipart] > fDownComPIDneg[ipart]) && (fComPID[kPos][ipart] > fDownComPIDpos[ipart])) return 1;
933 if(!(TPCdEdxCuts(ipart, track))){//momentary solution for direct TPC dE/dx cut
940 //Armenteros-Polanski cut
941 if(!(fArmenteros[iDecay])) return -13;
943 //Kalman filter Chi2/NDF cut
944 if(fChi2ndf[iDecay] > fUpChi2ndf[iDecay]) return -14;
946 //Invariant mass cut for K0s and photons, assuming two pions/two electrons as daughters:
948 if((fInvMass[iDecay] > fUpInvMass[iDecay][iPSlot]) || (fInvMass[iDecay] < fDownInvMass[iDecay])) {
957 //_________________________________________________
958 void AliTRDv0Info::Print(Option_t *opt) const
960 //prints text for debugging etc.
961 printf("V0 P[%d] N[%d]\n", fPindex, fNindex);
962 printf(" DCA[%5.3f] Radius[%5.3f]\n", fDCA, fRadius);
963 printf(" Angles : Pointing[%5.3f] Open[%5.3f] Psi[%5.3f]\n", fPointingAngle, fOpenAngle, fPsiPair);
964 if(strcmp(opt, "a")!=0) return;
965 printf(" Reconstructed PID\n"
966 " sgn spec ITS TPC TOF COM\n");
967 for(Int_t idt=0; idt<kNDaughters; idt++){
968 printf(" %c", idt?'-':'+');
969 for(Int_t is(0); is<AliPID::kSPECIES; is++){
970 printf("%s%s%s", is==0?" ":" ", AliPID::ParticleShortName(is), (is==1||is==2)?" ":" ");
971 for(Int_t id(0); id<kNDetectors; id++){
972 printf("%5.1f ", 1.e2*fDetPID[idt][id][is]);
974 printf("%5.1f\n", 1.e2*fComPID[idt][is]);
979 //_________________________________________________
980 void AliTRDv0Info::SetV0tracks(AliESDtrack *p, AliESDtrack *n)
982 //sets the two daughter trex and their indices
983 fTrackP = p; fPindex = p->GetID();
984 fTrackN = n; fNindex = n->GetID();
986 //_________________________________________________
988 AliESDtrack *AliTRDv0Info::GetV0Daughter(Int_t sign)
990 //Gets positive of negative daughter of decay