]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliAODRecoDecay.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
The description of changes:
[u/mrichter/AliRoot.git] / STEER / AliAODRecoDecay.cxx
CommitLineData
7de7497b 1/**************************************************************************
2 * Copyright(c) 1998-2006, 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//
18// Base class for AOD reconstructed decay
19//
20// Author: A.Dainese, andrea.dainese@lnl.infn.it
21/////////////////////////////////////////////////////////////
22
23#include <TDatabasePDG.h>
24#include <TVector3.h>
bfd68286 25#include <TClonesArray.h>
e045ffda 26#include "AliLog.h"
27#include "AliVTrack.h"
bfd68286 28#include "AliAODMCParticle.h"
7de7497b 29#include "AliAODRecoDecay.h"
30
31ClassImp(AliAODRecoDecay)
32
33//--------------------------------------------------------------------------
34AliAODRecoDecay::AliAODRecoDecay() :
e045ffda 35 AliVTrack(),
7de7497b 36 fSecondaryVtx(0x0),
a11de4a0 37 fOwnSecondaryVtx(0x0),
7de7497b 38 fCharge(0),
39 fNProngs(0), fNDCA(0), fNPID(0),
40 fPx(0x0), fPy(0x0), fPz(0x0),
41 fd0(0x0),
42 fDCA(0x0),
43 fPID(0x0),
44 fEventNumber(-1),fRunNumber(-1)
45{
46 //
47 // Default Constructor
48 //
7de7497b 49}
50//--------------------------------------------------------------------------
51AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
52 Short_t charge,
53 Double_t *px,Double_t *py,Double_t *pz,
54 Double_t *d0) :
e045ffda 55 AliVTrack(),
7de7497b 56 fSecondaryVtx(vtx2),
a11de4a0 57 fOwnSecondaryVtx(0x0),
7de7497b 58 fCharge(charge),
59 fNProngs(nprongs), fNDCA(0), fNPID(0),
60 fPx(0x0), fPy(0x0), fPz(0x0),
61 fd0(0x0),
62 fDCA(0x0),
63 fPID(0x0),
64 fEventNumber(-1),fRunNumber(-1)
65{
66 //
67 // Constructor with AliAODVertex for decay vertex
68 //
69
70 fPx = new Double_t[GetNProngs()];
71 fPy = new Double_t[GetNProngs()];
72 fPz = new Double_t[GetNProngs()];
73 fd0 = new Double_t[GetNProngs()];
74 for(Int_t i=0; i<GetNProngs(); i++) {
75 fPx[i] = px[i];
76 fPy[i] = py[i];
77 fPz[i] = pz[i];
78 fd0[i] = d0[i];
79 }
80}
81//--------------------------------------------------------------------------
82AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
83 Short_t charge,
84 Double_t *d0) :
e045ffda 85 AliVTrack(),
7de7497b 86 fSecondaryVtx(vtx2),
a11de4a0 87 fOwnSecondaryVtx(0x0),
7de7497b 88 fCharge(charge),
89 fNProngs(nprongs), fNDCA(0), fNPID(0),
90 fPx(0x0), fPy(0x0), fPz(0x0),
91 fd0(0x0),
92 fDCA(0x0),
93 fPID(0x0),
94 fEventNumber(-1),fRunNumber(-1)
95{
96 //
97 // Constructor with AliAODVertex for decay vertex and without prongs momenta
98 //
99
100 fd0 = new Double_t[GetNProngs()];
101 for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
102}
103//--------------------------------------------------------------------------
104AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
e045ffda 105 AliVTrack(source),
7de7497b 106 fSecondaryVtx(source.fSecondaryVtx),
a11de4a0 107 fOwnSecondaryVtx(source.fOwnSecondaryVtx),
7de7497b 108 fCharge(source.fCharge),
109 fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
110 fPx(0x0), fPy(0x0), fPz(0x0),
111 fd0(0x0),
112 fDCA(0x0),
113 fPID(0x0),
114 fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
115{
116 //
117 // Copy constructor
118 //
119 if(source.GetNProngs()>0) {
58b0186f 120 fd0 = new Double32_t[GetNProngs()];
121 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
7de7497b 122 if(source.fPx) {
58b0186f 123 fPx = new Double32_t[GetNProngs()];
124 fPy = new Double32_t[GetNProngs()];
125 fPz = new Double32_t[GetNProngs()];
126 memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
127 memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
128 memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
7de7497b 129 }
130 if(source.fPID) {
58b0186f 131 fPID = new Double32_t[fNPID];
132 memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
7de7497b 133 }
134 if(source.fDCA) {
58b0186f 135 fDCA = new Double32_t[fNDCA];
136 memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
7de7497b 137 }
138 }
139}
140//--------------------------------------------------------------------------
141AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
142{
143 //
144 // assignment operator
145 //
146 if(&source == this) return *this;
147 fSecondaryVtx = source.fSecondaryVtx;
a11de4a0 148 fOwnSecondaryVtx = source.fOwnSecondaryVtx;
7de7497b 149 fCharge = source.fCharge;
150 fNProngs = source.fNProngs;
151 fNDCA = source.fNDCA;
152 fNPID = source.fNPID;
153 fEventNumber = source.fEventNumber;
154 fRunNumber = source.fRunNumber;
155 if(source.GetNProngs()>0) {
341952b8 156 if(fd0)delete [] fd0;
58b0186f 157 fd0 = new Double32_t[GetNProngs()];
158 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
7de7497b 159 if(source.fPx) {
341952b8 160 if(fPx) delete [] fPx;
58b0186f 161 fPx = new Double32_t[GetNProngs()];
341952b8 162 if(fPy) delete [] fPy;
58b0186f 163 fPy = new Double32_t[GetNProngs()];
341952b8 164 if(fPz) delete [] fPz;
58b0186f 165 fPz = new Double32_t[GetNProngs()];
166 memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
167 memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
168 memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
7de7497b 169 }
170 if(source.fPID) {
341952b8 171 if(fPID) delete [] fPID;
58b0186f 172 fPID = new Double32_t[fNPID];
173 memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
7de7497b 174 }
175 if(source.fDCA) {
341952b8 176 if(fDCA) delete [] fDCA;
58b0186f 177 fDCA = new Double32_t[fNDCA];
178 memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
7de7497b 179 }
180 }
181 return *this;
182}
183//--------------------------------------------------------------------------
184AliAODRecoDecay::~AliAODRecoDecay() {
185 //
186 // Default Destructor
187 //
a11de4a0 188 if(fPx) { delete [] fPx; fPx=NULL; }
189 if(fPy) { delete [] fPy; fPy=NULL; }
190 if(fPz) { delete [] fPz; fPz=NULL; }
191 if(fd0) { delete [] fd0; fd0=NULL; }
192 if(fPID) { delete [] fPID; fPID=NULL; }
193 if(fDCA) { delete [] fDCA; fDCA=NULL; }
194 if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
7de7497b 195}
196//----------------------------------------------------------------------------
197Double_t AliAODRecoDecay::Alpha() const
198{
199 //
200 // Armenteros-Podolanski alpha for 2-prong decays
201 //
202 if(GetNProngs()!=2) {
203 printf("Can be called only for 2-prong decays");
204 return (Double_t)-99999.;
205 }
206 return 1.-2./(1.+QlProng(0)/QlProng(1));
207}
208//----------------------------------------------------------------------------
209Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
210{
211 //
212 // Decay length assuming it is produced at "point" [cm]
213 //
214 return TMath::Sqrt((point[0]-GetSecVtxX())
215 *(point[0]-GetSecVtxX())
216 +(point[1]-GetSecVtxY())
217 *(point[1]-GetSecVtxY())
218 +(point[2]-GetSecVtxZ())
219 *(point[2]-GetSecVtxZ()));
220}
221//----------------------------------------------------------------------------
222Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
223{
224 //
225 // Decay length in XY assuming it is produced at "point" [cm]
226 //
227 return TMath::Sqrt((point[0]-GetSecVtxX())
228 *(point[0]-GetSecVtxX())
229 +(point[1]-GetSecVtxY())
230 *(point[1]-GetSecVtxY()));
231}
232//----------------------------------------------------------------------------
233Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
234{
235 //
236 // Cosine of pointing angle in space assuming it is produced at "point"
237 //
238 TVector3 mom(Px(),Py(),Pz());
239 TVector3 fline(GetSecVtxX()-point[0],
240 GetSecVtxY()-point[1],
241 GetSecVtxZ()-point[2]);
242
243 Double_t pta = mom.Angle(fline);
244
245 return TMath::Cos(pta);
246}
247//----------------------------------------------------------------------------
248Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
249{
250 //
251 // Cosine of pointing angle in transverse plane assuming it is produced
252 // at "point"
253 //
254 TVector3 momXY(Px(),Py(),0.);
255 TVector3 flineXY(GetSecVtxX()-point[0],
256 GetSecVtxY()-point[1],
257 0.);
258
259 Double_t ptaXY = momXY.Angle(flineXY);
260
261 return TMath::Cos(ptaXY);
262}
263//----------------------------------------------------------------------------
264Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
265{
266 //
267 // Only for 2-prong decays:
268 // Cosine of decay angle (theta*) in the rest frame of the mother particle
269 // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
270 // pdgprong0 for prong 0 and pdgprong1 for prong1
271 //
272 if(GetNProngs()!=2) {
273 printf("Can be called only for 2-prong decays");
274 return (Double_t)-99999.;
275 }
276 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
277 Double_t massp[2];
278 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
279 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
280
281 Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
282
283 Double_t beta = P()/E(pdgvtx);
284 Double_t gamma = E(pdgvtx)/massvtx;
285
286 Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
287
288 return cts;
289}
290//---------------------------------------------------------------------------
291Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
292{
293 //
294 // Decay time * c assuming it is produced at "point" [cm]
295 //
296 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
297 return DecayLength(point)*mass/P();
298}
299//---------------------------------------------------------------------------
300Double_t AliAODRecoDecay::E(UInt_t pdg) const
301{
302 //
303 // Energy
304 //
305 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
306 return TMath::Sqrt(mass*mass+P()*P());
307}
308//---------------------------------------------------------------------------
309Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
310{
311 //
312 // Energy of ip-th prong
313 //
314 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
315 return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
316}
e045ffda 317//--------------------------------------------------------------------------
318Bool_t AliAODRecoDecay::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
319 //
320 // This function returns the global covariance matrix of the track params
321 //
322 // Cov(x,x) ... : cv[0]
323 // Cov(y,x) ... : cv[1] cv[2]
324 // Cov(z,x) ... : cv[3] cv[4] cv[5]
325 // Cov(px,x)... : cv[6] cv[7] cv[8] cv[9]
326 // Cov(py,x)... : cv[10] cv[11] cv[12] cv[13] cv[14]
327 // Cov(pz,x)... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
7de7497b 328 //
e045ffda 329 // For XYZ we take the cov of the vertex, for PxPyPz we take the
330 // sum of the covs of PxPyPz from the daughters, for the moment
331 // we set the cov between position and momentum as the sum of
332 // the same cov from the daughters.
7de7497b 333 //
e045ffda 334
335 Int_t j;
336 for(j=0;j<21;j++) cv[j]=0.;
337
338 if(!GetNDaughters()) {
339 AliError("No daughters available");
340 return kFALSE;
341 }
342
343 Double_t v[6];
344 AliAODVertex *secv=GetSecondaryVtx();
345 if(!secv) {
346 AliError("Vertex covariance matrix not available");
347 return kFALSE;
348 }
349 if(!secv->GetCovMatrix(v)) {
350 AliError("Vertex covariance matrix not available");
351 return kFALSE;
352 }
353
354 Double_t p[21]; for(j=0;j<21;j++) p[j]=0.;
355 Bool_t error=kFALSE;
356 for(Int_t i=1; i<GetNDaughters(); i++) {
357 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
358 Double_t dcov[21];
359 if(!daugh->GetCovarianceXYZPxPyPz(dcov)) error=kTRUE;
360 for(j=0;j<21;j++) p[j] += dcov[j];
361 }
362 if(error) {
363 AliError("No covariance for at least one daughter")
364 return kFALSE;
365 }
366
367 for(j=0; j<21; j++) {
368 if(j<6) {
369 cv[j] = v[j];
370 } else {
371 cv[j] = p[j];
372 }
373 }
374
375 return kTRUE;
376}
7de7497b 377//----------------------------------------------------------------------------
e045ffda 378UChar_t AliAODRecoDecay::GetITSClusterMap() const {
379 //
380 // We take the logical AND of the daughters cluster maps
381 // (only if all daughters have the bit for given layer, we set the bit)
382 //
383 UChar_t map=0;
384
385 if(!GetNDaughters()) {
386 AliError("No daughters available");
387 return map;
388 }
389
390 for(Int_t l=0; l<12; l++) { // loop on ITS layers (from here we cannot know how many they are; let's put 12 to be conservative)
391 Int_t bit = 1;
392 for(Int_t i=0; i<GetNDaughters(); i++) {
393 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
394 if(!TESTBIT(daugh->GetITSClusterMap(),l)) bit=0;
395 }
396 if(bit) SETBIT(map,l);
397 }
398
399 return map;
400}
401//--------------------------------------------------------------------------
402ULong_t AliAODRecoDecay::GetStatus() const {
403 //
404 // Same as for ITSClusterMap
405 //
406 ULong_t status=0;
407
408 if(!GetNDaughters()) {
409 AliError("No daughters available");
410 return status;
411 }
412
413 AliVTrack *daugh0 = (AliVTrack*)GetDaughter(0);
414 status = status&(daugh0->GetStatus());
415
416 for(Int_t i=1; i<GetNDaughters(); i++) {
417 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
418 status = status&(daugh->GetStatus());
419 }
420
421 return status;
422}
423//--------------------------------------------------------------------------
7de7497b 424Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
425{
426 //
427 // Impact parameter in the bending plane of the particle
428 // w.r.t. to "point"
429 //
430 Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
431 k /= Pt()*Pt();
432 Double_t dx = GetSecVtxX()-point[0]+k*Px();
433 Double_t dy = GetSecVtxY()-point[1]+k*Py();
434 Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
435 TVector3 mom(Px(),Py(),Pz());
436 TVector3 fline(GetSecVtxX()-point[0],
437 GetSecVtxY()-point[1],
438 GetSecVtxZ()-point[2]);
439 TVector3 cross = mom.Cross(fline);
440 return (cross.Z()>0. ? absImpPar : -absImpPar);
441}
442//----------------------------------------------------------------------------
443Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
444{
445 //
446 // Invariant mass for prongs mass hypotheses in pdg array
447 //
448 if(GetNProngs()!=npdg) {
449 printf("npdg != GetNProngs()");
450 return (Double_t)-99999.;
451 }
452 Double_t energysum = 0.;
453
454 for(Int_t i=0; i<GetNProngs(); i++) {
455 energysum += EProng(i,pdg[i]);
456 }
457
458 Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
459
460 return mass;
461}
462//----------------------------------------------------------------------------
463Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
464 UInt_t pdg1,UInt_t pdg2) const
465{
466 //
467 // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
468 //
469 Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
470 Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
471 +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
472 +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
473 Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
474
475 return mass;
476}
bfd68286 477//----------------------------------------------------------------------------
9c4d0454 478Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
479 Int_t ndgCk,Int_t *pdgDg) const
bfd68286 480{
481 //
482 // Check if this candidate is matched to a MC signal
483 // If no, return -1
484 // If yes, return label (>=0) of the AliAODMCParticle
485 //
9c4d0454 486 // if ndgCk>0, checks also daughters PDGs
487 //
8862e6bf 488 Int_t ndg=GetNDaughters();
489 if(!ndg) {
bfd68286 490 AliError("No daughters available");
491 return -1;
492 }
8862e6bf 493 if(ndg>10) {
494 AliError("Only decays with <10 daughters supported");
495 return -1;
496 }
9c4d0454 497 if(ndgCk>0 && ndgCk!=ndg) {
498 AliError("Wrong number of daughter PDGs passed");
499 return -1;
500 }
bfd68286 501
8862e6bf 502 Int_t dgLabels[10];
bfd68286 503
504 // loop on daughters and write the labels
8862e6bf 505 for(Int_t i=0; i<ndg; i++) {
506 AliAODTrack *trk = (AliAODTrack*)GetDaughter(i);
bfd68286 507 dgLabels[i] = trk->GetLabel();
508 }
509
9c4d0454 510 return MatchToMC(pdgabs,mcArray,dgLabels,ndg,ndgCk,pdgDg);
bfd68286 511}
512//----------------------------------------------------------------------------
513Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
9c4d0454 514 Int_t dgLabels[10],Int_t ndg,
515 Int_t ndgCk,Int_t *pdgDg) const
bfd68286 516{
517 //
518 // Check if this candidate is matched to a MC signal
519 // If no, return -1
520 // If yes, return label (>=0) of the AliAODMCParticle
521 //
522
0f0c06de 523 Int_t labMom[10]={0,0,0,0,0,0,0,0,0,0};
9c4d0454 524 Int_t i,j,lab,labMother,pdgMother,pdgPart;
bfd68286 525 AliAODMCParticle *part=0;
8862e6bf 526 AliAODMCParticle *mother=0;
54ea2d63 527 Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.;
9c4d0454 528 Bool_t pdgUsed[10]={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
bfd68286 529
530 // loop on daughter labels
8862e6bf 531 for(i=0; i<ndg; i++) {
bfd68286 532 labMom[i]=-1;
533 lab = dgLabels[i];
534 if(lab<0) {
8862e6bf 535 printf("daughter with negative label %d\n",lab);
536 return -1;
bfd68286 537 }
538 part = (AliAODMCParticle*)mcArray->At(lab);
539 if(!part) {
540 printf("no MC particle\n");
8862e6bf 541 return -1;
bfd68286 542 }
8862e6bf 543
9c4d0454 544 // check the PDG of the daughter, if requested
545 if(ndgCk>0) {
546 pdgPart=TMath::Abs(part->GetPdgCode());
547 for(j=0; j<ndg; j++) {
548 if(!pdgUsed[j] && pdgPart==pdgDg[j]) {
549 pdgUsed[j]=kTRUE;
550 break;
551 }
552 }
553 }
554
54ea2d63 555 // for the J/psi, check that the daughters are electrons
8862e6bf 556 if(pdgabs==443) {
557 if(TMath::Abs(part->GetPdgCode())!=11) return -1;
558 }
54ea2d63 559
8862e6bf 560 mother = part;
561 while(mother->GetMother()>=0) {
562 labMother=mother->GetMother();
563 mother = (AliAODMCParticle*)mcArray->At(labMother);
564 if(!mother) {
bfd68286 565 printf("no MC mother particle\n");
566 break;
567 }
8862e6bf 568 pdgMother = TMath::Abs(mother->GetPdgCode());
bfd68286 569 if(pdgMother==pdgabs) {
570 labMom[i]=labMother;
8862e6bf 571 // keep sum of daughters' momenta, to check for mom conservation
572 pxSumDgs += part->Px();
573 pySumDgs += part->Py();
574 pzSumDgs += part->Pz();
575 break;
576 } else if(pdgMother>pdgabs || pdgMother<10) {
bfd68286 577 break;
578 }
579 }
8862e6bf 580 if(labMom[i]==-1) return -1; // mother PDG not ok for this daughter
581 } // end loop on daughters
bfd68286 582
583 // check if the candidate is signal
bfd68286 584 labMother=labMom[0];
585 // all labels have to be the same and !=-1
8862e6bf 586 for(i=0; i<ndg; i++) {
587 if(labMom[i]==-1) return -1;
588 if(labMom[i]!=labMother) return -1;
bfd68286 589 }
590
9c4d0454 591 // check that all daughter PDGs are matched
592 if(ndgCk>0) {
593 for(i=0; i<ndg; i++) {
594 if(pdgUsed[i]==kFALSE) return -1;
595 }
596 }
8862e6bf 597
598 /*
bfd68286 599 // check that the mother decayed in <GetNDaughters()> prongs
8862e6bf 600 Int_t ndg2 = TMath::Abs(mother->GetDaughter(1)-mother->GetDaughter(0))+1;
601 printf(" MC daughters %d\n",ndg2);
602 //if(ndg!=GetNDaughters()) return -1;
603 AliAODMCParticle* p1=(AliAODMCParticle*)(mcArray->At(mother->GetDaughter(1)));
604 AliAODMCParticle* p0=(AliAODMCParticle*)(mcArray->At(mother->GetDaughter(0)));
605 printf(">>>>>>>> pdg %d %d %d %d %d %d\n",mother->GetDaughter(0),mother->GetDaughter(1),dgLabels[0],dgLabels[1],p0->GetPdgCode(),p1->GetPdgCode());
54ea2d63 606 */
bfd68286 607
54ea2d63 608 // the above works only for non-resonant decays,
609 // it's better to check for mom conservation
8862e6bf 610 mother = (AliAODMCParticle*)mcArray->At(labMother);
611 Double_t pxMother = mother->Px();
612 Double_t pyMother = mother->Py();
613 Double_t pzMother = mother->Pz();
a364c44a 614 // within 0.1%
9c4d0454 615 if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.00001 &&
616 (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.00001 &&
617 (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.00001)
a364c44a 618 return -1;
54ea2d63 619
bfd68286 620 return labMother;
621}
7de7497b 622//---------------------------------------------------------------------------
623void AliAODRecoDecay::Print(Option_t* /*option*/) const
624{
625 //
626 // Print some information
627 //
628 printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
629 printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
630
631 return;
632}
633//----------------------------------------------------------------------------
634Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
635{
636 //
637 // Relative angle between two prongs
638 //
639 TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
640 TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
641
642 Double_t angle = momA.Angle(momB);
643
644 return angle;
645}
646//----------------------------------------------------------------------------
647Double_t AliAODRecoDecay::QlProng(Int_t ip) const
648{
649 //
650 // Longitudinal momentum of prong w.r.t. to total momentum
651 //
652 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
653 TVector3 momTot(Px(),Py(),Pz());
654
655 return mom.Dot(momTot)/momTot.Mag();
656}
657//----------------------------------------------------------------------------
658Double_t AliAODRecoDecay::QtProng(Int_t ip) const
659{
660 //
661 // Transverse momentum of prong w.r.t. to total momentum
662 //
663 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
664 TVector3 momTot(Px(),Py(),Pz());
665
666 return mom.Perp(momTot);
667}
668//----------------------------------------------------------------------------
669Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
670{
671 //
672 // Longitudinal momentum of prong w.r.t. to flight line between "point"
673 // and fSecondaryVtx
674 //
675 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
676 TVector3 fline(GetSecVtxX()-point[0],
677 GetSecVtxY()-point[1],
678 GetSecVtxZ()-point[2]);
679
680 return mom.Dot(fline)/fline.Mag();
681}
682//----------------------------------------------------------------------------
683Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
684{
685 //
686 // Transverse momentum of prong w.r.t. to flight line between "point" and
687 // fSecondaryVtx
688 //
689 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
690 TVector3 fline(GetSecVtxX()-point[0],
691 GetSecVtxY()-point[1],
692 GetSecVtxZ()-point[2]);
693
694 return mom.Perp(fline);
695}
696//--------------------------------------------------------------------------