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