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