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