]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG/CaloTrackCorrBase/AliMCAnalysisUtils.cxx
Converting PWGCaloTrackCorrBase to native cmake
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / AliMCAnalysisUtils.cxx
index daf8cb6232e18ff4e5b32587942bfbf36f4e1d5a..0f45b63757d2a176d5fedfe989946deea069cbd9 100755 (executable)
@@ -37,6 +37,7 @@
 #include "AliStack.h"
 #include "AliGenPythiaEventHeader.h"
 #include "AliAODMCParticle.h"
+#include "AliLog.h"
 
 ClassImp(AliMCAnalysisUtils)
 
@@ -44,9 +45,12 @@ ClassImp(AliMCAnalysisUtils)
 AliMCAnalysisUtils::AliMCAnalysisUtils() : 
 TObject(), 
 fCurrentEvent(-1), 
-fDebug(-1), 
+fDebug(0),
 fJetsList(new TList), 
-fMCGenerator("PYTHIA")
+fMCGenerator(kPythia),
+fMCGeneratorString("PYTHIA"),
+fDaughMom(),  fDaughMom2(),
+fMotherMom(), fGMotherMom()
 {
   //Ctor
 }
@@ -56,19 +60,21 @@ AliMCAnalysisUtils::~AliMCAnalysisUtils()
 {
   // Remove all pointers.
   
-  if (fJetsList) {
+  if (fJetsList)
+  {
     fJetsList->Clear();
     delete fJetsList ;
   }     
 }
 
 //_____________________________________________________________________________________________
-Int_t AliMCAnalysisUtils::CheckCommonAncestor(const Int_t index1, const Int_t index2, 
+Int_t AliMCAnalysisUtils::CheckCommonAncestor(Int_t index1, Int_t index2, 
                                               const AliCaloTrackReader* reader, 
                                               Int_t & ancPDG, Int_t & ancStatus, 
                                               TLorentzVector & momentum, TVector3 & prodVertex) 
 {
-  //Check the first common ancestor of 2 clusters, given the most likely labels of the primaries generating such clusters.
+  // Check the first common ancestor of 2 clusters, given the most likely labels of the primaries generating such clusters.
+  
   Int_t label1[100];
   Int_t label2[100];
   label1[0]= index1;
@@ -91,46 +97,60 @@ Int_t AliMCAnalysisUtils::CheckCommonAncestor(const Int_t index1, const Int_t in
       Int_t label=label1[0];
       if(label < 0) return -1;
       
-      while(label > -1 && counter1 < 99){
+      while(label > -1 && counter1 < 99)
+      {
         counter1++;
         AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
-        if(mom){
+        if(mom)
+        {
           label  = mom->GetMother() ;
           label1[counter1]=label;
         }
         //printf("\t counter %d, label %d\n", counter1,label);
       }
+     
       //printf("Org label2=%d,\n",label2[0]);
       label=label2[0];
       if(label < 0) return -1;
-      while(label > -1 && counter2 < 99){
+      
+      while(label > -1 && counter2 < 99)
+      {
         counter2++;
         AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
-        if(mom){
+        if(mom)
+        {
           label  = mom->GetMother() ;
           label2[counter2]=label;
         }
         //printf("\t counter %d, label %d\n", counter2,label);
       }
     }//AOD MC
-    else { //Kine stack from ESDs 
+    else
+    { //Kine stack from ESDs
       AliStack * stack = reader->GetStack();
+      
       Int_t label=label1[0];
-      while(label > -1 && counter1 < 99){
+      while(label > -1 && counter1 < 99)
+      {
         counter1++;
         TParticle * mom = stack->Particle(label);
-        if(mom){
+        if(mom)
+        {
           label  = mom->GetFirstMother() ;
           label1[counter1]=label;
         }
         //printf("\t counter %d, label %d\n", counter1,label);
       }
+      
       //printf("Org label2=%d,\n",label2[0]);
+      
       label=label2[0];
-      while(label > -1 && counter2 < 99){
+      while(label > -1 && counter2 < 99)
+      {
         counter2++;
         TParticle * mom = stack->Particle(label);
-        if(mom){
+        if(mom)
+        {
           label  = mom->GetFirstMother() ;
           label2[counter2]=label;
         }
@@ -139,34 +159,47 @@ Int_t AliMCAnalysisUtils::CheckCommonAncestor(const Int_t index1, const Int_t in
     }// Kine stack from ESDs
   }//First labels not the same
   
-  if((counter1==99 || counter2==99) && fDebug >=0) printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2);
+//  if((counter1==99 || counter2==99) && fDebug >=0)
+//    printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2);
   //printf("CheckAncestor:\n");
+  
   Int_t commonparents = 0;
   Int_t ancLabel = -1;
   //printf("counters %d %d \n",counter1, counter2);
-  for (Int_t c1 = 0; c1 < counter1; c1++) {
-    for (Int_t c2 = 0; c2 < counter2; c2++) {
-      if(label1[c1]==label2[c2] && label1[c1]>-1) {
+  for (Int_t c1 = 0; c1 < counter1; c1++)
+  {
+    for (Int_t c2 = 0; c2 < counter2; c2++)
+    {
+      if(label1[c1]==label2[c2] && label1[c1]>-1)
+      {
         ancLabel = label1[c1];
         commonparents++;
-        if(reader->ReadAODMCParticles()){
+        
+        if(reader->ReadAODMCParticles())
+        {
           AliAODMCParticle * mom = (AliAODMCParticle *) reader->GetAODMCParticles()->At(label1[c1]);
-          if (mom) {
+          
+          if (mom)
+          {
             ancPDG    = mom->GetPdgCode();
             ancStatus = mom->GetStatus();
             momentum.SetPxPyPzE(mom->Px(),mom->Py(),mom->Pz(),mom->E());
             prodVertex.SetXYZ(mom->Xv(),mom->Yv(),mom->Zv());
           }
         }
-        else {
+        else
+        {
           TParticle * mom = (reader->GetStack())->Particle(label1[c1]);
-          if (mom) {
+          
+          if (mom)
+          {
             ancPDG    = mom->GetPdgCode();
             ancStatus = mom->GetStatusCode();
             mom->Momentum(momentum);
             prodVertex.SetXYZ(mom->Vx(),mom->Vy(),mom->Vz());
           }
         }
+        
         //First ancestor found, end the loops
         counter1=0;
         counter2=0;
@@ -174,62 +207,79 @@ Int_t AliMCAnalysisUtils::CheckCommonAncestor(const Int_t index1, const Int_t in
     }//second cluster loop
   }//first cluster loop
   
+  if(ancLabel < 0)
+  {
+    ancPDG    = -10000;
+    ancStatus = -10000;
+    momentum.SetXYZT(0,0,0,0);
+    prodVertex.SetXYZ(-10,-10,-10);
+  }
+  
   return ancLabel;
 }
 
-//_____________________________________________________________________
-Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, 
-                                      const Int_t nlabels, 
-                                      const AliCaloTrackReader* reader)
+//________________________________________________________________________________________
+Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, Int_t nlabels,
+                                      const AliCaloTrackReader* reader, Int_t calorimeter)
 {
-  //Play with the montecarlo particles if available
+  // Play with the montecarlo particles if available.
+  
   Int_t tag = 0;
   
-  if(nlabels<=0) {
-    printf("AliMCAnalysisUtils::CheckOrigin(nlabel<=0) - No MC labels available, please check!!!\n");
+  if( nlabels <= 0 )
+  {
+    AliWarning("No MC labels available, please check!!!");
     return kMCBadLabel;
   }
   
+  TObjArray* arrayCluster = 0;
+  if      ( calorimeter == AliCaloTrackReader::kEMCAL ) arrayCluster = reader->GetEMCALClusters();
+  else if ( calorimeter == AliCaloTrackReader::kPHOS  ) arrayCluster = reader->GetPHOSClusters ();
+  
   //Select where the information is, ESD-galice stack or AOD mcparticles branch
   if(reader->ReadStack()){
-    tag = CheckOriginInStack(label, nlabels, reader->GetStack());
+    tag = CheckOriginInStack(label, nlabels, reader->GetStack(), arrayCluster);
   }
   else if(reader->ReadAODMCParticles()){
-    tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles());
+    tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles(),arrayCluster);
   }
   
   return tag ;
 }
 
-//_____________________________________________________________________
-Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t label, 
-                                      const AliCaloTrackReader* reader)
+//____________________________________________________________________________________________________
+Int_t AliMCAnalysisUtils::CheckOrigin(Int_t label, const AliCaloTrackReader* reader, Int_t calorimeter)
 {
-  //Play with the montecarlo particles if available
+  // Play with the montecarlo particles if available.
+  
   Int_t tag = 0;
   
-  if(label<0) {
-    printf("AliMCAnalysisUtils::CheckOrigin(label<0) - No MC labels available, please check!!!\n");
+  if( label < 0 )
+  {
+    AliWarning("No MC labels available, please check!!!");
     return kMCBadLabel;
   }
   
+  TObjArray* arrayCluster = 0;
+  if      ( calorimeter == AliCaloTrackReader::kEMCAL ) arrayCluster = reader->GetEMCALClusters();
+  else if ( calorimeter == AliCaloTrackReader::kPHOS  ) arrayCluster = reader->GetPHOSClusters();
+  
   Int_t labels[]={label};
   
   //Select where the information is, ESD-galice stack or AOD mcparticles branch
   if(reader->ReadStack()){
-    tag = CheckOriginInStack(labels, 1,reader->GetStack());
+    tag = CheckOriginInStack(labels, 1,reader->GetStack(),arrayCluster);
   }
   else if(reader->ReadAODMCParticles()){
-    tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles());
+    tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles(),arrayCluster);
   }
   
   return tag ;
 }      
 
-//_________________________________________________________________
-Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, 
-                                             const Int_t nlabels, 
-                                             AliStack* stack) 
+//__________________________________________________________________________________________
+Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, Int_t nlabels,
+                                             AliStack* stack, const TObjArray* arrayCluster)
 {
   // Play with the MC stack if available. Tag particles depending on their origin.
   // Do same things as in CheckOriginInAOD but different input.
@@ -239,16 +289,17 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
   //about its heritage, but one can also use it directly with stack 
   //particles not connected to reconstructed entities
   
-  if(!stack) {
-    if (fDebug >=0) 
-      printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
+  if(!stack)
+  {
+    AliDebug(1,"Stack is not available, check analysis settings in configuration file, STOP!!");
     return -1;
   }
   
   Int_t tag = 0;
   Int_t label=labels[0];//Most significant particle contributing to the cluster
   
-  if(label >= 0 && label < stack->GetNtrack()){
+  if(label >= 0 && label < stack->GetNtrack())
+  {
     //MC particle of interest is the "mom" of the entity
     TParticle * mom = stack->Particle(label);
     Int_t iMom     = label;
@@ -256,32 +307,37 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
     Int_t mPdg     = TMath::Abs(mPdgSign);
     Int_t mStatus  = mom->GetStatusCode() ;
     Int_t iParent  = mom->GetFirstMother() ;
-    if(fDebug > 0 && label < 8 && fMCGenerator!="") printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
+    
+    //if( label < 8 && fMCGenerator != kBoxLike ) AliDebug(1,Form("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d",iParent));
     
     //GrandParent of the entity
     TParticle * parent = NULL;
     Int_t pPdg = -1;
     Int_t pStatus =-1;
-    if(iParent >= 0){
+    if(iParent >= 0)
+    {
       parent = stack->Particle(iParent);
-      if(parent){
+      
+      if(parent)
+      {
         pPdg = TMath::Abs(parent->GetPdgCode());
         pStatus = parent->GetStatusCode();  
       }
     }
-    else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent);
+    else AliDebug(1,Form("Parent with label %d",iParent));
+    
+    AliDebug(2,"Cluster most contributing mother and its parent:");
+    AliDebug(2,Form("\t Mother label %d, pdg %d, status %d",iMom, mPdg, mStatus));
+    AliDebug(2,Form("\t Parent label %d, pdg %d, status %d",iParent, pPdg, pStatus));
     
-    if(fDebug > 2 ) {
-      printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n");
-      printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
-      printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
-    }
-         
     //Check if "mother" of entity is converted, if not, get the first non converted mother
-    if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0){
+    if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0)
+    {
       SetTagBit(tag,kMCConversion);
+      
       //Check if the mother is photon or electron with status not stable
-      while ((pPdg == 22 || pPdg == 11) && mStatus != 1) {
+      while ((pPdg == 22 || pPdg == 11) && mStatus != 1)
+      {
         //Mother
         iMom     = mom->GetFirstMother();
         mom      = stack->Particle(iMom);
@@ -289,12 +345,15 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
         mPdg     = TMath::Abs(mPdgSign);
         mStatus  = mom->GetStatusCode() ;
         iParent  = mom->GetFirstMother() ;
-        if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent);
+       
+        //if(label < 8 ) AliDebug(1,Form("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent));
         
         //GrandParent
-        if(iParent >= 0){
+        if(iParent >= 0)
+        {
           parent = stack->Particle(iParent);
-          if(parent){
+          if(parent)
+          {
             pPdg = TMath::Abs(parent->GetPdgCode());
             pStatus = parent->GetStatusCode();  
           }
@@ -304,28 +363,28 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
           pStatus = 0;
           break;
         }
-      }//while   
-      if(fDebug > 2 ) {
-        printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n");
-        printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
-        printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus);
-      }
+      }//while
+      
+      AliDebug(2,"Converted photon/electron:");
+      AliDebug(2,Form("\t Mother label %d, pdg %d, status %d",iMom, mPdg, mStatus));
+      AliDebug(2,Form("\t Parent label %d, pdg %d, status %d",iParent, pPdg, pStatus));
       
     }//mother and parent are electron or photon and have status 0
-    else if((mPdg == 22 || mPdg == 11) && mStatus == 0){       
+    else if((mPdg == 22 || mPdg == 11) && mStatus == 0)
+    {
       //Still a conversion but only one electron/photon generated. Just from hadrons but not decays.
       if(pPdg == 2112 ||  pPdg == 211  ||  pPdg == 321 ||
-         pPdg == 2212 ||  pPdg == 130  ||  pPdg == 13 ) {
+         pPdg == 2212 ||  pPdg == 130  ||  pPdg == 13 )
+      {
         SetTagBit(tag,kMCConversion);
         iMom     = mom->GetFirstMother();
         mom      = stack->Particle(iMom);
         mPdgSign = mom->GetPdgCode();
         mPdg     = TMath::Abs(mPdgSign);
         
-        if(fDebug > 2 ) {
-          printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n");
-          printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus);
-        }
+        AliDebug(2,"Converted hadron:");
+        AliDebug(2,Form("\t Mother label %d, pdg %d, status %d",iMom, mPdg, mStatus));
+        
       }//hadron converted
       
       //Comment for the next lines, we do not check the parent of the hadron for the moment.
@@ -350,62 +409,77 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
     else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
     
     //check for pi0 and eta (shouldn't happen unless their decays were turned off)
-    else if(mPdg == 111)  {
+    else if(mPdg == 111)
+    {
+      
       SetTagBit(tag,kMCPi0Decay);
-      if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n");
+      
+      AliDebug(2,"First mother is directly pi0, not decayed by generator");
+      
       CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
+      
     }
-    else if(mPdg == 221) {
+    else if(mPdg == 221)
+    {
       SetTagBit(tag,kMCEtaDecay);
-      if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n");
+      
+      AliDebug(2,"First mother is directly eta, not decayed by generator");
+      
       CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster
     }
     //Photons  
-    else if(mPdg == 22){
+    else if(mPdg == 22)
+    {
       SetTagBit(tag,kMCPhoton);
-      if(mStatus == 1){ //undecayed particle
-        if(fMCGenerator == "PYTHIA"){
-          if(iParent < 8 && iParent > 5) {//outgoing partons
+      
+      if(pPdg == 111)
+      {
+        SetTagBit(tag,kMCPi0Decay);
+        
+        AliDebug(2,"PYTHIA pi0 decay photon,  parent pi0 with status 11");
+        
+        CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
+        // In case it did not merge, check if the decay companion is lost
+        if(!CheckTagBit(tag, kMCPi0) && !CheckTagBit(tag,kMCDecayPairInCalo))
+          CheckLostDecayPair(arrayCluster,iMom, iParent, stack, tag);
+
+        //printf("Bit set is Merged %d, Pair in calo %d, Lost %d\n",CheckTagBit(tag, kMCPi0),CheckTagBit(tag,kMCDecayPairInCalo),CheckTagBit(tag,kMCDecayPairLost));
+      }
+      else if (pPdg == 221)
+      {
+        SetTagBit(tag, kMCEtaDecay);
+        
+        AliDebug(2,"PYTHIA eta decay photon,  parent pi0 with status 11");
+        
+        CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
+        // In case it did not merge, check if the decay companion is lost
+        if(!CheckTagBit(tag, kMCEta) && !CheckTagBit(tag,kMCDecayPairInCalo))
+          CheckLostDecayPair(arrayCluster,iMom, iParent, stack, tag);
+      }
+      else if(mStatus == 1)
+      { //undecayed particle
+        if(fMCGenerator == kPythia)
+        {
+          if(iParent < 8 && iParent > 5)
+          {//outgoing partons
             if(pPdg == 22) SetTagBit(tag,kMCPrompt);
-            else SetTagBit(tag,kMCFragmentation);
+            else           SetTagBit(tag,kMCFragmentation);
           }//Outgoing partons 
-          else  if(iParent <= 5) {
+          else  if(iParent <= 5)
+          {
             SetTagBit(tag, kMCISR); //Initial state radiation
           }
-          else if(pStatus == 11){//Decay
-            if(pPdg == 111) {
-              SetTagBit(tag,kMCPi0Decay);
-              if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon,  parent pi0 with status 11 \n");
-              CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
-            }
-            else if (pPdg == 221) {
-              SetTagBit(tag, kMCEtaDecay);
-              if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon,  parent pi0 with status 11 \n");
-              CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
-            }
-            else SetTagBit(tag,kMCOtherDecay);
-          }//Decay
-          else {
-            if(fDebug > 1 && parent) printf("AliMCAnalysisUtils::CheckOrigingInStack() - what is it in PYTHIA? Wrong generator setting? Mother mPdg %d, status %d \n    Parent  iParent %d, pPdg %d %s, status %d\n",
-                                            mPdg, mStatus,iParent, pPdg, parent->GetName(),pStatus);
-            if(pPdg == 111) {
-              SetTagBit(tag,kMCPi0Decay);
-              if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon,  parent pi0 with status 11 \n");
-              CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
-            }
-            else if (pPdg == 221) {
-              SetTagBit(tag, kMCEtaDecay);
-              if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon,  parent pi0 with status 11 \n");
-              CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster
-            }
-            else SetTagBit(tag,kMCOtherDecay);
-          }
-        }//PYTHIA
-        
-        else if(fMCGenerator == "HERWIG"){       
-          if(pStatus < 197){//Not decay
-            while(1){
-              if(parent){
+          else  SetTagBit(tag,kMCUnknown);
+         }//PYTHIA
+      
+        else if(fMCGenerator == kHerwig)
+        {
+          if(pStatus < 197)
+          {//Not decay
+            while(1)
+            {
+              if(parent)
+              {
                 if(parent->GetFirstMother()<=5) break;
                 iParent = parent->GetFirstMother();
                 parent=stack->Particle(iParent);
@@ -414,93 +488,84 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
               } else break;
             }//Look for the parton
             
-            if(iParent < 8 && iParent > 5) {
+            if(iParent < 8 && iParent > 5)
+            {
               if(pPdg == 22) SetTagBit(tag,kMCPrompt);
-              else SetTagBit(tag,kMCFragmentation);
+              else           SetTagBit(tag,kMCFragmentation);
             }
             else SetTagBit(tag,kMCISR);//Initial state radiation
           }//Not decay
-          else{//Decay
-            if(pPdg == 111) {
-              SetTagBit(tag,kMCPi0Decay);
-              if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG pi0 decay photon \n");
-              CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
-            }
-            else if (pPdg == 221) {
-              SetTagBit(tag,kMCEtaDecay);
-              if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - HERWIG eta decay photon \n");
-              CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
-            }
-            else SetTagBit(tag,kMCOtherDecay);
-          }//Decay
+          else  SetTagBit(tag,kMCUnknown);
         }//HERWIG
-        
-        else SetTagBit(tag,kMCUnknown);
-        
-      }//Status 1 : created by event generator
-      
-      else if(mStatus == 0){ // geant
-        if(pPdg == 111) {
-          SetTagBit(tag,kMCPi0Decay);
-          if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC pi0 decay photon \n");
-          CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster
-        }
-        else if (pPdg == 221) {
-          SetTagBit(tag,kMCEtaDecay);
-          if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Transport MC eta decay photon \n");
-          CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCEta if 2 gammas in same cluster
-        }
-        else  SetTagBit(tag,kMCOtherDecay);    
-      }//status 0 : geant generated
-      
+      }
+      else  SetTagBit(tag,kMCOtherDecay);
+                  
     }//Mother Photon
     
     //Electron check.  Where did that electron come from?
     else if(mPdg == 11){ //electron
-      if(pPdg == 11 && parent){
+      if(pPdg == 11 && parent)
+      {
         Int_t iGrandma = parent->GetFirstMother();
-        if(iGrandma >= 0) {
+        if(iGrandma >= 0)
+        {
           TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother
           Int_t gPdg = TMath::Abs(gma->GetPdgCode());
           
-          if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
+          if      (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
           else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
         }
       }
-      SetTagBit(tag,kMCElectron);      
-      if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
-      if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
-      else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
+      
+      SetTagBit(tag,kMCElectron);
+      
+      AliDebug(1,"Checking ancestors of electrons");
+     
+      if      (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag, kMCDecayDalitz) ; } //Pi0 Dalitz decay
+      else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag, kMCDecayDalitz) ; } //Eta Dalitz decay
       else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay
-      else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay
-        if(parent){
+      else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000))
+      { //check charm decay
+        if(parent)
+        {
           Int_t iGrandma = parent->GetFirstMother();
-          if(iGrandma >= 0) {
+          if(iGrandma >= 0)
+          {
             TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm
             Int_t gPdg = TMath::Abs(gma->GetPdgCode());
             if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e
-            else SetTagBit(tag,kMCEFromC); //c-->e 
-          } else SetTagBit(tag,kMCEFromC); //c-->e 
+            else SetTagBit(tag,kMCEFromC); //c-->e
+          }
+          else SetTagBit(tag,kMCEFromC); //c-->e
         }//parent
-      } else {
+      }
+      else
+      {
         //if it is not from any of the above, where is it from?
         if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
+        
         else SetTagBit(tag,kMCOtherDecay);
-        if(fDebug > 0 && parent) printf("AliMCAnalysisUtils::CheckOriginInStack() - Status %d Electron from other origin: %s (pPdg = %d) %s (mpdg = %d)\n",mStatus,parent->GetName(),pPdg,mom->GetName(),mPdg);
+        
+        //if(parent) AliDebug(1,Form("Status %d Electron from other origin: %s (pPdg = %d) %s (mpdg = %d)",mStatus,parent->GetName(),pPdg,mom->GetName(),mPdg));
       }
     }//electron check
     //Cluster was made by something else
-    else {
-      if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n",mom->GetName(),mPdg,pPdg);
+    else
+    {
+      AliDebug(2,Form("\t Setting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)",
+                      mom->GetName(),mPdg,pPdg));
+      
       SetTagBit(tag,kMCUnknown);
     }
   }//Good label value
-  else{// Bad label 
-         
-    if(label < 0 && (fDebug >= 0)) 
-      printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***:  label %d \n", label);
-    if(label >=  stack->GetNtrack() &&  (fDebug >= 0)) 
-      printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
+  else
+  {// Bad label
+    if(label < 0)
+      AliWarning(Form("*** bad label or no stack ***:  label %d ", label));
+    
+    if(label >=  stack->GetNtrack())
+      AliWarning(Form("*** large label ***:  label %d, n tracks %d", label, stack->GetNtrack()));
+    
     SetTagBit(tag,kMCUnknown);
   }//Bad label
   
@@ -509,16 +574,16 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
 }
 
 
-//_________________________________________________________________________
-Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, 
-                                           const Int_t nlabels, 
-                                           const TClonesArray *mcparticles) 
+//________________________________________________________________________________________________________
+Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, Int_t nlabels,
+                                           const TClonesArray *mcparticles, const TObjArray* arrayCluster)
 {
   // Play with the MCParticles in AOD if available. Tag particles depending on their origin.
   // Do same things as in CheckOriginInStack but different input.
-  if(!mcparticles) {
-    if(fDebug >= 0)
-      printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
+  
+  if(!mcparticles)
+  {
+    AliDebug(1,"AODMCParticles is not available, check analysis settings in configuration file!!");
     return -1;
   }
        
@@ -526,46 +591,50 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
   Int_t label=labels[0];//Most significant particle contributing to the cluster
   
   Int_t nprimaries = mcparticles->GetEntriesFast();
-  if(label >= 0 && label < nprimaries){
+  if(label >= 0 && label < nprimaries)
+  {
     //Mother
     AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label);
     Int_t iMom     = label;
     Int_t mPdgSign = mom->GetPdgCode();
     Int_t mPdg     = TMath::Abs(mPdgSign);
     Int_t iParent  = mom->GetMother() ;
-    if(fDebug > 0 && label < 8 && fMCGenerator!="") printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
+    
+    //if(label < 8 && fMCGenerator != kBoxLike) AliDebug(1,Form("Mother is parton %d\n",iParent));
     
     //GrandParent
     AliAODMCParticle * parent = NULL ;
     Int_t pPdg = -1;
-    if(iParent >= 0){
+    if(iParent >= 0)
+    {
       parent = (AliAODMCParticle *) mcparticles->At(iParent);
       pPdg = TMath::Abs(parent->GetPdgCode());
     }
-    else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent);
+    else AliDebug(1,Form("Parent with label %d",iParent));
+    
+    AliDebug(2,"Cluster most contributing mother and its parent:");
+    AliDebug(2,Form("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()));
+    AliDebug(2,Form("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d",iParent, pPdg, parent?parent->IsPrimary():-1, parent?parent->IsPhysicalPrimary():-1));
     
-    if(fDebug > 2 ) {
-      printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n");
-      printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
-      if(parent)
-        printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
-    }
-         
     //Check if mother is converted, if not, get the first non converted mother
-    if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()){
+    if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary())
+    {
       SetTagBit(tag,kMCConversion);
+      
       //Check if the mother is photon or electron with status not stable
-      while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) {
+      while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary())
+      {
         //Mother
         iMom     = mom->GetMother();
         mom      = (AliAODMCParticle *) mcparticles->At(iMom);
         mPdgSign = mom->GetPdgCode();
         mPdg     = TMath::Abs(mPdgSign);
         iParent  = mom->GetMother() ;
-        if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent);
+        //if(label < 8 ) AliDebug(1, Form("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent));
         
         //GrandParent
-        if(iParent >= 0 && parent){
+        if(iParent >= 0 && parent)
+        {
           parent = (AliAODMCParticle *) mcparticles->At(iParent);
           pPdg = TMath::Abs(parent->GetPdgCode());
         }
@@ -574,28 +643,26 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
         
       }//while 
       
-      if(fDebug > 2 ) {
-        printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n");
-        printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
-        if(parent)
-          printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
-      }
+      AliDebug(2,"AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron:");
+      AliDebug(2,Form("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()));
+      AliDebug(2,Form("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d",iParent, pPdg, parent?parent->IsPrimary():-1, parent?parent->IsPhysicalPrimary():-1));
       
     }//mother and parent are electron or photon and have status 0 and parent is photon or electron
-    else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary()){  
+    else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary())
+    {
       //Still a conversion but only one electron/photon generated. Just from hadrons
       if(pPdg == 2112 ||  pPdg == 211 ||  pPdg == 321 ||  
-         pPdg == 2212 ||  pPdg == 130 ||  pPdg == 13 ) {
+         pPdg == 2212 ||  pPdg == 130 ||  pPdg == 13 )
+      {
         SetTagBit(tag,kMCConversion);
         iMom     = mom->GetMother();
         mom      = (AliAODMCParticle *) mcparticles->At(iMom);
         mPdgSign = mom->GetPdgCode();
         mPdg     = TMath::Abs(mPdgSign);
         
-        if(fDebug > 2 ) {
-          printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n");
-          printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary());
-        }
+       AliDebug(2,"AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron:");
+       AliDebug(2,Form("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()));
+      
       }//hadron converted
       
       //Comment for next lines, we do not check the parent of the hadron for the moment.
@@ -623,115 +690,133 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
     else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron);
     
     //check for pi0 and eta (shouldn't happen unless their decays were turned off)
-    else if(mPdg == 111)  {
+    else if(mPdg == 111)
+    {
       SetTagBit(tag,kMCPi0Decay);
-      if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n");
+      
+      AliDebug(2,"First mother is directly pi0, not decayed by generator");
+      
       CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
     }
-    else if(mPdg == 221)  {
+    else if(mPdg == 221)
+    {
       SetTagBit(tag,kMCEtaDecay);   
-      if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n");
+      
+      AliDebug(2,"First mother is directly eta, not decayed by generator");
+      
       CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
     }
     //Photons  
-    else if(mPdg == 22){
+    else if(mPdg == 22)
+    {
       SetTagBit(tag,kMCPhoton);
-      if(mom->IsPhysicalPrimary() && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG")) //undecayed particle
+      
+      if(pPdg == 111)
+      {
+        SetTagBit(tag,kMCPi0Decay);
+        
+        AliDebug(2,"Generator pi0 decay photon");
+        
+        CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
+        // In case it did not merge, check if the decay companion is lost
+        if(!CheckTagBit(tag, kMCPi0) && !CheckTagBit(tag,kMCDecayPairInCalo) && !CheckTagBit(tag,kMCDecayPairLost))
+        {
+          CheckLostDecayPair(arrayCluster,iMom, iParent, mcparticles, tag);
+        }
+        
+        //printf("Bit set is Merged %d, Pair in calo %d, Lost %d\n",CheckTagBit(tag, kMCPi0),CheckTagBit(tag,kMCDecayPairInCalo),CheckTagBit(tag,kMCDecayPairLost));
+      }
+      else if (pPdg == 221)
+      {
+        SetTagBit(tag, kMCEtaDecay);
+        
+        AliDebug(2,"Generator eta decay photon");
+        
+        CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
+        // In case it did not merge, check if the decay companion is lost
+        if(!CheckTagBit(tag, kMCEta) && !CheckTagBit(tag,kMCDecayPairInCalo))
+          CheckLostDecayPair(arrayCluster,iMom, iParent, mcparticles, tag);
+      }
+      else if( mom->IsPhysicalPrimary() && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) ) //undecayed particle
       {
-        if(iParent < 8 && iParent > 5 ) {//outgoing partons
+        if(iParent < 8 && iParent > 5 )
+        {//outgoing partons
           if(pPdg == 22) SetTagBit(tag,kMCPrompt);
           else SetTagBit(tag,kMCFragmentation);
         }//Outgoing partons
-        else if(iParent <= 5 && (fMCGenerator=="PYTHIA" || fMCGenerator=="HERWIG")) {
+        else if( iParent <= 5 && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) )
+        {
           SetTagBit(tag, kMCISR); //Initial state radiation
         }
-        else if(parent && parent->IsPrimary() && !parent->IsPhysicalPrimary()){//Decay
-          if(pPdg == 111){
-            SetTagBit(tag,kMCPi0Decay);
-            if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
-            CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster
-          }
-          else if (pPdg == 221) {
-            SetTagBit(tag, kMCEtaDecay);
-            if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
-            CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
-          }
-          else SetTagBit(tag,kMCOtherDecay);
-        }//Decay
-        else {
-          if(parent)printf("AliMCAnalysisUtils::CheckOriginInAOD() - what is it? Mother mPdg %d, is primary? %d, is physical %d \n    Parent  iParent %d, pPdg %d, is primary? %d, is physical? %d\n",
-                           mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(),iParent, pPdg,parent->IsPrimary(), parent->IsPhysicalPrimary());
-          SetTagBit(tag,kMCOtherDecay);//Check
-        }
+        else SetTagBit(tag,kMCUnknown);
       }//Physical primary
-      else if(!mom->IsPrimary()){      //Decays  
-        if(pPdg == 111){ 
-          SetTagBit(tag,kMCPi0Decay); 
-          if(fDebug > 2 ) 
-            printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC pi0 decay photon \n");
-          CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster          
-        }
-        else if (pPdg == 221) {
-          SetTagBit(tag,kMCEtaDecay);
-          if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Transport MC eta decay photon \n");
-          CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster
-        }
-        else  SetTagBit(tag,kMCOtherDecay);
-      }//not primary : geant generated, decays
-      else  {
-        //printf("UNKNOWN 1, mom  pdg %d, primary %d, physical primary %d; parent %d, pdg %d, primary %d, physical primary %d \n",
-        //mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary(), iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary());
-        SetTagBit(tag,kMCUnknown);
-      }
+      else SetTagBit(tag,kMCOtherDecay);
+
     }//Mother Photon
     
     //Electron check.  Where did that electron come from?
-    else if(mPdg == 11){ //electron
-      if(pPdg == 11 && parent){
+    else if(mPdg == 11)
+    { //electron
+      if(pPdg == 11 && parent)
+      {
         Int_t iGrandma = parent->GetMother();
-        if(iGrandma >= 0) {
+        if(iGrandma >= 0)
+        {
           AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma);
           Int_t gPdg = TMath::Abs(gma->GetPdgCode());
           
-          if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
+          if      (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson
           else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson
         }
       }
-      SetTagBit(tag,kMCElectron);      
-      if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons");
-      if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
-      else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
+      
+      SetTagBit(tag,kMCElectron);
+      
+      AliDebug(1,"Checking ancestors of electrons");
+      
+      if      (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag,kMCDecayDalitz);} //Pi0 Dalitz decay
+      else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag,kMCDecayDalitz);} //Eta Dalitz decay
       else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay
-      else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //c-hadron decay check
-        if(parent){
+      else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000))
+      { //c-hadron decay check
+        if(parent)
+        {
           Int_t iGrandma = parent->GetMother();
-          if(iGrandma >= 0) {
+          if(iGrandma >= 0)
+          {
             AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother
             Int_t gPdg = TMath::Abs(gma->GetPdgCode());
             if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay
             else SetTagBit(tag,kMCEFromC); //c-hadron decay
-          } else SetTagBit(tag,kMCEFromC); //c-hadron decay
+          }
+          else SetTagBit(tag,kMCEFromC); //c-hadron decay
         }//parent
-      } else { //prompt or other decay
+      } else
+      { //prompt or other decay
         TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg);
         TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg);
-        if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n",foo->GetName(), pPdg,foo1->GetName(),mPdg);
+        
+        AliDebug(1,Form("Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)",foo->GetName(), pPdg,foo1->GetName(),mPdg));
+        
         if(pPdg > 10000) SetTagBit(tag,kMCUnknown);
-        else SetTagBit(tag,kMCOtherDecay);
+        else             SetTagBit(tag,kMCOtherDecay);
       }      
     }//electron check
     //cluster was made by something else
-    else {
-      if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n",mPdg,pPdg);
+    else
+    {
+      AliDebug(1,Form("\t Setting kMCUnknown for cluster with pdg = %d, Parent pdg = %d",mPdg,pPdg));
       SetTagBit(tag,kMCUnknown);
     }
   }//Good label value
-  else{//Bad label
-         
-    if(label < 0 && (fDebug >= 0) ) 
-      printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***:  label %d \n", label);
-    if(label >=  mcparticles->GetEntriesFast() &&  (fDebug >= 0) ) 
-      printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***:  label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
+  else
+  {//Bad label
+    if(label < 0)
+      AliWarning(Form("*** bad label or no mcparticles ***:  label %d", label));
+    
+    if(label >=  mcparticles->GetEntriesFast())
+      AliWarning(Form("*** large label ***:  label %d, n tracks %d", label, mcparticles->GetEntriesFast()));
+  
     SetTagBit(tag,kMCUnknown);
     
   }//Bad label
@@ -740,34 +825,33 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
   
 }
 
-//_________________________________________________________________________
-void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, 
-                                                    const Int_t nlabels, 
-                                                    const Int_t mesonIndex, 
-                                                    AliStack *stack, 
+//_________________________________________________________________________________________
+void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,    Int_t nlabels,
+                                                    Int_t mesonIndex, AliStack *stack,
                                                     Int_t &tag)
 {
-  //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack
+  // Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack.
   
-  if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) {
-    if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
-                          labels[0],stack->GetNtrack(), nlabels);
+  if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1)
+  {
+    AliDebug(2,Form("Exit : label[0] %d, n primaries %d, nlabels %d",labels[0],stack->GetNtrack(), nlabels));
     return;
   }
   
   TParticle * meson = stack->Particle(mesonIndex);
   Int_t mesonPdg    = meson->GetPdgCode();
-  if(mesonPdg!=111 && mesonPdg!=221){ 
-    printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Wrong pi0/eta PDG : %d \n",mesonPdg);
+  if(mesonPdg!=111 && mesonPdg!=221)
+  {
+    AliWarning(Form("Wrong pi0/eta PDG : %d",mesonPdg));
     return;
   }
   
-  if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s, label %d\n",meson->GetName(), mesonIndex);
+  AliDebug(2,Form("%s, label %d",meson->GetName(), mesonIndex));
   
   //Check if meson decayed into 2 daughters or if both were kept.
-  if(meson->GetNDaughters() != 2){
-    if(fDebug > 2) 
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
+  if(meson->GetNDaughters() != 2)
+  {
+    AliDebug(2,Form("Not overalapped. Number of daughters is %d, not 2",meson->GetNDaughters()));
     return;
   }
   
@@ -778,92 +862,103 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
   TParticle *photon1 = stack->Particle(iPhoton1);
   
   //Check if both daughters are photons
-  if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22){
-    if(fDebug > 2) 
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
+  if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22)
+  {
+    AliDebug(2,Form("Not overalapped. PDG:  daughter 1 = %d, of daughter 2 = %d",photon0->GetPdgCode(),photon1->GetPdgCode()));
     return;
   }
   
-  if(fDebug > 2) 
-    printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
+  AliDebug(2,Form("Daughter labels : photon0 = %d, photon1 = %d",iPhoton0,iPhoton1));
   
   //Check if both photons contribute to the cluster
   Bool_t okPhoton0 = kFALSE;
   Bool_t okPhoton1 = kFALSE;
   
-  if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Labels loop:\n");
+  AliDebug(3,"Labels loop:");
+  
+  Bool_t conversion = kFALSE;
   
-  for(Int_t i = 0; i < nlabels; i++){
-    if(fDebug > 3) printf("\t  at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1);
+  for(Int_t i = 0; i < nlabels; i++)
+  {
+    AliDebug(3,Form("\t  at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d", i+1, nlabels, labels[i], okPhoton0, okPhoton1));
     
     //If we already found both, break the loop
     if(okPhoton0 && okPhoton1) break;
     
-    Int_t index =      labels[i];
-    if      (iPhoton0 == index) {
+    Int_t index = labels[i];
+    if      (iPhoton0 == index)
+    {
       okPhoton0 = kTRUE;
       continue;
     }
-    else if (iPhoton1 == index) {
+    else if (iPhoton1 == index)
+    {
       okPhoton1 = kTRUE;
       continue;
     }
     
     //Trace back the mother in case it was a conversion
     
-    if(index >= stack->GetNtrack()){
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) Particle index %d larger than size of list %d !!\n",index,stack->GetNtrack());
+    if(index >= stack->GetNtrack())
+    {
+      AliWarning(Form("Particle index %d larger than size of list %d!!",index,stack->GetNtrack()));
       continue;
     }
     
     TParticle * daught = stack->Particle(index);
     Int_t tmpindex = daught->GetFirstMother();         
-    if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex);
-    while(tmpindex>=0){
+    
+    AliDebug(3,Form("\t Conversion? : mother %d",tmpindex));
+    
+    while(tmpindex>=0)
+    {
       //MC particle of interest is the mother
-      if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex);
+      AliDebug(3,Form("\t \t parent index %d",tmpindex));
       daught   = stack->Particle(tmpindex);
-      if      (iPhoton0 == tmpindex) {
-        okPhoton0 = kTRUE;
+      if      (iPhoton0 == tmpindex)
+      {
+        conversion = kTRUE;
+        okPhoton0  = kTRUE;
         break;
       }
-      else if (iPhoton1 == tmpindex) {
-        okPhoton1 = kTRUE;
+      else if (iPhoton1 == tmpindex)
+      {
+        conversion = kTRUE;
+        okPhoton1  = kTRUE;
         break;
       }
+      
       tmpindex = daught->GetFirstMother();
+      
     }//While to check if pi0/eta daughter was one of these contributors to the cluster
     
-    if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0) 
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Something happens, first label should be from a photon decay!\n");
+    //if(i == 0 && (!okPhoton0 && !okPhoton1))
+    //  AliDebug(1,Form("Something happens, first label should be from a photon decay!"));
     
   }//loop on list of labels
   
   //If both photons contribute tag as the corresponding meson.
-  if(okPhoton0 && okPhoton1){
-    if(fDebug > 2) 
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - %s OVERLAPPED DECAY \n", meson->GetName());
+  if(okPhoton0 && okPhoton1)
+  {
+    AliDebug(2,Form("%s OVERLAPPED DECAY", meson->GetName()));
+    
+    if(!CheckTagBit(tag,kMCConversion) && conversion) SetTagBit(tag,kMCConversion) ;
     
     if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
-    else  SetTagBit(tag,kMCEta);
+    else                SetTagBit(tag,kMCEta);
   }
   
 }      
 
-//__________________________________________________________________________________
-void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
-                                                    const Int_t nlabels, 
-                                                    const Int_t mesonIndex, 
-                                                    const TClonesArray *mcparticles, 
-                                                    Int_t & tag )
+//________________________________________________________________________________________________________
+void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels, Int_t mesonIndex,
+                                                    const TClonesArray *mcparticles, Int_t & tag   )
 {
-  //Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles
+  // Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles.
   
   if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1)
   {
-    if(fDebug > 2) 
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n",
-             labels[0],mcparticles->GetEntriesFast(), nlabels);
+    AliDebug(2,Form("Exit : label[0] %d, n primaries %d, nlabels %d",labels[0],mcparticles->GetEntriesFast(), nlabels));
     return;
   }
   
@@ -871,18 +966,16 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
   Int_t mesonPdg = meson->GetPdgCode();
   if(mesonPdg != 111 && mesonPdg != 221)
   {
-    printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
+    AliWarning(Form("Wrong pi0/eta PDG : %d",mesonPdg));
     return;
   }
   
-  if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters());
-  
+  AliDebug(2,Form("pdg %d, label %d, ndaughters %d", mesonPdg, mesonIndex, meson->GetNDaughters()));
   
   //Get the daughters
   if(meson->GetNDaughters() != 2)
   {
-    if(fDebug > 2) 
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(stack) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters());
+    AliDebug(2,Form("Not overalapped. Number of daughters is %d, not 2",meson->GetNDaughters()));
     return;
   }
   
@@ -899,24 +992,23 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
   //Check if both daughters are photons
   if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22)
   {
-    printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overlapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
+    AliWarning(Form("Not overlapped. PDG:  daughter 1 = %d, of daughter 2 = %d",photon0->GetPdgCode(),photon1->GetPdgCode()));
     return;
   }
   
-  if(fDebug > 2) 
-    printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1);
+  AliDebug(2,Form("Daughter labels : photon0 = %d, photon1 = %d",iPhoton0,iPhoton1));
   
   //Check if both photons contribute to the cluster
   Bool_t okPhoton0 = kFALSE;
   Bool_t okPhoton1 = kFALSE;
   
-  if(fDebug > 3) 
-    printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n");
+  AliDebug(3,"Labels loop:");
+  
+  Bool_t conversion = kFALSE;
   
   for(Int_t i = 0; i < nlabels; i++)
   {
-    if(fDebug > 3)
-      printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1);
+    AliDebug(3, Form("\t label %d/%d: %d, ok? %d, %d", i, nlabels, labels[i], okPhoton0, okPhoton1));
     
     if(labels[i]<0) continue;
     
@@ -924,54 +1016,66 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
     if(okPhoton0 && okPhoton1) break;
     
     Int_t index =      labels[i];
-    if      (iPhoton0 == index) {
+    if      (iPhoton0 == index)
+    {
       okPhoton0 = kTRUE;
       continue;
     }
-    else if (iPhoton1 == index) {
+    else if (iPhoton1 == index)
+    {
       okPhoton1 = kTRUE;
       continue;
     }
     
     //Trace back the mother in case it was a conversion
     
-    if(index >= mcparticles->GetEntriesFast()){
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Particle index %d larger than size of list %d !!\n",index,mcparticles->GetEntriesFast());
+    if(index >= mcparticles->GetEntriesFast())
+    {
+      AliWarning(Form("Particle index %d larger than size of list %d!!",index,mcparticles->GetEntriesFast()));
       continue;
     }
     
     AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index);
     Int_t tmpindex = daught->GetMother();
-    if(fDebug > 3) 
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex);
+    AliDebug(3,Form("Conversion? : mother %d",tmpindex));
     
     while(tmpindex>=0){
       
       //MC particle of interest is the mother
-      if(fDebug > 3) 
-        printf("\t parent index %d\n",tmpindex);
+      AliDebug(3,Form("\t parent index %d",tmpindex));
       daught   = (AliAODMCParticle*) mcparticles->At(tmpindex);
       //printf("tmpindex %d\n",tmpindex);
-      if      (iPhoton0 == tmpindex) {
-        okPhoton0 = kTRUE;
+      if      (iPhoton0 == tmpindex)
+      {
+        conversion = kTRUE;
+        okPhoton0  = kTRUE;
         break;
       }
-      else if (iPhoton1 == tmpindex) {
-        okPhoton1 = kTRUE;
+      else if (iPhoton1 == tmpindex)
+      {
+        conversion = kTRUE;
+        okPhoton1  = kTRUE;
         break;
-      }                
+      }
+      
       tmpindex = daught->GetMother();
+      
     }//While to check if pi0/eta daughter was one of these contributors to the cluster
     
-    if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 )
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n");
+    //if(i == 0 && (!okPhoton0 && !okPhoton1)) AliDebug(1,"Something happens, first label should be from a photon decay!");
     
   }//loop on list of labels
   
   //If both photons contribute tag as the corresponding meson.
-  if(okPhoton0 && okPhoton1){
-    if(fDebug > 2) 
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName());
+  if(okPhoton0 && okPhoton1)
+  {
+    AliDebug(2,Form("%s OVERLAPPED DECAY",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName()));
+    
+    if(!CheckTagBit(tag,kMCConversion) && conversion)
+    {
+      AliDebug(2,"Second decay photon produced a conversion");
+      SetTagBit(tag,kMCConversion) ;
+    }
     
     if(mesonPdg == 111) SetTagBit(tag,kMCPi0);
     else                SetTagBit(tag,kMCEta);
@@ -979,10 +1083,224 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
   
 }
 
-//_________________________________________________________________________
+//______________________________________________________________________________________________________
+void    AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray* arrayCluster,   Int_t iMom, Int_t iParent,
+                                               AliStack * stack,                Int_t & tag)
+{
+  // Check on ESDs if the current decay photon has the second photon companion lost.
+  
+  if(!arrayCluster || iMom < 0 || iParent < 0|| !stack) return;
+  
+  TParticle * parent= stack->Particle(iParent);
+  
+  if(parent->GetNDaughters()!=2)
+  {
+    SetTagBit(tag, kMCDecayPairLost);
+    return ;
+  }
+  
+  Int_t pairLabel = -1;
+  if     ( iMom != parent->GetDaughter(0) ) pairLabel = parent->GetDaughter(0);
+  else if( iMom != parent->GetDaughter(1) ) pairLabel = parent->GetDaughter(1);
+  
+  if(pairLabel<0)
+  {
+    SetTagBit(tag, kMCDecayPairLost);
+    return ;
+  }
+  
+  for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++)
+  {
+    AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus);
+    for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++)
+    {
+      Int_t label = cluster->GetLabels()[ilab];
+      if(label==pairLabel)
+      {
+        SetTagBit(tag, kMCDecayPairInCalo);
+        return ;
+      }
+      else if(label== iParent || label== iMom)
+      {
+        continue;
+      }
+      else // check the ancestry
+      {
+        TParticle * mother = stack->Particle(label);
+        Int_t momPDG = TMath::Abs(mother->GetPdgCode());
+        if(momPDG!=11 && momPDG!=22) continue;
+        
+        //Check if "mother" of entity is converted, if not, get the first non converted mother
+        Int_t iParentClus = mother->GetFirstMother();
+        if(iParentClus < 0) continue;
+        
+        TParticle * parentClus = stack->Particle(iParentClus);
+        if(!parentClus) continue;
+        
+        Int_t parentClusPDG    = TMath::Abs(parentClus->GetPdgCode());
+        Int_t parentClusStatus = parentClus->GetStatusCode();
+        
+        if( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0) continue;
+        
+        //printf("Conversion\n");
+        
+        //Check if the mother is photon or electron with status not stable
+        while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1)
+        {
+          //New Mother
+          label            = iParentClus;
+          momPDG           = parentClusPDG;
+          
+          iParentClus      = parentClus->GetFirstMother();
+          if(iParentClus < 0) break;
+          
+          parentClus       = stack->Particle(iParentClus);
+          if(!parentClus) break;
+          
+          parentClusPDG    = TMath::Abs(parentClus->GetPdgCode());
+          parentClusStatus = parentClus->GetStatusCode() ;
+        }//while
+    
+        if((momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel))
+        {
+          SetTagBit(tag, kMCDecayPairInCalo);
+          //printf("Conversion is paired\n");
+          return ;
+        }
+        else continue;
+        
+      }
+    }
+  } // cluster loop
+  
+  SetTagBit(tag, kMCDecayPairLost);
+
+}
+
+//______________________________________________________________________________________________________
+void    AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray   * arrayCluster,Int_t iMom, Int_t iParent,
+                                               const TClonesArray* mcparticles, Int_t & tag)
+{
+  // Check on AODs if the current decay photon has the second photon companion lost.
+
+  if(!arrayCluster || iMom < 0 || iParent < 0|| !mcparticles) return;
+
+  AliAODMCParticle * parent = (AliAODMCParticle*) mcparticles->At(iParent);
+  
+  //printf("*** Check label %d with parent %d\n",iMom, iParent);
+  
+  if(parent->GetNDaughters()!=2)
+  {
+    SetTagBit(tag, kMCDecayPairLost);
+    //printf("\t ndaugh = %d\n",parent->GetNDaughters());
+    return ;
+  }
+  
+  Int_t pairLabel = -1;
+  if     ( iMom != parent->GetDaughter(0) ) pairLabel = parent->GetDaughter(0);
+  else if( iMom != parent->GetDaughter(1) ) pairLabel = parent->GetDaughter(1);
+  
+  if(pairLabel<0)
+  {
+    //printf("\t pair Label not found = %d\n",pairLabel);
+    SetTagBit(tag, kMCDecayPairLost);
+    return ;
+  }
+  
+  //printf("\t *** find pair %d\n",pairLabel);
+  
+  for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++)
+  {
+    AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus);
+    //printf("\t \t ** Cluster %d, nlabels %d\n",iclus,cluster->GetNLabels());
+    for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++)
+    {
+      Int_t label = cluster->GetLabels()[ilab];
+      
+      //printf("\t \t label %d\n",label);
+
+      if(label==pairLabel)
+      {
+        //printf("\t \t Pair found\n");
+        SetTagBit(tag, kMCDecayPairInCalo);
+        return ;
+      }
+      else if(label== iParent || label== iMom)
+      {
+        //printf("\t \t skip\n");
+        continue;
+      }
+      else // check the ancestry
+      {
+        AliAODMCParticle * mother = (AliAODMCParticle*) mcparticles->At(label);
+        Int_t momPDG = TMath::Abs(mother->GetPdgCode());
+        if(momPDG!=11 && momPDG!=22)
+        {
+          //printf("\t \t skip, pdg %d\n",momPDG);
+          continue;
+        }
+        
+        //Check if "mother" of entity is converted, if not, get the first non converted mother
+        Int_t iParentClus = mother->GetMother();
+        if(iParentClus < 0) continue;
+        
+        AliAODMCParticle * parentClus =  (AliAODMCParticle*) mcparticles->At(iParentClus);
+        if(!parentClus) continue;
+        
+        Int_t parentClusPDG    = TMath::Abs(parentClus->GetPdgCode());
+        Int_t parentClusStatus = parentClus->GetStatus();
+        
+        if( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0)
+        {
+          //printf("\t \t skip, not a conversion, parent: pdg %d, status %d\n",parentClusPDG,parentClusStatus);
+          continue;
+        }
+        
+        //printf("\t \t Conversion\n");
+        
+        //Check if the mother is photon or electron with status not stable
+        while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1)
+        {
+          //New Mother
+          label            = iParentClus;
+          momPDG           = parentClusPDG;
+          
+          iParentClus      = parentClus->GetMother();
+          if(iParentClus < 0) break;
+          
+          parentClus       =  (AliAODMCParticle*) mcparticles->At(iParentClus);
+          if(!parentClus) break;
+          
+          parentClusPDG    = TMath::Abs(parentClus->GetPdgCode());
+          parentClusStatus = parentClus->GetStatus() ;
+        }//while
+        
+        if((momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel))
+        {
+          SetTagBit(tag, kMCDecayPairInCalo);
+          //printf("\t \t Conversion is paired: mom %d, parent %d\n",label,iParentClus);
+          return ;
+        }
+        else
+        {
+          //printf("\t \t Skip, finally label %d, pdg %d, parent label %d, pdg %d, status %d\n",label,momPDG,iParentClus,parentClusPDG,parentClusStatus);
+          continue;
+        }
+        
+      }
+    }
+  } // cluster loop
+
+  
+  SetTagBit(tag, kMCDecayPairLost);
+  
+}
+
+//_____________________________________________________________________
 TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
 {
-  //Return list of jets (TParticles) and index of most likely parton that originated it.
+  // Return list of jets (TParticles) and index of most likely parton that originated it.
+  
   AliStack * stack = reader->GetStack();
   Int_t iEvent = reader->GetEventNumber();     
   AliGenEventHeader * geh = reader->GetGenEventHeader();
@@ -997,12 +1315,12 @@ TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
     if(stack->GetNtrack() < 8) return fJetsList;
     TParticle * parton1 =  stack->Particle(6);
     TParticle * parton2 =  stack->Particle(7);
-    if(fDebug > 2){
-      printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
-             parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta());
-      printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
-             parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta());
-               }
+    
+    AliDebug(2,Form("Parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
+                    parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta()));
+    AliDebug(2,Form("Parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
+                    parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta()));
+    
     //                 //Trace the jet from the mother parton
     //                 Float_t pt  = 0;
     //                 Float_t pt1 = 0;
@@ -1042,16 +1360,15 @@ TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
                
                //Get the jet, different way for different generator
                //PYTHIA
-    if(fMCGenerator == "PYTHIA"){
+    if(fMCGenerator == kPythia)
+    {
       TParticle * jet =  0x0;
       AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh;
       nTriggerJets =  pygeh->NTriggerJets();
-      if(fDebug > 1)
-        printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets);
+      AliDebug(2,Form("PythiaEventHeader: Njets: %d",nTriggerJets));
       
-      Int_t iparton = -1;
-      for(Int_t i = 0; i< nTriggerJets; i++){
-        iparton=-1;
+      for(Int_t i = 0; i< nTriggerJets; i++)
+      {
         pygeh->TriggerJet(i, tmpjet);
         jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0);
         //Assign an outgoing parton as mother
@@ -1060,18 +1377,19 @@ TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
         if(phidiff1 > phidiff2) jet->SetFirstMother(7);
         else  jet->SetFirstMother(6);
         //jet->Print();
-        if(fDebug > 1)
-          printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
-                 i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta());
+        AliDebug(1,Form("PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
+                        i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta()));
         fJetsList->Add(jet);                   
       }
     }//Pythia triggered jets
     //HERWIG
-    else if (fMCGenerator=="HERWIG"){
+    else if (fMCGenerator == kHerwig)
+    {
       Int_t pdg = -1;          
       //Check parton 1
       TParticle * tmp = parton1;
-      if(parton1->GetPdgCode()!=22){
+      if(parton1->GetPdgCode()!=22)
+      {
         while(pdg != 94){
           if(tmp->GetFirstDaughter()==-1) return fJetsList;
           tmp = stack->Particle(tmp->GetFirstDaughter());
@@ -1086,30 +1404,29 @@ TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
         //tmp = stack->Particle(tmp->GetFirstDaughter());
         //tmp->Print();
         //jet1->Print();
-        if(fDebug > 1)                 
-          printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
-                 jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta());
+        AliDebug(1,Form("HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
+                        jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta()));
       }//not photon
       
       //Check parton 2
       pdg = -1;
       tmp = parton2;
-      Int_t i = -1;
-      if(parton2->GetPdgCode()!=22){
-        while(pdg != 94){
+      if(parton2->GetPdgCode()!=22)
+      {
+        while(pdg != 94)
+        {
           if(tmp->GetFirstDaughter()==-1) return fJetsList;
-          i = tmp->GetFirstDaughter();
           tmp = stack->Particle(tmp->GetFirstDaughter());
           pdg = tmp->GetPdgCode();
         }//while
+        
         //Add found jet to list
         TParticle *jet2 = new TParticle(*tmp);
         jet2->SetFirstMother(7);
         fJetsList->Add(jet2);
         //jet2->Print();
-        if(fDebug > 1)
-          printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n",
-                 jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta());
+        AliDebug(2,Form("HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f",
+                        jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta()));
         //Int_t first =  tmp->GetFirstDaughter();
         //Int_t last  =  tmp->GetLastDaughter();
         //printf("jet 2:  first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode());
@@ -1128,45 +1445,135 @@ TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
 }
 
 
-//_______________________________________________________________________________________________
-TLorentzVector AliMCAnalysisUtils::GetMother(const Int_t label, const AliCaloTrackReader* reader,
-                                             Bool_t & ok) 
+//__________________________________________________________________________________________________________
+TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label,
+                                               const AliCaloTrackReader* reader,
+                                               Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & daughlabel)
+{
+  // Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother.
+  fDaughMom.SetPxPyPzE(0,0,0,0);
+  
+  if(reader->ReadStack())
+  {
+    if(!reader->GetStack())
+    {
+      AliWarning("Stack is not available, check analysis settings in configuration file!!");
+      
+      ok=kFALSE;
+      return fDaughMom;
+    }
+    
+    Int_t nprimaries = reader->GetStack()->GetNtrack();
+    if(label >= 0 && label < nprimaries)
+    {
+      TParticle * momP = reader->GetStack()->Particle(label);
+      daughlabel       = momP->GetDaughter(idaugh);
+      
+      if(daughlabel < 0 || daughlabel >= nprimaries)
+      {
+        ok = kFALSE;
+        return fDaughMom;
+      }
+      
+      TParticle * daughP = reader->GetStack()->Particle(daughlabel);
+      daughP->Momentum(fDaughMom);
+      pdg    = daughP->GetPdgCode();
+      status = daughP->GetStatusCode();
+    }
+    else
+    {
+      ok = kFALSE;
+      return fDaughMom;
+    }
+  }
+  else if(reader->ReadAODMCParticles())
+  {
+    TClonesArray* mcparticles = reader->GetAODMCParticles();
+    if(!mcparticles)
+    {
+      AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
+      
+      ok=kFALSE;
+      return fDaughMom;
+    }
+    
+    Int_t nprimaries = mcparticles->GetEntriesFast();
+    if(label >= 0 && label < nprimaries)
+    {
+      AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
+      daughlabel              = momP->GetDaughter(idaugh);
+      
+      if(daughlabel < 0 || daughlabel >= nprimaries)
+      {
+        ok = kFALSE;
+        return fDaughMom;
+      }
+      
+      AliAODMCParticle * daughP = (AliAODMCParticle *) mcparticles->At(daughlabel);
+      fDaughMom.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
+      pdg    = daughP->GetPdgCode();
+      status = daughP->GetStatus();
+    }
+    else
+    {
+      ok = kFALSE;
+      return fDaughMom;
+    }
+  }
+  
+  ok = kTRUE;
+  
+  return fDaughMom;
+}
+
+//______________________________________________________________________________________________________
+TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
+{
+  // Return the kinematics of the particle that generated the signal.
+  
+  Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1;
+  return GetMother(label,reader,pdg,status, ok,momlabel);
+}
+
+//_________________________________________________________________________________________
+TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader,
+                                             Int_t & pdg, Int_t & status, Bool_t & ok)
 {
-  //Return the kinematics of the particle that generated the signal
+  // Return the kinematics of the particle that generated the signal.
   
-  Int_t pdg = -1; Int_t status = -1;
-  return GetMother(label,reader,pdg,status, ok);
+  Int_t momlabel = -1;
+  return GetMother(label,reader,pdg,status, ok,momlabel);
 }
 
-//_______________________________________________________________________________________________
-TLorentzVector AliMCAnalysisUtils::GetMother(const Int_t label, const AliCaloTrackReader* reader, 
-                                             Int_t & pdg, Int_t & status, Bool_t & ok
+//______________________________________________________________________________________________________
+TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, 
+                                             Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel)
 {
-  //Return the kinematics of the particle that generated the signal, its pdg and its status
+  // Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother.
   
-  TLorentzVector mom(0,0,0,0);
+  fMotherMom.SetPxPyPzE(0,0,0,0);
   
   if(reader->ReadStack())
   {
     if(!reader->GetStack()) 
     {
-      if (fDebug >=0) 
-        printf("AliMCAnalysisUtils::GetMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
+      AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
      
       ok=kFALSE;
-      return mom;
+      return fMotherMom;
     }
     if(label >= 0 && label < reader->GetStack()->GetNtrack())
     {
       TParticle * momP = reader->GetStack()->Particle(label);
-      momP->Momentum(mom);
-      pdg    = momP->GetPdgCode();
-      status = momP->GetStatusCode();
+      momP->Momentum(fMotherMom);
+      pdg      = momP->GetPdgCode();
+      status   = momP->GetStatusCode();
+      momlabel = momP->GetFirstMother();
     } 
     else 
     {
       ok = kFALSE;
-      return mom;
+      return fMotherMom;
     }
   }
   else if(reader->ReadAODMCParticles())
@@ -1174,52 +1581,52 @@ TLorentzVector AliMCAnalysisUtils::GetMother(const Int_t label, const AliCaloTra
     TClonesArray* mcparticles = reader->GetAODMCParticles();
     if(!mcparticles) 
     {
-      if(fDebug >= 0)
-        printf("AliMCAnalysisUtils::GetMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
+      AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
       
       ok=kFALSE;
-      return mom;
+      return fMotherMom;
     }
     
     Int_t nprimaries = mcparticles->GetEntriesFast();
     if(label >= 0 && label < nprimaries)
     {
       AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
-      mom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
-      pdg    = momP->GetPdgCode();
-      status = momP->GetStatus();
-    }     
+      fMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E());
+      pdg      = momP->GetPdgCode();
+      status   = momP->GetStatus();
+      momlabel = momP->GetMother();
+    }
     else 
     {
       ok = kFALSE;
-      return mom;
+      return fMotherMom;
     }
   }
   
   ok = kTRUE;
   
-  return mom;
+  return fMotherMom;
 }
 
-
-//_____________________________________________________________________________________
-TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int_t pdg, const AliCaloTrackReader* reader, Bool_t & ok) 
+//___________________________________________________________________________________
+TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(Int_t label, Int_t pdg,
+                                                    const AliCaloTrackReader* reader,
+                                                    Bool_t & ok, Int_t & momlabel)
 {
-  //Return the kinematics of the particle that generated the signal
-  
-  TLorentzVector grandmom(0,0,0,0);
+  // Return the kinematics of the particle that generated the signal.
   
+  fGMotherMom.SetPxPyPzE(0,0,0,0);
   
   if(reader->ReadStack())
   {
     if(!reader->GetStack())
     {
-      if (fDebug >=0) 
-        printf("AliMCAnalysisUtils::GetMotherWithPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
+     AliWarning("Stack is not available, check analysis settings in configuration file!!");
       
       ok = kFALSE;
-      return grandmom;
+      return fGMotherMom;
     }
+    
     if(label >= 0 && label < reader->GetStack()->GetNtrack())
     {
       TParticle * momP = reader->GetStack()->Particle(label);
@@ -1233,7 +1640,8 @@ TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int
         grandmomPDG = grandmomP->GetPdgCode();
         if(grandmomPDG==pdg)
         {
-          grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
+          momlabel = grandmomLabel;
+          fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
           break;
         }
         
@@ -1241,7 +1649,7 @@ TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int
         
       }
       
-      if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(ESD) - mother with PDG %d, not found! \n",pdg);
+      if(grandmomPDG!=pdg) AliInfo(Form("Mother with PDG %d, not found! \n",pdg));
     }
   }
   else if(reader->ReadAODMCParticles())
@@ -1249,11 +1657,10 @@ TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int
     TClonesArray* mcparticles = reader->GetAODMCParticles();
     if(!mcparticles) 
     {
-      if(fDebug >= 0)
-        printf("AliMCAnalysisUtils::GetMotherWithPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
+      AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
       
       ok=kFALSE;
-      return grandmom;
+      return fGMotherMom;
     }
     
     Int_t nprimaries = mcparticles->GetEntriesFast();
@@ -1271,8 +1678,8 @@ TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int
         if(grandmomPDG==pdg)
         {
           //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
-
-          grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
+          momlabel = grandmomLabel;
+          fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
           break;
         }
         
@@ -1280,32 +1687,106 @@ TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int
         
       }
       
-      if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d, NOT found! \n",pdg);
+      if(grandmomPDG!=pdg) AliInfo(Form("Mother with PDG %d, NOT found!",pdg));
             
     }
   }
   
   ok = kTRUE;
   
-  return grandmom;
+  return fGMotherMom;
 }
 
-//_____________________________________________________________________________________
-Float_t AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(const Int_t label, const Int_t pdg, const AliCaloTrackReader* reader, Bool_t & ok) 
+//______________________________________________________________________________________________
+TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliCaloTrackReader* reader,
+                                                  Int_t & pdg, Int_t & status, Bool_t & ok,
+                                                  Int_t & grandMomLabel, Int_t & greatMomLabel)
 {
-  //In case of an eta or pi0 decay into 2 photons, get the asymmetry  in the energy of the photons
+  // Return the kinematics of the particle that generated the signal.
   
-  Float_t asym = -2;
+  fGMotherMom.SetPxPyPzE(0,0,0,0);
   
   if(reader->ReadStack())
   {
     if(!reader->GetStack())
     {
-      if (fDebug >=0) 
-        printf("AliMCAnalysisUtils::GetMotherWithPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
+      AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
+      
+      ok = kFALSE;
+      return fGMotherMom;
+    }
+    
+    if(label >= 0 && label < reader->GetStack()->GetNtrack())
+    {
+      TParticle * momP = reader->GetStack()->Particle(label);
+      
+      grandMomLabel = momP->GetFirstMother();
+      
+      TParticle * grandmomP = 0x0;
+      
+      if (grandMomLabel >=0 )
+      {
+        grandmomP   = reader->GetStack()->Particle(grandMomLabel);
+        pdg    = grandmomP->GetPdgCode();
+        status = grandmomP->GetStatusCode();
+       
+        fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
+        greatMomLabel =  grandmomP->GetFirstMother();
+
+      }
+    }
+  }
+  else if(reader->ReadAODMCParticles())
+  {
+    TClonesArray* mcparticles = reader->GetAODMCParticles();
+    if(!mcparticles)
+    {
+      AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
+      
+      ok=kFALSE;
+      return fGMotherMom;
+    }
+    
+    Int_t nprimaries = mcparticles->GetEntriesFast();
+    if(label >= 0 && label < nprimaries)
+    {
+      AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
+      
+      grandMomLabel = momP->GetMother();
+      
+      AliAODMCParticle * grandmomP = 0x0;
+      
+      if(grandMomLabel >=0 )
+      {
+        grandmomP   = (AliAODMCParticle *) mcparticles->At(grandMomLabel);
+        pdg    = grandmomP->GetPdgCode();
+        status = grandmomP->GetStatus();
+      
+        fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
+        greatMomLabel =  grandmomP->GetMother();
+        
+      }
+    }
+  }
+  
+  ok = kTRUE;
+  
+  return fGMotherMom;
+}
+
+//_______________________________________________________________________________________________________________
+void AliMCAnalysisUtils::GetMCDecayAsymmetryAngleForPDG(Int_t label, Int_t pdg, const AliCaloTrackReader* reader,
+                                                        Float_t & asym, Float_t & angle, Bool_t & ok)
+{
+  // In case of an eta or pi0 decay into 2 photons, get the asymmetry  in the energy of the photons.
+  
+  if(reader->ReadStack())
+  {
+    if(!reader->GetStack())
+    {
+      AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
       
       ok = kFALSE;
-      return asym;
     }
     if(label >= 0 && label < reader->GetStack()->GetNtrack())
     {
@@ -1332,12 +1813,15 @@ Float_t AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(const Int_t label, const I
         if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
         {
           asym = (d1->Energy()-d2->Energy())/grandmomP->Energy();
+          d1->Momentum(fDaughMom );
+          d2->Momentum(fDaughMom2);
+          angle = fDaughMom.Angle(fDaughMom2.Vect());
         }
       }
       else 
       {
         ok=kFALSE;
-        printf("AliMCAnalysisUtils::GetMotherWithPDG(ESD) - mother with PDG %d, not found! \n",pdg);
+        AliInfo(Form("Mother with PDG %d, not found!",pdg));
       }
       
       } // good label
@@ -1347,11 +1831,10 @@ Float_t AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(const Int_t label, const I
     TClonesArray* mcparticles = reader->GetAODMCParticles();
     if(!mcparticles) 
     {
-      if(fDebug >= 0)
-        printf("AliMCAnalysisUtils::GetMotherWithPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n");
+      AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
       
       ok=kFALSE;
-      return asym;
+      return;
     }
     
     Int_t nprimaries = mcparticles->GetEntriesFast();
@@ -1380,12 +1863,15 @@ Float_t AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(const Int_t label, const I
         if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22)
         {
           asym = (d1->E()-d2->E())/grandmomP->E();
+          fDaughMom .SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E());
+          fDaughMom2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E());
+          angle = fDaughMom.Angle(fDaughMom2.Vect());
         }
       }
       else 
       {
         ok=kFALSE;
-        printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d, not found! \n",pdg);
+        AliInfo(Form("Mother with PDG %d, not found! \n",pdg));
       }
       
     } // good label
@@ -1393,14 +1879,148 @@ Float_t AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(const Int_t label, const I
   
   ok = kTRUE;
   
-  return asym;
 }
 
+//_________________________________________________________________________________________________
+Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
+{
+  // Return the the number of daughters of a given MC particle.
+  
+  if(reader->ReadStack())
+  {
+    if(!reader->GetStack())
+    {
+      AliWarning("Stack is not available, check analysis settings in configuration file, STOP!!");
+      
+      ok=kFALSE;
+      return -1;
+    }
+    if(label >= 0 && label < reader->GetStack()->GetNtrack())
+    {
+      TParticle * momP = reader->GetStack()->Particle(label);
+      ok=kTRUE;
+      return momP->GetNDaughters();
+    }
+    else
+    {
+      ok = kFALSE;
+      return -1;
+    }
+  }
+  else if(reader->ReadAODMCParticles())
+  {
+    TClonesArray* mcparticles = reader->GetAODMCParticles();
+    if(!mcparticles)
+    {
+      AliWarning("AODMCParticles is not available, check analysis settings in configuration file!!");
+      
+      ok=kFALSE;
+      return -1;
+    }
+    
+    Int_t nprimaries = mcparticles->GetEntriesFast();
+    if(label >= 0 && label < nprimaries)
+    {
+      AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label);
+      ok = kTRUE;
+      return momP->GetNDaughters();
+    }
+    else
+    {
+      ok = kFALSE;
+      return -1;
+    }
+  }
+  
+  ok = kFALSE;
+  
+  return -1;
+}
+
+//_________________________________________________________________________________
+Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, UInt_t nlabels,
+                                       Int_t mctag, Int_t mesonLabel,
+                                       AliCaloTrackReader * reader, Int_t *overpdg)
+{
+  // Compare the primary depositing more energy with the rest,
+  // if no photon/electron (conversion) or neutral meson as comon ancestor, consider it as other particle contributing.
+  // Give as input the meson label in case it was a pi0 or eta merged cluster.
+  // Init overpdg with nlabels.
+  
+  Int_t ancPDG = 0, ancStatus = -1;
+  TVector3 prodVertex;
+  Int_t ancLabel = 0;
+  Int_t noverlaps = 0;
+  Bool_t ok = kFALSE;
+  
+  for (UInt_t ilab = 1; ilab < nlabels; ilab++ )
+  {
+    ancLabel = CheckCommonAncestor(label[0],label[ilab],reader,ancPDG,ancStatus,fMotherMom,prodVertex);
+    
+    //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n",
+    //       ilab,label[0],label[ilab],ancLabel,ancPDG, mctag);
+    
+    Bool_t overlap = kFALSE;
+    
+    if     ( ancLabel < 0 )
+    {
+      overlap = kTRUE;
+      //printf("\t \t \t No Label = %d\n",ancLabel);
+    }
+    else if( ( ancPDG==111 || ancPDG==221 ) && ( CheckTagBit(mctag,kMCPi0) ||  CheckTagBit(mctag,kMCEta)) && mesonLabel != ancLabel)
+    {
+      //printf("\t \t  meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel);
+      overlap = kTRUE;
+    }
+    else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 )
+    {
+      //printf("\t \t \t Non EM PDG = %d\n",ancPDG);
+      overlap = kTRUE ;
+    }
+    
+    if( !overlap ) continue ;
+    
+    // We have at least one overlap
+    
+    //printf("Overlap!!!!!!!!!!!!!!\n");
+    
+    noverlaps++;
+    
+    // What is the origin of the overlap?
+    Bool_t  mOK = 0,      gOK = 0;
+    Int_t   mpdg = -999999,  gpdg = -1;
+    Int_t   mstatus = -1, gstatus = -1;
+    Int_t   gLabel = -1, ggLabel = -1;
+    
+    GetMother     (label[ilab],reader,mpdg,mstatus,mOK);
+    fGMotherMom =
+    GetGrandMother(label[ilab],reader,gpdg,gstatus,gOK, gLabel,ggLabel);
+    
+    //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg);
+    
+    if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) &&
+        ( gpdg == 22 || TMath::Abs(gpdg==11) ) &&
+       gLabel >=0 )
+    {
+      Int_t labeltmp = gLabel;
+      while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
+      {
+        mpdg=gpdg;
+        fGMotherMom = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel);
+        labeltmp=gLabel;
+      }
+    }
+    overpdg[noverlaps-1] = mpdg;
+  }
+  
+  return noverlaps ;
+  
+}
 
 //________________________________________________________
 void AliMCAnalysisUtils::Print(const Option_t * opt) const
 {
-  //Print some relevant parameters set for the analysis
+  // Print some relevant parameters set for the analysis.
   
   if(! opt)
     return;
@@ -1408,17 +2028,17 @@ void AliMCAnalysisUtils::Print(const Option_t * opt) const
   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
   
   printf("Debug level    = %d\n",fDebug);
-  printf("MC Generator   = %s\n",fMCGenerator.Data());
+  printf("MC Generator   = %s\n",fMCGeneratorString.Data());
   printf(" \n");
   
 } 
 
-//________________________________________________________
-void AliMCAnalysisUtils::PrintMCTag(const Int_t tag) const
+//__________________________________________________
+void AliMCAnalysisUtils::PrintMCTag(Int_t tag) const
 {
-  // print the assigned origins to this particle
+  // Print the assigned origins to this particle.
   
-  printf("AliMCAnalysisUtils::PrintMCTag() - tag %d \n    photon %d, conv %d, prompt %d, frag %d, isr %d, \n    pi0 decay %d, eta decay %d, other decay %d  pi0 %d,  eta %d \n    electron %d, muon %d,pion %d, proton %d, neutron %d, \n    kaon %d, a-proton %d, a-neutron %d, other %d, unk %d, bad %d\n",
+  printf("AliMCAnalysisUtils::PrintMCTag() - tag %d \n    photon %d, conv %d, prompt %d, frag %d, isr %d, \n    pi0 decay %d, eta decay %d, other decay %d  pi0 %d,  eta %d \n    electron %d, muon %d,pion %d, proton %d, neutron %d, \n    kaon %d, a-proton %d, a-neutron %d, unk %d, bad %d\n",
          tag,
          CheckTagBit(tag,kMCPhoton),
          CheckTagBit(tag,kMCConversion),
@@ -1438,11 +2058,44 @@ void AliMCAnalysisUtils::PrintMCTag(const Int_t tag) const
          CheckTagBit(tag,kMCKaon), 
          CheckTagBit(tag,kMCAntiProton), 
          CheckTagBit(tag,kMCAntiNeutron),
-         CheckTagBit(tag,kMCOther),
          CheckTagBit(tag,kMCUnknown),
          CheckTagBit(tag,kMCBadLabel)
          );
-  
 } 
 
+//__________________________________________________
+void AliMCAnalysisUtils::SetMCGenerator(Int_t mcgen)
+{
+  // Set the generator type.
+  
+  fMCGenerator = mcgen ;
+  if     (mcgen == kPythia) fMCGeneratorString = "PYTHIA";
+  else if(mcgen == kHerwig) fMCGeneratorString = "HERWIG";
+  else if(mcgen == kHijing) fMCGeneratorString = "HIJING";
+  else
+  {
+    fMCGeneratorString = "";
+    fMCGenerator       = kBoxLike ;
+  }
+  
+}
+
+//__________________________________________________
+void AliMCAnalysisUtils::SetMCGenerator(TString mcgen)
+{
+  // Set the generator type.
+  
+  fMCGeneratorString = mcgen ;
+  
+  if     (mcgen == "PYTHIA") fMCGenerator = kPythia;
+  else if(mcgen == "HERWIG") fMCGenerator = kHerwig;
+  else if(mcgen == "HIJING") fMCGenerator = kHijing;
+  else
+  {
+    fMCGenerator       = kBoxLike;
+    fMCGeneratorString = "" ;
+  }
+}
+
+