#include <AliAnalysisManager.h>
#include <AliAODHandler.h>
#include <AliESDInputHandler.h>
#include <AliAODInputHandler.h>
#include <AliMCEventHandler.h>
#include <AliMCEvent.h>
#include <AliMCParticle.h>
#include <AliAODMCParticle.h>
#include <AliAODMCHeader.h>
#include <AliStack.h>
#include <AliESDEvent.h>
#include <AliAODEvent.h>
#include <AliESDtrack.h>
#include <AliAODTrack.h>
#include <AliLog.h>
#include <TClonesArray.h>
#include <TParticle.h>
#include <TMCProcess.h>
#include "AliDielectronSignalMC.h"
#include "AliDielectronMC.h"
ClassImp(AliDielectronMC)
AliDielectronMC* AliDielectronMC::fgInstance=0x0;
AliDielectronMC* AliDielectronMC::Instance()
{
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()==AliAODInputHandler::Class()) type=kAOD;
AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
if(type == kESD) hasMC=mcHandler!=0x0;
}
fgInstance=new AliDielectronMC(type);
fgInstance->SetHasMC(hasMC);
return fgInstance;
}
AliDielectronMC::AliDielectronMC(AnalysisType type):
fMCEvent(0x0),
fStack(0x0),
fAnaType(type),
fHasMC(kTRUE),
fMcArray(0x0)
{
}
AliDielectronMC::~AliDielectronMC()
{
}
void AliDielectronMC::Initialize()
{
if (!ConnectMCEvent()) AliError("Initialization of MC object failed!");
}
Int_t AliDielectronMC::GetNMCTracks()
{
if(fAnaType == kESD){
if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
return fMCEvent->GetNumberOfTracks();}
else if(fAnaType == kAOD){
if(!fMcArray) { AliError("No fMcArray"); return 0; }
return fMcArray->GetEntriesFast();
}
return 0;
}
Int_t AliDielectronMC::GetNMCTracksFromStack()
{
if (!fStack){ AliError("No fStack"); return -999; }
return fStack->GetNtrack();
}
Int_t AliDielectronMC::GetNPrimary() const
{
if (!fMCEvent){ AliError("No fMCEvent"); return 0; }
return fMCEvent->GetNumberOfPrimaries();
}
Int_t AliDielectronMC::GetNPrimaryFromStack()
{
if (!fStack){ AliError("No fStack"); return -999; }
return fStack->GetNprimary();
}
AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t label) const
{
if (label<0) return NULL;
AliVParticle * track=0x0;
if(fAnaType == kESD){
if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
track = fMCEvent->GetTrack(label);
} else if(fAnaType == kAOD) {
if (!fMcArray){ AliError("No fMcArray"); return NULL;}
if (label>fMcArray->GetEntriesFast()) { AliDebug(10,Form("track %d out of array size %d",label,fMcArray->GetEntriesFast())); return NULL;}
track = (AliVParticle*)fMcArray->At(label);
}
return track;
}
Bool_t AliDielectronMC::ConnectMCEvent()
{
fMcArray = 0x0;
fMCEvent = 0x0;
if(fAnaType == kESD){
AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
if (!mcHandler){ return kFALSE; }
if (!mcHandler->InitOk() ) return kFALSE;
if (!mcHandler->TreeK() ) return kFALSE;
if (!mcHandler->TreeTR() ) return kFALSE;
AliMCEvent* mcEvent = mcHandler->MCEvent();
if (!mcEvent){ return kFALSE; }
fMCEvent = mcEvent;
if (!UpdateStack()) return kFALSE;
}
else if(fAnaType == kAOD)
{
AliAODInputHandler* aodHandler=(AliAODInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
if (!aodHandler) return kFALSE;
AliAODEvent *aod=aodHandler->GetEvent();
if (!aod) return kFALSE;
fMcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
if (!fMcArray){ return kFALSE; }
else fHasMC=kTRUE;
}
return kTRUE;
}
Bool_t AliDielectronMC::UpdateStack()
{
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)
{
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;
}
AliAODMCParticle* AliDielectronMC::GetMCTrack( const AliAODTrack* _track)
{
if(!fMcArray) { AliError("No fMcArray"); return NULL;}
Int_t nStack = fMcArray->GetEntriesFast();
Int_t label = TMath::Abs(_track->GetLabel());
if(label > nStack) return NULL;
AliAODMCParticle *mctrack = (AliAODMCParticle*)fMcArray->At(label);
return mctrack;
}
TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
{
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)
{
AliMCParticle* mcpart = GetMCTrack(_track);
if (!mcpart) return NULL;
if(mcpart->GetMother()<0) return NULL;
AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(mcpart->GetMother()));
if (!mcmother) return NULL;
return mcmother;
}
AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODTrack* _track)
{
AliAODMCParticle* mcpart = GetMCTrack(_track);
if (!mcpart) return NULL;
if(mcpart->GetMother() < 0) return NULL;
AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(mcpart->GetMother()));
if (!mcmother) return NULL;
return mcmother;
}
AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){
if(_particle->GetMother() < 0) return NULL;
AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
return mcmother;
}
AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
if( _particle->GetMother() < 0) return NULL;
AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMcArray->At(_particle->GetMother()));
return mcmother;
}
TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
{
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)
{
AliMCParticle* mcpart = GetMCTrack(_track);
if (!mcpart) return -999;
return mcpart->PdgCode();
}
Int_t AliDielectronMC::GetMCPID(const AliAODTrack* _track)
{
AliAODMCParticle* mcpart = GetMCTrack(_track);
if (!mcpart) return -999;
return mcpart->PdgCode();
}
Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
{
TParticle* mcpart = GetMCTrackFromStack(_track);
if (!mcpart) return -999;
return mcpart->GetPdgCode();
}
Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track)
{
AliMCParticle* mcmother = GetMCTrackMother(_track);
if (!mcmother) return -999;
return mcmother->PdgCode();
}
Int_t AliDielectronMC::GetMotherPDG( const AliAODTrack* _track)
{
AliAODMCParticle* mcmother = GetMCTrackMother(_track);
if (!mcmother) return -999;
return mcmother->PdgCode();
}
Int_t AliDielectronMC::GetMotherPDG( const AliMCParticle* _track)
{
AliMCParticle* mcmother = GetMCTrackMother(_track);
if (!mcmother) return -999;
return mcmother->PdgCode();
}
Int_t AliDielectronMC::GetMotherPDG( const AliAODMCParticle* _track)
{
AliAODMCParticle* mcmother = GetMCTrackMother(_track);
if (!mcmother) return -999;
return mcmother->PdgCode();
}
Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
{
TParticle* mcmother = GetMCTrackMotherFromStack(_track);
if (!mcmother) return -999;
return mcmother->GetPdgCode();
}
Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
{
AliMCParticle* mcpart = GetMCTrack(_track);
if (!mcpart) return -999;
return 0;
}
Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
{
TParticle* mcpart = GetMCTrackFromStack(_track);
if (!mcpart) return -999;
return mcpart->GetUniqueID();
}
Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
{
AliMCParticle *mcmother=GetMCTrackMother(track);
if(!mcmother||!mcmother->Particle()) return -999;
return mcmother->Particle()->GetNDaughters();
}
Int_t AliDielectronMC::NumberOfDaughters(const AliAODTrack* track)
{
AliAODMCParticle *mcmother=GetMCTrackMother(track);
if(!mcmother) return -999;
return NumberOfDaughters(mcmother);
}
Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle)
{
AliMCParticle *mcmother=GetMCTrackMother(particle);
if(!mcmother||!mcmother->Particle()) return -999;
return mcmother->Particle()->GetNDaughters();
}
Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle)
{
AliAODMCParticle *mcmother=GetMCTrackMother(particle);
if(!mcmother) return -999;
return mcmother->GetNDaughters();
}
Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track)
{
AliMCParticle* mcmother = GetMCTrackMother(_track);
if (!mcmother) return -999;
return 0;
}
Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
{
TParticle* mcmother = GetMCTrackMotherFromStack(_track);
if (!mcmother) return -999;
return mcmother->GetUniqueID();
}
Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother)
{
if (fAnaType==kESD && !fMCEvent) return kFALSE;
if (fAnaType==kAOD && !fMcArray) 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)
{
if (particle->PdgCode()!=pdgMother) return kFALSE;
Int_t ifirst = particle->GetFirstDaughter();
Int_t ilast = particle->GetLastDaughter();
if ((ilast-ifirst)!=1) return kFALSE;
AliMCParticle *firstD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ifirst));
AliMCParticle *secondD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ilast));
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)
{
if (particle->GetPdgCode()!=pdgMother) return kFALSE;
if (particle->GetNDaughters()!=2) return kFALSE;
Int_t ifirst = particle->GetDaughter(0);
Int_t ilast = particle->GetDaughter(1);
if ((ilast-ifirst)!=1) return kFALSE;
AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ilast));
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)
{
if (fAnaType==kESD){
if (!fMCEvent) return -1;
return GetLabelMotherWithPdgESD(particle1, particle2, pdgMother);
}
else if (fAnaType==kAOD)
{
if (!fMcArray) return -1;
return GetLabelMotherWithPdgAOD(particle1, particle2, pdgMother);
}
return -1;
}
Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother)
{
Int_t lblPart1 = TMath::Abs(particle1->GetLabel());
Int_t lblPart2 = TMath::Abs(particle2->GetLabel());
AliMCParticle *mcPart1=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblPart1));
AliMCParticle *mcPart2=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(lblPart2));
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)
{
Int_t lblPart1 = TMath::Abs(particle1->GetLabel());
Int_t lblPart2 = TMath::Abs(particle2->GetLabel());
AliAODMCParticle *mcPart1=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblPart1));
AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(lblPart2));
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)
{
Int_t lblD1=-1;
Int_t lblD2=-1;
d1=0;
d2=0;
if (fAnaType==kAOD){
if(!fMcArray) return;
const AliAODMCParticle *aodMother=static_cast<const AliAODMCParticle*>(mother);
lblD1=aodMother->GetDaughter(0);
lblD2=aodMother->GetDaughter(1);
d1 = (AliVParticle*)fMcArray->At(lblD1);
d2 = (AliVParticle*)fMcArray->At(lblD2);
} else if (fAnaType==kESD){
if (!fMCEvent) return;
const AliMCParticle *aodMother=static_cast<const AliMCParticle*>(mother);
lblD1=aodMother->GetFirstDaughter();
lblD2=aodMother->GetLastDaughter();
d1=fMCEvent->GetTrack(lblD1);
d2=fMCEvent->GetTrack(lblD2);
}
}
Int_t AliDielectronMC::GetMothersLabel(Int_t daughterLabel) const {
if(daughterLabel<0) return -1;
if (fAnaType==kAOD) {
if(!fMcArray) return -1;
if(GetMCTrackFromMCEvent(daughterLabel))
return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
} else if(fAnaType==kESD) {
if (!fMCEvent) return -1;
if(GetMCTrackFromMCEvent(daughterLabel))
return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
}
return -1;
}
Int_t AliDielectronMC::GetPdgFromLabel(Int_t label) const {
if(label<0) return 0;
if(fAnaType==kAOD) {
if(!fMcArray) return 0;
return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
} else if(fAnaType==kESD) {
if (!fMCEvent) return 0;
return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(label)))->PdgCode();
}
return 0;
}
Bool_t AliDielectronMC::ComparePDG(Int_t particlePDG, Int_t requiredPDG, Bool_t pdgExclusion, Bool_t checkBothCharges) const {
Bool_t result = kTRUE;
Int_t absRequiredPDG = TMath::Abs(requiredPDG);
switch(absRequiredPDG) {
case 0:
result = kTRUE;
break;
case 100:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=100 && TMath::Abs(particlePDG)<=199;
else {
if(requiredPDG>0) result = particlePDG>=100 && particlePDG<=199;
if(requiredPDG<0) result = particlePDG>=-199 && particlePDG<=-100;
}
break;
case 1000:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=1000 && TMath::Abs(particlePDG)<=1999;
else {
if(requiredPDG>0) result = particlePDG>=1000 && particlePDG<=1999;
if(requiredPDG<0) result = particlePDG>=-1999 && particlePDG<=-1000;
}
break;
case 200:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=200 && TMath::Abs(particlePDG)<=299;
else {
if(requiredPDG>0)result = particlePDG>=200 && particlePDG<=299;
if(requiredPDG<0)result = particlePDG>=-299 && particlePDG<=-200;
}
break;
case 2000:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=2000 && TMath::Abs(particlePDG)<=2999;
else {
if(requiredPDG>0) result = particlePDG>=2000 && particlePDG<=2999;
if(requiredPDG<0) result = particlePDG>=-2999 && particlePDG<=-2000;
}
break;
case 300:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=300 && TMath::Abs(particlePDG)<=399;
else {
if(requiredPDG>0) result = particlePDG>=300 && particlePDG<=399;
if(requiredPDG<0) result = particlePDG>=-399 && particlePDG<=-300;
}
break;
case 3000:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=3000 && TMath::Abs(particlePDG)<=3999;
else {
if(requiredPDG>0) result = particlePDG>=3000 && particlePDG<=3999;
if(requiredPDG<0) result = particlePDG>=-3999 && particlePDG<=-3000;
}
break;
case 400:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=499;
else {
if(requiredPDG>0) result = particlePDG>=400 && particlePDG<=499;
if(requiredPDG<0) result = particlePDG>=-499 && particlePDG<=-400;
}
break;
case 401:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439;
else {
if(requiredPDG>0) result = particlePDG>=400 && particlePDG<=439;
if(requiredPDG<0) result = particlePDG>=-439 && particlePDG<=-400;
}
break;
case 402:
if(checkBothCharges)
result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439) ||
(TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399);
else {
if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=439) ||
(particlePDG>=4000 && particlePDG<=4399);
if(requiredPDG<0) result = (particlePDG>=-439 && particlePDG<=-400) ||
(particlePDG>=-4399 && particlePDG<=-4000);
}
break;
case 403:
if(checkBothCharges)
result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=499) ||
(TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4999);
else {
if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=499) ||
(particlePDG>=4000 && particlePDG<=4999);
if(requiredPDG<0) result = (particlePDG>=-499 && particlePDG<=-400) ||
(particlePDG>=-4999 && particlePDG<=-4000);
}
break;
case 4000:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4999;
else {
if(requiredPDG>0) result = particlePDG>=4000 && particlePDG<=4999;
if(requiredPDG<0) result = particlePDG>=-4999 && particlePDG<=-4000;
}
break;
case 4001:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399;
else {
if(requiredPDG>0) result = particlePDG>=4000 && particlePDG<=4399;
if(requiredPDG<0) result = particlePDG>=-4399 && particlePDG<=-4000;
}
break;
case 500:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=599;
else {
if(requiredPDG>0) result = particlePDG>=500 && particlePDG<=599;
if(requiredPDG<0) result = particlePDG>=-599 && particlePDG<=-500;
}
break;
case 501:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549;
else {
if(requiredPDG>0) result = particlePDG>=500 && particlePDG<=549;
if(requiredPDG<0) result = particlePDG>=-549 && particlePDG<=-500;
}
break;
case 502:
if(checkBothCharges)
result = (TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549) ||
(TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499);
else {
if(requiredPDG>0) result = (particlePDG>=500 && particlePDG<=549) ||
(particlePDG>=5000 && particlePDG<=5499);
if(requiredPDG<0) result = (particlePDG>=-549 && particlePDG<=-500) ||
(particlePDG>=-5499 && particlePDG<=-5000);
}
break;
case 503:
if(checkBothCharges)
result = (TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=599) ||
(TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5999);
else {
if(requiredPDG>0) result = (particlePDG>=500 && particlePDG<=599) ||
(particlePDG>=5000 && particlePDG<=5999);
if(requiredPDG<0) result = (particlePDG>=-599 && particlePDG<=-500) ||
(particlePDG>=-5999 && particlePDG<=-5000);
}
break;
case 5000:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5999;
else {
if(requiredPDG>0) result = particlePDG>=5000 && particlePDG<=5999;
if(requiredPDG<0) result = particlePDG>=-5999 && particlePDG<=-5000;
}
break;
case 5001:
if(checkBothCharges)
result = TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499;
else {
if(requiredPDG>0) result = particlePDG>=5000 && particlePDG<=5499;
if(requiredPDG<0) result = particlePDG>=-5499 && particlePDG<=-5000;
}
break;
case 902:
if(checkBothCharges)
result = (TMath::Abs(particlePDG)>=400 && TMath::Abs(particlePDG)<=439) ||
(TMath::Abs(particlePDG)>=4000 && TMath::Abs(particlePDG)<=4399) ||
(TMath::Abs(particlePDG)>=500 && TMath::Abs(particlePDG)<=549) ||
(TMath::Abs(particlePDG)>=5000 && TMath::Abs(particlePDG)<=5499);
else {
if(requiredPDG>0) result = (particlePDG>=400 && particlePDG<=439) ||
(particlePDG>=4000 && particlePDG<=4399) ||
(particlePDG>=500 && particlePDG<=549) ||
(particlePDG>=5000 && particlePDG<=5499);
if(requiredPDG<0) result = (particlePDG>=-439 && particlePDG<=-400) ||
(particlePDG>=-4399 && particlePDG<=-4000) ||
(particlePDG>=-549 && particlePDG<=-500) ||
(particlePDG>=-5499 && particlePDG<=-5000);
}
break;
default:
if(checkBothCharges)
result = (absRequiredPDG==TMath::Abs(particlePDG));
else
result = (requiredPDG==particlePDG);
}
if(absRequiredPDG!=0 && pdgExclusion) result = !result;
return result;
}
Bool_t AliDielectronMC::IsPhysicalPrimary(Int_t label) const {
if(label<0) return kFALSE;
if(fAnaType==kAOD) {
if(!fMcArray) return kFALSE;
return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->IsPhysicalPrimary();
} else if(fAnaType==kESD) {
if (!fMCEvent) return kFALSE;
return fStack->IsPhysicalPrimary(label);
}
return kFALSE;
}
Bool_t AliDielectronMC::IsSecondaryFromWeakDecay(Int_t label) const {
if(label<0) return kFALSE;
if(fAnaType==kAOD) {
if(!fMcArray) return kFALSE;
return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->IsSecondaryFromWeakDecay();
} else if(fAnaType==kESD) {
if (!fMCEvent) return kFALSE;
return fStack->IsSecondaryFromWeakDecay(label);
}
return kFALSE;
}
Bool_t AliDielectronMC::IsSecondaryFromMaterial(Int_t label) const {
if(label<0) return kFALSE;
if(fAnaType==kAOD) {
if(!fMcArray) return kFALSE;
return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label)))->IsSecondaryFromMaterial();
} else if(fAnaType==kESD) {
if (!fMCEvent) return kFALSE;
return fStack->IsSecondaryFromMaterial(label);
}
return kFALSE;
}
Bool_t AliDielectronMC::CheckGEANTProcess(Int_t label, TMCProcess process) const {
if(label<0) return kFALSE;
if(fAnaType==kAOD) {
if(!fMcArray) return kFALSE;
UInt_t processID = static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label))->GetMCProcessCode();
return (process==processID);
} else if(fAnaType==kESD) {
if (!fMCEvent) return kFALSE;
AliError(Form("return of GEANT process not implemented for ESD "));
return kFALSE;
}
return kFALSE;
}
Bool_t AliDielectronMC::CheckParticleSource(Int_t label, AliDielectronSignalMC::ESource source) const {
switch (source) {
case AliDielectronSignalMC::kDontCare :
return kTRUE;
break;
case AliDielectronSignalMC::kPrimary :
return (label>=0 && label<=GetNPrimary());
break;
case AliDielectronSignalMC::kFinalState :
return IsPhysicalPrimary(label);
break;
case AliDielectronSignalMC::kDirect :
return (label>=0 && GetMothersLabel(label)<0);
break;
case AliDielectronSignalMC::kNoCocktail :
return (label>=0 && GetMothersLabel(label)>=0);
break;
case AliDielectronSignalMC::kSecondary :
return (label>=GetNPrimary() && !IsPhysicalPrimary(label));
break;
case AliDielectronSignalMC::kSecondaryFromWeakDecay :
return (IsSecondaryFromWeakDecay(label));
break;
case AliDielectronSignalMC::kSecondaryFromMaterial :
return (IsSecondaryFromMaterial(label));
break;
default :
return kFALSE;
}
return kFALSE;
}
Bool_t AliDielectronMC::CheckIsRadiative(Int_t label) const
{
if(label<0) return kFALSE;
if(fAnaType==kAOD) {
if(!fMcArray) return kFALSE;
AliAODMCParticle *mother=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label));
if (!mother) return kFALSE;
const Int_t nd=mother->GetNDaughters();
if (nd==2) return kFALSE;
for (Int_t i=2; i<nd; ++i)
if (GetMCTrackFromMCEvent(mother->GetDaughter(0)+i)->PdgCode()!=22) return kFALSE;
} else if(fAnaType==kESD) {
if (!fMCEvent) return kFALSE;
AliMCParticle *mother=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(label));
if (!mother) return kFALSE;
const Int_t nd=(mother->GetLastDaughter()-mother->GetFirstDaughter()+1);
if (nd==2) return kFALSE;
for (Int_t i=2; i<nd; ++i)
if (GetMCTrackFromMCEvent(mother->GetFirstDaughter()+i)->PdgCode()!=22) return kFALSE;
}
return kTRUE;
}
Bool_t AliDielectronMC::CheckRadiativeDecision(Int_t mLabel, const AliDielectronSignalMC * const signalMC) const
{
if (!signalMC) return kFALSE;
if (signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kAll) return kTRUE;
Bool_t isRadiative=CheckIsRadiative(mLabel);
if ((signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kIsRadiative) && !isRadiative) return kFALSE;
if ((signalMC->GetJpsiRadiative()==AliDielectronSignalMC::kIsNotRadiative) && isRadiative) return kFALSE;
return kTRUE;
}
Bool_t AliDielectronMC::IsMCTruth(Int_t label, AliDielectronSignalMC* signalMC, Int_t branch) const {
if(label<0) label *= -1;
AliVParticle* part = GetMCTrackFromMCEvent(label);
if (!part) {
AliError(Form("Could not find MC particle with label %d",label));
return kFALSE;
}
if(signalMC->GetCheckGEANTProcess() && !CheckGEANTProcess(label,signalMC->GetGEANTProcess())) return kFALSE;
if(!ComparePDG(part->PdgCode(),signalMC->GetLegPDG(branch),signalMC->GetLegPDGexclude(branch),signalMC->GetCheckBothChargesLegs(branch))) return kFALSE;
if(!CheckParticleSource(label, signalMC->GetLegSource(branch))) return kFALSE;
AliVParticle* mcMother=0x0;
Int_t mLabel = -1;
if(signalMC->GetMotherPDG(branch)!=0 || signalMC->GetMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
if(part) {
mLabel = GetMothersLabel(label);
mcMother = GetMCTrackFromMCEvent(mLabel);
}
if(!mcMother && !signalMC->GetMotherPDGexclude(branch)) return kFALSE;
if(!ComparePDG((mcMother ? mcMother->PdgCode() : 0),signalMC->GetMotherPDG(branch),signalMC->GetMotherPDGexclude(branch),signalMC->GetCheckBothChargesMothers(branch))) return kFALSE;
if(!CheckParticleSource(mLabel, signalMC->GetMotherSource(branch))) return kFALSE;
if (!CheckRadiativeDecision(mLabel, signalMC)) return kFALSE;
}
AliVParticle* mcGrandMother=0x0;
if(signalMC->GetGrandMotherPDG(branch)!=0 || signalMC->GetGrandMotherSource(branch)!=AliDielectronSignalMC::kDontCare) {
Int_t gmLabel = -1;
if(mcMother) {
gmLabel = GetMothersLabel(mLabel);
mcGrandMother = static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(gmLabel));
}
if(!mcGrandMother && !signalMC->GetGrandMotherPDGexclude(branch)) return kFALSE;
if(!ComparePDG((mcGrandMother ? mcGrandMother->PdgCode() : 0),signalMC->GetGrandMotherPDG(branch),signalMC->GetGrandMotherPDGexclude(branch),signalMC->GetCheckBothChargesGrandMothers(branch))) return kFALSE;
if(!CheckParticleSource(gmLabel, signalMC->GetGrandMotherSource(branch))) return kFALSE;
}
return kTRUE;
}
Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielectronSignalMC* signalMC) const {
const AliVParticle * mcD1 = pair->GetFirstDaughterP();
const AliVParticle * mcD2 = pair->GetSecondDaughterP();
Int_t labelD1 = (mcD1 ? TMath::Abs(mcD1->GetLabel()) : -1);
Int_t labelD2 = (mcD2 ? TMath::Abs(mcD2->GetLabel()) : -1);
Int_t d1Pdg = GetPdgFromLabel(labelD1);
Int_t d2Pdg = GetPdgFromLabel(labelD2);
AliVParticle* mcM1=0x0;
AliVParticle* mcM2=0x0;
AliVParticle* mcG1 = 0x0;
AliVParticle* mcG2 = 0x0;
Bool_t directTerm = kTRUE;
directTerm = directTerm && mcD1 && ComparePDG(d1Pdg,signalMC->GetLegPDG(1),signalMC->GetLegPDGexclude(1),signalMC->GetCheckBothChargesLegs(1))
&& CheckParticleSource(labelD1, signalMC->GetLegSource(1));
directTerm = directTerm && mcD2 && ComparePDG(d2Pdg,signalMC->GetLegPDG(2),signalMC->GetLegPDGexclude(2),signalMC->GetCheckBothChargesLegs(2))
&& CheckParticleSource(labelD2, signalMC->GetLegSource(2));
Int_t labelM1 = -1;
if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
labelM1 = GetMothersLabel(labelD1);
if(labelD1>-1 && labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
directTerm = directTerm && (mcM1 || signalMC->GetMotherPDGexclude(1))
&& ComparePDG((mcM1 ? mcM1->PdgCode() : 0),signalMC->GetMotherPDG(1),signalMC->GetMotherPDGexclude(1),signalMC->GetCheckBothChargesMothers(1))
&& CheckParticleSource(labelM1, signalMC->GetMotherSource(1))
&& CheckRadiativeDecision(labelM1,signalMC);
}
Int_t labelM2 = -1;
if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
labelM2 = GetMothersLabel(labelD2);
if(labelD2>-1 && labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
directTerm = directTerm && (mcM2 || signalMC->GetMotherPDGexclude(2))
&& ComparePDG((mcM2 ? mcM2->PdgCode() : 0),signalMC->GetMotherPDG(2),signalMC->GetMotherPDGexclude(2),signalMC->GetCheckBothChargesMothers(2))
&& CheckParticleSource(labelM2, signalMC->GetMotherSource(2))
&& CheckRadiativeDecision(labelM2,signalMC);
}
Int_t labelG1 = -1;
if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
labelG1 = GetMothersLabel(labelM1);
if(mcM1 && labelG1>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
directTerm = directTerm && (mcG1 || signalMC->GetGrandMotherPDGexclude(1))
&& ComparePDG((mcG1 ? mcG1->PdgCode() : 0),signalMC->GetGrandMotherPDG(1),signalMC->GetGrandMotherPDGexclude(1),signalMC->GetCheckBothChargesGrandMothers(1))
&& CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(1));
}
Int_t labelG2 = -1;
if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
labelG2 = GetMothersLabel(labelM2);
if(mcM2 && labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
directTerm = directTerm && (mcG2 || signalMC->GetGrandMotherPDGexclude(2))
&& ComparePDG((mcG2 ? mcG2->PdgCode() : 0),signalMC->GetGrandMotherPDG(2),signalMC->GetGrandMotherPDGexclude(2),signalMC->GetCheckBothChargesGrandMothers(2))
&& CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(2));
}
Bool_t crossTerm = kTRUE;
crossTerm = crossTerm && mcD2
&& ComparePDG(d2Pdg,signalMC->GetLegPDG(1),signalMC->GetLegPDGexclude(1),signalMC->GetCheckBothChargesLegs(1))
&& CheckParticleSource(labelD2, signalMC->GetLegSource(1));
crossTerm = crossTerm && mcD1
&& ComparePDG(d1Pdg,signalMC->GetLegPDG(2),signalMC->GetLegPDGexclude(2),signalMC->GetCheckBothChargesLegs(2))
&& CheckParticleSource(labelD1, signalMC->GetLegSource(2));
if(signalMC->GetMotherPDG(1)!=0 || signalMC->GetMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
if(!mcM2 && labelD2>-1) {
labelM2 = GetMothersLabel(labelD2);
if(labelM2>-1) mcM2 = GetMCTrackFromMCEvent(labelM2);
}
crossTerm = crossTerm && (mcM2 || signalMC->GetMotherPDGexclude(1))
&& ComparePDG((mcM2 ? mcM2->PdgCode() : 0),signalMC->GetMotherPDG(1),signalMC->GetMotherPDGexclude(1),signalMC->GetCheckBothChargesMothers(1))
&& CheckParticleSource(labelM2, signalMC->GetMotherSource(1))
&& CheckRadiativeDecision(labelM2,signalMC);
}
if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
if(!mcM1 && labelD1>-1) {
labelM1 = GetMothersLabel(labelD1);
if(labelM1>-1) mcM1 = GetMCTrackFromMCEvent(labelM1);
}
crossTerm = crossTerm && (mcM1 || signalMC->GetMotherPDGexclude(2))
&& ComparePDG((mcM1 ? mcM1->PdgCode() : 0),signalMC->GetMotherPDG(2),signalMC->GetMotherPDGexclude(2),signalMC->GetCheckBothChargesMothers(2))
&& CheckParticleSource(labelM1, signalMC->GetMotherSource(2))
&& CheckRadiativeDecision(labelM1,signalMC);
}
if(signalMC->GetGrandMotherPDG(1)!=0 || signalMC->GetGrandMotherSource(1)!=AliDielectronSignalMC::kDontCare) {
if(!mcG2 && mcM2) {
labelG2 = GetMothersLabel(labelM2);
if(labelG2>-1) mcG2 = GetMCTrackFromMCEvent(labelG2);
}
crossTerm = crossTerm && (mcG2 || signalMC->GetGrandMotherPDGexclude(1))
&& ComparePDG((mcG2 ? mcG2->PdgCode() : 0),signalMC->GetGrandMotherPDG(1),signalMC->GetGrandMotherPDGexclude(1),signalMC->GetCheckBothChargesGrandMothers(1))
&& CheckParticleSource(labelG2, signalMC->GetGrandMotherSource(1));
}
if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
if(!mcG1 && mcM1) {
labelG1 = GetMothersLabel(labelM1);
if(labelG1>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
}
crossTerm = crossTerm && (mcG1 || signalMC->GetGrandMotherPDGexclude(2))
&& ComparePDG((mcG1 ? mcG1->PdgCode() : 0),signalMC->GetGrandMotherPDG(2),signalMC->GetGrandMotherPDGexclude(2),signalMC->GetCheckBothChargesGrandMothers(2))
&& CheckParticleSource(labelG1, signalMC->GetGrandMotherSource(2));
}
Bool_t motherRelation = kTRUE;
if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kSame) {
motherRelation = motherRelation && HaveSameMother(pair);
}
if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent) {
motherRelation = motherRelation && !HaveSameMother(pair);
}
Bool_t processGEANT = kTRUE;
if(signalMC->GetCheckGEANTProcess()) {
if(!CheckGEANTProcess(labelD1,signalMC->GetGEANTProcess()) &&
!CheckGEANTProcess(labelD2,signalMC->GetGEANTProcess()) ) processGEANT= kFALSE;
}
return ((directTerm || crossTerm) && motherRelation && processGEANT);
}
Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair) const
{
const AliVParticle * daughter1 = pair->GetFirstDaughterP();
const AliVParticle * daughter2 = pair->GetSecondDaughterP();
if (!daughter1 || !daughter2) return 0;
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;
}
Int_t AliDielectronMC::IsJpsiPrimary(const AliDielectronPair * pair)
{
if(!HaveSameMother(pair)) return 2;
AliVParticle *mcDaughter1=GetMCTrackFromMCEvent((pair->GetFirstDaughterP())->GetLabel());
Int_t labelMother=-1;
if (mcDaughter1->IsA()==AliMCParticle::Class()){
labelMother=(static_cast<AliMCParticle*>(mcDaughter1))->GetMother();
} else if (mcDaughter1->IsA()==AliAODMCParticle::Class()) {
labelMother=(static_cast<AliAODMCParticle*>(mcDaughter1))->GetMother();
}
AliVParticle* mcMother=GetMCTrackFromMCEvent(labelMother);
if(!IsMCMotherToEE(mcMother,443)) return 2;
return IsJpsiPrimary(mcMother);
}
Int_t AliDielectronMC::IsJpsiPrimary(const AliVParticle * particle)
{
Int_t labelMoth=-1;
Int_t pdgCode;
if (particle->IsA()==AliMCParticle::Class()){
labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
while(labelMoth>0){
particle = GetMCTrackFromMCEvent(labelMoth);
pdgCode = TMath::Abs((static_cast<const AliMCParticle*>(particle))->PdgCode());
if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
labelMoth = (static_cast<const AliMCParticle*>(particle))->GetMother();
}
}
else if (particle->IsA()==AliAODMCParticle::Class()){
labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
while(labelMoth>0){
particle = GetMCTrackFromMCEvent(labelMoth);
pdgCode = TMath::Abs((static_cast<const AliAODMCParticle*>(particle))->PdgCode());
if((pdgCode>500 && pdgCode<600) || (pdgCode>5000 && pdgCode<6000)) return 1;
labelMoth = (static_cast<const AliAODMCParticle*>(particle))->GetMother();
}
}
return 0;
}
Bool_t AliDielectronMC::GetPrimaryVertex(Double_t &primVtxX, Double_t &primVtxY, Double_t &primVtxZ){
if(fAnaType == kESD){
const AliVVertex* mcVtx = fMCEvent->GetPrimaryVertex();
if(!mcVtx) return kFALSE;
primVtxX = mcVtx->GetX();
primVtxY = mcVtx->GetY();
primVtxZ = mcVtx->GetZ();
}else if(fAnaType == kAOD){
AliAODEvent *aod=((AliAODInputHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()))->GetEvent();
if(!aod) return kFALSE;
AliAODMCHeader *mcHead = dynamic_cast<AliAODMCHeader*>(aod->FindListObject(AliAODMCHeader::StdBranchName()));
if(!mcHead) return kFALSE;
primVtxX = mcHead->GetVtxX();
primVtxY = mcHead->GetVtxY();
primVtxZ = mcHead->GetVtxZ();
}
return kTRUE;
}