]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEV0pid.cxx
Fixes to cure warnings
[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//
27#include <TDatabasePDG.h>
28#include <TObjArray.h>
29#include <TPDGCode.h>
30#include <TString.h>
31
32#include <TDatabasePDG.h>
33
34#include "AliAODEvent.h"
35#include "AliAODTrack.h"
36#include "AliAODv0.h"
37#include "AliAODVertex.h"
38#include "AliESDEvent.h"
39#include "AliESDtrack.h"
40#include "AliESDv0.h"
41#include "AliESDVertex.h"
42#include "AliKFParticle.h"
43#include "AliKFVertex.h"
44#include "AliVEvent.h"
45#include "AliVTrack.h"
46
47#include "AliHFEV0pid.h"
48ClassImp(AliHFEV0pid)
49
50//____________________________________________________________
51AliHFEV0pid::AliHFEV0pid():
52 TObject()
53 , fInputEvent(NULL)
54 , fPrimaryVertex(NULL)
55 , fElectrons(NULL)
56 , fPionsK0(NULL)
57 , fPionsL(NULL)
58 , fKaons(NULL)
59 , fProtons(NULL)
60 , fIndices(NULL)
61 , fQA(NULL)
62{
63 //
64 // Default constructor
65 //
66 fElectrons = new TObjArray();
67 fPionsK0 = new TObjArray();
68 fPionsL = new TObjArray();
69 fKaons = new TObjArray();
70 fProtons = new TObjArray();
71 fIndices = new AliHFEV0pidTrackIndex;
72}
73
74//____________________________________________________________
75AliHFEV0pid::~AliHFEV0pid(){
76 //
77 // Destructor
78 // Remove Containers
79 //
80 if(fInputEvent) delete fInputEvent;
81 //if(fPrimaryVertex) delete fPrimaryVertex;
82 if(fElectrons) delete fElectrons;
83 if(fPionsK0) delete fPionsK0;
84 if(fPionsL) delete fPionsL;
85 if(fKaons) delete fKaons;
86 if(fProtons) delete fProtons;
87 if(fIndices) delete fIndices;
88 if(fQA) delete fQA;
89}
90
91//____________________________________________________________
92void AliHFEV0pid::InitQA(){
93 //
94 // Initialize QA histograms
95 //
96 if(!fQA){
97 fQA = new AliHFEcollection("v0pidQA", "QA histograms for V0 PID");
98
99 // QA histograms for cut statistics
100 fQA->CreateTH1F("h_cutEfficiencyGamma", "Cut Efficiency for Gammas", 10, 0, 10);
101 fQA->CreateTH1F("h_cutEfficiencyK0s", "Cut Efficiency for K0s", 10, 0, 10);
102 fQA->CreateTH1F("h_cutEfficiencyPhi", "Cut Efficiency for Phi", 10, 0, 10);
103 fQA->CreateTH1F("h_cutEfficiencyLambda", "Cut Efficiency for Lambdas", 10, 0, 10);
104
105 // QA histograms for invariant mass
91c7e1ec 106 fQA->CreateTH1F("h_InvMassGamma", "Gamma invariant mass; inv mass [GeV/c^{2}]; counts", 100, 0, 0.25);
70da6c5a 107 fQA->CreateTH1F("h_InvMassK0s", "K0s invariant mass; inv mass [GeV/c^{2}]; counts", 100, 0.4, 0.65);
108 fQA->CreateTH1F("h_InvMassPhi", "Phi invariant mass; inv mass [GeV/c^{2}]; counts", 100, 0.4, 0.65);
109 fQA->CreateTH1F("h_InvMassLambda", "Lambda invariant mass; inv mass [GeV/c^{2}]; counts", 100, 1.05, 1.15);
110
91c7e1ec 111 // QA histograms for p distribution (of the daughters)
112 fQA->CreateTH1F("h_P_electron", "P distribution of the gamma electrons; p (GeV/c); counts", 100, 0.1, 10);
113 fQA->CreateTH1F("h_P_K0pion", "P distribution of the K0 pions; p (GeV/c); counts", 100, 0.1, 10);
114 fQA->CreateTH1F("h_P_Lpion", "P distribution of the Lambda pions; p (GeV/c); counts", 100, 0.1, 10);
115 fQA->CreateTH1F("h_P_Lproton", "P distribution of the Lambda protons; p (GeV/c); counts", 100, 0.1, 10);
116
70da6c5a 117 // QA invariant mass as a functin of pt
118 fQA->CreateTH1Fvector1(20, "h_InvMassGamma_pt", "Gamma invarinat mass in pt bins; inv mass [GeV/c^{2}]; counts", 250, 0, 2);
119 fQA->CreateTH1Fvector1(20, "h_InvMassK0_pt", "K0 invarinat mass in pt bins; inv mass [GeV/c^{2}]; counts", 250, 0, 2);
120 fQA->CreateTH1Fvector1(20, "h_InvMassPhi_pt", "Phi invarinat mass in pt bins; inv mass [GeV/c^{2}]; counts", 250, 0, 2);
121 fQA->CreateTH1Fvector1(20, "h_InvMassLambda_pt", "Lambda invarinat mass in pt bins; inv mass [GeV/c^{2}]; counts", 250, 0, 2);
122
123 // QA pt of the V0
124 fQA->CreateTH1F("h_Pt_Gamma", "Pt of the gamma conversion; p_{T} (GeV/c); counts", 100, 0, 10);
125 fQA->CreateTH1F("h_Pt_K0", "Pt of the K0; p_{T} (GeV/c); counts", 100, 0, 10);
126 fQA->CreateTH1F("h_Pt_Phi", "Pt of the Phi; p_{T} (GeV/c); counts", 100, 0, 10);
127 fQA->CreateTH1F("h_Pt_Lambda", "Pt of the Lambda; p_{T} (GeV/c); counts", 100, 0, 10);
128 //fQA->CreateTH1F("h_Pt_electrons", "Pt of the conversion electrons; p_{T} (GeV/c); counts");
129 //fQA->CreateTH1F("h_Pt_pionsK0", "Pt of the K0 pions; p_{T} (GeV/c); counts");
130 //fQA->CreateTH1F("h_Pt_pionsL", "Pt of the Lambda pions; p_{T} (GeV/c); counts");
131 //fQA->CreateTH1F("h_Pt_protons", "Pt of the Lambda protons; p_{T} (GeV/c); counts");
132
133
134 // QA histogram for both Lambda candidate combinations -
135 fQA->CreateTH2F("h_L0_dca_v_dMass", "L0 dca verus dMass; dMass [GeV/c^{2}]; dDCA [cm]; ", 100, -1., 1., 100, 0., 5.);
136
137 // Chi2 histograms
138 fQA->CreateTH1F("h_chi2_gamma", "Chi2 for gammas", 10000, 0, 1000);
139 fQA->CreateTH1F("h_chi2_K0s", "Chi2 for K0s", 10000, 0, 500);
140 fQA->CreateTH1F("h_chi2_Phi", "Chi2 for K0s", 10000, 0, 500);
141 fQA->CreateTH1F("h_chi2_Lambda", "Chi2 for Lambdas", 10000, 0, 1000);
142 }
143}
144
145//____________________________________________________________
146void AliHFEV0pid::Process(AliVEvent * const inputEvent){
147
148 //
149 // Find protons, pions and electrons using V0 decays and
150 // store the pointers in the TObjArray
151 //
152 Int_t nGamma = 0, nK0s = 0, nLambda = 0, nPhi = 0;
153 fInputEvent = inputEvent;
154
155 fIndices->Init(fInputEvent->GetNumberOfV0s() * 2);
156 fPrimaryVertex = new AliKFVertex(*(fInputEvent->GetPrimaryVertex()));
157 Int_t v0status = 0;
158 for(Int_t iv0 = 0; iv0 < fInputEvent->GetNumberOfV0s(); iv0++){
159 if(!TString(fInputEvent->IsA()->GetName()).CompareTo("AliESDEvent")){
160 // case ESD
161 SetESDanalysis();
162 AliESDv0 *esdV0 = (dynamic_cast<AliESDEvent *>(fInputEvent))->GetV0(iv0);
163 if(!esdV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
164 v0status = ProcessV0(esdV0);
165 } else {
166 // case AOD
167 SetAODanalysis();
168 AliAODv0 *aodV0 = (dynamic_cast<AliAODEvent *>(fInputEvent))->GetV0(iv0);
169 if(aodV0->GetOnFlyStatus()) continue; // Take only V0s from the On-the-fly v0 finder
170 v0status = ProcessV0(aodV0);
171 }
172 switch(v0status){
173 case kRecoGamma: nGamma++; break;
174 case kRecoK0s: nK0s++; break;
175 case kRecoPhi: nPhi++; break;
176 case kRecoLambda: nLambda++; break;
177 };
178 }
179
180 AliDebug(1, Form("Number of gammas : %d", nGamma));
181 AliDebug(1, Form("Number of K0s : %d", nK0s));
182 AliDebug(1, Form("Number of Phis : %d", nPhi));
183 AliDebug(1, Form("Number of Lambdas : %d", nLambda));
184
185 AliDebug(1, "Number of stored tracks:");
186 AliDebug(1, Form("Number of electrons : %d", fElectrons->GetEntries()));
187 AliDebug(1, Form("Number of K0 pions : %d", fPionsK0->GetEntries()));
188 AliDebug(1, Form("Number of Lambda pions : %d", fPionsL->GetEntries()));
189 AliDebug(1, Form("Number of Phi kaons : %d", fKaons->GetEntries()));
190 AliDebug(1, Form("Number of protons : %d", fProtons->GetEntries()));
191
192 delete fPrimaryVertex;
193}
194
195//____________________________________________________________
196Int_t AliHFEV0pid::ProcessV0(TObject *v0){
197 //
198 // Process single V0
199 // Apply general cut and special cuts for gamma, K0s, Lambda
200 //
201 AliVTrack* daughter[2];
202 daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack((dynamic_cast<AliESDv0 *>(v0))->GetPindex()));
203 daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack((dynamic_cast<AliESDv0 *>(v0))->GetPindex()));
204 if(!daughter[0] || !daughter[1]) return kUndef;
205
206 if(IsESDanalysis()){
207 for(Int_t i=0; i<2; ++i){
208 if(!CutESDtrack(dynamic_cast<AliESDtrack*>(daughter[i]))) return kUndef;
209 }
210 }
211
212 if(IsGammaConv(v0)) return kRecoGamma;
213 else if(IsK0s(v0)) return kRecoK0s;
214 else if(IsLambda(v0)) return kRecoLambda;
215 else return kUndef;
216
217
218}
219//____________________________________________________________
220void AliHFEV0pid::Flush(){
221 //
222 // Clear the Lists
223 //
224 AliDebug(1, "Flushing containers");
225 fProtons->Clear();
226 fPionsK0->Clear();
227 fPionsL->Clear();
228 fElectrons->Clear();
229 fIndices->Flush();
230}
231
232//____________________________________________________________
233Bool_t AliHFEV0pid::IsGammaConv(TObject *v0){
234 //
235 // Identify Gamma
236 //
237 AliVTrack* daughter[2];
238 Int_t pIndex = 0, nIndex = 0;
239 Double_t invMass = 0.;
240 if(IsESDanalysis()){
241 // ESD - cut V0
242 AliESDv0 *esdV0 = dynamic_cast<AliESDv0 *>(v0);
243 if(!CutV0(esdV0, kRecoGamma)) return kFALSE;
244 if(LooseRejectK0(esdV0) || LooseRejectLambda(esdV0)) return kFALSE;
245 // DEBUG
246 //invMass = esdV0->GetEffMass(AliPID::kElectron, AliPID::kElectron);
247 invMass = GetEffMass(esdV0, AliPID::kElectron, AliPID::kElectron);
248 //..
249
250 pIndex = esdV0->GetPindex();
251 nIndex = esdV0->GetNindex();
252 } else {
253 // AOD Analysis - not possible to cut
254 AliAODv0 *aodV0 = dynamic_cast<AliAODv0 *>(v0);
255 pIndex = aodV0->GetPosID();
256 nIndex = aodV0->GetNegID();
257 invMass = aodV0->InvMass2Prongs(0, 1, kElectron, kElectron);
258 }
259 daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
260 daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
261 if(!daughter[0] || !daughter[1]) return kFALSE;
262
263 // Get Invariant mass and chi2/ndf
264 AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kElectron), TMath::Abs(kElectron));
265 kfMother->SetProductionVertex(*fPrimaryVertex);
266 kfMother->SetMassConstraint(0, 1.);
267 Double_t ndf = kfMother->GetNDF();
268 Double_t chi2 = kfMother->GetChi2();
269 delete kfMother;
270
271 if(fQA) fQA->Fill("h_chi2_gamma", chi2/ndf);
272 if(chi2/ndf > 7) return kFALSE;
273 Double_t mPt = kfMother->GetPt();
274 fQA->Fill("h_Pt_Gamma", mPt);
275 Int_t ptBin = (int)(mPt*10.0);
70da6c5a 276
277 if(fQA) fQA->Fill("h_InvMassGamma", invMass);
278 if(invMass > 0.05) return kFALSE;
91c7e1ec 279 fQA->Fill("h_InvMassGamma_pt", ptBin+1, invMass);
280
70da6c5a 281 AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
282 // Identified gamma - store tracks in the electron containers
283 if(!fIndices->Find(daughter[0]->GetID())){
284 AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
285 fElectrons->Add(daughter[0]);
286 fIndices->Add(daughter[0]->GetID(), AliHFEV0pid::kRecoElectron);
287 }
288 if(!fIndices->Find(daughter[1]->GetID())){
289 AliDebug(1, Form("Gamma identified, daughter IDs: %d,%d", daughter[1]->GetID(), daughter[1]->GetID()));
290 fElectrons->Add(daughter[1]);
291 fIndices->Add(daughter[1]->GetID(), AliHFEV0pid::kRecoElectron);
292 }
293 return kTRUE;
294}
295
296
297
298//____________________________________________________________
299Bool_t AliHFEV0pid::IsK0s(TObject *v0){
300 //
301 // Identify K0s
302 //
303 const Double_t kK0smass=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass(); // PDG K0s mass
304 AliVTrack* daughter[2];
305 Int_t pIndex = 0, nIndex = 0;
306 Double_t invMass = 0.;
307 if(IsESDanalysis()){
308 // ESD - cut V0
309 AliESDv0 *esdV0 = dynamic_cast<AliESDv0 *>(v0);
310 if(!CutV0(esdV0, kRecoK0s)) return kFALSE;
311 invMass = esdV0->GetEffMass(AliPID::kPion, AliPID::kPion);
312 pIndex = esdV0->GetPindex();
313 nIndex = esdV0->GetNindex();
314 } else {
315 // AOD Analysis - not possible to cut
316 AliAODv0 *aodV0 = dynamic_cast<AliAODv0 *>(v0);
317 pIndex = aodV0->GetPosID();
318 nIndex = aodV0->GetNegID();
319 invMass = aodV0->MassK0Short();
320 }
321 daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
322 daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
323 if(!daughter[0] || !daughter[1]) return kFALSE;
324
325 // Get Invariant mass and chi2/ndf
326 AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kPiPlus));
327 kfMother->SetProductionVertex(*fPrimaryVertex);
328 kfMother->SetMassConstraint(kK0smass, 0.);
329 Double_t ndf = kfMother->GetNDF();
330 Double_t chi2 = kfMother->GetChi2();
331 delete kfMother;
332
333 if(fQA) fQA->Fill("h_chi2_K0s", chi2/ndf);
334 if(chi2/ndf > 5) return kFALSE;
335 Double_t mPt = kfMother->GetPt();
336 fQA->Fill("h_Pt_K0", mPt);
337 Int_t ptBin = (int)(mPt*10.0);
338 fQA->Fill("h_InvMassK0_pt", ptBin+1, invMass);
339
340 if(fQA) fQA->Fill("h_InvMassK0s", invMass);
341
342 if(invMass < 0.485 || invMass > 0.51) return kFALSE;
343 AliDebug(1, Form("K0 identified, daughter IDs: %d,%d", daughter[0]->GetID(), daughter[1]->GetID()));
344
345 // Identified gamma - store tracks in the electron containers
346 if(!fIndices->Find(daughter[0]->GetID())){
347 AliDebug(1, Form("Adding K0 Pion track with ID %d", daughter[0]->GetID()));
348 fPionsK0->Add(daughter[0]);
349 fIndices->Add(daughter[0]->GetID(), AliHFEV0pid::kRecoPionK0);
350 }
351 if(!fIndices->Find(daughter[1]->GetID())){
352 AliDebug(1, Form("Adding K0 Pion track with ID %d", daughter[1]->GetID()));
353 fPionsK0->Add(daughter[1]);
354 fIndices->Add(daughter[1]->GetID(), AliHFEV0pid::kRecoPionK0);
355 }
356 return kTRUE;
357}
358
359//____________________________________________________________
360Bool_t AliHFEV0pid::IsPhi(TObject *v0){
361 //
362 // Identify Phi - very preliminary - requires diffrent approach (V0 fnder is not effective)
363 //
364
365 //const Double_t kPhiMass=TDatabasePDG::Instance()->GetParticle(333)->Mass(); // PDG phi mass
366 //AliVTrack* daughter[2];
367 //AliKFParticle *mother = NULL;
368 //Double_t invMass = 0.;
369
370 Int_t pIndex = 0, nIndex = 0;
371 if(IsESDanalysis()){
372 // ESD - cut V0
373 AliESDv0 *esdV0 = dynamic_cast<AliESDv0 *>(v0);
374 if(!CutV0(esdV0, kRecoPhi)) return kFALSE;
375 if(LooseRejectGamma(esdV0) || LooseRejectK0(esdV0)) return kFALSE;
376 pIndex = esdV0->GetPindex();
377 nIndex = esdV0->GetNindex();
378 } else {
379 // PRELIMINARY - !!!
380 // AOD Analysis - not possible to cut
381 }
382
383 return kTRUE;
384}
385
386//____________________________________________________________
387Bool_t AliHFEV0pid::IsLambda(TObject *v0){
388 //
389 // Identify Lambda
390 //
391 const Double_t kL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(); // PDG lambda mass
392 AliVTrack* daughter[2];
393 Int_t pIndex = 0, nIndex = 0;
394 AliKFParticle *mother[2] = {NULL, NULL};
395 Float_t dMass[2]; // lambda mass difference for the two hypoteses
396 //Int_t lambda = 0; // [1] for lambda and [-1] for anti-lambda
397 Double_t chi2 = 0.;
398
399 Double_t invMass = 0.;
400 Double_t invMassEOD = 0;
401 if(IsESDanalysis()){
402 // ESD - cut V0
403 AliESDv0 *esdV0 = dynamic_cast<AliESDv0 *>(v0);
404 if(!CutV0(esdV0, kRecoLambda)) return kFALSE;
405 if(LooseRejectK0(esdV0) || LooseRejectGamma(esdV0)) return kFALSE;
406 pIndex = esdV0->GetPindex();
407 nIndex = esdV0->GetNindex();
408
409 //invMass = esdV0->GetEffMass(AliPID::kPion, AliPID::kPion);
410
411 } else {
412 // PRELIMINARY - !!!
413 // AOD Analysis - not possible to cut
414
415 // again - two cases as above
416 AliAODv0 *aodV0 = dynamic_cast<AliAODv0 *>(v0);
417 pIndex = aodV0->GetPosID();
418 nIndex = aodV0->GetNegID();
419 invMassEOD = aodV0->MassLambda();
420 }
421
422 daughter[0] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(pIndex));
423 daughter[1] = dynamic_cast<AliVTrack *>(fInputEvent->GetTrack(nIndex));
424 if(!daughter[0] || !daughter[1]) return kFALSE;
425
426 //
427 // now - go over two case - Lambda and AntiLambda
428 // choose the on ewhere the resulting lambda mass is close to the
429 // expected value and at the same time points closer to the primary vertex in XY plane
430 //
431 // A) lambda -> proton + negative-pion
432 // B) anti-lambda -> anti-proton + positive-pion
433 //
434
435 //
436 // A) :: proton
437 //
438 mother[0] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kProton), TMath::Abs(kPiPlus));
439 AliHFELambdaInf lambda0(mother[0], fPrimaryVertex);
440
441 // Undo changes
442 *fPrimaryVertex -= *mother[0];
443 AddTrackToKFVertex(daughter[0], TMath::Abs(kProton));
444 AddTrackToKFVertex(daughter[1], TMath::Abs(kPiPlus));
445
446 //
447 // B) :: anti-proton
448 //
449 mother[1] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kProton));
450 AliHFELambdaInf lambda1(mother[1], fPrimaryVertex);
451
452 // qa histograms
453 dMass[0] = lambda0.GetInvariantMass() - kL0mass;
454 dMass[1] = lambda1.GetInvariantMass() - kL0mass;
455 fQA->Fill("h_L0_dca_v_dMass", dMass[0], lambda0.GetDistanceFromPrimaryVertex());
456 fQA->Fill("h_L0_dca_v_dMass", dMass[1], lambda1.GetDistanceFromPrimaryVertex());
457
458 //
459 // decide on the true Lambda candidate bsaed on
460 // - mass difference
461 // - vertex dca (testing needed first)
462 //
463 AliHFELambdaInf *lambdaInf[2] = {&lambda0, &lambda1};
464 Bool_t isLambda = kFALSE;
465 Int_t index = -1;
466 index = (dMass[0] <dMass[1]) ? 0 : 1;
467 invMass = lambdaInf[index]->GetInvariantMass();
468 chi2 = lambdaInf[index]->GetChi2NDF();
469 //if(!IsESDanalysis()){
470 // invMass = invMassEOD;
471 //}
472 //else{
473 // invMass = mother[index]->GetMass();
474 //}
475 if(fQA) fQA->Fill("h_chi2_Lambda", chi2);
476 if(chi2 < 3){
477 if(fQA) fQA->Fill("h_InvMassLambda", invMass);
478 Double_t mPt = mother[index]->GetPt();
479 fQA->Fill("h_Pt_Lambda", mPt);
480 Int_t ptBin = (int)(mPt*10.0);
481 fQA->Fill("h_InvMassLambda_pt", ptBin+1, invMass);
482
483 // cut on the invariant mass for the proton and pion selection
484 if(invMass > 1.11 || invMass < 1.12){
485 // Identified lambdas - store the protons and pions and update primary vertex
486 *fPrimaryVertex += *mother[index];
487
488 // lambda
489 if(0 == index){
490 if(!fIndices->Find(daughter[0]->GetID())){
491 fProtons->Add(daughter[0]);
492 fIndices->Add(daughter[0]->GetID(), AliHFEV0pid::kRecoProton);
493 }
494 if(!fIndices->Find(daughter[1]->GetID())){
495 fPionsL->Add(daughter[1]);
496 fIndices->Add(daughter[1]->GetID(), AliHFEV0pid::kRecoPionL);
497 }
498 }
499 // antilambda
500 if(1 == index){
501 if(!fIndices->Find(daughter [1]->GetID())){
502 fProtons->Add(daughter[1]);
503 fIndices->Add(daughter[1]->GetID(), AliHFEV0pid::kRecoProton);
504 }
505 if(!fIndices->Find(daughter [0]->GetID())){
506 fPionsL->Add(daughter[0]);
507 fIndices->Add(daughter [0]->GetID(), AliHFEV0pid::kRecoPionL);
508 }
509 }
510 isLambda = kTRUE;
511 }
512 }
513
514 if(!isLambda){
515 // Add the daughters again to the primary vertex
516 AddTrackToKFVertex(daughter[0], TMath::Abs(kPiPlus));
517 AddTrackToKFVertex(daughter[1], TMath::Abs(kProton));
518 }
519 // remove the objects
520 for(Int_t i=0; i<2; ++i){
521 if (mother[i]) delete mother[i];
522 }
523
524 return isLambda;
525}
526
527//____________________________________________________________
528AliKFParticle *AliHFEV0pid::CreateMotherParticle(AliVTrack *pdaughter, AliVTrack *ndaughter, Int_t pspec, Int_t nspec){
529 //
530 // Creates a mother particle
531 //
532 AliKFParticle pkfdaughter(*pdaughter, pspec);
533 AliKFParticle nkfdaughter(*ndaughter, nspec);
534
535 // check if the daughter particles are coming from the primary vertex
536 if(IsESDanalysis()){
537 // ESD Analyis
538 const AliESDVertex *esdvertex = dynamic_cast<const AliESDVertex *>(fInputEvent->GetPrimaryVertex());
539 UShort_t *contrib = esdvertex->GetIndices();
540
541 Int_t nfound = 0;
542 for(Int_t id = 0; id < esdvertex->GetNIndices(); id++){
543 if(contrib[id] == pdaughter->GetID()){
544 *fPrimaryVertex -= pkfdaughter;
545 nfound++;
546 }
547 if(contrib[id] == ndaughter->GetID()){
548 *fPrimaryVertex -= nkfdaughter;
549 nfound++;
550 }
551 if(nfound == 2) break;
552 }
553 } else {
554 // AOD Analysis: AOD Vertex
555 const AliAODVertex *aodvertex = dynamic_cast<const AliAODVertex *>(fInputEvent->GetPrimaryVertex());
556 if(aodvertex->HasDaughter(pdaughter))
557 *fPrimaryVertex -= pkfdaughter;
558 if(aodvertex->HasDaughter(ndaughter))
559 *fPrimaryVertex -= nkfdaughter;
560 }
561
562 // Create the mother particle and add them to the primary vertex
563 AliKFParticle *mother = new AliKFParticle(pkfdaughter, nkfdaughter);
564 *fPrimaryVertex += *mother;
565
566 return mother;
567}
568
569//____________________________________________________________
570void AliHFEV0pid::AddTrackToKFVertex(AliVTrack *track, Int_t species){
571 //
572 // Add track to the primary vertex (if it was used in the vertex
573 // calculation)
574 //
575 Bool_t isAdd = kFALSE;
576 if(IsESDanalysis()){
577 const AliESDVertex *esdvertex = dynamic_cast<const AliESDVertex *>(fInputEvent->GetPrimaryVertex());
578 UShort_t *contrib = esdvertex->GetIndices();
579 for(Int_t id = 0; id < esdvertex->GetNIndices(); id++){
580 if(contrib[id] == track->GetID()){
581 isAdd = kTRUE;
582 break;
583 }
584 }
585 } else {
586 const AliAODVertex *aodvertex = dynamic_cast<const AliAODVertex *>(fInputEvent->GetPrimaryVertex());
587 if(aodvertex->HasDaughter(track)) isAdd = kTRUE;
588 }
589 if(isAdd){
590 AliKFParticle kftrack(*track, species);
591 *fPrimaryVertex += kftrack;
592 }
593}
594
595//____________________________________________________________
596Bool_t AliHFEV0pid::CutV0(AliESDv0 *v0, Int_t V0species){
597 //
598 // Cut the V0
599 //
600 Int_t cutRequired = 0;
601 // For cut always take min and max
602 Double_t cutCosPoint[2] = {0., 0.},
603 cutDCA[2] = {0., 0.},
604 cutProdVtx[2] = {0., 0.},
605 cutOpAng[2] = {0., 0.},
606 cutPsiPair[2] = {0., 0.};
607
608 switch(V0species){
609 case kRecoGamma:
610 cutCosPoint[1] = 0.03;
611 SETBIT(cutRequired, 1);
612 cutDCA[1] = 0.25;
613 SETBIT(cutRequired, 3);
614 cutProdVtx[0] = 6;
615 SETBIT(cutRequired, 4);
616 cutOpAng[1] = 0.1;
617 SETBIT(cutRequired, 7);
618 cutPsiPair[1] = 0.05;
619 SETBIT(cutRequired, 9);
620 break;
621 case kRecoK0s:
622 cutCosPoint[1] = 0.03;
623 SETBIT(cutRequired, 1);
624 cutDCA[1] = 0.1;
625 SETBIT(cutRequired, 3);
626 cutProdVtx[1] = 8.1;
627 SETBIT(cutRequired, 5);
628 break;
629 case kRecoPhi:
630 break;
631 case kRecoLambda:
632 cutCosPoint[1] = 0.03;
633 SETBIT(cutRequired, 1);
634 cutDCA[1] = 0.2;
635 SETBIT(cutRequired, 3);
636 cutProdVtx[1] = 24;
637 SETBIT(cutRequired, 5);
638 break;
639 default:
640 // unidentified, return
641 return kFALSE;
642 };
643
644 Char_t hname[256];
645 const Char_t *specname[4] = {"Gamma", "K0s", ", Phi", "Lambda"};
646 sprintf(hname, "h_cutEfficiency%s", specname[V0species-1]);
647
648 // Cut on pointing angle
649 Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
650 if(TESTBIT(cutRequired, 0) && TMath::ACos(cosPoint) < cutCosPoint[0]) return kFALSE;
651 if(fQA) fQA->Fill(hname, 0);
652 if(TESTBIT(cutRequired, 1) && TMath::ACos(cosPoint) > cutCosPoint[1]) return kFALSE;
653 if(fQA) fQA->Fill(hname, 1);
654
655 // Cut on DCA between daughters
656 Double_t dca = v0->GetDcaV0Daughters();
657 if(TESTBIT(cutRequired, 2) && dca < cutDCA[0]) return kFALSE;
658 if(fQA) fQA->Fill(hname, 2);
659 if(TESTBIT(cutRequired, 3) && dca > cutDCA[1]) return kFALSE;
660 if(fQA) fQA->Fill(hname, 3);
661
662 // Cut on reconstructed verted position
663 Double_t x, y, z;
664 v0->GetXYZ(x,y,z);
665 Double_t r = TMath::Sqrt(x*x + y*y);
666 if(TESTBIT(cutRequired, 4) && r < cutProdVtx[0]) return kFALSE;
667 if(fQA) fQA->Fill(hname, 4);
668 if(TESTBIT(cutRequired, 5) && r > cutProdVtx[1]) return kFALSE;
669 if(fQA) fQA->Fill(hname, 5);
670
671 //Cut on Opening angle (conversions only)
672 if(TESTBIT(cutRequired, 6) && OpenAngle(v0) < cutOpAng[0]) return kFALSE;
673 if(fQA) fQA->Fill(hname, 6);
674 if(TESTBIT(cutRequired, 7) && OpenAngle(v0) > cutOpAng[1]) return kFALSE;
675 if(fQA) fQA->Fill(hname, 7);
676
677 //Cut on PsiPair angle (conversons only)
678 if(TESTBIT(cutRequired, 8) && PsiPair(v0) < cutPsiPair[0]) return kFALSE;
679 if(fQA) fQA->Fill(hname, 8);
680 if(TESTBIT(cutRequired, 9) && PsiPair(v0) > cutPsiPair[1]) return kFALSE;
681 if(fQA) fQA->Fill(hname, 9);
682 return kTRUE;
683}
684
685//____________________________________________________________
686Bool_t AliHFEV0pid::CutESDtrack(AliESDtrack *track){
687 //
688 // Single track cuts
689 //
690 // Hard coaded cut values for the beginning
691 //
692
693 if(!track) return kFALSE;
694
695 // status word
696 ULong_t status = track->GetStatus();
697
698 // DCA - to Vertex R & Z
699 Float_t dcaR = -1.;
700 Float_t dcaZ = -1.;
701 track->GetImpactParameters(dcaR, dcaZ);
702 //if(dcaR > 4.0) return kFALSE;
703 //if(dcaZ > 10.0) return kFALSE;
704
705 // No. of TPC clusters
706 if(track->GetTPCNcls() < 80) return kFALSE;
707
708 // TPC refit
709 if(!(status & AliESDtrack::kTPCrefit)) return kFALSE;
710
711 // Chi2 per TPC cluster
712 Int_t nTPCclusters = track->GetTPCclusters(0);
713 Float_t chi2perTPCcluster = track->GetTPCchi2()/Float_t(nTPCclusters);
714 if(chi2perTPCcluster > 3.5) return kFALSE;
715
716 // TPC cluster ratio
717 Float_t cRatioTPC = track->GetTPCNclsF() > 0. ? static_cast<Float_t>(track->GetTPCNcls())/static_cast<Float_t> (track->GetTPCNclsF()) : 1.;
718 if(cRatioTPC < 0.6) return kFALSE;
719
720 // kinks
721 if(track->GetKinkIndex(0) != 0) return kFALSE;
722
723 // Covariance matrix - TO BE RECONSIDERED
724 Double_t extCov[15];
725 track->GetExternalCovariance(extCov);
726 //if(extCov[0] > 2. ) return kFALSE;
727 //if(extCov[2] > 2. ) return kFALSE;
728 //if(extCov[5] > 0.5) return kFALSE;
729 //if(extCov[9] > 0.5) return kFALSE;
730 //if(extCov[14] > 2. ) return kFALSE;
731
732 // pt
733 if(track->Pt() < 0.1 || track->Pt() > 100) return kFALSE;
734
735 // eta
736 if(TMath::Abs(track->Eta()) > 0.9) return kFALSE;
737
738 // the track made it through! :-)
739 return kTRUE;
740}
741
742//_________________________________________________
743Bool_t AliHFEV0pid::LooseRejectK0(AliESDv0 * const v0) const {
744 //
745 // Reject K0 based on loose cuts
746 //
747 Double_t mass = v0->GetEffMass(AliPID::kPion, AliPID::kPion);
748 if(mass > 0.494 && mass < 0.501) return kTRUE;
749 return kFALSE;
750}
751
752//_________________________________________________
753Bool_t AliHFEV0pid::LooseRejectLambda(AliESDv0 * const v0) const {
754 //
755 // Reject Lambda based on loose cuts
756 //
757 Double_t mass1 = v0->GetEffMass(AliPID::kPion, AliPID::kProton);
758 Double_t mass2 = v0->GetEffMass(AliPID::kProton, AliPID::kPion);
759
760 if(mass1 > 1.1 && mass1 < 1.12) return kTRUE;
761 if(mass2 > 1.1 && mass2 < 1.12) return kTRUE;
762 return kFALSE;
763}
764
765//_________________________________________________
766Bool_t AliHFEV0pid::LooseRejectGamma(AliESDv0 * const v0) const {
767 //
768 // Reject Lambda based on loose cuts
769 //
770
771 //Double_t mass = v0->GetEffMass(AliPID::kElectron, AliPID::kElectron);
772 // DEBUG temporary solution, see the comment in GetEffMass
773 Double_t mass = GetEffMass(v0, AliPID::kElectron, AliPID::kElectron);
774 //..
775 if(mass < 0.02) return kTRUE;
776 return kFALSE;
777}
778
779//_________________________________________________
780Float_t AliHFEV0pid::OpenAngle(AliESDv0 *v0) const {
781 //
782 // Opening angle between two daughter tracks
783 //
784 Double_t mn[3] = {0,0,0};
785 Double_t mp[3] = {0,0,0};
786
787
788 v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
789 v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
790
791
792 Float_t openAngle = TMath::ACos((mp[0]*mn[0] + mp[1]*mn[1] + mp[2]*mn[2])/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1] + mp[2]*mp[2])*TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1] + mn[2]*mn[2])));
793
794 return TMath::Abs(openAngle);
795}
796
797//_________________________________________________
798Float_t AliHFEV0pid::PsiPair(AliESDv0 *v0) {
799 //
800 // Angle between daughter momentum plane and plane
801 //
802 Float_t magField = fInputEvent->GetMagneticField();
803
804 Int_t pIndex = v0->GetPindex();
805 Int_t nIndex = v0->GetNindex();
806
807 AliESDtrack* daughter[2];
808
809 daughter[0] = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(pIndex));
810 daughter[1] = dynamic_cast<AliESDtrack *>(fInputEvent->GetTrack(nIndex));
811
812 Double_t x, y, z;
813 v0->GetXYZ(x,y,z);//Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
814
815 Double_t mn[3] = {0,0,0};
816 Double_t mp[3] = {0,0,0};
817
818
819 v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
820 v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
821
822
823 Double_t deltat = 1.;
824 deltat = TMath::ATan(mp[2]/(TMath::Sqrt(mp[0]*mp[0] + mp[1]*mp[1])+1.e-13)) - TMath::ATan(mn[2]/(TMath::Sqrt(mn[0]*mn[0] + mn[1]*mn[1])+1.e-13));//difference of angles of the two daughter tracks with z-axis
825
826 Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
827
828 Double_t momPosProp[3];
829 Double_t momNegProp[3];
830
831 AliExternalTrackParam pt(*daughter[0]), nt(*daughter[1]);
832
833 Float_t psiPair = 4.;
834
835 if(nt.PropagateTo(radiussum,magField) == 0)//propagate tracks to the outside
836 psiPair = -5.;
837 if(pt.PropagateTo(radiussum,magField) == 0)
838 psiPair = -5.;
839 pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
840 nt.GetPxPyPz(momNegProp);
841
842 Double_t pEle =
843 TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
844 Double_t pPos =
845 TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
846
847 Double_t scalarproduct =
848 momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
849
850 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
851
852 psiPair = TMath::Abs(TMath::ASin(deltat/chipair));
853
854 return psiPair;
855}
856
857//____________________________________________________________
858AliHFEV0pid::AliHFEV0pidTrackIndex::AliHFEV0pidTrackIndex():
859 fNElectrons(0)
860 , fNPionsK0(0)
861 , fNPionsL(0)
862 , fNKaons(0)
863 , fNProtons(0)
864 , fIndexElectron(NULL)
865 , fIndexPionK0(NULL)
866 , fIndexPionL(NULL)
867 , fIndexKaon(NULL)
868 , fIndexProton(NULL)
869{
870 //
871 // Default Constructor
872 //
873}
874
875//____________________________________________________________
876AliHFEV0pid::AliHFEV0pidTrackIndex::~AliHFEV0pidTrackIndex(){
877 //
878 // Destructor
879 //
880 if(fIndexElectron) delete[] fIndexElectron;
881 if(fIndexPionK0) delete[] fIndexPionK0;
882 if(fIndexPionL) delete[] fIndexPionL;
883 if(fIndexProton) delete[] fIndexProton;
884}
885
886//____________________________________________________________
887void AliHFEV0pid::AliHFEV0pidTrackIndex::Flush(){
888 //
889 // Reset containers
890 //
891
892 if(fIndexElectron) delete[] fIndexElectron;
893 fIndexElectron = NULL;
894 if(fIndexPionK0) delete[] fIndexPionK0;
895 fIndexPionK0 = NULL;
896 if(fIndexPionL) delete[] fIndexPionL;
897 fIndexPionL = NULL;
898 if(fIndexKaon) delete[] fIndexKaon;
899 fIndexKaon = NULL;
900 if(fIndexProton) delete[] fIndexProton;
901 fIndexProton = NULL;
902
903 fNElectrons = 0;
904 fNPionsK0 = 0;
905 fNPionsL = 0;
906 fNKaons = 0;
907 fNProtons = 0;
908}
909
910//____________________________________________________________
911void AliHFEV0pid::AliHFEV0pidTrackIndex::Init(Int_t capacity){
912 //
913 // Initialize container
914 //
915 fIndexElectron = new Int_t[capacity];
916 fIndexPionK0 = new Int_t[capacity];
917 fIndexPionL = new Int_t[capacity];
918 fIndexProton = new Int_t[capacity];
919}
920
921//____________________________________________________________
922void AliHFEV0pid::AliHFEV0pidTrackIndex::Add(Int_t index, Int_t species){
923 //
924 // Add new index to the list of identified particles
925 //
926 switch(species){
927 case AliHFEV0pid::kRecoElectron:
928 fIndexElectron[fNElectrons++] = index;
929 break;
930 case AliHFEV0pid::kRecoPionK0:
931 fIndexPionK0[fNPionsK0++] = index;
932 break;
933 case AliHFEV0pid::kRecoPionL:
934 fIndexPionL[fNPionsL++] = index;
935 break;
936 case AliHFEV0pid::kRecoProton:
937 fIndexProton[fNProtons++] = index;
938 break;
939 };
940}
941
942//____________________________________________________________
943Bool_t AliHFEV0pid::AliHFEV0pidTrackIndex::Find(Int_t index, Int_t species) const {
944 //
945 // Find track index in the specific sample of particles
946 //
947
948 Int_t *container = NULL; Int_t n = 0;
949 switch(species){
950 case AliHFEV0pid::kRecoElectron:
951 container = fIndexElectron;
952 n = fNElectrons;
953 break;
954 case AliHFEV0pid::kRecoPionK0:
955 container = fIndexPionK0;
956 n = fNPionsK0;
957 break;
958 case AliHFEV0pid::kRecoPionL:
959 container = fIndexPionL;
960 n = fNPionsL;
961 break;
962 case AliHFEV0pid::kRecoProton:
963 container = fIndexProton;
964 n = fNProtons;
965 break;
966 }
967 if(!container) return kFALSE;
968 if(n == 0) return kFALSE;
969 Bool_t found = kFALSE;
970 for(Int_t i = 0; i < n; i++){
971 if(container[i] == index){
972 found = kTRUE;
973 break;
974 }
975 }
976 return found;
977}
978
979//____________________________________________________________
980Bool_t AliHFEV0pid::AliHFEV0pidTrackIndex::Find(Int_t index) const {
981 //
982 // Find index in all samples
983 //
984 if(Find(index, AliHFEV0pid::kRecoElectron)) return kTRUE;
985 else if(Find(index, AliHFEV0pid::kRecoPionK0)) return kTRUE;
986 else if(Find(index, AliHFEV0pid::kRecoPionL)) return kTRUE;
987 else return Find(index, AliHFEV0pid::kRecoProton);
988}
989
990//____________________________________________________________
991AliHFEV0pid::AliHFELambdaInf::AliHFELambdaInf(AliKFParticle *mother, AliKFVertex * const primaryVertex):
992 fInvariantMass(0.),
993 fChi2NDF(0.),
994 fDistVtx(0.)
995{
996 //
997 // Constructor
998 // Fill infos
999 //
1000 const Double_t kL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(); // PDG lambda mass
1001 fInvariantMass = mother->GetMass();
1002 // distance frm the privary vertex and mass difference for the (A) case
1003 fDistVtx = mother->GetDistanceFromVertex(*primaryVertex);
1004 // apply constraints for the fit
1005 mother->SetMassConstraint(kL0mass, 0.);
1006 mother->SetProductionVertex(*primaryVertex);
1007 fChi2NDF = mother->GetChi2()/mother->GetNDF();
1008}
1009//____________________________________________________________
1010AliHFEV0pid::AliHFELambdaInf::~AliHFELambdaInf(){
1011 //
1012 // destructor
1013 //
1014}
1015//____________________________________________________________
1016// DEBUG
1017//____________________________________________________________
1018Double_t AliHFEV0pid::GetEffMass(AliESDv0 *v0, UInt_t p1, UInt_t p2) const{
1019 //
1020 // TEMPORARY - this function should become obsolete with v4-18-Rev-10 or 11
1021 // calculate effective mass
1022 //
1023 const Double_t kpmass[5] = {TDatabasePDG::Instance()->GetParticle(kElectron)->Mass(),
1024 TDatabasePDG::Instance()->GetParticle(kMuonMinus)->Mass(),
1025 TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass(),
1026 TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass(),
1027 TDatabasePDG::Instance()->GetParticle(kProton)->Mass()}; // float
1028 if (p1>4) return -1;
1029 if (p2>4) return -1;
1030 Double_t mass1 = kpmass[p1]; // float
1031 Double_t mass2 = kpmass[p2]; // float
1032
1033 Double_t pMom[3];
1034 Double_t nMom[3];
1035
1036 v0->GetPPxPyPz(pMom[0], pMom[1], pMom[2]);
1037 v0->GetNPxPyPz(nMom[0], nMom[1], nMom[2]);
1038
1039
1040 const Double_t *m1 = pMom;
1041 const Double_t *m2 = nMom;
1042 //
1043 //if (fRP[p1]+fRM[p2]<fRP[p2]+fRM[p1]){
1044 // m1 = fPM;
1045 // m2 = fPP;
1046 //}
1047 //
1048 Double_t e1 = TMath::Sqrt(mass1*mass1+
1049 m1[0]*m1[0]+
1050 m1[1]*m1[1]+
1051 m1[2]*m1[2]); // float
1052 Double_t e2 = TMath::Sqrt(mass2*mass2+
1053 m2[0]*m2[0]+
1054 m2[1]*m2[1]+
1055 m2[2]*m2[2]); // float
1056 Double_t mass =
1057 (m2[0]+m1[0])*(m2[0]+m1[0])+
1058 (m2[1]+m1[1])*(m2[1]+m1[1])+
1059 (m2[2]+m1[2])*(m2[2]+m1[2]); // float
1060
1061
1062 mass = (e1+e2)*(e1+e2)-mass;
1063 //if(mass < 0.00001){
1064 // printf("-D: mass: %f\n", mass);
1065 // }
1066 mass = TMath::Sqrt(mass);
1067 return mass;
1068}