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 "AliAODMCParticle.h"
29 #include "AliAODRecoDecay.h"
31 ClassImp(AliAODRecoDecay)
33 //--------------------------------------------------------------------------
34 AliAODRecoDecay::AliAODRecoDecay() :
37 fOwnSecondaryVtx(0x0),
39 fNProngs(0), fNDCA(0), fNPID(0),
40 fPx(0x0), fPy(0x0), fPz(0x0),
44 fEventNumber(-1),fRunNumber(-1)
47 // Default Constructor
50 //--------------------------------------------------------------------------
51 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
53 Double_t *px,Double_t *py,Double_t *pz,
57 fOwnSecondaryVtx(0x0),
59 fNProngs(nprongs), fNDCA(0), fNPID(0),
60 fPx(0x0), fPy(0x0), fPz(0x0),
64 fEventNumber(-1),fRunNumber(-1)
67 // Constructor with AliAODVertex for decay vertex
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++) {
81 //--------------------------------------------------------------------------
82 AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
87 fOwnSecondaryVtx(0x0),
89 fNProngs(nprongs), fNDCA(0), fNPID(0),
90 fPx(0x0), fPy(0x0), fPz(0x0),
94 fEventNumber(-1),fRunNumber(-1)
97 // Constructor with AliAODVertex for decay vertex and without prongs momenta
100 fd0 = new Double_t[GetNProngs()];
101 for(Int_t i=0; i<GetNProngs(); i++) fd0[i] = d0[i];
103 //--------------------------------------------------------------------------
104 AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
106 fSecondaryVtx(source.fSecondaryVtx),
107 fOwnSecondaryVtx(source.fOwnSecondaryVtx),
108 fCharge(source.fCharge),
109 fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
110 fPx(0x0), fPy(0x0), fPz(0x0),
114 fEventNumber(source.fEventNumber),fRunNumber(source.fRunNumber)
119 if(source.GetNProngs()>0) {
120 fd0 = new Double32_t[GetNProngs()];
121 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
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));
131 fPID = new Double32_t[fNPID];
132 memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
135 fDCA = new Double32_t[fNDCA];
136 memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
140 //--------------------------------------------------------------------------
141 AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
144 // assignment operator
146 if(&source == this) return *this;
147 fSecondaryVtx = source.fSecondaryVtx;
148 fOwnSecondaryVtx = source.fOwnSecondaryVtx;
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) {
156 if(fd0)delete [] fd0;
157 fd0 = new Double32_t[GetNProngs()];
158 memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
160 if(fPx) delete [] fPx;
161 fPx = new Double32_t[GetNProngs()];
162 if(fPy) delete [] fPy;
163 fPy = new Double32_t[GetNProngs()];
164 if(fPz) delete [] fPz;
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));
171 if(fPID) delete [] fPID;
172 fPID = new Double32_t[fNPID];
173 memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
176 if(fDCA) delete [] fDCA;
177 fDCA = new Double32_t[fNDCA];
178 memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
183 //--------------------------------------------------------------------------
184 AliAODRecoDecay::~AliAODRecoDecay() {
186 // Default Destructor
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; }
196 //----------------------------------------------------------------------------
197 Double_t AliAODRecoDecay::Alpha() const
200 // Armenteros-Podolanski alpha for 2-prong decays
202 if(GetNProngs()!=2) {
203 printf("Can be called only for 2-prong decays");
204 return (Double_t)-99999.;
206 return 1.-2./(1.+QlProng(0)/QlProng(1));
208 //----------------------------------------------------------------------------
209 Double_t AliAODRecoDecay::DecayLength(Double_t point[3]) const
212 // Decay length assuming it is produced at "point" [cm]
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()));
221 //----------------------------------------------------------------------------
222 Double_t AliAODRecoDecay::DecayLengthXY(Double_t point[3]) const
225 // Decay length in XY assuming it is produced at "point" [cm]
227 return TMath::Sqrt((point[0]-GetSecVtxX())
228 *(point[0]-GetSecVtxX())
229 +(point[1]-GetSecVtxY())
230 *(point[1]-GetSecVtxY()));
232 //----------------------------------------------------------------------------
233 Double_t AliAODRecoDecay::CosPointingAngle(Double_t point[3]) const
236 // Cosine of pointing angle in space assuming it is produced at "point"
238 TVector3 mom(Px(),Py(),Pz());
239 TVector3 fline(GetSecVtxX()-point[0],
240 GetSecVtxY()-point[1],
241 GetSecVtxZ()-point[2]);
243 Double_t pta = mom.Angle(fline);
245 return TMath::Cos(pta);
247 //----------------------------------------------------------------------------
248 Double_t AliAODRecoDecay::CosPointingAngleXY(Double_t point[3]) const
251 // Cosine of pointing angle in transverse plane assuming it is produced
254 TVector3 momXY(Px(),Py(),0.);
255 TVector3 flineXY(GetSecVtxX()-point[0],
256 GetSecVtxY()-point[1],
259 Double_t ptaXY = momXY.Angle(flineXY);
261 return TMath::Cos(ptaXY);
263 //----------------------------------------------------------------------------
264 Double_t AliAODRecoDecay::CosThetaStar(Int_t ip,UInt_t pdgvtx,UInt_t pdgprong0,UInt_t pdgprong1) const
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
272 if(GetNProngs()!=2) {
273 printf("Can be called only for 2-prong decays");
274 return (Double_t)-99999.;
276 Double_t massvtx = TDatabasePDG::Instance()->GetParticle(pdgvtx)->Mass();
278 massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
279 massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
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);
283 Double_t beta = P()/E(pdgvtx);
284 Double_t gamma = E(pdgvtx)/massvtx;
286 Double_t cts = (QlProng(ip)/gamma-beta*TMath::Sqrt(pStar*pStar+massp[ip]*massp[ip]))/pStar;
290 //---------------------------------------------------------------------------
291 Double_t AliAODRecoDecay::Ct(UInt_t pdg,Double_t point[3]) const
294 // Decay time * c assuming it is produced at "point" [cm]
296 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
297 return DecayLength(point)*mass/P();
299 //---------------------------------------------------------------------------
300 Double_t AliAODRecoDecay::E(UInt_t pdg) const
305 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
306 return TMath::Sqrt(mass*mass+P()*P());
308 //---------------------------------------------------------------------------
309 Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) const
312 // Energy of ip-th prong
314 Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
315 return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
317 //--------------------------------------------------------------------------
318 Bool_t AliAODRecoDecay::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
320 // This function returns the global covariance matrix of the track params
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]
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.
336 for(j=0;j<21;j++) cv[j]=0.;
338 if(!GetNDaughters()) {
339 AliError("No daughters available");
344 AliAODVertex *secv=GetSecondaryVtx();
346 AliError("Vertex covariance matrix not available");
349 if(!secv->GetCovMatrix(v)) {
350 AliError("Vertex covariance matrix not available");
354 Double_t p[21]; for(j=0;j<21;j++) p[j]=0.;
356 for(Int_t i=1; i<GetNDaughters(); i++) {
357 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
359 if(!daugh->GetCovarianceXYZPxPyPz(dcov)) error=kTRUE;
360 for(j=0;j<21;j++) p[j] += dcov[j];
363 AliError("No covariance for at least one daughter")
367 for(j=0; j<21; j++) {
377 //----------------------------------------------------------------------------
378 UChar_t AliAODRecoDecay::GetITSClusterMap() const {
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)
385 if(!GetNDaughters()) {
386 AliError("No daughters available");
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)
392 for(Int_t i=0; i<GetNDaughters(); i++) {
393 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
394 if(!TESTBIT(daugh->GetITSClusterMap(),l)) bit=0;
396 if(bit) SETBIT(map,l);
401 //--------------------------------------------------------------------------
402 ULong_t AliAODRecoDecay::GetStatus() const {
404 // Same as for ITSClusterMap
408 if(!GetNDaughters()) {
409 AliError("No daughters available");
413 AliVTrack *daugh0 = (AliVTrack*)GetDaughter(0);
414 status = status&(daugh0->GetStatus());
416 for(Int_t i=1; i<GetNDaughters(); i++) {
417 AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
418 status = status&(daugh->GetStatus());
423 //--------------------------------------------------------------------------
424 Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
427 // Impact parameter in the bending plane of the particle
430 Double_t k = -(GetSecVtxX()-point[0])*Px()-(GetSecVtxY()-point[1])*Py();
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);
442 //----------------------------------------------------------------------------
443 Double_t AliAODRecoDecay::InvMass(Int_t npdg,UInt_t *pdg) const
446 // Invariant mass for prongs mass hypotheses in pdg array
448 if(GetNProngs()!=npdg) {
449 printf("npdg != GetNProngs()");
450 return (Double_t)-99999.;
452 Double_t energysum = 0.;
454 for(Int_t i=0; i<GetNProngs(); i++) {
455 energysum += EProng(i,pdg[i]);
458 Double_t mass = TMath::Sqrt(energysum*energysum-P()*P());
462 //----------------------------------------------------------------------------
463 Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
464 UInt_t pdg1,UInt_t pdg2) const
467 // 2-prong(ip1,ip2) invariant mass for prongs mass hypotheses in pdg1,2
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);
477 //----------------------------------------------------------------------------
478 Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
479 Int_t ndgCk,Int_t *pdgDg) const
482 // Check if this candidate is matched to a MC signal
484 // If yes, return label (>=0) of the AliAODMCParticle
486 // if ndgCk>0, checks also daughters PDGs
488 Int_t ndg=GetNDaughters();
490 AliError("No daughters available");
494 AliError("Only decays with <10 daughters supported");
497 if(ndgCk>0 && ndgCk!=ndg) {
498 AliError("Wrong number of daughter PDGs passed");
504 // loop on daughters and write the labels
505 for(Int_t i=0; i<ndg; i++) {
506 AliAODTrack *trk = (AliAODTrack*)GetDaughter(i);
507 dgLabels[i] = trk->GetLabel();
510 return MatchToMC(pdgabs,mcArray,dgLabels,ndg,ndgCk,pdgDg);
512 //----------------------------------------------------------------------------
513 Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
514 Int_t dgLabels[10],Int_t ndg,
515 Int_t ndgCk,Int_t *pdgDg) const
518 // Check if this candidate is matched to a MC signal
520 // If yes, return label (>=0) of the AliAODMCParticle
523 Int_t labMom[10]={0,0,0,0,0,0,0,0,0,0};
524 Int_t i,j,lab,labMother,pdgMother,pdgPart;
525 AliAODMCParticle *part=0;
526 AliAODMCParticle *mother=0;
527 Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.;
528 Bool_t pdgUsed[10]={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
530 // loop on daughter labels
531 for(i=0; i<ndg; i++) {
535 printf("daughter with negative label %d\n",lab);
538 part = (AliAODMCParticle*)mcArray->At(lab);
540 printf("no MC particle\n");
544 // check the PDG of the daughter, if requested
546 pdgPart=TMath::Abs(part->GetPdgCode());
547 for(j=0; j<ndg; j++) {
548 if(!pdgUsed[j] && pdgPart==pdgDg[j]) {
555 // for the J/psi, check that the daughters are electrons
557 if(TMath::Abs(part->GetPdgCode())!=11) return -1;
561 while(mother->GetMother()>=0) {
562 labMother=mother->GetMother();
563 mother = (AliAODMCParticle*)mcArray->At(labMother);
565 printf("no MC mother particle\n");
568 pdgMother = TMath::Abs(mother->GetPdgCode());
569 if(pdgMother==pdgabs) {
571 // keep sum of daughters' momenta, to check for mom conservation
572 pxSumDgs += part->Px();
573 pySumDgs += part->Py();
574 pzSumDgs += part->Pz();
576 } else if(pdgMother>pdgabs || pdgMother<10) {
580 if(labMom[i]==-1) return -1; // mother PDG not ok for this daughter
581 } // end loop on daughters
583 // check if the candidate is signal
585 // all labels have to be the same and !=-1
586 for(i=0; i<ndg; i++) {
587 if(labMom[i]==-1) return -1;
588 if(labMom[i]!=labMother) return -1;
591 // check that all daughter PDGs are matched
593 for(i=0; i<ndg; i++) {
594 if(pdgUsed[i]==kFALSE) return -1;
599 // check that the mother decayed in <GetNDaughters()> prongs
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());
608 // the above works only for non-resonant decays,
609 // it's better to check for mom conservation
610 mother = (AliAODMCParticle*)mcArray->At(labMother);
611 Double_t pxMother = mother->Px();
612 Double_t pyMother = mother->Py();
613 Double_t pzMother = mother->Pz();
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)
622 //---------------------------------------------------------------------------
623 void AliAODRecoDecay::Print(Option_t* /*option*/) const
626 // Print some information
628 printf("AliAODRecoDecay with %d prongs\n",GetNProngs());
629 printf("Secondary Vertex: (%f, %f, %f)\n",GetSecVtxX(),GetSecVtxY(),GetSecVtxZ());
633 //----------------------------------------------------------------------------
634 Double_t AliAODRecoDecay::ProngsRelAngle(Int_t ip1,Int_t ip2) const
637 // Relative angle between two prongs
639 TVector3 momA(PxProng(ip1),PyProng(ip1),PzProng(ip1));
640 TVector3 momB(PxProng(ip2),PyProng(ip2),PzProng(ip2));
642 Double_t angle = momA.Angle(momB);
646 //----------------------------------------------------------------------------
647 Double_t AliAODRecoDecay::QlProng(Int_t ip) const
650 // Longitudinal momentum of prong w.r.t. to total momentum
652 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
653 TVector3 momTot(Px(),Py(),Pz());
655 return mom.Dot(momTot)/momTot.Mag();
657 //----------------------------------------------------------------------------
658 Double_t AliAODRecoDecay::QtProng(Int_t ip) const
661 // Transverse momentum of prong w.r.t. to total momentum
663 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
664 TVector3 momTot(Px(),Py(),Pz());
666 return mom.Perp(momTot);
668 //----------------------------------------------------------------------------
669 Double_t AliAODRecoDecay::QlProngFlightLine(Int_t ip,Double_t point[3]) const
672 // Longitudinal momentum of prong w.r.t. to flight line between "point"
675 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
676 TVector3 fline(GetSecVtxX()-point[0],
677 GetSecVtxY()-point[1],
678 GetSecVtxZ()-point[2]);
680 return mom.Dot(fline)/fline.Mag();
682 //----------------------------------------------------------------------------
683 Double_t AliAODRecoDecay::QtProngFlightLine(Int_t ip,Double_t point[3]) const
686 // Transverse momentum of prong w.r.t. to flight line between "point" and
689 TVector3 mom(PxProng(ip),PyProng(ip),PzProng(ip));
690 TVector3 fline(GetSecVtxX()-point[0],
691 GetSecVtxY()-point[1],
692 GetSecVtxZ()-point[2]);
694 return mom.Perp(fline);
696 //--------------------------------------------------------------------------