]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/vertexingHF/AliVertexingHFUtils.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliVertexingHFUtils.cxx
index 4e1814ae2974c82e817eb5a8b92ccbdf540384f0..685a62055bcc3ab62c650424e405a24127f63a2b 100644 (file)
@@ -239,6 +239,37 @@ Int_t AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(TClonesArray*
   }
   return nChargedMC;
 }
+
+
+//______________________________________________________________________
+Double_t AliVertexingHFUtils::GetVZEROAEqualizedMultiplicity(AliAODEvent* ev){
+  //
+  // Method to get VZERO-A equalized multiplicty as done in AliCentralitySelectionTask
+  //  getting the equalized VZERO factors from the tender or AOD
+  // http://git.cern.ch/pubweb/AliRoot.git/blob/HEAD:/ANALYSIS/AliCentralitySelectionTask.cxx#l1345
+
+  Double_t multV0AEq=0;
+  for(Int_t iCh = 32; iCh < 64; ++iCh) {
+    Double_t mult = ev->GetVZEROEqMultiplicity(iCh);
+    multV0AEq += mult;
+  }
+  return multV0AEq;
+}
+
+//______________________________________________________________________
+Double_t AliVertexingHFUtils::GetVZEROCEqualizedMultiplicity(AliAODEvent* ev){
+  // Method to get VZERO-C equalized multiplicty as done in AliCentralitySelectionTask
+  //  getting the equalized VZERO factors from the tender or AOD
+  // http://git.cern.ch/pubweb/AliRoot.git/blob/HEAD:/ANALYSIS/AliCentralitySelectionTask.cxx#l1345
+
+  Double_t multV0CEq=0;
+  for(Int_t iCh = 0; iCh < 32; ++iCh) {
+    Double_t mult = ev->GetVZEROEqMultiplicity(iCh);
+    multV0CEq += mult;
+  }
+  return multV0CEq;
+}
+
 //______________________________________________________________________
 void AliVertexingHFUtils::AveragePt(Float_t& averagePt, Float_t& errorPt,Float_t ptmin,Float_t ptmax, TH2F* hMassD, Float_t massFromFit, Float_t sigmaFromFit, TF1* funcB2, Float_t sigmaRangeForSig,Float_t sigmaRangeForBkg, Float_t minMass, Float_t maxMass, Int_t rebin){
 
@@ -1020,6 +1051,173 @@ Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCPartic
   return -1;
 
 }
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
+  // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
+
+  if(label<0) return -1;
+  TParticle* mcPart = stack->Particle(label);
+  Int_t pdgD=mcPart->GetPdgCode();
+  if(TMath::Abs(pdgD)!=411) return -1;
+
+  Int_t nDau=mcPart->GetNDaughters();
+  Int_t labelFirstDau = mcPart->GetDaughter(0);
+  Int_t nKaons=0;
+  Int_t nPions=0;
+  Double_t sumPxDau=0.;
+  Double_t sumPyDau=0.;
+  Double_t sumPzDau=0.;
+  Int_t nFoundKpi=0;
+  Bool_t isPhi=kFALSE;
+  Bool_t isk0st=kFALSE;
+  
+  if(nDau==3 || nDau==2){
+    for(Int_t iDau=0; iDau<nDau; iDau++){
+      Int_t indDau = labelFirstDau+iDau;
+      if(indDau<0) return -1;
+      TParticle* dau=stack->Particle(indDau);
+      if(!dau) return -1;
+      Int_t pdgdau=dau->GetPdgCode();
+      if(TMath::Abs(pdgdau)==321){
+       nKaons++;
+       sumPxDau+=dau->Px();
+       sumPyDau+=dau->Py();
+       sumPzDau+=dau->Pz();
+       arrayDauLab[nFoundKpi++]=indDau;
+       if(nFoundKpi>3) return -1;
+      }else if(TMath::Abs(pdgdau)==211){
+       nPions++;
+       sumPxDau+=dau->Px();
+       sumPyDau+=dau->Py();
+       sumPzDau+=dau->Pz();
+       arrayDauLab[nFoundKpi++]=indDau;
+       if(nFoundKpi>3) return -1;
+      }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){
+       if(TMath::Abs(pdgdau)==313) isk0st=kTRUE;
+       if(TMath::Abs(pdgdau)==333) isPhi=kTRUE;
+       Int_t nResDau=dau->GetNDaughters();
+       if(nResDau!=2) return -1;
+       Int_t indFirstResDau=dau->GetDaughter(0);
+       for(Int_t resDau=0; resDau<2; resDau++){
+         Int_t indResDau=indFirstResDau+resDau;
+         if(indResDau<0) return -1;
+         TParticle* resdau=stack->Particle(indResDau);
+         if(!resdau) return -1;
+         Int_t pdgresdau=resdau->GetPdgCode();
+         if(TMath::Abs(pdgresdau)==321){
+           sumPxDau+=resdau->Px();
+           sumPyDau+=resdau->Py();
+           sumPzDau+=resdau->Pz();
+           nKaons++;
+           arrayDauLab[nFoundKpi++]=indResDau;
+           if(nFoundKpi>3) return -1;
+         }
+         if(TMath::Abs(pdgresdau)==211){
+           sumPxDau+=resdau->Px();
+           sumPyDau+=resdau->Py();
+           sumPzDau+=resdau->Pz();
+           nPions++;
+           arrayDauLab[nFoundKpi++]=indResDau;
+           if(nFoundKpi>3) return -1;
+         }
+         }
+      }else{
+       return -1;
+      }
+    }
+    if(nPions!=1) return -1;
+    if(nKaons!=2) return -1;
+    if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
+    if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
+    if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
+    if(nDau==3) return 3;
+    else if(nDau==2){
+      if(isk0st) return 2;
+      if(isPhi) return 1;
+    }
+  }
+  
+  return -1;
+  
+}
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckDplusK0spiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
+  // Checks the Dplus->V0+pion decay channel. Returns 1 if success, -1 otherwise
+
+  if(label<0) return -1;
+  TParticle* mcPart = stack->Particle(label);
+  Int_t pdgD=mcPart->GetPdgCode();
+  if(TMath::Abs(pdgD)!=411) return -1;
+
+  Int_t nDau=mcPart->GetNDaughters();
+  Int_t labelFirstDau = mcPart->GetDaughter(0);
+  Int_t nPions=0;
+  Double_t sumPxDau=0.;
+  Double_t sumPyDau=0.;
+  Double_t sumPzDau=0.;
+  Int_t nFoundpi=0;
+
+  Int_t codeV0=-1;
+  if(nDau==2){
+    for(Int_t iDau=0; iDau<nDau; iDau++){
+      Int_t indDau = labelFirstDau+iDau;
+      if(indDau<0) return -1;
+      TParticle* dau=stack->Particle(indDau);
+      if(!dau) return -1;
+      Int_t pdgdau=dau->GetPdgCode();
+      if(TMath::Abs(pdgdau)==211){
+       nPions++;
+       sumPxDau+=dau->Px();
+       sumPyDau+=dau->Py();
+       sumPzDau+=dau->Pz();
+       arrayDauLab[nFoundpi++]=indDau;
+       if(nFoundpi>3) return -1;
+      }else if(TMath::Abs(pdgdau)==311){
+       codeV0=TMath::Abs(pdgdau);
+       TParticle* v0=dau;
+       if(codeV0==311){
+         Int_t nK0Dau=dau->GetNDaughters();
+         if(nK0Dau!=1) return -1;
+         Int_t indK0s=dau->GetDaughter(0);
+         if(indK0s<0) return -1;
+         v0=stack->Particle(indK0s);
+         if(!v0) return -1;
+         Int_t pdgK0sl=v0->GetPdgCode();
+         codeV0=TMath::Abs(pdgK0sl);
+       }
+       Int_t nV0Dau=v0->GetNDaughters();
+       if(nV0Dau!=2) return -1;
+       Int_t indFirstV0Dau=v0->GetDaughter(0);
+       for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
+         Int_t indV0Dau=indFirstV0Dau+v0Dau;
+         if(indV0Dau<0) return -1;
+         TParticle* v0dau=stack->Particle(indV0Dau);
+         if(!v0dau) return -1;
+         Int_t pdgv0dau=v0dau->GetPdgCode();
+         if(TMath::Abs(pdgv0dau)==211){
+           sumPxDau+=v0dau->Px();
+           sumPyDau+=v0dau->Py();
+           sumPzDau+=v0dau->Pz();
+           nPions++;
+           arrayDauLab[nFoundpi++]=indV0Dau;
+           if(nFoundpi>3) return -1;
+         }
+       }
+      }else{
+       return -1;
+      }
+    }
+    if(nPions!=3) return -1;
+    if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
+    if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
+    if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
+    if(codeV0==310) return 1;
+  }
+  return -1;
+  
+}
+
 //____________________________________________________________________________
 Int_t AliVertexingHFUtils::CheckDsDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
   // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
