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