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