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