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