PWGHF/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 <TObjArray.h>
  28
  29 #include "AliAODEvent.h"
  30 #include "AliAODv0.h"
  31 #include "AliESDEvent.h"
  32 #include "AliESDtrack.h"
  33 #include "AliESDv0.h"
  34 #include "AliKFVertex.h"
  35 #include "AliVEvent.h"
  36 #include "AliVTrack.h"
  37 #include "AliMCParticle.h"
  38 #include "AliStack.h"
  39 #include "AliMCEvent.h"
  40
  41 #include "AliHFEV0cuts.h"
  42 #include "AliHFEV0info.h"
  43 #include "AliHFEcollection.h"
  44
  45 #include "AliHFEV0pid.h"
  46 ClassImp(AliHFEV0pid)
  47
  48 AliHFEV0pid::AliHFEV0pid():
  49   TNamed()
  50   , fInputEvent(NULL)
  51   , fNtracks(0)
  52   , fMCEvent(NULL)
  53   , fMCon(kFALSE)
  54   , fPrimaryVertex(NULL)
  55   , fElectrons(NULL)
  56   , fPionsK0(NULL)
  57   , fPionsL(NULL)
  58   , fKaons(NULL)
  59   , fProtons(NULL)
  60   , fGammas(NULL)
  61   , fK0s(NULL)
  62   , fLambdas(NULL)
  63   , fAntiLambdas(NULL)
  64   , fIndices(NULL)
  65   , fQA(NULL)
  66   , fV0cuts(NULL)
  67   , fOutput(NULL)
  68   , fDestBits(0)
  69
  70 {
  71   //
  72   // Default constructor
  73   //
  74 }
  75 //____________________________________________________________
  76 AliHFEV0pid::AliHFEV0pid(const char *name):
  77   TNamed(name, "")
  78   , fInputEvent(NULL)
  79   , fNtracks(0)
  80   , fMCEvent(NULL)
  81   , fMCon(kFALSE)
  82   , fPrimaryVertex(NULL)
  83   , fElectrons(NULL)
  84   , fPionsK0(NULL)
  85   , fPionsL(NULL)
  86   , fKaons(NULL)
  87   , fProtons(NULL)
  88   , fGammas(NULL)
  89   , fK0s(NULL)
  90   , fLambdas(NULL)
  91   , fAntiLambdas(NULL)
  92   , fIndices(NULL)
  93   , fQA(NULL)
  94   , fV0cuts(NULL)
  95   , fOutput(NULL)
  96   , fDestBits(0)
  97 {
  98   //
  99   // Standard constructor
 100   //
 101   fElectrons = new TObjArray();
 102   fPionsK0 = new TObjArray();
 103   fPionsL = new TObjArray();
 104   fKaons = new TObjArray();
 105   fProtons = new TObjArray();
 106
 107   fElectrons->SetOwner();
 108   fPionsK0->SetOwner();
 109   fProtons->SetOwner();
 110   fPionsL->SetOwner();
 111   fKaons->SetOwner();
 112
 113   fGammas = new TObjArray();
 114   fK0s = new TObjArray();
 115   fLambdas = new TObjArray();
 116   fAntiLambdas = new TObjArray();
 117
 118   fIndices = new AliHFEV0pidTrackIndex();
 119
 120 }
 121 //____________________________________________________________
 122
 123 AliHFEV0pid::~AliHFEV0pid(){
 124   //
 125   // Destructor
 126   // Remove Containers
 127   //
 128   if(fElectrons) delete fElectrons;
 129   if(fPionsK0) delete fPionsK0;
 130   if(fPionsL) delete fPionsL;
 131   if(fKaons) delete fKaons;
 132   if(fProtons) delete fProtons;
 133
 134   if(fGammas) delete fGammas;
 135   if(fK0s) delete fK0s;
 136   if(fLambdas) delete fLambdas;
 137   if(fAntiLambdas) delete fAntiLambdas;
 138
 139   if(fIndices) delete fIndices;
 140   if(fV0cuts) delete fV0cuts;
 141
 142   if(TESTBIT(fDestBits, 1)){
 143     if(fQA) delete fQA;
 144     if(fOutput) delete fOutput;
 145   }
 146 }
 147
 148 //____________________________________________________________
 149 void AliHFEV0pid::InitQA(){
 150   //
 151   // Initialize QA histograms
 152   //
 153
 154   memset(&fDestBits, 0, sizeof(fDestBits));
 155   SETBIT(fDestBits, 1);
 156
 157   fOutput = new TList();
 158   fOutput->SetOwner();
 159
 160   fV0cuts = new AliHFEV0cuts();
 161   fV0cuts->Init("V0cuts");
 162
 163   if(!fQA){
 164     fQA = new AliHFEcollection("v0pidQA", "QA histograms for V0 selection");
 165
 166     fQA->CreateTH1F("h_nV0s", "No. of found and accepted V0s", 5, -0.5, 4.5);
 167
 168     // QA histograms for invariant mass
 169     fQA->CreateTH1F("h_InvMassGamma", "Gamma invariant mass; inv mass [GeV/c^{2}]; counts", 100, 0, 0.25);
 170     fQA->CreateTH1F("h_InvMassK0s", "K0s invariant mass; inv mass [GeV/c^{2}]; counts", 200, 0.4, 0.65);
 171     fQA->CreateTH1F("h_InvMassL", "Lambda invariant mass; inv mass [GeV/c^{2}]; counts", 100, 1.05, 1.15);
 172     fQA->CreateTH1F("h_InvMassAL", "Lambda invariant mass; inv mass [GeV/c^{2}]; counts", 100, 1.05, 1.15);
 173
 174     // QA histograms for p distribution (of the daughters)
 175     fQA->CreateTH1F("h_P_electron", "P distribution of the gamma electrons; p (GeV/c); counts", 100, 0.1, 10);
 176     fQA->CreateTH1F("h_P_K0pion", "P distribution of the K0 pions; p (GeV/c); counts", 100, 0.1, 10);
 177     fQA->CreateTH1F("h_P_Lpion", "P distribution of the Lambda pions; p (GeV/c); counts", 100, 0.1, 10);
 178     fQA->CreateTH1F("h_P_Lproton", "P distribution of the Lambda protons; p (GeV/c); counts", 100, 0.1, 10);
 179
 180     // QA pt of the V0
 181     fQA->CreateTH1F("h_Pt_Gamma", "Pt of the gamma conversion; p_{T} (GeV/c); counts", 100, 0, 10);
 182     fQA->CreateTH1F("h_Pt_K0", "Pt of the K0; p_{T} (GeV/c); counts", 100, 0, 10);
 183     fQA->CreateTH1F("h_Pt_L", "Pt of the Lambda; p_{T} (GeV/c); counts", 100, 0, 10);
 184     fQA->CreateTH1F("h_Pt_AL", "Pt of the Lambda; p_{T} (GeV/c); counts", 100, 0, 10);
 185
 186     // Armenteros plot V0 preselection
 187     fQA->CreateTH2F("h_AP_all_V0s", "armenteros plot for all V0 candidates; #alpha; Q_{T}", 200, -1, 1, 200, 0, 0.25);
 188     fQA->CreateTH2F("h_AP_selected_V0s", "armenteros plot for all V0 candidates; #alpha; Q_{T}", 200, -1, 1, 200, 0, 0.25);
 189
 190     //
 191     // !!! MC plots !!!
 192     //
 193     fQA->CreateTH2F("h_AP_MC_all_V0s", "armenteros plot for all MC tagged V0s; #alpha; Q_{T}", 200, -1, 1, 200, 0, 0.25);
 194     fQA->CreateTH2F("h_AP_MC_Gamma", "armenteros plot for MC tagged Gammas; #alpha; Q_{T}",  200, -1, 1, 200, 0, 0.25);
 195     fQA->CreateTH2F("h_AP_MC_K0", "armenteros plot for MC tagged K0s; #alpha; Q_{T}",  200, -1, 1, 200, 0, 0.25);
 196     fQA->CreateTH2F("h_AP_MC_Lambda", "armenteros plot for MC tagged Lambdas; #alpha; Q_{T}",  200, -1, 1, 200, 0, 0.25);
 197     fQA->CreateTH2F("h_AP_MC_ALambda", "armenteros plot for MC tagged A-Lambdass; #alpha; Q_{T}",  200, -1, 1, 200, 0, 0.25);
 198
 199
 200
 201     // armenteros plots for different V0 momenta - MC signal only
 202     fQA->CreateTH2Fvector1(12, "h_AP_MC_Gamma_p", "armenteros plot for MC tagged Gammas; #alpha; Q_{T}",  200, -1., 1., 200, 0., 0.25);
 203     fQA->CreateTH2Fvector1(12, "h_AP_MC_K0_p", "armenteros plot for MC tagged K0s; #alpha; Q_{T}",  200, -1., 1., 200, 0., 0.25);
 204     fQA->CreateTH2Fvector1(12, "h_AP_MC_Lambda_p", "armenteros plot for MC tagged Lambdas; #alpha; Q_{T}",  200, -1., 1., 200, 0., 0.25);
 205    //
 206
 207
 208   }
 209 }
 210
 211 //____________________________________________________________
 212 void AliHFEV0pid::Process(AliVEvent * const inputEvent){
 213
 214   //
 215   // Find protons, pions and electrons using V0 decays and
 216   // store the pointers in the TObjArray
 217   //
 218
 219
 220
 221   Int_t nGamma = 0, nK0s = 0, nLambda = 0, nPhi = 0;
 222   fInputEvent = inputEvent;
 223   fNtracks = fInputEvent->GetNumberOfTracks();
 224   fIndices->Init(fInputEvent->GetNumberOfV0s() * 2);
 225   fPrimaryVertex = new AliKFVertex(*(fInputEvent->GetPrimaryVertex()));
 226   if(!fPrimaryVertex) return;
 227   fV0cuts->SetInputEvent(fInputEvent);
 228   fV0cuts->SetPrimaryVertex(fPrimaryVertex);
 229   if(fMCEvent) fV0cuts->SetMCEvent(fMCEvent);
 230   Int_t check[fNtracks];
 231   memset(check, 0, sizeof(Int_t)*fNtracks);
 232   Int_t v0status = 0;
 233
 234   //BenchmarkV0finder();
 235
 236   for(Int_t iv0 = 0; iv0 < fInputEvent->GetNumberOfV0s(); iv0++){
 237     if(!TString(fInputEvent->IsA()->GetName()).CompareTo("AliESDEvent")){
 238       // case ESD
 239       SetESDanalysis();
 240       AliESDv0 *esdV0 = (static_cast<AliESDEvent *>(fInputEvent))->GetV0(iv0);
 241       if(!esdV0) continue;
 242       if(!esdV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
 243       v0status = ProcessV0(esdV0);
 244     } else {
 245       // case AOD
 246       SetAODanalysis();
 247       AliAODv0 *aodV0 = (static_cast<AliAODEvent *>(fInputEvent))->GetV0(iv0);
 248       if(aodV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
 249       v0status = ProcessV0(aodV0);
 250       if(AliHFEV0cuts::kUndef != v0status){
 251       }
 252     }
 253     switch(v0status){
 254     case AliHFEV0cuts::kRecoGamma: nGamma++; break;
 255     case AliHFEV0cuts::kRecoK0: nK0s++; break;
 256     case AliHFEV0cuts::kRecoPhi: nPhi++; break;
 257     case AliHFEV0cuts::kRecoLambda: nLambda++; break;
 258     };
 259   }
 260
 261   AliDebug(1, Form("Number of gammas  : %d", nGamma));
 262   AliDebug(1, Form("Number of K0s     : %d", nK0s));
 263   AliDebug(1, Form("Number of Phis    : %d", nPhi));
 264   AliDebug(1, Form("Number of Lambdas : %d", nLambda));
 265
 266   AliDebug(1, "Number of stored tracks:");
 267   AliDebug(1, Form("Number of electrons      : %d", fElectrons->GetEntries()));
 268   AliDebug(1, Form("Number of K0 pions       : %d", fPionsK0->GetEntries()));
 269   AliDebug(1, Form("Number of Lambda pions   : %d", fPionsL->GetEntries()));
 270   AliDebug(1, Form("Number of Phi kaons      : %d", fKaons->GetEntries()));
 271   AliDebug(1, Form("Number of protons        : %d", fProtons->GetEntries()));
 272
 273   delete  fPrimaryVertex;
 274
 275 }
 276
 277 //____________________________________________________________
 278 Int_t AliHFEV0pid::ProcessV0(TObject *v0){
 279   //
 280   // Process single V0
 281   // Apply general cut and special cuts for gamma, K0s, Lambda
 282   //
 283   if(!v0)  return AliHFEV0cuts::kUndef;
 284   // CHECK
 285   AliESDv0* esdV0 =  dynamic_cast<AliESDv0 *>(v0);
 286   if( ! esdV0 ) {
 287     AliError("Unexpected v0 type.");
 288     return AliHFEV0cuts::kUndef;
 289   }
 290   AliVTrack* daughter[2];
 291   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(esdV0->GetPindex()));
 292   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(esdV0->GetNindex()));
 293   if(!daughter[0] || !daughter[1]) return AliHFEV0cuts::kUndef;
 294
 295   if(fMCEvent != NULL) fMCon = kTRUE;
 296   //printf("-D: fMCEvent %x, fMCon: %i\n", fMCEvent, fMCon);
 297
 298
 299   Int_t dMC[2] = {-1, -1};
 300   Int_t idMC = AliHFEV0cuts::kUndef;
 301
 302   if(IsESDanalysis()){
 303     if(fMCon){
 304       idMC = IdentifyV0(v0, dMC);
 305       //printf("--D: V0 pid: %i, P: %i, N: %i\n", id, d[0], d[1]);
 306       fV0cuts->SetCurrentV0id(idMC);
 307       fV0cuts->SetDaughtersID(dMC);
 308     }
 309     // check common single track cuts
 310     for(Int_t i=0; i<2; ++i){
 311       if(!fV0cuts->TrackCutsCommon(static_cast<AliESDtrack*>(daughter[i]))) return AliHFEV0cuts::kUndef;
 312     }
 313     // check commom V0 cuts
 314     // CHECK
 315     if(!fV0cuts->V0CutsCommon(esdV0)) return AliHFEV0cuts::kUndef;
 316   }
 317
 318   // preselect the V0 candidates based on the Armenteros plot
 319   Int_t id = PreselectV0(esdV0, idMC);
 320   // store the resutls
 321   if(AliHFEV0cuts::kRecoGamma == id && IsGammaConv(v0)){
 322     fQA->Fill("h_nV0s", AliHFEV0cuts::kRecoGamma);
 323     return AliHFEV0cuts::kRecoGamma;
 324   }
 325   else if(AliHFEV0cuts::kRecoK0 == id && IsK0s(v0)){
 326     fQA->Fill("h_nV0s", AliHFEV0cuts::kRecoK0);
 327     return AliHFEV0cuts::kRecoK0;
 328   }
 329   else if(AliHFEV0cuts::kRecoLambda == TMath::Abs(id) && IsLambda(v0)){
 330     fQA->Fill("h_nV0s", AliHFEV0cuts::kRecoLambda);
 331     return AliHFEV0cuts::kRecoLambda;
 332   }
 333   else return AliHFEV0cuts::kUndef;
 334
 335
 336
 337 }
 338 //____________________________________________________________
 339 void AliHFEV0pid::Flush(){
 340   //
 341   // Clear the Lists
 342   //
 343   AliDebug(1, "Flushing containers");
 344   fProtons->Delete();
 345   fPionsK0->Delete();
 346   fPionsL->Delete();
 347   fElectrons->Delete();
 348   fIndices->Flush();
 349 }
 350 //____________________________________________________________
 351 Int_t AliHFEV0pid::PreselectV0(AliESDv0 * const v0, Int_t idMC){
 352   //
 353   // Based on Armenteros plot preselet the possible identity of the V0 candidate
 354   //
 355
 356   if(!v0) return -1;
 357
 358   // momentum dependent armenteros plots
 359   ArmenterosPlotMC(v0, idMC);
 360
 361   // comupte the armenteros variables
 362   Float_t ar[2];
 363   fV0cuts->Armenteros(v0, ar);
 364   // for clarity
 365   const Float_t alpha = ar[0];
 366   const Float_t qt = ar[1];
 367   //printf(" -D: Alpha: %f, QT: %f \n", alpha, qt);
 368
 369   if(TMath::Abs(alpha) > 1) return AliHFEV0cuts::kUndef;
 370
 371   fQA->Fill("h_AP_all_V0s", alpha, qt);
 372
 373   // fill a few MC tagged histograms - AP plots
 374   if(fMCEvent){
 375     switch(idMC){
 376     case AliHFEV0cuts::kRecoGamma :
 377       fQA->Fill("h_AP_MC_all_V0s", alpha, qt);
 378       fQA->Fill("h_AP_MC_Gamma", alpha, qt);
 379       break;
 380     case AliHFEV0cuts::kRecoK0 :
 381       fQA->Fill("h_AP_MC_all_V0s", alpha, qt);
 382     fQA->Fill("h_AP_MC_K0", alpha, qt);
 383       break;
 384     case AliHFEV0cuts::kRecoLambda :
 385       fQA->Fill("h_AP_MC_all_V0s", alpha, qt);
 386       fQA->Fill("h_AP_MC_Lambda", alpha, qt);
 387       break;
 388     case AliHFEV0cuts::kRecoALambda :
 389       fQA->Fill("h_AP_MC_all_V0s", alpha, qt);
 390       fQA->Fill("h_AP_MC_ALambda", alpha, qt);
 391       break;
 392     }
 393   }
 394
 395   // Gamma cuts
 396   const Double_t cutAlphaG = 0.35;
 397   const Double_t cutQTG = 0.05;
 398   const Double_t cutAlphaG2[2] = {0.6, 0.8};
 399   const Double_t cutQTG2 = 0.04;
 400
 401   // K0 cuts
 402   const Float_t cutQTK0[2] = {0.1075, 0.215};
 403   const Float_t cutAPK0[2] = {0.199, 0.8};   // parameters for curved QT cut
 404
 405   // Lambda & A-Lambda cuts
 406   const Float_t cutQTL = 0.03;
 407   const Float_t cutAlphaL[2] = {0.35, 0.7};
 408   const Float_t cutAlphaAL[2] = {-0.7,  -0.35};
 409   const Float_t cutAPL[3] = {0.107, -0.69, 0.5};  // parameters fir curved QT cut
 410
 411
 412   // Check for Gamma candidates
 413   if(qt < cutQTG){
 414     if( (TMath::Abs(alpha) < cutAlphaG) ){
 415       fQA->Fill("h_AP_selected_V0s", alpha, qt);
 416       return  AliHFEV0cuts::kRecoGamma;
 417     }
 418   }
 419   // additional region - should help high pT gammas
 420   if(qt < cutQTG2){
 421     if( (TMath::Abs(alpha) > cutAlphaG2[0]) &&  (TMath::Abs(alpha) < cutAlphaG2[1]) ){
 422       fQA->Fill("h_AP_selected_V0s", alpha, qt);
 423       return  AliHFEV0cuts::kRecoGamma;
 424     }
 425   }
 426
 427
 428   // Check for K0 candidates
 429   Float_t q = cutAPK0[0] * TMath::Sqrt(TMath::Abs(1 - alpha*alpha/(cutAPK0[1]*cutAPK0[1])));
 430   if( (qt > cutQTK0[0]) && (qt < cutQTK0[1]) && (qt > q) ){
 431     fQA->Fill("h_AP_selected_V0s", alpha, qt);
 432     return AliHFEV0cuts::kRecoK0;
 433   }
 434
 435   if( (alpha > 0) && (alpha > cutAlphaL[0])  && (alpha < cutAlphaL[1]) && (qt > cutQTL)){
 436     q = cutAPL[0] * TMath::Sqrt(1 - ( (alpha + cutAPL[1]) * (alpha + cutAPL[1]))  / (cutAPL[2]*cutAPL[2]) );
 437     if( qt < q  ){
 438       fQA->Fill("h_AP_selected_V0s", alpha, qt);
 439       return AliHFEV0cuts::kRecoLambda;
 440     }
 441   }
 442
 443   // Check for A-Lambda candidates
 444   if( (alpha < 0) && (alpha > cutAlphaAL[0]) && (alpha < cutAlphaAL[1]) && (qt > cutQTL)){
 445     q = cutAPL[0] * TMath::Sqrt(1 - ( (alpha - cutAPL[1]) * (alpha - cutAPL[1]) ) / (cutAPL[2]*cutAPL[2]) );
 446     if( qt < q ){
 447       fQA->Fill("h_AP_selected_V0s", alpha, qt);
 448       return AliHFEV0cuts::kRecoLambda;
 449     }
 450   }
 451
 452   return AliHFEV0cuts::kUndef;
 453
 454 }
 455 //____________________________________________________________
 456 Bool_t AliHFEV0pid::IsGammaConv(TObject *v0){
 457   //
 458   // Identify Gamma
 459   //
 460
 461   if(!v0) return kFALSE;
 462
 463   AliVTrack* daughter[2];
 464   Int_t pIndex = 0, nIndex = 0;
 465   Double_t invMass = 0.;
 466   Double_t mPt = 0.;
 467   Int_t v0id = -1;
 468   if(IsESDanalysis()){
 469     // ESD - cut V0
 470     AliESDv0 *esdV0 = static_cast<AliESDv0 *>(v0);
 471     v0id = esdV0->GetLabel();
 472     // apply FULL gamma cuts
 473     if(!fV0cuts->GammaCuts(esdV0)) return kFALSE;
 474     invMass = esdV0->GetEffMass(AliPID::kElectron, AliPID::kElectron);
 475     pIndex = esdV0->GetPindex();
 476     nIndex = esdV0->GetNindex();
 477     mPt = esdV0->Pt();
 478   } else {
 479     // AOD Analysis - not possible to cut
 480     AliAODv0 *aodV0 = static_cast<AliAODv0 *>(v0);
 481     v0id = aodV0->GetID();
 482     pIndex = aodV0->GetPosID();
 483     nIndex = aodV0->GetNegID();
 484     invMass = aodV0->InvMass2Prongs(0, 1, kElectron, kElectron);
 485   }
 486   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
 487   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
 488   if(!daughter[0] || !daughter[1]) return kFALSE;
 489
 490   // DEBUG
 491   AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
 492
 493   // AFTER all gamma cuts
 494   fQA->Fill("h_Pt_Gamma", mPt);
 495   fQA->Fill("h_InvMassGamma", invMass);
 496
 497   // Identified gamma - store tracks in the electron containers
 498   if(!fIndices->Find(daughter[0]->GetID())){
 499     AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
 500     fElectrons->Add(new AliHFEV0info(daughter[0], daughter[1]->GetID(), v0id));
 501     fIndices->Add(daughter[0]->GetID(), AliHFEV0cuts::kRecoElectron);
 502   }
 503   if(!fIndices->Find(daughter[1]->GetID())){
 504     AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[1]->GetID(), daughter[1]->GetID()));
 505     fElectrons->Add(new AliHFEV0info(daughter[1], daughter[0]->GetID(), v0id));
 506     fIndices->Add(daughter[1]->GetID(), AliHFEV0cuts::kRecoElectron);
 507   }
 508   fGammas->Add(v0);
 509
 510   return kTRUE;
 511 }
 512 //____________________________________________________________
 513 Bool_t AliHFEV0pid::IsK0s(TObject *v0){
 514   //
 515   // Identify K0s
 516   //
 517
 518   if(!v0) return kFALSE;
 519
 520   AliVTrack* daughter[2];
 521   Int_t pIndex = 0, nIndex = 0;
 522   Int_t v0id = -1;
 523   Double_t invMass = 0.;
 524   Double_t mPt = 0.;
 525   if(IsESDanalysis()){
 526     // ESD - cut V0
 527     AliESDv0 *esdV0 = static_cast<AliESDv0 *>(v0);
 528     if(!fV0cuts->K0Cuts(esdV0)) return kFALSE;
 529     v0id = esdV0->GetLabel();
 530     pIndex = esdV0->GetPindex();
 531     nIndex = esdV0->GetNindex();
 532     invMass = esdV0->GetEffMass(AliPID::kPion, AliPID::kPion);
 533     mPt = esdV0->Pt();
 534   } else {
 535     // AOD Analysis - not possible to cut
 536     AliAODv0 *aodV0 = static_cast<AliAODv0 *>(v0);
 537     aodV0->GetID();
 538     pIndex = aodV0->GetPosID();
 539     nIndex = aodV0->GetNegID();
 540     invMass = aodV0->MassK0Short();
 541   }
 542   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
 543   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
 544   if(!daughter[0] || !daughter[1]) return kFALSE;
 545
 546   fQA->Fill("h_Pt_K0", mPt);
 547   fQA->Fill("h_InvMassK0s", invMass);
 548
 549   AliDebug(1, Form("K0 identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
 550
 551   // AFTER all K0 cuts
 552   // Identified gamma - store tracks in the electron containers
 553   if(!fIndices->Find(daughter[0]->GetID())){
 554     AliDebug(1, Form("Adding K0 Pion track with ID %d", daughter[0]->GetID()));
 555     fPionsK0->Add(new AliHFEV0info(daughter[0], daughter[1]->GetID(), v0id));
 556     fIndices->Add(daughter[0]->GetID(), AliHFEV0cuts::kRecoPionK0);
 557   }
 558   if(!fIndices->Find(daughter[1]->GetID())){
 559     AliDebug(1, Form("Adding K0 Pion track with ID %d", daughter[1]->GetID()));
 560     fPionsK0->Add(new AliHFEV0info(daughter[1], daughter[0]->GetID(), v0id));
 561     fIndices->Add(daughter[1]->GetID(), AliHFEV0cuts::kRecoPionK0);
 562   }
 563   fK0s->Add(v0);
 564   return kTRUE;
 565 }
 566
 567 //____________________________________________________________
 568 Bool_t AliHFEV0pid::IsPhi(const TObject *v0) const {
 569   //
 570   // Identify Phi - very preliminary - requires diffrent approach (V0 fnder is not effective)
 571   //
 572
 573   //const Double_t kPhiMass=TDatabasePDG::Instance()->GetParticle(333)->Mass();  // PDG phi mass
 574   //AliVTrack* daughter[2];
 575   //Double_t invMass = 0.;
 576
 577   if(!v0) return kFALSE;
 578
 579   //Int_t pIndex = 0, nIndex = 0;
 580   if(IsESDanalysis()){
 581     // ESD - cut V0
 582     //AliESDv0 *esdV0 = static_cast<AliESDv0 *>(v0);
 583     //pIndex = esdV0->GetPindex();
 584     //nIndex = esdV0->GetNindex();
 585   } else {
 586     // PRELIMINARY - !!!
 587     // AOD Analysis - not possible to cut
 588   }
 589
 590   return kTRUE;
 591 }
 592
 593 //____________________________________________________________
 594 Bool_t AliHFEV0pid::IsLambda(TObject *v0){
 595   //
 596   // Identify Lambda
 597   //
 598
 599   if(!v0) return kFALSE;
 600
 601   AliVTrack* daughter[2];
 602   Int_t pIndex = 0, nIndex = 0;
 603   Double_t invMass = 0.;
 604   Bool_t isLambda = kTRUE; // Lambda - kTRUE, Anti Lambda - kFALSE
 605   Double_t mPt = 0.;
 606   Int_t v0id = -1;
 607   if(IsESDanalysis()){
 608     // ESD - cut V0
 609     AliESDv0 *esdV0 = static_cast<AliESDv0 *>(v0);
 610     v0id = esdV0->GetLabel();
 611     if(!fV0cuts->LambdaCuts(esdV0,isLambda)) return kFALSE;
 612     mPt = esdV0->Pt();
 613     if(fV0cuts->CheckSigns(esdV0)){
 614       pIndex = esdV0->GetPindex();
 615       nIndex = esdV0->GetNindex();
 616     }
 617     else{
 618       pIndex = esdV0->GetNindex();
 619       nIndex = esdV0->GetPindex();
 620     }
 621   } else {
 622     // PRELIMINARY - !!!
 623     // AOD Analysis - not possible to cut
 624
 625     // again - two cases as above
 626     AliAODv0 *aodV0 = static_cast<AliAODv0 *>(v0);
 627     v0id = aodV0->GetID();
 628     pIndex = aodV0->GetPosID();
 629     nIndex = aodV0->GetNegID();
 630     invMass = aodV0->MassLambda();
 631   }
 632
 633   daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
 634   daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
 635   if(!daughter[0] || !daughter[1]) return kFALSE;
 636
 637   // lambda
 638   if(isLambda){
 639     fQA->Fill("h_Pt_L", mPt);
 640     fQA->Fill("h_InvMassL", invMass);
 641
 642     if(!fIndices->Find(daughter[0]->GetID())){
 643       fProtons->Add(new AliHFEV0info(daughter[0], daughter[1]->GetID(), v0id));
 644       fIndices->Add(daughter[0]->GetID(), AliHFEV0cuts::kRecoProton);
 645     }
 646     if(!fIndices->Find(daughter[1]->GetID())){
 647       fPionsL->Add(new AliHFEV0info(daughter[1], daughter[0]->GetID(), v0id));
 648       fIndices->Add(daughter[1]->GetID(), AliHFEV0cuts::kRecoPionL);
 649     }
 650   }
 651   // antilambda
 652   else{
 653     fQA->Fill("h_Pt_AL", mPt);
 654     fQA->Fill("h_InvMassAL", invMass);
 655     if(!fIndices->Find(daughter [1]->GetID())){
 656       fProtons->Add(new AliHFEV0info(daughter[1], daughter[0]->GetID(), v0id));
 657       fIndices->Add(daughter[1]->GetID(), AliHFEV0cuts::kRecoProton);
 658     }
 659     if(!fIndices->Find(daughter [0]->GetID())){
 660       fPionsL->Add(new AliHFEV0info(daughter[0], daughter[1]->GetID(), v0id));
 661       fIndices->Add(daughter [0]->GetID(), AliHFEV0cuts::kRecoPionL);
 662     }
 663   }
 664   if(isLambda) fLambdas->Add(v0);
 665   else fAntiLambdas->Add(v0);
 666
 667   return kTRUE;
 668 }
 669
 670 //____________________________________________________________
 671 Int_t AliHFEV0pid::IdentifyV0(TObject *esdV0, Int_t d[2]){
 672   //
 673   // for MC only, returns the V0 Id
 674   //
 675
 676   //
 677   // be carefull about changing the return values - they are used later selectively
 678   // In particulra "-2" means that identity of either of the daughters could not be
 679   // estimated
 680   //
 681
 682   AliESDv0 *v0 = dynamic_cast<AliESDv0 *>(esdV0);
 683
 684   if(!v0) return -1;
 685   AliESDtrack* dN, *dP;
 686   Int_t iN, iP;
 687   iN = iP = -1;
 688   iP = v0->GetPindex();
 689   iN = v0->GetNindex();
 690   if(iN < 0 || iP < 0) return -1;
 691   if(iN >= fNtracks || iP >= fNtracks) return -1;
 692   dP = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(iP));
 693   dN = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(iN));
 694   if(!dN || !dP) return -1;
 695
 696   // as of 26/10/2010
 697   // there is still a problem with wrong assignment of positive and negative
 698   // V0 daughter in V0 finder - a check is necessary
 699   // if the V0 daughters are miss-assigned - swap their labels
 700   Bool_t sign = fV0cuts->CheckSigns(v0);
 701
 702   // get the MC labels
 703   Int_t lN, lP;
 704   if(sign){
 705     lN = dN->GetLabel();
 706     lP = dP->GetLabel();
 707   }
 708   else{
 709     lP = dN->GetLabel();
 710     lN = dP->GetLabel();
 711   }
 712   if(lN < 0 || lP < 0) return -2;
 713   // get the associated MC particles
 714   AliMCParticle *mcP, *mcN;
 715   mcP = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(lP));
 716   mcN = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(lN));
 717   if(!mcP || !mcN) return -2;
 718
 719   // identify the daughter tracks and their mothers
 720   Int_t pdgP, pdgN;
 721   pdgP = TMath::Abs(mcP->PdgCode());
 722   pdgN = TMath::Abs(mcN->PdgCode());
 723   // store the daughter ID for later use
 724   d[0] = pdgP;
 725   d[1] = pdgN;
 726   //printf(" -D: pdgP: %i, pdgN: %i\n", pdgP, pdgN);
 727   // lablel of the mother particle
 728   // -1 may mean it was a primary particle
 729   Int_t lPm, lNm;
 730   lPm = mcP->GetMother();
 731   lNm = mcN->GetMother();
 732   if(-1==lPm || -1==lNm) return -3;
 733
 734   // return if mothers are not the same particle
 735   if(lPm != lNm) return -3;
 736   // get MC mother particle - now we need only 1
 737   AliMCParticle *m = dynamic_cast<AliMCParticle*>(fMCEvent->GetTrack(lPm));
 738   if(!m) return -2;
 739   // mother PDG
 740   Int_t pdgM = m->PdgCode();
 741
 742   //   if(3122 == TMath::Abs(pdgM)){
 743   //     printf("-D: v0 signs: %i\n", fV0cuts->CheckSigns(v0));
 744   //     printf("-D: pdgM: %i, pdgN: %i, pdgP: %i \n", pdgM, pdgN, pdgP);
 745   //   }
 746
 747
 748   // now check the mother and daughters identity
 749   if(22 == TMath::Abs(pdgM) && 11 == pdgN && 11 == pdgP) return AliHFEV0cuts::kRecoGamma;
 750   if(310 == TMath::Abs(pdgM) && 211 == pdgN && 211 == pdgP) return AliHFEV0cuts::kRecoK0;
 751   if(-3122 == pdgM && 2212 == pdgN && 211 == pdgP) return AliHFEV0cuts::kRecoALambda;
 752   if(3122 == pdgM && 211 == pdgN && 2212 == pdgP) return AliHFEV0cuts::kRecoLambda;
 753
 754   return AliHFEV0cuts::kUndef;
 755
 756 }
 757 //____________________________________________________________
 758 void AliHFEV0pid::BenchmarkV0finder(){
 759   //
 760   // produce histograms for all findable V0s that are
 761   // were selected byt the (oline) V0 finder and can
 762   // be used to estimate the efficiency of teh V0 cuts
 763   //
 764
 765   for(Int_t iv0 = 0; iv0 < fInputEvent->GetNumberOfV0s(); iv0++){
 766     AliESDv0 *esdV0 = (static_cast<AliESDEvent *>(fInputEvent))->GetV0(iv0);
 767     if(!esdV0) continue;
 768     if(!esdV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
 769     // indetify the V0 candidate
 770     Int_t idV0 = AliHFEV0cuts::kUndef;
 771     Int_t idD[2] = {-1, -1};
 772     idV0 = IdentifyV0(esdV0, idD);
 773   }
 774 }
 775 //____________________________________________________________
 776 void   AliHFEV0pid::ArmenterosPlotMC(AliESDv0 * const v0, Int_t idMC){
 777   //
 778   // Armenteros plots as a function of Mohter Momentum
 779   //
 780   //const Float_t minP = 0.1;
 781   //const Float_t maxP = 10.;
 782   // approx log bins - over the 0.1 - 10 GeV/c
 783   const Float_t bins[13] = {0.1, 0.1468, 0.2154, 0.3162, 0.4642, 0.6813, 1.0, 1.4678, 2.1544, 3.1623, 4.6416, 6.8129, 10.0};
 784
 785   Float_t ar[2];
 786   fV0cuts->Armenteros(v0, ar);
 787   Float_t p = v0->P();
 788
 789   if( (p <=  bins[0]) || (p >= bins[12])) return;
 790
 791   Int_t pBin = 0;
 792   Float_t tmp = bins[0];
 793   while(tmp < p){
 794     ++pBin;
 795     tmp = bins[pBin];
 796   }
 797   pBin--;
 798
 799   if(AliHFEV0cuts::kRecoGamma == idMC) fQA->Fill("h_AP_MC_Gamma_p", pBin, ar[0], ar[1]);
 800   if(AliHFEV0cuts::kRecoK0 == idMC) fQA->Fill("h_AP_MC_K0_p", pBin, ar[0], ar[1]);
 801   if(AliHFEV0cuts::kRecoLambda == TMath::Abs(idMC)) fQA->Fill("h_AP_MC_Lambda_p", pBin, ar[0], ar[1]);
 802
 803
 804
 805 }
 806 //____________________________________________________________
 807 AliHFEV0pid::AliHFEV0pidTrackIndex::AliHFEV0pidTrackIndex():
 808     fNElectrons(0)
 809   , fNPionsK0(0)
 810   , fNPionsL(0)
 811   , fNKaons(0)
 812   , fNProtons(0)
 813   , fIndexElectron(NULL)
 814   , fIndexPionK0(NULL)
 815   , fIndexPionL(NULL)
 816   , fIndexKaon(NULL)
 817   , fIndexProton(NULL)
 818 {
 819   //
 820   // Default Constructor
 821   //
 822 }
 823
 824 //____________________________________________________________
 825 AliHFEV0pid::AliHFEV0pidTrackIndex::~AliHFEV0pidTrackIndex(){
 826   //
 827   // Destructor
 828   //
 829   if(fIndexElectron) delete[] fIndexElectron;
 830   if(fIndexPionK0) delete[] fIndexPionK0;
 831   if(fIndexPionL) delete[] fIndexPionL;
 832   if(fIndexProton) delete[] fIndexProton;
 833 }
 834
 835 //____________________________________________________________
 836 void AliHFEV0pid::AliHFEV0pidTrackIndex::Flush(){
 837   //
 838   // Reset containers
 839   //
 840
 841   if(fIndexElectron) delete[] fIndexElectron;
 842   fIndexElectron = NULL;
 843   if(fIndexPionK0) delete[] fIndexPionK0;
 844   fIndexPionK0 = NULL;
 845   if(fIndexPionL) delete[] fIndexPionL;
 846   fIndexPionL = NULL;
 847   if(fIndexKaon) delete[] fIndexKaon;
 848   fIndexKaon = NULL;
 849   if(fIndexProton) delete[] fIndexProton;
 850   fIndexProton = NULL;
 851
 852   fNElectrons = 0;
 853   fNPionsK0 = 0;
 854   fNPionsL = 0;
 855   fNKaons = 0;
 856   fNProtons = 0;
 857 }
 858
 859 //____________________________________________________________
 860 void AliHFEV0pid::AliHFEV0pidTrackIndex::Init(Int_t capacity){
 861   //
 862   // Initialize container
 863   //
 864   fIndexElectron = new Int_t[capacity];
 865   fIndexPionK0 = new Int_t[capacity];
 866   fIndexPionL = new Int_t[capacity];
 867   fIndexProton = new Int_t[capacity];
 868 }
 869
 870 //____________________________________________________________
 871 void AliHFEV0pid::AliHFEV0pidTrackIndex::Add(Int_t index, Int_t species){
 872   //
 873   // Add new index to the list of identified particles
 874   //
 875   switch(species){
 876     case AliHFEV0cuts::kRecoElectron:
 877       fIndexElectron[fNElectrons++] = index;
 878       break;
 879     case AliHFEV0cuts::kRecoPionK0:
 880       fIndexPionK0[fNPionsK0++] = index;
 881       break;
 882     case AliHFEV0cuts::kRecoPionL:
 883       fIndexPionL[fNPionsL++] = index;
 884       break;
 885     case AliHFEV0cuts::kRecoProton:
 886       fIndexProton[fNProtons++] = index;
 887       break;
 888   };
 889 }
 890
 891 //____________________________________________________________
 892 Bool_t AliHFEV0pid::AliHFEV0pidTrackIndex::Find(Int_t index, Int_t species) const {
 893   //
 894   // Find track index in the specific sample of particles
 895   //
 896
 897   Int_t *container = NULL; Int_t n = 0;
 898   switch(species){
 899   case AliHFEV0cuts::kRecoElectron:
 900     container = fIndexElectron;
 901     n = fNElectrons;
 902     break;
 903   case AliHFEV0cuts::kRecoPionK0:
 904     container = fIndexPionK0;
 905     n = fNPionsK0;
 906     break;
 907   case AliHFEV0cuts::kRecoPionL:
 908     container = fIndexPionL;
 909     n = fNPionsL;
 910     break;
 911   case AliHFEV0cuts::kRecoProton:
 912     container = fIndexProton;
 913     n = fNProtons;
 914     break;
 915   }
 916   if(!container) return kFALSE;
 917   if(n == 0) return kFALSE;
 918   Bool_t found = kFALSE;
 919   for(Int_t i = 0; i < n; i++){
 920     if(container[i] == index){
 921       found = kTRUE;
 922       break;
 923     }
 924   }
 925   return found;
 926 }
 927
 928 //____________________________________________________________
 929 Bool_t AliHFEV0pid::AliHFEV0pidTrackIndex::Find(Int_t index) const {
 930   //
 931   // Find index in all samples
 932   //
 933   if(Find(index, AliHFEV0cuts::kRecoElectron)) return kTRUE;
 934   else if(Find(index, AliHFEV0cuts::kRecoPionK0)) return kTRUE;
 935   else if(Find(index, AliHFEV0cuts::kRecoPionL)) return kTRUE;
 936   else return Find(index, AliHFEV0cuts::kRecoProton);
 937 }
 938
 939 //____________________________________________________________
 940 TList *AliHFEV0pid::GetListOfQAhistograms(){
 941   //
 942   // Getter for V0 PID QA histograms
 943   //
 944
 945   CLRBIT(fDestBits, 1);
 946
 947   TList *tmp = fV0cuts->GetList();
 948   tmp->SetName("V0cuts");
 949   fOutput->Add(tmp);
 950   if(fQA){
 951     tmp = 0x0;
 952     tmp = fQA->GetList();
 953     tmp->SetName("V0pid");
 954     fOutput->Add(tmp);
 955   }
 956   tmp = 0x0;
 957   tmp = fV0cuts->GetListMC();
 958   tmp->SetName("V0cutsMC");
 959   //printf(" -D: adding MC V0 cuts stuff\n");
 960   fOutput->Add(tmp);
 961
 962   return fOutput;
 963 }