[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / AliRsnLoopEffPair.cxx

//
// Class AliRsnLoopEffPair
//
// Inherits from basic AliRsnLoopEff for efficiency,
// and computed efficiencies for single-tracks
//
// author: Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
//

#include <Riostream.h>

#include "AliStack.h"
#include "AliGenEventHeader.h"
#include "AliAODMCHeader.h"

#include "AliRsnEvent.h"
#include "AliRsnCutManager.h"
#include "AliRsnDaughterDef.h"
#include "AliRsnPairDef.h"
#include "AliRsnLoopEffPair.h"

ClassImp(AliRsnLoopEffPair)

//_____________________________________________________________________________
AliRsnLoopEffPair::AliRsnLoopEffPair(const char *name, AliRsnPairDef *def) :
   AliRsnLoopEff(name, 1),
   fDef(def),
   fMother()
{
//
// Default constructor.
// Do not repeat 'DefineOutput' since it is done in base class and we don't add new ones.
//
}

//_____________________________________________________________________________
AliRsnLoopEffPair::AliRsnLoopEffPair(const AliRsnLoopEffPair &copy) :
   AliRsnLoopEff(copy),
   fDef(copy.fDef),
   fMother(copy.fMother)
{
//
// Copy constrtuctor.
//
}

//_____________________________________________________________________________
AliRsnLoopEffPair &AliRsnLoopEffPair::operator=(const AliRsnLoopEffPair &copy)
{
//
// Assignment operator.
// Owned data members are meaningless for this operator.
//

   AliRsnLoopEff::operator=(copy);
   if (this == &copy)
      return *this;
   fDef = copy.fDef;

   return (*this);
}

//__________________________________________________________________________________________________
Bool_t AliRsnLoopEffPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
{
//
// Calls the appropriate assignment method
//

   // setup pointers
   fMother.SetDaughter(0, &fDaughter[0]);
   fMother.SetDaughter(1, &fDaughter[1]);
   fMother.SetRefEvent(rsnEvent);
   fDaughter[0].SetOwnerEvent(rsnEvent);
   fDaughter[1].SetOwnerEvent(rsnEvent);

   if (rsnEvent->IsESD())
      return AssignMotherAndDaughtersESD(rsnEvent, ipart);
   else if (rsnEvent->IsAOD())
      return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
   else {
      AliError("Unrecognized input event");
      return kFALSE;
   }
}

//__________________________________________________________________________________________________
Bool_t AliRsnLoopEffPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
{
//
// Gets a particle in the MC event and try to assign it to the mother.
// If it has two daughters, retrieve them and assign also them.
// NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
//       for ESD and AOD, if available
// ---
// Implementation for ESD inputs
//

   AliMCEvent    *mc      = rsnEvent->GetRefMCESD();
   AliStack      *stack   = mc->Stack();
   AliMCParticle *mother  = (AliMCParticle *)mc->GetTrack(ipart);
   TParticle     *motherP = mother->Particle();
   Int_t          ntracks = stack->GetNtrack();

   // check PDG code and exit if it is wrong
   if (TMath::Abs(motherP->GetPdgCode()) != fDef->GetMotherPDG()) return kFALSE;

   // check number of daughters and exit if it is not 2
   if (motherP->GetNDaughters() < 2) return kFALSE;

   // check distance from primary vertex
   TLorentzVector vprod;
   motherP->ProductionVertex(vprod);
   if (!CheckDistanceFromPV(vprod.X(), vprod.Y(), vprod.Z())) {
      AliDebugClass(1, "Distant production vertex");
      return kFALSE;
   }

   // get the daughters and check their PDG code and charge:
   // if they match one of the pair daughter definitions,
   // assign them as MC reference of the 'fDaughter' objects
   fDaughter[0].Reset();
   fDaughter[1].Reset();
   Int_t index[2] = {motherP->GetFirstDaughter(), motherP->GetLastDaughter()};
   Int_t i, pdg;
   Short_t charge;
   AliMCParticle *daughter = 0x0;
   for (i = 0; i < 2; i++) {
      // check index for stack
      if (index[i] < 0 || index[i] > ntracks) {
         AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
         return kFALSE;
      }
      // get daughter and its PDG and charge
      daughter = (AliMCParticle *)mc->GetTrack(index[i]);
      pdg      = TMath::Abs(daughter->Particle()->GetPdgCode());
      charge   = (Short_t)(daughter->Particle()->GetPDG()->Charge() / 3);
      // check if it matches one definition
      if (fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge)) {
         fDaughter[0].SetGood();
         fDaughter[0].SetRefMC(daughter);
         fDaughter[0].SetLabel(index[i]);
      } else if (fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge)) {
         fDaughter[1].SetGood();
         fDaughter[1].SetRefMC(daughter);
         fDaughter[1].SetLabel(index[i]);
      }
   }

   // return success if both daughters were assigned
   if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
      return kTRUE;
   } else {
      fDaughter[0].Reset();
      fDaughter[1].Reset();
      return kFALSE;
   }
}

