PWG3/hfe/AliHFEV0pid.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 // Utility class for V0 PID
  17 // Identifies Electrons, Pions and Protons using gamma conversions and
  18 // the decays of K0s and Lambdas
  19 // Containers with samples of Electrons, Pions and Protons can be accessed
  20 // via GetListOfElectrons() etc.
  21 //
  22 // Authors:
  23 //    Matus Kalisky <matus.kalisky@cern.ch>
  24 //    Markus Heide <mheide@uni-muenster.de>
  25 //    Markus Fasel <M.Fasel@gsi.de>
  26 //
  27 #include <TDatabasePDG.h>
  28 #include <TObjArray.h>
  29 #include <TPDGCode.h>
  30 #include <TString.h>
  31
  32 #include <TDatabasePDG.h>
  33
  34 #include "AliAODEvent.h"
  35 #include "AliAODTrack.h"
  36 #include "AliAODv0.h"
  37 #include "AliAODVertex.h"
  38 #include "AliESDEvent.h"
  39 #include "AliESDtrack.h"
  40 #include "AliESDv0.h"
  41 #include "AliESDVertex.h"
  42 #include "AliKFParticle.h"
  43 #include "AliKFVertex.h"
  44 #include "AliVEvent.h"
  45 #include "AliVTrack.h"
  46
  47 #include "AliHFEV0pid.h"
  48 ClassImp(AliHFEV0pid)
  49
  50 //____________________________________________________________
  51 AliHFEV0pid::AliHFEV0pid():
  52   TObject()
  53   , fInputEvent(NULL)
  54   , fPrimaryVertex(NULL)
  55   , fElectrons(NULL)
  56   , fPionsK0(NULL)
  57   , fPionsL(NULL)
  58   , fKaons(NULL)
  59   , fProtons(NULL)
  60   , fIndices(NULL)
  61   , fQA(NULL)
  62 {
  63   //
  64   // Default constructor
  65   //
  66   fElectrons = new TObjArray();
  67   fPionsK0 = new TObjArray();
  68   fPionsL = new TObjArray();
  69   fKaons = new TObjArray();
  70   fProtons = new TObjArray();
  71   fIndices = new AliHFEV0pidTrackIndex;
  72 }
  73
  74 //____________________________________________________________
  75 AliHFEV0pid::~AliHFEV0pid(){
  76   //
  77   // Destructor
  78   // Remove Containers
  79   //
  80   if(fInputEvent) delete fInputEvent;
  81   //if(fPrimaryVertex) delete fPrimaryVertex;
  82   if(fElectrons) delete fElectrons;
  83   if(fPionsK0) delete fPionsK0;
  84   if(fPionsL) delete fPionsL;
  85   if(fKaons) delete fKaons;
  86   if(fProtons) delete fProtons;
  87   if(fIndices) delete fIndices;
  88   if(fQA) delete fQA;
  89 }
  90
  91 //____________________________________________________________
  92 void AliHFEV0pid::InitQA(){
  93   //
  94   // Initialize QA histograms
  95   //
  96   if(!fQA){
  97     fQA = new AliHFEcollection("v0pidQA", "QA histograms for V0 PID");
  98
  99     // QA histograms for cut statistics
 100     fQA->CreateTH1F("h_cutEfficiencyGamma", "Cut Efficiency for Gammas", 10, 0, 10);
 101     fQA->CreateTH1F("h_cutEfficiencyK0s", "Cut Efficiency for K0s", 10, 0, 10);
 102     fQA->CreateTH1F("h_cutEfficiencyPhi", "Cut Efficiency for Phi", 10, 0, 10);
 103     fQA->CreateTH1F("h_cutEfficiencyLambda", "Cut Efficiency for Lambdas", 10, 0, 10);
 104
 105     // QA histograms for invariant mass
 106     fQA->CreateTH1F("h_InvMassGamma", "Gamma invariant mass; inv mass [GeV/c^{2}]; counts", 100, 0, 0.25);
 107     fQA->CreateTH1F("h_InvMassK0s", "K0s invariant mass; inv mass [GeV/c^{2}]; counts", 100, 0.4, 0.65);
 108     fQA->CreateTH1F("h_InvMassPhi", "Phi invariant mass; inv mass [GeV/c^{2}]; counts", 100, 0.4, 0.65);
 109     fQA->CreateTH1F("h_InvMassLambda", "Lambda invariant mass; inv mass [GeV/c^{2}]; counts", 100, 1.05, 1.15);
 110
 111     // QA histograms for p distribution (of the daughters)
 112     fQA->CreateTH1F("h_P_electron", "P distribution of the gamma electrons; p (GeV/c); counts", 100, 0.1, 10);
 113     fQA->CreateTH1F("h_P_K0pion", "P distribution of the K0 pions; p (GeV/c); counts", 100, 0.1, 10);
 114     fQA->CreateTH1F("h_P_Lpion", "P distribution of the Lambda pions; p (GeV/c); counts", 100, 0.1, 10);
 115     fQA->CreateTH1F("h_P_Lproton", "P distribution of the Lambda protons; p (GeV/c); counts", 100, 0.1, 10);
 116
 117     // QA invariant mass as a functin of pt
 118     fQA->CreateTH1Fvector1(20, "h_InvMassGamma_pt", "Gamma invarinat mass in pt bins; inv mass [GeV/c^{2}]; counts", 250, 0, 2);
 119     fQA->CreateTH1Fvector1(20, "h_InvMassK0_pt", "K0 invarinat mass in pt bins; inv mass [GeV/c^{2}]; counts", 250, 0, 2);
 120     fQA->CreateTH1Fvector1(20, "h_InvMassPhi_pt", "Phi invarinat mass in pt bins; inv mass [GeV/c^{2}]; counts", 250, 0, 2);
 121     fQA->CreateTH1Fvector1(20, "h_InvMassLambda_pt", "Lambda invarinat mass in pt bins; inv mass [GeV/c^{2}]; counts", 250, 0, 2);
 122
 123     // QA pt of the V0
 124     fQA->CreateTH1F("h_Pt_Gamma", "Pt of the gamma conversion; p_{T} (GeV/c); counts", 100, 0, 10);
 125     fQA->CreateTH1F("h_Pt_K0", "Pt of the K0; p_{T} (GeV/c); counts", 100, 0, 10);
 126     fQA->CreateTH1F("h_Pt_Phi", "Pt of the Phi; p_{T} (GeV/c); counts", 100, 0, 10);
 127     fQA->CreateTH1F("h_Pt_Lambda", "Pt of the Lambda; p_{T} (GeV/c); counts", 100, 0, 10);
 128     //fQA->CreateTH1F("h_Pt_electrons", "Pt of the conversion electrons; p_{T} (GeV/c); counts");
 129     //fQA->CreateTH1F("h_Pt_pionsK0", "Pt of the K0 pions; p_{T} (GeV/c); counts");
 130     //fQA->CreateTH1F("h_Pt_pionsL", "Pt of the Lambda pions; p_{T} (GeV/c); counts");
 131     //fQA->CreateTH1F("h_Pt_protons", "Pt of the Lambda protons; p_{T} (GeV/c); counts");
 132
 133
 134     // QA histogram for both Lambda candidate combinations -
 135     fQA->CreateTH2F("h_L0_dca_v_dMass", "L0 dca verus dMass; dMass [GeV/c^{2}]; dDCA [cm]; ", 100, -1., 1., 100, 0., 5.);
 136
 137     // Chi2 histograms
 138     fQA->CreateTH1F("h_chi2_gamma", "Chi2 for gammas", 10000, 0, 1000);
 139     fQA->CreateTH1F("h_chi2_K0s", "Chi2 for K0s", 10000, 0, 500);
 140     fQA->CreateTH1F("h_chi2_Phi", "Chi2 for K0s", 10000, 0, 500);
 141     fQA->CreateTH1F("h_chi2_Lambda", "Chi2 for Lambdas", 10000, 0, 1000);
 142   }
 143 }
 144
 145 //____________________________________________________________
 146 void AliHFEV0pid::Process(AliVEvent * const inputEvent){
 147
 148   //
 149   // Find protons, pions and electrons using V0 decays and
 150   // store the pointers in the TObjArray
 151   //
 152   Int_t nGamma = 0, nK0s = 0, nLambda = 0, nPhi = 0;
 153   fInputEvent = inputEvent;
 154
 155   fIndices->Init(fInputEvent->GetNumberOfV0s() * 2);
 156   fPrimaryVertex = new AliKFVertex(*(fInputEvent->GetPrimaryVertex()));
 157   Int_t v0status = 0;
 158   for(Int_t iv0 = 0; iv0 < fInputEvent->GetNumberOfV0s(); iv0++){
 159     if(!TString(fInputEvent->IsA()->GetName()).CompareTo("AliESDEvent")){
 160       // case ESD
 161       SetESDanalysis();
 162       AliESDv0 *esdV0 = (dynamic_cast<AliESDEvent *>(fInputEvent))->GetV0(iv0);
 163       if(!esdV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
 164       v0status = ProcessV0(esdV0);
 165     } else {
 166       // case AOD
 167       SetAODanalysis();
 168       AliAODv0 *aodV0 = (dynamic_cast<AliAODEvent *>(fInputEvent))->GetV0(iv0);
 169       if(aodV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
 170       v0status = ProcessV0(aodV0);
 171     }
 172     switch(v0status){
 173     case kRecoGamma: nGamma++; break;
 174     case kRecoK0s: nK0s++; break;
 175     case kRecoPhi: nPhi++; break;
 176     case kRecoLambda: nLambda++; break;
 177     };
 178   }
 179
 180   AliDebug(1, Form("Number of gammas  : %d", nGamma));
 181   AliDebug(1, Form("Number of K0s     : %d", nK0s));
 182   AliDebug(1, Form("Number of Phis    : %d", nPhi));
 183   AliDebug(1, Form("Number of Lambdas : %d", nLambda));
 184
 185   AliDebug(1, "Number of stored tracks:");
 186   AliDebug(1, Form("Number of electrons      : %d", fElectrons->GetEntries()));
 187   AliDebug(1, Form("Number of K0 pions       : %d", fPionsK0->GetEntries()));
 188   AliDebug(1, Form("Number of Lambda pions   : %d", fPionsL->GetEntries()));
 189   AliDebug(1, Form("Number of Phi kaons      : %d", fKaons->GetEntries()));
 190   AliDebug(1, Form("Number of protons        : %d", fProtons->GetEntries()));
 191
 192   delete  fPrimaryVertex;
 193 }
 194
 195 //____________________________________________________________
 196 Int_t AliHFEV0pid::ProcessV0(TObject *v0){
 197   //
 198   // Process single V0
 199   // Apply general cut and special cuts for gamma, K0s, Lambda
 200   //
 201   AliVTrack* daughter[2];
 202   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack((dynamic_cast<AliESDv0 *>(v0))->GetPindex()));
 203   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack((dynamic_cast<AliESDv0 *>(v0))->GetPindex()));
 204   if(!daughter[0] || !daughter[1]) return kUndef;
 205
 206   if(IsESDanalysis()){
 207     for(Int_t i=0; i<2; ++i){
 208       if(!CutESDtrack(dynamic_cast<AliESDtrack*>(daughter[i]))) return kUndef;
 209     }
 210   }
 211
 212   if(IsGammaConv(v0)) return kRecoGamma;
 213   else if(IsK0s(v0)) return kRecoK0s;
 214   else if(IsLambda(v0)) return kRecoLambda;
 215   else return kUndef;
 216
 217
 218 }
 219 //____________________________________________________________
 220 void AliHFEV0pid::Flush(){
 221   //
 222   // Clear the Lists
 223   //
 224   AliDebug(1, "Flushing containers");
 225   fProtons->Clear();
 226   fPionsK0->Clear();
 227   fPionsL->Clear();
 228   fElectrons->Clear();
 229   fIndices->Flush();
 230 }
 231
 232 //____________________________________________________________
 233 Bool_t AliHFEV0pid::IsGammaConv(TObject *v0){
 234   //
 235   // Identify Gamma
 236   //
 237   AliVTrack* daughter[2];
 238   Int_t pIndex = 0, nIndex = 0;
 239   Double_t invMass = 0.;
 240   if(IsESDanalysis()){
 241     // ESD - cut V0
 242     AliESDv0 *esdV0 = dynamic_cast<AliESDv0 *>(v0);
 243     if(!CutV0(esdV0, kRecoGamma)) return kFALSE;
 244     if(LooseRejectK0(esdV0) || LooseRejectLambda(esdV0)) return kFALSE;
 245     // DEBUG
 246     //invMass = esdV0->GetEffMass(AliPID::kElectron, AliPID::kElectron);
 247     invMass = GetEffMass(esdV0, AliPID::kElectron, AliPID::kElectron);
 248     //..
 249
 250     pIndex = esdV0->GetPindex();
 251     nIndex = esdV0->GetNindex();
 252   } else {
 253     // AOD Analysis - not possible to cut
 254     AliAODv0 *aodV0 = dynamic_cast<AliAODv0 *>(v0);
 255     pIndex = aodV0->GetPosID();
 256     nIndex = aodV0->GetNegID();
 257     invMass = aodV0->InvMass2Prongs(0, 1, kElectron, kElectron);
 258   }
 259   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
 260   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
 261   if(!daughter[0] || !daughter[1]) return kFALSE;
 262
 263   // Get Invariant mass and chi2/ndf
 264   AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kElectron), TMath::Abs(kElectron));
 265   kfMother->SetProductionVertex(*fPrimaryVertex);
 266   kfMother->SetMassConstraint(0, 1.);
 267   Double_t ndf = kfMother->GetNDF();
 268   Double_t chi2 = kfMother->GetChi2();
 269   delete kfMother;
 270
 271   if(fQA) fQA->Fill("h_chi2_gamma", chi2/ndf);
 272   if(chi2/ndf > 7) return kFALSE;
 273   Double_t mPt = kfMother->GetPt();
 274   fQA->Fill("h_Pt_Gamma", mPt);
 275   Int_t ptBin = (int)(mPt*10.0);
 276
 277   if(fQA) fQA->Fill("h_InvMassGamma", invMass);
 278   if(invMass > 0.05) return kFALSE;
 279   fQA->Fill("h_InvMassGamma_pt", ptBin+1, invMass);
 280
 281   AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
 282   // Identified gamma - store tracks in the electron containers
 283   if(!fIndices->Find(daughter[0]->GetID())){
 284     AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
 285     fElectrons->Add(daughter[0]);
 286     fIndices->Add(daughter[0]->GetID(), AliHFEV0pid::kRecoElectron);
 287   }
 288   if(!fIndices->Find(daughter[1]->GetID())){
 289     AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[1]->GetID(), daughter[1]->GetID()));
 290     fElectrons->Add(daughter[1]);
 291     fIndices->Add(daughter[1]->GetID(), AliHFEV0pid::kRecoElectron);
 292   }
 293   return kTRUE;
 294 }
 295
 296
 297
 298 //____________________________________________________________
 299 Bool_t AliHFEV0pid::IsK0s(TObject *v0){
 300   //
 301   // Identify K0s
 302   //
 303   const Double_t kK0smass=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();  // PDG K0s mass
 304   AliVTrack* daughter[2];
 305   Int_t pIndex = 0, nIndex = 0;
 306   Double_t invMass = 0.;
 307   if(IsESDanalysis()){
 308     // ESD - cut V0
 309     AliESDv0 *esdV0 = dynamic_cast<AliESDv0 *>(v0);
 310     if(!CutV0(esdV0, kRecoK0s)) return kFALSE;
 311     invMass = esdV0->GetEffMass(AliPID::kPion, AliPID::kPion);
 312     pIndex = esdV0->GetPindex();
 313     nIndex = esdV0->GetNindex();
 314   } else {
 315     // AOD Analysis - not possible to cut
 316     AliAODv0 *aodV0 = dynamic_cast<AliAODv0 *>(v0);
 317     pIndex = aodV0->GetPosID();
 318     nIndex = aodV0->GetNegID();
 319     invMass = aodV0->MassK0Short();
 320   }
 321   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
 322   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
 323   if(!daughter[0] || !daughter[1]) return kFALSE;
 324
 325   // Get Invariant mass and chi2/ndf
 326   AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kPiPlus));
 327   kfMother->SetProductionVertex(*fPrimaryVertex);
 328   kfMother->SetMassConstraint(kK0smass, 0.);
 329   Double_t ndf = kfMother->GetNDF();
 330   Double_t chi2 = kfMother->GetChi2();
 331   delete kfMother;
 332
 333   if(fQA) fQA->Fill("h_chi2_K0s", chi2/ndf);
 334   if(chi2/ndf > 5) return kFALSE;
 335   Double_t mPt = kfMother->GetPt();
 336   fQA->Fill("h_Pt_K0", mPt);
 337   Int_t ptBin = (int)(mPt*10.0);
 338   fQA->Fill("h_InvMassK0_pt", ptBin+1, invMass);
 339
 340   if(fQA) fQA->Fill("h_InvMassK0s", invMass);
 341
 342   if(invMass < 0.485 || invMass > 0.51) return kFALSE;
 343   AliDebug(1, Form("K0 identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
 344
 345   // Identified gamma - store tracks in the electron containers
 346   if(!fIndices->Find(daughter[0]->GetID())){
 347     AliDebug(1, Form("Adding K0 Pion track with ID %d", daughter[0]->GetID()));
 348     fPionsK0->Add(daughter[0]);
 349     fIndices->Add(daughter[0]->GetID(), AliHFEV0pid::kRecoPionK0);
 350   }
 351   if(!fIndices->Find(daughter[1]->GetID())){
 352     AliDebug(1, Form("Adding K0 Pion track with ID %d", daughter[1]->GetID()));
 353     fPionsK0->Add(daughter[1]);
 354     fIndices->Add(daughter[1]->GetID(), AliHFEV0pid::kRecoPionK0);
 355   }
 356   return kTRUE;
 357 }
 358
 359 //____________________________________________________________
 360 Bool_t AliHFEV0pid::IsPhi(TObject *v0){
 361   //
 362   // Identify Phi - very preliminary - requires diffrent approach (V0 fnder is not effective)
 363   //
 364
 365   //const Double_t kPhiMass=TDatabasePDG::Instance()->GetParticle(333)->Mass();  // PDG phi mass
 366   //AliVTrack* daughter[2];
 367   //AliKFParticle *mother = NULL;
 368   //Double_t invMass = 0.;
 369
 370   Int_t pIndex = 0, nIndex = 0;
 371   if(IsESDanalysis()){
 372     // ESD - cut V0
 373     AliESDv0 *esdV0 = dynamic_cast<AliESDv0 *>(v0);
 374     if(!CutV0(esdV0, kRecoPhi)) return kFALSE;
 375     if(LooseRejectGamma(esdV0) || LooseRejectK0(esdV0)) return kFALSE;
 376     pIndex = esdV0->GetPindex();
 377     nIndex = esdV0->GetNindex();
 378   } else {
 379     // PRELIMINARY - !!!
 380     // AOD Analysis - not possible to cut
 381   }
 382
 383   return kTRUE;
 384 }
 385
 386 //____________________________________________________________
 387 Bool_t AliHFEV0pid::IsLambda(TObject *v0){
 388   //
 389   // Identify Lambda
 390   //
 391   const Double_t kL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();  // PDG lambda mass
 392   AliVTrack* daughter[2];
 393   Int_t pIndex = 0, nIndex = 0;
 394   AliKFParticle *mother[2] = {NULL, NULL};
 395   Float_t dMass[2]; // lambda mass difference for the two hypoteses
 396   //Int_t   lambda = 0; // [1] for lambda and [-1] for anti-lambda
 397   Double_t chi2 = 0.;
 398
 399   Double_t invMass = 0.;
 400   Double_t invMassEOD = 0;
 401   if(IsESDanalysis()){
 402     // ESD - cut V0
 403     AliESDv0 *esdV0 = dynamic_cast<AliESDv0 *>(v0);
 404     if(!CutV0(esdV0, kRecoLambda)) return kFALSE;
 405     if(LooseRejectK0(esdV0) || LooseRejectGamma(esdV0)) return kFALSE;
 406     pIndex = esdV0->GetPindex();
 407     nIndex = esdV0->GetNindex();
 408
 409     //invMass = esdV0->GetEffMass(AliPID::kPion, AliPID::kPion);
 410
 411   } else {
 412     // PRELIMINARY - !!!
 413     // AOD Analysis - not possible to cut
 414
 415     // again - two cases as above
 416     AliAODv0 *aodV0 = dynamic_cast<AliAODv0 *>(v0);
 417     pIndex = aodV0->GetPosID();
 418     nIndex = aodV0->GetNegID();
 419     invMassEOD = aodV0->MassLambda();
 420   }
 421
 422   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
 423   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
 424   if(!daughter[0] || !daughter[1]) return kFALSE;
 425
 426   //
 427   // now - go over two case - Lambda and AntiLambda
 428   // choose the on ewhere the resulting lambda mass is close to the
 429   // expected value and at the same time points closer to the primary vertex in XY plane
 430   //
 431   // A)      lambda -> proton + negative-pion
 432   // B) anti-lambda -> anti-proton + positive-pion
 433   //
 434
 435   //
 436   // A) :: proton
 437   //
 438   mother[0] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kProton), TMath::Abs(kPiPlus));
 439   AliHFELambdaInf lambda0(mother[0], fPrimaryVertex);
 440
 441   // Undo changes
 442   *fPrimaryVertex -= *mother[0];
 443   AddTrackToKFVertex(daughter[0], TMath::Abs(kProton));
 444   AddTrackToKFVertex(daughter[1], TMath::Abs(kPiPlus));
 445
 446   //
 447   // B) :: anti-proton
 448   //
 449   mother[1] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kProton));
 450   AliHFELambdaInf lambda1(mother[1], fPrimaryVertex);
 451
 452   // qa histograms
 453   dMass[0] = lambda0.GetInvariantMass() - kL0mass;
 454   dMass[1] = lambda1.GetInvariantMass() - kL0mass;
 455   fQA->Fill("h_L0_dca_v_dMass", dMass[0], lambda0.GetDistanceFromPrimaryVertex());
 456   fQA->Fill("h_L0_dca_v_dMass", dMass[1], lambda1.GetDistanceFromPrimaryVertex());
 457
 458   //
 459   // decide on the true Lambda candidate bsaed on
 460   // - mass difference
 461   // - vertex dca  (testing needed first)
 462   //
 463   AliHFELambdaInf *lambdaInf[2] = {&lambda0, &lambda1};
 464   Bool_t isLambda = kFALSE;
 465   Int_t index = -1;
 466   index = (dMass[0] <dMass[1]) ? 0 : 1;
 467   invMass = lambdaInf[index]->GetInvariantMass();
 468   chi2 = lambdaInf[index]->GetChi2NDF();
 469   //if(!IsESDanalysis()){
 470   //  invMass = invMassEOD;
 471   //}
 472   //else{
 473   //  invMass = mother[index]->GetMass();
 474   //}
 475   if(fQA) fQA->Fill("h_chi2_Lambda", chi2);
 476   if(chi2 < 3){
 477     if(fQA) fQA->Fill("h_InvMassLambda", invMass);
 478     Double_t mPt = mother[index]->GetPt();
 479     fQA->Fill("h_Pt_Lambda", mPt);
 480     Int_t ptBin = (int)(mPt*10.0);
 481     fQA->Fill("h_InvMassLambda_pt", ptBin+1, invMass);
 482
 483     // cut on the invariant mass for the proton and pion selection
 484     if(invMass > 1.11 || invMass < 1.12){
 485       // Identified lambdas - store the protons and pions and update primary vertex
 486       *fPrimaryVertex += *mother[index];
 487
 488       // lambda
 489       if(0 == index){
 490         if(!fIndices->Find(daughter[0]->GetID())){
 491                 fProtons->Add(daughter[0]);
 492                 fIndices->Add(daughter[0]->GetID(), AliHFEV0pid::kRecoProton);
 493         }
 494         if(!fIndices->Find(daughter[1]->GetID())){
 495                 fPionsL->Add(daughter[1]);
 496                 fIndices->Add(daughter[1]->GetID(), AliHFEV0pid::kRecoPionL);
 497         }
 498       }
 499       // antilambda
 500       if(1 == index){
 501         if(!fIndices->Find(daughter [1]->GetID())){
 502                 fProtons->Add(daughter[1]);
 503                 fIndices->Add(daughter[1]->GetID(), AliHFEV0pid::kRecoProton);
 504         }
 505         if(!fIndices->Find(daughter [0]->GetID())){
 506                 fPionsL->Add(daughter[0]);
 507                 fIndices->Add(daughter [0]->GetID(), AliHFEV0pid::kRecoPionL);
 508         }
 509       }
 510       isLambda = kTRUE;
 511     }
 512   }
 513
 514   if(!isLambda){
 515     // Add the daughters again to the primary vertex
 516     AddTrackToKFVertex(daughter[0], TMath::Abs(kPiPlus));
 517     AddTrackToKFVertex(daughter[1], TMath::Abs(kProton));
 518   }
 519   // remove the objects
 520   for(Int_t i=0; i<2; ++i){
 521     if (mother[i]) delete mother[i];
 522   }
 523
 524   return isLambda;
 525 }
 526
 527 //____________________________________________________________
 528 AliKFParticle *AliHFEV0pid::CreateMotherParticle(AliVTrack *pdaughter, AliVTrack *ndaughter, Int_t pspec, Int_t nspec){
 529   //
 530   // Creates a mother particle
 531   //
 532   AliKFParticle pkfdaughter(*pdaughter, pspec);
 533   AliKFParticle nkfdaughter(*ndaughter, nspec);
 534
 535   // check if the daughter particles are coming from the primary vertex
 536   if(IsESDanalysis()){
 537     // ESD Analyis
 538     const AliESDVertex *esdvertex = dynamic_cast<const AliESDVertex *>(fInputEvent->GetPrimaryVertex());
 539     UShort_t *contrib = esdvertex->GetIndices();
 540
 541     Int_t nfound = 0;
 542     for(Int_t id = 0; id < esdvertex->GetNIndices(); id++){
 543       if(contrib[id] == pdaughter->GetID()){
 544         *fPrimaryVertex -= pkfdaughter;
 545         nfound++;
 546       }
 547       if(contrib[id] == ndaughter->GetID()){
 548         *fPrimaryVertex -= nkfdaughter;
 549         nfound++;
 550       }
 551       if(nfound == 2) break;
 552     }
 553   } else {
 554     // AOD Analysis: AOD Vertex
 555     const AliAODVertex *aodvertex = dynamic_cast<const AliAODVertex *>(fInputEvent->GetPrimaryVertex());
 556     if(aodvertex->HasDaughter(pdaughter))
 557       *fPrimaryVertex -= pkfdaughter;
 558     if(aodvertex->HasDaughter(ndaughter))
 559       *fPrimaryVertex -= nkfdaughter;
 560   }
 561
 562   // Create the mother particle and add them to the primary vertex
 563   AliKFParticle *mother = new AliKFParticle(pkfdaughter, nkfdaughter);
 564   *fPrimaryVertex += *mother;
 565
 566   return mother;
 567 }
 568
 569 //____________________________________________________________
 570 void AliHFEV0pid::AddTrackToKFVertex(AliVTrack *track, Int_t species){
 571   //
 572   // Add track to the primary vertex (if it was used in the vertex
 573   // calculation)
 574   //
 575   Bool_t isAdd = kFALSE;
 576   if(IsESDanalysis()){
 577      const AliESDVertex *esdvertex = dynamic_cast<const AliESDVertex *>(fInputEvent->GetPrimaryVertex());
 578     UShort_t *contrib = esdvertex->GetIndices();
 579     for(Int_t id = 0; id < esdvertex->GetNIndices(); id++){
 580       if(contrib[id] == track->GetID()){
 581         isAdd = kTRUE;
 582         break;
 583       }
 584     }
 585   } else {
 586     const AliAODVertex *aodvertex = dynamic_cast<const AliAODVertex *>(fInputEvent->GetPrimaryVertex());
 587     if(aodvertex->HasDaughter(track)) isAdd = kTRUE;
 588   }
 589   if(isAdd){
 590     AliKFParticle kftrack(*track, species);
 591     *fPrimaryVertex += kftrack;
 592   }
 593 }
 594
 595 //____________________________________________________________
 596 Bool_t AliHFEV0pid::CutV0(AliESDv0 *v0, Int_t V0species){
 597   //
 598   // Cut the V0
 599   //
 600   Int_t cutRequired = 0;
 601   // For cut always take min and max
 602   Double_t  cutCosPoint[2] = {0., 0.},
 603   cutDCA[2] = {0., 0.},
 604   cutProdVtx[2] = {0., 0.},
 605         cutOpAng[2] = {0., 0.},
 606         cutPsiPair[2] = {0., 0.};
 607
 608   switch(V0species){
 609   case kRecoGamma:
 610     cutCosPoint[1] = 0.03;
 611     SETBIT(cutRequired, 1);
 612     cutDCA[1] = 0.25;
 613     SETBIT(cutRequired, 3);
 614     cutProdVtx[0] = 6;
 615     SETBIT(cutRequired, 4);
 616     cutOpAng[1] = 0.1;
 617     SETBIT(cutRequired, 7);
 618     cutPsiPair[1] = 0.05;
 619     SETBIT(cutRequired, 9);
 620     break;
 621   case kRecoK0s:
 622     cutCosPoint[1] = 0.03;
 623     SETBIT(cutRequired, 1);
 624     cutDCA[1] = 0.1;
 625     SETBIT(cutRequired, 3);
 626     cutProdVtx[1] = 8.1;
 627     SETBIT(cutRequired, 5);
 628     break;
 629   case kRecoPhi:
 630     break;
 631   case kRecoLambda:
 632     cutCosPoint[1] = 0.03;
 633     SETBIT(cutRequired, 1);
 634     cutDCA[1] = 0.2;
 635     SETBIT(cutRequired, 3);
 636     cutProdVtx[1] = 24;
 637     SETBIT(cutRequired, 5);
 638     break;
 639   default:
 640     // unidentified, return
 641     return kFALSE;
 642   };
 643
 644   Char_t hname[256];
 645   const Char_t *specname[4] = {"Gamma", "K0s", ", Phi", "Lambda"};
 646   sprintf(hname, "h_cutEfficiency%s", specname[V0species-1]);
 647
 648   // Cut on pointing angle
 649   Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
 650   if(TESTBIT(cutRequired, 0) && TMath::ACos(cosPoint) < cutCosPoint[0]) return kFALSE;
 651   if(fQA) fQA->Fill(hname, 0);
 652   if(TESTBIT(cutRequired, 1) && TMath::ACos(cosPoint) > cutCosPoint[1]) return kFALSE;
 653   if(fQA) fQA->Fill(hname, 1);
 654
 655   // Cut on DCA between daughters
 656   Double_t dca = v0->GetDcaV0Daughters();
 657   if(TESTBIT(cutRequired, 2) && dca < cutDCA[0]) return kFALSE;
 658   if(fQA) fQA->Fill(hname, 2);
 659   if(TESTBIT(cutRequired, 3) && dca > cutDCA[1]) return kFALSE;
 660   if(fQA) fQA->Fill(hname, 3);
 661
 662   // Cut on reconstructed verted position
 663   Double_t x, y, z;
 664   v0->GetXYZ(x,y,z);
 665   Double_t r = TMath::Sqrt(x*x + y*y);
 666   if(TESTBIT(cutRequired, 4) && r < cutProdVtx[0]) return kFALSE;
 667   if(fQA) fQA->Fill(hname, 4);
 668   if(TESTBIT(cutRequired, 5) && r > cutProdVtx[1]) return kFALSE;
 669   if(fQA) fQA->Fill(hname, 5);
 670
 671   //Cut on Opening angle (conversions only)
 672   if(TESTBIT(cutRequired, 6) && OpenAngle(v0) < cutOpAng[0]) return kFALSE;
 673   if(fQA) fQA->Fill(hname, 6);
 674   if(TESTBIT(cutRequired, 7) && OpenAngle(v0) > cutOpAng[1]) return kFALSE;
 675   if(fQA) fQA->Fill(hname, 7);
 676
 677   //Cut on PsiPair angle (conversons only)
 678   if(TESTBIT(cutRequired, 8) && PsiPair(v0) < cutPsiPair[0]) return kFALSE;
 679   if(fQA) fQA->Fill(hname, 8);
 680   if(TESTBIT(cutRequired, 9) && PsiPair(v0) > cutPsiPair[1]) return kFALSE;
 681   if(fQA) fQA->Fill(hname, 9);
 682   return kTRUE;
 683 }
 684
 685 //____________________________________________________________
 686 Bool_t AliHFEV0pid::CutESDtrack(AliESDtrack *track){
 687   //
 688   // Single track cuts
 689   //
 690   // Hard coaded cut values for the beginning
 691   //
 692
 693   if(!track) return kFALSE;
 694
 695   // status word
 696   ULong_t status = track->GetStatus();
 697
 698   // DCA - to Vertex R & Z
 699   Float_t dcaR = -1.;
 700   Float_t dcaZ = -1.;
 701   track->GetImpactParameters(dcaR, dcaZ);
 702   //if(dcaR > 4.0) return kFALSE;
 703   //if(dcaZ > 10.0) return kFALSE;
 704
 705   // No. of TPC clusters
 706   if(track->GetTPCNcls() < 80) return kFALSE;
 707
 708   // TPC refit
 709   if(!(status & AliESDtrack::kTPCrefit)) return kFALSE;
 710
 711   // Chi2 per TPC cluster
 712   Int_t nTPCclusters = track->GetTPCclusters(0);
 713   Float_t chi2perTPCcluster = track->GetTPCchi2()/Float_t(nTPCclusters);
 714   if(chi2perTPCcluster > 3.5) return kFALSE;
 715
 716   // TPC cluster ratio
 717   Float_t cRatioTPC = track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t> (track->GetTPCNclsF()) : 1.;
 718   if(cRatioTPC < 0.6) return kFALSE;
 719
 720   // kinks
 721   if(track->GetKinkIndex(0) != 0) return kFALSE;
 722
 723   // Covariance matrix - TO BE RECONSIDERED
 724   Double_t extCov[15];
 725   track->GetExternalCovariance(extCov);
 726   //if(extCov[0]  > 2. ) return kFALSE;
 727   //if(extCov[2]  > 2. ) return kFALSE;
 728   //if(extCov[5]  > 0.5) return kFALSE;
 729   //if(extCov[9]  > 0.5) return kFALSE;
 730   //if(extCov[14] > 2. ) return kFALSE;
 731
 732   // pt
 733   if(track->Pt() < 0.1 || track->Pt() > 100) return kFALSE;
 734
 735   // eta
 736   if(TMath::Abs(track->Eta()) > 0.9) return kFALSE;
 737
 738   // the track made it through! :-)
 739   return kTRUE;
 740 }
 741
 742 //_________________________________________________
 743 Bool_t AliHFEV0pid::LooseRejectK0(AliESDv0 * const v0) const {
 744   //
 745   // Reject K0 based on loose cuts
 746   //
 747   Double_t mass = v0->GetEffMass(AliPID::kPion, AliPID::kPion);
 748   if(mass > 0.494 && mass < 0.501) return kTRUE;
 749   return kFALSE;
 750 }
 751
 752 //_________________________________________________
 753 Bool_t AliHFEV0pid::LooseRejectLambda(AliESDv0 * const v0) const {
 754   //
 755   // Reject Lambda based on loose cuts
 756   //
 757   Double_t mass1 = v0->GetEffMass(AliPID::kPion, AliPID::kProton);
 758   Double_t mass2 = v0->GetEffMass(AliPID::kProton, AliPID::kPion);
 759
 760   if(mass1 > 1.1 && mass1 < 1.12) return kTRUE;
 761   if(mass2 > 1.1 && mass2 < 1.12) return kTRUE;
 762   return kFALSE;
 763 }
 764
 765 //_________________________________________________
 766 Bool_t AliHFEV0pid::LooseRejectGamma(AliESDv0 * const v0) const {
 767   //
 768   // Reject Lambda based on loose cuts
 769   //
 770
 771   //Double_t mass = v0->GetEffMass(AliPID::kElectron, AliPID::kElectron);
 772   // DEBUG temporary solution, see the comment in GetEffMass
 773   Double_t mass = GetEffMass(v0, AliPID::kElectron, AliPID::kElectron);
 774   //..
 775   if(mass < 0.02) return kTRUE;
 776   return kFALSE;
 777 }
 778
 779 //_________________________________________________
 780 Float_t AliHFEV0pid::OpenAngle(AliESDv0 *v0) const {
 781   //
 782   // Opening angle between two daughter tracks
 783   //
 784   Double_t mn[3] = {0,0,0};
 785   Double_t mp[3] = {0,0,0};
 786
 787
 788   v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
 789   v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
 790
 791
 792   Float_t openAngle = TMath::ACos((mp[0]*mn[0] + mp[1]*mn[1] + mp[2]*mn[2])/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1] + mp[2]*mp[2])*TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1] + mn[2]*mn[2])));
 793
 794   return TMath::Abs(openAngle);
 795 }
 796
 797 //_________________________________________________
 798 Float_t AliHFEV0pid::PsiPair(AliESDv0 *v0) {
 799   //
 800   // Angle between daughter momentum plane and plane
 801   //
 802   Float_t magField = fInputEvent->GetMagneticField();
 803
 804   Int_t pIndex = v0->GetPindex();
 805   Int_t nIndex = v0->GetNindex();
 806
 807   AliESDtrack* daughter[2];
 808
 809   daughter[0] = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(pIndex));
 810   daughter[1] = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(nIndex));
 811
 812   Double_t x, y, z;
 813   v0->GetXYZ(x,y,z);//Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
 814
 815   Double_t mn[3] = {0,0,0};
 816   Double_t mp[3] = {0,0,0};
 817
 818
 819   v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
 820   v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
 821
 822
 823   Double_t deltat = 1.;
 824   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
 825
 826   Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
 827
 828   Double_t momPosProp[3];
 829   Double_t momNegProp[3];
 830
 831   AliExternalTrackParam pt(*daughter[0]), nt(*daughter[1]);
 832
 833   Float_t psiPair = 4.;
 834
 835   if(nt.PropagateTo(radiussum,magField) == 0)//propagate tracks to the outside
 836     psiPair =  -5.;
 837   if(pt.PropagateTo(radiussum,magField) == 0)
 838     psiPair = -5.;
 839   pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
 840   nt.GetPxPyPz(momNegProp);
 841
 842   Double_t pEle =
 843     TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
 844   Double_t pPos =
 845     TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
 846
 847   Double_t scalarproduct =
 848     momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
 849
 850   Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
 851
 852   psiPair =  TMath::Abs(TMath::ASin(deltat/chipair));
 853
 854   return psiPair;
 855 }
 856
 857 //____________________________________________________________
 858 AliHFEV0pid::AliHFEV0pidTrackIndex::AliHFEV0pidTrackIndex():
 859     fNElectrons(0)
 860   , fNPionsK0(0)
 861   , fNPionsL(0)
 862   , fNKaons(0)
 863   , fNProtons(0)
 864   , fIndexElectron(NULL)
 865   , fIndexPionK0(NULL)
 866   , fIndexPionL(NULL)
 867   , fIndexKaon(NULL)
 868   , fIndexProton(NULL)
 869 {
 870   //
 871   // Default Constructor
 872   //
 873 }
 874
 875 //____________________________________________________________
 876 AliHFEV0pid::AliHFEV0pidTrackIndex::~AliHFEV0pidTrackIndex(){
 877   //
 878   // Destructor
 879   //
 880   if(fIndexElectron) delete[] fIndexElectron;
 881   if(fIndexPionK0) delete[] fIndexPionK0;
 882   if(fIndexPionL) delete[] fIndexPionL;
 883   if(fIndexProton) delete[] fIndexProton;
 884 }
 885
 886 //____________________________________________________________
 887 void AliHFEV0pid::AliHFEV0pidTrackIndex::Flush(){
 888   //
 889   // Reset containers
 890   //
 891
 892   if(fIndexElectron) delete[] fIndexElectron;
 893   fIndexElectron = NULL;
 894   if(fIndexPionK0) delete[] fIndexPionK0;
 895   fIndexPionK0 = NULL;
 896   if(fIndexPionL) delete[] fIndexPionL;
 897   fIndexPionL = NULL;
 898   if(fIndexKaon) delete[] fIndexKaon;
 899   fIndexKaon = NULL;
 900   if(fIndexProton) delete[] fIndexProton;
 901   fIndexProton = NULL;
 902
 903   fNElectrons = 0;
 904   fNPionsK0 = 0;
 905   fNPionsL = 0;
 906   fNKaons = 0;
 907   fNProtons = 0;
 908 }
 909
 910 //____________________________________________________________
 911 void AliHFEV0pid::AliHFEV0pidTrackIndex::Init(Int_t capacity){
 912   //
 913   // Initialize container
 914   //
 915   fIndexElectron = new Int_t[capacity];
 916   fIndexPionK0 = new Int_t[capacity];
 917   fIndexPionL = new Int_t[capacity];
 918   fIndexProton = new Int_t[capacity];
 919 }
 920
 921 //____________________________________________________________
 922 void AliHFEV0pid::AliHFEV0pidTrackIndex::Add(Int_t index, Int_t species){
 923   //
 924   // Add new index to the list of identified particles
 925   //
 926   switch(species){
 927     case AliHFEV0pid::kRecoElectron:
 928       fIndexElectron[fNElectrons++] = index;
 929       break;
 930     case AliHFEV0pid::kRecoPionK0:
 931       fIndexPionK0[fNPionsK0++] = index;
 932       break;
 933     case AliHFEV0pid::kRecoPionL:
 934       fIndexPionL[fNPionsL++] = index;
 935       break;
 936     case AliHFEV0pid::kRecoProton:
 937       fIndexProton[fNProtons++] = index;
 938       break;
 939   };
 940 }
 941
 942 //____________________________________________________________
 943 Bool_t AliHFEV0pid::AliHFEV0pidTrackIndex::Find(Int_t index, Int_t species) const {
 944   //
 945   // Find track index in the specific sample of particles
 946   //
 947
 948   Int_t *container = NULL; Int_t n = 0;
 949   switch(species){
 950   case AliHFEV0pid::kRecoElectron:
 951     container = fIndexElectron;
 952     n = fNElectrons;
 953     break;
 954   case AliHFEV0pid::kRecoPionK0:
 955     container = fIndexPionK0;
 956     n = fNPionsK0;
 957     break;
 958   case AliHFEV0pid::kRecoPionL:
 959     container = fIndexPionL;
 960     n = fNPionsL;
 961     break;
 962   case AliHFEV0pid::kRecoProton:
 963     container = fIndexProton;
 964     n = fNProtons;
 965     break;
 966   }
 967   if(!container) return kFALSE;
 968   if(n == 0) return kFALSE;
 969   Bool_t found = kFALSE;
 970   for(Int_t i = 0; i < n; i++){
 971     if(container[i] == index){
 972       found = kTRUE;
 973       break;
 974     }
 975   }
 976   return found;
 977 }
 978
 979 //____________________________________________________________
 980 Bool_t AliHFEV0pid::AliHFEV0pidTrackIndex::Find(Int_t index) const {
 981   //
 982   // Find index in all samples
 983   //
 984   if(Find(index, AliHFEV0pid::kRecoElectron)) return kTRUE;
 985   else if(Find(index, AliHFEV0pid::kRecoPionK0)) return kTRUE;
 986   else if(Find(index, AliHFEV0pid::kRecoPionL)) return kTRUE;
 987   else return Find(index, AliHFEV0pid::kRecoProton);
 988 }
 989
 990 //____________________________________________________________
 991 AliHFEV0pid::AliHFELambdaInf::AliHFELambdaInf(AliKFParticle *mother, AliKFVertex * const primaryVertex):
 992   fInvariantMass(0.),
 993   fChi2NDF(0.),
 994   fDistVtx(0.)
 995 {
 996   //
 997   // Constructor
 998   // Fill infos
 999   //
1000   const Double_t kL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();  // PDG lambda mass
1001   fInvariantMass = mother->GetMass();
1002   // distance frm the privary vertex and mass difference for the (A) case
1003   fDistVtx = mother->GetDistanceFromVertex(*primaryVertex);
1004   // apply constraints for the fit
1005   mother->SetMassConstraint(kL0mass, 0.);
1006   mother->SetProductionVertex(*primaryVertex);
1007   fChi2NDF = mother->GetChi2()/mother->GetNDF();
1008 }
1009 //____________________________________________________________
1010 AliHFEV0pid::AliHFELambdaInf::~AliHFELambdaInf(){
1011   //
1012   // destructor
1013   //
1014 }
1015 //____________________________________________________________
1016 // DEBUG
1017 //____________________________________________________________
1018 Double_t AliHFEV0pid::GetEffMass(AliESDv0 *v0, UInt_t p1, UInt_t p2) const{
1019   //
1020   // TEMPORARY - this function should become obsolete with v4-18-Rev-10 or 11
1021   // calculate effective mass
1022   //
1023   const Double_t kpmass[5] = {TDatabasePDG::Instance()->GetParticle(kElectron)->Mass(),
1024                              TDatabasePDG::Instance()->GetParticle(kMuonMinus)->Mass(),
1025                              TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass(),
1026                              TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass(),
1027                               TDatabasePDG::Instance()->GetParticle(kProton)->Mass()}; // float
1028   if (p1>4) return -1;
1029   if (p2>4) return -1;
1030   Double_t mass1 = kpmass[p1]; // float
1031   Double_t mass2 = kpmass[p2]; // float
1032
1033   Double_t pMom[3];
1034   Double_t nMom[3];
1035
1036   v0->GetPPxPyPz(pMom[0], pMom[1], pMom[2]);
1037   v0->GetNPxPyPz(nMom[0], nMom[1], nMom[2]);
1038
1039
1040   const Double_t *m1 = pMom;
1041   const Double_t *m2 = nMom;
1042   //
1043   //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
1044   //  m1 = fPM;
1045   //  m2 = fPP;
1046   //}
1047   //
1048   Double_t e1    = TMath::Sqrt(mass1*mass1+
1049                               m1[0]*m1[0]+
1050                               m1[1]*m1[1]+
1051                               m1[2]*m1[2]); // float
1052   Double_t e2    = TMath::Sqrt(mass2*mass2+
1053                               m2[0]*m2[0]+
1054                               m2[1]*m2[1]+
1055                               m2[2]*m2[2]); // float
1056   Double_t mass =
1057     (m2[0]+m1[0])*(m2[0]+m1[0])+
1058     (m2[1]+m1[1])*(m2[1]+m1[1])+
1059     (m2[2]+m1[2])*(m2[2]+m1[2]); // float
1060
1061
1062   mass = (e1+e2)*(e1+e2)-mass;
1063   //if(mass < 0.00001){
1064   //  printf("-D: mass: %f\n", mass);
1065   // }
1066   mass = TMath::Sqrt(mass);
1067   return mass;
1068 }