I improved (by a factor 2.5) the speed of the MatchToMC method
[u/mrichter/AliRoot.git] / STEER / AliAODRecoDecay.cxx
index 34c4feacbbd0521833f2609f8209db1e0683dda4..0c5c4fc17fb1d75f451eb9c57485577be4e9776c 100644 (file)
 
 #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),
@@ -42,15 +46,15 @@ AliAODRecoDecay::AliAODRecoDecay() :
   //
   // 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),
@@ -78,8 +82,9 @@ 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),
   fNProngs(nprongs), fNDCA(0), fNPID(0),
   fPx(0x0), fPy(0x0), fPz(0x0),
@@ -97,8 +102,9 @@ 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),
   fNProngs(source.fNProngs), fNDCA(source.fNDCA), fNPID(source.fNPID),
   fPx(0x0), fPy(0x0), fPz(0x0),
@@ -111,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 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));
     }
   }
 }
@@ -139,6 +145,7 @@ AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
   //
   if(&source == this) return *this;
   fSecondaryVtx = source.fSecondaryVtx;
+  fOwnSecondaryVtx = source.fOwnSecondaryVtx;
   fCharge = source.fCharge;
   fNProngs = source.fNProngs;
   fNDCA = source.fNDCA;
@@ -146,23 +153,29 @@ AliAODRecoDecay &AliAODRecoDecay::operator=(const AliAODRecoDecay &source)
   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;
@@ -172,13 +185,13 @@ AliAODRecoDecay::~AliAODRecoDecay() {
   //  
   // 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 
@@ -301,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 = 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 
 {
   //
@@ -365,6 +474,126 @@ Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
 
   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 
 {