@@ -1109,6 +1307,93 @@ Int_t AliVertexingHFUtils::CheckDsDecay(AliStack* stack, Int_t label, Int_t* arr
   return -1;
   
 }
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckDsK0sKDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
+  // Checks the Ds->K0s+S decay channel. Returns 1 in case of success, otherwise -1
+
+  if(label<0) return -1;
+  TParticle* mcPart = stack->Particle(label);
+  Int_t pdgD=mcPart->GetPdgCode();
+  if(TMath::Abs(pdgD)!=431) return -1;
+
+  Int_t nDau=mcPart->GetNDaughters();
+  Int_t labelFirstDau = mcPart->GetDaughter(0);
+  Int_t nPions=0;
+  Int_t nKaons=0;
+  Double_t sumPxDau=0.;
+  Double_t sumPyDau=0.;
+  Double_t sumPzDau=0.;
+  Int_t nFoundKpi=0;
+
+  Int_t codeV0=-1;
+  if(nDau==2){
+    for(Int_t iDau=0; iDau<nDau; iDau++){
+      Int_t indDau = labelFirstDau+iDau;
+      if(indDau<0) return -1;
+      TParticle* dau=stack->Particle(indDau);
+      if(!dau) return -1;
+      Int_t pdgdau=dau->GetPdgCode();
+      if(TMath::Abs(pdgdau)==211){
+       nPions++;
+       sumPxDau+=dau->Px();
+       sumPyDau+=dau->Py();
+       sumPzDau+=dau->Pz();
+       arrayDauLab[nFoundKpi++]=indDau;
+       if(nFoundKpi>3) return -1;
+      }else if(TMath::Abs(pdgdau)==321){
+       nKaons++;
+       sumPxDau+=dau->Px();
+       sumPyDau+=dau->Py();
+       sumPzDau+=dau->Pz();
+       arrayDauLab[nFoundKpi++]=indDau;
+       if(nFoundKpi>3) return -1;
+      }else if(TMath::Abs(pdgdau)==311){
+       codeV0=TMath::Abs(pdgdau);
+       TParticle* v0=dau;
+       if(codeV0==311){
+         Int_t nK0Dau=dau->GetNDaughters();
+         if(nK0Dau!=1) return -1;
+         Int_t indK0s=dau->GetDaughter(0);
+         if(indK0s<0) return -1;
+         v0=stack->Particle(indK0s);
+         if(!v0) return -1;
+         Int_t pdgK0sl=v0->GetPdgCode();
+         codeV0=TMath::Abs(pdgK0sl);
+       }
+       Int_t nV0Dau=v0->GetNDaughters();
+       if(nV0Dau!=2) return -1;
+       Int_t indFirstV0Dau=v0->GetDaughter(0);
+       for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
+         Int_t indV0Dau=indFirstV0Dau+v0Dau;
+         if(indV0Dau<0) return -1;
+         TParticle* v0dau=stack->Particle(indV0Dau);
+         if(!v0dau) return -1;
+         Int_t pdgv0dau=v0dau->GetPdgCode();
+         if(TMath::Abs(pdgv0dau)==211){
+           sumPxDau+=v0dau->Px();
+           sumPyDau+=v0dau->Py();
+           sumPzDau+=v0dau->Pz();
+           nPions++;
+           arrayDauLab[nFoundKpi++]=indV0Dau;
+           if(nFoundKpi>3) return -1;
+         }
+       }
+      }else{
+       return -1;
+      }
+    }
+    if(nPions!=2) return -1;
+    if(nKaons!=1) return -1;
+    if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
+    if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
+    if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
+    if(codeV0==310) return 1;
+  }
+  return -1;
+  
+}
+
 //____________________________________________________________________________
 Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
   // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case