//__________________________________________________________________________________________________
Bool_t AliRsnLoopEffPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
{
//
// Gets a particle in the MC event and try to assign it to the mother.
// If it has two daughters, retrieve them and assign also them.
// NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
//       for ESD and AOD, if available
// ---
// Implementation for AOD inputs
//

   AliAODEvent      *aod     = rsnEvent->GetRefAOD();
   TClonesArray     *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
   AliAODMCParticle *mother  = (AliAODMCParticle *)listAOD->At(ipart);
   Int_t             ntracks = listAOD->GetEntries();

   // check PDG code and exit if it is wrong
   if (TMath::Abs(mother->GetPdgCode()) != fDef->GetMotherPDG()) return kFALSE;

   // check number of daughters and exit if it is not 2
   if (mother->GetNDaughters() < 2) return kFALSE;

   // check distance from primary vertex
   Double_t vprod[3] = {(Double_t)mother->Xv(), (Double_t)mother->Yv(), (Double_t)mother->Zv()};
   Double_t dv = DistanceFromPV(vprod[0], vprod[1], vprod[2]);
   if (dv > fMaxDistPV) {
      AliDebugClass(1, "Distant production vertex");
      return kFALSE;
   }

   // get the daughters and check their PDG code and charge:
   // if they match one of the pair daughter definitions,
   // assign them as MC reference of the 'fDaughter' objects
   fDaughter[0].Reset();
   fDaughter[1].Reset();
   Int_t index[2] = {(Int_t)mother->GetDaughter(0), (Int_t)mother->GetDaughter(1)};
   Int_t i, pdg;
   Short_t charge;
   AliAODMCParticle *daughter = 0x0;
   for (i = 0; i < 2; i++) {
      // check index for stack
      if (index[i] < 0 || index[i] > ntracks) {
         AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
         return kFALSE;
      }
      // get daughter and its PDG and charge
      daughter = (AliAODMCParticle *)listAOD->At(index[i]);
      pdg      = TMath::Abs(daughter->GetPdgCode());
      charge   = (Short_t)(daughter->Charge() / 3);
      // check if it matches one definition
      if (fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge)) {
         fDaughter[0].SetGood();
         fDaughter[0].SetRefMC(daughter);
         fDaughter[0].SetLabel(index[i]);
      } else if (fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge)) {
         fDaughter[1].SetGood();
         fDaughter[1].SetRefMC(daughter);
         fDaughter[1].SetLabel(index[i]);
      }
   }

   // return success if both daughters were assigned
   if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
      return kTRUE;
   } else {
      fDaughter[0].Reset();
      fDaughter[1].Reset();
      return kFALSE;
   }
}

