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

//
// Class AliRsnAnalysisTaskEffPair
//
// Inherits from basic AliRsnAnalysisTaskEff for efficiency,
// and computed efficiencies for pairs
//
// author: Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
//

#include "AliStack.h"
#include "AliMCEvent.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"

#include "AliRsnPairDef.h"
#include "AliRsnCutManager.h"
#include "AliRsnAnalysisTaskEffPair.h"

ClassImp(AliRsnAnalysisTaskEffPair)

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

//_____________________________________________________________________________
AliRsnAnalysisTaskEffPair::AliRsnAnalysisTaskEffPair(const AliRsnAnalysisTaskEffPair& copy) :
   AliRsnAnalysisTaskEff(copy),
   fMother()
{
//
// Copy constrtuctor.
//
}

//_____________________________________________________________________________
AliRsnAnalysisTaskEffPair& AliRsnAnalysisTaskEffPair::operator=(const AliRsnAnalysisTaskEffPair& copy)
{
//
// Assignment operator.
// Owned data members are meaningless for this operator.
//
   
   AliRsnAnalysisTaskEff::operator=(copy);
   return (*this);
}

//_____________________________________________________________________________
Bool_t AliRsnAnalysisTaskEffPair::RsnEventProcess()
{
//
// Preprocessing to the event
// Calls the same function in base class plus somethin specific
//

   if (!AliRsnAnalysisTaskEff::RsnEventProcess()) return kFALSE;
   
   // assign this event to the mother reference
   fMother.SetRefEvent(&fRsnEvent[0]);
   fDaughter[0].SetOwnerEvent(&fRsnEvent[0]);
   fDaughter[1].SetOwnerEvent(&fRsnEvent[0]);
   
   return kTRUE;
}

//_____________________________________________________________________________
Int_t AliRsnAnalysisTaskEffPair::NGoodSteps()
{
//
// Checks how many 'reconstruction' steps are passed by current daughter
//

   Int_t istep, count = 0;
   Int_t nSteps = fStepsRec.GetEntries();
   
   for (istep = 0; istep < nSteps; istep++) {
      AliRsnCutManager *cutMgr = (AliRsnCutManager*)fStepsRec[istep];
      
      if (!cutMgr->PassCommonDaughterCuts(&fDaughter[0])) break;
      if (!cutMgr->PassCommonDaughterCuts(&fDaughter[1])) break;
      if (!cutMgr->PassDaughter1Cuts(&fDaughter[0])) break;
      if (!cutMgr->PassDaughter2Cuts(&fDaughter[1])) break;
      if (!cutMgr->PassMotherCuts(&fMother)) break;
      
      count++;
   }
   
   return count;
}

