[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;
  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());

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