]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGDQ/dielectron/AliDielectronMC.cxx
add proptections (marcel), add GEANT process for MC
[u/mrichter/AliRoot.git] / PWGDQ / dielectron / AliDielectronMC.cxx
index eb56285c92d535e780b12035187fa5e76cf0af16..f34857bb1af66168c117b1c9a9705336b443e16c 100644 (file)
@@ -30,6 +30,7 @@
 #include <AliMCEvent.h>
 #include <AliMCParticle.h>
 #include <AliAODMCParticle.h>
+#include <AliAODMCHeader.h>
 #include <AliStack.h>
 #include <AliESDEvent.h>
 #include <AliAODEvent.h>
 
 #include <TClonesArray.h>
 #include <TParticle.h>
+#include <TMCProcess.h>
 
 #include "AliDielectronSignalMC.h"
 #include "AliDielectronMC.h"
 
+
+ClassImp(AliDielectronMC)
+
 AliDielectronMC* AliDielectronMC::fgInstance=0x0;
 
 //____________________________________________________________
@@ -149,20 +154,21 @@ Int_t AliDielectronMC::GetNPrimaryFromStack()
 }
 
 //____________________________________________________________
-AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t itrk) const
+AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t label) const
 {
   //
   // return MC track directly from MC event
+  // used not only for tracks but for mothers as well, therefore do not use abs(label)
   //
-  if (itrk<0) return NULL;
+  if (label<0) return NULL;
   AliVParticle * track=0x0;
   if(fAnaType == kESD){
     if (!fMCEvent){ AliError("No fMCEvent"); return NULL;}
-    track = fMCEvent->GetTrack(itrk); //  tracks from MC event (ESD)
+    track = fMCEvent->GetTrack(label); //  tracks from MC event (ESD)
   } else if(fAnaType == kAOD) {
     if (!fMcArray){ AliError("No fMcArray"); return NULL;}
-    if (itrk>fMcArray->GetEntriesFast()) { AliWarning(Form("track %d out of array size %d",itrk,fMcArray->GetEntriesFast())); return NULL;}
-    track = (AliVParticle*)fMcArray->At(itrk); //  tracks from MC event (AOD)
+    if (label>fMcArray->GetEntriesFast()) { AliDebug(10,Form("track %d out of array size %d",label,fMcArray->GetEntriesFast())); return NULL;}
+    track = (AliVParticle*)fMcArray->At(label); //  tracks from MC event (AOD)
   }
   return track;
 }
@@ -173,8 +179,11 @@ Bool_t AliDielectronMC::ConnectMCEvent()
   //
   // connect stack object from the mc handler
   //
+
+  fMcArray = 0x0;
+  fMCEvent = 0x0;
+  
   if(fAnaType == kESD){
-    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;
@@ -189,8 +198,11 @@ Bool_t AliDielectronMC::ConnectMCEvent()
   }
   else if(fAnaType == kAOD)
   {
-    fMcArray = 0x0;
-    AliAODEvent *aod=((AliAODInputHandler*)((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()))->GetEvent();
+    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){ /*AliError("Could not retrieve MC array!");*/ return kFALSE; }
     else fHasMC=kTRUE;
@@ -218,11 +230,9 @@ 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());
+  Int_t label  = TMath::Abs(_track->GetLabel()); // negative label indicate poor matching quality
   if(label>nStack)return NULL;
-
   AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(label));
   return mctrack;
 }
@@ -235,8 +245,9 @@ AliAODMCParticle* AliDielectronMC::GetMCTrack( const AliAODTrack* _track)
   // return MC track
   //
  if(!fMcArray) { AliError("No fMcArray"); return NULL;}
- Int_t label = _track->GetLabel();
- if(label < 0 || label > fMcArray->GetEntriesFast()) return NULL;
+ Int_t nStack = fMcArray->GetEntriesFast();
+ Int_t label  = TMath::Abs(_track->GetLabel()); // negative label indicate poor matching quality
+ if(label > nStack) return NULL;
  AliAODMCParticle *mctrack = (AliAODMCParticle*)fMcArray->At(label);
  return mctrack; 
 }
