Bug corrected.
[u/mrichter/AliRoot.git] / STEER / AliAODRecoDecay.cxx
index 6037574..d7d192b 100644 (file)
@@ -25,6 +25,8 @@
 #include <TClonesArray.h>
 #include "AliLog.h"
 #include "AliVTrack.h"
+#include "AliExternalTrackParam.h"
+#include "AliNeutralTrackParam.h"
 #include "AliAODMCParticle.h"
 #include "AliAODRecoDecay.h"
 
@@ -421,6 +423,27 @@ ULong_t AliAODRecoDecay::GetStatus() const {
   return status;
 }
 //--------------------------------------------------------------------------
+Bool_t AliAODRecoDecay::PropagateToDCA(const AliVVertex* vtx,Double_t b,Double_t maxd,Double_t dz[2],Double_t covar[3]) 
+{
+  // compute impact parameters to the vertex vtx and their covariance matrix
+  // b is the Bz, needed to propagate correctly the track to vertex 
+  // return kFALSE is something went wrong
+
+  AliWarning("The AliAODRecoDecay momentum is not updated to the DCA");
+
+  Bool_t retval=kTRUE;
+  if(Charge()==0) {
+    // convert to AliNeutralTrackParam
+    AliNeutralTrackParam ntp(this);  
+    retval = ntp.PropagateToDCA(vtx,b,maxd,dz,covar);
+  } else {
+    // convert to AliExternalTrackParam
+    AliExternalTrackParam etp(this);  
+    retval = etp.PropagateToDCA(vtx,b,maxd,dz,covar);
+  }
+  return retval;
+}
+//--------------------------------------------------------------------------
 Double_t AliAODRecoDecay::ImpParXY(Double_t point[3]) const 
 {
   //
@@ -475,36 +498,44 @@ Double_t AliAODRecoDecay::InvMass2Prongs(Int_t ip1,Int_t ip2,
   return mass;
 }
 //----------------------------------------------------------------------------
-Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray) const
+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(!GetNDaughters()) {
+  // 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 = new Int_t[GetNDaughters()];
+  Int_t dgLabels[10];
 
   // loop on daughters and write the labels
-  for(Int_t i=0; i<GetNDaughters(); i++) {
-    AliVTrack *trk = (AliVTrack*)GetDaughter(i);
+  for(Int_t i=0; i<ndg; i++) {
+    AliAODTrack *trk = (AliAODTrack*)GetDaughter(i);
     dgLabels[i] = trk->GetLabel();
   }
 
-  Int_t labMother = MatchToMC(pdgabs,mcArray,dgLabels);
-
-  delete [] dgLabels; dgLabels=NULL;
-
-  return labMother;
+  return MatchToMC(pdgabs,mcArray,dgLabels,ndg,ndgCk,pdgDg);
 }
 //----------------------------------------------------------------------------
 Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
-                                Int_t *dgLabels) const
+                                Int_t dgLabels[10],Int_t ndg,
+                                Int_t ndgCk,Int_t *pdgDg) const
 {
   //
   // Check if this candidate is matched to a MC signal
@@ -512,82 +543,101 @@ Int_t AliAODRecoDecay::MatchToMC(Int_t pdgabs,TClonesArray *mcArray,
   // If yes, return label (>=0) of the AliAODMCParticle
   // 
 
-  if(!GetNDaughters()) {
-    AliError("No daughters available");
-    return -1;
-  }
-  
-  Int_t *labMom = new Int_t[GetNDaughters()];
-  Int_t i,lab,labMother,pdgMother;
+  Int_t labMom[10]={0,0,0,0,0,0,0,0,0,0};
+  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; i<GetNDaughters(); i++) {
+  for(i=0; i<ndg; i++) {
     labMom[i]=-1;
     lab = dgLabels[i];
     if(lab<0) {
-      printf("daughter with negative label\n");
-      continue;
+      printf("daughter with negative label %d\n",lab);
+      return -1;
     }
     part = (AliAODMCParticle*)mcArray->At(lab);
     if(!part) { 
       printf("no MC particle\n");
-      continue;
+      return -1;
     }
-    // keep sum of daughters' momenta, to check for mom conservation
-    pxSumDgs += part->Px();
-    pySumDgs += part->Py();
-    pzSumDgs += part->Pz();
+
+    // check the PDG of the daughter, if requested
+    if(ndgCk>0) {
+      pdgPart=TMath::Abs(part->GetPdgCode());
+      for(j=0; j<ndg; j++) {
+       if(!pdgUsed[j] && pdgPart==pdgDg[j]) {
+         pdgUsed[j]=kTRUE;
+         break;
+       }
+      }
+    }
+
     // for the J/psi, check that the daughters are electrons
-    if(pdgabs==443 && TMath::Abs(part->GetPdgCode())!=11) continue;
+    if(pdgabs==443) {
+      if(TMath::Abs(part->GetPdgCode())!=11) return -1;
+    }
 
-    while(part->GetMother()>=0) {
-      labMother=part->GetMother();
-      part = (AliAODMCParticle*)mcArray->At(labMother);
-      if(!part) {
+    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(part->GetPdgCode());
+      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
-  Bool_t isSignal=kTRUE;
   labMother=labMom[0];
   // all labels have to be the same and !=-1
-  for(i=0; i<GetNDaughters(); i++) {
-    if(labMom[i]==-1 || labMom[i]!=labMother) isSignal=kFALSE;
+  for(i=0; i<ndg; i++) {
+    if(labMom[i]==-1)        return -1;
+    if(labMom[i]!=labMother) return -1;
   }
 
-  delete [] labMom; labMom=NULL;
-
-  if(!isSignal) return -1;
-
-  /*  
+  // check that all daughter PDGs are matched
+  if(ndgCk>0) {
+    for(i=0; i<ndg; i++) {
+      if(pdgUsed[i]==kFALSE) return -1;
+    }
+  }
+  
+  /*
   // check that the mother decayed in <GetNDaughters()> prongs
-  Int_t ndg = TMath::Abs(part->GetDaughter(1)-part->GetDaughter(0))+1;
-  if(ndg!=GetNDaughters()) return -1;
-  AliAODMCParticle* p1=(AliAODMCParticle*)(mcArray->At(part->GetDaughter(1)));
-  AliAODMCParticle* p0=(AliAODMCParticle*)(mcArray->At(part->GetDaughter(0)));
-  printf("pdg %d  %d %d %d   %d %d\n",part->GetDaughter(1),part->GetDaughter(0),dgLabels[0],dgLabels[1],p0->GetPdgCode(),p1->GetPdgCode());
+  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
-  part = (AliAODMCParticle*)mcArray->At(labMother);
-  Double_t pxMother = part->Px();
-  Double_t pyMother = part->Py();
-  Double_t pzMother = part->Pz();
+  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) 
+  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;