//_____________________________________________________________________________
Int_t AliRsnLoopEffPair::DoLoop(AliRsnEvent *rsn, AliRsnDaughterSelector *, AliRsnEvent *, AliRsnDaughterSelector *)
{
//
// Loop on event and fill containers
//

   // check event cuts
   if (!OkEvent(rsn)) return 0;

   // retrieve output
   fOutput = (AliRsnListOutput *)fOutputs[0];

   // check presence of MC reference
   if (!rsn->GetRefMC()) {
      AliError("Need a MC to compute efficiency");
      return 0;
   }

   // check presence of event
   if (!rsn->GetRef()) {
      AliError("Need an event to compute efficiency");
      return 0;
   }

   // check event type:
   // must be ESD or AOD, and then use a bool to know in the rest
   if (!rsn->IsESD() && !rsn->IsAOD()) {
      AliError("Need to process ESD or AOD input");
      return 0;
   }

   // retrieve the MC primary vertex position
   // and do some additional coherence checks
   Int_t i, npart = 0;
   if (rsn->IsESD()) {
      TArrayF mcVertex(3);
      AliGenEventHeader *genEH = rsn->GetRefMCESD()->GenEventHeader();
      genEH->PrimaryVertex(mcVertex);
      for (i = 0; i < 3; i++) fVertex[i] = (Double_t)mcVertex[i];
      npart = rsn->GetRefMCESD()->GetNumberOfTracks();
   } else {
      for (i = 0; i < 3; i++) fVertex[i] = 0.0;
      AliAODEvent *aod = rsn->GetRefMCAOD();
      TClonesArray *listAOD = (TClonesArray *)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
      if (listAOD) npart = listAOD->GetEntries();
      AliAODMCHeader *mcH = static_cast<AliAODMCHeader *>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
      if (mcH) mcH->GetVertex(fVertex);
   }

   // check number of particles
   if (!npart) {
      AliInfo("Empty event");
      return 0;
   }

   // utility variables
   Int_t ipart, istep, count = 0, nsteps = fSteps.GetEntries();
   Int_t ntracks = rsn->GetAbsoluteSum();
   AliRsnDaughter check;

   // loop over particles
   for (ipart = 0; ipart < npart; ipart++) {
      // check i-th particle
      if (!AssignMotherAndDaughters(rsn, ipart)) continue;
      // if assignment was successful, for first step, use MC info
      fDaughter[0].SetRef(fDaughter[0].GetRefMC());
      fDaughter[1].SetRef(fDaughter[1].GetRefMC());
      fMother.SetRefEvent(rsn);
      fMother.ComputeSum(fDef->GetDef1().GetMass(), fDef->GetDef2().GetMass(), fDef->GetMotherMass());
      fOutput->Fill(&fMother, 0);
      count++;
      // for each further step, try to find two tracks which pass the related cuts
      for (istep = 0; istep < nsteps; istep++) {
         AliRsnCutManager *cuts = (AliRsnCutManager *)fSteps[istep];
         fDaughter[0].SetBad();
         fDaughter[1].SetBad();
         for (i = 0; i < ntracks; i++) {
            rsn->SetDaughter(check, i);
            if (!cuts->PassCommonDaughterCuts(&check)) continue;
            if (TMath::Abs(check.GetLabel()) == fDaughter[0].GetLabel()) {
               if (!cuts->PassDaughter1Cuts(&check)) continue;
               fDaughter[0].SetRef(check.GetRef());
               fDaughter[0].SetGood();
            } else if (TMath::Abs(check.GetLabel()) == fDaughter[1].GetLabel()) {
               if (!cuts->PassDaughter2Cuts(&check)) continue;
               fDaughter[1].SetRef(check.GetRef());
               fDaughter[1].SetGood();
            }
            if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
               fMother.ComputeSum(fDef->GetDef1().GetMass(), fDef->GetDef2().GetMass(), fDef->GetMotherMass());
               if (cuts->PassMotherCuts(&fMother)) {
                  fOutput->Fill(&fMother, 1 + istep);
                  break;
               }
            }
         }
      }
   }

   return count;
}
Commit	Line	Data
c865cb1d	1	//
	2	// Class AliRsnLoopEffPair
	3	//
	4	// Inherits from basic AliRsnLoopEff for efficiency,
	5	// and computed efficiencies for single-tracks
	6	//
	7	// author: Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
	8	//
	9
f34f960b	10	#include <Riostream.h>
f34f960b	11
c865cb1d	12	#include "AliStack.h"
f34f960b	13	#include "AliGenEventHeader.h"
f34f960b	14	#include "AliAODMCHeader.h"
c865cb1d	15
e187bd70	16	#include "AliRsnEvent.h"
c865cb1d	17	#include "AliRsnCutManager.h"
	18	#include "AliRsnDaughterDef.h"
	19	#include "AliRsnPairDef.h"
	20	#include "AliRsnLoopEffPair.h"
	21
	22	ClassImp(AliRsnLoopEffPair)
	23
	24	//_____________________________________________________________________________
	25	AliRsnLoopEffPair::AliRsnLoopEffPair(const char name, AliRsnPairDef def) :
f34f960b	26	AliRsnLoopEff(name, 1),
c865cb1d	27	fDef(def),
	28	fMother()
	29	{
	30	//
	31	// Default constructor.
	32	// Do not repeat 'DefineOutput' since it is done in base class and we don't add new ones.
	33	//
	34	}
	35
	36	//_____________________________________________________________________________
61f275d1	37	AliRsnLoopEffPair::AliRsnLoopEffPair(const AliRsnLoopEffPair &copy) :
c865cb1d	38	AliRsnLoopEff(copy),
	39	fDef(copy.fDef),
	40	fMother(copy.fMother)
	41	{
	42	//
	43	// Copy constrtuctor.
	44	//
	45	}
	46
	47	//_____________________________________________________________________________