//_____________________________________________________________________________
void AliRsnAnalysisTaskEffPair::ProcessEventESD()
{
//
// Process current event with the definitions of the specified step in MC list
// and store results in the container slot defined in second argument.
// It is associated with the AliCFContainer with the name of the pair.
//

   AliESDEvent   *esd   = fRsnEvent[0].GetRefESD();
   AliMCEvent    *mc    = fRsnEvent[0].GetRefMCESD();
   AliStack      *stack = mc->Stack();
   TArrayI        indexes[2];
   Int_t          i, j, istep, imax, icheck, itrack[2], ipart;
   Int_t          pdg, label[2];
   Short_t        charge, pairDefMatch[2];
   TParticle     *part = 0x0;
   AliMCParticle *mother = 0x0;
   
   // set pointers
   fMother.SetDaughter(0, &fDaughter[0]);
   fMother.SetDaughter(1, &fDaughter[1]);
   
   // loop on definitions
   AliRsnPairDef *def = 0x0;
   TObjArrayIter nextDef(&fDefs);
   while ( (def = (AliRsnPairDef*)nextDef()) ) {
      
      // loop on the MC list of particles
      for (ipart = 0; ipart < stack->GetNprimary(); ipart++) {
         
         // reset daughters
         fDaughter[0].Reset();
         fDaughter[1].Reset();
         
         // taks MC particle
         mother = (AliMCParticle*)mc->GetTrack(ipart);
         
         // check that it is a binary decay and the PDG code matches
         if (mother->Particle()->GetNDaughters() != 2) continue;
         if (mother->Particle()->GetPdgCode() != def->GetMotherPDG()) continue;
         
         // store the labels of the two daughters
         // and check that they are in the stack
         label[0] = mother->Particle()->GetFirstDaughter();
         label[1] = mother->Particle()->GetLastDaughter();
         if (label[0] < 0 || label[0] > stack->GetNtrack()) continue;
         if (label[1] < 0 || label[1] > stack->GetNtrack()) continue;
         
         // for each daughter, check what slot in the pair definition it matches
         // or if it does not match any of them
         for (i = 0; i < 2; i++) {
            pairDefMatch[i] = -1;
            part = stack->Particle(label[i]);
            if (part) {
               pdg    = TMath::Abs(part->GetPdgCode());
               charge = (Short_t)(part->GetPDG()->Charge() / 3);
               if (def->GetDef1()->MatchesPDG(pdg) && def->GetDef1()->MatchesCharge(charge)) pairDefMatch[i] = 0;
               if (def->GetDef2()->MatchesPDG(pdg) && def->GetDef2()->MatchesCharge(charge)) pairDefMatch[i] = 1;
            }
         }
         
         // if the two label match the two definitions for the pair
         // and if they are in the wrong order, swap them,
         // otherwise, if they don't match the definition in any order, skip
         if (pairDefMatch[0] == 1 && pairDefMatch[1] == 0) {
            icheck = label[0];
            label[0] = label[1];
            label[1] = icheck;
         }
         else if (pairDefMatch[0] < 0 || pairDefMatch[1] < 0) continue;
         
         // from now on, we are sure that label[0] refers to the particle
         // that matches definitions of first daughter, and label[1] to
         // the particle that matches definitions of second daughter
         fDaughter[0].SetRefMC(mc->GetTrack(label[0]));
         fDaughter[1].SetRefMC(mc->GetTrack(label[1]));
         
         // assign masses and fill the MC steps,
         // where reconstruction is not taken into account
         fMother.ComputeSum(def->GetMass1(), def->GetMass2());
         FillContainer(kTRUE, def);
         
         // search for all reconstructed tracks which have these labels
         for (i = 0; i < 2; i++) indexes[i] = FindTracks(label[i], esd);
         
         // if not both tracks have been reconstructed, stop here
         if (indexes[0].GetSize() < 1 || indexes[1].GetSize() < 1) continue;
         
         // if both daughters were reconstructed
         // search for the best combination of indexes for this pair
         imax = itrack[0] = itrack[1] = 0;
         for (i = 0; i < indexes[0].GetSize(); i++) {
            for (j = 0; j < indexes[1].GetSize(); j++) {
               fDaughter[0].SetRef(esd->GetTrack(indexes[0][i]));
               fDaughter[1].SetRef(esd->GetTrack(indexes[1][j]));
               fMother.ComputeSum(def->GetMass1(), def->GetMass2());
               istep = NGoodSteps();
               if (istep > imax) {
                  itrack[0] = indexes[0][i];
                  itrack[1] = indexes[1][j];
                  imax = istep;
               }
            }
         }
         
         // then assign definitely the best combination and fill rec container
         fDaughter[0].SetRef(esd->GetTrack(itrack[0]));
         fDaughter[1].SetRef(esd->GetTrack(itrack[1]));
         fMother.ComputeSum(def->GetMass1(), def->GetMass2());
         FillContainer(kFALSE, def);
      }
   }
}

//_____________________________________________________________________________
void AliRsnAnalysisTaskEffPair::ProcessEventAOD()
{
//
// Process current event with the definitions of the specified step in MC list
// and store results in the container slot defined in second argument.
// It is associated with the AliCFContainer with the name of the pair.
//

   AliAODEvent  *aod     = fRsnEvent[0].GetRefAOD();
   TClonesArray *mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
   if (!mcArray) return;
   
   TArrayI           indexes[2];
   Int_t             i, j, pdg, imax, istep, icheck, itrack[2], label[2];
   AliAODMCParticle *part = 0x0, *mother = 0x0;
   Short_t           charge = 0, pairDefMatch[2];
   
   // set pointers
   fMother.SetDaughter(0, &fDaughter[0]);
   fMother.SetDaughter(1, &fDaughter[1]);
   
   // loop on definitions
   AliRsnPairDef *def = 0x0;
   TObjArrayIter nextDef(&fDefs);
   while ( (def = (AliRsnPairDef*)nextDef()) ) {
      
      // loop on the MC list of particles
      TObjArrayIter next(mcArray);
      while ((mother = (AliAODMCParticle*)next())) {
         
         // check that it is a binary decay and the PDG code matches
         if (mother->GetNDaughters() != 2) continue;
         if (mother->GetPdgCode() != def->GetMotherPDG()) continue;
         
         // store the labels of the two daughters
         label[0] = mother->GetDaughter(0);
         label[1] = mother->GetDaughter(1);
         
         // for each daughter, check what slot in the pair definition it matches
         // or if it does not match any of them
         for (i = 0; i < 2; i++) {
            pairDefMatch[i] = -1;
            part = (AliAODMCParticle*)mcArray->At(label[i]);
            if (part) {
               pdg    = TMath::Abs(part->GetPdgCode());
               charge = (Short_t)part->Charge();
               if (def->GetDef1()->MatchesPDG(pdg) && def->GetDef1()->MatchesCharge(charge)) pairDefMatch[i] = 0;
               if (def->GetDef2()->MatchesPDG(pdg) && def->GetDef2()->MatchesCharge(charge)) pairDefMatch[i] = 1;
            }
         }
         
         // if the two label match the two definitions for the pair
         // and if they are in the wrong order, swap them,
         // otherwise, if they don't match the definition in any order, skip
         if (pairDefMatch[0] == 1 && pairDefMatch[1] == 0) {
            icheck = label[0];
            label[0] = label[1];
            label[1] = icheck;
         }
         else if (pairDefMatch[0] < 0 || pairDefMatch[1] < 0) continue;
         
         // from now on, we are sure that label[0] refers to the particle
         // that matches definitions of first daughter, and label[1] to
         // the particle that matches definitions of second daughter
         fDaughter[0].SetRefMC((AliAODMCParticle*)mcArray->At(label[0]));
         fDaughter[1].SetRefMC((AliAODMCParticle*)mcArray->At(label[1]));
         
         // assign masses and fill the MC steps,
         // where reconstruction is not taken into account
         fMother.ComputeSum(def->GetMass1(), def->GetMass2());
         FillContainer(kTRUE, def);
         
         // search for all reconstructed tracks which have these labels
         for (i = 0; i < 2; i++) indexes[i] = FindTracks(label[i], aod);
         
         // if not both tracks have been reconstructed, stop here
         if (indexes[0].GetSize() < 1 || indexes[1].GetSize() < 1) continue;
         
         // if both daughters were reconstructed
         // search for the best combination of indexes for this pair
         imax = itrack[0] = itrack[1] = 0;
         for (i = 0; i < indexes[0].GetSize(); i++) {
            for (j = 0; j < indexes[1].GetSize(); j++) {
               fDaughter[0].SetRef(aod->GetTrack(indexes[0][i]));
               fDaughter[1].SetRef(aod->GetTrack(indexes[1][j]));
               fMother.ComputeSum(def->GetMass1(), def->GetMass2());
               istep = NGoodSteps();
               if (istep > imax) {
                  itrack[0] = indexes[0][i];
                  itrack[1] = indexes[1][j];
               }
            }
         }
         
         // then assign definitely the best combination and fill rec container
         fDaughter[0].SetRef(aod->GetTrack(itrack[0]));
         fDaughter[1].SetRef(aod->GetTrack(itrack[1]));
         fMother.ComputeSum(def->GetMass1(), def->GetMass2());
         FillContainer(kFALSE, def);
      }
   }
}

//_____________________________________________________________________________
void AliRsnAnalysisTaskEffPair::FillContainer(Bool_t mcList, TObject *defObj)
{
//
// Fill the container corresponding to current definition.
//
   
   // cast the object into the def
   if (!defObj->InheritsFrom(AliRsnPairDef::Class())) {
      AliError("Def object does not inherit from AliRsnPairDef!");
      return;
   }
   AliRsnPairDef *def = static_cast<AliRsnPairDef*>(defObj);
   
   // retrieve container
   AliCFContainer *cont = (AliCFContainer*)fOutList->FindObject(def->GetName());
   if (!cont) return;

   TObjArray &stepList =  (mcList ? fStepsMC : fStepsRec);
   Int_t      firstStep = (mcList ? 0 : ((Int_t)fStepsMC.GetEntries()));
   Int_t      iaxis, nAxes  = fAxes.GetEntries();
   Int_t      istep, nSteps = stepList.GetEntries();
   Bool_t     computeOK;
   
   // compute values for all axes
   for (iaxis = 0; iaxis < nAxes; iaxis++) {
      AliRsnValue *fcnAxis = (AliRsnValue*)fAxes.At(iaxis);
      fVar[iaxis] = -1E10;
      switch (fcnAxis->GetTargetType()) {
         case AliRsnTarget::kMother:
            fcnAxis->SetSupportObject(def);
            computeOK = fcnAxis->Eval(&fMother, mcList);
            break;
         case AliRsnTarget::kEvent:
            computeOK = fcnAxis->Eval(&fRsnEvent[0]);
            break;
         default:
            AliError(Form("Allowed targets are mothers and events; cannot use axis '%s' which has target '%s'", fcnAxis->GetName(), fcnAxis->GetTargetTypeName()));
            computeOK = kFALSE;
      }
      if (computeOK) fVar[iaxis] = ((Float_t)fcnAxis->GetComputedValue());
   }

   // fill all successful steps
   for (istep = 0; istep < nSteps; istep++) {
      AliRsnCutManager *cutMgr = (AliRsnCutManager*)stepList[istep];
      
      if (!cutMgr->PassCommonDaughterCuts(&fDaughter[0])) break;
      if (!cutMgr->PassCommonDaughterCuts(&fDaughter[1])) break;
      if (!cutMgr->PassDaughter1Cuts(&fDaughter[0])) break;
      if (!cutMgr->PassDaughter2Cuts(&fDaughter[1])) break;
      if (!cutMgr->PassMotherCuts(&fMother)) break;
      
      AliDebug(AliLog::kDebug + 2, Form("DEF: %s --> filling step %d", def->GetName(), istep));
      cont->Fill(fVar.GetArray(), istep + firstStep);
   }
}
Commit	Line	Data
99261e24	1	//
	2	// Class AliRsnAnalysisTaskEffPair
	3	//
	4	// Inherits from basic AliRsnAnalysisTaskEff for efficiency,
	5	// and computed efficiencies for pairs
	6	//
	7	// author: Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
	8	//
	9
	10	#include "AliStack.h"
	11	#include "AliMCEvent.h"
	12	#include "AliESDEvent.h"
	13	#include "AliAODEvent.h"
	14
	15	#include "AliRsnPairDef.h"
	16	#include "AliRsnCutManager.h"
	17	#include "AliRsnAnalysisTaskEffPair.h"
	18
	19	ClassImp(AliRsnAnalysisTaskEffPair)
	20
	21	//_____________________________________________________________________________
	22	AliRsnAnalysisTaskEffPair::AliRsnAnalysisTaskEffPair(const char *name) :
	23	AliRsnAnalysisTaskEff(name),
d7712d44	24	fMother()
99261e24	25	{
	26	//
	27	// Default constructor.
	28	// Do not repeat 'DefineOutput' since it is done in base class and we don't add new ones.
	29	//
	30	}
	31
	32	//_____________________________________________________________________________
	33	AliRsnAnalysisTaskEffPair::AliRsnAnalysisTaskEffPair(const AliRsnAnalysisTaskEffPair& copy) :
	34	AliRsnAnalysisTaskEff(copy),
d7712d44	35	fMother()
99261e24	36	{
	37	//
	38	// Copy constrtuctor.
	39	//
	40	}
	41
	42	//_____________________________________________________________________________
	43	AliRsnAnalysisTaskEffPair& AliRsnAnalysisTaskEffPair::operator=(const AliRsnAnalysisTaskEffPair& copy)
	44	{
	45	//
	46	// Assignment operator.
	47	// Owned data members are meaningless for this operator.
	48	//
	49
	50	AliRsnAnalysisTaskEff::operator=(copy);
	51	return (*this);
	52	}
	53
f4459279	54	//_____________________________________________________________________________
	55	Bool_t AliRsnAnalysisTaskEffPair::RsnEventProcess()
	56	{
	57	//
	58	// Preprocessing to the event
	59	// Calls the same function in base class plus somethin specific
	60	//
	61
	62	if (!AliRsnAnalysisTaskEff::RsnEventProcess()) return kFALSE;
	63
	64	// assign this event to the mother reference
	65	fMother.SetRefEvent(&fRsnEvent[0]);
	66	fDaughter[0].SetOwnerEvent(&fRsnEvent[0]);
	67	fDaughter[1].SetOwnerEvent(&fRsnEvent[0]);
	68
	69	return kTRUE;
	70	}
	71