@@ -1176,7 +1461,7 @@ Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle
            arrayDauLab[nFoundKpi++]=indResDau;
            if(nFoundKpi>3) return -1;
          }
-         }
+       }
       }else{
        return -1;
       }
@@ -1344,3 +1629,257 @@ Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCPartic
   
 }
 
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckLcpKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
+  // Checks the Lc->pKpi decay channel. Returns 1 for non-resonant decays and 2, 3 or 4 for resonant ones, -1 in other cases
+
+  if(label<0) return -1;
+  TParticle* mcPart = stack->Particle(label);
+  Int_t pdgD=mcPart->GetPdgCode();
+  if(TMath::Abs(pdgD)!=4122) return -1;
+
+  Int_t nDau=mcPart->GetNDaughters();
+  Int_t labelFirstDau = mcPart->GetDaughter(0);
+  Int_t nKaons=0;
+  Int_t nPions=0;
+  Int_t nProtons=0;
+  Double_t sumPxDau=0.;
+  Double_t sumPyDau=0.;
+  Double_t sumPzDau=0.;
+  Int_t nFoundpKpi=0;
+
+  Int_t codeRes=-1;
+  if(nDau==3 || nDau==2){
+    for(Int_t iDau=0; iDau<nDau; iDau++){
+      Int_t indDau = labelFirstDau+iDau;
+      if(indDau<0) return -1;
+      TParticle* dau=stack->Particle(indDau);
+      if(!dau) return -1;
+      Int_t pdgdau=dau->GetPdgCode();
+      if(TMath::Abs(pdgdau)==321){
+       nKaons++;
+       sumPxDau+=dau->Px();
+       sumPyDau+=dau->Py();
+       sumPzDau+=dau->Pz();
+       arrayDauLab[nFoundpKpi++]=indDau;
+       if(nFoundpKpi>3) return -1;
+      }else if(TMath::Abs(pdgdau)==211){
+       nPions++;
+       sumPxDau+=dau->Px();
+       sumPyDau+=dau->Py();
+       sumPzDau+=dau->Pz();
+       arrayDauLab[nFoundpKpi++]=indDau;
+       if(nFoundpKpi>3) return -1;
+      }else if(TMath::Abs(pdgdau)==2212){
+       nProtons++;
+       sumPxDau+=dau->Px();
+       sumPyDau+=dau->Py();
+       sumPzDau+=dau->Pz();
+       arrayDauLab[nFoundpKpi++]=indDau;
+       if(nFoundpKpi>3) return -1;
+      }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 || 
+              TMath::Abs(pdgdau)==2224){
+       codeRes=TMath::Abs(pdgdau);
+       Int_t nResDau=dau->GetNDaughters();
+       if(nResDau!=2) return -1;
+       Int_t indFirstResDau=dau->GetDaughter(0);
+       for(Int_t resDau=0; resDau<2; resDau++){
+         Int_t indResDau=indFirstResDau+resDau;
+         if(indResDau<0) return -1;
+         TParticle* resdau=stack->Particle(indResDau);
+         if(!resdau) return -1;
+         Int_t pdgresdau=resdau->GetPdgCode();
+         if(TMath::Abs(pdgresdau)==321){
+           sumPxDau+=resdau->Px();
+           sumPyDau+=resdau->Py();
+           sumPzDau+=resdau->Pz();
+           nKaons++;
+           arrayDauLab[nFoundpKpi++]=indResDau;
+           if(nFoundpKpi>3) return -1;
+         }else if(TMath::Abs(pdgresdau)==211){
+           sumPxDau+=resdau->Px();
+           sumPyDau+=resdau->Py();
+           sumPzDau+=resdau->Pz();
+           nPions++;
+           arrayDauLab[nFoundpKpi++]=indResDau;
+           if(nFoundpKpi>3) return -1;
+         }else if(TMath::Abs(pdgresdau)==2212){
+           sumPxDau+=resdau->Px();
+           sumPyDau+=resdau->Py();
+           sumPzDau+=resdau->Pz();
+           nProtons++;
+           arrayDauLab[nFoundpKpi++]=indResDau;
+           if(nFoundpKpi>3) return -1;
+         }
+       }
+      }else{
+       return -1;
+      }
+    }
+    if(nPions!=1) return -1;
+    if(nKaons!=1) return -1;
+    if(nProtons!=1) return -1;
+    if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
+    if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
+    if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
+    if(nDau==3) return 1;
+    else if(nDau==2){
+      if(codeRes==313) return 2;
+      else if(codeRes==2224) return 3;
+      else if(codeRes==3124) return 4;
+    }  
+  }
+  return -1;
+  
+}
+
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
+  // Checks the Lc->V0+bachelor decay channel. Returns 1 for pK0s, 2 for piLambda, 3 for pK0l -1 in other cases
+
+  if(label<0) return -1;
+  TParticle* mcPart = stack->Particle(label);
+  Int_t pdgD=mcPart->GetPdgCode();
+  if(TMath::Abs(pdgD)!=4122) return -1;
+
+  Int_t nDau=mcPart->GetNDaughters();
+  Int_t labelFirstDau = mcPart->GetDaughter(0);
+  Int_t nPions=0;
+  Int_t nProtons=0;
+  Double_t sumPxDau=0.;
+  Double_t sumPyDau=0.;
+  Double_t sumPzDau=0.;
+  Int_t nFoundppi=0;
+
+  Int_t codeV0=-1;
+  if(nDau==2){
+    for(Int_t iDau=0; iDau<nDau; iDau++){
+      Int_t indDau = labelFirstDau+iDau;
+      if(indDau<0) return -1;
+      TParticle* dau=stack->Particle(indDau);
+      if(!dau) return -1;
+      Int_t pdgdau=dau->GetPdgCode();
+      if(TMath::Abs(pdgdau)==211){
+       nPions++;
+       sumPxDau+=dau->Px();
+       sumPyDau+=dau->Py();
+       sumPzDau+=dau->Pz();
+       arrayDauLab[nFoundppi++]=indDau;
+       if(nFoundppi>3) return -1;
+      }else if(TMath::Abs(pdgdau)==2212){
+       nProtons++;
+       sumPxDau+=dau->Px();
+       sumPyDau+=dau->Py();
+       sumPzDau+=dau->Pz();
+       arrayDauLab[nFoundppi++]=indDau;
+       if(nFoundppi>3) return -1;
+      }else if(TMath::Abs(pdgdau)==311 ||  TMath::Abs(pdgdau)==3122){
+       codeV0=TMath::Abs(pdgdau);
+       TParticle* v0=dau;
+       if(codeV0==311){
+         Int_t nK0Dau=dau->GetNDaughters();
+         if(nK0Dau!=1) return -1;
+         Int_t indK0s=dau->GetDaughter(0);
+         if(indK0s<0) return -1;
+         v0=stack->Particle(indK0s);
+         if(!v0) return -1;
+         Int_t pdgK0sl=v0->GetPdgCode();
+         codeV0=TMath::Abs(pdgK0sl);
+       }
+       Int_t nV0Dau=v0->GetNDaughters();
+       if(nV0Dau!=2) return -1;
+       Int_t indFirstV0Dau=v0->GetDaughter(0);
+       for(Int_t v0Dau=0; v0Dau<2; v0Dau++){
+         Int_t indV0Dau=indFirstV0Dau+v0Dau;
+         if(indV0Dau<0) return -1;
+         TParticle* v0dau=stack->Particle(indV0Dau);
+         if(!v0dau) return -1;
+         Int_t pdgv0dau=v0dau->GetPdgCode();
+         if(TMath::Abs(pdgv0dau)==211){
+           sumPxDau+=v0dau->Px();
+           sumPyDau+=v0dau->Py();
+           sumPzDau+=v0dau->Pz();
+           nPions++;
+           arrayDauLab[nFoundppi++]=indV0Dau;
+           if(nFoundppi>3) return -1;
+         }else if(TMath::Abs(pdgv0dau)==2212){
+           sumPxDau+=v0dau->Px();
+           sumPyDau+=v0dau->Py();
+           sumPzDau+=v0dau->Pz();
+           nProtons++;
+           arrayDauLab[nFoundppi++]=indV0Dau;
+           if(nFoundppi>3) return -1;
+         }
+       }
+      }else{
+       return -1;
+      }
+    }
+    if(nPions!=2) return -1;
+    if(nProtons!=1) return -1;
+    if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
+    if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
+    if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
+    if(codeV0==310) return 1;
+    else if(codeV0==3122) return 2;
+  }
+  return -1;
+}
+
+//__________________________________xic______________________________________
+Int_t AliVertexingHFUtils::CheckXicXipipiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
+  // Checks the Xic decay channel. Returns 1 for Xic->Xipipi, -1 in other cases
+
+  if(label<0) return -1;
+  TParticle* mcPart = stack->Particle(label);
+  Int_t pdgD=mcPart->GetPdgCode();
+  if(TMath::Abs(pdgD)!=4232) return -1;
+
+  Int_t nDau=mcPart->GetNDaughters();
+  if(nDau!=3) return -1;
+
+  Int_t labelFirstDau = mcPart->GetDaughter(0);
+  Int_t nXi=0;
+  Int_t nPions=0;
+  Double_t sumPxDau=0.;
+  Double_t sumPyDau=0.;
+  Double_t sumPzDau=0.;
+  Int_t nFoundXi=0;
+
+  for(Int_t iDau=0; iDau<nDau; iDau++){
+    Int_t indDau = labelFirstDau+iDau;
+    if(indDau<0) return -1;
+    TParticle* dau=stack->Particle(indDau);
+    if(!dau) return -1;
+    Int_t pdgdau=dau->GetPdgCode();
+    if(TMath::Abs(pdgdau)==3312){
+      if(pdgD*pdgdau<0) return -1;
+      sumPxDau+=dau->Px();
+      sumPyDau+=dau->Py();
+      sumPzDau+=dau->Pz();
+      nXi++;
+      arrayDauLab[nFoundXi++]=indDau;
+      
+    }
+    if(TMath::Abs(pdgdau)==211){
+      if(pdgD*pdgdau<0) return -1;
+      nPions++;
+      sumPxDau+=dau->Px();
+      sumPyDau+=dau->Py();
+      sumPzDau+=dau->Pz();
+      arrayDauLab[nFoundXi++]=indDau;
+      if(nFoundXi>3) return -1;
+    }
+  }
+  
+  if(nPions!=2) return -1;
+  if(nXi!=1) return -1;
+  if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2;
+  if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2;
+  if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2;
+  return 1;
+  
+}