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