[u/mrichter/AliRoot.git] / PWG3 / dielectron / AliDielectronMC.cxx

/*************************************************************************
* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

//#####################################################
//#                                                   # 
//#              Class AliDielectronMC                #
//#       Cut Class for Jpsi->e+e- analysis           #
//#                                                   #
//#   by WooJin J. Park, GSI / W.J.Park@gsi.de        #
//#                                                   #
//#####################################################

#include <AliAnalysisManager.h>
#include <AliAODHandler.h>
#include <AliESDInputHandler.h>
#include <AliMCEventHandler.h>
#include <AliMCEvent.h>
#include <AliMCParticle.h>
#include <AliAODMCParticle.h>
#include <AliStack.h>
#include <AliESDEvent.h>
#include <AliESDtrack.h>
#include <AliLog.h>

#include "AliDielectronMC.h"

AliDielectronMC* AliDielectronMC::fgInstance=0x0;

//____________________________________________________________
AliDielectronMC* AliDielectronMC::Instance()
{
  //
  // return pointer to singleton implementation
  //
  if (fgInstance) return fgInstance;
  
  AnalysisType type=kUNSET;
  Bool_t hasMC=kFALSE;
  if (AliAnalysisManager::GetAnalysisManager()){
    if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD;
    else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODHandler::Class()) type=kAOD;

    AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
    hasMC=mcHandler!=0x0;
  }
  
  fgInstance=new AliDielectronMC(type);
  fgInstance->SetHasMC(hasMC);
  
  return fgInstance;
}

//____________________________________________________________
AliDielectronMC::AliDielectronMC(AnalysisType type):
  fMCEvent(0x0),
  fStack(0x0),
  fAnaType(type),
  fHasMC(kTRUE)
{
  //
  // default constructor
  //
}


//____________________________________________________________
AliDielectronMC::~AliDielectronMC()
{
  //
  // default destructor
  //
  
}

//____________________________________________________________
void AliDielectronMC::Initialize()
{
  //
  // initialize MC class
  //
  if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
}

//____________________________________________________________
Int_t AliDielectronMC::GetNMCTracks()
{
  //
  //  return the number of generated tracks from MC event
  //
  if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
  return fMCEvent->GetNumberOfTracks();
}

//____________________________________________________________
Int_t AliDielectronMC::GetNMCTracksFromStack()
{
  //
  //  return the number of generated tracks from stack
  //
  if (!fStack){ AliError("No fStack"); return -999; }
  return fStack->GetNtrack();
}

//____________________________________________________________
AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t itrk)
{
  //
  // return MC track directly from MC event
  //
  if (itrk<0) return NULL;
  if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
  AliVParticle * track = fMCEvent->GetTrack(itrk); //  tracks from MC event
  return track;
}

//____________________________________________________________
Bool_t AliDielectronMC::ConnectMCEvent()
{
  //
  // connect stack object from the mc handler
  //
  fMCEvent=0x0;
  AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
  if (!mcHandler){ AliError("Could not retrive MC event handler!"); return kFALSE; }
  if (!mcHandler->InitOk() ) return kFALSE;
  if (!mcHandler->TreeK() )  return kFALSE;
  if (!mcHandler->TreeTR() ) return kFALSE;
  
  AliMCEvent* mcEvent = mcHandler->MCEvent();
  if (!mcEvent){ AliError("Could not retrieve MC event!"); return kFALSE; }
  fMCEvent = mcEvent;
  
  if (!UpdateStack()) return kFALSE;
  return kTRUE;
}

//____________________________________________________________
Bool_t AliDielectronMC::UpdateStack()
{
  //
  // update stack with new event
  //
  if (!fMCEvent){ AliError("No fMCEvent"); return kFALSE;}
  AliStack* stack = fMCEvent->Stack();
  if (!stack){ AliError("Could not retrive stack!"); return kFALSE; }
  fStack = stack;
  return kTRUE;
}

//____________________________________________________________
AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
{
  //
  // return MC track
  //
  if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}

  Int_t nStack = fMCEvent->GetNumberOfTracks();
  Int_t label = TMath::Abs(_track->GetLabel());
  if(label>nStack)return NULL;

  AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
  return mctrack;
}

//____________________________________________________________
TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
{
  //
  // return MC track from stack
  //
  Int_t label = TMath::Abs(_track->GetLabel());
  if (!fStack) AliWarning("fStack is not available. Update stack first.");
  TParticle* mcpart = fStack->Particle(label);
  if (!mcpart) return NULL;
  return mcpart;
}

//____________________________________________________________
AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track)
{
  //
  // return MC track mother
  //
  AliMCParticle* mcpart = GetMCTrack(_track);
  if (!mcpart) return NULL;
  AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
  if (!mcmother) return NULL;
  return mcmother;
}

//____________________________________________________________
AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
  //
  // return MC track mother
  //
  AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
  return mcmother;
}

//____________________________________________________________
AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
  //
  // return MC track mother
  //
  AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
  return mcmother;
}

//____________________________________________________________
TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
{
  //
  // return MC track mother from stack
  //
  TParticle* mcpart = GetMCTrackFromStack(_track);
  if ( !mcpart || mcpart->GetFirstMother()<=0 ) return NULL; 
  TParticle* mcmother = fStack->Particle(mcpart->GetFirstMother());
  if (!mcmother) return NULL;
  return mcmother;
}

//____________________________________________________________
Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
{
  //
  // return PDG code of the track from the MC truth info
  //
  AliMCParticle* mcpart = GetMCTrack(_track);
  if (!mcpart) return -999;
  return mcpart->PdgCode();
}

//____________________________________________________________
Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
{
  // 
  // return MC PDG code from stack
  //
  TParticle* mcpart = GetMCTrackFromStack(_track);
  if (!mcpart) return -999;
  return mcpart->GetPdgCode();
}

//____________________________________________________________
Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
{
  //
  // return PDG code of the mother track from the MC truth info
  //
  AliMCParticle* mcmother = GetMCTrackMother(_track);
  if (!mcmother) return -999;
  return mcmother->PdgCode();
}

//____________________________________________________________
Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
{
  //
  // return PDG code of the mother track from stack
  //
  TParticle* mcmother = GetMCTrackMotherFromStack(_track);
  if (!mcmother) return -999;
  return mcmother->GetPdgCode();
}

//____________________________________________________________
Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
{
  //
  // return process number of the track
  //
  AliMCParticle* mcpart = GetMCTrack(_track);
  if (!mcpart) return -999;
  return 0;
}

//____________________________________________________________
Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
{
  //
  // return process number of the track
  //
  TParticle* mcpart = GetMCTrackFromStack(_track);
  if (!mcpart) return -999;
  return mcpart->GetUniqueID();
}

//____________________________________________________________
Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
{
  //
  // returns the number of daughters
  //
  AliMCParticle *mcmother=GetMCTrackMother(track);
  if(!mcmother||!mcmother->Particle()) return -999;
//   return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
  return mcmother->Particle()->GetNDaughters();
}

//____________________________________________________________
Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle)
{
  //
  // returns the number of daughters
  //
  AliMCParticle *mcmother=GetMCTrackMother(particle);
  if(!mcmother||!mcmother->Particle()) return -999;
//   return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
  return mcmother->Particle()->GetNDaughters();
}

//____________________________________________________________
Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle)
{
  //
  // returns the number of daughters
  //
  AliAODMCParticle *mcmother=GetMCTrackMother(particle);
  if(!mcmother) return -999;
  return mcmother->GetNDaughters();
}

//____________________________________________________________
Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
{
  //
  // return process number of the mother of the track
  //
  AliMCParticle* mcmother = GetMCTrackMother(_track);
  if (!mcmother) return -999;
  return 0;
}

//____________________________________________________________
Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
{
  //
  // return process number of the mother of the track
  //
  TParticle* mcmother = GetMCTrackMotherFromStack(_track);
  if (!mcmother) return -999;
  return mcmother->GetUniqueID();
}

//____________________________________________________________
Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother)
{
  //
  // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
  //
  
  if (!fMCEvent) return kFALSE;
  if (!particle) return kFALSE;
  
  if (particle->IsA()==AliMCParticle::Class()){
    return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
  } else if (particle->IsA()==AliAODMCParticle::Class()){
    return IsMCMotherToEEaod(static_cast<const AliAODMCParticle*>(particle),pdgMother);
  } else {
    AliError("Unknown particle type");
  }
  return kFALSE;
  
}

//____________________________________________________________
Bool_t AliDielectronMC::IsMCMotherToEEesd(const AliMCParticle *particle, Int_t pdgMother)
{
  //
  // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
  // ESD case
  //
  
  //check pdg code
  if (particle->PdgCode()!=pdgMother) return kFALSE;
  Int_t ifirst = particle->GetFirstDaughter();
  Int_t ilast  = particle->GetLastDaughter();
  
  //check number of daughters
  if ((ilast-ifirst)!=1) return kFALSE;
  AliMCParticle *firstD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ifirst));
  AliMCParticle *secondD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ilast));

  //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
  if (firstD->Charge()>0){
    if (firstD->PdgCode()!=-11) return kFALSE;
    if (secondD->PdgCode()!=11) return kFALSE;
  }else{
    if (firstD->PdgCode()!=11) return kFALSE;
    if (secondD->PdgCode()!=-11) return kFALSE;
  }
  
  return kTRUE;
}

//____________________________________________________________
Bool_t AliDielectronMC::IsMCMotherToEEaod(const AliAODMCParticle *particle, Int_t pdgMother)
{
  //
  // Check if the Mother 'particle' is of type pdgMother and decays to e+e-
  // AOD case
  //
  if (particle->GetPdgCode()!=pdgMother) return kFALSE;
  if (particle->GetNDaughters()!=2) return kFALSE;
  
  Int_t ifirst = particle->GetDaughter(0);
  Int_t ilast  = particle->GetDaughter(1);
  
  //check number of daughters
  if ((ilast-ifirst)!=1) return kFALSE;
  AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
  AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ilast));
  
  //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
  if (firstD->Charge()>0){
    if (firstD->GetPdgCode()!=11) return kFALSE;
    if (secondD->GetPdgCode()!=-11) return kFALSE;
  }else{
    if (firstD->GetPdgCode()!=-11) return kFALSE;
    if (secondD->GetPdgCode()!=11) return kFALSE;
  }
  return kTRUE;
}

//____________________________________________________________
Int_t AliDielectronMC::GetLabelMotherWithPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
{
  //
  // test if mother of particle 1 and 2 has pdgCode pdgMother and is the same;
  //
  if (!fMCEvent) return -1;
  
  if (fAnaType==kESD) return GetLabelMotherWithPdgESD(particle1, particle2, pdgMother);
  else if (fAnaType==kAOD) return GetLabelMotherWithPdgAOD(particle1, particle2, pdgMother);
  
  return -1;
}

//____________________________________________________________
Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
{
  //
  // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
  //    have the same mother and the mother had pdg code pdgMother
  // ESD case
  //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
  //
  
  AliMCParticle *mcPart1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
  AliMCParticle *mcPart2=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
  
  if (!mcPart1||!mcPart2) return -1;
  
  Int_t lblMother1=mcPart1->GetMother();
  Int_t lblMother2=mcPart2->GetMother();
  
  AliMCParticle *mcMother1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
  if (!mcMother1) return -1;
  if (lblMother1!=lblMother2) return -1;
  if (TMath::Abs(mcPart1->PdgCode())!=11) return -1;
  if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1;
  if (mcMother1->PdgCode()!=pdgMother) return -1;
  
  return lblMother1;
}

//____________________________________________________________
Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
{
  //
  // test if mother of particle 1 and 2 has pdgCode +-11 (electron),
  //    have the same mother and the mother had pdg code pdgMother
  // AOD case
  //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
  //
  AliAODMCParticle *mcPart1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
  AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
  
  if (!mcPart1||!mcPart2) return -1;
  
  Int_t lblMother1=mcPart1->GetMother();
  Int_t lblMother2=mcPart2->GetMother();
  
  AliAODMCParticle *mcMother1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblMother1));
  
  if (!mcMother1) return -1;
  if (lblMother1!=lblMother2) return -1;
  if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1;
  if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1;
  if (mcMother1->GetPdgCode()!=pdgMother) return -1;
  
  return lblMother1;
}

//____________________________________________________________
void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, AliVParticle* &d2)
{
  //
  // Get First two daughters of the mother
  //

  Int_t lblD1=-1;
  Int_t lblD2=-1;
  d1=0;
  d2=0;
  if (!fMCEvent) return;
  if (fAnaType==kAOD){
    const AliAODMCParticle *aodMother=static_cast<const AliAODMCParticle*>(mother);
    lblD1=aodMother->GetDaughter(0);
    lblD2=aodMother->GetDaughter(1);
  } else if (fAnaType==kESD){
    const AliMCParticle *aodMother=static_cast<const AliMCParticle*>(mother);
    lblD1=aodMother->GetFirstDaughter();
    lblD2=aodMother->GetLastDaughter();
  }
  d1=fMCEvent->GetTrack(lblD1);
  d2=fMCEvent->GetTrack(lblD2);
}

//____________________________________________________________
Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair)
{
  //
  // Check whether two particles have the same mother
  //

  const AliVParticle * daughter1 = pair->GetFirstDaughter();
  const AliVParticle * daughter2 = pair->GetSecondDaughter();

  AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
  AliVParticle *mcDaughter2=GetMCTrackFromMCEvent(daughter2->GetLabel());
  if (!mcDaughter1 || !mcDaughter2) return 0;

  Int_t labelMother1=-1;
  Int_t labelMother2=-1;

  if (mcDaughter1->IsA()==AliMCParticle::Class()){
    labelMother1=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
    labelMother2=(static_cast<AliMCParticle*>(mcDaughter2))->GetMother();
  } else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
    labelMother1=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
    labelMother2=(static_cast<AliAODMCParticle*>(mcDaughter2))->GetMother();
  }

  Bool_t sameMother=(labelMother1>-1)&&(labelMother2>-1)&&(labelMother1==labelMother2);

  return sameMother;
}
Commit	Line	Data
b2a297fa	1	/*************************************************************************
	2	* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//#####################################################
	17	//# #
	18	//# Class AliDielectronMC #
	19	//# Cut Class for Jpsi->e+e- analysis #
	20	//# #
	21	//# by WooJin J. Park, GSI / W.J.Park@gsi.de #
	22	//# #
	23	//#####################################################
	24
	25	#include <AliAnalysisManager.h>
	26	#include <AliAODHandler.h>
	27	#include <AliESDInputHandler.h>
	28	#include <AliMCEventHandler.h>
	29	#include <AliMCEvent.h>
	30	#include <AliMCParticle.h>
	31	#include <AliAODMCParticle.h>
	32	#include <AliStack.h>
	33	#include <AliESDEvent.h>
	34	#include <AliESDtrack.h>
	35	#include <AliLog.h>
	36
	37	#include "AliDielectronMC.h"
	38
	39	AliDielectronMC* AliDielectronMC::fgInstance=0x0;
	40
	41	//____________________________________________________________
	42	AliDielectronMC* AliDielectronMC::Instance()
	43	{
	44	//
	45	// return pointer to singleton implementation
	46	//
	47	if (fgInstance) return fgInstance;
	48
	49	AnalysisType type=kUNSET;
3505bfad	50	Bool_t hasMC=kFALSE;
	51	if (AliAnalysisManager::GetAnalysisManager()){
	52	if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD;
	53	else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODHandler::Class()) type=kAOD;
8df8e382	54
3505bfad	55	AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
	56	hasMC=mcHandler!=0x0;
	57	}
	58
b2a297fa	59	fgInstance=new AliDielectronMC(type);
3505bfad	60	fgInstance->SetHasMC(hasMC);
b2a297fa	61
	62	return fgInstance;
	63	}
	64
	65	//____________________________________________________________
	66	AliDielectronMC::AliDielectronMC(AnalysisType type):
	67	fMCEvent(0x0),
	68	fStack(0x0),
572b0139	69	fAnaType(type),
572b0139	70	fHasMC(kTRUE)
b2a297fa	71	{
	72	//
	73	// default constructor
	74	//
	75	}
	76
	77
	78	//____________________________________________________________
	79	AliDielectronMC::~AliDielectronMC()
	80	{
	81	//
	82	// default destructor
	83	//
	84
	85	}
	86
	87	//____________________________________________________________
	88	void AliDielectronMC::Initialize()
	89	{
	90	//
	91	// initialize MC class
	92	//
	93	if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
	94	}
	95
	96	//____________________________________________________________
	97	Int_t AliDielectronMC::GetNMCTracks()
	98	{
	99	//
	100	// return the number of generated tracks from MC event
	101	//
	102	if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
	103	return fMCEvent->GetNumberOfTracks();
	104	}
	105
	106	//____________________________________________________________
	107	Int_t AliDielectronMC::GetNMCTracksFromStack()
	108	{
	109	//
	110	// return the number of generated tracks from stack
	111	//
	112	if (!fStack){ AliError("No fStack"); return -999; }
	113	return fStack->GetNtrack();
	114	}
	115
	116	//____________________________________________________________
83e37742	117	AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t itrk)
b2a297fa	118	{
	119	//
	120	// return MC track directly from MC event
	121	//
83e37742	122	if (itrk<0) return NULL;
b2a297fa	123	if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
83e37742	124	AliVParticle * track = fMCEvent->GetTrack(itrk); // tracks from MC event
b2a297fa	125	return track;
	126	}
	127
	128	//____________________________________________________________
	129	Bool_t AliDielectronMC::ConnectMCEvent()
	130	{
	131	//
	132	// connect stack object from the mc handler
	133	//
8df8e382	134	fMCEvent=0x0;
b2a297fa	135	AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
b2a297fa	136	if (!mcHandler){ AliError("Could not retrive MC event handler!"); return kFALSE; }
8df8e382	137	if (!mcHandler->InitOk() ) return kFALSE;
	138	if (!mcHandler->TreeK() ) return kFALSE;
	139	if (!mcHandler->TreeTR() ) return kFALSE;
b2a297fa	140
	141	AliMCEvent* mcEvent = mcHandler->MCEvent();
	142	if (!mcEvent){ AliError("Could not retrieve MC event!"); return kFALSE; }
	143	fMCEvent = mcEvent;
	144
	145	if (!UpdateStack()) return kFALSE;
	146	return kTRUE;
	147	}
	148
	149	//____________________________________________________________
	150	Bool_t AliDielectronMC::UpdateStack()
	151	{
	152	//
	153	// update stack with new event
	154	//
	155	if (!fMCEvent){ AliError("No fMCEvent"); return kFALSE;}
	156	AliStack* stack = fMCEvent->Stack();
	157	if (!stack){ AliError("Could not retrive stack!"); return kFALSE; }
	158	fStack = stack;
	159	return kTRUE;
	160	}
	161
	162	//____________________________________________________________
8df8e382	163	AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
b2a297fa	164	{
	165	//
	166	// return MC track
	167	//
9143d69f	168	if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
	169
	170	Int_t nStack = fMCEvent->GetNumberOfTracks();
b2a297fa	171	Int_t label = TMath::Abs(_track->GetLabel());
9143d69f	172	if(label>nStack)return NULL;
9143d69f	173
b2a297fa	174	AliMCParticle mctrack = dynamic_cast<AliMCParticle >(fMCEvent->GetTrack(label));
	175	return mctrack;
	176	}
	177
	178	//____________________________________________________________
8df8e382	179	TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
b2a297fa	180	{
	181	//
	182	// return MC track from stack
	183	//
	184	Int_t label = TMath::Abs(_track->GetLabel());
	185	if (!fStack) AliWarning("fStack is not available. Update stack first.");
	186	TParticle* mcpart = fStack->Particle(label);
	187	if (!mcpart) return NULL;
	188	return mcpart;
	189	}
	190
	191	//____________________________________________________________
8df8e382	192	AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track)
b2a297fa	193	{
	194	//
	195	// return MC track mother
	196	//
	197	AliMCParticle* mcpart = GetMCTrack(_track);
	198	if (!mcpart) return NULL;
b2a297fa	199	AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
	200	if (!mcmother) return NULL;
	201	return mcmother;
	202	}
	203
	204	//____________________________________________________________
8df8e382	205	AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
	206	//
	207	// return MC track mother
	208	//
	209	AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
	210	return mcmother;
	211	}
	212
	213	//____________________________________________________________
	214	AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
	215	//
	216	// return MC track mother
	217	//
	218	AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
	219	return mcmother;
	220	}
	221
	222	//____________________________________________________________
	223	TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
b2a297fa	224	{
	225	//
	226	// return MC track mother from stack
	227	//
	228	TParticle* mcpart = GetMCTrackFromStack(_track);
	229	if ( !mcpart \|\| mcpart->GetFirstMother()<=0 ) return NULL;
	230	TParticle* mcmother = fStack->Particle(mcpart->GetFirstMother());
	231	if (!mcmother) return NULL;
	232	return mcmother;
	233	}
	234
	235	//____________________________________________________________
8df8e382	236	Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
b2a297fa	237	{
	238	//
	239	// return PDG code of the track from the MC truth info
	240	//
	241	AliMCParticle* mcpart = GetMCTrack(_track);
	242	if (!mcpart) return -999;
	243	return mcpart->PdgCode();
	244	}
	245
	246	//____________________________________________________________
8df8e382	247	Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
b2a297fa	248	{
	249	//
	250	// return MC PDG code from stack
	251	//
	252	TParticle* mcpart = GetMCTrackFromStack(_track);
	253	if (!mcpart) return -999;
	254	return mcpart->GetPdgCode();
	255	}
	256
	257	//____________________________________________________________
8df8e382	258	Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
b2a297fa	259	{
	260	//
	261	// return PDG code of the mother track from the MC truth info
	262	//
	263	AliMCParticle* mcmother = GetMCTrackMother(_track);
	264	if (!mcmother) return -999;
	265	return mcmother->PdgCode();
	266	}
	267
	268	//____________________________________________________________
8df8e382	269	Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
b2a297fa	270	{
	271	//
	272	// return PDG code of the mother track from stack
	273	//
	274	TParticle* mcmother = GetMCTrackMotherFromStack(_track);
	275	if (!mcmother) return -999;
	276	return mcmother->GetPdgCode();
	277	}
	278
	279	//____________________________________________________________
8df8e382	280	Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
b2a297fa	281	{
	282	//
	283	// return process number of the track
	284	//
	285	AliMCParticle* mcpart = GetMCTrack(_track);
	286	if (!mcpart) return -999;
	287	return 0;
	288	}
	289
	290	//____________________________________________________________
8df8e382	291	Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
b2a297fa	292	{
	293	//
	294	// return process number of the track
	295	//
	296	TParticle* mcpart = GetMCTrackFromStack(_track);
	297	if (!mcpart) return -999;
	298	return mcpart->GetUniqueID();
	299	}
	300
	301	//____________________________________________________________
8df8e382	302	Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
	303	{
	304	//
	305	// returns the number of daughters
	306	//
	307	AliMCParticle *mcmother=GetMCTrackMother(track);
	308	if(!mcmother\|\|!mcmother->Particle()) return -999;
	309	// return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
	310	return mcmother->Particle()->GetNDaughters();
	311	}
	312
	313	//____________________________________________________________
	314	Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle)
	315	{
	316	//
	317	// returns the number of daughters
	318	//
	319	AliMCParticle *mcmother=GetMCTrackMother(particle);
	320	if(!mcmother\|\|!mcmother->Particle()) return -999;
	321	// return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
	322	return mcmother->Particle()->GetNDaughters();
	323	}
	324
	325	//____________________________________________________________
	326	Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle)
	327	{
	328	//
	329	// returns the number of daughters
	330	//
	331	AliAODMCParticle *mcmother=GetMCTrackMother(particle);
	332	if(!mcmother) return -999;
	333	return mcmother->GetNDaughters();
	334	}
	335
	336	//____________________________________________________________
	337	Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
b2a297fa	338	{
	339	//
	340	// return process number of the mother of the track
	341	//
	342	AliMCParticle* mcmother = GetMCTrackMother(_track);
	343	if (!mcmother) return -999;
	344	return 0;
	345	}
	346
	347	//____________________________________________________________
8df8e382	348	Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
b2a297fa	349	{
	350	//
	351	// return process number of the mother of the track
	352	//
	353	TParticle* mcmother = GetMCTrackMotherFromStack(_track);
	354	if (!mcmother) return -999;
	355	return mcmother->GetUniqueID();
	356	}
	357
	358	//____________________________________________________________
	359	Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother)
	360	{
	361	//
	362	// Check if the Mother 'particle' is of type pdgMother and decays to e+e-
	363	//
	364
	365	if (!fMCEvent) return kFALSE;
8df8e382	366	if (!particle) return kFALSE;
b2a297fa	367
	368	if (particle->IsA()==AliMCParticle::Class()){
	369	return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
	370	} else if (particle->IsA()==AliAODMCParticle::Class()){
	371	return IsMCMotherToEEaod(static_cast<const AliAODMCParticle*>(particle),pdgMother);
	372	} else {
	373	AliError("Unknown particle type");
	374	}
	375	return kFALSE;
	376
	377	}
	378
	379	//____________________________________________________________
	380	Bool_t AliDielectronMC::IsMCMotherToEEesd(const AliMCParticle *particle, Int_t pdgMother)
	381	{
	382	//
	383	// Check if the Mother 'particle' is of type pdgMother and decays to e+e-
	384	// ESD case
	385	//
	386
	387	//check pdg code
	388	if (particle->PdgCode()!=pdgMother) return kFALSE;
	389	Int_t ifirst = particle->GetFirstDaughter();
	390	Int_t ilast = particle->GetLastDaughter();
	391
	392	//check number of daughters
	393	if ((ilast-ifirst)!=1) return kFALSE;
	394	AliMCParticle firstD=static_cast<AliMCParticle>(GetMCTrackFromMCEvent(ifirst));
	395	AliMCParticle secondD=static_cast<AliMCParticle>(GetMCTrackFromMCEvent(ilast));
8df8e382	396
8df8e382	397	//TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
b2a297fa	398	if (firstD->Charge()>0){
	399	if (firstD->PdgCode()!=-11) return kFALSE;
	400	if (secondD->PdgCode()!=11) return kFALSE;
	401	}else{
	402	if (firstD->PdgCode()!=11) return kFALSE;
	403	if (secondD->PdgCode()!=-11) return kFALSE;
	404	}
	405
	406	return kTRUE;
	407	}
	408
	409	//____________________________________________________________
	410	Bool_t AliDielectronMC::IsMCMotherToEEaod(const AliAODMCParticle *particle, Int_t pdgMother)
	411	{
	412	//
	413	// Check if the Mother 'particle' is of type pdgMother and decays to e+e-
	414	// AOD case
	415	//
	416	if (particle->GetPdgCode()!=pdgMother) return kFALSE;
	417	if (particle->GetNDaughters()!=2) return kFALSE;
	418
	419	Int_t ifirst = particle->GetDaughter(0);
	420	Int_t ilast = particle->GetDaughter(1);
	421
	422	//check number of daughters
	423	if ((ilast-ifirst)!=1) return kFALSE;
	424	AliAODMCParticle firstD=static_cast<AliAODMCParticle>(GetMCTrackFromMCEvent(ifirst));
	425	AliAODMCParticle secondD=static_cast<AliAODMCParticle>(GetMCTrackFromMCEvent(ilast));
	426
8df8e382	427	//TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
b2a297fa	428	if (firstD->Charge()>0){
	429	if (firstD->GetPdgCode()!=11) return kFALSE;
	430	if (secondD->GetPdgCode()!=-11) return kFALSE;
	431	}else{
	432	if (firstD->GetPdgCode()!=-11) return kFALSE;
	433	if (secondD->GetPdgCode()!=11) return kFALSE;
	434	}
	435	return kTRUE;
	436	}
	437
	438	//____________________________________________________________
a655b716	439	Int_t AliDielectronMC::GetLabelMotherWithPdg(const AliVParticle particle1, const AliVParticle particle2, Int_t pdgMother)
b2a297fa	440	{
	441	//
	442	// test if mother of particle 1 and 2 has pdgCode pdgMother and is the same;
	443	//
a655b716	444	if (!fMCEvent) return -1;
b2a297fa	445
a655b716	446	if (fAnaType==kESD) return GetLabelMotherWithPdgESD(particle1, particle2, pdgMother);
a655b716	447	else if (fAnaType==kAOD) return GetLabelMotherWithPdgAOD(particle1, particle2, pdgMother);
b2a297fa	448
a655b716	449	return -1;
b2a297fa	450	}
	451
	452	//____________________________________________________________
a655b716	453	Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle particle1, const AliVParticle particle2, Int_t pdgMother)
b2a297fa	454	{
b2a297fa	455	//
a655b716	456	// test if mother of particle 1 and 2 has pdgCode +-11 (electron),
a655b716	457	// have the same mother and the mother had pdg code pdgMother
b2a297fa	458	// ESD case
8df8e382	459	//TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
b2a297fa	460	//
	461
	462	AliMCParticle mcPart1=static_cast<AliMCParticle>(GetMCTrackFromMCEvent(particle1->GetLabel()));
	463	AliMCParticle mcPart2=static_cast<AliMCParticle>(GetMCTrackFromMCEvent(particle2->GetLabel()));
	464
a655b716	465	if (!mcPart1\|\|!mcPart2) return -1;
b2a297fa	466
	467	Int_t lblMother1=mcPart1->GetMother();
	468	Int_t lblMother2=mcPart2->GetMother();
	469
	470	AliMCParticle mcMother1=static_cast<AliMCParticle>(GetMCTrackFromMCEvent(lblMother1));
a655b716	471	if (!mcMother1) return -1;
	472	if (lblMother1!=lblMother2) return -1;
	473	if (TMath::Abs(mcPart1->PdgCode())!=11) return -1;
572b0139	474	if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1;
a655b716	475	if (mcMother1->PdgCode()!=pdgMother) return -1;
b2a297fa	476
a655b716	477	return lblMother1;
b2a297fa	478	}
	479
	480	//____________________________________________________________
a655b716	481	Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle particle1, const AliVParticle particle2, Int_t pdgMother)
b2a297fa	482	{
b2a297fa	483	//
a655b716	484	// test if mother of particle 1 and 2 has pdgCode +-11 (electron),
a655b716	485	// have the same mother and the mother had pdg code pdgMother
b2a297fa	486	// AOD case
8df8e382	487	//TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!!
b2a297fa	488	//
	489	AliAODMCParticle mcPart1=static_cast<AliAODMCParticle>(GetMCTrackFromMCEvent(particle1->GetLabel()));
	490	AliAODMCParticle mcPart2=static_cast<AliAODMCParticle>(GetMCTrackFromMCEvent(particle2->GetLabel()));
	491
a655b716	492	if (!mcPart1\|\|!mcPart2) return -1;
b2a297fa	493
	494	Int_t lblMother1=mcPart1->GetMother();
	495	Int_t lblMother2=mcPart2->GetMother();
	496
	497	AliAODMCParticle mcMother1=static_cast<AliAODMCParticle>(GetMCTrackFromMCEvent(lblMother1));
	498
a655b716	499	if (!mcMother1) return -1;
	500	if (lblMother1!=lblMother2) return -1;
	501	if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1;
572b0139	502	if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1;
a655b716	503	if (mcMother1->GetPdgCode()!=pdgMother) return -1;
b2a297fa	504
a655b716	505	return lblMother1;
b2a297fa	506	}
b2a297fa	507
6551594b	508	//____________________________________________________________
	509	void AliDielectronMC::GetDaughters(const TObject mother, AliVParticle &d1, AliVParticle* &d2)
	510	{
	511	//
	512	// Get First two daughters of the mother
	513	//
	514
	515	Int_t lblD1=-1;
	516	Int_t lblD2=-1;
	517	d1=0;
	518	d2=0;
	519	if (!fMCEvent) return;
	520	if (fAnaType==kAOD){
	521	const AliAODMCParticle aodMother=static_cast<const AliAODMCParticle>(mother);
	522	lblD1=aodMother->GetDaughter(0);
	523	lblD2=aodMother->GetDaughter(1);
	524	} else if (fAnaType==kESD){
	525	const AliMCParticle aodMother=static_cast<const AliMCParticle>(mother);
	526	lblD1=aodMother->GetFirstDaughter();
	527	lblD2=aodMother->GetLastDaughter();
	528	}
	529	d1=fMCEvent->GetTrack(lblD1);
	530	d2=fMCEvent->GetTrack(lblD2);
	531	}
8df8e382	532
	533	//____________________________________________________________
	534	Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair)
	535	{
	536	//
	537	// Check whether two particles have the same mother
	538	//
	539
	540	const AliVParticle * daughter1 = pair->GetFirstDaughter();
	541	const AliVParticle * daughter2 = pair->GetSecondDaughter();
	542
83e37742	543	AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
83e37742	544	AliVParticle *mcDaughter2=GetMCTrackFromMCEvent(daughter2->GetLabel());
8df8e382	545	if (!mcDaughter1 \|\| !mcDaughter2) return 0;
8df8e382	546
83e37742	547	Int_t labelMother1=-1;
83e37742	548	Int_t labelMother2=-1;
8df8e382	549
83e37742	550	if (mcDaughter1->IsA()==AliMCParticle::Class()){
	551	labelMother1=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
	552	labelMother2=(static_cast<AliMCParticle*>(mcDaughter2))->GetMother();
	553	} else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
	554	labelMother1=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
	555	labelMother2=(static_cast<AliAODMCParticle*>(mcDaughter2))->GetMother();
	556	}
8df8e382	557
83e37742	558	Bool_t sameMother=(labelMother1>-1)&&(labelMother2>-1)&&(labelMother1==labelMother2);
8df8e382	559
83e37742	560	return sameMother;
8df8e382	561	}
	562
	563
	564
	565
	566
	567
	568
	569