]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG/CaloTrackCorrBase/AliMCAnalysisUtils.cxx
remove commented line
[u/mrichter/AliRoot.git] / PWG / CaloTrackCorrBase / AliMCAnalysisUtils.cxx
index 0171460dd35cd04f45c1c051f997a23b044f0d79..61b37022be18daba3a3c3448149e1a7b9a47214f 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;
@@ -153,8 +159,8 @@ 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;
@@ -212,59 +218,68 @@ Int_t AliMCAnalysisUtils::CheckCommonAncestor(const Int_t index1, const Int_t in
   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.
@@ -276,8 +291,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
   
   if(!stack)
   {
-    if (fDebug >=0) 
-      printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n");
+    AliDebug(1,"Stack is not available, check analysis settings in configuration file, STOP!!");
     return -1;
   }
   
@@ -294,7 +308,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
     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;
@@ -310,15 +324,12 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
         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)
     {
@@ -335,7 +346,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
         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)
@@ -354,12 +365,9 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
         }
       }//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);
-      }
+      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)
@@ -374,11 +382,9 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
         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.
@@ -408,17 +414,16 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
       
       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)
     {
       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
     }
@@ -431,21 +436,29 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
       {
         SetTagBit(tag,kMCPi0Decay);
         
-        if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon,  parent pi0 with status 11 \n");
+        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);
         
-        if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon,  parent pi0 with status 11 \n");
+        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 == "PYTHIA")
+        if(fMCGenerator == kPythia)
         {
           if(iParent < 8 && iParent > 5)
           {//outgoing partons
@@ -459,7 +472,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
           else  SetTagBit(tag,kMCUnknown);
          }//PYTHIA
       
-        else if(fMCGenerator == "HERWIG")
+        else if(fMCGenerator == kHerwig)
         {
           if(pStatus < 197)
           {//Not decay
@@ -506,12 +519,13 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
       
       SetTagBit(tag,kMCElectron);
       
-      if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n");
+      AliDebug(1,"Checking ancestors of electrons");
      
-      if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); } //Pi0 Dalitz decay
-      else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); } //Eta Dalitz decay
+      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
+      else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000))
+      { //check charm decay
         if(parent)
         {
           Int_t iGrandma = parent->GetFirstMother();
@@ -520,7 +534,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
             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
         }//parent
@@ -531,27 +545,26 @@ Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels,
         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);
+      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 < 0)
+      AliWarning(Form("*** bad label or no stack ***:  label %d ", label));
     
-    if(label >=  stack->GetNtrack() &&  (fDebug >= 0))
-      printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
+    if(label >=  stack->GetNtrack())
+      AliWarning(Form("*** large label ***:  label %d, n tracks %d", label, stack->GetNtrack()));
     
     SetTagBit(tag,kMCUnknown);
   }//Bad label
@@ -561,17 +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");
+    AliDebug(1,"AODMCParticles is not available, check analysis settings in configuration file!!");
     return -1;
   }
        
@@ -588,7 +600,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
     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 ;
@@ -598,20 +610,12 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
       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->IsPrimary(), parent->IsPhysicalPrimary()));
     
-    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())
     {
@@ -626,7 +630,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
         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)
@@ -639,17 +643,9 @@ 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->IsPrimary(), parent->IsPhysicalPrimary()));
       
     }//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())
@@ -664,11 +660,9 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
         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.
@@ -700,7 +694,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
     {
       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
     }
@@ -708,7 +702,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
     {
       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
     }
@@ -721,26 +715,36 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
       {
         SetTagBit(tag,kMCPi0Decay);
         
-        if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n");
+        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);
         
-        if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n");
+        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=="PYTHIA" || fMCGenerator=="HERWIG")) //undecayed particle
+      else if( mom->IsPhysicalPrimary() && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) ) //undecayed particle
       {
         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
         }
@@ -768,9 +772,10 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
       
       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
+      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
@@ -791,8 +796,7 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
         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);
@@ -801,18 +805,17 @@ Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels,
     //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);
+      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 < 0)
+      AliWarning(Form("*** bad label or no mcparticles ***:  label %d", label));
     