@@ -262,6 +273,7 @@ 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;
@@ -285,6 +297,7 @@ AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle)
   //
   // return MC track mother
   //
+  if(_particle->GetMother() < 0) return NULL;
   AliMCParticle* mcmother = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
   return mcmother;
 }
@@ -367,6 +380,28 @@ Int_t AliDielectronMC::GetMotherPDG( const AliAODTrack* _track)
   return mcmother->PdgCode();
 }
 
+//________________________________________________________
+Int_t AliDielectronMC::GetMotherPDG( const AliMCParticle* _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::GetMotherPDG( const AliAODMCParticle* _track)
+{
+  //
+  // return PDG code of the mother track from the MC truth info
+  //
+  AliAODMCParticle* mcmother = GetMCTrackMother(_track);
+  if (!mcmother) return -999;
+  return mcmother->PdgCode();
+}
+
 //____________________________________________________________
 Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
 {
@@ -408,7 +443,7 @@ Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track)
   //
   AliMCParticle *mcmother=GetMCTrackMother(track);
   if(!mcmother||!mcmother->Particle()) return -999;
-//   return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
+  //   return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0;
   return mcmother->Particle()->GetNDaughters();
 }
 
@@ -579,16 +614,18 @@ Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, c
   // 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()));
+  // negative label indicate poor matching quality
+  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;
@@ -607,14 +644,16 @@ Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, c
   // 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()));
+  // negative label indicate poor matching quality
+  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;
@@ -658,14 +697,17 @@ void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, Ali
 Int_t AliDielectronMC::GetMothersLabel(Int_t daughterLabel) const {
   //
   //  Get the label of the mother for particle with label daughterLabel
+  //  NOTE: for tracks, the absolute label should be passed
   //
   if(daughterLabel<0) return -1;
   if (fAnaType==kAOD) {
     if(!fMcArray) return -1;
-    return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
+    if(GetMCTrackFromMCEvent(daughterLabel))
+      return (static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
   } else if(fAnaType==kESD) {
     if (!fMCEvent) return -1;
-    return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
+    if(GetMCTrackFromMCEvent(daughterLabel))
+      return (static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(daughterLabel)))->GetMother();
   }
   return -1;
 }
