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