-    if(label >=  mcparticles->GetEntriesFast() &&  (fDebug >= 0) )
-      printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***:  label %d, n tracks %d \n", label, mcparticles->GetEntriesFast());
+    if(label >=  mcparticles->GetEntriesFast())
+      AliWarning(Form("*** large label ***:  label %d, n tracks %d", label, mcparticles->GetEntriesFast()));
   
     SetTagBit(tag,kMCUnknown);
     
@@ -822,18 +825,16 @@ 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(ESD) - 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;
   }
   
@@ -841,17 +842,16 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
   Int_t mesonPdg    = meson->GetPdgCode();
   if(mesonPdg!=111 && mesonPdg!=221)
   {
-    printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Wrong pi0/eta PDG : %d \n",mesonPdg);
+    AliWarning(Form("Wrong pi0/eta PDG : %d",mesonPdg));
     return;
   }
   
-  if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %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(ESD) - 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;
   }
   
@@ -864,25 +864,23 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
   //Check if both daughters are photons
   if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22)
   {
-    if(fDebug > 2) 
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. PDG:  daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode());
+    AliDebug(2,Form("Not overalapped. PDG:  daughter 1 = %d, of daughter 2 = %d",photon0->GetPdgCode(),photon1->GetPdgCode()));
     return;
   }
   
-  if(fDebug > 2) 
-    printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - 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(ESD) - 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);
+    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;
@@ -903,18 +901,19 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
     
     if(index >= stack->GetNtrack())
     {
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) Particle index %d larger than size of list %d !!\n",
-             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);
+    
+    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)
       {
@@ -933,16 +932,15 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
       
     }//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(ESD) - 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(ESD) - %s OVERLAPPED DECAY \n", meson->GetName());
+    AliDebug(2,Form("%s OVERLAPPED DECAY", meson->GetName()));
     
     if(!CheckTagBit(tag,kMCConversion) && conversion) SetTagBit(tag,kMCConversion) ;
     
@@ -952,20 +950,15 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
   
 }      
 
-//__________________________________________________________________________________
-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;
   }
   
@@ -973,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(ESD) - 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;
   }
   
@@ -1001,26 +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;
     
@@ -1043,20 +1031,18 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
     
     if(index >= mcparticles->GetEntriesFast())
     {
-      printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Particle index %d larger than size of list %d !!\n",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)
@@ -1076,22 +1062,18 @@ void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels,
       
     }//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());
+    AliDebug(2,Form("%s OVERLAPPED DECAY",(TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName()));
     
     if(!CheckTagBit(tag,kMCConversion) && conversion)
     {
-      if(fDebug > 2)
-        printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Second decay photon produced a conversion \n");
+      AliDebug(2,"Second decay photon produced a conversion");
       SetTagBit(tag,kMCConversion) ;
     }
     
@@ -1101,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();
@@ -1119,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;
@@ -1164,17 +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
@@ -1183,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());
@@ -1209,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());
@@ -1251,39 +1445,45 @@ TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader)
 }
 
 
-//________________________________________________________________________________________________________
-TLorentzVector AliMCAnalysisUtils::GetDaughter(const Int_t idaugh, const Int_t label,
+//__________________________________________________________________________________________________________
+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
-  
-  TLorentzVector daughter(0,0,0,0);
+  // 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())
     {
-      if (fDebug >=0)
-        printf("AliMCAnalysisUtils::GetDaughter() - 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 daughter;
+      return fDaughMom;
     }
-    if(label >= 0 && label < reader->GetStack()->GetNtrack())
+    
+    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(daughter);
+      daughP->Momentum(fDaughMom);
       pdg    = daughP->GetPdgCode();
       status = daughP->GetStatusCode();
     }
     else
     {
       ok = kFALSE;
-      return daughter;
+      return fDaughMom;
     }
   }
   else if(reader->ReadAODMCParticles())
