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] = {static_cast<Float_t>(TMath::Abs(mMass[0] - cL0mass)), static_cast<Float_t>(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] = {static_cast<Float_t>(d[0]->GetP()), static_cast<Float_t>(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 }