Adding a reminder for coders
[u/mrichter/AliRoot.git] / ANALYSIS / AliESDv0KineCuts.cxx
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
40 ClassImp(AliESDv0KineCuts)\r
41 \r
42 //____________________________________________________________________\r
43 AliESDv0KineCuts::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
104 AliESDv0KineCuts::~AliESDv0KineCuts(){\r
105   //\r
106   // Destructor\r
107   //\r
108 \r
109 \r
110 }\r
111 //____________________________________________________________________\r
112 AliESDv0KineCuts::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
135 AliESDv0KineCuts &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
144 void 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
183 Bool_t  AliESDv0KineCuts::ProcessV0(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN) const \r
184 {\r
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
219 Bool_t  AliESDv0KineCuts::ProcessV0(AliESDv0* const v0, Int_t &pdgP, Int_t &pdgN) const \r
220 {\r
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
236 Int_t AliESDv0KineCuts::PreselectV0(AliESDv0* const v0) const \r
237 {\r
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
308 Bool_t  AliESDv0KineCuts::SingleTrackCuts(AliESDv0 * const v0) const \r
309 {\r
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
353 Bool_t AliESDv0KineCuts::CaseGamma(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN) const \r
354 {\r
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
441 Bool_t  AliESDv0KineCuts::CaseK0(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN) const {\r
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
518 Bool_t  AliESDv0KineCuts::CaseLambda(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN, Int_t id) const {\r
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
661 Bool_t  AliESDv0KineCuts::V0CutsCommon(AliESDv0 * const v0) const \r
662 {\r
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
682 void AliESDv0KineCuts::Armenteros(AliESDv0* const v0, Float_t val[2]) const \r
683 {\r
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
719 Bool_t AliESDv0KineCuts::CheckSigns(AliESDv0* const v0) const \r
720 {\r
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
750 Double_t AliESDv0KineCuts::PsiPair(AliESDv0* const v0) const \r
751 {\r
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
822 Bool_t  AliESDv0KineCuts::GetConvPosXY(AliESDtrack * const ptrack, AliESDtrack * const ntrack, Double_t convpos[2]) const\r
823 {\r
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
855 Bool_t  AliESDv0KineCuts::GetHelixCenter(AliESDtrack * const track, Double_t b,Int_t charge, Double_t center[2]) const\r
856 {\r
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
907 AliKFParticle *AliESDv0KineCuts::CreateMotherParticle(const AliVTrack* const pdaughter, const AliVTrack* const ndaughter, Int_t pspec, Int_t nspec) const\r
908 {\r
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
918   m->SetField(fEvent->GetMagneticField());\r
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
941 void  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
953 void  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
966 void 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
977 void 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
1014 void 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
1032 void 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