1 /**************************************************************************
2 * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 /////////////////////////////////////////////////////////////
18 // Base class for AOD reconstructed decay
20 // Author: A.Dainese, andrea.dainese@lnl.infn.it
21 /////////////////////////////////////////////////////////////
23 #include <TDatabasePDG.h>
25 #include <TClonesArray.h>
27 #include "AliVTrack.h"
28 #include "AliExternalTrackParam.h"
29 #include "AliNeutralTrackParam.h"
30 #include "AliAODMCParticle.h"
31 #include "AliAODRecoDecay.h"
33 ClassImp(AliAODRecoDecay)
35 //--------------------------------------------------------------------------
36 AliAODRecoDecay::AliAODRecoDecay() :
39 fOwnSecondaryVtx(0x0),
41 fNProngs(0), fNDCA(0), fNPID(0),
42 fPx(0x0), fPy(0x0), fPz(0x0),
46 fEventNumber(-1),fRunNumber(-1)
49 // Default Constructor
52 //--------------------------------------------------------------------------
53 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
55 Double_t *px,Double_t *py,Double_t *pz,
59 fOwnSecondaryVtx(0x0),
61 fNProngs(nprongs), fNDCA(0), fNPID(0),
62 fPx(0x0), fPy(0x0), fPz(0x0),
66 fEventNumber(-1),fRunNumber(-1)
69 // Constructor with AliAODVertex for decay vertex
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++) {
83 //--------------------------------------------------------------------------
84 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
89 fOwnSecondaryVtx(0x0),
91 fNProngs(nprongs), fNDCA(0), fNPID(0),
92 fPx(0x0), fPy(0x0), fPz(0x0),
96 fEventNumber(-1),fRunNumber(-1)
99 // Constructor with AliAODVertex for decay vertex and without prongs momenta
102 fd0 = new Double_t[GetNProngs()];
103 for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
105 //--------------------------------------------------------------------------
106 AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
108 fSecondaryVtx(source.fSecondaryVtx),
109 fOwnSecondaryVtx(source.fOwnSecondaryVtx),
110 fCharge(source.fCharge),
111 fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
112 fPx(0x0), fPy(0x0), fPz(0x0),
116 fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
121 if(source.GetNProngs()>0) {
122 fd0 = new Double32_t[GetNProngs()];
123 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
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));
133 fPID = new Double32_t[fNPID];
134 memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
137 fDCA = new Double32_t[fNDCA];
138 memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
142 //--------------------------------------------------------------------------
143 AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
146 // assignment operator
148 if(&source == this) return *this;
149 fSecondaryVtx = source.fSecondaryVtx;
150 fOwnSecondaryVtx = source.fOwnSecondaryVtx;
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) {
158 if(fd0)delete [] fd0;
159 fd0 = new Double32_t[GetNProngs()];
160 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
162 if(fPx) delete [] fPx;
163 fPx = new Double32_t[GetNProngs()];
164 if(fPy) delete [] fPy;
165 fPy = new Double32_t[GetNProngs()];
166 if(fPz) delete [] fPz;
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));
173 if(fPID) delete [] fPID;
174 fPID = new Double32_t[fNPID];
175 memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
178 if(fDCA) delete [] fDCA;
179 fDCA = new Double32_t[fNDCA];
180 memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
185 //--------------------------------------------------------------------------
186 AliAODRecoDecay::~AliAODRecoDecay() {
188 // Default Destructor
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; }
198 //----------------------------------------------------------------------------
199 Double_t AliAODRecoDecay::Alpha() const
202 // Armenteros-Podolanski alpha for 2-prong decays
204 if(GetNProngs()!=2) {
205 printf("Can be called only for 2-prong decays");
206 return (Double_t)-99999.;
208 return 1.-2./(1.+QlProng(0)/QlProng(1));
210 //----------------------------------------------------------------------------
211 Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
214 // Decay length assuming it is produced at "point" [cm]
216 return TMath::Sqrt((point[0]-GetSecVtxX())
217 *(point[0]-GetSecVtxX())
218 +(point[1]-GetSecVtxY())
219 *(point[1]-GetSecVtxY())
220 +(point[2]-GetSecVtxZ())
221 *(point[2]-GetSecVtxZ()));
223 //----------------------------------------------------------------------------
224 Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
227 // Decay length in XY assuming it is produced at "point" [cm]
229 return TMath::Sqrt((point[0]-GetSecVtxX())
230 *(point[0]-GetSecVtxX())
231 +(point[1]-GetSecVtxY())
232 *(point[1]-GetSecVtxY()));
234 //----------------------------------------------------------------------------
235 Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
238 // Cosine of pointing angle in space assuming it is produced at "point"
240 TVector3 mom(Px(),Py(),Pz());
241 TVector3 fline(GetSecVtxX()-point[0],
242 GetSecVtxY()-point[1],
243 GetSecVtxZ()-point[2]);
245 Double_t pta = mom.Angle(fline);
247 return TMath::Cos(pta);
249 //----------------------------------------------------------------------------
250 Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
253 // Cosine of pointing angle in transverse plane assuming it is produced
256 TVector3 momXY(Px(),Py(),0.);
257 TVector3 flineXY(GetSecVtxX()-point[0],
258 GetSecVtxY()-point[1],
261 Double_t ptaXY = momXY.Angle(flineXY);
263 return TMath::Cos(ptaXY);
265 //----------------------------------------------------------------------------
266 Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
269 // Only for 2-prong decays:
270 // Cosine of decay angle (theta*) in the rest frame of the mother particle
271 // for prong ip (0 or 1) with mass hypotheses pdgvtx for mother particle,
272 // pdgprong0 for prong 0 and pdgprong1 for prong1
274 if(GetNProngs()!=2) {
275 printf("Can be called only for 2-prong decays");
276 return (Double_t)-99999.;
278 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
280 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
281 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
283 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);
285 Double_t beta = P()/E(pdgvtx);
286 Double_t gamma = E(pdgvtx)/massvtx;
288 Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
292 //---------------------------------------------------------------------------
293 Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
296 // Decay time * c assuming it is produced at "point" [cm]
298 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
299 return DecayLength(point)*mass/P();
301 //---------------------------------------------------------------------------
302 Double_t AliAODRecoDecay::E(UInt_t pdg) const
307 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
308 return TMath::Sqrt(mass*mass+P()*P());
310 //---------------------------------------------------------------------------
311 Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
314 // Energy of ip-th prong
316 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
317 return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
319 //--------------------------------------------------------------------------
320 Bool_t AliAODRecoDecay::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
322 // This function returns the global covariance matrix of the track params
324 // Cov(x,x) ... : cv[0]
325 // Cov(y,x) ... : cv[1] cv[2]
326 // Cov(z,x) ... : cv[3] cv[4] cv[5]
327 // Cov(px,x)... : cv[6] cv[7] cv[8] cv[9]
328 // Cov(py,x)... : cv[10] cv[11] cv[12] cv[13] cv[14]
329 // Cov(pz,x)... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
331 // For XYZ we take the cov of the vertex, for PxPyPz we take the
332 // sum of the covs of PxPyPz from the daughters, for the moment
333 // we set the cov between position and momentum as the sum of
334 // the same cov from the daughters.
338 for(j=0;j<21;j++) cv[j]=0.;
340 if(!GetNDaughters()) {
341 AliError("No daughters available");
346 AliAODVertex *secv=GetSecondaryVtx();
348 AliError("Vertex covariance matrix not available");
351 if(!secv->GetCovMatrix(v)) {
352 AliError("Vertex covariance matrix not available");
356 Double_t p[21]; for(j=0;j<21;j++) p[j]=0.;
358 for(Int_t i=1; i<GetNDaughters(); i++) {
359 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
361 if(!daugh->GetCovarianceXYZPxPyPz(dcov)) error=kTRUE;
362 for(j=0;j<21;j++) p[j] += dcov[j];
365 AliError("No covariance for at least one daughter")
369 for(j=0; j<21; j++) {
379 //----------------------------------------------------------------------------
380 UChar_t AliAODRecoDecay::GetITSClusterMap() const {
382 // We take the logical AND of the daughters cluster maps
383 // (only if all daughters have the bit for given layer, we set the bit)
387 if(!GetNDaughters()) {
388 AliError("No daughters available");
392 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)
394 for(Int_t i=0; i<GetNDaughters(); i++) {
395 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
396 if(!TESTBIT(daugh->GetITSClusterMap(),l)) bit=0;
398 if(bit) SETBIT(map,l);
403 //--------------------------------------------------------------------------
404 ULong_t AliAODRecoDecay::GetStatus() const {
406 // Same as for ITSClusterMap
410 if(!GetNDaughters()) {
411 AliError("No daughters available");
415 AliVTrack *daugh0 = (AliVTrack*)GetDaughter(0);
416 status = status&(daugh0->GetStatus());
418 for(Int_t i=1; i<GetNDaughters(); i++) {
419 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
420 status = status&(daugh->GetStatus());
425 //--------------------------------------------------------------------------
426 Bool_t AliAODRecoDecay::PropagateToDCA(const AliVVertex* vtx,Double_t b,Double_t maxd,Double_t dz[2],Double_t covar[3])
428 // compute impact parameters to the vertex vtx and their covariance matrix
429 // b is the Bz, needed to propagate correctly the track to vertex
430 // return kFALSE is something went wrong
432 AliWarning("The AliAODRecoDecay momentum is not updated to the DCA");
436 // convert to AliNeutralTrackParam
437 AliNeutralTrackParam ntp(this);
438 retval = ntp.PropagateToDCA(vtx,b,maxd,dz,covar);
440 // convert to AliExternalTrackParam
441 AliExternalTrackParam etp(this);
442 retval = etp.PropagateToDCA(vtx,b,maxd,dz,covar);
446 //--------------------------------------------------------------------------
447 Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
450 // Impact parameter in the bending plane of the particle
453 Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
455 Double_t dx = GetSecVtxX()-point[0]+k*Px();
456 Double_t dy = GetSecVtxY()-point[1]+k*Py();
457 Double_t absImpPar = TMath::Sqrt(dx*dx+dy*dy);
458 TVector3 mom(Px(),Py(),Pz());
459 TVector3 fline(GetSecVtxX()-point[0],
460 GetSecVtxY()-point[1],
461 GetSecVtxZ()-point[2]);
462 TVector3 cross = mom.Cross(fline);
463 return (cross.Z()>0. ? absImpPar : -absImpPar);
465 //----------------------------------------------------------------------------
466 Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
469 // Invariant mass for prongs mass hypotheses in pdg array
471 if(GetNProngs()!=npdg) {
472 printf("npdg != GetNProngs()");
473 return (Double_t)-99999.;
475 Double_t energysum = 0.;
477 for(Int_t i=0; i<GetNProngs(); i++) {
478 energysum += EProng(i,pdg[i]);
481 Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
485 //----------------------------------------------------------------------------
486 Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
487 UInt_t pdg1,UInt_t pdg2) const
490 // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
492 Double_t energysum = EProng(ip1,pdg1) + EProng(ip2,pdg2);
493 Double_t psum2 = (PxProng(ip1)+PxProng(ip2))*(PxProng(ip1)+PxProng(ip2))
494 +(PyProng(ip1)+PyProng(ip2))*(PyProng(ip1)+PyProng(ip2))
495 +(PzProng(ip1)+PzProng(ip2))*(PzProng(ip1)+PzProng(ip2));
496 Double_t mass = TMath::Sqrt(energysum*energysum-psum2);
500 //----------------------------------------------------------------------------
501 Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
502 Int_t ndgCk,Int_t *pdgDg) const
505 // Check if this candidate is matched to a MC signal
507 // If yes, return label (>=0) of the AliAODMCParticle
509 // if ndgCk>0, checks also daughters PDGs
511 Int_t ndg=GetNDaughters();
513 AliError("No daughters available");
517 AliError("Only decays with <10 daughters supported");
520 if(ndgCk>0 && ndgCk!=ndg) {
521 AliError("Wrong number of daughter PDGs passed");
527 // loop on daughters and write the labels
528 for(Int_t i=0; i<ndg; i++) {
529 AliAODTrack *trk = (AliAODTrack*)GetDaughter(i);
530 dgLabels[i] = trk->GetLabel();
533 return MatchToMC(pdgabs,mcArray,dgLabels,ndg,ndgCk,pdgDg);
535 //----------------------------------------------------------------------------
536 Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
537 Int_t dgLabels[10],Int_t ndg,
538 Int_t ndgCk,Int_t *pdgDg) const
541 // Check if this candidate is matched to a MC signal
543 // If yes, return label (>=0) of the AliAODMCParticle
546 Int_t labMom[10]={0,0,0,0,0,0,0,0,0,0};
547 Int_t i,j,lab,labMother,pdgMother,pdgPart;
548 AliAODMCParticle *part=0;
549 AliAODMCParticle *mother=0;
550 Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.;
551 Bool_t pdgUsed[10]={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
553 // loop on daughter labels
554 for(i=0; i<ndg; i++) {
558 printf("daughter with negative label %d\n",lab);
561 part = (AliAODMCParticle*)mcArray->At(lab);
563 printf("no MC particle\n");
567 // check the PDG of the daughter, if requested
569 pdgPart=TMath::Abs(part->GetPdgCode());
570 for(j=0; j<ndg; j++) {
571 if(!pdgUsed[j] && pdgPart==pdgDg[j]) {
578 // for the J/psi, check that the daughters are electrons
580 if(TMath::Abs(part->GetPdgCode())!=11) return -1;
584 while(mother->GetMother()>=0) {
585 labMother=mother->GetMother();
586 mother = (AliAODMCParticle*)mcArray->At(labMother);
588 printf("no MC mother particle\n");
591 pdgMother = TMath::Abs(mother->GetPdgCode());
592 if(pdgMother==pdgabs) {
594 // keep sum of daughters' momenta, to check for mom conservation
595 pxSumDgs += part->Px();
596 pySumDgs += part->Py();
597 pzSumDgs += part->Pz();
599 } else if(pdgMother>pdgabs || pdgMother<10) {
603 if(labMom[i]==-1) return -1; // mother PDG not ok for this daughter
604 } // end loop on daughters
606 // check if the candidate is signal
608 // all labels have to be the same and !=-1
609 for(i=0; i<ndg; i++) {
610 if(labMom[i]==-1) return -1;
611 if(labMom[i]!=labMother) return -1;
614 // check that all daughter PDGs are matched
616 for(i=0; i<ndg; i++) {
617 if(pdgUsed[i]==kFALSE) return -1;
622 // check that the mother decayed in <GetNDaughters()> prongs
623 Int_t ndg2 = TMath::Abs(mother->GetDaughter(1)-mother->GetDaughter(0))+1;
624 printf(" MC daughters %d\n",ndg2);
625 //if(ndg!=GetNDaughters()) return -1;
626 AliAODMCParticle* p1=(AliAODMCParticle*)(mcArray->At(mother->GetDaughter(1)));
627 AliAODMCParticle* p0=(AliAODMCParticle*)(mcArray->At(mother->GetDaughter(0)));
628 printf(">>>>>>>> pdg %d %d %d %d %d %d\n",mother->GetDaughter(0),mother->GetDaughter(1),dgLabels[0],dgLabels[1],p0->GetPdgCode(),p1->GetPdgCode());
631 // the above works only for non-resonant decays,
632 // it's better to check for mom conservation
633 mother = (AliAODMCParticle*)mcArray->At(labMother);
634 Double_t pxMother = mother->Px();
635 Double_t pyMother = mother->Py();
636 Double_t pzMother = mother->Pz();
638 if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.00001 &&
639 (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.00001 &&
640 (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.00001)
645 //---------------------------------------------------------------------------
646 void AliAODRecoDecay::Print(Option_t* /*option*/) const
649 // Print some information
651 printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
652 printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
656 //----------------------------------------------------------------------------
657 Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
660 // Relative angle between two prongs
662 TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
663 TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
665 Double_t angle = momA.Angle(momB);
669 //----------------------------------------------------------------------------
670 Double_t AliAODRecoDecay::QlProng(Int_t ip) const
673 // Longitudinal momentum of prong w.r.t. to total momentum
675 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
676 TVector3 momTot(Px(),Py(),Pz());
678 return mom.Dot(momTot)/momTot.Mag();
680 //----------------------------------------------------------------------------
681 Double_t AliAODRecoDecay::QtProng(Int_t ip) const
684 // Transverse momentum of prong w.r.t. to total momentum
686 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
687 TVector3 momTot(Px(),Py(),Pz());
689 return mom.Perp(momTot);
691 //----------------------------------------------------------------------------
692 Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
695 // Longitudinal momentum of prong w.r.t. to flight line between "point"
698 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
699 TVector3 fline(GetSecVtxX()-point[0],
700 GetSecVtxY()-point[1],
701 GetSecVtxZ()-point[2]);
703 return mom.Dot(fline)/fline.Mag();
705 //----------------------------------------------------------------------------
706 Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
709 // Transverse momentum of prong w.r.t. to flight line between "point" and
712 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
713 TVector3 fline(GetSecVtxX()-point[0],
714 GetSecVtxY()-point[1],
715 GetSecVtxZ()-point[2]);
717 return mom.Perp(fline);
719 //--------------------------------------------------------------------------