Update of the HFE package
[u/mrichter/AliRoot.git] / 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 <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(TObject *v0){
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 const 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 }