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
// 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;
+ }
+
+ // 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;
+ }
+ }
}
- // keep sum of daughters' momenta, to check for mom conservation
- pxSumDgs += part->Px();
- pySumDgs += part->Py();
- pzSumDgs += part->Pz();
+
// 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();
- // within 1 MeV
- if(TMath::Abs(pxMother-pxSumDgs) > 0.001 ||
- TMath::Abs(pyMother-pySumDgs) > 0.001 ||
- TMath::Abs(pzMother-pzSumDgs) > 0.001) return -1;
+ 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;
}