@@ -1291,11 +1491,10 @@ TLorentzVector AliMCAnalysisUtils::GetDaughter(const Int_t idaugh, const Int_t l
     TClonesArray* mcparticles = reader->GetAODMCParticles();
     if(!mcparticles)
     {
-      if(fDebug >= 0)
-        printf("AliMCAnalysisUtils::GetDaughter() - 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 daughter;
+      return fDaughMom;
     }
     
     Int_t nprimaries = mcparticles->GetEntriesFast();
@@ -1304,73 +1503,77 @@ TLorentzVector AliMCAnalysisUtils::GetDaughter(const Int_t idaugh, const Int_t l
       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);
-      daughter.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
+      fDaughMom.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E());
       pdg    = daughP->GetPdgCode();
       status = daughP->GetStatus();
     }
     else
     {
       ok = kFALSE;
-      return daughter;
+      return fDaughMom;
     }
   }
   
   ok = kTRUE;
   
-  return daughter;
+  return fDaughMom;
 }
 
-//_______________________________________________________________________________________________
-TLorentzVector AliMCAnalysisUtils::GetMother(const Int_t label, const AliCaloTrackReader* reader,
-                                             Bool_t & ok) 
+//______________________________________________________________________________________________________
+TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, 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; Int_t momlabel = -1;
   return GetMother(label,reader,pdg,status, ok,momlabel);
 }
 
-//_______________________________________________________________________________________________
-TLorentzVector AliMCAnalysisUtils::GetMother(const Int_t label, const AliCaloTrackReader* reader,
+//_________________________________________________________________________________________
+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 momlabel = -1;
   return GetMother(label,reader,pdg,status, ok,momlabel);
 }
 
-//_______________________________________________________________________________________________
-TLorentzVector AliMCAnalysisUtils::GetMother(const Int_t label, const AliCaloTrackReader* reader, 
+//______________________________________________________________________________________________________
+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 and its label mother
+  // 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())
@@ -1378,52 +1581,50 @@ 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,
+//___________________________________________________________________________________
+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())
@@ -1440,7 +1641,7 @@ TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int
         if(grandmomPDG==pdg)
         {
           momlabel = grandmomLabel;
-          grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
+          fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
           break;
         }
         
@@ -1448,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())
@@ -1456,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();
@@ -1479,7 +1679,7 @@ TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(const Int_t label, const Int
         {
           //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg);
           momlabel = grandmomLabel;
-          grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
+          fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
           break;
         }
         
@@ -1487,34 +1687,33 @@ 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;
 }
 
