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