]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ANALYSIS/AliESDv0KineCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
No need for TAlienFile.h
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDv0KineCuts.cxx
CommitLineData
7cd657f1 1/**************************************************************************\r
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3 * *\r
4 * Author: The ALICE Off-line Project. *\r
5 * Contributors are mentioned in the code where appropriate. *\r
6 * *\r
7 * Permission to use, copy, modify and distribute this software and its *\r
8 * documentation strictly for non-commercial purposes is hereby granted *\r
9 * without fee, provided that the above copyright notice appears in all *\r
10 * copies and that both the copyright notice and this permission notice *\r
11 * appear in the supporting documentation. The authors make no claims *\r
12 * about the suitability of this software for any purpose. It is *\r
13 * provided "as is" without express or implied warranty. *\r
14 **************************************************************************/\r
15\r
16/*\r
17 * author: M.Kalisky@gsi.de\r
18 * 08/Dec/2010\r
19 *\r
20 * Description: This class allows with purely kinematical cuts\r
21 * to select clean samples of electrons, pions and protons from the\r
22 * V0 online finder ESD V0 candidates for PID and dectector resonse\r
23 * studies.\r
24 */\r
25\r
26#include <TVector3.h>\r
27#include <TDatabasePDG.h>\r
28\r
29#include "AliESDv0.h"\r
30#include "AliESDtrack.h"\r
31#include "AliESDEvent.h"\r
32#include "AliVEvent.h"\r
33#include "AliLog.h"\r
34#include "AliKFParticle.h"\r
35#include "AliVTrack.h"\r
36#include "AliKFVertex.h"\r
37\r
38#include "AliESDv0KineCuts.h"\r
39\r
40ClassImp(AliESDv0KineCuts)\r
41\r
42//____________________________________________________________________\r
43AliESDv0KineCuts::AliESDv0KineCuts() :\r
44 fEvent(0x0)\r
45 , fPrimaryVertex(0x0)\r
46 , fType(0)\r
47 , fMode(0)\r
48 , fTPCNcls(1)\r
49 , fTPCrefit(kTRUE)\r
50 , fTPCchi2perCls(4.0)\r
51 , fTPCclsRatio(0.6)\r
52 , fNoKinks(kTRUE)\r
53 , fGcutChi2NDF(10)\r
54 , fGcutInvMass(0.05)\r
55 , fK0cutChi2NDF(10)\r
56 , fLcutChi2NDF(10)\r
57{\r
58 //\r
59 // Default constructor\r
60 //\r
61\r
62 // default single track cuts\r
63 fTPCNcls = 1; // minimal number of the TPC clusters\r
64 fTPCrefit = kTRUE; // TPC refit\r
65 fTPCchi2perCls = 4.0; // chi2 per TPC cluster\r
66 fTPCclsRatio = 0.6; // minimal foun/findable TPC cluster ratio\r
67 fNoKinks = kTRUE; // kinks - no [kTRUE] or do not care [kFalse]\r
68\r
69\r
70 // default gamma cuts values\r
71 fGcutChi2NDF = 10; // Chi2NF cut value for the AliKFparticle gamma\r
72 fGcutCosPoint[0] = 0; // cos of the pointing angle [min, max]\r
73 fGcutCosPoint[1] = 0.02; // cos of the pointing angle [min, max]\r
74 fGcutDCA[0] = 0.; // DCA between the daughter tracks [min, max]\r
75 fGcutDCA[1] = 0.25; // DCA between the daughter tracks [min, max]\r
76 fGcutVertexR[0] = 3.; // radius of the conversion point [min, max]\r
77 fGcutVertexR[1] = 90.; // radius of the conversion point [min, max]\r
78 fGcutPsiPair[0] = 0.; // value of the psi pair cut [min, max]\r
79 fGcutPsiPair[1] = 0.05; // value of the psi pair cut [min, max]\r
80 fGcutInvMass = 0.05; // upper value on the gamma invariant mass\r
81 // default K0 cuts\r
82 fK0cutChi2NDF = 10; // Chi2NF cut value for the AliKFparticle K0\r
83 fK0cutCosPoint[0] = 0.; // cos of the pointing angle [min, max]\r
84 fK0cutCosPoint[1] = 0.02; // cos of the pointing angle [min, max]\r
85 fK0cutDCA[0] = 0.; // DCA between the daughter tracks [min, max]\r
86 fK0cutDCA[1] = 0.2; // DCA between the daughter tracks [min, max]\r
87 fK0cutVertexR[0] = 2.0; // radius of the decay point [min, max]\r
88 fK0cutVertexR[1] = 30.0; // radius of the decay point [min, max]\r
89 fK0cutInvMass[0] = 0.486; // invariant mass window\r
90 fK0cutInvMass[1] = 0.508; // invariant mass window\r
91 // Lambda & anti-Lambda cut values\r
92 fLcutChi2NDF = 10; // Chi2NF cut value for the AliKFparticle K0\r
93 fLcutCosPoint[0] = 0.; // cos of the pointing angle [min, max]\r
94 fLcutCosPoint[1] = 0.02; // cos of the pointing angle [min, max]\r
95 fLcutDCA[0] = 0.; // DCA between the daughter tracks [min, max]\r
96 fLcutDCA[1] = 0.2; // DCA between the daughter tracks [min, max]\r
97 fLcutVertexR[0] = 2.0; // radius of the decay point [min, max]\r
98 fLcutVertexR[1] = 40.0; // radius of the decay point [min, max]\r
99 fLcutInvMass[0] = 1.11; // invariant mass window\r
100 fLcutInvMass[1] = 1.12; // invariant mass window\r
101 \r
102}\r
103//____________________________________________________________________\r
104AliESDv0KineCuts::~AliESDv0KineCuts(){\r
105 //\r
106 // Destructor\r
107 //\r
108\r
109\r
110}\r
111//____________________________________________________________________\r
112AliESDv0KineCuts::AliESDv0KineCuts(const AliESDv0KineCuts &ref):\r
113 TObject(ref)\r
114 , fEvent(0x0)\r
115 , fPrimaryVertex(0x0)\r
116 , fType(0)\r
117 , fMode(0)\r
118 , fTPCNcls(1)\r
119 , fTPCrefit(kTRUE)\r
120 , fTPCchi2perCls(4.0)\r
121 , fTPCclsRatio(0.6)\r
122 , fNoKinks(kTRUE)\r
123 , fGcutChi2NDF(10)\r
124 , fGcutInvMass(0.05)\r
125 , fK0cutChi2NDF(10)\r
126 , fLcutChi2NDF(10)\r
127{\r
128 //\r
129 // Copy operator\r
130 //\r
131\r
132 ref.Copy(*this);\r
133}\r
134//____________________________________________________________________\r
135AliESDv0KineCuts &AliESDv0KineCuts::operator=(const AliESDv0KineCuts &ref){\r
136 //\r
137 // assignment operator\r
138 //\r
139 if(this != &ref)\r
140 ref.Copy(*this);\r
141 return *this; \r
142}\r
143//____________________________________________________________________\r
144void AliESDv0KineCuts::Copy(TObject &ref) const {\r
145 //\r
146 // Performs the copying of the object\r
147 //\r
148\r
149 TObject::Copy(ref);\r
150\r
151 AliESDv0KineCuts &target = dynamic_cast<AliESDv0KineCuts &>(ref);\r
152\r
153 // default single track cuts\r
154 target.fTPCNcls = fTPCNcls;\r
155 target.fTPCrefit = fTPCrefit;\r
156 target.fTPCchi2perCls = fTPCchi2perCls;\r
157 target.fTPCclsRatio = fTPCclsRatio;\r
158 target.fNoKinks = fNoKinks;\r
159\r
160\r
161 // default gamma cuts values\r
162 target.fGcutChi2NDF = fGcutChi2NDF;\r
163 memcpy(target.fGcutCosPoint, fGcutCosPoint, sizeof(Float_t) * 2);\r
164 memcpy(target.fGcutDCA, fGcutDCA, sizeof(Float_t) * 2); \r
165 memcpy(target.fGcutVertexR, fGcutVertexR, sizeof(Float_t) * 2);\r
166 memcpy(target.fGcutPsiPair, fGcutPsiPair, sizeof(Float_t) * 2);\r
167 target.fGcutInvMass = fGcutInvMass;\r
168 // default K0 cuts\r
169 target.fK0cutChi2NDF = fK0cutChi2NDF;\r
170 memcpy(target.fK0cutCosPoint, fK0cutCosPoint, sizeof(Float_t) * 2);\r
171 memcpy(target.fK0cutDCA, fK0cutDCA, sizeof(Float_t) * 2);\r
172 memcpy(target.fK0cutVertexR, fK0cutVertexR, sizeof(Float_t) * 2);\r
173 memcpy(target.fK0cutInvMass, fK0cutInvMass, sizeof(Float_t) * 2);\r
174 // Lambda & anti-Lambda cut values\r
175 target.fLcutChi2NDF = fLcutChi2NDF;\r
176 memcpy(target.fLcutCosPoint, fLcutCosPoint, sizeof(Float_t) * 2);\r
177 memcpy(target.fLcutDCA, fLcutDCA, sizeof(Float_t) * 2);\r
178 memcpy(target.fLcutVertexR, fLcutVertexR, sizeof(Float_t) * 2);\r
179 memcpy(target.fLcutInvMass, fLcutInvMass, sizeof(Float_t) * 2);\r
180 \r
181}\r
182//____________________________________________________________________\r
5abc3dfe 183Bool_t AliESDv0KineCuts::ProcessV0(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN) const \r
184{\r
7cd657f1 185 //\r
186 // main user function\r
187 //\r
188\r
189 if(!v0) return kFALSE;\r
190 if(!fEvent){\r
191 AliErrorClass("No valid Event pointer available, provide it first");\r
192 return kFALSE;\r
193 }\r
194\r
195 if(!V0CutsCommon(v0)) return kFALSE;\r
196\r
197 const Int_t id = PreselectV0(v0);\r
198\r
199 if(!SingleTrackCuts(v0)) return kFALSE;\r
200\r
201 switch(id){\r
202 case kUndef:\r
203 return kFALSE;\r
204 case kGamma:\r
205 return CaseGamma(v0, pdgV0, pdgP, pdgN);\r
206 case kK0:\r
207 return CaseK0(v0, pdgV0, pdgP, pdgN);\r
208 case kLambda:\r
209 return CaseLambda(v0, pdgV0, pdgP, pdgN, 0);\r
210 case kALambda:\r
211 return CaseLambda(v0, pdgV0, pdgP, pdgN, 1);\r
212 default:\r
213 return kFALSE; \r
214 }\r
215\r
216 return kFALSE;\r
217}\r
218//____________________________________________________________________\r
5abc3dfe 219Bool_t AliESDv0KineCuts::ProcessV0(AliESDv0* const v0, Int_t &pdgP, Int_t &pdgN) const \r
220{\r
7cd657f1 221 //\r
222 // main user function, simplified if the V0 identity is not necessary\r
223 //\r
224\r
225 if(!v0) return kFALSE;\r
226 if(!fEvent){\r
227 AliErrorClass("No valid Event pointer available, provide it first");\r
228 return kFALSE;\r
229 }\r
230\r
231 Int_t idV0 = -1;\r
232 return ProcessV0(v0, idV0, pdgP, pdgN);\r
233\r
234}\r
235//____________________________________________________________________\r
5abc3dfe 236Int_t AliESDv0KineCuts::PreselectV0(AliESDv0* const v0) const \r
237{\r
7cd657f1 238 //\r
239 // Make a preselection (exclusive) of the V0 cadidates based on\r
240 // Armenteros plot\r
241 // the armenteros cut values are currently fixed and user is not able to set them via\r
242 // set funcions. The reason is that these cuts are optimized and furneter changes should \r
243 // not be necessary. To prove otherwise please study in detail before changing the values\r
244 //\r
245 \r
246 Float_t ap[2] = {-1., -1.};\r
247 Armenteros(v0, ap);\r
248 // for clarity\r
249 const Float_t alpha = ap[0];\r
250 const Float_t qt = ap[1];\r
251\r
252 // selection cuts \r
253 // - the reagions for different candidates must not overlap \r
254\r
255 // Gamma cuts\r
256 const Double_t cutAlphaG = 0.35; \r
257 const Double_t cutQTG = 0.05;\r
258 const Double_t cutAlphaG2[2] = {0.6, 0.8};\r
259 const Double_t cutQTG2 = 0.04;\r
260\r
261 // K0 cuts\r
262 const Float_t cutQTK0[2] = {0.1075, 0.215};\r
263 const Float_t cutAPK0[2] = {0.199, 0.8}; // parameters for curved QT cut\r
264 \r
265 // Lambda & A-Lambda cuts\r
266 const Float_t cutQTL = 0.03;\r
267 const Float_t cutAlphaL[2] = {0.35, 0.7};\r
268 const Float_t cutAlphaAL[2] = {-0.7, -0.35};\r
269 const Float_t cutAPL[3] = {0.107, -0.69, 0.5}; // parameters fir curved QT cut\r
270\r
271\r
272 if(kPurity == fMode){\r
273 // Check for Gamma candidates\r
274 if(qt < cutQTG){\r
275 if( (TMath::Abs(alpha) < cutAlphaG) ) return kGamma;\r
276 }\r
277 // additional region - should help high pT gammas\r
278 if(qt < cutQTG2){\r
279 if( (TMath::Abs(alpha) > cutAlphaG2[0]) && (TMath::Abs(alpha) < cutAlphaG2[1]) ) return kGamma;\r
280 }\r
281 }\r
282 if(kEffGamma == fMode){\r
283 if(qt < cutQTG) return kGamma;\r
284 }\r
285\r
286 \r
287 // Check for K0 candidates\r
288 Float_t q = cutAPK0[0] * TMath::Sqrt(TMath::Abs(1 - alpha*alpha/(cutAPK0[1]*cutAPK0[1])));\r
289 if( (qt > cutQTK0[0]) && (qt < cutQTK0[1]) && (qt > q) ){\r
290 return kK0;\r
291 }\r
292\r
293 // Check for Lambda candidates\r
294 q = cutAPL[0] * TMath::Sqrt(TMath::Abs(1 - ( (alpha + cutAPL[1]) * (alpha + cutAPL[1]) ) / (cutAPL[2]*cutAPL[2]) ));\r
295 if( (alpha > cutAlphaL[0]) && (alpha < cutAlphaL[1]) && (qt > cutQTL) && (qt < q) ){\r
296 return kLambda;\r
297 }\r
298\r
299 // Check for A-Lambda candidates\r
300 q = cutAPL[0] * TMath::Sqrt(TMath::Abs(1 - ( (alpha - cutAPL[1]) * (alpha - cutAPL[1]) ) / (cutAPL[2]*cutAPL[2]) ));\r
301 if( (alpha > cutAlphaAL[0]) && (alpha < cutAlphaAL[1]) && (qt > cutQTL) && (qt < q) ){\r
302 return kALambda;\r
303 }\r
304 \r
305 return kUndef;\r
306}\r
307//____________________________________________________________________\r
5abc3dfe 308Bool_t AliESDv0KineCuts::SingleTrackCuts(AliESDv0 * const v0) const \r
309{\r
7cd657f1 310 //\r
311 // apply single track cuts\r
312 // correct sign not relevat here\r
313 //\r
314\r
315 if(!v0) return kFALSE;\r
316 \r
317 Int_t pIndex = 0, nIndex = 0;\r
318 pIndex = v0->GetPindex();\r
319 nIndex = v0->GetNindex();\r
320 AliESDtrack* d[2];\r
321 d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));\r
322 d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));\r
323 \r
324 for(Int_t i=0; i<2; ++i){\r
325 if(!d[i]) return kFALSE;\r
326 \r
327 // status word\r
328 ULong_t status = d[i]->GetStatus();\r
329\r
330 // No. of TPC clusters leave to the users\r
331 if(d[i]->GetTPCNcls() < 1) return kFALSE;\r
332\r
333 // TPC refit\r
334 if(!(status & AliESDtrack::kTPCrefit)) return kFALSE;\r
335 \r
336 // Chi2 per TPC cluster\r
337 Int_t nTPCclusters = d[i]->GetTPCNcls();\r
338 Float_t chi2perTPCcluster = d[i]->GetTPCchi2()/Float_t(nTPCclusters);\r
339 if(chi2perTPCcluster > 4) return kFALSE;\r
340\r
341 // TPC cluster ratio\r
342 Float_t cRatioTPC = d[i]->GetTPCNclsF() > 0. ? static_cast<Float_t>(d[i]->GetTPCNcls())/static_cast<Float_t> (d[i]->GetTPCNclsF()) : 1.;\r
343 if(cRatioTPC < 0.6) return kFALSE;\r
344 \r
345 // kinks\r
346 if(d[i]->GetKinkIndex(0) != 0) return kFALSE;\r
347 \r
348 }\r
349\r
350 return kTRUE;\r
351}\r
352//____________________________________________________________________\r
5abc3dfe 353Bool_t AliESDv0KineCuts::CaseGamma(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN) const \r
354{\r
7cd657f1 355 //\r
356 // process the gamma conversion candidate\r
357 //\r
358\r
359 if(!v0) return kFALSE;\r
360\r
361 AliVTrack* daughter[2];\r
362 Int_t pIndex = 0, nIndex = 0;\r
363\r
364 Bool_t sign = CheckSigns(v0);\r
365 if(sign){\r
366 pIndex = v0->GetPindex();\r
367 nIndex = v0->GetNindex();\r
368 }\r
369 else{\r
370 pIndex = v0->GetNindex();\r
371 nIndex = v0->GetPindex(); \r
372 }\r
373 daughter[0] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(pIndex));\r
374 daughter[1] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(nIndex));\r
375 if(!daughter[0] || !daughter[1]) return kFALSE;\r
376\r
377 AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kElectron), TMath::Abs(kElectron));\r
378 if(!kfMother) return kFALSE;\r
379\r
380 AliESDtrack* d[2];\r
381 d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));\r
382 d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));\r
383\r
384 Float_t iMass = v0->GetEffMass(0, 0);\r
385\r
386 // cos pointing angle\r
387 Double_t cosPoint = v0->GetV0CosineOfPointingAngle();\r
388 cosPoint = TMath::ACos(cosPoint);\r
389\r
390 // DCA between daughters\r
391 Double_t dca = v0->GetDcaV0Daughters();\r
392\r
393 // Production vertex\r
394 Double_t x, y, z; \r
395 v0->GetXYZ(x,y,z);\r
396 Double_t r = TMath::Sqrt(x*x + y*y);\r
397\r
398 Double_t xy[2];\r
399 Double_t r2 = -1.;\r
400 if ( GetConvPosXY(d[0], d[1], xy) ){\r
401 r2 = TMath::Sqrt(xy[0]*xy[0] + xy[1]*xy[1]);\r
402 }\r
403\r
404 // psi pair \r
405 Double_t psiPair = PsiPair(v0);\r
406 \r
407 // V0 chi2/ndf\r
408 Double_t chi2ndf = kfMother->GetChi2()/kfMother->GetNDF();\r
409\r
410 if(kfMother) delete kfMother; \r
411 \r
412 // apply the cuts\r
413\r
414 if(iMass > fGcutInvMass) return kFALSE;\r
415\r
416 if(chi2ndf > fGcutChi2NDF) return kFALSE;\r
417\r
418 if(cosPoint < fGcutCosPoint[0] || cosPoint > fGcutCosPoint[1]) return kFALSE;\r
419\r
420 if(dca < fGcutDCA[0] || dca > fGcutDCA[1]) return kFALSE;\r
421\r
422 if(r < fGcutVertexR[0] || r > fGcutVertexR[1]) return kFALSE;\r
423\r
424 if(psiPair < fGcutPsiPair[0] || psiPair > fGcutPsiPair[1]) return kFALSE;\r
425 \r
426 // all cuts passed\r
427\r
428 pdgV0 = 22;\r
429 if(sign){\r
430 pdgP = -11;\r
431 pdgN = 11;\r
432 }\r
433 else{\r
434 pdgP = 11;\r
435 pdgN = -11;\r
436 }\r
437\r
438 return kTRUE;\r
439}\r
440//____________________________________________________________________\r
5abc3dfe 441Bool_t AliESDv0KineCuts::CaseK0(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN) const {\r
7cd657f1 442 //\r
443 // process the K0 candidate\r
444 //\r
445\r
446 if(!v0) return kFALSE;\r
447 \r
448 AliVTrack* daughter[2];\r
449 Int_t pIndex = 0, nIndex = 0;\r
450 Bool_t sign = CheckSigns(v0);\r
451 if(sign){\r
452 pIndex = v0->GetPindex();\r
453 nIndex = v0->GetNindex();\r
454 }\r
455 else{\r
456 pIndex = v0->GetNindex();\r
457 nIndex = v0->GetPindex(); \r
458 }\r
459 \r
460 daughter[0] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(pIndex));\r
461 daughter[1] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(nIndex));\r
462 if(!daughter[0] || !daughter[1]) return kFALSE;\r
463\r
464 AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kPiPlus));\r
465 if(!kfMother) return kFALSE;\r
466\r
467 AliESDtrack* d[2];\r
468 d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));\r
469 d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));\r
470\r
471 Float_t iMass = v0->GetEffMass(2, 2);\r
472\r
473 // cos pointing angle\r
474 Double_t cosPoint = v0->GetV0CosineOfPointingAngle();\r
475 cosPoint = TMath::ACos(cosPoint);\r
476\r
477 // DCA between daughters\r
478 Double_t dca = v0->GetDcaV0Daughters();\r
479\r
480 // Production vertex\r
481 Double_t x, y, z; \r
482 v0->GetXYZ(x,y,z);\r
483\r
484 Double_t r = TMath::Sqrt(x*x + y*y); \r
485\r
486 // V0 chi2/ndf\r
487 Double_t chi2ndf = kfMother->GetChi2()/kfMother->GetNDF();\r
488 \r
489 if(kfMother) delete kfMother; \r
490\r
491 //\r
492 // apply the cuts\r
493 //\r
494 if(iMass < fK0cutInvMass[0] || iMass > fK0cutInvMass[1]) return kFALSE;\r
495\r
496 if(chi2ndf > fK0cutChi2NDF) return kFALSE;\r
497\r
498 if(cosPoint < fK0cutCosPoint[0] || cosPoint > fK0cutCosPoint[1]) return kFALSE;\r
499\r
500 if(dca < fK0cutDCA[0] || dca > fK0cutDCA[1]) return kFALSE;\r
501\r
502 if(r < fK0cutVertexR[0] || r > fK0cutVertexR[1]) return kFALSE;\r
503\r
504 // all cuts passed\r
505 pdgV0 = 310;\r
506 if(sign){\r
507 pdgP = 211;\r
508 pdgN = -211;\r
509 }\r
510 else{\r
511 pdgP = -211;\r
512 pdgN = 211;\r
513 }\r
514\r
515 return kTRUE;\r
516}\r
517//____________________________________________________________________\r
5abc3dfe 518Bool_t AliESDv0KineCuts::CaseLambda(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN, Int_t id) const {\r
7cd657f1 519 //\r
520 // process teh Lambda and Anti-Lambda candidate\r
521 //\r
522 \r
523 if(!v0) return kFALSE;\r
524\r
525 const Double_t cL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(); // PDG lambda mass\r
526\r
527 AliVTrack* daughter[2];\r
528 Int_t pIndex = 0, nIndex = 0;\r
529 Float_t mMass[2] = {-1., -1.};\r
530 Bool_t sign = CheckSigns(v0);\r
531 if(sign){\r
532 pIndex = v0->GetPindex();\r
533 nIndex = v0->GetNindex();\r
534 mMass[0] = v0->GetEffMass(4, 2);\r
535 mMass[1] = v0->GetEffMass(2, 4);\r
536 }\r
537 else{\r
538 pIndex = v0->GetNindex();\r
539 nIndex = v0->GetPindex(); \r
540 mMass[0] = v0->GetEffMass(2, 4);\r
541 mMass[1] = v0->GetEffMass(4, 2);\r
542 }\r
543 \r
544 daughter[0] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(pIndex));\r
545 daughter[1] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(nIndex));\r
546 if(!daughter[0] || !daughter[1]) return kFALSE;\r
547\r
548 AliKFParticle *kfMother[2] = {0x0, 0x0};\r
549 // Lambda\r
550 kfMother[0] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kProton), TMath::Abs(kPiPlus));\r
551 if(!kfMother[0]) return kFALSE;\r
552 \r
553 // Anti-Lambda\r
554 kfMother[1] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kProton));\r
555 if(!kfMother[1]) return kFALSE;\r
556\r
557 Float_t dMass[2] = {TMath::Abs(mMass[0] - cL0mass), TMath::Abs(mMass[1] - cL0mass)};\r
558 \r
559 AliESDtrack* d[2];\r
560 d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));\r
561 d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));\r
562 if(!d[0] || !d[1]) return kFALSE;\r
563 \r
564 Float_t p[2] = {d[0]->GetP(), d[1]->GetP()}; \r
565\r
566 // check the 3 lambda - antilambda variables\r
567 Int_t check[2] = {-1, -1}; // 0 : lambda, 1 : antilambda\r
568 // 1) momentum of the daughter particles - proton is expected to have higher momentum than pion\r
569 check[0] = (p[0] > p[1]) ? 0 : 1;\r
570 // 2) mass of the mother particle\r
571 check[1] = (dMass[0] < dMass[1]) ? 0 : 1;\r
572 \r
573 // require positive correlation of (1) and (2)\r
574 if(check[0] != check[1]){\r
575 if(kfMother[0]) delete kfMother[0]; \r
576 if(kfMother[1]) delete kfMother[1]; \r
577 return kFALSE;\r
578 }\r
579\r
580 // now that the check[0] == check[1]\r
581 const Int_t type = check[0];\r
582\r
583 // require that the input armenteros preselection agree:\r
584 if(type != id) return kFALSE;\r
585\r
586 Float_t iMass =0.;\r
587 if(sign){\r
588 iMass = (type == 0) ? v0->GetEffMass(4, 2) : v0->GetEffMass(2, 4);\r
589 }\r
590 else{\r
591 iMass = (type == 0) ? v0->GetEffMass(2, 4) : v0->GetEffMass(4, 2);\r
592 }\r
593\r
594 // cos pointing angle\r
595 Double_t cosPoint = v0->GetV0CosineOfPointingAngle();\r
596 cosPoint = TMath::ACos(cosPoint);\r
597\r
598 // DCA between daughters\r
599 Double_t dca = v0->GetDcaV0Daughters();\r
600 \r
601 // Production vertex\r
602 Double_t x, y, z; \r
603 v0->GetXYZ(x,y,z);\r
604 Double_t r = TMath::Sqrt(x*x + y*y);\r
605\r
606 // proton - pion indices\r
607 Int_t ix[2] = {0, 1};\r
608 if(1 == type){\r
609 ix[0] = 1;\r
610 ix[1] = 0;\r
611 }\r
612\r
613 // V0 chi2/ndf\r
614 Double_t chi2ndf = kfMother[type]->GetChi2()/kfMother[type]->GetNDF();\r
615\r
616 if(kfMother[0]) delete kfMother[0]; \r
617 if(kfMother[1]) delete kfMother[1]; \r
618\r
619 //\r
620 // apply the cuts\r
621 //\r
622\r
623 if(iMass < fLcutInvMass[0] || iMass > fLcutInvMass[1]) return kFALSE;\r
624\r
625 if(chi2ndf > fLcutChi2NDF) return kFALSE;\r
626\r
627 if(cosPoint < fLcutCosPoint[0] || cosPoint > fLcutCosPoint[1]) return kFALSE;\r
628\r
629 if(dca < fLcutDCA[0] || dca > fLcutDCA[1]) return kFALSE;\r
630\r
631 if(r < fLcutVertexR[0] || r > fLcutVertexR[1]) return kFALSE;\r
632\r
633 // all cuts passed\r
634\r
635 if(0 == type){\r
636 pdgV0 = 3122;\r
637 if(sign){\r
638 pdgP = 2212;\r
639 pdgN = -211;\r
640 }\r
641 else{\r
642 pdgP = -211;\r
643 pdgN = 2212;\r
644 }\r
645 }\r
646 else{\r
647 pdgV0 = -3122;\r
648 if(sign){\r
649 pdgP = 211;\r
650 pdgN = -2212;\r
651 }\r
652 else{\r
653 pdgP = -2212;\r
654 pdgN = 211;\r
655 }\r
656 }\r
657\r
658 return kTRUE;\r
659}\r
660//____________________________________________________________________\r
5abc3dfe 661Bool_t AliESDv0KineCuts::V0CutsCommon(AliESDv0 * const v0) const \r
662{\r
7cd657f1 663 //\r
664 // V0 cuts common to all V0s\r
665 //\r
666\r
667 AliESDtrack* dN, *dP; \r
668 \r
669 dP = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(v0->GetPindex()));\r
670 dN = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(v0->GetNindex())); \r
671 \r
672 if(!dN || !dP) return kFALSE;\r
673\r
674 Int_t qP = dP->Charge();\r
675 Int_t qN = dN->Charge();\r
676\r
677 if((qP*qN) != -1) return kFALSE;\r
678\r
679 return kTRUE;\r
680}\r
681//____________________________________________________________________\r
5abc3dfe 682void AliESDv0KineCuts::Armenteros(AliESDv0* const v0, Float_t val[2]) const \r
683{\r
7cd657f1 684 //\r
685 // computes the Armenteros variables for given V0\r
686 // fills the histogram\r
687 // returns the values via "val"\r
688 //\r
689 \r
690 Double_t mn[3] = {0,0,0};\r
691 Double_t mp[3] = {0,0,0}; \r
692 Double_t mm[3] = {0,0,0}; \r
693\r
694 if(CheckSigns(v0)){\r
695 v0->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter\r
696 v0->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter\r
697 }\r
698 else{\r
699 v0->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter\r
700 v0->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter\r
701 }\r
702 v0->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother\r
703\r
704 TVector3 vecN(mn[0],mn[1],mn[2]);\r
705 TVector3 vecP(mp[0],mp[1],mp[2]);\r
706 TVector3 vecM(mm[0],mm[1],mm[2]);\r
707 \r
708 Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));\r
709 Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));\r
710 \r
711 Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/\r
712 ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;\r
713 Double_t qt = vecP.Mag()*sin(thetaP);\r
714\r
715 val[0] = alfa;\r
716 val[1] = qt;\r
717}\r
718//____________________________________________________________________\r
5abc3dfe 719Bool_t AliESDv0KineCuts::CheckSigns(AliESDv0* const v0) const \r
720{\r
7cd657f1 721 //\r
722 // check wheter the sign was correctly applied to \r
723 // V0 daughter tracks\r
724 //\r
725 \r
726 Bool_t correct = kFALSE;\r
727\r
728 Int_t pIndex = 0, nIndex = 0;\r
729 pIndex = v0->GetPindex();\r
730 nIndex = v0->GetNindex();\r
731 \r
732 AliESDtrack* d[2];\r
733 d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));\r
734 d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));\r
735\r
736 Int_t sign[2];\r
737 sign[0] = (int)d[0]->GetSign();\r
738 sign[1] = (int)d[1]->GetSign();\r
739 \r
740 if(-1 == sign[0] && 1 == sign[1]){\r
741 correct = kFALSE;\r
742 }\r
743 else{\r
744 correct = kTRUE;\r
745 }\r
746 \r
747 return correct;\r
748}\r
749//________________________________________________________________\r
5abc3dfe 750Double_t AliESDv0KineCuts::PsiPair(AliESDv0* const v0) const \r
751{\r
7cd657f1 752 //\r
753 // Angle between daughter momentum plane and plane \r
754 // \r
755\r
756 if(!fEvent) return -1.;\r
757\r
758 Float_t magField = fEvent->GetMagneticField();\r
759\r
760 Int_t pIndex = -1;\r
761 Int_t nIndex = -1;\r
762 if(CheckSigns(v0)){\r
763 pIndex = v0->GetPindex();\r
764 nIndex = v0->GetNindex();\r
765 }\r
766 else{\r
767 pIndex = v0->GetNindex();\r
768 nIndex = v0->GetPindex(); \r
769 }\r
770 \r
771\r
772 AliESDtrack* daughter[2];\r
773\r
774 daughter[0] = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(pIndex));\r
775 daughter[1] = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(nIndex));\r
776\r
777 Double_t x, y, z;\r
778 v0->GetXYZ(x,y,z);//Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!\r
779 \r
780 Double_t mn[3] = {0,0,0};\r
781 Double_t mp[3] = {0,0,0};\r
782 \r
783\r
784 v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;\r
785 v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter; \r
786\r
787\r
788 Double_t deltat = 1.;\r
789 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\r
790\r
791 Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated\r
792\r
793 Double_t momPosProp[3];\r
794 Double_t momNegProp[3];\r
795 \r
796 AliExternalTrackParam pt(*daughter[0]), nt(*daughter[1]);\r
797 \r
798 Double_t psiPair = 4.;\r
799\r
800 if(nt.PropagateTo(radiussum,magField) == 0)//propagate tracks to the outside\r
801 psiPair = -5.;\r
802 if(pt.PropagateTo(radiussum,magField) == 0)\r
803 psiPair = -5.;\r
804 pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation\r
805 nt.GetPxPyPz(momNegProp);\r
806 \r
807 Double_t pEle =\r
808 TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter\r
809 Double_t pPos =\r
810 TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter\r
811 \r
812 Double_t scalarproduct =\r
813 momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta\r
814 \r
815 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks\r
816\r
817 psiPair = TMath::Abs(TMath::ASin(deltat/chipair)); \r
818\r
819 return psiPair; \r
820}\r
821//___________________________________________________________________\r
5abc3dfe 822Bool_t AliESDv0KineCuts::GetConvPosXY(AliESDtrack * const ptrack, AliESDtrack * const ntrack, Double_t convpos[2]) const\r
823{\r
7cd657f1 824 //\r
825 // recalculate the gamma conversion XY postition\r
826 //\r
827\r
828 const Double_t b = fEvent->GetMagneticField();\r
829\r
830 Double_t helixcenterpos[2];\r
831 GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);\r
832\r
833 Double_t helixcenterneg[2];\r
834 GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);\r
835\r
836 Double_t poshelix[6];\r
837 ptrack->GetHelixParameters(poshelix,b);\r
838 Double_t posradius = TMath::Abs(1./poshelix[4]);\r
839\r
840 Double_t neghelix[6];\r
841 ntrack->GetHelixParameters(neghelix,b);\r
842 Double_t negradius = TMath::Abs(1./neghelix[4]);\r
843\r
844 Double_t xpos = helixcenterpos[0];\r
845 Double_t ypos = helixcenterpos[1];\r
846 Double_t xneg = helixcenterneg[0];\r
847 Double_t yneg = helixcenterneg[1];\r
848\r
849 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);\r
850 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);\r
851\r
852 return 1;\r
853}\r
854//___________________________________________________________________\r
5abc3dfe 855Bool_t AliESDv0KineCuts::GetHelixCenter(AliESDtrack * const track, Double_t b,Int_t charge, Double_t center[2]) const\r
856{\r
7cd657f1 857 //\r
858 // computes the center of the track helix\r
859 //\r
860 \r
861 Double_t pi = TMath::Pi();\r
862 \r
863 Double_t helix[6];\r
864 track->GetHelixParameters(helix,b);\r
865 \r
866 Double_t xpos = helix[5];\r
867 Double_t ypos = helix[0];\r
868 Double_t radius = TMath::Abs(1./helix[4]);\r
869 Double_t phi = helix[2];\r
870\r
871 if(phi < 0){\r
872 phi = phi + 2*pi;\r
873 }\r
874\r
875 phi -= pi/2.;\r
876 Double_t xpoint = radius * TMath::Cos(phi);\r
877 Double_t ypoint = radius * TMath::Sin(phi);\r
878\r
879 if(b<0){\r
880 if(charge > 0){\r
881 xpoint = - xpoint;\r
882 ypoint = - ypoint;\r
883 }\r
884\r
885 if(charge < 0){\r
886 xpoint = xpoint;\r
887 ypoint = ypoint;\r
888 }\r
889 }\r
890 if(b>0){\r
891 if(charge > 0){\r
892 xpoint = xpoint;\r
893 ypoint = ypoint;\r
894 }\r
895\r
896 if(charge < 0){\r
897 xpoint = - xpoint;\r
898 ypoint = - ypoint;\r
899 }\r
900 }\r
901 center[0] = xpos + xpoint;\r
902 center[1] = ypos + ypoint;\r
903\r
904 return 1;\r
905}\r
906//___________________________________________________________________\r
5abc3dfe 907AliKFParticle *AliESDv0KineCuts::CreateMotherParticle(const AliVTrack* const pdaughter, const AliVTrack* const ndaughter, Int_t pspec, Int_t nspec) const\r
908{\r
7cd657f1 909 //\r
910 // Creates a mother particle\r
911 //\r
912 AliKFParticle pkfdaughter(*pdaughter, pspec);\r
913 AliKFParticle nkfdaughter(*ndaughter, nspec);\r
914 \r
915 \r
916 // Create the mother particle \r
917 AliKFParticle *m = new AliKFParticle(pkfdaughter, nkfdaughter);\r
5abc3dfe 918 m->SetField(fEvent->GetMagneticField());\r
7cd657f1 919 if(TMath::Abs(kElectron) == pspec && TMath::Abs(kElectron) == nspec) m->SetMassConstraint(0, 0.001);\r
920 else if(TMath::Abs(kPiPlus) == pspec && TMath::Abs(kPiPlus) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass(), 0.);\r
921 else if(TMath::Abs(kProton) == pspec && TMath::Abs(kPiPlus) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(), 0.);\r
922 else if(TMath::Abs(kPiPlus) == pspec && TMath::Abs(kProton) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(), 0.);\r
923 else{\r
924 AliErrorClass("Wrong daughter ID - mass constraint can not be set");\r
925 }\r
926\r
927 AliKFVertex improvedVertex = *fPrimaryVertex;\r
928 improvedVertex += *m;\r
929 m->SetProductionVertex(improvedVertex);\r
930 \r
931 // update 15/06/2010\r
932 // mother particle will not be added to primary vertex but only to its copy \r
933 // as this confilcts with calling\r
934 // m->SetPrimaryVertex() function and\r
935 // subsequently removing the mother particle afterwards\r
936 // Source: Sergey Gorbunov\r
937\r
938 return m;\r
939}\r
940//____________________________________________________________________\r
941void AliESDv0KineCuts::SetEvent(AliESDEvent* const event){\r
942 //\r
943 // direct setter of ESD event\r
944 //\r
945 fEvent = event;\r
946 if(!fEvent){\r
947 AliErrorClass("Invalid input event pointer");\r
948 return;\r
949 }\r
950\r
951}\r
952//____________________________________________________________________\r
953void AliESDv0KineCuts::SetEvent(AliVEvent* const event){\r
954 //\r
955 // direct setter of ESD event\r
956 //\r
957 if(event)\r
958 fEvent = static_cast<AliESDEvent*>(event);\r
959 if(!fEvent){\r
960 AliErrorClass("Invalid input event pointer");\r
961 return;\r
962 }\r
963\r
964}\r
965//________________________________________________________________\r
966void AliESDv0KineCuts::SetPrimaryVertex(AliKFVertex* const v){\r
967 //\r
968 // set the primary vertex of the event\r
969 //\r
970 fPrimaryVertex = v;\r
971 if(!fPrimaryVertex){\r
972 AliErrorClass("Failed to initialize the primary vertex");\r
973 return;\r
974 }\r
975}\r
976//___________________________________________________________________\r
977void AliESDv0KineCuts::SetMode(Int_t mode, Int_t type){\r
978 //\r
979 // this function allows the user to select (prior running the 'ProcessV0' function)\r
980 // to select different approaches to V0 selection - the 'mode'\r
981 // - and -\r
982 // different systems (pp, PbPb) - 'type' \r
983 //\r
984 // To see the cut values for different modes please refer to the\r
985 // function SetCuts()\r
986 //\r
987 // Important notice: based on the parameters particular sets of cuts will\r
988 // be activated for teh V0 selection. If some additional changes to single\r
989 // cuts are needed please us the SetXXXcut function (see the header file)\r
990 // \r
991\r
992 switch(mode){\r
993 case kPurity:\r
994 fMode = kPurity; // used to obtain highest purity possible - the efficiency may be low\r
995 case kEffGamma:\r
996 fMode = kEffGamma; // used to obtain highes efficiency possible - the purity may be worse\r
997 default:\r
998 AliError("V0 selection mode not recognozed, setting 'kPurity'");\r
999 fMode = kPurity;\r
1000 }\r
1001\r
1002 switch(type){\r
1003 case kPP:\r
1004 fType = kPP; // cuts optimized for low multiplicity \r
1005 case kPbPb:\r
1006 fType = kPbPb; // cuts optimized for high multiplicity\r
1007 }\r
1008 \r
1009 // setup the cut values for selected mode & type\r
1010 SetCuts();\r
1011\r
1012}\r
1013//___________________________________________________________________\r
1014void AliESDv0KineCuts::SetMode(Int_t mode, const char* type){\r
1015 //\r
1016 // overloaded function - please see above\r
1017 // \r
1018 \r
1019 Int_t t = -1;\r
1020\r
1021 if(!strcmp("pp", type)) t = kPP;\r
1022 else if(!(strcmp("PbPb", type))) t = kPbPb;\r
1023 else{\r
1024 AliError("data type not recognized, setting 'pp'");\r
1025 t = kPP; \r
1026 }\r
1027\r
1028 SetMode(mode, t);\r
1029\r
1030}\r
1031//___________________________________________________________________\r
1032void AliESDv0KineCuts::SetCuts(){\r
1033 //\r
1034 // this funciton sets the default cut values based on the selected\r
1035 // fMode and fType.\r
1036 // please note that only the cuts that have different values than the default\r
1037 // cuts are updated here\r
1038 //\r
1039 \r
1040 // last update: 14/02/2011\r
1041 // as a very preliminary - the only change to default cuts is to apply\r
1042 // less restricting gamma conversion selection in PreselectV0() function\r
1043 \r
1044\r
1045 \r
1046}\r