61f275d1	48	AliRsnLoopEffPair &AliRsnLoopEffPair::operator=(const AliRsnLoopEffPair &copy)
c865cb1d	49	{
	50	//
	51	// Assignment operator.
	52	// Owned data members are meaningless for this operator.
	53	//
	54
	55	AliRsnLoopEff::operator=(copy);
e6f3a909	56	if (this == &copy)
61f275d1	57	return *this;
c865cb1d	58	fDef = copy.fDef;
	59
	60	return (*this);
	61	}
	62
f34f960b	63	//__________________________________________________________________________________________________
	64	Bool_t AliRsnLoopEffPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
	65	{
	66	//
	67	// Calls the appropriate assignment method
	68	//
	69
	70	// setup pointers
	71	fMother.SetDaughter(0, &fDaughter[0]);
	72	fMother.SetDaughter(1, &fDaughter[1]);
	73	fMother.SetRefEvent(rsnEvent);
	74	fDaughter[0].SetOwnerEvent(rsnEvent);
	75	fDaughter[1].SetOwnerEvent(rsnEvent);
	76
	77	if (rsnEvent->IsESD())
	78	return AssignMotherAndDaughtersESD(rsnEvent, ipart);
	79	else if (rsnEvent->IsAOD())
	80	return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
	81	else {
	82	AliError("Unrecognized input event");
	83	return kFALSE;
	84	}
	85	}
	86
	87	//__________________________________________________________________________________________________
	88	Bool_t AliRsnLoopEffPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
	89	{
	90	//
	91	// Gets a particle in the MC event and try to assign it to the mother.
	92	// If it has two daughters, retrieve them and assign also them.
	93	// NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
	94	// for ESD and AOD, if available
	95	// ---
	96	// Implementation for ESD inputs
	97	//
	98
	99	AliMCEvent *mc = rsnEvent->GetRefMCESD();
	100	AliStack *stack = mc->Stack();
61f275d1	101	AliMCParticle mother = (AliMCParticle )mc->GetTrack(ipart);
f34f960b	102	TParticle *motherP = mother->Particle();
f34f960b	103	Int_t ntracks = stack->GetNtrack();
61f275d1	104
f34f960b	105	// check PDG code and exit if it is wrong
f34f960b	106	if (TMath::Abs(motherP->GetPdgCode()) != fDef->GetMotherPDG()) return kFALSE;
61f275d1	107
f34f960b	108	// check number of daughters and exit if it is not 2
f34f960b	109	if (motherP->GetNDaughters() < 2) return kFALSE;
61f275d1	110
f34f960b	111	// check distance from primary vertex
	112	TLorentzVector vprod;
	113	motherP->ProductionVertex(vprod);
	114	if (!CheckDistanceFromPV(vprod.X(), vprod.Y(), vprod.Z())) {
	115	AliDebugClass(1, "Distant production vertex");
	116	return kFALSE;
	117	}
61f275d1	118
f34f960b	119	// get the daughters and check their PDG code and charge:
61f275d1	120	// if they match one of the pair daughter definitions,
f34f960b	121	// assign them as MC reference of the 'fDaughter' objects
	122	fDaughter[0].Reset();
	123	fDaughter[1].Reset();
	124	Int_t index[2] = {motherP->GetFirstDaughter(), motherP->GetLastDaughter()};
	125	Int_t i, pdg;
	126	Short_t charge;
	127	AliMCParticle *daughter = 0x0;
	128	for (i = 0; i < 2; i++) {
	129	// check index for stack
	130	if (index[i] < 0 \|\| index[i] > ntracks) {
	131	AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
	132	return kFALSE;
	133	}
	134	// get daughter and its PDG and charge
61f275d1	135	daughter = (AliMCParticle *)mc->GetTrack(index[i]);
f34f960b	136	pdg = TMath::Abs(daughter->Particle()->GetPdgCode());
	137	charge = (Short_t)(daughter->Particle()->GetPDG()->Charge() / 3);
	138	// check if it matches one definition
	139	if (fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge)) {
	140	fDaughter[0].SetGood();
	141	fDaughter[0].SetRefMC(daughter);
	142	fDaughter[0].SetLabel(index[i]);
	143	} else if (fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge)) {
	144	fDaughter[1].SetGood();
	145	fDaughter[1].SetRefMC(daughter);
	146	fDaughter[1].SetLabel(index[i]);
	147	}
	148	}
61f275d1	149
f34f960b	150	// return success if both daughters were assigned
	151	if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
	152	return kTRUE;
	153	} else {
	154	fDaughter[0].Reset();
	155	fDaughter[1].Reset();
	156	return kFALSE;
	157	}
	158	}
	159
	160	//__________________________________________________________________________________________________
	161	Bool_t AliRsnLoopEffPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
	162	{
b63357a0	163	//
	164	// Gets a particle in the MC event and try to assign it to the mother.
	165	// If it has two daughters, retrieve them and assign also them.
	166	// NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
	167	// for ESD and AOD, if available
	168	// ---
	169	// Implementation for AOD inputs
	170	//
	171
f34f960b	172	AliAODEvent *aod = rsnEvent->GetRefAOD();
61f275d1	173	TClonesArray listAOD = (TClonesArray )(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
61f275d1	174	AliAODMCParticle mother = (AliAODMCParticle )listAOD->At(ipart);
f34f960b	175	Int_t ntracks = listAOD->GetEntries();
61f275d1	176
f34f960b	177	// check PDG code and exit if it is wrong
f34f960b	178	if (TMath::Abs(mother->GetPdgCode()) != fDef->GetMotherPDG()) return kFALSE;
61f275d1	179
f34f960b	180	// check number of daughters and exit if it is not 2
f34f960b	181	if (mother->GetNDaughters() < 2) return kFALSE;
61f275d1	182
f34f960b	183	// check distance from primary vertex
	184	Double_t vprod[3] = {(Double_t)mother->Xv(), (Double_t)mother->Yv(), (Double_t)mother->Zv()};
	185	Double_t dv = DistanceFromPV(vprod[0], vprod[1], vprod[2]);
	186	if (dv > fMaxDistPV) {
	187	AliDebugClass(1, "Distant production vertex");
	188	return kFALSE;
	189	}
61f275d1	190
f34f960b	191	// get the daughters and check their PDG code and charge:
61f275d1	192	// if they match one of the pair daughter definitions,
f34f960b	193	// assign them as MC reference of the 'fDaughter' objects
	194	fDaughter[0].Reset();
	195	fDaughter[1].Reset();
	196	Int_t index[2] = {(Int_t)mother->GetDaughter(0), (Int_t)mother->GetDaughter(1)};
	197	Int_t i, pdg;
	198	Short_t charge;
	199	AliAODMCParticle *daughter = 0x0;
	200	for (i = 0; i < 2; i++) {
	201	// check index for stack
	202	if (index[i] < 0 \|\| index[i] > ntracks) {
	203	AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
	204	return kFALSE;
	205	}
	206	// get daughter and its PDG and charge
61f275d1	207	daughter = (AliAODMCParticle *)listAOD->At(index[i]);
f34f960b	208	pdg = TMath::Abs(daughter->GetPdgCode());
	209	charge = (Short_t)(daughter->Charge() / 3);
	210	// check if it matches one definition
	211	if (fDef->GetDef1().MatchesPDG(pdg) && fDef->GetDef1().MatchesChargeS(charge)) {
	212	fDaughter[0].SetGood();
	213	fDaughter[0].SetRefMC(daughter);
	214	fDaughter[0].SetLabel(index[i]);
	215	} else if (fDef->GetDef2().MatchesPDG(pdg) && fDef->GetDef2().MatchesChargeS(charge)) {
	216	fDaughter[1].SetGood();
	217	fDaughter[1].SetRefMC(daughter);
	218	fDaughter[1].SetLabel(index[i]);
	219	}
	220	}
61f275d1	221
f34f960b	222	// return success if both daughters were assigned
	223	if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
	224	return kTRUE;
	225	} else {
	226	fDaughter[0].Reset();
	227	fDaughter[1].Reset();
	228	return kFALSE;
	229	}
	230	}
	231
c865cb1d	232	//_____________________________________________________________________________
61f275d1	233	Int_t AliRsnLoopEffPair::DoLoop(AliRsnEvent rsn, AliRsnDaughterSelector , AliRsnEvent , AliRsnDaughterSelector )
c865cb1d	234	{
	235	//
	236	// Loop on event and fill containers
	237	//
	238
	239	// check event cuts
	240	if (!OkEvent(rsn)) return 0;
61f275d1	241
f34f960b	242	// retrieve output
61f275d1	243	fOutput = (AliRsnListOutput *)fOutputs[0];
61f275d1	244
c865cb1d	245	// check presence of MC reference
	246	if (!rsn->GetRefMC()) {
	247	AliError("Need a MC to compute efficiency");
	248	return 0;
	249	}
61f275d1	250
f34f960b	251	// check presence of event
	252	if (!rsn->GetRef()) {
	253	AliError("Need an event to compute efficiency");
	254	return 0;
	255	}
61f275d1	256
c865cb1d	257	// check event type:
	258	// must be ESD or AOD, and then use a bool to know in the rest
	259	if (!rsn->IsESD() && !rsn->IsAOD()) {
	260	AliError("Need to process ESD or AOD input");
	261	return 0;
	262	}
61f275d1	263
f34f960b	264	// retrieve the MC primary vertex position
	265	// and do some additional coherence checks
	266	Int_t i, npart = 0;
e187bd70	267	if (rsn->IsESD()) {
f34f960b	268	TArrayF mcVertex(3);
	269	AliGenEventHeader *genEH = rsn->GetRefMCESD()->GenEventHeader();
	270	genEH->PrimaryVertex(mcVertex);
	271	for (i = 0; i < 3; i++) fVertex[i] = (Double_t)mcVertex[i];
	272	npart = rsn->GetRefMCESD()->GetNumberOfTracks();
c865cb1d	273	} else {
f34f960b	274	for (i = 0; i < 3; i++) fVertex[i] = 0.0;
f34f960b	275	AliAODEvent *aod = rsn->GetRefMCAOD();
61f275d1	276	TClonesArray listAOD = (TClonesArray )(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
f34f960b	277	if (listAOD) npart = listAOD->GetEntries();
61f275d1	278	AliAODMCHeader mcH = static_cast<AliAODMCHeader >(aod->FindListObject(AliAODMCHeader::StdBranchName()));
f34f960b	279	if (mcH) mcH->GetVertex(fVertex);
c865cb1d	280	}
61f275d1	281
c865cb1d	282	// check number of particles
	283	if (!npart) {
	284	AliInfo("Empty event");
	285	return 0;
	286	}
61f275d1	287
c865cb1d	288	// utility variables
f34f960b	289	Int_t ipart, istep, count = 0, nsteps = fSteps.GetEntries();
	290	Int_t ntracks = rsn->GetAbsoluteSum();
	291	AliRsnDaughter check;
61f275d1	292
c865cb1d	293	// loop over particles
c865cb1d	294	for (ipart = 0; ipart < npart; ipart++) {
f34f960b	295	// check i-th particle
	296	if (!AssignMotherAndDaughters(rsn, ipart)) continue;
	297	// if assignment was successful, for first step, use MC info
	298	fDaughter[0].SetRef(fDaughter[0].GetRefMC());
	299	fDaughter[1].SetRef(fDaughter[1].GetRefMC());
	300	fMother.SetRefEvent(rsn);
	301	fMother.ComputeSum(fDef->GetDef1().GetMass(), fDef->GetDef2().GetMass(), fDef->GetMotherMass());
	302	fOutput->Fill(&fMother, 0);
c865cb1d	303	count++;
f34f960b	304	// for each further step, try to find two tracks which pass the related cuts
f34f960b	305	for (istep = 0; istep < nsteps; istep++) {
61f275d1	306	AliRsnCutManager cuts = (AliRsnCutManager )fSteps[istep];
f34f960b	307	fDaughter[0].SetBad();
	308	fDaughter[1].SetBad();
	309	for (i = 0; i < ntracks; i++) {
	310	rsn->SetDaughter(check, i);
	311	if (!cuts->PassCommonDaughterCuts(&check)) continue;
	312	if (TMath::Abs(check.GetLabel()) == fDaughter[0].GetLabel()) {
	313	if (!cuts->PassDaughter1Cuts(&check)) continue;
	314	fDaughter[0].SetRef(check.GetRef());
	315	fDaughter[0].SetGood();
	316	} else if (TMath::Abs(check.GetLabel()) == fDaughter[1].GetLabel()) {
	317	if (!cuts->PassDaughter2Cuts(&check)) continue;
	318	fDaughter[1].SetRef(check.GetRef());
	319	fDaughter[1].SetGood();
	320	}
	321	if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
	322	fMother.ComputeSum(fDef->GetDef1().GetMass(), fDef->GetDef2().GetMass(), fDef->GetMotherMass());
	323	if (cuts->PassMotherCuts(&fMother)) {
	324	fOutput->Fill(&fMother, 1 + istep);
	325	break;
	326	}
	327	}
c865cb1d	328	}
	329	}
	330	}
61f275d1	331
c865cb1d	332	return count;
c865cb1d	333	}