#include <TMath.h>
#include <TRandom.h>
+#include <TProfile.h>
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TF1.h>
+#include <TParticle.h>
+#include "AliStack.h"
#include "AliAODEvent.h"
+#include "AliAODMCHeader.h"
+#include "AliGenEventHeader.h"
#include "AliAODMCParticle.h"
+#include "AliAODRecoDecayHF.h"
#include "AliVertexingHFUtils.h"
/* $Id$ */
}
+//______________________________________________________________________
+void AliVertexingHFUtils::ComputeSignificance(Double_t signal, Double_t errsignal, Double_t background, Double_t errbackground, Double_t &significance,Double_t &errsignificance){
+ // calculate significance from S, B and errors
+
+
+ Double_t errSigSq=errsignal*errsignal;
+ Double_t errBkgSq=errbackground*errbackground;
+ Double_t sigPlusBkg=signal+background;
+ if(sigPlusBkg>0. && signal>0.){
+ significance = signal/TMath::Sqrt(signal+background);
+ errsignificance = significance*TMath::Sqrt((errSigSq+errBkgSq)/(4.*sigPlusBkg*sigPlusBkg)+(background/sigPlusBkg)*errSigSq/signal/signal);
+ }else{
+ significance=0.;
+ errsignificance=0.;
+ }
+ return;
+
+}
//______________________________________________________________________
Double_t AliVertexingHFUtils::Pol(Double_t x, Int_t k){
// compute chi from polynomial approximation
}
return nChargedMC;
}
+
+
//______________________________________________________________________
-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, Int_t rebin){
+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){
// Compute <pt> from 2D histogram M vs pt
// printf("Signal Fit Limits = %f %f\n",minMassSigBin,maxMassSigBin);
Double_t maxMassBkgLow=massFromFit-sigmaRangeForBkg*sigmaFromFit;
- Int_t minBinBkgLow=2;
+ Int_t minBinBkgLow=TMath::Max(hMassD->GetYaxis()->FindBin(minMass),2);
Int_t maxBinBkgLow=hMassD->GetYaxis()->FindBin(maxMassBkgLow);
Double_t minMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgLow);
Double_t maxMassBkgLowBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgLow)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgLow);
Double_t minMassBkgHi=massFromFit+sigmaRangeForBkg*sigmaFromFit;
Int_t minBinBkgHi=hMassD->GetYaxis()->FindBin(minMassBkgHi);
- Int_t maxBinBkgHi=hMassD->GetNbinsY()-1;
+ Int_t maxBinBkgHi=TMath::Min(hMassD->GetYaxis()->FindBin(maxMass),hMassD->GetNbinsY()-1);
Double_t minMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(minBinBkgHi);
Double_t maxMassBkgHiBin=hMassD->GetYaxis()->GetBinLowEdge(maxBinBkgHi)+hMassD->GetYaxis()->GetBinWidth(maxBinBkgHi);
// printf("BKG Fit Limits = %f %f && %f %f\n",minMassBkgLowBin,maxMassBkgLowBin,minMassBkgHiBin,maxMassBkgHiBin);
delete hMptSig;
}
+//____________________________________________________________________________
+Bool_t AliVertexingHFUtils::CheckT0TriggerFired(AliAODEvent* aodEv){
+ // check if T0VTX trigger was fired, based on a workaround suggested by Alla
+ const Double32_t *mean = aodEv->GetT0TOF();
+ if(mean && mean[0]<9999.) return kTRUE;
+ else return kFALSE;
+}
+//____________________________________________________________________________
+Double_t AliVertexingHFUtils::GetTrueImpactParameterDzero(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
+ // true impact parameter calculation for Dzero
+
+ if(!partD || !arrayMC || !mcHeader) return 99999.;
+ Int_t code=partD->GetPdgCode();
+ if(TMath::Abs(code)!=421) return 99999.;
+
+ Double_t vtxTrue[3];
+ mcHeader->GetVertex(vtxTrue);
+ Double_t origD[3];
+ partD->XvYvZv(origD);
+ Short_t charge=partD->Charge();
+ Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
+ for(Int_t iDau=0; iDau<2; iDau++){
+ pXdauTrue[iDau]=0.;
+ pYdauTrue[iDau]=0.;
+ pZdauTrue[iDau]=0.;
+ }
+
+ Int_t nDau=partD->GetNDaughters();
+ Int_t labelFirstDau = partD->GetDaughter(0);
+ if(nDau==2){
+ for(Int_t iDau=0; iDau<2; iDau++){
+ Int_t ind = labelFirstDau+iDau;
+ AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
+ if(!part){
+ printf("Daughter particle not found in MC array");
+ return 99999.;
+ }
+ pXdauTrue[iDau]=part->Px();
+ pYdauTrue[iDau]=part->Py();
+ pZdauTrue[iDau]=part->Pz();
+ }
+ }else{
+ return 99999.;
+ }
+
+ Double_t d0dummy[2]={0.,0.};
+ AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
+ return aodDvsMC.ImpParXY();
+
+}
+//____________________________________________________________________________
+Double_t AliVertexingHFUtils::GetTrueImpactParameterDplus(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD) {
+ // true impact parameter calculation for Dplus
+
+ if(!partD || !arrayMC || !mcHeader) return 99999.;
+ Int_t code=partD->GetPdgCode();
+ if(TMath::Abs(code)!=411) return 99999.;
+
+ Double_t vtxTrue[3];
+ mcHeader->GetVertex(vtxTrue);
+ Double_t origD[3];
+ partD->XvYvZv(origD);
+ Short_t charge=partD->Charge();
+ Double_t pXdauTrue[3],pYdauTrue[3],pZdauTrue[3];
+ for(Int_t iDau=0; iDau<3; iDau++){
+ pXdauTrue[iDau]=0.;
+ pYdauTrue[iDau]=0.;
+ pZdauTrue[iDau]=0.;
+ }
+
+ Int_t nDau=partD->GetNDaughters();
+ Int_t labelFirstDau = partD->GetDaughter(0);
+ if(nDau==3){
+ for(Int_t iDau=0; iDau<3; iDau++){
+ Int_t ind = labelFirstDau+iDau;
+ AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
+ if(!part){
+ printf("Daughter particle not found in MC array");
+ return 99999.;
+ }
+ pXdauTrue[iDau]=part->Px();
+ pYdauTrue[iDau]=part->Py();
+ pZdauTrue[iDau]=part->Pz();
+ }
+ }else if(nDau==2){
+ Int_t theDau=0;
+ for(Int_t iDau=0; iDau<2; iDau++){
+ Int_t ind = labelFirstDau+iDau;
+ AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
+ if(!part){
+ printf("Daughter particle not found in MC array");
+ return 99999.;
+ }
+ Int_t pdgCode=TMath::Abs(part->GetPdgCode());
+ if(pdgCode==211 || pdgCode==321){
+ pXdauTrue[theDau]=part->Px();
+ pYdauTrue[theDau]=part->Py();
+ pZdauTrue[theDau]=part->Pz();
+ ++theDau;
+ }else{
+ Int_t nDauRes=part->GetNDaughters();
+ if(nDauRes==2){
+ Int_t labelFirstDauRes = part->GetDaughter(0);
+ for(Int_t iDauRes=0; iDauRes<2; iDauRes++){
+ Int_t indDR = labelFirstDauRes+iDauRes;
+ AliAODMCParticle* partDR = dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDR));
+ if(!partDR){
+ printf("Daughter particle not found in MC array");
+ return 99999.;
+ }
+
+ Int_t pdgCodeDR=TMath::Abs(partDR->GetPdgCode());
+ if(pdgCodeDR==211 || pdgCodeDR==321){
+ pXdauTrue[theDau]=partDR->Px();
+ pYdauTrue[theDau]=partDR->Py();
+ pZdauTrue[theDau]=partDR->Pz();
+ ++theDau;
+ }
+ }
+ }
+ }
+ }
+ if(theDau!=3){
+ printf("Wrong number of decay prongs");
+ return 99999.;
+ }
+ }
+
+ Double_t d0dummy[3]={0.,0.,0.};
+ AliAODRecoDecayHF aodDvsMC(vtxTrue,origD,3,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
+ return aodDvsMC.ImpParXY();
+
+}
+
+
+
//____________________________________________________________________________
Double_t AliVertexingHFUtils::GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult) {
//
//
if(TMath::Abs(vtxZ)>10.0){
- printf("ERROR: Z vertex out of range for correction of multiplicity\n");
+ // printf("ERROR: Z vertex out of range for correction of multiplicity\n");
return uncorrectedNacc;
}
return correctedNacc;
}
+//______________________________________________________________________
+TString AliVertexingHFUtils::GetGenerator(Int_t label, AliAODMCHeader* header){
+ // get the name of the generator that produced a given particle
+
+ Int_t nsumpart=0;
+ TList *lh=header->GetCocktailHeaders();
+ Int_t nh=lh->GetEntries();
+ for(Int_t i=0;i<nh;i++){
+ AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
+ TString genname=gh->GetName();
+ Int_t npart=gh->NProduced();
+ if(label>=nsumpart && label<(nsumpart+npart)) return genname;
+ nsumpart+=npart;
+ }
+ TString empty="";
+ return empty;
+}
+//_____________________________________________________________________
+void AliVertexingHFUtils::GetTrackPrimaryGenerator(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
+
+ // method to check if a track comes from a given generator
+
+ Int_t lab=TMath::Abs(track->GetLabel());
+ nameGen=GetGenerator(lab,header);
+
+ // Int_t countControl=0;
+
+ while(nameGen.IsWhitespace()){
+ AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
+ if(!mcpart){
+ printf("AliVertexingHFUtils::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
+ break;
+ }
+ Int_t mother = mcpart->GetMother();
+ if(mother<0){
+ printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
+ break;
+ }
+ lab=mother;
+ nameGen=GetGenerator(mother,header);
+ // countControl++;
+ // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
+ // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
+ // break;
+ // }
+ }
+
+ return;
+}
+//----------------------------------------------------------------------
+Bool_t AliVertexingHFUtils::IsTrackInjected(AliAODTrack *track,AliAODMCHeader *header,TClonesArray *arrayMC){
+ // method to check if a track comes from the signal event or from the underlying Hijing event
+ TString nameGen;
+ GetTrackPrimaryGenerator(track,header,arrayMC,nameGen);
+
+ if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;
+
+ return kTRUE;
+}
+//____________________________________________________________________________
+Bool_t AliVertexingHFUtils::IsCandidateInjected(AliAODRecoDecayHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
+ // method to check if a D meson candidate comes from the signal event or from the underlying Hijing event
+
+ Int_t nprongs=cand->GetNProngs();
+ for(Int_t i=0;i<nprongs;i++){
+ AliAODTrack *daugh=(AliAODTrack*)cand->GetDaughter(i);
+ if(IsTrackInjected(daugh,header,arrayMC)) return kTRUE;
+ }
+ return kFALSE;
+}
+//____________________________________________________________________________
+Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascadeHF *cand, AliAODMCHeader *header,TClonesArray *arrayMC){
+ // method to check if a cascade candidate comes from the signal event or from the underlying Hijing event
+
+ AliAODTrack* bach = cand->GetBachelor();
+ if(IsTrackInjected(bach, header, arrayMC)) {
+ AliDebug(2, "Bachelor is injected, the whole candidate is then injected");
+ return kTRUE;
+ }
+ AliAODv0* v0 = cand->Getv0();
+ Int_t nprongs = v0->GetNProngs();
+ for(Int_t i = 0; i < nprongs; i++){
+ AliAODTrack *daugh = (AliAODTrack*)v0->GetDaughter(i);
+ if(IsTrackInjected(daugh,header,arrayMC)) {
+ AliDebug(2, Form("V0 daughter number %d is injected, the whole candidate is then injected", i));
+ return kTRUE;
+ }
+ }
+ return kFALSE;
+}
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckOrigin(AliStack* stack, TParticle *mcPart, Bool_t searchUpToQuark){
+ // checking whether the mother of the particles come from a charm or a bottom quark
+
+ Int_t pdgGranma = 0;
+ Int_t mother = 0;
+ mother = mcPart->GetFirstMother();
+ Int_t istep = 0;
+ Int_t abspdgGranma =0;
+ Bool_t isFromB=kFALSE;
+ Bool_t isQuarkFound=kFALSE;
+ while (mother >0 ){
+ istep++;
+ TParticle* mcGranma = stack->Particle(mother);
+ if (mcGranma){
+ pdgGranma = mcGranma->GetPdgCode();
+ abspdgGranma = TMath::Abs(pdgGranma);
+ if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
+ isFromB=kTRUE;
+ }
+ if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
+ mother = mcGranma->GetFirstMother();
+ }else{
+ printf("CheckOrigin: Failed casting the mother particle!");
+ break;
+ }
+ }
+ if(searchUpToQuark && !isQuarkFound) return 0;
+ if(isFromB) return 5;
+ else return 4;
+
+}
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){
+ // checking whether the mother of the particles come from a charm or a bottom quark
+
+ Int_t pdgGranma = 0;
+ Int_t mother = 0;
+ mother = mcPart->GetMother();
+ Int_t istep = 0;
+ Int_t abspdgGranma =0;
+ Bool_t isFromB=kFALSE;
+ Bool_t isQuarkFound=kFALSE;
+ while (mother >0 ){
+ istep++;
+ AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
+ if (mcGranma){
+ pdgGranma = mcGranma->GetPdgCode();
+ abspdgGranma = TMath::Abs(pdgGranma);
+ if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
+ isFromB=kTRUE;
+ }
+ if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
+ mother = mcGranma->GetMother();
+ }else{
+ printf("AliVertexingHFUtils::CheckOrigin: Failed casting the mother particle!");
+ break;
+ }
+ }
+ if(searchUpToQuark && !isQuarkFound) return 0;
+ if(isFromB) return 5;
+ else return 4;
+
+}
+
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckD0Decay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
+ // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
+
+ if(label<0) return -1;
+ TParticle* mcPart = stack->Particle(label);
+ Int_t pdgD=mcPart->GetPdgCode();
+ if(TMath::Abs(pdgD)!=421) return -1;
+
+ Int_t nDau=mcPart->GetNDaughters();
+
+ if(nDau==2){
+ Int_t daughter0 = mcPart->GetDaughter(0);
+ Int_t daughter1 = mcPart->GetDaughter(1);
+ TParticle* mcPartDaughter0 = stack->Particle(daughter0);
+ TParticle* mcPartDaughter1 = stack->Particle(daughter1);
+ if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
+ arrayDauLab[0]=daughter0;
+ arrayDauLab[1]=daughter1;
+ Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
+ Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
+ if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
+ !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
+ return -1;
+ }
+ if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
+ if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
+ if((pdgCode0*pdgCode1)>0) return -1;
+ Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
+ Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
+ Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
+ 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;
+ }
+
+ if(nDau==3 || nDau==4){
+ Int_t nKaons=0;
+ Int_t nPions=0;
+ Double_t sumPxDau=0.;
+ Double_t sumPyDau=0.;
+ Double_t sumPzDau=0.;
+ Int_t nFoundKpi=0;
+ Int_t labelFirstDau = mcPart->GetDaughter(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)==321){
+ if(pdgD>0 && pdgdau>0) return -1;
+ if(pdgD<0 && pdgdau<0) return -1;
+ nKaons++;
+ sumPxDau+=dau->Px();
+ sumPyDau+=dau->Py();
+ sumPzDau+=dau->Pz();
+ arrayDauLab[nFoundKpi++]=indDau;
+ if(nFoundKpi>4) return -1;
+ }else if(TMath::Abs(pdgdau)==211){
+ nPions++;
+ sumPxDau+=dau->Px();
+ sumPyDau+=dau->Py();
+ sumPzDau+=dau->Pz();
+ arrayDauLab[nFoundKpi++]=indDau;
+ if(nFoundKpi>4) return -1;
+ }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
+ 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){
+ if(pdgD>0 && pdgresdau>0) return -1;
+ if(pdgD<0 && pdgresdau<0) return -1;
+ nKaons++;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ arrayDauLab[nFoundKpi++]=indResDau;
+ if(nFoundKpi>4) return -1;
+ }
+ if(TMath::Abs(pdgresdau)==211){
+ nPions++;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ arrayDauLab[nFoundKpi++]=indResDau;
+ if(nFoundKpi>4) return -1;
+ }
+ }
+ }else{
+ return -1;
+ }
+ }
+ if(nPions!=3) return -1;
+ if(nKaons!=1) return -1;
+ if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
+ if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
+ if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
+ return 2;
+ }
+ return -1;
+}
+
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
+ // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases
+
+ Int_t pdgD=mcPart->GetPdgCode();
+ if(TMath::Abs(pdgD)!=421) return -1;
+
+ Int_t nDau=mcPart->GetNDaughters();
+
+ if(nDau==2){
+ Int_t daughter0 = mcPart->GetDaughter(0);
+ Int_t daughter1 = mcPart->GetDaughter(1);
+ AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter0));
+ AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(daughter1));
+ if(!mcPartDaughter0 || !mcPartDaughter1) return -1;
+ arrayDauLab[0]=daughter0;
+ arrayDauLab[1]=daughter1;
+ Int_t pdgCode0=mcPartDaughter0->GetPdgCode();
+ Int_t pdgCode1=mcPartDaughter1->GetPdgCode();
+ if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) &&
+ !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){
+ return -1;
+ }
+ if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1;
+ if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1;
+ if((pdgCode0*pdgCode1)>0) return -1;
+ Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px();
+ Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py();
+ Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
+ 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;
+ }
+
+ if(nDau==3 || nDau==4){
+ Int_t nKaons=0;
+ Int_t nPions=0;
+ Double_t sumPxDau=0.;
+ Double_t sumPyDau=0.;
+ Double_t sumPzDau=0.;
+ Int_t nFoundKpi=0;
+ Int_t labelFirstDau = mcPart->GetDaughter(0);
+ for(Int_t iDau=0; iDau<nDau; iDau++){
+ Int_t indDau = labelFirstDau+iDau;
+ if(indDau<0) return -1;
+ AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
+ if(!dau) return -1;
+ Int_t pdgdau=dau->GetPdgCode();
+ if(TMath::Abs(pdgdau)==321){
+ if(pdgD>0 && pdgdau>0) return -1;
+ if(pdgD<0 && pdgdau<0) return -1;
+ nKaons++;
+ sumPxDau+=dau->Px();
+ sumPyDau+=dau->Py();
+ sumPzDau+=dau->Pz();
+ arrayDauLab[nFoundKpi++]=indDau;
+ if(nFoundKpi>4) return -1;
+ }else if(TMath::Abs(pdgdau)==211){
+ nPions++;
+ sumPxDau+=dau->Px();
+ sumPyDau+=dau->Py();
+ sumPzDau+=dau->Pz();
+ arrayDauLab[nFoundKpi++]=indDau;
+ if(nFoundKpi>4) return -1;
+ }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){
+ 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;
+ AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
+ if(!resdau) return -1;
+ Int_t pdgresdau=resdau->GetPdgCode();
+ if(TMath::Abs(pdgresdau)==321){
+ if(pdgD>0 && pdgresdau>0) return -1;
+ if(pdgD<0 && pdgresdau<0) return -1;
+ nKaons++;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ arrayDauLab[nFoundKpi++]=indResDau;
+ if(nFoundKpi>4) return -1;
+ }
+ if(TMath::Abs(pdgresdau)==211){
+ nPions++;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ arrayDauLab[nFoundKpi++]=indResDau;
+ if(nFoundKpi>4) return -1;
+ }
+ }
+ }else{
+ return -1;
+ }
+ }
+ if(nPions!=3) return -1;
+ if(nKaons!=1) return -1;
+ if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1;
+ if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1;
+ if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1;
+ return 2;
+ }
+ return -1;
+}
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckDplusDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
+ // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
+
+ 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;
+
+ 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){
+ if(pdgD*pdgdau>0) return -1;
+ nKaons++;
+ sumPxDau+=dau->Px();
+ sumPyDau+=dau->Py();
+ sumPzDau+=dau->Pz();
+ arrayDauLab[nFoundKpi++]=indDau;
+ if(nFoundKpi>3) return -1;
+ }else if(TMath::Abs(pdgdau)==211){
+ if(pdgD*pdgdau<0) return -1;
+ nPions++;
+ sumPxDau+=dau->Px();
+ sumPyDau+=dau->Py();
+ sumPzDau+=dau->Pz();
+ arrayDauLab[nFoundKpi++]=indDau;
+ if(nFoundKpi>3) return -1;
+ }else if(TMath::Abs(pdgdau)==313){
+ 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){
+ if(pdgD*pdgresdau>0) return -1;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ nKaons++;
+ arrayDauLab[nFoundKpi++]=indResDau;
+ if(nFoundKpi>3) return -1;
+ }
+ if(TMath::Abs(pdgresdau)==211){
+ if(pdgD*pdgresdau<0) return -1;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ nPions++;
+ arrayDauLab[nFoundKpi++]=indResDau;
+ 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(nDau==3) return 1;
+ else if(nDau==2) return 2;
+ }
+
+ return -1;
+
+}
+
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
+ // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases
+
+ 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;
+
+ if(nDau==3 || nDau==2){
+ for(Int_t iDau=0; iDau<nDau; iDau++){
+ Int_t indDau = labelFirstDau+iDau;
+ if(indDau<0) return -1;
+ AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
+ if(!dau) return -1;
+ Int_t pdgdau=dau->GetPdgCode();
+ if(TMath::Abs(pdgdau)==321){
+ if(pdgD*pdgdau>0) return -1;
+ nKaons++;
+ sumPxDau+=dau->Px();
+ sumPyDau+=dau->Py();
+ sumPzDau+=dau->Pz();
+ arrayDauLab[nFoundKpi++]=indDau;
+ if(nFoundKpi>3) return -1;
+ }else if(TMath::Abs(pdgdau)==211){
+ if(pdgD*pdgdau<0) return -1;
+ nPions++;
+ sumPxDau+=dau->Px();
+ sumPyDau+=dau->Py();
+ sumPzDau+=dau->Pz();
+ arrayDauLab[nFoundKpi++]=indDau;
+ if(nFoundKpi>3) return -1;
+ }else if(TMath::Abs(pdgdau)==313){
+ 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;
+ AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
+ if(!resdau) return -1;
+ Int_t pdgresdau=resdau->GetPdgCode();
+ if(TMath::Abs(pdgresdau)==321){
+ if(pdgD*pdgresdau>0) return -1;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ nKaons++;
+ arrayDauLab[nFoundKpi++]=indResDau;
+ if(nFoundKpi>3) return -1;
+ }
+ if(TMath::Abs(pdgresdau)==211){
+ if(pdgD*pdgresdau<0) return -1;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ nPions++;
+ arrayDauLab[nFoundKpi++]=indResDau;
+ 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(nDau==3) return 1;
+ else if(nDau==2) return 2;
+ }
+
+ 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
+
+ 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 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::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
+
+ 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 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;
+ AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(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;
+ AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(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::CheckDstarDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){
+ // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases
+
+ if(label<0) return -1;
+ TParticle* mcPart = stack->Particle(label);
+ Int_t pdgD=mcPart->GetPdgCode();
+ if(TMath::Abs(pdgD)!=413) return -1;
+
+ Int_t nDau=mcPart->GetNDaughters();
+ if(nDau!=2) return -1;
+
+ 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;
+
+ 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)==421){
+ 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){
+ if(pdgD*pdgresdau>0) return -1;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ nKaons++;
+ arrayDauLab[nFoundKpi++]=indResDau;
+ if(nFoundKpi>3) return -1;
+ }
+ if(TMath::Abs(pdgresdau)==211){
+ if(pdgD*pdgresdau<0) return -1;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ nPions++;
+ arrayDauLab[nFoundKpi++]=indResDau;
+ if(nFoundKpi>3) return -1;
+ }
+ }
+ }else if(TMath::Abs(pdgdau)==211){
+ if(pdgD*pdgdau<0) return -1;
+ nPions++;
+ sumPxDau+=dau->Px();
+ sumPyDau+=dau->Py();
+ sumPzDau+=dau->Pz();
+ arrayDauLab[nFoundKpi++]=indDau;
+ if(nFoundKpi>3) 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;
+ return 1;
+
+}
+
+//____________________________________________________________________________
+Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){
+ // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases
+
+ Int_t pdgD=mcPart->GetPdgCode();
+ if(TMath::Abs(pdgD)!=413) return -1;
+
+ Int_t nDau=mcPart->GetNDaughters();
+ if(nDau!=2) return -1;
+
+ 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;
+
+ for(Int_t iDau=0; iDau<nDau; iDau++){
+ Int_t indDau = labelFirstDau+iDau;
+ if(indDau<0) return -1;
+ AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indDau));
+ if(!dau) return -1;
+ Int_t pdgdau=dau->GetPdgCode();
+ if(TMath::Abs(pdgdau)==421){
+ 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;
+ AliAODMCParticle* resdau=dynamic_cast<AliAODMCParticle*>(arrayMC->At(indResDau));
+ if(!resdau) return -1;
+ Int_t pdgresdau=resdau->GetPdgCode();
+ if(TMath::Abs(pdgresdau)==321){
+ if(pdgD*pdgresdau>0) return -1;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ nKaons++;
+ arrayDauLab[nFoundKpi++]=indResDau;
+ if(nFoundKpi>3) return -1;
+ }
+ if(TMath::Abs(pdgresdau)==211){
+ if(pdgD*pdgresdau<0) return -1;
+ sumPxDau+=resdau->Px();
+ sumPyDau+=resdau->Py();
+ sumPzDau+=resdau->Pz();
+ nPions++;
+ arrayDauLab[nFoundKpi++]=indResDau;
+ if(nFoundKpi>3) return -1;
+ }
+ }
+ }else if(TMath::Abs(pdgdau)==211){
+ if(pdgD*pdgdau<0) return -1;
+ nPions++;
+ sumPxDau+=dau->Px();
+ sumPyDau+=dau->Py();
+ sumPzDau+=dau->Pz();
+ arrayDauLab[nFoundKpi++]=indDau;
+ if(nFoundKpi>3) 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;
+ return 1;
+
+}
+
+//____________________________________________________________________________
+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;
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff