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