@@ -675,6 +717,7 @@ Int_t AliDielectronMC::GetMothersLabel(Int_t daughterLabel) const {
 Int_t AliDielectronMC::GetPdgFromLabel(Int_t label) const {
   //
   //  Get particle code using the label from stack
+  //  NOTE: for tracks, the absolute label should be passed
   //
   if(label<0) return 0;
   if(fAnaType==kAOD) {
@@ -856,6 +899,23 @@ Bool_t AliDielectronMC::ComparePDG(Int_t particlePDG, Int_t requiredPDG, Bool_t
       if(requiredPDG<0) result = particlePDG>=-5499 && particlePDG<=-5000;
     }
     break;
+  case 902:      // // open charm,beauty  mesons and baryons together
+    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:          // all specific cases
     if(checkBothCharges)
       result = (absRequiredPDG==TMath::Abs(particlePDG));
@@ -867,7 +927,6 @@ Bool_t AliDielectronMC::ComparePDG(Int_t particlePDG, Int_t requiredPDG, Bool_t
   return result;
 }
 
-
 //________________________________________________________________________________
 Bool_t AliDielectronMC::IsPhysicalPrimary(Int_t label) const {
   //
@@ -892,11 +951,67 @@ Bool_t AliDielectronMC::IsPhysicalPrimary(Int_t label) const {
   return kFALSE;
 }
 
+//________________________________________________________________________________
+Bool_t AliDielectronMC::IsSecondaryFromWeakDecay(Int_t label) const {
+  //
+  // Check if the particle with label "label" is a physical secondary from weak decay according to the
+  // definition in AliStack::IsSecondaryFromWeakDecay(Int_t label)
+  //
+  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 {
+  //
+  // Check if the particle with label "label" is a physical secondary from weak decay according to the
+  // definition in AliStack::IsSecondaryFromMaterial(Int_t label)
+  //
+  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 {
+  //
+  //  Check the GEANT process for the particle 
+  //  NOTE: for tracks the absolute label should be passed
+  //
+  if(label<0) return kFALSE;
+  if(fAnaType==kAOD) {
+    if(!fMcArray) return kFALSE;
+    UInt_t processID = static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(label))->GetMCProcessCode();
+    //    printf("process: id %d --> %s \n",processID,TMCProcessName[processID]);
+    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 {
   //
   //  Check the source for the particle 
+  //  NOTE: for tracks the absolute label should be passed
   //
 
   switch (source) {
@@ -930,20 +1045,79 @@ Bool_t AliDielectronMC::CheckParticleSource(Int_t label, AliDielectronSignalMC::
       // 3.) Certain particles added via MC generator cocktails (e.g. J/psi added to pythia MB events)
       return (label>=0 && GetMothersLabel(label)<0);
       break;
+    case AliDielectronSignalMC::kNoCocktail :
+      // Particles from the HIJING event and NOT from the AliGenCocktail
+      return (label>=0 && GetMothersLabel(label)>=0);
+      break;
     case AliDielectronSignalMC::kSecondary :          
       // particles which are created by the interaction of final state primaries with the detector
       // or particles from strange weakly decaying particles (e.g. lambda, kaons, etc.)
       return (label>=GetNPrimary() && !IsPhysicalPrimary(label));
     break;
+    case AliDielectronSignalMC::kSecondaryFromWeakDecay :          
+      // secondary particle from weak decay
+      // or particles from strange weakly decaying particles (e.g. lambda, kaons, etc.)
+      return (IsSecondaryFromWeakDecay(label));
+    break;
+    case AliDielectronSignalMC::kSecondaryFromMaterial :
+      // secondary particle from material
+      return (IsSecondaryFromMaterial(label));
+    break;
     default :
       return kFALSE;
   }
   return kFALSE;
 }
 
+//________________________________________________________________________________
+Bool_t AliDielectronMC::CheckIsRadiative(Int_t label) const
+{
+  //
+  // Check if the particle has a three body decay, one being a photon
+  //
+  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; //last daughter is photon
+  } 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; //last daughters are photons
+  }
+  return kTRUE;
+}
+
+//________________________________________________________________________________
+Bool_t AliDielectronMC::CheckRadiativeDecision(Int_t mLabel, const AliDielectronSignalMC * const signalMC) const
+{
+  //
+  // Check for the decision of the radiative type request
+  //
+  
+  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) {
+Bool_t AliDielectronMC::IsMCTruth(Int_t label, AliDielectronSignalMC* signalMC, Int_t branch) const {
   //
   // Check if the particle corresponds to the MC truth in signalMC in the branch specified
   //
@@ -958,7 +1132,10 @@ Bool_t AliDielectronMC::IsMCTruth(Int_t label, AliDielectronSignalMC* signalMC,
     AliError(Form("Could not find MC particle with label %d",label));
     return kFALSE;
   }
-  
+
+  // check geant process if set
+  if(signalMC->GetCheckGEANTProcess() && !CheckGEANTProcess(label,signalMC->GetGEANTProcess())) return kFALSE;
+
   // check the leg
   if(!ComparePDG(part->PdgCode(),signalMC->GetLegPDG(branch),signalMC->GetLegPDGexclude(branch),signalMC->GetCheckBothChargesLegs(branch))) return kFALSE;
   if(!CheckParticleSource(label, signalMC->GetLegSource(branch))) return kFALSE;
@@ -972,9 +1149,12 @@ Bool_t AliDielectronMC::IsMCTruth(Int_t label, AliDielectronSignalMC* signalMC,
       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;
+
+    //check for radiative deday
+    if (!CheckRadiativeDecision(mLabel, signalMC)) return kFALSE;
   }
   
   // check the grandmother
@@ -1003,16 +1183,12 @@ Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielec
   //
  
   // legs (daughters)
-  const AliVParticle * mcD1 = pair->GetFirstDaughter();
-  const AliVParticle * mcD2 = pair->GetSecondDaughter();
-  Int_t labelD1 = (mcD1 ? mcD1->GetLabel() : -1);
-  Int_t labelD2 = (mcD2 ? mcD2->GetLabel() : -1);
-  if(labelD1<0) labelD1 *= -1;
-  if(labelD2<0) labelD2 *= -1;
-  Int_t d1Pdg = 0;
-  Int_t d2Pdg = 0;
-  d1Pdg=GetPdgFromLabel(labelD1);
-  d2Pdg=GetPdgFromLabel(labelD2);
+  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);
   
   // mothers
   AliVParticle* mcM1=0x0;
@@ -1036,9 +1212,10 @@ Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielec
   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)) 
+    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));
+                 && CheckParticleSource(labelM1, signalMC->GetMotherSource(1))
+                 && CheckRadiativeDecision(labelM1,signalMC);
   }
   
   Int_t labelM2 = -1;
@@ -1047,7 +1224,8 @@ Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielec
     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));
