]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/hfe/AliHFEV0pid.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates for TRD HFE analysis
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEV0pid.cxx
CommitLineData
70da6c5a 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//
70da6c5a 27#include <TObjArray.h>
70da6c5a 28
29#include "AliAODEvent.h"
70da6c5a 30#include "AliAODv0.h"
70da6c5a 31#include "AliESDEvent.h"
32#include "AliESDtrack.h"
33#include "AliESDv0.h"
70da6c5a 34#include "AliKFVertex.h"
35#include "AliVEvent.h"
36#include "AliVTrack.h"
3a72645a 37#include "AliMCParticle.h"
38#include "AliStack.h"
39#include "AliMCEvent.h"
70da6c5a 40
faee3b18 41#include "AliHFEV0cuts.h"
42#include "AliHFEV0info.h"
43#include "AliHFEcollection.h"
44
70da6c5a 45#include "AliHFEV0pid.h"
46ClassImp(AliHFEV0pid)
47
70da6c5a 48AliHFEV0pid::AliHFEV0pid():
e156c3bb 49 TNamed()
70da6c5a 50 , fInputEvent(NULL)
3a72645a 51 , fNtracks(0)
52 , fMCEvent(NULL)
53 , fMCon(kFALSE)
70da6c5a 54 , fPrimaryVertex(NULL)
55 , fElectrons(NULL)
56 , fPionsK0(NULL)
57 , fPionsL(NULL)
58 , fKaons(NULL)
59 , fProtons(NULL)
faee3b18 60 , fGammas(NULL)
61 , fK0s(NULL)
62 , fLambdas(NULL)
63 , fAntiLambdas(NULL)
70da6c5a 64 , fIndices(NULL)
65 , fQA(NULL)
faee3b18 66 , fV0cuts(NULL)
67 , fOutput(NULL)
e156c3bb 68 , fDestBits(0)
69
70da6c5a 70{
71 //
72 // Default constructor
73 //
e156c3bb 74}
75//____________________________________________________________
76AliHFEV0pid::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 //
70da6c5a 101 fElectrons = new TObjArray();
102 fPionsK0 = new TObjArray();
103 fPionsL = new TObjArray();
104 fKaons = new TObjArray();
105 fProtons = new TObjArray();
faee3b18 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
70da6c5a 120}
70da6c5a 121//____________________________________________________________
e156c3bb 122
70da6c5a 123AliHFEV0pid::~AliHFEV0pid(){
124 //
125 // Destructor
126 // Remove Containers
127 //
70da6c5a 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;
faee3b18 133
134 if(fGammas) delete fGammas;
135 if(fK0s) delete fK0s;
136 if(fLambdas) delete fLambdas;
137 if(fAntiLambdas) delete fAntiLambdas;
138
70da6c5a 139 if(fIndices) delete fIndices;
faee3b18 140 if(fV0cuts) delete fV0cuts;
e156c3bb 141
142 if(TESTBIT(fDestBits, 1)){
143 if(fQA) delete fQA;
144 if(fOutput) delete fOutput;
145 }
70da6c5a 146}
147
148//____________________________________________________________
149void AliHFEV0pid::InitQA(){
150 //
151 // Initialize QA histograms
152 //
faee3b18 153
e156c3bb 154 memset(&fDestBits, 0, sizeof(fDestBits));
155 SETBIT(fDestBits, 1);
156
faee3b18 157 fOutput = new TList();
e156c3bb 158 fOutput->SetOwner();
faee3b18 159
160 fV0cuts = new AliHFEV0cuts();
161 fV0cuts->Init("V0cuts");
162
70da6c5a 163 if(!fQA){
3a72645a 164 fQA = new AliHFEcollection("v0pidQA", "QA histograms for V0 selection");
70da6c5a 165
faee3b18 166 fQA->CreateTH1F("h_nV0s", "No. of found and accepted V0s", 5, -0.5, 4.5);
70da6c5a 167
168 // QA histograms for invariant mass
91c7e1ec 169 fQA->CreateTH1F("h_InvMassGamma", "Gamma invariant mass; inv mass [GeV/c^{2}]; counts", 100, 0, 0.25);
faee3b18 170 fQA->CreateTH1F("h_InvMassK0s", "K0s invariant mass; inv mass [GeV/c^{2}]; counts", 200, 0.4, 0.65);
3a72645a 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);
70da6c5a 173
91c7e1ec 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
70da6c5a 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);
3a72645a 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);
faee3b18 185
3a72645a 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
70da6c5a 207
70da6c5a 208 }
209}
210
211//____________________________________________________________
212void 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 //
faee3b18 218
e97c2edf 219
220
70da6c5a 221 Int_t nGamma = 0, nK0s = 0, nLambda = 0, nPhi = 0;
222 fInputEvent = inputEvent;
3a72645a 223 fNtracks = fInputEvent->GetNumberOfTracks();
70da6c5a 224 fIndices->Init(fInputEvent->GetNumberOfV0s() * 2);
225 fPrimaryVertex = new AliKFVertex(*(fInputEvent->GetPrimaryVertex()));
faee3b18 226 if(!fPrimaryVertex) return;
227 fV0cuts->SetInputEvent(fInputEvent);
228 fV0cuts->SetPrimaryVertex(fPrimaryVertex);
3a72645a 229 if(fMCEvent) fV0cuts->SetMCEvent(fMCEvent);
e97c2edf 230 Int_t check[fNtracks];
231 memset(check, 0, sizeof(Int_t)*fNtracks);
70da6c5a 232 Int_t v0status = 0;
e97c2edf 233
234 //BenchmarkV0finder();
235
70da6c5a 236 for(Int_t iv0 = 0; iv0 < fInputEvent->GetNumberOfV0s(); iv0++){
237 if(!TString(fInputEvent->IsA()->GetName()).CompareTo("AliESDEvent")){
238 // case ESD
239 SetESDanalysis();
bf892a6a 240 AliESDv0 *esdV0 = (static_cast<AliESDEvent *>(fInputEvent))->GetV0(iv0);
faee3b18 241 if(!esdV0) continue;
70da6c5a 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();
bf892a6a 247 AliAODv0 *aodV0 = (static_cast<AliAODEvent *>(fInputEvent))->GetV0(iv0);
70da6c5a 248 if(aodV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
249 v0status = ProcessV0(aodV0);
3a72645a 250 if(AliHFEV0cuts::kUndef != v0status){
faee3b18 251 }
70da6c5a 252 }
253 switch(v0status){
3a72645a 254 case AliHFEV0cuts::kRecoGamma: nGamma++; break;
255 case AliHFEV0cuts::kRecoK0: nK0s++; break;
256 case AliHFEV0cuts::kRecoPhi: nPhi++; break;
257 case AliHFEV0cuts::kRecoLambda: nLambda++; break;
70da6c5a 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;
c2690925 274
70da6c5a 275}
276
277//____________________________________________________________
278Int_t AliHFEV0pid::ProcessV0(TObject *v0){
279 //
280 // Process single V0
281 // Apply general cut and special cuts for gamma, K0s, Lambda
282 //
bf892a6a 283 if(!v0) return AliHFEV0cuts::kUndef;
215ffe88 284 // CHECK
285 AliESDv0* esdV0 = dynamic_cast<AliESDv0 *>(v0);
286 if( ! esdV0 ) {
287 AliError("Unexpected v0 type.");
288 return AliHFEV0cuts::kUndef;
289 }
70da6c5a 290 AliVTrack* daughter[2];
215ffe88 291 daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(esdV0->GetPindex()));
292 daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(esdV0->GetNindex()));
3a72645a 293 if(!daughter[0] || !daughter[1]) return AliHFEV0cuts::kUndef;
294
bf892a6a 295 if(fMCEvent != NULL) fMCon = kTRUE;
3a72645a 296 //printf("-D: fMCEvent %x, fMCon: %i\n", fMCEvent, fMCon);
297
e97c2edf 298
3a72645a 299 Int_t dMC[2] = {-1, -1};
300 Int_t idMC = AliHFEV0cuts::kUndef;
70da6c5a 301
302 if(IsESDanalysis()){
3a72645a 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
70da6c5a 310 for(Int_t i=0; i<2; ++i){
bf892a6a 311 if(!fV0cuts->TrackCutsCommon(static_cast<AliESDtrack*>(daughter[i]))) return AliHFEV0cuts::kUndef;
3a72645a 312 }
faee3b18 313 // check commom V0 cuts
215ffe88 314 // CHECK
315 if(!fV0cuts->V0CutsCommon(esdV0)) return AliHFEV0cuts::kUndef;
70da6c5a 316 }
317
3a72645a 318 // preselect the V0 candidates based on the Armenteros plot
215ffe88 319 Int_t id = PreselectV0(esdV0, idMC);
faee3b18 320 // store the resutls
3a72645a 321 if(AliHFEV0cuts::kRecoGamma == id && IsGammaConv(v0)){
322 fQA->Fill("h_nV0s", AliHFEV0cuts::kRecoGamma);
323 return AliHFEV0cuts::kRecoGamma;
faee3b18 324 }
3a72645a 325 else if(AliHFEV0cuts::kRecoK0 == id && IsK0s(v0)){
326 fQA->Fill("h_nV0s", AliHFEV0cuts::kRecoK0);
327 return AliHFEV0cuts::kRecoK0;
faee3b18 328 }
3a72645a 329 else if(AliHFEV0cuts::kRecoLambda == TMath::Abs(id) && IsLambda(v0)){
330 fQA->Fill("h_nV0s", AliHFEV0cuts::kRecoLambda);
331 return AliHFEV0cuts::kRecoLambda;
faee3b18 332 }
c2690925 333 else return AliHFEV0cuts::kUndef;
334
335
faee3b18 336
70da6c5a 337}
338//____________________________________________________________
339void AliHFEV0pid::Flush(){
340 //
341 // Clear the Lists
342 //
343 AliDebug(1, "Flushing containers");
faee3b18 344 fProtons->Delete();
345 fPionsK0->Delete();
346 fPionsL->Delete();
347 fElectrons->Delete();
70da6c5a 348 fIndices->Flush();
349}
70da6c5a 350//____________________________________________________________
3a72645a 351Int_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//____________________________________________________________
70da6c5a 456Bool_t AliHFEV0pid::IsGammaConv(TObject *v0){
457 //
458 // Identify Gamma
459 //
bf892a6a 460
461 if(!v0) return kFALSE;
462
70da6c5a 463 AliVTrack* daughter[2];
464 Int_t pIndex = 0, nIndex = 0;
465 Double_t invMass = 0.;
faee3b18 466 Double_t mPt = 0.;
467 Int_t v0id = -1;
70da6c5a 468 if(IsESDanalysis()){
469 // ESD - cut V0
bf892a6a 470 AliESDv0 *esdV0 = static_cast<AliESDv0 *>(v0);
faee3b18 471 v0id = esdV0->GetLabel();
472 // apply FULL gamma cuts
473 if(!fV0cuts->GammaCuts(esdV0)) return kFALSE;
474 invMass = esdV0->GetEffMass(AliPID::kElectron, AliPID::kElectron);
70da6c5a 475 pIndex = esdV0->GetPindex();
476 nIndex = esdV0->GetNindex();
faee3b18 477 mPt = esdV0->Pt();
70da6c5a 478 } else {
479 // AOD Analysis - not possible to cut
bf892a6a 480 AliAODv0 *aodV0 = static_cast<AliAODv0 *>(v0);
faee3b18 481 v0id = aodV0->GetID();
70da6c5a 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;
faee3b18 489
490 // DEBUG
491 AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
70da6c5a 492
faee3b18 493 // AFTER all gamma cuts
70da6c5a 494 fQA->Fill("h_Pt_Gamma", mPt);
faee3b18 495 fQA->Fill("h_InvMassGamma", invMass);
70da6c5a 496
70da6c5a 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()));
faee3b18 500 fElectrons->Add(new AliHFEV0info(daughter[0], daughter[1]->GetID(), v0id));
3a72645a 501 fIndices->Add(daughter[0]->GetID(), AliHFEV0cuts::kRecoElectron);
70da6c5a 502 }
503 if(!fIndices->Find(daughter[1]->GetID())){
504 AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[1]->GetID(), daughter[1]->GetID()));
faee3b18 505 fElectrons->Add(new AliHFEV0info(daughter[1], daughter[0]->GetID(), v0id));
3a72645a 506 fIndices->Add(daughter[1]->GetID(), AliHFEV0cuts::kRecoElectron);
70da6c5a 507 }
faee3b18 508 fGammas->Add(v0);
509
70da6c5a 510 return kTRUE;
511}
70da6c5a 512//____________________________________________________________
513Bool_t AliHFEV0pid::IsK0s(TObject *v0){
514 //
515 // Identify K0s
516 //
bf892a6a 517
518 if(!v0) return kFALSE;
519
70da6c5a 520 AliVTrack* daughter[2];
521 Int_t pIndex = 0, nIndex = 0;
faee3b18 522 Int_t v0id = -1;
70da6c5a 523 Double_t invMass = 0.;
faee3b18 524 Double_t mPt = 0.;
70da6c5a 525 if(IsESDanalysis()){
526 // ESD - cut V0
bf892a6a 527 AliESDv0 *esdV0 = static_cast<AliESDv0 *>(v0);
faee3b18 528 if(!fV0cuts->K0Cuts(esdV0)) return kFALSE;
529 v0id = esdV0->GetLabel();
70da6c5a 530 pIndex = esdV0->GetPindex();
531 nIndex = esdV0->GetNindex();
faee3b18 532 invMass = esdV0->GetEffMass(AliPID::kPion, AliPID::kPion);
533 mPt = esdV0->Pt();
70da6c5a 534 } else {
535 // AOD Analysis - not possible to cut
bf892a6a 536 AliAODv0 *aodV0 = static_cast<AliAODv0 *>(v0);
faee3b18 537 aodV0->GetID();
70da6c5a 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
70da6c5a 546 fQA->Fill("h_Pt_K0", mPt);
faee3b18 547 fQA->Fill("h_InvMassK0s", invMass);
70da6c5a 548
70da6c5a 549 AliDebug(1, Form("K0 identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
550
faee3b18 551 // AFTER all K0 cuts
70da6c5a 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()));
faee3b18 555 fPionsK0->Add(new AliHFEV0info(daughter[0], daughter[1]->GetID(), v0id));
3a72645a 556 fIndices->Add(daughter[0]->GetID(), AliHFEV0cuts::kRecoPionK0);
70da6c5a 557 }
558 if(!fIndices->Find(daughter[1]->GetID())){
559 AliDebug(1, Form("Adding K0 Pion track with ID %d", daughter[1]->GetID()));
faee3b18 560 fPionsK0->Add(new AliHFEV0info(daughter[1], daughter[0]->GetID(), v0id));
3a72645a 561 fIndices->Add(daughter[1]->GetID(), AliHFEV0cuts::kRecoPionK0);
70da6c5a 562 }
faee3b18 563 fK0s->Add(v0);
70da6c5a 564 return kTRUE;
565}
566
567//____________________________________________________________
8c1c76e9 568Bool_t AliHFEV0pid::IsPhi(const TObject *v0) const {
70da6c5a 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];
70da6c5a 575 //Double_t invMass = 0.;
bf892a6a 576
577 if(!v0) return kFALSE;
70da6c5a 578
8c1c76e9 579 //Int_t pIndex = 0, nIndex = 0;
70da6c5a 580 if(IsESDanalysis()){
581 // ESD - cut V0
8c1c76e9 582 //AliESDv0 *esdV0 = static_cast<AliESDv0 *>(v0);
583 //pIndex = esdV0->GetPindex();
584 //nIndex = esdV0->GetNindex();
70da6c5a 585 } else {
586 // PRELIMINARY - !!!
587 // AOD Analysis - not possible to cut
588 }
589
590 return kTRUE;
591}
592
593//____________________________________________________________
594Bool_t AliHFEV0pid::IsLambda(TObject *v0){
595 //
596 // Identify Lambda
597 //
bf892a6a 598
599 if(!v0) return kFALSE;
600
70da6c5a 601 AliVTrack* daughter[2];
602 Int_t pIndex = 0, nIndex = 0;
70da6c5a 603 Double_t invMass = 0.;
faee3b18 604 Bool_t isLambda = kTRUE; // Lambda - kTRUE, Anti Lambda - kFALSE
3a72645a 605 Double_t mPt = 0.;
faee3b18 606 Int_t v0id = -1;
70da6c5a 607 if(IsESDanalysis()){
608 // ESD - cut V0
bf892a6a 609 AliESDv0 *esdV0 = static_cast<AliESDv0 *>(v0);
faee3b18 610 v0id = esdV0->GetLabel();
611 if(!fV0cuts->LambdaCuts(esdV0,isLambda)) return kFALSE;
3a72645a 612 mPt = esdV0->Pt();
faee3b18 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 }
70da6c5a 621 } else {
622 // PRELIMINARY - !!!
623 // AOD Analysis - not possible to cut
624
625 // again - two cases as above
bf892a6a 626 AliAODv0 *aodV0 = static_cast<AliAODv0 *>(v0);
faee3b18 627 v0id = aodV0->GetID();
70da6c5a 628 pIndex = aodV0->GetPosID();
629 nIndex = aodV0->GetNegID();
faee3b18 630 invMass = aodV0->MassLambda();
631 }
632
70da6c5a 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;
faee3b18 636
637 // lambda
638 if(isLambda){
3a72645a 639 fQA->Fill("h_Pt_L", mPt);
640 fQA->Fill("h_InvMassL", invMass);
641
faee3b18 642 if(!fIndices->Find(daughter[0]->GetID())){
643 fProtons->Add(new AliHFEV0info(daughter[0], daughter[1]->GetID(), v0id));
3a72645a 644 fIndices->Add(daughter[0]->GetID(), AliHFEV0cuts::kRecoProton);
70da6c5a 645 }
faee3b18 646 if(!fIndices->Find(daughter[1]->GetID())){
647 fPionsL->Add(new AliHFEV0info(daughter[1], daughter[0]->GetID(), v0id));
3a72645a 648 fIndices->Add(daughter[1]->GetID(), AliHFEV0cuts::kRecoPionL);
70da6c5a 649 }
70da6c5a 650 }
faee3b18 651 // antilambda
652 else{
3a72645a 653 fQA->Fill("h_Pt_AL", mPt);
654 fQA->Fill("h_InvMassAL", invMass);
faee3b18 655 if(!fIndices->Find(daughter [1]->GetID())){
656 fProtons->Add(new AliHFEV0info(daughter[1], daughter[0]->GetID(), v0id));
3a72645a 657 fIndices->Add(daughter[1]->GetID(), AliHFEV0cuts::kRecoProton);
faee3b18 658 }
659 if(!fIndices->Find(daughter [0]->GetID())){
660 fPionsL->Add(new AliHFEV0info(daughter[0], daughter[1]->GetID(), v0id));
3a72645a 661 fIndices->Add(daughter [0]->GetID(), AliHFEV0cuts::kRecoPionL);
70da6c5a 662 }
70da6c5a 663 }
faee3b18 664 if(isLambda) fLambdas->Add(v0);
665 else fAntiLambdas->Add(v0);
70da6c5a 666
70da6c5a 667 return kTRUE;
668}
669
3a72645a 670//____________________________________________________________
671Int_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//____________________________________________________________
e97c2edf 758void 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//____________________________________________________________
3a72645a 776void 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}
70da6c5a 806//____________________________________________________________
807AliHFEV0pid::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//____________________________________________________________
825AliHFEV0pid::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//____________________________________________________________
63bffcf1 836void AliHFEV0pid::AliHFEV0pidTrackIndex::Flush(){
70da6c5a 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//____________________________________________________________
860void 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//____________________________________________________________
871void AliHFEV0pid::AliHFEV0pidTrackIndex::Add(Int_t index, Int_t species){
872 //
873 // Add new index to the list of identified particles
874 //
875 switch(species){
3a72645a 876 case AliHFEV0cuts::kRecoElectron:
70da6c5a 877 fIndexElectron[fNElectrons++] = index;
878 break;
3a72645a 879 case AliHFEV0cuts::kRecoPionK0:
70da6c5a 880 fIndexPionK0[fNPionsK0++] = index;
881 break;
3a72645a 882 case AliHFEV0cuts::kRecoPionL:
70da6c5a 883 fIndexPionL[fNPionsL++] = index;
884 break;
3a72645a 885 case AliHFEV0cuts::kRecoProton:
70da6c5a 886 fIndexProton[fNProtons++] = index;
887 break;
888 };
889}
890
891//____________________________________________________________
892Bool_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){
3a72645a 899 case AliHFEV0cuts::kRecoElectron:
70da6c5a 900 container = fIndexElectron;
901 n = fNElectrons;
902 break;
3a72645a 903 case AliHFEV0cuts::kRecoPionK0:
70da6c5a 904 container = fIndexPionK0;
905 n = fNPionsK0;
906 break;
3a72645a 907 case AliHFEV0cuts::kRecoPionL:
70da6c5a 908 container = fIndexPionL;
909 n = fNPionsL;
910 break;
3a72645a 911 case AliHFEV0cuts::kRecoProton:
70da6c5a 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//____________________________________________________________
929Bool_t AliHFEV0pid::AliHFEV0pidTrackIndex::Find(Int_t index) const {
930 //
931 // Find index in all samples
932 //
3a72645a 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);
70da6c5a 937}
938
939//____________________________________________________________
faee3b18 940TList *AliHFEV0pid::GetListOfQAhistograms(){
70da6c5a 941 //
faee3b18 942 // Getter for V0 PID QA histograms
70da6c5a 943 //
3a72645a 944
e156c3bb 945 CLRBIT(fDestBits, 1);
946
faee3b18 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 }
3a72645a 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
faee3b18 962 return fOutput;
70da6c5a 963}