]>
Commit | Line | Data |
---|---|---|
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" | |
48 | ClassImp(AliHFEV0pid) | |
49 | ||
50 | //____________________________________________________________ | |
51 | AliHFEV0pid::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 | //____________________________________________________________ | |
75 | AliHFEV0pid::~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 | //____________________________________________________________ | |
92 | void 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 | //____________________________________________________________ | |
146 | void 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 | //____________________________________________________________ | |
196 | Int_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 | //____________________________________________________________ | |
220 | void 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 | //____________________________________________________________ | |
233 | Bool_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 | //____________________________________________________________ | |
299 | Bool_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 | //____________________________________________________________ | |
360 | Bool_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 | //____________________________________________________________ | |
387 | Bool_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 | //____________________________________________________________ | |
528 | AliKFParticle *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 | //____________________________________________________________ | |
570 | void 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 | //____________________________________________________________ | |
596 | Bool_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 | //____________________________________________________________ | |
686 | Bool_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 | //_________________________________________________ | |
743 | Bool_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 | //_________________________________________________ | |
753 | Bool_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 | //_________________________________________________ | |
766 | Bool_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 | //_________________________________________________ | |
780 | Float_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 | //_________________________________________________ | |
798 | Float_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 | //____________________________________________________________ | |
858 | AliHFEV0pid::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 | //____________________________________________________________ | |
876 | AliHFEV0pid::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 | //____________________________________________________________ | |
887 | void 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 | //____________________________________________________________ | |
911 | void 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 | //____________________________________________________________ | |
922 | void 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 | //____________________________________________________________ | |
943 | Bool_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 | //____________________________________________________________ | |
980 | Bool_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 | //____________________________________________________________ | |
991 | AliHFEV0pid::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 | //____________________________________________________________ | |
1010 | AliHFEV0pid::AliHFELambdaInf::~AliHFELambdaInf(){ | |
1011 | // | |
1012 | // destructor | |
1013 | // | |
1014 | } | |
1015 | //____________________________________________________________ | |
1016 | // DEBUG | |
1017 | //____________________________________________________________ | |
1018 | Double_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 | } |