#include <TDatabasePDG.h>
#include <TVector3.h>
-#include "AliVParticle.h"
+#include <TClonesArray.h>
+#include "AliLog.h"
+#include "AliVTrack.h"
+#include "AliAODMCParticle.h"
#include "AliAODRecoDecay.h"
ClassImp(AliAODRecoDecay)
//--------------------------------------------------------------------------
AliAODRecoDecay::AliAODRecoDecay() :
- AliVParticle(),
+ AliVTrack(),
fSecondaryVtx(0x0),
+ fOwnSecondaryVtx(0x0),
fCharge(0),
fNProngs(0), fNDCA(0), fNPID(0),
fPx(0x0), fPy(0x0), fPz(0x0),
//
// Default Constructor
//
- fSecondaryVtx = new AliAODVertex();
}
//--------------------------------------------------------------------------
AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
Short_t charge,
Double_t *px,Double_t *py,Double_t *pz,
Double_t *d0) :
- AliVParticle(),
+ AliVTrack(),
fSecondaryVtx(vtx2),
+ fOwnSecondaryVtx(0x0),
fCharge(charge),
fNProngs(nprongs), fNDCA(0), fNPID(0),
fPx(0x0), fPy(0x0), fPz(0x0),
AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs,
Short_t charge,
Double_t *d0) :
- AliVParticle(),
+ AliVTrack(),
fSecondaryVtx(vtx2),
+ fOwnSecondaryVtx(0x0),
fCharge(charge),
fNProngs(nprongs), fNDCA(0), fNPID(0),
fPx(0x0), fPy(0x0), fPz(0x0),
}
//--------------------------------------------------------------------------
AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) :
- AliVParticle(source),
+ AliVTrack(source),
fSecondaryVtx(source.fSecondaryVtx),
+ fOwnSecondaryVtx(source.fOwnSecondaryVtx),
fCharge(source.fCharge),
fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
fPx(0x0), fPy(0x0), fPz(0x0),
// Copy constructor
//
if(source.GetNProngs()>0) {
- fd0 = new Double_t[GetNProngs()];
- memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double_t));
+ fd0 = new Double32_t[GetNProngs()];
+ memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
if(source.fPx) {
- fPx = new Double_t[GetNProngs()];
- fPy = new Double_t[GetNProngs()];
- fPz = new Double_t[GetNProngs()];
- memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double_t));
- memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double_t));
- memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double_t));
+ fPx = new Double32_t[GetNProngs()];
+ fPy = new Double32_t[GetNProngs()];
+ fPz = new Double32_t[GetNProngs()];
+ memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
+ memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
+ memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
}
if(source.fPID) {
- fPID = new Double_t[5*GetNProngs()];
- memcpy(fPID,source.fPID,GetNProngs()*sizeof(Double_t));
+ fPID = new Double32_t[fNPID];
+ memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
}
if(source.fDCA) {
- fDCA = new Float_t[GetNProngs()*(GetNProngs()-1)/2];
- memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Float_t));
+ fDCA = new Double32_t[fNDCA];
+ memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
}
}
}
//
if(&source == this) return *this;
fSecondaryVtx = source.fSecondaryVtx;
+ fOwnSecondaryVtx = source.fOwnSecondaryVtx;
fCharge = source.fCharge;
fNProngs = source.fNProngs;
fNDCA = source.fNDCA;
fEventNumber = source.fEventNumber;
fRunNumber = source.fRunNumber;
if(source.GetNProngs()>0) {
- fd0 = new Double_t[GetNProngs()];
- memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double_t));
+ if(fd0)delete [] fd0;
+ fd0 = new Double32_t[GetNProngs()];
+ memcpy(fd0,source.fd0,GetNProngs()*sizeof(Double32_t));
if(source.fPx) {
- fPx = new Double_t[GetNProngs()];
- fPy = new Double_t[GetNProngs()];
- fPz = new Double_t[GetNProngs()];
- memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double_t));
- memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double_t));
- memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double_t));
+ if(fPx) delete [] fPx;
+ fPx = new Double32_t[GetNProngs()];
+ if(fPy) delete [] fPy;
+ fPy = new Double32_t[GetNProngs()];
+ if(fPz) delete [] fPz;
+ fPz = new Double32_t[GetNProngs()];
+ memcpy(fPx,source.fPx,GetNProngs()*sizeof(Double32_t));
+ memcpy(fPy,source.fPy,GetNProngs()*sizeof(Double32_t));
+ memcpy(fPz,source.fPz,GetNProngs()*sizeof(Double32_t));
}
if(source.fPID) {
- fPID = new Double_t[5*GetNProngs()];
- memcpy(fPID,source.fPID,GetNProngs()*sizeof(Double_t));
+ if(fPID) delete [] fPID;
+ fPID = new Double32_t[fNPID];
+ memcpy(fPID,source.fPID,fNPID*sizeof(Double32_t));
}
if(source.fDCA) {
- fDCA = new Float_t[GetNProngs()*(GetNProngs()-1)/2];
- memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Float_t));
+ if(fDCA) delete [] fDCA;
+ fDCA = new Double32_t[fNDCA];
+ memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t));
}
}
return *this;
//
// Default Destructor
//
- if(fSecondaryVtx) delete fSecondaryVtx;
- if(fPx) delete [] fPx;
- if(fPy) delete [] fPy;
- if(fPz) delete [] fPz;
- if(fd0) delete [] fd0;
- if(fPID) delete [] fPID;
- if(fDCA) delete [] fDCA;
+ if(fPx) { delete [] fPx; fPx=NULL; }
+ if(fPy) { delete [] fPy; fPy=NULL; }
+ if(fPz) { delete [] fPz; fPz=NULL; }
+ if(fd0) { delete [] fd0; fd0=NULL; }
+ if(fPID) { delete [] fPID; fPID=NULL; }
+ if(fDCA) { delete [] fDCA; fDCA=NULL; }
+ if(fOwnSecondaryVtx) { delete fOwnSecondaryVtx; fOwnSecondaryVtx=NULL; }
}
//----------------------------------------------------------------------------
Double_t AliAODRecoDecay::Alpha() const
Double_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
return TMath::Sqrt(mass*mass+PProng(ip)*PProng(ip));
}
-//---------------------------------------------------------------------------
-/*Int_t AliAODRecoDecay::GetIndexProng(Int_t ip) const
-{
- //
- // Index of prong ip
+//--------------------------------------------------------------------------
+Bool_t AliAODRecoDecay::GetCovarianceXYZPxPyPz(Double_t cv[21]) const {
+ //
+ // This function returns the global covariance matrix of the track params
+ //
+ // Cov(x,x) ... : cv[0]
+ // Cov(y,x) ... : cv[1] cv[2]
+ // Cov(z,x) ... : cv[3] cv[4] cv[5]
+ // Cov(px,x)... : cv[6] cv[7] cv[8] cv[9]
+ // Cov(py,x)... : cv[10] cv[11] cv[12] cv[13] cv[14]
+ // Cov(pz,x)... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
+ //
+ // For XYZ we take the cov of the vertex, for PxPyPz we take the
+ // sum of the covs of PxPyPz from the daughters, for the moment
+ // we set the cov between position and momentum as the sum of
+ // the same cov from the daughters.
//
- if(!GetNProngs()) return 999999;
- UShort_t *indices = fSecondaryVtx->GetIndices();
- return indices[ip];
-}*/
+
+ Int_t j;
+ for(j=0;j<21;j++) cv[j]=0.;
+
+ if(!GetNDaughters()) {
+ AliError("No daughters available");
+ return kFALSE;
+ }
+
+ Double_t v[6];
+ AliAODVertex *secv=GetSecondaryVtx();
+ if(!secv) {
+ AliError("Vertex covariance matrix not available");
+ return kFALSE;
+ }
+ if(!secv->GetCovMatrix(v)) {
+ AliError("Vertex covariance matrix not available");
+ return kFALSE;
+ }
+
+ Double_t p[21]; for(j=0;j<21;j++) p[j]=0.;
+ Bool_t error=kFALSE;
+ for(Int_t i=1; i<GetNDaughters(); i++) {
+ AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
+ Double_t dcov[21];
+ if(!daugh->GetCovarianceXYZPxPyPz(dcov)) error=kTRUE;
+ for(j=0;j<21;j++) p[j] += dcov[j];
+ }
+ if(error) {
+ AliError("No covariance for at least one daughter")
+ return kFALSE;
+ }
+
+ for(j=0; j<21; j++) {
+ if(j<6) {
+ cv[j] = v[j];
+ } else {
+ cv[j] = p[j];
+ }
+ }
+
+ return kTRUE;
+}
//----------------------------------------------------------------------------
+UChar_t AliAODRecoDecay::GetITSClusterMap() const {
+ //
+ // We take the logical AND of the daughters cluster maps
+ // (only if all daughters have the bit for given layer, we set the bit)
+ //
+ UChar_t map=0;
+
+ if(!GetNDaughters()) {
+ AliError("No daughters available");
+ return map;
+ }
+
+ 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)
+ Int_t bit = 1;
+ for(Int_t i=0; i<GetNDaughters(); i++) {
+ AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
+ if(!TESTBIT(daugh->GetITSClusterMap(),l)) bit=0;
+ }
+ if(bit) SETBIT(map,l);
+ }
+
+ return map;
+}
+//--------------------------------------------------------------------------
+ULong_t AliAODRecoDecay::GetStatus() const {
+ //
+ // Same as for ITSClusterMap
+ //
+ ULong_t status=0;
+
+ if(!GetNDaughters()) {
+ AliError("No daughters available");
+ return status;
+ }
+
+ AliVTrack *daugh0 = (AliVTrack*)GetDaughter(0);
+ status = status&(daugh0->GetStatus());
+
+ for(Int_t i=1; i<GetNDaughters(); i++) {
+ AliVTrack *daugh = (AliVTrack*)GetDaughter(i);
+ status = status&(daugh->GetStatus());
+ }
+
+ return status;
+}
+//--------------------------------------------------------------------------
Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const
{
//
return mass;
}
+//----------------------------------------------------------------------------
+Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray) const
+{
+ //
+ // Check if this candidate is matched to a MC signal
+ // If no, return -1
+ // If yes, return label (>=0) of the AliAODMCParticle
+ //
+
+ Int_t ndg=GetNDaughters();
+ if(!ndg) {
+ AliError("No daughters available");
+ return -1;
+ }
+ if(ndg>10) {
+ AliError("Only decays with <10 daughters supported");
+ return -1;
+ }
+
+ Int_t dgLabels[10];
+
+ // loop on daughters and write the labels
+ for(Int_t i=0; i<ndg; i++) {
+ AliAODTrack *trk = (AliAODTrack*)GetDaughter(i);
+ dgLabels[i] = trk->GetLabel();
+ }
+
+ return MatchToMC(pdgabs,mcArray,dgLabels,ndg);
+}
+//----------------------------------------------------------------------------
+Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
+ Int_t dgLabels[10],Int_t ndg) const
+{
+ //
+ // Check if this candidate is matched to a MC signal
+ // If no, return -1
+ // If yes, return label (>=0) of the AliAODMCParticle
+ //
+
+ Int_t labMom[10];
+ Int_t i,lab,labMother,pdgMother;
+ AliAODMCParticle *part=0;
+ AliAODMCParticle *mother=0;
+ Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.;
+
+ // loop on daughter labels
+ for(i=0; i<ndg; i++) {
+ labMom[i]=-1;
+ lab = dgLabels[i];
+ if(lab<0) {
+ printf("daughter with negative label %d\n",lab);
+ return -1;
+ }
+ part = (AliAODMCParticle*)mcArray->At(lab);
+ if(!part) {
+ printf("no MC particle\n");
+ return -1;
+ }
+
+ // for the J/psi, check that the daughters are electrons
+ if(pdgabs==443) {
+ if(TMath::Abs(part->GetPdgCode())!=11) return -1;
+ }
+
+ mother = part;
+ while(mother->GetMother()>=0) {
+ labMother=mother->GetMother();
+ mother = (AliAODMCParticle*)mcArray->At(labMother);
+ if(!mother) {
+ printf("no MC mother particle\n");
+ break;
+ }
+ pdgMother = TMath::Abs(mother->GetPdgCode());
+ if(pdgMother==pdgabs) {
+ labMom[i]=labMother;
+ // keep sum of daughters' momenta, to check for mom conservation
+ pxSumDgs += part->Px();
+ pySumDgs += part->Py();
+ pzSumDgs += part->Pz();
+ break;
+ } else if(pdgMother>pdgabs || pdgMother<10) {
+ break;
+ }
+ }
+ if(labMom[i]==-1) return -1; // mother PDG not ok for this daughter
+ } // end loop on daughters
+
+ // check if the candidate is signal
+ labMother=labMom[0];
+ // all labels have to be the same and !=-1
+ for(i=0; i<ndg; i++) {
+ if(labMom[i]==-1) return -1;
+ if(labMom[i]!=labMother) return -1;
+ }
+
+
+ /*
+ // check that the mother decayed in <GetNDaughters()> prongs
+ Int_t ndg2 = TMath::Abs(mother->GetDaughter(1)-mother->GetDaughter(0))+1;
+ printf(" MC daughters %d\n",ndg2);
+ //if(ndg!=GetNDaughters()) return -1;
+ AliAODMCParticle* p1=(AliAODMCParticle*)(mcArray->At(mother->GetDaughter(1)));
+ AliAODMCParticle* p0=(AliAODMCParticle*)(mcArray->At(mother->GetDaughter(0)));
+ printf(">>>>>>>> pdg %d %d %d %d %d %d\n",mother->GetDaughter(0),mother->GetDaughter(1),dgLabels[0],dgLabels[1],p0->GetPdgCode(),p1->GetPdgCode());
+ */
+
+ // the above works only for non-resonant decays,
+ // it's better to check for mom conservation
+ mother = (AliAODMCParticle*)mcArray->At(labMother);
+ Double_t pxMother = mother->Px();
+ Double_t pyMother = mother->Py();
+ Double_t pzMother = mother->Pz();
+ // within 0.1%
+ if((TMath::Abs(pxMother-pxSumDgs)/(TMath::Abs(pxMother)+1.e-13)) > 0.001 &&
+ (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.001 &&
+ (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.001)
+ return -1;
+
+ return labMother;
+}
//---------------------------------------------------------------------------
void AliAODRecoDecay::Print(Option_t* /*option*/) const
{
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff