[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnLoopPair.cxx

//
// *** Class AliRsnLoopPair ***
//
// "Core" method for defining the work on a pari of particles.
// For one analysis, one must setup one of this for each pair he wants to analyze,
// adding to it all analysis which he desires to do.
// Here he defines the cuts, and the particle types and charges, and can add
// functions which do different operations on the same pair, and some binning
// with respect to some kinematic variables (eta, momentum)
//
// authors: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
//          M. Vala (email: martin.vala@cern.ch)
//

#include <Riostream.h>
#include <TList.h>
#include <TEntryList.h>

#include "AliLog.h"

#include "AliRsnEvent.h"
#include "AliRsnPairDef.h"
#include "AliRsnMother.h"
#include "AliRsnCutSet.h"
#include "AliRsnDaughterSelector.h"

#include "AliRsnLoopPair.h"

ClassImp(AliRsnLoopPair)

//_____________________________________________________________________________
AliRsnLoopPair::AliRsnLoopPair(const char *name, AliRsnPairDef *def, Bool_t isMixed) :
   AliRsnLoop(name, isMixed),
   fTrueMC(kFALSE),
   fOnlyTrue(kFALSE),
   fCheckDecay(kFALSE),
   fRangeY(1E20),
   fPairDef(def),
   fPairCuts(0x0),
   fMother()
{
//
// Default constructor
//

   fListID[0] = -1;
   fListID[1] = -1;
}

//_____________________________________________________________________________
AliRsnLoopPair::AliRsnLoopPair(const AliRsnLoopPair& copy) :
   AliRsnLoop(copy),
   fTrueMC(copy.fTrueMC),
   fOnlyTrue(copy.fOnlyTrue),
   fCheckDecay(copy.fCheckDecay),
   fRangeY(copy.fRangeY),
   fPairDef(copy.fPairDef),
   fPairCuts(copy.fPairCuts),
   fMother(copy.fMother)
{
//
// Copy constructor
//

   fListID[0] = copy.fListID[0];
   fListID[1] = copy.fListID[1];
}

//_____________________________________________________________________________
AliRsnLoopPair& AliRsnLoopPair::operator=(const AliRsnLoopPair& copy)
{
//
// Assignment operator
//

   AliRsnLoop::operator=(copy);
   
   fTrueMC = copy.fTrueMC;
   fOnlyTrue = copy.fOnlyTrue;
   fCheckDecay = copy.fCheckDecay;
   fRangeY = copy.fRangeY;
   fPairDef = copy.fPairDef;
   fPairCuts = copy.fPairCuts;
   fMother = copy.fMother;
   fListID[0] = copy.fListID[0];
   fListID[1] = copy.fListID[1];

   return (*this);
}

//_____________________________________________________________________________
AliRsnLoopPair::~AliRsnLoopPair()
{
//
// Destructor
//
}

//_____________________________________________________________________________
void AliRsnLoopPair::Print(Option_t* /*option*/) const
{
//
// Prints info about pair
//

   AliRsnLoop::Print();
}

//_____________________________________________________________________________
Bool_t AliRsnLoopPair::Init(const char *prefix, TList *list)
{
//
// Initialization function.
// Loops on all functions and eventual the ntuple, to initialize output objects.
//

   // assign data members relations
   fMother.SetDaughter(0, &fDaughter[0]);
   fMother.SetDaughter(1, &fDaughter[1]);
   AliInfo(Form("[%s] Initialization", GetName()));
   
   TString name(prefix);
   name += '_';
   name += GetName();
   if (IsMixed()) name.Prepend("mix_");
   
   return AliRsnLoop::Init(name.Data(), list);
}

//_____________________________________________________________________________
Int_t AliRsnLoopPair::DoLoop
(AliRsnEvent *evMain, AliRsnDaughterSelector *selMain, AliRsnEvent *evMix, AliRsnDaughterSelector *selMix)
{
//
// Loop function.
// Computes what is needed from passed events.
// Returns the number of pairs successfully processed.
//

   if (fIsMixed) {
      AliDebugClass(3, Form("[%s]: event-mixing loop", GetName()));
      if (!evMix || !selMix) {
         AliError(Form("[%s] NULL mixed event when mixing is required: cannot process", GetName()));
         return 0;
      }
   } else {
      AliDebugClass(3, Form("[%s]: single-event loop", GetName()));
      evMix = evMain;
      selMix = selMain;
   }
   fMother.SetRefEvent(evMain);
   
   // check cuts
   if (!OkEvent(evMain)) {
      AliDebugClass(3, Form("[%s]: main event not accepted", GetName()));
      return 0;
   }
   if (!OkEvent(evMix)) {
      AliDebugClass(3, Form("[%s]: mixed event not accepted", GetName()));
      return 0;
   }
   
   // if it is required to loop over True MC, do this here and skip the rest of the method
   if (fTrueMC) return LoopTrueMC(evMain);
   
   Int_t i0, i1, start, npairs = 0;
   
   TEntryList *list0 = selMain->GetSelected(fListID[0], fPairDef->GetDef1().GetChargeC());
   TEntryList *list1 = selMix ->GetSelected(fListID[1], fPairDef->GetDef2().GetChargeC());
   if (!list0 || !list1) {
      AliError("Can't process NULL lists");
      return 0;
   }
   AliDebugClass(3, Form("[%s]: list counts: %d, %d", GetName(), (Int_t)list0->GetN(), (Int_t)list1->GetN()));
   if (!list0->GetN() || !list1->GetN()) {
      AliDebugClass(3, Form("[%s]: at least one list is empty", GetName()));
      return 0;
   }
   
   TObjArrayIter next(&fOutputs);
   AliRsnListOutput *out = 0x0;
   
   for (i0 = 0; i0 < list0->GetN(); i0++) {
      evMain->SetDaughter(fDaughter[0], (Int_t)list0->GetEntry(i0));
      fDaughter[0].FillP(fPairDef->GetDef1().GetMass());
      start = 0;
      if (!fIsMixed && list0 == list1) start = i0 + 1;
      for (i1 = start; i1 < list1->GetN(); i1++) {
         AliDebugClass(4, Form("Checking entries pair: %d (%d) with %d (%d)", (Int_t)i0, (Int_t)list0->GetEntry(i0), (Int_t)i1, (Int_t)list1->GetEntry(i1)));
         evMix->SetDaughter(fDaughter[1], (Int_t)list1->GetEntry(i1));
         fDaughter[1].FillP(fPairDef->GetDef2().GetMass());
         fMother.Sum(0) = fDaughter[0].Prec() + fDaughter[1].Prec();
         fMother.Sum(1) = fDaughter[0].Psim() + fDaughter[1].Psim();
         fMother.Ref(0).SetXYZM(fMother.Sum(0).X(), fMother.Sum(0).Y(), fMother.Sum(0).Z(), fPairDef->GetMotherMass());
         fMother.Ref(1).SetXYZM(fMother.Sum(1).X(), fMother.Sum(1).Y(), fMother.Sum(1).Z(), fPairDef->GetMotherMass());
         // check rapidity range
         if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
            AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
            continue;
         }
         // check mother
         if (fOnlyTrue) {
            if (!IsTrueMother()) {
               AliDebugClass(2, Form("[%s]: candidate mother is not true", GetName()));
               continue;
            }
         }
         // check cuts
         if (fPairCuts) {
            if (!fPairCuts->IsSelected(&fMother)) {
               AliDebugClass(2, Form("[%s]: candidate mother didn't pass the cuts", GetName()));
               continue;
            }
         }
         // fill outputs
         next.Reset();
         while ( (out = (AliRsnListOutput*)next()) ) {
            if (out->Fill(&fMother)) npairs++;
            else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
         }
      }
   }

   return npairs;
}

//_____________________________________________________________________________
Bool_t AliRsnLoopPair::IsTrueMother()
{
//
// Checks to see if the mother comes from a true resonance.
// It is triggered by the 'SetOnlyTrue()' function
//
      
   // check #1:
   // daughters have same mother with the right PDG code
   Int_t commonPDG = fMother.CommonMother();
   if (commonPDG != fPairDef->GetMotherPDG()) return kFALSE;
   AliDebugClass(1, "Found a true mother");
   
   // check #2:
   // checks if daughter have the right particle type
   // (activated by fCheckDecay)
   if (fCheckDecay) {
      AliRsnDaughterDef &def1 = fPairDef->GetDef1();
      AliRsnDaughterDef &def2 = fPairDef->GetDef2();
      if (!def1.MatchesPID(&fDaughter[0])) return kFALSE;
      if (!def2.MatchesPID(&fDaughter[1])) return kFALSE;
   }
   AliDebugClass(1, "Decay products match");
   
   return kTRUE;
}

//__________________________________________________________________________________________________
Bool_t AliRsnLoopPair::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 AliRsnLoopPair::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()) != fPairDef->GetMotherPDG()) return kFALSE;
   
   // check number of daughters and exit if it is not 2
   if (motherP->GetNDaughters() < 2) return kFALSE;
   //if (!stack->IsPhysicalPrimary(ipart)) 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->GetDaughter(0), motherP->GetDaughter(1)};
   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 (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
         fDaughter[0].SetGood();
         fDaughter[0].SetRefMC(daughter);
         fDaughter[0].SetLabel(index[i]);
      } else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->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 AliRsnLoopPair::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);
   
   // check PDG code and exit if it is wrong
   if (TMath::Abs(mother->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
   
   // check number of daughters and exit if it is not 2
   if (mother->GetNDaughters() < 2) return kFALSE;
   if (!mother->IsPrimary()) 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 ntracks = listAOD->GetEntriesFast();
   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 (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
         fDaughter[0].SetGood();
         fDaughter[0].SetRefMC(daughter);
         fDaughter[0].SetLabel(index[i]);
      } else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->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 AliRsnLoopPair::LoopTrueMC(AliRsnEvent *rsn)
{
//
// Loop on event and fill containers
//

   // check presence of MC reference
   if (!rsn->GetRefMC()) {
      AliError("Need a MC 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
   //Double_t fVertex[3] = {0.0, 0.0, 0.0};
   Int_t 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, count = 0;
   AliRsnDaughter check;
   
   TObjArrayIter next(&fOutputs);
   AliRsnListOutput *out = 0x0;
   
   // 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(fPairDef->GetDef1().GetMass(), fPairDef->GetDef2().GetMass(), fPairDef->GetMotherMass());
      // check rapidity range
      if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
         AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
         continue;
      }
      // fill outputs
      next.Reset();
      while ( (out = (AliRsnListOutput*)next()) ) {
         if (out->Fill(&fMother)) count++;
         else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
      }
   }
   
   return count;
}
Commit	Line	Data
c865cb1d	1	//
	2	// * Class AliRsnLoopPair *
	3	//
	4	// "Core" method for defining the work on a pari of particles.
	5	// For one analysis, one must setup one of this for each pair he wants to analyze,
	6	// adding to it all analysis which he desires to do.
	7	// Here he defines the cuts, and the particle types and charges, and can add
	8	// functions which do different operations on the same pair, and some binning
	9	// with respect to some kinematic variables (eta, momentum)
	10	//
	11	// authors: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
	12	// M. Vala (email: martin.vala@cern.ch)
	13	//
	14
f34f960b	15	#include <Riostream.h>
c865cb1d	16	#include <TList.h>
	17	#include <TEntryList.h>
	18
	19	#include "AliLog.h"
	20
b63357a0	21	#include "AliRsnEvent.h"
b63357a0	22	#include "AliRsnPairDef.h"
c865cb1d	23	#include "AliRsnMother.h"
	24	#include "AliRsnCutSet.h"
	25	#include "AliRsnDaughterSelector.h"
	26
	27	#include "AliRsnLoopPair.h"
	28
	29	ClassImp(AliRsnLoopPair)
	30
	31	//_____________________________________________________________________________
	32	AliRsnLoopPair::AliRsnLoopPair(const char name, AliRsnPairDef def, Bool_t isMixed) :
	33	AliRsnLoop(name, isMixed),
b63357a0	34	fTrueMC(kFALSE),
c865cb1d	35	fOnlyTrue(kFALSE),
c865cb1d	36	fCheckDecay(kFALSE),
b63357a0	37	fRangeY(1E20),
c865cb1d	38	fPairDef(def),
	39	fPairCuts(0x0),
	40	fMother()
	41	{
	42	//
	43	// Default constructor
	44	//
	45
	46	fListID[0] = -1;
	47	fListID[1] = -1;
	48	}
	49
	50	//_____________________________________________________________________________
	51	AliRsnLoopPair::AliRsnLoopPair(const AliRsnLoopPair& copy) :
	52	AliRsnLoop(copy),
b63357a0	53	fTrueMC(copy.fTrueMC),
c865cb1d	54	fOnlyTrue(copy.fOnlyTrue),
c865cb1d	55	fCheckDecay(copy.fCheckDecay),
b63357a0	56	fRangeY(copy.fRangeY),
c865cb1d	57	fPairDef(copy.fPairDef),
	58	fPairCuts(copy.fPairCuts),
	59	fMother(copy.fMother)
	60	{
	61	//
	62	// Copy constructor
	63	//
	64
	65	fListID[0] = copy.fListID[0];
	66	fListID[1] = copy.fListID[1];
	67	}
	68
	69	//_____________________________________________________________________________
	70	AliRsnLoopPair& AliRsnLoopPair::operator=(const AliRsnLoopPair& copy)
	71	{
b63357a0	72	//
	73	// Assignment operator
	74	//
	75
c865cb1d	76	AliRsnLoop::operator=(copy);
c865cb1d	77
b63357a0	78	fTrueMC = copy.fTrueMC;
c865cb1d	79	fOnlyTrue = copy.fOnlyTrue;
c865cb1d	80	fCheckDecay = copy.fCheckDecay;
b63357a0	81	fRangeY = copy.fRangeY;
c865cb1d	82	fPairDef = copy.fPairDef;
	83	fPairCuts = copy.fPairCuts;
	84	fMother = copy.fMother;
	85	fListID[0] = copy.fListID[0];
	86	fListID[1] = copy.fListID[1];
	87
	88	return (*this);
	89	}
	90
	91	//_____________________________________________________________________________
	92	AliRsnLoopPair::~AliRsnLoopPair()
	93	{
	94	//
	95	// Destructor
	96	//
	97	}
	98
	99	//_____________________________________________________________________________
	100	void AliRsnLoopPair::Print(Option_t* /option/) const
	101	{
	102	//
	103	// Prints info about pair
	104	//
	105
	106	AliRsnLoop::Print();
	107	}
	108
	109	//_____________________________________________________________________________
	110	Bool_t AliRsnLoopPair::Init(const char prefix, TList list)
	111	{
	112	//
	113	// Initialization function.
	114	// Loops on all functions and eventual the ntuple, to initialize output objects.
	115	//
	116
	117	// assign data members relations
	118	fMother.SetDaughter(0, &fDaughter[0]);
	119	fMother.SetDaughter(1, &fDaughter[1]);
	120	AliInfo(Form("[%s] Initialization", GetName()));
	121
	122	TString name(prefix);
	123	name += '_';
	124	name += GetName();
	125	if (IsMixed()) name.Prepend("mix_");
	126
	127	return AliRsnLoop::Init(name.Data(), list);
	128	}
	129
	130	//_____________________________________________________________________________
	131	Int_t AliRsnLoopPair::DoLoop
	132	(AliRsnEvent evMain, AliRsnDaughterSelector selMain, AliRsnEvent evMix, AliRsnDaughterSelector selMix)
	133	{
	134	//
	135	// Loop function.
	136	// Computes what is needed from passed events.
	137	// Returns the number of pairs successfully processed.
	138	//
	139
	140	if (fIsMixed) {
b63357a0	141	AliDebugClass(3, Form("[%s]: event-mixing loop", GetName()));
c865cb1d	142	if (!evMix \|\| !selMix) {
	143	AliError(Form("[%s] NULL mixed event when mixing is required: cannot process", GetName()));
	144	return 0;
	145	}
	146	} else {
b63357a0	147	AliDebugClass(3, Form("[%s]: single-event loop", GetName()));
c865cb1d	148	evMix = evMain;
	149	selMix = selMain;
	150	}
	151	fMother.SetRefEvent(evMain);
	152
	153	// check cuts
	154	if (!OkEvent(evMain)) {
b63357a0	155	AliDebugClass(3, Form("[%s]: main event not accepted", GetName()));
c865cb1d	156	return 0;
	157	}
	158	if (!OkEvent(evMix)) {
b63357a0	159	AliDebugClass(3, Form("[%s]: mixed event not accepted", GetName()));
c865cb1d	160	return 0;
	161	}
	162
b63357a0	163	// if it is required to loop over True MC, do this here and skip the rest of the method
	164	if (fTrueMC) return LoopTrueMC(evMain);
	165
c865cb1d	166	Int_t i0, i1, start, npairs = 0;
	167
	168	TEntryList *list0 = selMain->GetSelected(fListID[0], fPairDef->GetDef1().GetChargeC());
	169	TEntryList *list1 = selMix ->GetSelected(fListID[1], fPairDef->GetDef2().GetChargeC());
	170	if (!list0 \|\| !list1) {
	171	AliError("Can't process NULL lists");
	172	return 0;
	173	}
b63357a0	174	AliDebugClass(3, Form("[%s]: list counts: %d, %d", GetName(), (Int_t)list0->GetN(), (Int_t)list1->GetN()));
c865cb1d	175	if (!list0->GetN() \|\| !list1->GetN()) {
b63357a0	176	AliDebugClass(3, Form("[%s]: at least one list is empty", GetName()));
c865cb1d	177	return 0;
	178	}
	179
	180	TObjArrayIter next(&fOutputs);
	181	AliRsnListOutput *out = 0x0;
	182
	183	for (i0 = 0; i0 < list0->GetN(); i0++) {
f34f960b	184	evMain->SetDaughter(fDaughter[0], (Int_t)list0->GetEntry(i0));
f34f960b	185	fDaughter[0].FillP(fPairDef->GetDef1().GetMass());
c865cb1d	186	start = 0;
c865cb1d	187	if (!fIsMixed && list0 == list1) start = i0 + 1;
c865cb1d	188	for (i1 = start; i1 < list1->GetN(); i1++) {
f34f960b	189	AliDebugClass(4, Form("Checking entries pair: %d (%d) with %d (%d)", (Int_t)i0, (Int_t)list0->GetEntry(i0), (Int_t)i1, (Int_t)list1->GetEntry(i1)));
	190	evMix->SetDaughter(fDaughter[1], (Int_t)list1->GetEntry(i1));
	191	fDaughter[1].FillP(fPairDef->GetDef2().GetMass());
	192	fMother.Sum(0) = fDaughter[0].Prec() + fDaughter[1].Prec();
	193	fMother.Sum(1) = fDaughter[0].Psim() + fDaughter[1].Psim();
	194	fMother.Ref(0).SetXYZM(fMother.Sum(0).X(), fMother.Sum(0).Y(), fMother.Sum(0).Z(), fPairDef->GetMotherMass());
	195	fMother.Ref(1).SetXYZM(fMother.Sum(1).X(), fMother.Sum(1).Y(), fMother.Sum(1).Z(), fPairDef->GetMotherMass());
b63357a0	196	// check rapidity range
	197	if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
	198	AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
	199	continue;
c865cb1d	200	}
c865cb1d	201	// check mother
f34f960b	202	if (fOnlyTrue) {
	203	if (!IsTrueMother()) {
	204	AliDebugClass(2, Form("[%s]: candidate mother is not true", GetName()));
	205	continue;
	206	}
c865cb1d	207	}
b63357a0	208	// check cuts
	209	if (fPairCuts) {
	210	if (!fPairCuts->IsSelected(&fMother)) {
	211	AliDebugClass(2, Form("[%s]: candidate mother didn't pass the cuts", GetName()));
	212	continue;
	213	}
	214	}
c865cb1d	215	// fill outputs
	216	next.Reset();
	217	while ( (out = (AliRsnListOutput*)next()) ) {
	218	if (out->Fill(&fMother)) npairs++;
	219	else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
	220	}
	221	}
	222	}
f34f960b	223
c865cb1d	224	return npairs;
	225	}
	226
	227	//_____________________________________________________________________________
f34f960b	228	Bool_t AliRsnLoopPair::IsTrueMother()
c865cb1d	229	{
c865cb1d	230	//
f34f960b	231	// Checks to see if the mother comes from a true resonance.
f34f960b	232	// It is triggered by the 'SetOnlyTrue()' function
c865cb1d	233	//
f34f960b	234
	235	// check #1:
	236	// daughters have same mother with the right PDG code
	237	Int_t commonPDG = fMother.CommonMother();
	238	if (commonPDG != fPairDef->GetMotherPDG()) return kFALSE;
b63357a0	239	AliDebugClass(1, "Found a true mother");
c865cb1d	240
f34f960b	241	// check #2:
	242	// checks if daughter have the right particle type
	243	// (activated by fCheckDecay)
	244	if (fCheckDecay) {
	245	AliRsnDaughterDef &def1 = fPairDef->GetDef1();
	246	AliRsnDaughterDef &def2 = fPairDef->GetDef2();
	247	if (!def1.MatchesPID(&fDaughter[0])) return kFALSE;
	248	if (!def2.MatchesPID(&fDaughter[1])) return kFALSE;
c865cb1d	249	}
b63357a0	250	AliDebugClass(1, "Decay products match");
c865cb1d	251
f34f960b	252	return kTRUE;
c865cb1d	253	}
b63357a0	254
	255	//__________________________________________________________________________________________________
	256	Bool_t AliRsnLoopPair::AssignMotherAndDaughters(AliRsnEvent *rsnEvent, Int_t ipart)
	257	{
	258	//
	259	// Calls the appropriate assignment method
	260	//
	261
	262	// setup pointers
	263	fMother.SetDaughter(0, &fDaughter[0]);
	264	fMother.SetDaughter(1, &fDaughter[1]);
	265	fMother.SetRefEvent(rsnEvent);
	266	fDaughter[0].SetOwnerEvent(rsnEvent);
	267	fDaughter[1].SetOwnerEvent(rsnEvent);
	268
	269	if (rsnEvent->IsESD())
	270	return AssignMotherAndDaughtersESD(rsnEvent, ipart);
	271	else if (rsnEvent->IsAOD())
	272	return AssignMotherAndDaughtersAOD(rsnEvent, ipart);
	273	else {
	274	AliError("Unrecognized input event");
	275	return kFALSE;
	276	}
	277	}
	278
	279	//__________________________________________________________________________________________________
	280	Bool_t AliRsnLoopPair::AssignMotherAndDaughtersESD(AliRsnEvent *rsnEvent, Int_t ipart)
	281	{
	282	//
	283	// Gets a particle in the MC event and try to assign it to the mother.
	284	// If it has two daughters, retrieve them and assign also them.
	285	// NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
	286	// for ESD and AOD, if available
	287	// ---
	288	// Implementation for ESD inputs
	289	//
	290
	291	AliMCEvent *mc = rsnEvent->GetRefMCESD();
	292	AliStack *stack = mc->Stack();
	293	AliMCParticle mother = (AliMCParticle)mc->GetTrack(ipart);
	294	TParticle *motherP = mother->Particle();
	295	Int_t ntracks = stack->GetNtrack();
	296
	297	// check PDG code and exit if it is wrong
	298	if (TMath::Abs(motherP->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
	299
	300	// check number of daughters and exit if it is not 2
	301	if (motherP->GetNDaughters() < 2) return kFALSE;
	302	//if (!stack->IsPhysicalPrimary(ipart)) return kFALSE;
	303
	304	/*
	305	// check distance from primary vertex
	306	TLorentzVector vprod;
	307	motherP->ProductionVertex(vprod);
	308	if (!CheckDistanceFromPV(vprod.X(), vprod.Y(), vprod.Z())) {
	309	AliDebugClass(1, "Distant production vertex");
	310	return kFALSE;
	311	}
	312	*/
	313
	314	// get the daughters and check their PDG code and charge:
	315	// if they match one of the pair daughter definitions,
	316	// assign them as MC reference of the 'fDaughter' objects
	317	fDaughter[0].Reset();
318	fDaughter[1].Reset();
319	Int_t index[2] = {motherP->GetDaughter(0), motherP->GetDaughter(1)};
320	Int_t i, pdg;
321	Short_t charge;
322	AliMCParticle *daughter = 0x0;
323	for (i = 0; i < 2; i++) {
324	// check index for stack
325	if (index[i] < 0 \|\| index[i] > ntracks) {
326	AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
327	return kFALSE;
328	}
329	// get daughter and its PDG and charge
330	daughter = (AliMCParticle*)mc->GetTrack(index[i]);
331	pdg = TMath::Abs(daughter->Particle()->GetPdgCode());
332	charge = (Short_t)(daughter->Particle()->GetPDG()->Charge() / 3);
333	// check if it matches one definition
334	if (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
335	fDaughter[0].SetGood();
336	fDaughter[0].SetRefMC(daughter);
337	fDaughter[0].SetLabel(index[i]);
338	} else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->GetDef2().MatchesChargeS(charge)) {
339	fDaughter[1].SetGood();
340	fDaughter[1].SetRefMC(daughter);
341	fDaughter[1].SetLabel(index[i]);
342	}
343	}
344
345	// return success if both daughters were assigned
346	if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
347	return kTRUE;
348	} else {
349	fDaughter[0].Reset();
350	fDaughter[1].Reset();
351	return kFALSE;
352	}
353	}
354
355	//__________________________________________________________________________________________________
356	Bool_t AliRsnLoopPair::AssignMotherAndDaughtersAOD(AliRsnEvent *rsnEvent, Int_t ipart)
357	{
358	//
359	// Gets a particle in the MC event and try to assign it to the mother.
360	// If it has two daughters, retrieve them and assign also them.
361	// NOTE: assignment is done only for MC, since reconstructed match is assigned in the same way
362	// for ESD and AOD, if available
363	// ---
364	// Implementation for AOD inputs
365	//
366
367	AliAODEvent *aod = rsnEvent->GetRefAOD();
368	TClonesArray listAOD = (TClonesArray)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
369	AliAODMCParticle mother = (AliAODMCParticle)listAOD->At(ipart);
370
371	// check PDG code and exit if it is wrong
372	if (TMath::Abs(mother->GetPdgCode()) != fPairDef->GetMotherPDG()) return kFALSE;
373
374	// check number of daughters and exit if it is not 2
375	if (mother->GetNDaughters() < 2) return kFALSE;
376	if (!mother->IsPrimary()) return kFALSE;
377
378	/*
379	// check distance from primary vertex
380	Double_t vprod[3] = {(Double_t)mother->Xv(), (Double_t)mother->Yv(), (Double_t)mother->Zv()};
381	Double_t dv = DistanceFromPV(vprod[0], vprod[1], vprod[2]);
382	if (dv > fMaxDistPV) {
383	AliDebugClass(1, "Distant production vertex");
384	return kFALSE;
385	}
386	*/
387
388	// get the daughters and check their PDG code and charge:
389	// if they match one of the pair daughter definitions,
390	// assign them as MC reference of the 'fDaughter' objects
391	fDaughter[0].Reset();
392	fDaughter[1].Reset();
393	Int_t index[2] = {(Int_t)mother->GetDaughter(0), (Int_t)mother->GetDaughter(1)};
394	Int_t ntracks = listAOD->GetEntriesFast();
395	Int_t i, pdg;
396	Short_t charge;
397	AliAODMCParticle *daughter = 0x0;
398	for (i = 0; i < 2; i++) {
399	// check index for stack
400	if (index[i] < 0 \|\| index[i] > ntracks) {
401	AliError(Form("Index %d overflow: value = %d, max = %d", i, index[i], ntracks));
402	return kFALSE;
403	}
404	// get daughter and its PDG and charge
405	daughter = (AliAODMCParticle*)listAOD->At(index[i]);
406	pdg = TMath::Abs(daughter->GetPdgCode());
407	charge = (Short_t)(daughter->Charge() / 3);
408	// check if it matches one definition
409	if (fPairDef->GetDef1().MatchesPDG(pdg) && fPairDef->GetDef1().MatchesChargeS(charge)) {
410	fDaughter[0].SetGood();
411	fDaughter[0].SetRefMC(daughter);
412	fDaughter[0].SetLabel(index[i]);
413	} else if (fPairDef->GetDef2().MatchesPDG(pdg) && fPairDef->GetDef2().MatchesChargeS(charge)) {
414	fDaughter[1].SetGood();
415	fDaughter[1].SetRefMC(daughter);
416	fDaughter[1].SetLabel(index[i]);
417	}
418	}
419
420	// return success if both daughters were assigned
421	if (fDaughter[0].IsOK() && fDaughter[1].IsOK()) {
422	return kTRUE;
423	} else {
424	fDaughter[0].Reset();
425	fDaughter[1].Reset();
426	return kFALSE;
427	}
428	}
429
430	//_____________________________________________________________________________
431	Int_t AliRsnLoopPair::LoopTrueMC(AliRsnEvent *rsn)
432	{
433	//
434	// Loop on event and fill containers
435	//
436
437	// check presence of MC reference
438	if (!rsn->GetRefMC()) {
439	AliError("Need a MC to compute efficiency");
440	return 0;
441	}
442
443	// check event type:
444	// must be ESD or AOD, and then use a bool to know in the rest
445	if (!rsn->IsESD() && !rsn->IsAOD()) {
446	AliError("Need to process ESD or AOD input");
447	return 0;
448	}
449
450	// retrieve the MC primary vertex position
451	// and do some additional coherence checks
452	//Double_t fVertex[3] = {0.0, 0.0, 0.0};
453	Int_t npart = 0;
454	if (rsn->IsESD()) {
455	//TArrayF mcVertex(3);
456	//AliGenEventHeader *genEH = rsn->GetRefMCESD()->GenEventHeader();
457	//genEH->PrimaryVertex(mcVertex);
458	//for (i = 0; i < 3; i++) fVertex[i] = (Double_t)mcVertex[i];
459	npart = rsn->GetRefMCESD()->GetNumberOfTracks();
460	} else {
461	//for (i = 0; i < 3; i++) fVertex[i] = 0.0;
462	AliAODEvent *aod = rsn->GetRefMCAOD();
463	TClonesArray listAOD = (TClonesArray)(aod->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
464	if (listAOD) npart = listAOD->GetEntries();
465	//AliAODMCHeader mcH = static_cast<AliAODMCHeader>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
466	//if (mcH) mcH->GetVertex(fVertex);
467	}
468
469	// check number of particles
470	if (!npart) {
471	AliInfo("Empty event");
472	return 0;
473	}
474
475	// utility variables
476	Int_t ipart, count = 0;
477	AliRsnDaughter check;
478
479	TObjArrayIter next(&fOutputs);
480	AliRsnListOutput *out = 0x0;
481
482	// loop over particles
483	for (ipart = 0; ipart < npart; ipart++) {
484	// check i-th particle
485	if (!AssignMotherAndDaughters(rsn, ipart)) continue;
486	// if assignment was successful, for first step, use MC info
487	fDaughter[0].SetRef(fDaughter[0].GetRefMC());
488	fDaughter[1].SetRef(fDaughter[1].GetRefMC());
489	fMother.SetRefEvent(rsn);
490	fMother.ComputeSum(fPairDef->GetDef1().GetMass(), fPairDef->GetDef2().GetMass(), fPairDef->GetMotherMass());
491	// check rapidity range
492	if (TMath::Abs(fMother.Rapidity(0)) > fRangeY) {
493	AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName()));
494	continue;
495	}
496	// fill outputs
497	next.Reset();
498	while ( (out = (AliRsnListOutput*)next()) ) {
499	if (out->Fill(&fMother)) count++;
500	else AliDebugClass(3, Form("[%s]: failed computation", GetName()));
501	}
502	}
503
504	return count;
505	}