]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliESDv0KineCuts.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Major update of the HFE package (comments inside the code
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliESDv0KineCuts.cxx
CommitLineData
3a72645a 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 dectero 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 "AliLog.h"
33#include "AliKFParticle.h"
34#include "AliVTrack.h"
35#include "AliKFVertex.h"
36
37#include "AliESDv0KineCuts.h"
38
39ClassImp(AliESDv0KineCuts)
40
41//____________________________________________________________________
42AliESDv0KineCuts::AliESDv0KineCuts() :
43 fV0(0x0)
44 , fEvent(0x0)
45 , fPrimaryVertex(0x0)
46{
47 //
48 // Default constructor
49 //
50
51}
52//____________________________________________________________________
53AliESDv0KineCuts::~AliESDv0KineCuts(){
54 //
55 // Destructor
56 //
57
58 if (fV0) delete fV0;
59
60}
61//____________________________________________________________________
62AliESDv0KineCuts::AliESDv0KineCuts(const AliESDv0KineCuts &ref):
63 TObject(ref)
64 , fV0(0x0)
65 , fEvent(0x0)
66 , fPrimaryVertex(0x0)
67{
68 //
69 // Copy operator
70 //
71
72 ref.Copy(*this);
73}
74//____________________________________________________________________
75AliESDv0KineCuts &AliESDv0KineCuts::operator=(const AliESDv0KineCuts &ref){
76 //
77 // assignment operator
78 //
79 if(this != &ref)
80 ref.Copy(*this);
81 return *this;
82}
83//____________________________________________________________________
84void AliESDv0KineCuts::Copy(TObject &ref) const {
85 //
86 // Performs the copying of the object
87 //
88
89 TObject::Copy(ref);
90
91 AliESDv0KineCuts &target = dynamic_cast<AliESDv0KineCuts &>(ref);
92 if(fV0)
93 target.fV0 = dynamic_cast<AliESDv0 *>(fV0->Clone());
94 else
95 target.fV0 = NULL;
96
97}
98//____________________________________________________________________
99Bool_t AliESDv0KineCuts::ProcessV0(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN){
100 //
101 // main user function
102 //
103
104 if(!v0) return kFALSE;
105 if(!fEvent){
106 AliErrorClass("No valid Event pointer available, provide it first");
107 return kFALSE;
108 }
109
110 if(!V0CutsCommon(v0)) return kFALSE;
111
112 const Int_t id = PreselectV0(v0);
113
114 if(!SingleTrackCuts(v0)) return kFALSE;
115
116 switch(id){
117 case kUndef:
118 return kFALSE;
119 case kGamma:
120 return CaseGamma(v0, pdgV0, pdgP, pdgN);
121 case kK0:
122 return CaseK0(v0, pdgV0, pdgP, pdgN);
123 case kLambda:
124 return CaseLambda(v0, pdgV0, pdgP, pdgN, 0);
125 case kALambda:
126 return CaseLambda(v0, pdgV0, pdgP, pdgN, 1);
127 default:
128 return kFALSE;
129 }
130
131 return kFALSE;
132}
133//____________________________________________________________________
134Bool_t AliESDv0KineCuts::ProcessV0(AliESDv0* const v0, Int_t &pdgP, Int_t &pdgN){
135 //
136 // main user function, simplified if the V0 identity is not necessary
137 //
138
139 if(!v0) return kFALSE;
140 if(!fEvent){
141 AliErrorClass("No valid Event pointer available, provide it first");
142 return kFALSE;
143 }
144
145 Int_t idV0 = -1;
146 return ProcessV0(v0, idV0, pdgP, pdgN);
147
148}
149//____________________________________________________________________
150Int_t AliESDv0KineCuts::PreselectV0(AliESDv0* const v0){
151 //
152 // Make a preselection (exclusive) of the V0 cadidates based on
153 // Armenteros plot
154 //
155
156 Float_t ap[2] = {-1., -1.};
157 Armenteros(v0, ap);
158 // for clarity
159 const Float_t alpha = ap[0];
160 const Float_t qt = ap[1];
161
162 // selection cuts
163 // - the reagions for different candidates must not overlap
164
165 // Gamma cuts
166 const Double_t cutAlphaG = 0.35;
167 const Double_t cutAlphaG2[2] = {0.6, 0.8};
168 const Double_t cutQTG = 0.05;
169 const Double_t cutQTG2 = 0.04;
170
171 // K0 cuts
172 const Float_t cutQTK0[2] = {0.1075, 0.215};
173 const Float_t cutAPK0[2] = {0.199, 0.8}; // parameters for curved QT cut
174
175 // Lambda & A-Lambda cuts
176 const Float_t cutQTL = 0.03;
177 const Float_t cutAlphaL[2] = {0.35, 0.7};
178 const Float_t cutAlphaAL[2] = {-0.7, -0.35};
179 const Float_t cutAPL[3] = {0.107, -0.69, 0.5}; // parameters fir curved QT cut
180
181
182 // Check for Gamma candidates
183 if(qt < cutQTG){
184 if( (TMath::Abs(alpha) < cutAlphaG) ) return kGamma;
185 }
186 // additional region - should help high pT gammas
187 if(qt < cutQTG2){
188 if( (TMath::Abs(alpha) > cutAlphaG2[0]) && (TMath::Abs(alpha) < cutAlphaG2[1]) ) return kGamma;
189 }
190
191
192 // Check for K0 candidates
193 Float_t q = cutAPK0[0] * TMath::Sqrt(TMath::Abs(1 - alpha*alpha/(cutAPK0[1]*cutAPK0[1])));
194 if( (qt > cutQTK0[0]) && (qt < cutQTK0[1]) && (qt > q) ){
195 return kK0;
196 }
197
198 // Check for Lambda candidates
199 q = cutAPL[0] * TMath::Sqrt(TMath::Abs(1 - ( (alpha + cutAPL[1]) * (alpha + cutAPL[1]) ) / (cutAPL[2]*cutAPL[2]) ));
200 if( (alpha > cutAlphaL[0]) && (alpha < cutAlphaL[1]) && (qt > cutQTL) && (qt < q) ){
201 return kLambda;
202 }
203
204 // Check for A-Lambda candidates
205 q = cutAPL[0] * TMath::Sqrt(TMath::Abs(1 - ( (alpha - cutAPL[1]) * (alpha - cutAPL[1]) ) / (cutAPL[2]*cutAPL[2]) ));
206 if( (alpha > cutAlphaAL[0]) && (alpha < cutAlphaAL[1]) && (qt > cutQTL) && (qt < q) ){
207 return kALambda;
208 }
209
210 return kUndef;
211}
212//____________________________________________________________________
213Bool_t AliESDv0KineCuts::SingleTrackCuts(AliESDv0 * const v0){
214 //
215 // apply single track cuts
216 // correct sign not relevat here
217 //
218
219 if(!v0) return kFALSE;
220
221 Int_t pIndex = 0, nIndex = 0;
222 pIndex = v0->GetPindex();
223 nIndex = v0->GetNindex();
224 AliESDtrack* d[2];
225 d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));
226 d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));
227
228 for(Int_t i=0; i<2; ++i){
229 if(!d[i]) return kFALSE;
230
231 // status word
232 ULong_t status = d[i]->GetStatus();
233
234 // No. of TPC clusters leave to the users
235 if(d[i]->GetTPCNcls() < 1) return kFALSE;
236
237 // TPC refit
238 if(!(status & AliESDtrack::kTPCrefit)) return kFALSE;
239
240 // Chi2 per TPC cluster
241 Int_t nTPCclusters = d[i]->GetTPCNcls();
242 Float_t chi2perTPCcluster = d[i]->GetTPCchi2()/Float_t(nTPCclusters);
243 if(chi2perTPCcluster > 4) return kFALSE;
244
245 // TPC cluster ratio
246 Float_t cRatioTPC = d[i]->GetTPCNclsF() > 0. ? static_cast<Float_t>(d[i]->GetTPCNcls())/static_cast<Float_t> (d[i]->GetTPCNclsF()) : 1.;
247 if(cRatioTPC < 0.6) return kFALSE;
248
249 // kinks
250 if(d[i]->GetKinkIndex(0) != 0) return kFALSE;
251
252 }
253
254 return kTRUE;
255}
256//____________________________________________________________________
257Bool_t AliESDv0KineCuts::CaseGamma(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN){
258 //
259 // process the gamma conversion candidate
260 //
261
262 if(!v0) return kFALSE;
263
264 AliVTrack* daughter[2];
265 Int_t pIndex = 0, nIndex = 0;
266
267 Bool_t sign = CheckSigns(v0);
268 if(sign){
269 pIndex = v0->GetPindex();
270 nIndex = v0->GetNindex();
271 }
272 else{
273 pIndex = v0->GetNindex();
274 nIndex = v0->GetPindex();
275 }
276 daughter[0] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(pIndex));
277 daughter[1] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(nIndex));
278 if(!daughter[0] || !daughter[1]) return kFALSE;
279
280 AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kElectron), TMath::Abs(kElectron));
281 if(!kfMother) return kFALSE;
282
283 AliESDtrack* d[2];
284 d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));
285 d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));
286
287 Float_t iMass = v0->GetEffMass(0, 0);
288
289 // Cut values
290 const Double_t cutMass = 0.05; // old [0.05]
291 const Double_t cutChi2NDF = 40.; // old [7.]
292 const Double_t cutCosPoint[2] = {0., 0.02}; // old [0., 0.03]
293 const Double_t cutDCA[2] = {0., 0.25}; // old [0., 0.25]
294 const Double_t cutProdVtxR[2] = {8., 90.}; // old [6., 9999]
295 const Double_t cutPsiPair[2] = {0., 0.05}; // old [0. 0.05]
296 const Double_t cutOAngle[2] = {0, 0.1}; // old [0., 0.1]
297
298 // cos pointing angle
299 Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
300 cosPoint = TMath::ACos(cosPoint);
301
302 // DCA between daughters
303 Double_t dca = v0->GetDcaV0Daughters();
304
305 // Production vertex
306 Double_t x, y, z;
307 v0->GetXYZ(x,y,z);
308 Double_t r = TMath::Sqrt(x*x + y*y);
309
310 Double_t xy[2];
311 Double_t r2 = -1.;
312 if ( GetConvPosXY(d[0], d[1], xy) ){
313 r2 = TMath::Sqrt(xy[0]*xy[0] + xy[1]*xy[1]);
314 }
315
316 // Opening angle
317 Double_t oAngle = OpenAngle(v0);
318
319 // psi pair
320 Double_t psiPair = PsiPair(v0);
321
322 // V0 chi2/ndf
323 Double_t chi2ndf = kfMother->GetChi2()/kfMother->GetNDF();
324
325 if(kfMother) delete kfMother;
326
327 // apply the cuts
328
329 if(iMass > cutMass) return kFALSE;
330
331 if(chi2ndf > cutChi2NDF) return kFALSE;
332
333 if(cosPoint < cutCosPoint[0] || cosPoint > cutCosPoint[1]) return kFALSE;
334
335 if(dca < cutDCA[0] || dca > cutDCA[1]) return kFALSE;
336
337 if(r < cutProdVtxR[0] || r > cutProdVtxR[1]) return kFALSE;
338
339 if(psiPair < cutPsiPair[0] || psiPair > cutPsiPair[1]) return kFALSE;
340
341 if(oAngle < cutOAngle[0] || oAngle > cutOAngle[1]) return kFALSE;
342
343 // all cuts passed
344
345 pdgV0 = 22;
346 if(sign){
347 pdgP = -11;
348 pdgN = 11;
349 }
350 else{
351 pdgP = 11;
352 pdgN = -11;
353 }
354
355 return kTRUE;
356}
357//____________________________________________________________________
358Bool_t AliESDv0KineCuts::CaseK0(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN){
359 //
360 // process the K0 candidate
361 //
362
363 if(!v0) return kFALSE;
364
365 AliVTrack* daughter[2];
366 Int_t pIndex = 0, nIndex = 0;
367 Bool_t sign = CheckSigns(v0);
368 if(sign){
369 pIndex = v0->GetPindex();
370 nIndex = v0->GetNindex();
371 }
372 else{
373 pIndex = v0->GetNindex();
374 nIndex = v0->GetPindex();
375 }
376
377 daughter[0] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(pIndex));
378 daughter[1] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(nIndex));
379 if(!daughter[0] || !daughter[1]) return kFALSE;
380
381 AliKFParticle *kfMother = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kPiPlus));
382 if(!kfMother) return kFALSE;
383
384 AliESDtrack* d[2];
385 d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));
386 d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));
387
388 Float_t iMass = v0->GetEffMass(2, 2);
389
390 const Double_t cutMass[2] = {0.486, 0.508}; // ORG [0.485, 0.51]
391 const Double_t cutChi2NDF = 40.; // ORG [7.]
392 const Double_t cutCosPoint[2] = {0., 0.02}; // ORG [0., 0.03]
393 const Double_t cutDCA[2] = {0., 0.2}; // ORG [0., 0.1]
394 const Double_t cutProdVtxR[2] = {2.0, 30.}; // ORG [0., 8.1]
395
396 // cos pointing angle
397 Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
398 cosPoint = TMath::ACos(cosPoint);
399
400 // DCA between daughters
401 Double_t dca = v0->GetDcaV0Daughters();
402
403 // Production vertex
404 Double_t x, y, z;
405 v0->GetXYZ(x,y,z);
406
407 Double_t r = TMath::Sqrt(x*x + y*y);
408
409 // V0 chi2/ndf
410 Double_t chi2ndf = kfMother->GetChi2()/kfMother->GetNDF();
411
412 if(kfMother) delete kfMother;
413
414 //
415 // apply the cuts
416 //
417 if(iMass < cutMass[0] || iMass > cutMass[1]) return kFALSE;
418
419 if(chi2ndf > cutChi2NDF) return kFALSE;
420
421 if(cosPoint < cutCosPoint[0] || cosPoint > cutCosPoint[1]) return kFALSE;
422
423 if(dca < cutDCA[0] || dca > cutDCA[1]) return kFALSE;
424
425 if(r < cutProdVtxR[0] || r > cutProdVtxR[1]) return kFALSE;
426
427 // all cuts passed
428 pdgV0 = 310;
429 if(sign){
430 pdgP = 211;
431 pdgN = -211;
432 }
433 else{
434 pdgP = -211;
435 pdgN = 211;
436 }
437
438 return kTRUE;
439}
440//____________________________________________________________________
441Bool_t AliESDv0KineCuts::CaseLambda(AliESDv0* const v0, Int_t &pdgV0, Int_t &pdgP, Int_t &pdgN, Int_t id){
442 //
443 // process teh Lambda and Anti-Lambda candidate
444 //
445
446 if(!v0) return kFALSE;
447
448 const Double_t cL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(); // PDG lambda mass
449
450 AliVTrack* daughter[2];
451 Int_t pIndex = 0, nIndex = 0;
452 Float_t mMass[2] = {-1., -1.};
453 Bool_t sign = CheckSigns(v0);
454 if(sign){
455 pIndex = v0->GetPindex();
456 nIndex = v0->GetNindex();
457 mMass[0] = v0->GetEffMass(4, 2);
458 mMass[1] = v0->GetEffMass(2, 4);
459 }
460 else{
461 pIndex = v0->GetNindex();
462 nIndex = v0->GetPindex();
463 mMass[0] = v0->GetEffMass(2, 4);
464 mMass[1] = v0->GetEffMass(4, 2);
465 }
466
467 daughter[0] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(pIndex));
468 daughter[1] = dynamic_cast<AliVTrack *>(fEvent->GetTrack(nIndex));
469 if(!daughter[0] || !daughter[1]) return kFALSE;
470
471 AliKFParticle *kfMother[2] = {0x0, 0x0};
472 // Lambda
473 kfMother[0] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kProton), TMath::Abs(kPiPlus));
474 if(!kfMother[0]) return kFALSE;
475
476 // Anti-Lambda
477 kfMother[1] = CreateMotherParticle(daughter[0], daughter[1], TMath::Abs(kPiPlus), TMath::Abs(kProton));
478 if(!kfMother[1]) return kFALSE;
479
480 Float_t dMass[2] = {TMath::Abs(mMass[0] - cL0mass), TMath::Abs(mMass[1] - cL0mass)};
481
482 AliESDtrack* d[2];
483 d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));
484 d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));
485 if(!d[0] || !d[1]) return kFALSE;
486
487 Float_t p[2] = {d[0]->GetP(), d[1]->GetP()};
488
489 // check the 3 lambda - antilambda variables
490 Int_t check[2] = {-1, -1}; // 0 : lambda, 1 : antilambda
491 // 1) momentum of the daughter particles - proton is expected to have higher momentum than pion
492 check[0] = (p[0] > p[1]) ? 0 : 1;
493 // 2) mass of the mother particle
494 check[1] = (dMass[0] < dMass[1]) ? 0 : 1;
495
496 // require positive correlation of (1) and (2)
497 if(check[0] != check[1]){
498 if(kfMother[0]) delete kfMother[0];
499 if(kfMother[1]) delete kfMother[1];
500 return kFALSE;
501 }
502
503 // now that the check[0] == check[1]
504 const Int_t type = check[0];
505
506 // require that the input armenteros preselection agree:
507 if(type != id) return kFALSE;
508
509 Float_t iMass =0.;
510 if(sign){
511 iMass = (type == 0) ? v0->GetEffMass(4, 2) : v0->GetEffMass(2, 4);
512 }
513 else{
514 iMass = (type == 0) ? v0->GetEffMass(2, 4) : v0->GetEffMass(4, 2);
515 }
516
517 // Cuts
518 const Double_t cutMass[2] = {1.11, 1.12}; // ORG [1.11, 1.12]
519 const Double_t cutChi2NDF = 40.; // ORG [5.]
520 const Double_t cutCosPoint[2] = {0., 0.02}; // ORG [0., 0.03]
521 const Double_t cutDCA[2] = {0., 0.2}; // ORG [0., 0.2]
522 const Double_t cutProdVtxR[2] = {2., 40.}; // ORG [0., 24.]
523
524 // cos pointing angle
525 Double_t cosPoint = v0->GetV0CosineOfPointingAngle();
526 cosPoint = TMath::ACos(cosPoint);
527
528 // DCA between daughters
529 Double_t dca = v0->GetDcaV0Daughters();
530
531 // Production vertex
532 Double_t x, y, z;
533 v0->GetXYZ(x,y,z);
534 Double_t r = TMath::Sqrt(x*x + y*y);
535
536 // proton - pion indices
537 Int_t ix[2] = {0, 1};
538 if(1 == type){
539 ix[0] = 1;
540 ix[1] = 0;
541 }
542
543 // V0 chi2/ndf
544 Double_t chi2ndf = kfMother[type]->GetChi2()/kfMother[type]->GetNDF();
545
546 if(kfMother[0]) delete kfMother[0];
547 if(kfMother[1]) delete kfMother[1];
548
549 //
550 // apply the cuts
551 //
552
553 if(iMass < cutMass[0] || iMass > cutMass[1]) return kFALSE;
554
555 if(chi2ndf > cutChi2NDF) return kFALSE;
556
557 if(cosPoint < cutCosPoint[0] || cosPoint > cutCosPoint[1]) return kFALSE;
558
559 if(dca < cutDCA[0] || dca > cutDCA[1]) return kFALSE;
560
561 if(r < cutProdVtxR[0] || r > cutProdVtxR[1]) return kFALSE;
562
563 // all cuts passed
564
565 if(0 == type){
566 pdgV0 = 3122;
567 if(sign){
568 pdgP = 2212;
569 pdgN = -211;
570 }
571 else{
572 pdgP = -211;
573 pdgN = 2212;
574 }
575 }
576 else{
577 pdgV0 = -3122;
578 if(sign){
579 pdgP = 211;
580 pdgN = -2212;
581 }
582 else{
583 pdgP = -2212;
584 pdgN = 211;
585 }
586 }
587
588 return kTRUE;
589}
590//____________________________________________________________________
591Bool_t AliESDv0KineCuts::V0CutsCommon(const AliESDv0 * const v0){
592 //
593 // V0 cuts common to all V0s
594 //
595
596 AliESDtrack* dN, *dP;
597
598 dP = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(v0->GetPindex()));
599 dN = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(v0->GetNindex()));
600
601 if(!dN || !dP) return kFALSE;
602
603 Int_t qP = dP->Charge();
604 Int_t qN = dN->Charge();
605
606 if((qP*qN) != -1) return kFALSE;
607
608 return kTRUE;
609}
610//____________________________________________________________________
611void AliESDv0KineCuts::Armenteros(AliESDv0* const v0, Float_t val[2]){
612 //
613 // computes the Armenteros variables for given V0
614 // fills the histogram
615 // returns the values via "val"
616 //
617
618 Double_t mn[3] = {0,0,0};
619 Double_t mp[3] = {0,0,0};
620 Double_t mm[3] = {0,0,0};
621
622 if(CheckSigns(v0)){
623 v0->GetNPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
624 v0->GetPPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
625 }
626 else{
627 v0->GetPPxPyPz(mn[0],mn[1],mn[2]); //reconstructed cartesian momentum components of negative daughter
628 v0->GetNPxPyPz(mp[0],mp[1],mp[2]); //reconstructed cartesian momentum components of positive daughter
629 }
630 v0->GetPxPyPz(mm[0],mm[1],mm[2]); //reconstructed cartesian momentum components of mother
631
632 TVector3 vecN(mn[0],mn[1],mn[2]);
633 TVector3 vecP(mp[0],mp[1],mp[2]);
634 TVector3 vecM(mm[0],mm[1],mm[2]);
635
636 Double_t thetaP = acos((vecP * vecM)/(vecP.Mag() * vecM.Mag()));
637 Double_t thetaN = acos((vecN * vecM)/(vecN.Mag() * vecM.Mag()));
638
639 Double_t alfa = ((vecP.Mag())*cos(thetaP)-(vecN.Mag())*cos(thetaN))/
640 ((vecP.Mag())*cos(thetaP)+(vecN.Mag())*cos(thetaN)) ;
641 Double_t qt = vecP.Mag()*sin(thetaP);
642
643 val[0] = alfa;
644 val[1] = qt;
645}
646//____________________________________________________________________
647Bool_t AliESDv0KineCuts::CheckSigns(AliESDv0* const v0){
648 //
649 // check wheter the sign was correctly applied to
650 // V0 daughter tracks
651 //
652
653 Bool_t correct = kFALSE;
654
655 Int_t pIndex = 0, nIndex = 0;
656 pIndex = v0->GetPindex();
657 nIndex = v0->GetNindex();
658
659 AliESDtrack* d[2];
660 d[0] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(pIndex));
661 d[1] = dynamic_cast<AliESDtrack*>(fEvent->GetTrack(nIndex));
662
663 Int_t sign[2];
664 sign[0] = (int)d[0]->GetSign();
665 sign[1] = (int)d[1]->GetSign();
666
667 if(-1 == sign[0] && 1 == sign[1]){
668 correct = kFALSE;
669 }
670 else{
671 correct = kTRUE;
672 }
673
674 return correct;
675}
676//____________________________________________________________________
677void AliESDv0KineCuts::SetEvent(AliESDEvent* const event){
678 //
679 // direct setter of ESD event
680 //
681 if(!event){
682 AliErrorClass("Invalid input event pointer");
683 return;
684 }
685 fEvent = event;
686}
687//____________________________________________________________________
688void AliESDv0KineCuts::SetEvent(AliVEvent* const event){
689 //
690 // Set the current working ESD event
691 //
692 if(!event){
693 AliErrorClass("Invalid input event pointer");
694 return;
695 }
696
697 SetEvent(dynamic_cast<AliESDEvent*>(event));
698}
699//________________________________________________________________
700Double_t AliESDv0KineCuts::OpenAngle(AliESDv0 *v0) const {
701 //
702 // Opening angle between two daughter tracks
703 //
704 Double_t mn[3] = {0,0,0};
705 Double_t mp[3] = {0,0,0};
706
707
708 v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
709 v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
710
711
712 Double_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])));
713
714 return TMath::Abs(openAngle);
715}
716//________________________________________________________________
717Double_t AliESDv0KineCuts::PsiPair(AliESDv0* const v0) {
718 //
719 // Angle between daughter momentum plane and plane
720 //
721
722 if(!fEvent) return -1.;
723
724 Float_t magField = fEvent->GetMagneticField();
725
726 Int_t pIndex = -1;
727 Int_t nIndex = -1;
728 if(CheckSigns(v0)){
729 pIndex = v0->GetPindex();
730 nIndex = v0->GetNindex();
731 }
732 else{
733 pIndex = v0->GetNindex();
734 nIndex = v0->GetPindex();
735 }
736
737
738 AliESDtrack* daughter[2];
739
740 daughter[0] = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(pIndex));
741 daughter[1] = dynamic_cast<AliESDtrack *>(fEvent->GetTrack(nIndex));
742
743 Double_t x, y, z;
744 v0->GetXYZ(x,y,z);//Reconstructed coordinates of V0; to be replaced by Markus Rammler's method in case of conversions!
745
746 Double_t mn[3] = {0,0,0};
747 Double_t mp[3] = {0,0,0};
748
749
750 v0->GetNPxPyPz(mn[0],mn[1],mn[2]);//reconstructed cartesian momentum components of negative daughter;
751 v0->GetPPxPyPz(mp[0],mp[1],mp[2]);//reconstructed cartesian momentum components of positive daughter;
752
753
754 Double_t deltat = 1.;
755 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
756
757 Double_t radiussum = TMath::Sqrt(x*x + y*y) + 50;//radius to which tracks shall be propagated
758
759 Double_t momPosProp[3];
760 Double_t momNegProp[3];
761
762 AliExternalTrackParam pt(*daughter[0]), nt(*daughter[1]);
763
764 Double_t psiPair = 4.;
765
766 if(nt.PropagateTo(radiussum,magField) == 0)//propagate tracks to the outside
767 psiPair = -5.;
768 if(pt.PropagateTo(radiussum,magField) == 0)
769 psiPair = -5.;
770 pt.GetPxPyPz(momPosProp);//Get momentum vectors of tracks after propagation
771 nt.GetPxPyPz(momNegProp);
772
773 Double_t pEle =
774 TMath::Sqrt(momNegProp[0]*momNegProp[0]+momNegProp[1]*momNegProp[1]+momNegProp[2]*momNegProp[2]);//absolute momentum value of negative daughter
775 Double_t pPos =
776 TMath::Sqrt(momPosProp[0]*momPosProp[0]+momPosProp[1]*momPosProp[1]+momPosProp[2]*momPosProp[2]);//absolute momentum value of positive daughter
777
778 Double_t scalarproduct =
779 momPosProp[0]*momNegProp[0]+momPosProp[1]*momNegProp[1]+momPosProp[2]*momNegProp[2];//scalar product of propagated positive and negative daughters' momenta
780
781 Double_t chipair = TMath::ACos(scalarproduct/(pEle*pPos));//Angle between propagated daughter tracks
782
783 psiPair = TMath::Abs(TMath::ASin(deltat/chipair));
784
785 return psiPair;
786}
787//___________________________________________________________________
788Bool_t AliESDv0KineCuts::GetConvPosXY(AliESDtrack * const ptrack, AliESDtrack * const ntrack, Double_t convpos[2]){
789 //
790 // recalculate the gamma conversion XY postition
791 //
792
793 const Double_t b = fEvent->GetMagneticField();
794
795 Double_t helixcenterpos[2];
796 GetHelixCenter(ptrack,b,ptrack->Charge(),helixcenterpos);
797
798 Double_t helixcenterneg[2];
799 GetHelixCenter(ntrack,b,ntrack->Charge(),helixcenterneg);
800
801 Double_t poshelix[6];
802 ptrack->GetHelixParameters(poshelix,b);
803 Double_t posradius = TMath::Abs(1./poshelix[4]);
804
805 Double_t neghelix[6];
806 ntrack->GetHelixParameters(neghelix,b);
807 Double_t negradius = TMath::Abs(1./neghelix[4]);
808
809 Double_t xpos = helixcenterpos[0];
810 Double_t ypos = helixcenterpos[1];
811 Double_t xneg = helixcenterneg[0];
812 Double_t yneg = helixcenterneg[1];
813
814 convpos[0] = (xpos*negradius + xneg*posradius)/(negradius+posradius);
815 convpos[1] = (ypos*negradius+ yneg*posradius)/(negradius+posradius);
816
817 return 1;
818}
819//___________________________________________________________________
820Bool_t AliESDv0KineCuts::GetHelixCenter(AliESDtrack * const track, Double_t b,Int_t charge, Double_t center[2]){
821 //
822 // computes the center of the track helix
823 //
824
825 Double_t pi = TMath::Pi();
826
827 Double_t helix[6];
828 track->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
849 if(charge < 0){
850 xpoint = xpoint;
851 ypoint = ypoint;
852 }
853 }
854 if(b>0){
855 if(charge > 0){
856 xpoint = xpoint;
857 ypoint = ypoint;
858 }
859
860 if(charge < 0){
861 xpoint = - xpoint;
862 ypoint = - ypoint;
863 }
864 }
865 center[0] = xpos + xpoint;
866 center[1] = ypos + ypoint;
867
868 return 1;
869}
870//___________________________________________________________________
871AliKFParticle *AliESDv0KineCuts::CreateMotherParticle(const AliVTrack* const pdaughter, const AliVTrack* const ndaughter, Int_t pspec, Int_t nspec){
872 //
873 // Creates a mother particle
874 //
875 AliKFParticle pkfdaughter(*pdaughter, pspec);
876 AliKFParticle nkfdaughter(*ndaughter, nspec);
877
878
879 // Create the mother particle
880 AliKFParticle *m = new AliKFParticle(pkfdaughter, nkfdaughter);
881 // DEBUG - testing
882 if(TMath::Abs(kElectron) == pspec && TMath::Abs(kElectron) == nspec) m->SetMassConstraint(0, 0.001);
883 else if(TMath::Abs(kPiPlus) == pspec && TMath::Abs(kPiPlus) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass(), 0.);
884 else if(TMath::Abs(kProton) == pspec && TMath::Abs(kPiPlus) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(), 0.);
885 else if(TMath::Abs(kPiPlus) == pspec && TMath::Abs(kProton) == nspec) m->SetMassConstraint(TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass(), 0.);
886 else{
887 AliErrorClass("Wrong daughter ID - mass constraint can not be set");
888 }
889
890 AliKFVertex improvedVertex = *fPrimaryVertex;
891 improvedVertex += *m;
892 m->SetProductionVertex(improvedVertex);
893
894 // update 15/06/2010
895 // mother particle will not be added to primary vertex but only to its copy
896 // as this confilcts with calling
897 // m->SetPrimaryVertex() function and
898 // subsequently removing the mother particle afterwards
899 // Source: Sergey Gorbunov
900
901 return m;
902}