99261e24	72	//_____________________________________________________________________________
	73	Int_t AliRsnAnalysisTaskEffPair::NGoodSteps()
	74	{
	75	//
	76	// Checks how many 'reconstruction' steps are passed by current daughter
	77	//
	78
	79	Int_t istep, count = 0;
	80	Int_t nSteps = fStepsRec.GetEntries();
	81
	82	for (istep = 0; istep < nSteps; istep++) {
	83	AliRsnCutManager cutMgr = (AliRsnCutManager)fStepsRec[istep];
99261e24	84
	85	if (!cutMgr->PassCommonDaughterCuts(&fDaughter[0])) break;
	86	if (!cutMgr->PassCommonDaughterCuts(&fDaughter[1])) break;
	87	if (!cutMgr->PassDaughter1Cuts(&fDaughter[0])) break;
	88	if (!cutMgr->PassDaughter2Cuts(&fDaughter[1])) break;
	89	if (!cutMgr->PassMotherCuts(&fMother)) break;
	90
	91	count++;
	92	}
	93
	94	return count;
	95	}
	96
	97	//_____________________________________________________________________________
	98	void AliRsnAnalysisTaskEffPair::ProcessEventESD()
	99	{
	100	//
	101	// Process current event with the definitions of the specified step in MC list
	102	// and store results in the container slot defined in second argument.
	103	// It is associated with the AliCFContainer with the name of the pair.
	104	//
	105
	106	AliESDEvent *esd = fRsnEvent[0].GetRefESD();
	107	AliMCEvent *mc = fRsnEvent[0].GetRefMCESD();
	108	AliStack *stack = mc->Stack();
	109	TArrayI indexes[2];
	110	Int_t i, j, istep, imax, icheck, itrack[2], ipart;
	111	Int_t pdg, label[2];
	112	Short_t charge, pairDefMatch[2];
	113	TParticle *part = 0x0;
	114	AliMCParticle *mother = 0x0;
	115
d7712d44	116	// set pointers
	117	fMother.SetDaughter(0, &fDaughter[0]);
	118	fMother.SetDaughter(1, &fDaughter[1]);
	119
99261e24	120	// loop on definitions
d7712d44	121	AliRsnPairDef *def = 0x0;
99261e24	122	TObjArrayIter nextDef(&fDefs);
d7712d44	123	while ( (def = (AliRsnPairDef*)nextDef()) ) {
99261e24	124
	125	// loop on the MC list of particles
	126	for (ipart = 0; ipart < stack->GetNprimary(); ipart++) {
	127
d7712d44	128	// reset daughters
	129	fDaughter[0].Reset();
	130	fDaughter[1].Reset();
	131
99261e24	132	// taks MC particle
	133	mother = (AliMCParticle*)mc->GetTrack(ipart);
	134
	135	// check that it is a binary decay and the PDG code matches
	136	if (mother->Particle()->GetNDaughters() != 2) continue;
d7712d44	137	if (mother->Particle()->GetPdgCode() != def->GetMotherPDG()) continue;
99261e24	138
	139	// store the labels of the two daughters
	140	// and check that they are in the stack
	141	label[0] = mother->Particle()->GetFirstDaughter();
	142	label[1] = mother->Particle()->GetLastDaughter();
	143	if (label[0] < 0 \|\| label[0] > stack->GetNtrack()) continue;
	144	if (label[1] < 0 \|\| label[1] > stack->GetNtrack()) continue;
	145
	146	// for each daughter, check what slot in the pair definition it matches
	147	// or if it does not match any of them
	148	for (i = 0; i < 2; i++) {
	149	pairDefMatch[i] = -1;
	150	part = stack->Particle(label[i]);
	151	if (part) {
	152	pdg = TMath::Abs(part->GetPdgCode());
	153	charge = (Short_t)(part->GetPDG()->Charge() / 3);
d7712d44	154	if (def->GetDef1()->MatchesPDG(pdg) && def->GetDef1()->MatchesCharge(charge)) pairDefMatch[i] = 0;
d7712d44	155	if (def->GetDef2()->MatchesPDG(pdg) && def->GetDef2()->MatchesCharge(charge)) pairDefMatch[i] = 1;
99261e24	156	}
	157	}
	158
	159	// if the two label match the two definitions for the pair
	160	// and if they are in the wrong order, swap them,
	161	// otherwise, if they don't match the definition in any order, skip
	162	if (pairDefMatch[0] == 1 && pairDefMatch[1] == 0) {
	163	icheck = label[0];
	164	label[0] = label[1];
	165	label[1] = icheck;
	166	}
	167	else if (pairDefMatch[0] < 0 \|\| pairDefMatch[1] < 0) continue;
	168
	169	// from now on, we are sure that label[0] refers to the particle
	170	// that matches definitions of first daughter, and label[1] to
	171	// the particle that matches definitions of second daughter
	172	fDaughter[0].SetRefMC(mc->GetTrack(label[0]));
	173	fDaughter[1].SetRefMC(mc->GetTrack(label[1]));
	174
	175	// assign masses and fill the MC steps,
	176	// where reconstruction is not taken into account
d7712d44	177	fMother.ComputeSum(def->GetMass1(), def->GetMass2());
d7712d44	178	FillContainer(kTRUE, def);
99261e24	179
	180	// search for all reconstructed tracks which have these labels
	181	for (i = 0; i < 2; i++) indexes[i] = FindTracks(label[i], esd);
	182
	183	// if not both tracks have been reconstructed, stop here
	184	if (indexes[0].GetSize() < 1 \|\| indexes[1].GetSize() < 1) continue;
	185
	186	// if both daughters were reconstructed
	187	// search for the best combination of indexes for this pair
	188	imax = itrack[0] = itrack[1] = 0;
	189	for (i = 0; i < indexes[0].GetSize(); i++) {
	190	for (j = 0; j < indexes[1].GetSize(); j++) {
	191	fDaughter[0].SetRef(esd->GetTrack(indexes[0][i]));
	192	fDaughter[1].SetRef(esd->GetTrack(indexes[1][j]));
d7712d44	193	fMother.ComputeSum(def->GetMass1(), def->GetMass2());
99261e24	194	istep = NGoodSteps();
	195	if (istep > imax) {
	196	itrack[0] = indexes[0][i];
	197	itrack[1] = indexes[1][j];
	198	imax = istep;
	199	}
	200	}
	201	}
	202
	203	// then assign definitely the best combination and fill rec container
	204	fDaughter[0].SetRef(esd->GetTrack(itrack[0]));
	205	fDaughter[1].SetRef(esd->GetTrack(itrack[1]));
d7712d44	206	fMother.ComputeSum(def->GetMass1(), def->GetMass2());
d7712d44	207	FillContainer(kFALSE, def);
99261e24	208	}
	209	}
	210	}
	211
	212	//_____________________________________________________________________________
	213	void AliRsnAnalysisTaskEffPair::ProcessEventAOD()
	214	{
	215	//
	216	// Process current event with the definitions of the specified step in MC list
	217	// and store results in the container slot defined in second argument.
	218	// It is associated with the AliCFContainer with the name of the pair.
	219	//
	220
	221	AliAODEvent *aod = fRsnEvent[0].GetRefAOD();
	222	TClonesArray mcArray = (TClonesArray)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
	223	if (!mcArray) return;
	224
	225	TArrayI indexes[2];
	226	Int_t i, j, pdg, imax, istep, icheck, itrack[2], label[2];
	227	AliAODMCParticle part = 0x0, mother = 0x0;
	228	Short_t charge = 0, pairDefMatch[2];
	229
d7712d44	230	// set pointers
	231	fMother.SetDaughter(0, &fDaughter[0]);
	232	fMother.SetDaughter(1, &fDaughter[1]);
	233
99261e24	234	// loop on definitions
d7712d44	235	AliRsnPairDef *def = 0x0;
99261e24	236	TObjArrayIter nextDef(&fDefs);
d7712d44	237	while ( (def = (AliRsnPairDef*)nextDef()) ) {
99261e24	238
	239	// loop on the MC list of particles
	240	TObjArrayIter next(mcArray);
	241	while ((mother = (AliAODMCParticle*)next())) {
	242
	243	// check that it is a binary decay and the PDG code matches
	244	if (mother->GetNDaughters() != 2) continue;
d7712d44	245	if (mother->GetPdgCode() != def->GetMotherPDG()) continue;
99261e24	246
	247	// store the labels of the two daughters
	248	label[0] = mother->GetDaughter(0);
	249	label[1] = mother->GetDaughter(1);
	250
	251	// for each daughter, check what slot in the pair definition it matches
	252	// or if it does not match any of them
	253	for (i = 0; i < 2; i++) {
	254	pairDefMatch[i] = -1;
	255	part = (AliAODMCParticle*)mcArray->At(label[i]);
	256	if (part) {
	257	pdg = TMath::Abs(part->GetPdgCode());
	258	charge = (Short_t)part->Charge();
d7712d44	259	if (def->GetDef1()->MatchesPDG(pdg) && def->GetDef1()->MatchesCharge(charge)) pairDefMatch[i] = 0;
d7712d44	260	if (def->GetDef2()->MatchesPDG(pdg) && def->GetDef2()->MatchesCharge(charge)) pairDefMatch[i] = 1;
99261e24	261	}
	262	}
	263
	264	// if the two label match the two definitions for the pair
	265	// and if they are in the wrong order, swap them,
	266	// otherwise, if they don't match the definition in any order, skip
	267	if (pairDefMatch[0] == 1 && pairDefMatch[1] == 0) {
	268	icheck = label[0];
	269	label[0] = label[1];
	270	label[1] = icheck;
	271	}
	272	else if (pairDefMatch[0] < 0 \|\| pairDefMatch[1] < 0) continue;
	273
	274	// from now on, we are sure that label[0] refers to the particle
	275	// that matches definitions of first daughter, and label[1] to
	276	// the particle that matches definitions of second daughter
	277	fDaughter[0].SetRefMC((AliAODMCParticle*)mcArray->At(label[0]));
	278	fDaughter[1].SetRefMC((AliAODMCParticle*)mcArray->At(label[1]));
	279
	280	// assign masses and fill the MC steps,
	281	// where reconstruction is not taken into account
d7712d44	282	fMother.ComputeSum(def->GetMass1(), def->GetMass2());
d7712d44	283	FillContainer(kTRUE, def);
99261e24	284
	285	// search for all reconstructed tracks which have these labels
	286	for (i = 0; i < 2; i++) indexes[i] = FindTracks(label[i], aod);
	287
	288	// if not both tracks have been reconstructed, stop here
	289	if (indexes[0].GetSize() < 1 \|\| indexes[1].GetSize() < 1) continue;
	290
	291	// if both daughters were reconstructed
	292	// search for the best combination of indexes for this pair
	293	imax = itrack[0] = itrack[1] = 0;
	294	for (i = 0; i < indexes[0].GetSize(); i++) {
	295	for (j = 0; j < indexes[1].GetSize(); j++) {
	296	fDaughter[0].SetRef(aod->GetTrack(indexes[0][i]));
	297	fDaughter[1].SetRef(aod->GetTrack(indexes[1][j]));
d7712d44	298	fMother.ComputeSum(def->GetMass1(), def->GetMass2());
99261e24	299	istep = NGoodSteps();
	300	if (istep > imax) {
	301	itrack[0] = indexes[0][i];
	302	itrack[1] = indexes[1][j];
	303	}
	304	}
	305	}
	306
	307	// then assign definitely the best combination and fill rec container
	308	fDaughter[0].SetRef(aod->GetTrack(itrack[0]));
	309	fDaughter[1].SetRef(aod->GetTrack(itrack[1]));
d7712d44	310	fMother.ComputeSum(def->GetMass1(), def->GetMass2());
d7712d44	311	FillContainer(kFALSE, def);
99261e24	312	}
	313	}
	314	}
	315
	316	//_____________________________________________________________________________