-//____________________________________________________________________________________________________
-TLorentzVector AliMCAnalysisUtils::GetGrandMother(const Int_t label, const AliCaloTrackReader* reader,
+//______________________________________________________________________________________________
+TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliCaloTrackReader* reader,
                                                   Int_t & pdg, Int_t & status, Bool_t & ok,
                                                   Int_t & grandMomLabel, Int_t & greatMomLabel)
 {
-  //Return the kinematics of the particle that generated the signal
+  // Return the kinematics of the particle that generated the signal.
   
-  TLorentzVector grandmom(0,0,0,0);
+  fGMotherMom.SetPxPyPzE(0,0,0,0);
   
   if(reader->ReadStack())
   {
     if(!reader->GetStack())
     {
-      if (fDebug >=0)
-        printf("AliMCAnalysisUtils::GetGrandMother() - 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 grandmom;
+      return fGMotherMom;
     }
     
     if(label >= 0 && label < reader->GetStack()->GetNtrack())
@@ -1531,7 +1730,7 @@ TLorentzVector AliMCAnalysisUtils::GetGrandMother(const Int_t label, const AliCa
         pdg    = grandmomP->GetPdgCode();
         status = grandmomP->GetStatusCode();
        
-        grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
+        fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy());
         greatMomLabel =  grandmomP->GetFirstMother();
 
       }
@@ -1542,11 +1741,10 @@ TLorentzVector AliMCAnalysisUtils::GetGrandMother(const Int_t label, const AliCa
     TClonesArray* mcparticles = reader->GetAODMCParticles();
     if(!mcparticles)
     {
-      if(fDebug >= 0)
-        printf("AliMCAnalysisUtils::GetGrandMother() - 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();
@@ -1564,7 +1762,7 @@ TLorentzVector AliMCAnalysisUtils::GetGrandMother(const Int_t label, const AliCa
         pdg    = grandmomP->GetPdgCode();
         status = grandmomP->GetStatus();
       
-        grandmom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
+        fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E());
         greatMomLabel =  grandmomP->GetMother();
         
       }
@@ -1573,25 +1771,22 @@ TLorentzVector AliMCAnalysisUtils::GetGrandMother(const Int_t label, const AliCa
   
   ok = kTRUE;
   
-  return grandmom;
+  return fGMotherMom;
 }
 
-//_____________________________________________________________________________________
-Float_t AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(const Int_t label, const Int_t pdg, const AliCaloTrackReader* reader, Bool_t & ok) 
+//_______________________________________________________________________________________________________________
+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
-  
-  Float_t asym = -2;
+  // 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())
     {
-      if (fDebug >=0) 
-        printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - 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 asym;
     }
     if(label >= 0 && label < reader->GetStack()->GetNtrack())
     {
@@ -1618,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::GetMCDecayAsymmetryForPDG(ESD) - mother with PDG %d, not found! \n",pdg);
+        AliInfo(Form("Mother with PDG %d, not found!",pdg));
       }
       
       } // good label
@@ -1633,11 +1831,10 @@ Float_t AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(const Int_t label, const I
     TClonesArray* mcparticles = reader->GetAODMCParticles();
     if(!mcparticles) 
     {
-      if(fDebug >= 0)
-        printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - 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();
@@ -1666,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::GetMCDecayAsymmetryForPDG(AOD) - mother with PDG %d, not found! \n",pdg);
+        AliInfo(Form("Mother with PDG %d, not found! \n",pdg));
       }
       
     } // good label
@@ -1679,22 +1879,18 @@ Float_t AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(const Int_t label, const I
   
   ok = kTRUE;
   
-  return asym;
 }
 
-
-//_______________________________________________________________________________________________________
-Int_t AliMCAnalysisUtils::GetNDaughters(const Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
+//_________________________________________________________________________________________________
+Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok)
 {
-  // Return the the number of daughters of a given MC particle
-  
+  // Return the the number of daughters of a given MC particle.
   
   if(reader->ReadStack())
   {
     if(!reader->GetStack())
     {
-      if (fDebug >=0)
-        printf("AliMCAnalysisUtils::GetNDaughters() - 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 -1;
@@ -1716,8 +1912,7 @@ Int_t AliMCAnalysisUtils::GetNDaughters(const Int_t label, const AliCaloTrackRea
     TClonesArray* mcparticles = reader->GetAODMCParticles();
     if(!mcparticles)
     {
-      if(fDebug >= 0)
-        printf("AliMCAnalysisUtils::GetNDaughters() - 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 -1;
@@ -1742,25 +1937,25 @@ Int_t AliMCAnalysisUtils::GetNDaughters(const Int_t label, const AliCaloTrackRea
   return -1;
 }
 
-//_______________________________________________________________________________
-Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, const UInt_t nlabels,
-                                       const Int_t mctag, const Int_t mesonLabel,
+//_________________________________________________________________________________
+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
+  // 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;
-  TLorentzVector momentum; TVector3 prodVertex;
+  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,momentum,prodVertex);
+    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);
@@ -1796,8 +1991,10 @@ Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, const UInt_t nlabels
     Int_t   mpdg = -999999,  gpdg = -1;
     Int_t   mstatus = -1, gstatus = -1;
     Int_t   gLabel = -1, ggLabel = -1;
-    TLorentzVector mother      = GetMother     (label[ilab],reader,mpdg,mstatus,mOK);
-    TLorentzVector grandmother = GetGrandMother(label[ilab],reader,gpdg,gstatus,gOK, gLabel,ggLabel);
+    
+    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);
     
@@ -1809,7 +2006,7 @@ Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, const UInt_t nlabels
       while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 )
       {
         mpdg=gpdg;
-        grandmother = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel);
+        fGMotherMom = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel);
         labeltmp=gLabel;
       }
     }
@@ -1823,7 +2020,7 @@ Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, const UInt_t nlabels
 //________________________________________________________
 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;
@@ -1831,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),
@@ -1861,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 = "" ;
+  }
+}
+
+