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