X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliAODRecoDecay.cxx;h=9e4911da85120c3b321f3cc5706dcf827c3cfca4;hb=2aa9925fafff1e19cfe379b7e948ea12ef815241;hp=51326f35fe808e0f42ef02b10600f391b927594c;hpb=341952b8438301053cef475ba6066ea16b24a90b;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliAODRecoDecay.cxx b/STEER/AliAODRecoDecay.cxx index 51326f35fe8..9e4911da851 100644 --- a/STEER/AliAODRecoDecay.cxx +++ b/STEER/AliAODRecoDecay.cxx @@ -22,14 +22,17 @@ #include #include -#include "AliVParticle.h" +#include +#include "AliLog.h" +#include "AliVTrack.h" +#include "AliAODMCParticle.h" #include "AliAODRecoDecay.h" ClassImp(AliAODRecoDecay) //-------------------------------------------------------------------------- AliAODRecoDecay::AliAODRecoDecay() : - AliVParticle(), + AliVTrack(), fSecondaryVtx(0x0), fOwnSecondaryVtx(0x0), fCharge(0), @@ -49,7 +52,7 @@ 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), @@ -79,7 +82,7 @@ AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs, AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs, Short_t charge, Double_t *d0) : - AliVParticle(), + AliVTrack(), fSecondaryVtx(vtx2), fOwnSecondaryVtx(0x0), fCharge(charge), @@ -99,7 +102,7 @@ AliAODRecoDecay::AliAODRecoDecay(AliAODVertex *vtx2,Int_t nprongs, } //-------------------------------------------------------------------------- AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) : - AliVParticle(source), + AliVTrack(source), fSecondaryVtx(source.fSecondaryVtx), fOwnSecondaryVtx(source.fOwnSecondaryVtx), fCharge(source.fCharge), @@ -114,23 +117,23 @@ AliAODRecoDecay::AliAODRecoDecay(const AliAODRecoDecay &source) : // 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 Double_t[GetNProngs()*(GetNProngs()-1)/2]; - memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double_t)); + fDCA = new Double32_t[fNDCA]; + memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t)); } } } @@ -151,28 +154,28 @@ AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source) fRunNumber = source.fRunNumber; if(source.GetNProngs()>0) { if(fd0)delete [] fd0; - 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) { if(fPx) delete [] fPx; - fPx = new Double_t[GetNProngs()]; + fPx = new Double32_t[GetNProngs()]; if(fPy) delete [] fPy; - fPy = new Double_t[GetNProngs()]; + fPy = new Double32_t[GetNProngs()]; if(fPz) delete [] fPz; - 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)); + 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) { if(fPID) delete [] 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) { if(fDCA) delete [] fDCA; - fDCA = new Double_t[GetNProngs()*(GetNProngs()-1)/2]; - memcpy(fDCA,source.fDCA,(GetNProngs()*(GetNProngs()-1)/2)*sizeof(Double32_t)); + fDCA = new Double32_t[fNDCA]; + memcpy(fDCA,source.fDCA,fNDCA*sizeof(Double32_t)); } } return *this; @@ -311,17 +314,113 @@ Double_t AliAODRecoDecay::EProng(Int_t ip,UInt_t pdg) 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 = GetSecondaryVtx()->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; iGetCovarianceXYZPxPyPz(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; iGetITSClusterMap(),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; iGetStatus()); + } + + return status; +} +//-------------------------------------------------------------------------- Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const { // @@ -375,6 +474,151 @@ Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2, return mass; } +//---------------------------------------------------------------------------- +Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray, + Int_t ndgCk,Int_t *pdgDg) const +{ + // + // Check if this candidate is matched to a MC signal + // If no, return -1 + // If yes, return label (>=0) of the AliAODMCParticle + // + // if ndgCk>0, checks also daughters PDGs + // + Int_t ndg=GetNDaughters(); + if(!ndg) { + AliError("No daughters available"); + return -1; + } + if(ndg>10) { + AliError("Only decays with <10 daughters supported"); + return -1; + } + if(ndgCk>0 && ndgCk!=ndg) { + AliError("Wrong number of daughter PDGs passed"); + return -1; + } + + Int_t dgLabels[10]; + + // loop on daughters and write the labels + for(Int_t i=0; iGetLabel(); + } + + return MatchToMC(pdgabs,mcArray,dgLabels,ndg,ndgCk,pdgDg); +} +//---------------------------------------------------------------------------- +Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray, + Int_t dgLabels[10],Int_t ndg, + Int_t ndgCk,Int_t *pdgDg) 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,j,lab,labMother,pdgMother,pdgPart; + AliAODMCParticle *part=0; + AliAODMCParticle *mother=0; + Double_t pxSumDgs=0.,pySumDgs=0.,pzSumDgs=0.; + Bool_t pdgUsed[10]={kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE,kFALSE}; + + // loop on daughter labels + for(i=0; iAt(lab); + if(!part) { + printf("no MC particle\n"); + return -1; + } + + // check the PDG of the daughter, if requested + if(ndgCk>0) { + pdgPart=TMath::Abs(part->GetPdgCode()); + for(j=0; jGetPdgCode())!=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; i0) { + for(i=0; i 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.00001 && + (TMath::Abs(pyMother-pySumDgs)/(TMath::Abs(pyMother)+1.e-13)) > 0.00001 && + (TMath::Abs(pzMother-pzSumDgs)/(TMath::Abs(pzMother)+1.e-13)) > 0.00001) + return -1; + + return labMother; +} //--------------------------------------------------------------------------- void AliAODRecoDecay::Print(Option_t* /*option*/) const {