+                 && CheckParticleSource(labelM2, signalMC->GetMotherSource(2))
+                 && CheckRadiativeDecision(labelM2,signalMC);
   }
  
   // grand-mothers
@@ -1088,7 +1266,8 @@ Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielec
     }
     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));
+                && CheckParticleSource(labelM2, signalMC->GetMotherSource(1))
+                && CheckRadiativeDecision(labelM2,signalMC);
   }
   
   if(signalMC->GetMotherPDG(2)!=0 || signalMC->GetMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
@@ -1098,7 +1277,8 @@ Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielec
     }
     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));
+                && CheckParticleSource(labelM1, signalMC->GetMotherSource(2))
+                && CheckRadiativeDecision(labelM1,signalMC);
   }
 
   // grand-mothers
@@ -1115,7 +1295,7 @@ Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielec
   if(signalMC->GetGrandMotherPDG(2)!=0 || signalMC->GetGrandMotherSource(2)!=AliDielectronSignalMC::kDontCare) {
     if(!mcG1 && mcM1) {
       labelG1 = GetMothersLabel(labelM1);
-      if(labelG2>-1) mcG1 = GetMCTrackFromMCEvent(labelG1);
+      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))
@@ -1129,8 +1309,15 @@ Bool_t AliDielectronMC::IsMCTruth(const AliDielectronPair* pair, const AliDielec
   if(signalMC->GetMothersRelation()==AliDielectronSignalMC::kDifferent) {
     motherRelation = motherRelation && !HaveSameMother(pair);
   }
-  return ((directTerm || crossTerm) && motherRelation);
+
+  // check geant process if set
+  Bool_t processGEANT = kTRUE;
+  if(signalMC->GetCheckGEANTProcess()) {
+    if(!CheckGEANTProcess(labelD1,signalMC->GetGEANTProcess()) &&
+       !CheckGEANTProcess(labelD2,signalMC->GetGEANTProcess())   ) processGEANT= kFALSE;
+  }
+
+  return ((directTerm || crossTerm) && motherRelation && processGEANT);
 }
 
 
@@ -1142,8 +1329,8 @@ Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair) const
   // Check whether two particles have the same mother
   //
 
-  const AliVParticle * daughter1 = pair->GetFirstDaughter();
-  const AliVParticle * daughter2 = pair->GetSecondDaughter();
+  const AliVParticle * daughter1 = pair->GetFirstDaughterP();
+  const AliVParticle * daughter2 = pair->GetSecondDaughterP();
   if (!daughter1 || !daughter2) return 0;
 
   AliVParticle *mcDaughter1=GetMCTrackFromMCEvent(daughter1->GetLabel());
@@ -1173,7 +1360,7 @@ Int_t AliDielectronMC::IsJpsiPrimary(const AliDielectronPair * pair)
  //         "1" for secondary jpsi (from beauty)
  //         "2" for background  
  if(!HaveSameMother(pair)) return 2;
- AliVParticle *mcDaughter1=GetMCTrackFromMCEvent((pair->GetFirstDaughter())->GetLabel());
+ AliVParticle *mcDaughter1=GetMCTrackFromMCEvent((pair->GetFirstDaughterP())->GetLabel());
  Int_t labelMother=-1;
 
   if (mcDaughter1->IsA()==AliMCParticle::Class()){
@@ -1215,3 +1402,25 @@ Int_t AliDielectronMC::IsJpsiPrimary(const AliVParticle * particle)
   }
   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;
+}