d7712d44	317	void AliRsnAnalysisTaskEffPair::FillContainer(Bool_t mcList, TObject *defObj)
99261e24	318	{
	319	//
	320	// Fill the container corresponding to current definition.
	321	//
d7712d44	322
	323	// cast the object into the def
	324	if (!defObj->InheritsFrom(AliRsnPairDef::Class())) {
	325	AliError("Def object does not inherit from AliRsnPairDef!");
	326	return;
	327	}
	328	AliRsnPairDef def = static_cast<AliRsnPairDef>(defObj);
	329
99261e24	330	// retrieve container
d7712d44	331	AliCFContainer cont = (AliCFContainer)fOutList->FindObject(def->GetName());
99261e24	332	if (!cont) return;
	333
	334	TObjArray &stepList = (mcList ? fStepsMC : fStepsRec);
	335	Int_t firstStep = (mcList ? 0 : ((Int_t)fStepsMC.GetEntries()));
	336	Int_t iaxis, nAxes = fAxes.GetEntries();
	337	Int_t istep, nSteps = stepList.GetEntries();
	338	Bool_t computeOK;
	339
	340	// compute values for all axes
	341	for (iaxis = 0; iaxis < nAxes; iaxis++) {
	342	AliRsnValue fcnAxis = (AliRsnValue)fAxes.At(iaxis);
	343	fVar[iaxis] = -1E10;
	344	switch (fcnAxis->GetTargetType()) {
	345	case AliRsnTarget::kMother:
d7712d44	346	fcnAxis->SetSupportObject(def);
99261e24	347	computeOK = fcnAxis->Eval(&fMother, mcList);
	348	break;
	349	case AliRsnTarget::kEvent:
	350	computeOK = fcnAxis->Eval(&fRsnEvent[0]);
	351	break;
	352	default:
	353	AliError(Form("Allowed targets are mothers and events; cannot use axis '%s' which has target '%s'", fcnAxis->GetName(), fcnAxis->GetTargetTypeName()));
	354	computeOK = kFALSE;
	355	}
	356	if (computeOK) fVar[iaxis] = ((Float_t)fcnAxis->GetComputedValue());
	357	}
	358
	359	// fill all successful steps
	360	for (istep = 0; istep < nSteps; istep++) {
	361	AliRsnCutManager cutMgr = (AliRsnCutManager)stepList[istep];
99261e24	362
	363	if (!cutMgr->PassCommonDaughterCuts(&fDaughter[0])) break;
	364	if (!cutMgr->PassCommonDaughterCuts(&fDaughter[1])) break;
	365	if (!cutMgr->PassDaughter1Cuts(&fDaughter[0])) break;
	366	if (!cutMgr->PassDaughter2Cuts(&fDaughter[1])) break;
	367	if (!cutMgr->PassMotherCuts(&fMother)) break;
	368
d7712d44	369	AliDebug(AliLog::kDebug + 2, Form("DEF: %s --> filling step %d", def->GetName(), istep));
99261e24	370	cont->Fill(fVar.GetArray(), istep + firstStep);
	371	}
	372	}