Line data Source code
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 : /*
17 : * author: M.Kalisky@gsi.de
18 : * 08/Dec/2010
19 : *
20 : * Description: This class allows with purely kinematical cuts
21 : * to select clean samples of electrons, pions and protons from the
22 : * V0 online finder ESD V0 candidates for PID and dectector resonse
23 : * studies.
24 : */
25 :
26 : #include <TVector3.h>
27 : #include <TDatabasePDG.h>
28 :
29 : #include "AliESDv0.h"
30 : #include "AliESDtrack.h"
31 : #include "AliESDEvent.h"
32 : #include "AliVEvent.h"
33 : #include "AliLog.h"
34 : #include "AliKFParticle.h"
35 : #include "AliVTrack.h"
36 : #include "AliKFVertex.h"
37 :
38 : #include "AliESDv0KineCuts.h"
39 :
40 170 : ClassImp(AliESDv0KineCuts)
41 :
42 : //____________________________________________________________________
43 0 : AliESDv0KineCuts::AliESDv0KineCuts() :
44 0 : fEvent(0x0)
45 0 : , fPrimaryVertex(0x0)
46 0 : , fType(0)
47 0 : , fMode(0)
48 0 : , fTPCNcls(1)
49 0 : , fTPCrefit(kTRUE)
50 0 : , fTPCchi2perCls(4.0)
51 0 : , fTPCclsRatio(0.6)
52 0 : , fNoKinks(kTRUE)
53 0 : , fGcutChi2NDF(10)
54 0 : , fGcutInvMass(0.05)
55 0 : , fK0cutChi2NDF(10)
56 0 : , fLcutChi2NDF(10)
57 0 : , fUseExternalVertex(kFALSE)
58 0 : , fDeleteVertex(kFALSE)
59 0 : {
60 : //
61 : // Default constructor
62 : //
63 :
64 : // default single track cuts
65 0 : fTPCNcls = 1; // minimal number of the TPC clusters
66 0 : fTPCrefit = kTRUE; // TPC refit
67 0 : fTPCchi2perCls = 4.0; // chi2 per TPC cluster
68 0 : fTPCclsRatio = 0.6; // minimal foun/findable TPC cluster ratio
69 0 : fNoKinks = kTRUE; // kinks - no [kTRUE] or do not care [kFalse]
70 :
71 :
72 : // default gamma cuts values
73 0 : fGcutChi2NDF = 10; // Chi2NF cut value for the AliKFparticle gamma
74 0 : fGcutCosPoint[0] = 0; // cos of the pointing angle [min, max]
75 0 : fGcutCosPoint[1] = 0.02; // cos of the pointing angle [min, max]
76 0 : fGcutDCA[0] = 0.; // DCA between the daughter tracks [min, max]
77 0 : fGcutDCA[1] = 0.25; // DCA between the daughter tracks [min, max]
78 0 : fGcutVertexR[0] = 3.; // radius of the conversion point [min, max]
79 0 : fGcutVertexR[1] = 90.; // radius of the conversion point [min, max]
80 0 : fGcutPsiPair[0] = 0.; // value of the psi pair cut [min, max]
81 0 : fGcutPsiPair[1] = 0.05; // value of the psi pair cut [min, max]
82 0 : fGcutInvMass = 0.05; // upper value on the gamma invariant mass
83 : // default K0 cuts
84 0 : fK0cutChi2NDF = 10; // Chi2NF cut value for the AliKFparticle K0
85 0 : fK0cutCosPoint[0] = 0.; // cos of the pointing angle [min, max]
86 0 : fK0cutCosPoint[1] = 0.02; // cos of the pointing angle [min, max]
87 0 : fK0cutDCA[0] = 0.; // DCA between the daughter tracks [min, max]
88 0 : fK0cutDCA[1] = 0.2; // DCA between the daughter tracks [min, max]
89 0 : fK0cutVertexR[0] = 2.0; // radius of the decay point [min, max]
90 0 : fK0cutVertexR[1] = 30.0; // radius of the decay point [min, max]
91 0 : fK0cutInvMass[0] = 0.486; // invariant mass window
92 0 : fK0cutInvMass[1] = 0.508; // invariant mass window
93 : // Lambda & anti-Lambda cut values
94 0 : fLcutChi2NDF = 10; // Chi2NF cut value for the AliKFparticle K0
95 0 : fLcutCosPoint[0] = 0.; // cos of the pointing angle [min, max]
96 0 : fLcutCosPoint[1] = 0.02; // cos of the pointing angle [min, max]
97 0 : fLcutDCA[0] = 0.; // DCA between the daughter tracks [min, max]
98 0 : fLcutDCA[1] = 0.2; // DCA between the daughter tracks [min, max]
99 0 : fLcutVertexR[0] = 2.0; // radius of the decay point [min, max]
100 0 : fLcutVertexR[1] = 40.0; // radius of the decay point [min, max]
101 0 : fLcutInvMass[0] = 1.11; // invariant mass window
102 0 : fLcutInvMass[1] = 1.12; // invariant mass window
103 :
104 0 : }
105 : //____________________________________________________________________
106 0 : AliESDv0KineCuts::~AliESDv0KineCuts(){
107 : //
108 : // Destructor
109 : //
110 :
111 :
112 0 : }
113 : //____________________________________________________________________
114 : AliESDv0KineCuts::AliESDv0KineCuts(const AliESDv0KineCuts &ref):
115 0 : TObject(ref)
116 0 : , fEvent(0x0)
117 0 : , fPrimaryVertex(0x0)
118 0 : , fType(0)
119 0 : , fMode(0)
120 0 : , fTPCNcls(1)
121 0 : , fTPCrefit(kTRUE)
122 0 : , fTPCchi2perCls(4.0)
123 0 : , fTPCclsRatio(0.6)
124 0 : , fNoKinks(kTRUE)
125 0 : , fGcutChi2NDF(10)
126 0 : , fGcutInvMass(0.05)
127 0 : , fK0cutChi2NDF(10)
128 0 : , fLcutChi2NDF(10)
129 0 : , fUseExternalVertex(kFALSE)
130 0 : , fDeleteVertex(kFALSE)
131 0 : {
132 : //
133 : // Copy operator
134 : //
135 :
136 0 : ref.Copy(*this);
137 0 : }
138 : //____________________________________________________________________
139 : AliESDv0KineCuts &AliESDv0KineCuts::operator=(const AliESDv0KineCuts &ref){
140 : //
141 : // assignment operator
142 : //
143 0 : if(this != &ref)
144 0 : ref.Copy(*this);
145 0 : return *this;
146 : }
147 : //____________________________________________________________________
148 : void AliESDv0KineCuts::Copy(TObject &ref) const {
149 : //
150 : // Performs the copying of the object
151 : //
152 :
153 0 : TObject::Copy(ref);
154 :
155 0 : AliESDv0KineCuts &target = dynamic_cast<AliESDv0KineCuts &>(ref);
156 :
157 : // default single track cuts
158 0 : target.fTPCNcls = fTPCNcls;
159 0 : target.fTPCrefit = fTPCrefit;
160 0 : target.fTPCchi2perCls = fTPCchi2perCls;
161 0 : target.fTPCclsRatio = fTPCclsRatio;
162 0 : target.fNoKinks = fNoKinks;
163 0 : target.fUseExternalVertex = fUseExternalVertex; //added december 2nd 2011
164 0 : target.fDeleteVertex = fDeleteVertex; //added december 2nd 2011
165 :
166 : // default gamma cuts values
167 0 : target.fGcutChi2NDF = fGcutChi2NDF;
168 0 : memcpy(target.fGcutCosPoint, fGcutCosPoint, sizeof(Float_t) * 2);
169 0 : memcpy(target.fGcutDCA, fGcutDCA, sizeof(Float_t) * 2);
170 0 : memcpy(target.fGcutVertexR, fGcutVertexR, sizeof(Float_t) * 2);
171 0 : memcpy(target.fGcutPsiPair, fGcutPsiPair, sizeof(Float_t) * 2);
172 0 : target.fGcutInvMass = fGcutInvMass;
173 : // default K0 cuts
174 0 : target.fK0cutChi2NDF = fK0cutChi2NDF;
175 0 : memcpy(target.fK0cutCosPoint, fK0cutCosPoint, sizeof(Float_t) * 2);
176 0 : memcpy(target.fK0cutDCA, fK0cutDCA, sizeof(Float_t) * 2);
177 0 : memcpy(target.fK0cutVertexR, fK0cutVertexR, sizeof(Float_t) * 2);
178 0 : memcpy(target.fK0cutInvMass, fK0cutInvMass, sizeof(Float_t) * 2);
179 : // Lambda & anti-Lambda cut values
180 0 : target.fLcutChi2NDF = fLcutChi2NDF;
181 0 : memcpy(target.fLcutCosPoint, fLcutCosPoint, sizeof(Float_t) * 2);
182 0 : memcpy(target.fLcutDCA, fLcutDCA, sizeof(Float_t) * 2);
183 0 : memcpy(target.fLcutVertexR, fLcutVertexR, sizeof(Float_t) * 2);
184 0 : memcpy(target.fLcutInvMass, fLcutInvMass, sizeof(Float_t) * 2);
185 :
186 0 : }
187 : //____________________________________________________________________
188 : Bool_t AliESDv0KineCuts::ProcessV0(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN) const
189 : {
190 : //
191 : // main user function
192 : //
193 :
194 0 : if(!v0) return kFALSE;
195 0 : if(!fEvent){
196 0 : AliErrorClass("No valid Event pointer available, provide it first");
197 0 : return kFALSE;
198 : }
199 :
200 0 : if(!V0CutsCommon(v0)) return kFALSE;
201 :
202 0 : const Int_t id = PreselectV0(v0);
203 :
204 0 : if(!SingleTrackCuts(v0)) return kFALSE;
205 :
206 0 : switch(id){
207 : case kUndef:
208 0 : return kFALSE;
209 : case kGamma:
210 0 : return CaseGamma(v0, pdgV0, pdgP, pdgN);
211 : case kK0:
212 0 : return CaseK0(v0, pdgV0, pdgP, pdgN);
213 : case kLambda:
214 0 : return CaseLambda(v0, pdgV0, pdgP, pdgN, 0);
215 : case kALambda:
216 0 : return CaseLambda(v0, pdgV0, pdgP, pdgN, 1);
217 : default:
218 0 : return kFALSE;
219 : }
220 :
221 : return kFALSE;
222 0 : }
223 : //____________________________________________________________________
224 : Bool_t AliESDv0KineCuts::ProcessV0(AliESDv0* const v0, Int_t &pdgP, Int_t &pdgN) const
225 : {
226 : //
227 : // main user function, simplified if the V0 identity is not necessary
228 : //
229 :
230 0 : if(!v0) return kFALSE;
231 0 : if(!fEvent){
232 0 : AliErrorClass("No valid Event pointer available, provide it first");
233 0 : return kFALSE;
234 : }
235 :
236 0 : Int_t idV0 = -1;
237 0 : return ProcessV0(v0, idV0, pdgP, pdgN);
238 :
239 0 : }
240 : //____________________________________________________________________
241 : Int_t AliESDv0KineCuts::PreselectV0(AliESDv0* const v0) const
242 : {
243 : //
244 : // Make a preselection (exclusive) of the V0 cadidates based on
245 : // Armenteros plot
246 : // the armenteros cut values are currently fixed and user is not able to set them via
247 : // set funcions. The reason is that these cuts are optimized and furneter changes should
248 : // not be necessary. To prove otherwise please study in detail before changing the values
249 : //
250 :
251 0 : Float_t ap[2] = {-1., -1.};
252 0 : Armenteros(v0, ap);
253 : // for clarity
254 0 : const Float_t alpha = ap[0];
255 0 : const Float_t qt = ap[1];
256 :
257 : // selection cuts
258 : // - the reagions for different candidates must not overlap
259 :
260 : // Gamma cuts
261 : const Double_t cutAlphaG = 0.35;
262 : const Double_t cutQTG = 0.05;
263 : const Double_t cutAlphaG2[2] = {0.6, 0.8};
264 : const Double_t cutQTG2 = 0.04;
265 :
266 : // K0 cuts
267 : const Float_t cutQTK0[2] = {0.1075, 0.215};
268 : const Float_t cutAPK0[2] = {0.199, 0.8}; // parameters for curved QT cut
269 :
270 : // Lambda & A-Lambda cuts
271 : const Float_t cutQTL = 0.03;
272 : const Float_t cutAlphaL[2] = {0.35, 0.7};
273 : const Float_t cutAlphaAL[2] = {-0.7, -0.35};
274 : const Float_t cutAPL[3] = {0.107, -0.69, 0.5}; // parameters fir curved QT cut
275 :
276 :
277 0 : if(kPurity == fMode){
278 : // Check for Gamma candidates
279 0 : if(qt < cutQTG){
280 0 : if( (TMath::Abs(alpha) < cutAlphaG) ) return kGamma;
281 : }
282 : // additional region - should help high pT gammas
283 0 : if(qt < cutQTG2){
284 0 : if( (TMath::Abs(alpha) > cutAlphaG2[0]) && (TMath::Abs(alpha) < cutAlphaG2[1]) ) return kGamma;
285 : }
286 : }
287 0 : if(kEffGamma == fMode){
288 0 : if(qt < cutQTG) return kGamma;
289 : }
290 :
291 :
292 : // Check for K0 candidates
293 0 : Float_t q = cutAPK0[0] * TMath::Sqrt(TMath::Abs(1 - alpha*alpha/(cutAPK0[1]*cutAPK0[1])));
294 0 : if( (qt > cutQTK0[0]) && (qt < cutQTK0[1]) && (qt > q) ){
295 0 : return kK0;
296 : }
297 :
298 : // Check for Lambda candidates
299 0 : q = cutAPL[0] * TMath::Sqrt(TMath::Abs(1 - ( (alpha + cutAPL[1]) * (alpha + cutAPL[1]) ) / (cutAPL[2]*cutAPL[2]) ));
300 0 : if( (alpha > cutAlphaL[0]) && (alpha < cutAlphaL[1]) && (qt > cutQTL) && (qt < q) ){
301 0 : return kLambda;
302 : }
303 :
304 : // Check for A-Lambda candidates
305 0 : q = cutAPL[0] * TMath::Sqrt(TMath::Abs(1 - ( (alpha - cutAPL[1]) * (alpha - cutAPL[1]) ) / (cutAPL[2]*cutAPL[2]) ));
306 0 : if( (alpha > cutAlphaAL[0]) && (alpha < cutAlphaAL[1]) && (qt > cutQTL) && (qt < q) ){
307 0 : return kALambda;
308 : }
309 :
310 0 : return kUndef;
311 0 : }
312 : //____________________________________________________________________
313 : Bool_t AliESDv0KineCuts::SingleTrackCuts(AliESDv0 * const v0) const
314 : {
315 : //
316 : // apply single track cuts
317 : // correct sign not relevat here
318 : //
319 :
320 0 : if(!v0) return kFALSE;
321 :
322 : Int_t pIndex = 0, nIndex = 0;
323 0 : pIndex = v0->GetPindex();
324 0 : nIndex = v0->GetNindex();
325 0 : AliESDtrack* d[2];
326 0 : d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));
327 0 : d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));
328 :
329 0 : for(Int_t i=0; i<2; ++i){
330 0 : if(!d[i]) return kFALSE;
331 :
332 : // status word
333 0 : ULong_t status = d[i]->GetStatus();
334 :
335 : // No. of TPC clusters leave to the users
336 0 : if(d[i]->GetTPCNcls() < 1) return kFALSE;
337 :
338 : // TPC refit
339 0 : if(!(status & AliESDtrack::kTPCrefit)) return kFALSE;
340 :
341 : // Chi2 per TPC cluster
342 0 : Int_t nTPCclusters = d[i]->GetTPCNcls();
343 0 : Float_t chi2perTPCcluster = d[i]->GetTPCchi2()/Float_t(nTPCclusters);
344 0 : if(chi2perTPCcluster > 4) return kFALSE;
345 :
346 : // TPC cluster ratio
347 0 : Float_t cRatioTPC = d[i]->GetTPCNclsF() > 0. ? static_cast<Float_t>(d[i]->GetTPCNcls())/static_cast<Float_t> (d[i]->GetTPCNclsF()) : 1.;
348 0 : if(cRatioTPC < 0.6) return kFALSE;
349 :
350 : // kinks
351 0 : if(d[i]->GetKinkIndex(0) != 0) return kFALSE;
352 :
353 0 : }
354 :
355 0 : return kTRUE;
356 0 : }
357 : //____________________________________________________________________
358 : Bool_t AliESDv0KineCuts::CaseGamma(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN) const
359 : {
360 : //
361 : // process the gamma conversion candidate
362 : //
363 :
364 0 : if(!v0) return kFALSE;
365 :
366 : AliVTrack* daughter[2];
367 : Int_t pIndex = 0, nIndex = 0;
368 :
369 0 : Bool_t sign = CheckSigns(v0);
370 0 : if(sign){
371 0 : pIndex = v0->GetPindex();
372 0 : nIndex = v0->GetNindex();
373 0 : }
374 : else{
375 0 : pIndex = v0->GetNindex();
376 0 : nIndex = v0->GetPindex();
377 : }
378 0 : daughter[0] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(pIndex));
379 0 : daughter[1] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(nIndex));
380 0 : if(!daughter[0] || !daughter[1]) return kFALSE;
381 :
382 0 : AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kElectron), TMath::Abs(kElectron));
383 0 : if(!kfMother) return kFALSE;
384 :
385 : AliESDtrack* d[2];
386 0 : d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));
387 0 : d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));
388 :
389 0 : Float_t iMass = v0->GetEffMass(0, 0);
390 :
391 : // cos pointing angle
392 0 : Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
393 0 : cosPoint = TMath::ACos(cosPoint);
394 :
395 : // DCA between daughters
396 0 : Double_t dca = v0->GetDcaV0Daughters();
397 :
398 : // Production vertex
399 0 : Double_t x, y, z;
400 0 : v0->GetXYZ(x,y,z);
401 0 : Double_t r = TMath::Sqrt(x*x + y*y);
402 :
403 0 : Double_t xy[2];
404 : Double_t r2 = -1.;
405 0 : if ( GetConvPosXY(d[0], d[1], xy) ){
406 0 : r2 = TMath::Sqrt(xy[0]*xy[0] + xy[1]*xy[1]);
407 0 : }
408 :
409 : // psi pair
410 0 : Double_t psiPair = PsiPair(v0);
411 :
412 : // V0 chi2/ndf
413 0 : Double_t chi2ndf = kfMother->GetChi2()/kfMother->GetNDF();
414 :
415 0 : if(kfMother) delete kfMother;
416 :
417 : // apply the cuts
418 :
419 0 : if(iMass > fGcutInvMass) return kFALSE;
420 :
421 0 : if(chi2ndf > fGcutChi2NDF) return kFALSE;
422 :
423 0 : if(cosPoint < fGcutCosPoint[0] || cosPoint > fGcutCosPoint[1]) return kFALSE;
424 :
425 0 : if(dca < fGcutDCA[0] || dca > fGcutDCA[1]) return kFALSE;
426 :
427 0 : if(r < fGcutVertexR[0] || r > fGcutVertexR[1]) return kFALSE;
428 :
429 0 : if(psiPair < fGcutPsiPair[0] || psiPair > fGcutPsiPair[1]) return kFALSE;
430 :
431 : // all cuts passed
432 :
433 0 : pdgV0 = 22;
434 0 : if(sign){
435 0 : pdgP = -11;
436 0 : pdgN = 11;
437 0 : }
438 : else{
439 0 : pdgP = 11;
440 0 : pdgN = -11;
441 : }
442 :
443 0 : return kTRUE;
444 0 : }
445 : //____________________________________________________________________
446 : Bool_t AliESDv0KineCuts::CaseK0(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN) const {
447 : //
448 : // process the K0 candidate
449 : //
450 :
451 0 : if(!v0) return kFALSE;
452 :
453 : AliVTrack* daughter[2];
454 : Int_t pIndex = 0, nIndex = 0;
455 0 : Bool_t sign = CheckSigns(v0);
456 0 : if(sign){
457 0 : pIndex = v0->GetPindex();
458 0 : nIndex = v0->GetNindex();
459 0 : }
460 : else{
461 0 : pIndex = v0->GetNindex();
462 0 : nIndex = v0->GetPindex();
463 : }
464 :
465 0 : daughter[0] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(pIndex));
466 0 : daughter[1] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(nIndex));
467 0 : if(!daughter[0] || !daughter[1]) return kFALSE;
468 :
469 0 : AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kPiPlus));
470 0 : if(!kfMother) return kFALSE;
471 :
472 : AliESDtrack* d[2];
473 0 : d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));
474 0 : d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));
475 :
476 0 : Float_t iMass = v0->GetEffMass(2, 2);
477 :
478 : // cos pointing angle
479 0 : Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
480 0 : cosPoint = TMath::ACos(cosPoint);
481 :
482 : // DCA between daughters
483 0 : Double_t dca = v0->GetDcaV0Daughters();
484 :
485 : // Production vertex
486 0 : Double_t x, y, z;
487 0 : v0->GetXYZ(x,y,z);
488 :
489 0 : Double_t r = TMath::Sqrt(x*x + y*y);
490 :
491 : // V0 chi2/ndf
492 0 : Double_t chi2ndf = kfMother->GetChi2()/kfMother->GetNDF();
493 :
494 0 : if(kfMother) delete kfMother;
495 :
496 : //
497 : // apply the cuts
498 : //
499 0 : if(iMass < fK0cutInvMass[0] || iMass > fK0cutInvMass[1]) return kFALSE;
500 :
501 0 : if(chi2ndf > fK0cutChi2NDF) return kFALSE;
502 :
503 0 : if(cosPoint < fK0cutCosPoint[0] || cosPoint > fK0cutCosPoint[1]) return kFALSE;
504 :
505 0 : if(dca < fK0cutDCA[0] || dca > fK0cutDCA[1]) return kFALSE;
506 :
507 0 : if(r < fK0cutVertexR[0] || r > fK0cutVertexR[1]) return kFALSE;
508 :
509 : // all cuts passed
510 0 : pdgV0 = 310;
511 0 : if(sign){
512 0 : pdgP = 211;
513 0 : pdgN = -211;
514 0 : }
515 : else{
516 0 : pdgP = -211;
517 0 : pdgN = 211;
518 : }
519 :
520 0 : return kTRUE;
521 0 : }
522 : //____________________________________________________________________
523 : Bool_t AliESDv0KineCuts::CaseLambda(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN, Int_t id) const {
524 : //
525 : // process teh Lambda and Anti-Lambda candidate
526 : //
527 :
528 0 : if(!v0) return kFALSE;
529 :
530 0 : const Double_t cL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(); // PDG lambda mass
531 :
532 : AliVTrack* daughter[2];
533 : Int_t pIndex = 0, nIndex = 0;
534 : Float_t mMass[2] = {-1., -1.};
535 0 : Bool_t sign = CheckSigns(v0);
536 0 : if(sign){
537 0 : pIndex = v0->GetPindex();
538 0 : nIndex = v0->GetNindex();
539 0 : mMass[0] = v0->GetEffMass(4, 2);
540 0 : mMass[1] = v0->GetEffMass(2, 4);
541 0 : }
542 : else{
543 0 : pIndex = v0->GetNindex();
544 0 : nIndex = v0->GetPindex();
545 0 : mMass[0] = v0->GetEffMass(2, 4);
546 0 : mMass[1] = v0->GetEffMass(4, 2);
547 : }
548 :
549 0 : daughter[0] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(pIndex));
550 0 : daughter[1] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(nIndex));
551 0 : if(!daughter[0] || !daughter[1]) return kFALSE;
552 :
553 0 : AliKFParticle *kfMother[2] = {0x0, 0x0};
554 : // Lambda
555 0 : kfMother[0] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kProton), TMath::Abs(kPiPlus));
556 0 : if(!kfMother[0]) return kFALSE;
557 :
558 : // Anti-Lambda
559 0 : kfMother[1] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kProton));
560 0 : if(!kfMother[1]) return kFALSE;
561 :
562 0 : Float_t dMass[2] = {static_cast<Float_t>(TMath::Abs(mMass[0] - cL0mass)), static_cast<Float_t>(TMath::Abs(mMass[1] - cL0mass))};
563 :
564 : AliESDtrack* d[2];
565 0 : d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));
566 0 : d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));
567 0 : if(!d[0] || !d[1]) return kFALSE;
568 :
569 0 : Float_t p[2] = {static_cast<Float_t>(d[0]->GetP()), static_cast<Float_t>(d[1]->GetP())};
570 :
571 : // check the 3 lambda - antilambda variables
572 : Int_t check[2] = {-1, -1}; // 0 : lambda, 1 : antilambda
573 : // 1) momentum of the daughter particles - proton is expected to have higher momentum than pion
574 0 : check[0] = (p[0] > p[1]) ? 0 : 1;
575 : // 2) mass of the mother particle
576 0 : check[1] = (dMass[0] < dMass[1]) ? 0 : 1;
577 :
578 : // require positive correlation of (1) and (2)
579 0 : if(check[0] != check[1]){
580 0 : if(kfMother[0]) delete kfMother[0];
581 0 : if(kfMother[1]) delete kfMother[1];
582 0 : return kFALSE;
583 : }
584 :
585 : // now that the check[0] == check[1]
586 : const Int_t type = check[0];
587 :
588 : // require that the input armenteros preselection agree:
589 0 : if(type != id) return kFALSE;
590 :
591 : Float_t iMass =0.;
592 0 : if(sign){
593 0 : iMass = (type == 0) ? v0->GetEffMass(4, 2) : v0->GetEffMass(2, 4);
594 0 : }
595 : else{
596 0 : iMass = (type == 0) ? v0->GetEffMass(2, 4) : v0->GetEffMass(4, 2);
597 : }
598 :
599 : // cos pointing angle
600 0 : Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
601 0 : cosPoint = TMath::ACos(cosPoint);
602 :
603 : // DCA between daughters
604 0 : Double_t dca = v0->GetDcaV0Daughters();
605 :
606 : // Production vertex
607 0 : Double_t x, y, z;
608 0 : v0->GetXYZ(x,y,z);
609 0 : Double_t r = TMath::Sqrt(x*x + y*y);
610 :
611 : // proton - pion indices
612 : Int_t ix[2] = {0, 1};
613 0 : if(1 == type){
614 : ix[0] = 1;
615 : ix[1] = 0;
616 0 : }
617 :
618 : // V0 chi2/ndf
619 0 : Double_t chi2ndf = kfMother[type]->GetChi2()/kfMother[type]->GetNDF();
620 :
621 0 : if(kfMother[0]) delete kfMother[0];
622 0 : if(kfMother[1]) delete kfMother[1];
623 :
624 : //
625 : // apply the cuts
626 : //
627 :
628 0 : if(iMass < fLcutInvMass[0] || iMass > fLcutInvMass[1]) return kFALSE;
629 :
630 0 : if(chi2ndf > fLcutChi2NDF) return kFALSE;
631 :
632 0 : if(cosPoint < fLcutCosPoint[0] || cosPoint > fLcutCosPoint[1]) return kFALSE;
633 :
634 0 : if(dca < fLcutDCA[0] || dca > fLcutDCA[1]) return kFALSE;
635 :
636 0 : if(r < fLcutVertexR[0] || r > fLcutVertexR[1]) return kFALSE;
637 :
638 : // all cuts passed
639 :
640 0 : if(0 == type){
641 0 : pdgV0 = 3122;
642 0 : if(sign){
643 0 : pdgP = 2212;
644 0 : pdgN = -211;
645 0 : }
646 : else{
647 0 : pdgP = -211;
648 0 : pdgN = 2212;
649 : }
650 : }
651 : else{
652 0 : pdgV0 = -3122;
653 0 : if(sign){
654 0 : pdgP = 211;
655 0 : pdgN = -2212;
656 0 : }
657 : else{
658 0 : pdgP = -2212;
659 0 : pdgN = 211;
660 : }
661 : }
662 :
663 0 : return kTRUE;
664 0 : }
665 : //____________________________________________________________________
666 : Bool_t AliESDv0KineCuts::V0CutsCommon(const AliESDv0 * const v0) const
667 : {
668 : //
669 : // V0 cuts common to all V0s
670 : //
671 :
672 : AliESDtrack* dN, *dP;
673 :
674 0 : dP = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(v0->GetPindex()));
675 0 : dN = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(v0->GetNindex()));
676 :
677 0 : if(!dN || !dP) return kFALSE;
678 :
679 0 : Int_t qP = dP->Charge();
680 0 : Int_t qN = dN->Charge();
681 :
682 0 : if((qP*qN) != -1) return kFALSE;
683 :
684 0 : return kTRUE;
685 0 : }
686 : //____________________________________________________________________
687 : void AliESDv0KineCuts::Armenteros(AliESDv0* const v0, Float_t val[2]) const
688 : {
689 : //
690 : // computes the Armenteros variables for given V0
691 : // fills the histogram
692 : // returns the values via "val"
693 : //
694 :
695 0 : Double_t mn[3] = {0,0,0};
696 0 : Double_t mp[3] = {0,0,0};
697 0 : Double_t mm[3] = {0,0,0};
698 :
699 0 : if(CheckSigns(v0)){
700 0 : v0->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
701 0 : v0->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
702 0 : }
703 : else{
704 0 : v0->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
705 0 : v0->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
706 : }
707 0 : v0->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
708 :
709 0 : TVector3 vecN(mn[0],mn[1],mn[2]);
710 0 : TVector3 vecP(mp[0],mp[1],mp[2]);
711 0 : TVector3 vecM(mm[0],mm[1],mm[2]);
712 :
713 0 : Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
714 0 : Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
715 :
716 0 : Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
717 0 : ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
718 0 : Double_t qt = vecP.Mag()*sin(thetaP);
719 :
720 0 : val[0] = alfa;
721 0 : val[1] = qt;
722 0 : }
723 : //____________________________________________________________________
724 : Bool_t AliESDv0KineCuts::CheckSigns(AliESDv0* const v0) const
725 : {
726 : //
727 : // check wheter the sign was correctly applied to
728 : // V0 daughter tracks
729 : //
730 :
731 : Bool_t correct = kFALSE;
732 :
733 : Int_t pIndex = 0, nIndex = 0;
734 0 : pIndex = v0->GetPindex();
735 0 : nIndex = v0->GetNindex();
736 :
737 : AliESDtrack* d[2];
738 0 : d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));
739 0 : d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));
740 :
741 : Int_t sign[2];
742 0 : sign[0] = (int)d[0]->GetSign();
743 0 : sign[1] = (int)d[1]->GetSign();
744 :
745 0 : if(-1 == sign[0] && 1 == sign[1]){
746 : correct = kFALSE;
747 0 : }
748 : else{
749 : correct = kTRUE;
750 : }
751 :
752 0 : return correct;
753 : }
754 : //________________________________________________________________
755 : Double_t AliESDv0KineCuts::PsiPair(AliESDv0* const v0) const
756 : {
757 : //
758 : // Angle between daughter momentum plane and plane
759 : //
760 :
761 0 : if(!fEvent) return -1.;
762 :
763 0 : Float_t magField = fEvent->GetMagneticField();
764 :
765 : Int_t pIndex = -1;
766 : Int_t nIndex = -1;
767 0 : if(CheckSigns(v0)){
768 0 : pIndex = v0->GetPindex();
769 0 : nIndex = v0->GetNindex();
770 0 : }
771 : else{
772 0 : pIndex = v0->GetNindex();
773 0 : nIndex = v0->GetPindex();
774 : }
775 :
776 :
777 : AliESDtrack* daughter[2];
778 :
779 0 : daughter[0] = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(pIndex));
780 0 : daughter[1] = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(nIndex));
781 :
782 0 : Double_t x, y, z;
783 0 : v0->GetXYZ(x,y,z);//Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
784 :
785 0 : Double_t mn[3] = {0,0,0};
786 0 : Double_t mp[3] = {0,0,0};
787 :
788 :
789 0 : v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
790 0 : v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
791 :
792 :
793 : Double_t deltat = 1.;
794 0 : 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
795 :
796 0 : Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
797 :
798 0 : Double_t momPosProp[3];
799 0 : Double_t momNegProp[3];
800 :
801 0 : AliExternalTrackParam pt(*daughter[0]), nt(*daughter[1]);
802 :
803 : Double_t psiPair = 4.;
804 :
805 0 : if(nt.PropagateTo(radiussum,magField) == 0)//propagate tracks to the outside
806 0 : psiPair = -5.;
807 0 : if(pt.PropagateTo(radiussum,magField) == 0)
808 0 : psiPair = -5.;
809 0 : pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
810 0 : nt.GetPxPyPz(momNegProp);
811 :
812 : Double_t pEle =
813 0 : TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
814 : Double_t pPos =
815 0 : TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
816 :
817 : Double_t scalarproduct =
818 0 : momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
819 :
820 0 : Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
821 :
822 0 : psiPair = TMath::Abs(TMath::ASin(deltat/chipair));
823 :
824 : return psiPair;
825 0 : }
826 : //___________________________________________________________________
827 : Bool_t AliESDv0KineCuts::GetConvPosXY(AliESDtrack * const ptrack, AliESDtrack * const ntrack, Double_t convpos[2]) const
828 : {
829 : //
830 : // recalculate the gamma conversion XY postition
831 : //
832 :
833 0 : const Double_t b = fEvent->GetMagneticField();
834 :
835 0 : Double_t helixcenterpos[2];
836 0 : GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
837 :
838 0 : Double_t helixcenterneg[2];
839 0 : GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
840 :
841 0 : Double_t poshelix[6];
842 0 : ptrack->GetHelixParameters(poshelix,b);
843 0 : Double_t posradius = TMath::Abs(1./poshelix[4]);
844 :
845 0 : Double_t neghelix[6];
846 0 : ntrack->GetHelixParameters(neghelix,b);
847 0 : Double_t negradius = TMath::Abs(1./neghelix[4]);
848 :
849 0 : Double_t xpos = helixcenterpos[0];
850 0 : Double_t ypos = helixcenterpos[1];
851 0 : Double_t xneg = helixcenterneg[0];
852 0 : Double_t yneg = helixcenterneg[1];
853 :
854 0 : convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
855 0 : convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
856 :
857 0 : return 1;
858 0 : }
859 : //___________________________________________________________________
860 : Bool_t AliESDv0KineCuts::GetHelixCenter(AliESDtrack * const track, Double_t b,Int_t charge, Double_t center[2]) const
861 : {
862 : //
863 : // computes the center of the track helix
864 : //
865 :
866 0 : Double_t pi = TMath::Pi();
867 :
868 0 : Double_t helix[6];
869 0 : track->GetHelixParameters(helix,b);
870 :
871 0 : Double_t xpos = helix[5];
872 0 : Double_t ypos = helix[0];
873 0 : Double_t radius = TMath::Abs(1./helix[4]);
874 0 : Double_t phi = helix[2];
875 :
876 0 : if(phi < 0){
877 0 : phi = phi + 2*pi;
878 0 : }
879 :
880 0 : phi -= pi/2.;
881 0 : Double_t xpoint = radius * TMath::Cos(phi);
882 0 : Double_t ypoint = radius * TMath::Sin(phi);
883 :
884 0 : if(b<0){
885 0 : if(charge > 0){
886 0 : xpoint = - xpoint;
887 0 : ypoint = - ypoint;
888 0 : }
889 : /* avoid self assignment
890 : if(charge < 0){
891 : xpoint = xpoint;
892 : ypoint = ypoint;
893 : }
894 : */
895 : }
896 0 : if(b>0){
897 : /* avoid self assignment
898 : if(charge > 0){
899 : xpoint = xpoint;
900 : ypoint = ypoint;
901 : }
902 : */
903 0 : if(charge < 0){
904 0 : xpoint = - xpoint;
905 0 : ypoint = - ypoint;
906 0 : }
907 : }
908 0 : center[0] = xpos + xpoint;
909 0 : center[1] = ypos + ypoint;
910 :
911 0 : return 1;
912 0 : }
913 : //___________________________________________________________________
914 : AliKFParticle *AliESDv0KineCuts::CreateMotherParticle(const AliVTrack* const pdaughter, const AliVTrack* const ndaughter, Int_t pspec, Int_t nspec) const
915 : {
916 : //
917 : // Creates a mother particle
918 : //
919 0 : AliKFParticle pkfdaughter(*pdaughter, pspec);
920 0 : AliKFParticle nkfdaughter(*ndaughter, nspec);
921 :
922 :
923 : // Create the mother particle
924 0 : AliKFParticle *m = new AliKFParticle(pkfdaughter, nkfdaughter);
925 0 : m->SetField(fEvent->GetMagneticField());
926 0 : if(TMath::Abs(kElectron) == pspec && TMath::Abs(kElectron) == nspec) m->SetMassConstraint(0, 0.001);
927 0 : else if(TMath::Abs(kPiPlus) == pspec && TMath::Abs(kPiPlus) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass(), 0.);
928 0 : else if(TMath::Abs(kProton) == pspec && TMath::Abs(kPiPlus) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(), 0.);
929 0 : else if(TMath::Abs(kPiPlus) == pspec && TMath::Abs(kProton) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(), 0.);
930 : else{
931 0 : AliErrorClass("Wrong daughter ID - mass constraint can not be set");
932 : }
933 :
934 0 : AliKFVertex improvedVertex = *fPrimaryVertex;
935 0 : improvedVertex += *m;
936 0 : m->SetProductionVertex(improvedVertex);
937 :
938 : // update 15/06/2010
939 : // mother particle will not be added to primary vertex but only to its copy
940 : // as this confilcts with calling
941 : // m->SetPrimaryVertex() function and
942 : // subsequently removing the mother particle afterwards
943 : // Source: Sergey Gorbunov
944 :
945 : return m;
946 0 : }
947 : //____________________________________________________________________
948 : void AliESDv0KineCuts::SetEvent(AliESDEvent* const event){
949 : //
950 : // direct setter of ESD event
951 : //
952 0 : fEvent = event;
953 0 : if(!fEvent){
954 0 : AliErrorClass("Invalid input event pointer");
955 0 : return;
956 : }
957 0 : if (fUseExternalVertex) return;
958 : else{
959 0 : if(fPrimaryVertex && fDeleteVertex){
960 0 : delete fPrimaryVertex;
961 0 : fPrimaryVertex=0x0;
962 0 : }
963 0 : fPrimaryVertex = new AliKFVertex(*(fEvent->GetPrimaryVertex()));
964 0 : fDeleteVertex=kTRUE;
965 : }
966 :
967 :
968 :
969 0 : }
970 : //____________________________________________________________________
971 : void AliESDv0KineCuts::SetEvent(AliVEvent* const event){
972 : //
973 : // direct setter of ESD event
974 : //
975 :
976 0 : fEvent = dynamic_cast<AliESDEvent*>(event);
977 0 : if(!fEvent){
978 0 : AliErrorClass("Invalid input event pointer");
979 0 : return;
980 : }
981 :
982 0 : if (fUseExternalVertex) return;
983 : else{
984 0 : if(fPrimaryVertex && fDeleteVertex){
985 0 : delete fPrimaryVertex;
986 0 : fPrimaryVertex=0x0;
987 0 : }
988 0 : fPrimaryVertex = new AliKFVertex(*(fEvent->GetPrimaryVertex()));
989 0 : fDeleteVertex=kTRUE;
990 : }
991 0 : }
992 :
993 :
994 : //________________________________________________________________
995 : void AliESDv0KineCuts::UseExternalVertex(Bool_t use_external){
996 : //
997 : // Reenable primary Vertex from ESD event
998 : //
999 0 : if (use_external) fUseExternalVertex =kTRUE;
1000 0 : else fUseExternalVertex =kFALSE;
1001 0 : }
1002 :
1003 :
1004 :
1005 :
1006 : //________________________________________________________________
1007 : void AliESDv0KineCuts::SetPrimaryVertex(AliKFVertex* const v){
1008 : //
1009 : // set the primary vertex of the event
1010 : //
1011 0 : if(fPrimaryVertex && fDeleteVertex){
1012 0 : delete fPrimaryVertex;
1013 0 : fPrimaryVertex =0x0;
1014 0 : fDeleteVertex = kFALSE;
1015 0 : }
1016 0 : fUseExternalVertex=kTRUE;
1017 0 : fPrimaryVertex = v; // set primary Vertex
1018 0 : if(!fPrimaryVertex){
1019 0 : AliErrorClass("Failed to initialize the primary vertex");
1020 0 : return;
1021 : }
1022 0 : }
1023 : //___________________________________________________________________
1024 : void AliESDv0KineCuts::SetMode(Int_t mode, Int_t type){
1025 : //
1026 : // this function allows the user to select (prior running the 'ProcessV0' function)
1027 : // to select different approaches to V0 selection - the 'mode'
1028 : // - and -
1029 : // different systems (pp, PbPb) - 'type'
1030 : //
1031 : // To see the cut values for different modes please refer to the
1032 : // function SetCuts()
1033 : //
1034 : // Important notice: based on the parameters particular sets of cuts will
1035 : // be activated for teh V0 selection. If some additional changes to single
1036 : // cuts are needed please us the SetXXXcut function (see the header file)
1037 : //
1038 :
1039 0 : switch(mode){
1040 : case kPurity:
1041 0 : fMode = kPurity; // used to obtain highest purity possible - the efficiency may be low
1042 0 : break;
1043 : case kEffGamma:
1044 0 : fMode = kEffGamma; // used to obtain highes efficiency possible - the purity may be worse
1045 0 : break;
1046 : default:
1047 0 : AliError("V0 selection mode not recognozed, setting 'kPurity'");
1048 0 : fMode = kPurity;
1049 0 : }
1050 :
1051 0 : switch(type){
1052 : case kPP:
1053 0 : fType = kPP; // cuts optimized for low multiplicity
1054 0 : break;
1055 : case kPbPb:
1056 0 : fType = kPbPb; // cuts optimized for high multiplicity
1057 0 : break;
1058 : }
1059 :
1060 : // setup the cut values for selected mode & type
1061 0 : SetCuts();
1062 :
1063 0 : }
1064 : //___________________________________________________________________
1065 : void AliESDv0KineCuts::SetMode(Int_t mode, const char* type){
1066 : //
1067 : // overloaded function - please see above
1068 : //
1069 :
1070 : Int_t t = -1;
1071 :
1072 0 : if(!strcmp("pp", type)) t = kPP;
1073 0 : else if(!(strcmp("PbPb", type))) t = kPbPb;
1074 : else{
1075 0 : AliError("data type not recognized, setting 'pp'");
1076 : t = kPP;
1077 : }
1078 :
1079 0 : SetMode(mode, t);
1080 :
1081 0 : }
1082 : //___________________________________________________________________
1083 : void AliESDv0KineCuts::SetCuts(){
1084 : //
1085 : // this funciton sets the default cut values based on the selected
1086 : // fMode and fType.
1087 : // please note that only the cuts that have different values than the default
1088 : // cuts are updated here
1089 : //
1090 :
1091 : // last update: 14/02/2011
1092 : // as a very preliminary - the only change to default cuts is to apply
1093 : // less restricting gamma conversion selection in PreselectV0() function
1094 :
1095 :
1096 :
1097 0 : }
|
