ANALYSIS/AliAODv0KineCuts.cxx

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