]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTrackCorrelations/AliAnaInsideClusterInvariantMass.cxx
Add histograms related to asymmetry, add option of asymetry cut for splitting
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaInsideClusterInvariantMass.cxx
index aa2394d8e27526515c5895072e0e7f918b04b4ee..6dbbb2b118af530f070c390953f9262e7f4a3ae1 100755 (executable)
 #include <TList.h>
 #include <TClonesArray.h>
 #include <TObjString.h>
-#include <TH3F.h>
-#include <TCanvas.h>
-#include <TStyle.h>
-#include <TPad.h>
+#include <TH2F.h>
 
 // --- Analysis system --- 
 #include "AliAnaInsideClusterInvariantMass.h" 
@@ -57,11 +54,14 @@ AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
   AliAnaCaloTrackCorrBaseClass(),  
   fCalorimeter(""), 
   fM02MaxCut(0),    fM02MinCut(0),       
-  fMinNCells(0),  
-  fMassEtaMin(0),   fMassEtaMax(0),
-  fMassPi0Min(0),   fMassPi0Max(0),
-  fMassConMin(0),   fMassConMax(0),
-  fPlotCluster(0)
+  fMinNCells(0),    fMinBadDist(0),
+  fFillAngleHisto(kFALSE),
+  fFillTMResidualHisto(kFALSE),
+  fFillSSExtraHisto(kFALSE),
+  fFillMCFractionHisto(kFALSE),
+  fhMassM02CutNLocMax1(0),    fhMassM02CutNLocMax2(0),    fhMassM02CutNLocMaxN(0),
+  fhMassSplitECutNLocMax1(0), fhMassSplitECutNLocMax2(0), fhMassSplitECutNLocMaxN(0),
+  fhMassAsyCutNLocMax1(0),    fhMassAsyCutNLocMax2(0),    fhMassAsyCutNLocMaxN(0)
 {
   //default ctor
   
@@ -70,14 +70,10 @@ AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
   {
     for(Int_t j = 0; j < 2; j++)
     {
-      
       fhMassNLocMax1[i][j]  = 0;
       fhMassNLocMax2[i][j]  = 0;
       fhMassNLocMaxN[i][j]  = 0;
       fhNLocMax[i][j]       = 0;
-      fhNLocMaxNLabel[i][j] = 0;
-      fhNLocMaxEMax[i][j]   = 0;
-      fhNLocMaxEFrac[i][j]  = 0;
       fhNLocMaxM02Cut[i][j] = 0;
       fhM02NLocMax1[i][j]   = 0;
       fhM02NLocMax2[i][j]   = 0;
@@ -94,11 +90,81 @@ AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
       fhM02Pi0LocMaxN[i][j] = 0;
       fhM02EtaLocMaxN[i][j] = 0;
       fhM02ConLocMaxN[i][j] = 0;
+
+      fhMassPi0LocMax1[i][j] = 0;
+      fhMassEtaLocMax1[i][j] = 0;
+      fhMassConLocMax1[i][j] = 0;
+      fhMassPi0LocMax2[i][j] = 0;
+      fhMassEtaLocMax2[i][j] = 0;
+      fhMassConLocMax2[i][j] = 0;
+      fhMassPi0LocMaxN[i][j] = 0;
+      fhMassEtaLocMaxN[i][j] = 0;
+      fhMassConLocMaxN[i][j] = 0;
+      
+      fhAsyPi0LocMax1[i][j] = 0;
+      fhAsyEtaLocMax1[i][j] = 0;
+      fhAsyConLocMax1[i][j] = 0;
+      fhAsyPi0LocMax2[i][j] = 0;
+      fhAsyEtaLocMax2[i][j] = 0;
+      fhAsyConLocMax2[i][j] = 0;
+      fhAsyPi0LocMaxN[i][j] = 0;
+      fhAsyEtaLocMaxN[i][j] = 0;
+      fhAsyConLocMaxN[i][j] = 0;      
       
       fhMassM02NLocMax1[i][j]= 0;
       fhMassM02NLocMax2[i][j]= 0;
-      fhMassM02NLocMaxN[i][j]= 0;      
+      fhMassM02NLocMaxN[i][j]= 0;   
+      fhMassDispEtaNLocMax1[i][j]= 0;
+      fhMassDispEtaNLocMax2[i][j]= 0;
+      fhMassDispEtaNLocMaxN[i][j]= 0;      
+      fhMassDispPhiNLocMax1[i][j]= 0;
+      fhMassDispPhiNLocMax2[i][j]= 0;
+      fhMassDispPhiNLocMaxN[i][j]= 0;      
+      fhMassDispAsyNLocMax1[i][j]= 0;
+      fhMassDispAsyNLocMax2[i][j]= 0;
+      fhMassDispAsyNLocMaxN[i][j]= 0;      
+      
+      fhSplitEFractionNLocMax1[i][j]=0;
+      fhSplitEFractionNLocMax2[i][j]=0;
+      fhSplitEFractionNLocMaxN[i][j]=0;
+      
+      fhMCGenFracNLocMax1[i][j]= 0;
+      fhMCGenFracNLocMax2[i][j]= 0;
+      fhMCGenFracNLocMaxN[i][j]= 0;
+      
+      fhMCGenSplitEFracNLocMax1[i][j]= 0;
+      fhMCGenSplitEFracNLocMax2[i][j]= 0;
+      fhMCGenSplitEFracNLocMaxN[i][j]= 0;    
+      
+      fhMCGenEFracvsSplitEFracNLocMax1[i][j]= 0;
+      fhMCGenEFracvsSplitEFracNLocMax2[i][j]= 0;
+      fhMCGenEFracvsSplitEFracNLocMaxN[i][j]= 0;    
       
+      fhMCGenEvsSplitENLocMax1[i][j]= 0;
+      fhMCGenEvsSplitENLocMax2[i][j]= 0;
+      fhMCGenEvsSplitENLocMaxN[i][j]= 0;     
+      
+      fhAsymNLocMax1 [i][j] = 0;
+      fhAsymNLocMax2 [i][j] = 0;
+      fhAsymNLocMaxN [i][j] = 0;
+    }
+    
+    for(Int_t jj = 0; jj < 4; jj++)
+    {
+      fhM02MCGenFracNLocMax1Ebin[i][jj] = 0;
+      fhM02MCGenFracNLocMax2Ebin[i][jj] = 0;
+      fhM02MCGenFracNLocMaxNEbin[i][jj] = 0;
+      
+      fhMassMCGenFracNLocMax1Ebin[i][jj]= 0;
+      fhMassMCGenFracNLocMax2Ebin[i][jj]= 0;
+      fhMassMCGenFracNLocMaxNEbin[i][jj]= 0;
+      
+      fhMCGenFracNLocMaxEbin[i][jj]       = 0;
+      fhMCGenFracNLocMaxEbinMatched[i][jj]= 0;
+      
+      fhMassSplitEFractionNLocMax1Ebin[i][jj] = 0;
+      fhMassSplitEFractionNLocMax2Ebin[i][jj] = 0;
+      fhMassSplitEFractionNLocMaxNEbin[i][jj] = 0;
     }
     
     fhTrackMatchedDEtaLocMax1[i] = 0; 
@@ -118,6 +184,9 @@ AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
     fhAnglePairMassLocMax1[i] = 0;
     fhAnglePairMassLocMax2[i] = 0;
     fhAnglePairMassLocMaxN[i] = 0;
+    fhSplitEFractionvsAsyNLocMax1[i] = 0;
+    fhSplitEFractionvsAsyNLocMax2[i] = 0; 
+    fhSplitEFractionvsAsyNLocMaxN[i] = 0;    
   }
   
   for(Int_t i = 0; i < 4; i++)
@@ -125,6 +194,22 @@ AliAnaInsideClusterInvariantMass::AliAnaInsideClusterInvariantMass() :
     fhMassM02NLocMax1Ebin[i] = 0 ;
     fhMassM02NLocMax2Ebin[i] = 0 ;
     fhMassM02NLocMaxNEbin[i] = 0 ;
+    
+    fhMassDispEtaNLocMax1Ebin[i] = 0 ;
+    fhMassDispEtaNLocMax2Ebin[i] = 0 ;
+    fhMassDispEtaNLocMaxNEbin[i] = 0 ;
+    
+    fhMassDispPhiNLocMax1Ebin[i] = 0 ;
+    fhMassDispPhiNLocMax2Ebin[i] = 0 ;
+    fhMassDispPhiNLocMaxNEbin[i] = 0 ;    
+    
+    fhMassDispAsyNLocMax1Ebin[i] = 0 ;
+    fhMassDispAsyNLocMax2Ebin[i] = 0 ;
+    fhMassDispAsyNLocMaxNEbin[i] = 0 ;    
+
+    fhMCAsymM02NLocMax1MCPi0Ebin[i] = 0 ;
+    fhMCAsymM02NLocMax2MCPi0Ebin[i] = 0 ;
+    fhMCAsymM02NLocMaxNMCPi0Ebin[i] = 0 ;
   }
   
   InitParameters();
@@ -151,52 +236,14 @@ TObjString *  AliAnaInsideClusterInvariantMass::GetAnalysisCuts()
   snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n",    fM02MinCut, fM02MaxCut) ;
   parList+=onePar ;
   snprintf(onePar,buffersize,"fMinNCells =%d \n",        fMinNCells) ;
+  parList+=onePar ;    
+  snprintf(onePar,buffersize,"fMinBadDist =%1.1f \n",    fMinBadDist) ;
   parList+=onePar ;  
-  snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max);
-  parList+=onePar ;
-  snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax);
-  parList+=onePar ;
-  snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassConMin,fMassConMax);
-  parList+=onePar ;
 
   return new TObjString(parList) ;
   
 }
 
-
-//_____________________________________________________________________________________
-TLorentzVector AliAnaInsideClusterInvariantMass::GetCellMomentum(const Int_t absId,
-                                                                 Float_t en,
-                                                                 AliVCaloCells * cells)
-{
-
-  // Cell momentum calculation for invariant mass
-  
-  Double_t cellpos[] = {0, 0, 0};
-  GetEMCALGeometry()->GetGlobal(absId, cellpos);
-  
-  if(GetVertex(0)){//calculate direction from vertex
-    cellpos[0]-=GetVertex(0)[0];
-    cellpos[1]-=GetVertex(0)[1];
-    cellpos[2]-=GetVertex(0)[2];  
-  }
-  
-  Double_t r = TMath::Sqrt(cellpos[0]*cellpos[0]+cellpos[1]*cellpos[1]+cellpos[2]*cellpos[2] ) ;
-  
-  //If not calculated before, get the energy
-  if(en <=0 )
-  {
-    en = cells->GetCellAmplitude(absId);
-    GetCaloUtils()->RecalibrateCellAmplitude(en,fCalorimeter,absId);  
-  }
-  
-  TLorentzVector cellMom ;   
-  cellMom.SetPxPyPzE( en*cellpos[0]/r,  en*cellpos[1]/r, en*cellpos[2]/r,  en) ;   
-
-  return cellMom;
-  
-}
-
 //________________________________________________________________
 TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
 {  
@@ -228,86 +275,205 @@ TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
   
   TString sMatched[] = {"","Matched"};
   
+  
+  fhMassSplitECutNLocMax1  = new TH2F("hMassSplitECutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, (E1+E2)/E cut",
+                                   nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+  fhMassSplitECutNLocMax1->SetYTitle("M (GeV/c^{2})");
+  fhMassSplitECutNLocMax1->SetXTitle("E (GeV)");
+  outputContainer->Add(fhMassSplitECutNLocMax1) ;   
+  
+  fhMassSplitECutNLocMax2  = new TH2F("hMassSplitECutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, (E1+E2)/E cut",
+                                   nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+  fhMassSplitECutNLocMax2->SetYTitle("M (GeV/c^{2})");
+  fhMassSplitECutNLocMax2->SetXTitle("E (GeV)");
+  outputContainer->Add(fhMassSplitECutNLocMax2) ;   
+  
+  fhMassSplitECutNLocMaxN  = new TH2F("hMassSplitECutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, (E1+E2)/E cut",
+                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+  fhMassSplitECutNLocMaxN->SetYTitle("M (GeV/c^{2})");
+  fhMassSplitECutNLocMaxN->SetXTitle("E (GeV)");
+  outputContainer->Add(fhMassSplitECutNLocMaxN) ;   
+
+  fhMassM02CutNLocMax1  = new TH2F("hMassM02CutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, M02 cut",
+                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+  fhMassM02CutNLocMax1->SetYTitle("M (GeV/c^{2})");
+  fhMassM02CutNLocMax1->SetXTitle("E (GeV)");
+  outputContainer->Add(fhMassM02CutNLocMax1) ;   
+  
+  fhMassM02CutNLocMax2  = new TH2F("hMassM02CutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, M02 cut",
+                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+  fhMassM02CutNLocMax2->SetYTitle("M (GeV/c^{2})");
+  fhMassM02CutNLocMax2->SetXTitle("E (GeV)");
+  outputContainer->Add(fhMassM02CutNLocMax2) ;   
+  
+  fhMassM02CutNLocMaxN  = new TH2F("hMassM02CutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, M02 cut",
+                             nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+  fhMassM02CutNLocMaxN->SetYTitle("M (GeV/c^{2})");
+  fhMassM02CutNLocMaxN->SetXTitle("E (GeV)");
+  outputContainer->Add(fhMassM02CutNLocMaxN) ;   
+  
+  fhMassAsyCutNLocMax1  = new TH2F("hMassAsyCutNLocMax1","Invariant mass of splitted cluster with NLM=1 vs E, with |A|>0.8",
+                                   nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+  fhMassAsyCutNLocMax1->SetYTitle("M (GeV/c^{2})");
+  fhMassAsyCutNLocMax1->SetXTitle("E (GeV)");
+  outputContainer->Add(fhMassAsyCutNLocMax1) ;   
+  
+  fhMassAsyCutNLocMax2  = new TH2F("hMassAsyCutNLocMax2","Invariant mass of splitted cluster with NLM=2 vs E, with |A|>0.8",
+                                   nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+  fhMassAsyCutNLocMax2->SetYTitle("M (GeV/c^{2})");
+  fhMassAsyCutNLocMax2->SetXTitle("E (GeV)");
+  outputContainer->Add(fhMassAsyCutNLocMax2) ;   
+  
+  fhMassAsyCutNLocMaxN  = new TH2F("hMassAsyCutNLocMaxN","Invariant mass of splitted cluster with NLM>2 vs E, with |A|>0.8",
+                                   nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+  fhMassAsyCutNLocMaxN->SetYTitle("M (GeV/c^{2})");
+  fhMassAsyCutNLocMaxN->SetXTitle("E (GeV)");
+  outputContainer->Add(fhMassAsyCutNLocMaxN) ;   
+  
+  
   for(Int_t i = 0; i < n; i++)
   {  
-    
     for(Int_t j = 0; j < 2; j++)
     {  
       
       fhMassNLocMax1[i][j]  = new TH2F(Form("hMassNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
-                                       Form("Invariant mass of 2 highest energy cells vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                       Form("Invariant mass of splitted cluster with NLM=1 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
       fhMassNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
       fhMassNLocMax1[i][j]->SetXTitle("E (GeV)");
       outputContainer->Add(fhMassNLocMax1[i][j]) ;   
       
       fhMassNLocMax2[i][j]  = new TH2F(Form("hMassNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
-                                       Form("Invariant mass of 2 local maxima cells vs E,%s %s",ptype[i].Data(),sMatched[j].Data()),
+                                       Form("Invariant mass of splitted cluster with NLM=2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
       fhMassNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
       fhMassNLocMax2[i][j]->SetXTitle("E (GeV)");
       outputContainer->Add(fhMassNLocMax2[i][j]) ;   
       
       fhMassNLocMaxN[i][j]  = new TH2F(Form("hMassNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
-                                       Form("Invariant mass of N>2 local maxima cells vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                       Form("Invariant mass of splitted cluster with NLM>2 vs E, %s %s",ptype[i].Data(),sMatched[j].Data()),
                                        nptbins,ptmin,ptmax,mbins,mmin,mmax); 
       fhMassNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
       fhMassNLocMaxN[i][j]->SetXTitle("E (GeV)");
       outputContainer->Add(fhMassNLocMaxN[i][j]) ;   
       
       fhMassM02NLocMax1[i][j]  = new TH2F(Form("hMassM02NLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
-                                          Form("Invariant mass of 2 highest energy cells #lambda_{0}^{2}, E > 7 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
+                                          Form("Invariant mass of splitted cluster with NLM=1, #lambda_{0}^{2}, E > 8 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
       fhMassM02NLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
       fhMassM02NLocMax1[i][j]->SetXTitle("#lambda_{0}^{2}");
       outputContainer->Add(fhMassM02NLocMax1[i][j]) ;   
       
       fhMassM02NLocMax2[i][j]  = new TH2F(Form("hMassM02NLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
-                                          Form("Invariant mass of 2 local maxima cells #lambda_{0}^{2}, E > 7 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                          Form("Invariant mass of splitted cluster with NLM=2, #lambda_{0}^{2}, E > 8 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
       fhMassM02NLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
       fhMassM02NLocMax2[i][j]->SetXTitle("#lambda_{0}^{2}");
       outputContainer->Add(fhMassM02NLocMax2[i][j]) ;   
       
       fhMassM02NLocMaxN[i][j]  = new TH2F(Form("hMassM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
-                                          Form("Invariant mass of N>2 local maxima cells vs #lambda_{0}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                          Form("Invariant mass of splitted cluster with NLM>2, vs #lambda_{0}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
                                           ssbins,ssmin,ssmax,mbins,mmin,mmax); 
       fhMassM02NLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
       fhMassM02NLocMaxN[i][j]->SetXTitle("#lambda_{0}^{2}");
       outputContainer->Add(fhMassM02NLocMaxN[i][j]) ;   
       
       
+      fhAsymNLocMax1[i][j]  = new TH2F(Form("hAsymNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                    Form("Asymmetry of NLM=1  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                    nptbins,ptmin,ptmax,200,-1,1); 
+      fhAsymNLocMax1[i][j]->SetYTitle("#alpha (rad)");
+      fhAsymNLocMax1[i][j]->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsymNLocMax1[i][j]) ;   
+      
+      fhAsymNLocMax2[i][j]  = new TH2F(Form("hAsymNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                    Form("Asymmetry of NLM=2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                    nptbins,ptmin,ptmax,200,-1,1); 
+      fhAsymNLocMax2[i][j]->SetYTitle("#alpha (rad)");
+      fhAsymNLocMax2[i][j]->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsymNLocMax2[i][j]) ;   
+      
+      fhAsymNLocMaxN[i][j]  = new TH2F(Form("hAsymNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                    Form("Asymmetry of NLM>2  vs cluster Energy, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                    nptbins,ptmin,ptmax,200,-1,1); 
+      fhAsymNLocMaxN[i][j]->SetYTitle("(E_{1}-E_{2})/(E_{1}+E_{2})");
+      fhAsymNLocMaxN[i][j]->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsymNLocMaxN[i][j]) ;   
+      
+      
+      if(fFillSSExtraHisto)
+      {
+        fhMassDispEtaNLocMax1[i][j]  = new TH2F(Form("hMassDispEtaNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("Invariant mass of splitted cluster with NLM=1, #sigma_{#eta #eta}^{2}, E > 8 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+        fhMassDispEtaNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassDispEtaNLocMax1[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
+        outputContainer->Add(fhMassDispEtaNLocMax1[i][j]) ;   
+        
+        fhMassDispEtaNLocMax2[i][j]  = new TH2F(Form("hMassDispEtaNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("Invariant mass of splitted cluster with NLM=2 #sigma_{#eta #eta}^{2}, E > 8 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+        fhMassDispEtaNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassDispEtaNLocMax2[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
+        outputContainer->Add(fhMassDispEtaNLocMax2[i][j]) ;   
+        
+        fhMassDispEtaNLocMaxN[i][j]  = new TH2F(Form("hMassDispEtaNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("Invariant mass of splitted cluster with NLM>2, #sigma_{#eta #eta}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+        fhMassDispEtaNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassDispEtaNLocMaxN[i][j]->SetXTitle("#sigma_{#eta #eta}^{2}");
+        outputContainer->Add(fhMassDispEtaNLocMaxN[i][j]) ;   
+        
+        fhMassDispPhiNLocMax1[i][j]  = new TH2F(Form("hMassDispPhiNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E > 8 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+        fhMassDispPhiNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassDispPhiNLocMax1[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
+        outputContainer->Add(fhMassDispPhiNLocMax1[i][j]) ;   
+        
+        fhMassDispPhiNLocMax2[i][j]  = new TH2F(Form("hMassDispPhiNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E > 8 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+        fhMassDispPhiNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassDispPhiNLocMax2[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
+        outputContainer->Add(fhMassDispPhiNLocMax2[i][j]) ;   
+        
+        fhMassDispPhiNLocMaxN[i][j]  = new TH2F(Form("hMassDispPhiNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+        fhMassDispPhiNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassDispPhiNLocMaxN[i][j]->SetXTitle("#sigma_{#phi #phi}^{2}");
+        outputContainer->Add(fhMassDispPhiNLocMaxN[i][j]) ;   
+        
+        fhMassDispAsyNLocMax1[i][j]  = new TH2F(Form("hMassDispAsyNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("Invariant mass of 2 highest energy cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E > 8 GeV,%s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                200,-1,1,mbins,mmin,mmax); 
+        fhMassDispAsyNLocMax1[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassDispAsyNLocMax1[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
+        outputContainer->Add(fhMassDispAsyNLocMax1[i][j]) ;   
+        
+        fhMassDispAsyNLocMax2[i][j]  = new TH2F(Form("hMassDispAsyNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("Invariant mass of 2 local maxima cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E > 8 GeV, %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                200,-1,1,mbins,mmin,mmax); 
+        fhMassDispAsyNLocMax2[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassDispAsyNLocMax2[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
+        outputContainer->Add(fhMassDispAsyNLocMax2[i][j]) ;   
+        
+        fhMassDispAsyNLocMaxN[i][j]  = new TH2F(Form("hMassDispAsyNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("Invariant mass of N>2 local maxima cells vsA = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                200,-1,1,mbins,mmin,mmax); 
+        fhMassDispAsyNLocMaxN[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassDispAsyNLocMaxN[i][j]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
+        outputContainer->Add(fhMassDispAsyNLocMaxN[i][j]) ;   
+      }
+      
       fhNLocMax[i][j]     = new TH2F(Form("hNLocMax%s%s",pname[i].Data(),sMatched[j].Data()),
                                      Form("Number of local maxima in cluster %s %s",ptype[i].Data(),sMatched[j].Data()),
                                      nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins); 
       fhNLocMax[i][j]   ->SetYTitle("N maxima");
       fhNLocMax[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhNLocMax[i][j]) ; 
-      
-      if(IsDataMC())
-      {
-        fhNLocMaxNLabel[i][j]     = new TH2F(Form("hNLocMaxNLabel%s%s",pname[i].Data(),sMatched[j].Data()),
-                                             Form("Number of local maxima in cluster vs number of MC labels %s %s",ptype[i].Data(),sMatched[j].Data()),
-                                             nMaxBins,0,nMaxBins,nMaxBins,0,nMaxBins); 
-        fhNLocMaxNLabel[i][j]   ->SetYTitle("N maxima");
-        fhNLocMaxNLabel[i][j]   ->SetXTitle("N MC labels");
-        outputContainer->Add(fhNLocMaxNLabel[i][j]) ; 
-      }
-      
-      fhNLocMaxEMax[i][j]     = new TH2F(Form("hNLocMaxEMax%s%s",pname[i].Data(),sMatched[j].Data()),
-                                         Form("Number of local maxima in cluster vs energy of maxima %s %s",ptype[i].Data(),sMatched[j].Data()),
-                                         nptbins*10,ptmin,ptmax,nMaxBins,0,nMaxBins); 
-      fhNLocMaxEMax[i][j]   ->SetYTitle("N maxima");
-      fhNLocMaxEMax[i][j]   ->SetXTitle("E of maxima (GeV)");
-      outputContainer->Add(fhNLocMaxEMax[i][j]) ; 
-      
-      fhNLocMaxEFrac[i][j]     = new TH2F(Form("hNLocMaxEFrac%s%s",pname[i].Data(),sMatched[j].Data()),
-                                          Form("Number of local maxima in cluster vs fraction of cluster energy of maxima %s %s",ptype[i].Data(),sMatched[j].Data()),
-                                          100,0,1,nMaxBins,0,nMaxBins); 
-      fhNLocMaxEFrac[i][j]   ->SetYTitle("N maxima");
-      fhNLocMaxEFrac[i][j]   ->SetXTitle("E maxima / E cluster");
-      outputContainer->Add(fhNLocMaxEFrac[i][j]) ; 
-      
+            
       fhNLocMaxM02Cut[i][j] = new TH2F(Form("hNLocMaxM02Cut%s%s",pname[i].Data(),sMatched[j].Data()),
                                        Form("Number of local maxima in cluster %s for %2.2f < M02 < %2.2f",ptype[i].Data(),fM02MinCut,fM02MaxCut),
                                        nptbins,ptmin,ptmax,nMaxBins,0,nMaxBins); 
@@ -330,7 +496,6 @@ TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
       fhM02NLocMax2[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhM02NLocMax2[i][j]) ; 
       
-      
       fhM02NLocMaxN[i][j]    = new TH2F(Form("hM02NLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
                                         Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
@@ -339,99 +504,477 @@ TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
       outputContainer->Add(fhM02NLocMaxN[i][j]) ; 
       
       
-      fhNCellNLocMax1[i][j]  = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
-                                        Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
-                                        nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
-      fhNCellNLocMax1[i][j] ->SetYTitle("N cells");
-      fhNCellNLocMax1[i][j] ->SetXTitle("E (GeV)");
-      outputContainer->Add(fhNCellNLocMax1[i][j]) ; 
+      fhSplitEFractionNLocMax1[i][j]     = new TH2F(Form("hSplitEFractionNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                         Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                         nptbins,ptmin,ptmax,120,0,1.2); 
+      fhSplitEFractionNLocMax1[i][j]   ->SetXTitle("E_{cluster} (GeV)");
+      fhSplitEFractionNLocMax1[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
+      outputContainer->Add(fhSplitEFractionNLocMax1[i][j]) ; 
+      
+      fhSplitEFractionNLocMax2[i][j]     = new TH2F(Form("hSplitEFractionNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                         Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                         nptbins,ptmin,ptmax,120,0,1.2); 
+      fhSplitEFractionNLocMax2[i][j]   ->SetXTitle("E_{cluster} (GeV)");
+      fhSplitEFractionNLocMax2[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
+      outputContainer->Add(fhSplitEFractionNLocMax2[i][j]) ; 
+      
+      fhSplitEFractionNLocMaxN[i][j]    = new TH2F(Form("hSplitEFractionNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                        Form("(E1+E2)/E_{cluster} vs E_{cluster} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                        nptbins,ptmin,ptmax,120,0,1.2); 
+      fhSplitEFractionNLocMaxN[i][j]   ->SetXTitle("E_{cluster} (GeV)");
+      fhSplitEFractionNLocMaxN[i][j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
+      outputContainer->Add(fhSplitEFractionNLocMaxN[i][j]) ; 
       
-      fhNCellNLocMax2[i][j]     = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
-                                           Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
-                                           nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
-      fhNCellNLocMax2[i][j]   ->SetYTitle("N cells");
-      fhNCellNLocMax2[i][j]   ->SetXTitle("E (GeV)");
-      outputContainer->Add(fhNCellNLocMax2[i][j]) ; 
       
+      if(i > 0 && fFillMCFractionHisto) // skip first entry in array, general case not filled
+      {
+        fhMCGenFracNLocMax1[i][j]     = new TH2F(Form("hMCGenFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                 Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                 nptbins,ptmin,ptmax,200,0,2); 
+        fhMCGenFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
+        fhMCGenFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
+        outputContainer->Add(fhMCGenFracNLocMax1[i][j]) ; 
+        
+        fhMCGenFracNLocMax2[i][j]     = new TH2F(Form("hMCGenFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                 Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                 nptbins,ptmin,ptmax,200,0,2); 
+        fhMCGenFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
+        fhMCGenFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
+        outputContainer->Add(fhMCGenFracNLocMax2[i][j]) ; 
+        
+        
+        fhMCGenFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                nptbins,ptmin,ptmax,200,0,2); 
+        fhMCGenFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / E_{reco}");
+        fhMCGenFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
+        outputContainer->Add(fhMCGenFracNLocMaxN[i][j]) ; 
       
-      fhNCellNLocMaxN[i][j]     = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
-                                           Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
-                                           nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
-      fhNCellNLocMaxN[i][j]   ->SetYTitle("N cells");
-      fhNCellNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
-      outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
+        fhMCGenSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                 Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                 nptbins,ptmin,ptmax,200,0,2); 
+        fhMCGenSplitEFracNLocMax1[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
+        fhMCGenSplitEFracNLocMax1[i][j]   ->SetXTitle("E (GeV)");
+        outputContainer->Add(fhMCGenSplitEFracNLocMax1[i][j]) ; 
+        
+        fhMCGenSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                 Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                 nptbins,ptmin,ptmax,200,0,2); 
+        fhMCGenSplitEFracNLocMax2[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
+        fhMCGenSplitEFracNLocMax2[i][j]   ->SetXTitle("E (GeV)");
+        outputContainer->Add(fhMCGenSplitEFracNLocMax2[i][j]) ; 
+        
+        
+        fhMCGenSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                Form("E_{gen} / (E_{1 split}+E_{2 split}) vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                nptbins,ptmin,ptmax,200,0,2); 
+        fhMCGenSplitEFracNLocMaxN[i][j]   ->SetYTitle("E_{gen} / (E_{1 split}+E_{2 split})");
+        fhMCGenSplitEFracNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
+        outputContainer->Add(fhMCGenSplitEFracNLocMaxN[i][j]) ; 
+       
+        fhMCGenEFracvsSplitEFracNLocMax1[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                       Form("(E_{1 split}+E_{2 split})/E_{reco} vs E_{gen} / E_{reco} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                       200,0,2,200,0,2); 
+        fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
+        fhMCGenEFracvsSplitEFracNLocMax1[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
+        outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax1[i][j]) ; 
+        
+        fhMCGenEFracvsSplitEFracNLocMax2[i][j]     = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                       Form("(E_{1 split}+E_{2 split})/E_{reco} vs E_{gen} / E_{reco} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                       200,0,2,200,0,2); 
+        fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
+        fhMCGenEFracvsSplitEFracNLocMax2[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
+        outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMax2[i][j]) ; 
+        
+        
+        fhMCGenEFracvsSplitEFracNLocMaxN[i][j]    = new TH2F(Form("hMCGenEFracvsSplitEFracNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                      Form("(E_{1 split}+E_{2 split})/E_{reco} vs E_{gen} / E_{reco} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                      200,0,2,200,0,2); 
+        fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetYTitle("(E_{1 split}+E_{2 split})/E_{reco}");
+        fhMCGenEFracvsSplitEFracNLocMaxN[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
+        outputContainer->Add(fhMCGenEFracvsSplitEFracNLocMaxN[i][j]) ; 
+        
+        
+        fhMCGenEvsSplitENLocMax1[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                             Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
+        fhMCGenEvsSplitENLocMax1[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
+        fhMCGenEvsSplitENLocMax1[i][j]   ->SetXTitle("E_{gen} (GeV)");
+        outputContainer->Add(fhMCGenEvsSplitENLocMax1[i][j]) ; 
+        
+        fhMCGenEvsSplitENLocMax2[i][j]     = new TH2F(Form("hMCGenEvsSplitENLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                             Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
+        fhMCGenEvsSplitENLocMax2[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
+        fhMCGenEvsSplitENLocMax2[i][j]   ->SetXTitle("E_{gen} (GeV)");
+        outputContainer->Add(fhMCGenEvsSplitENLocMax2[i][j]) ; 
+        
+        
+        fhMCGenEvsSplitENLocMaxN[i][j]    = new TH2F(Form("hMCGenEvsSplitENLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                                            Form("E_{1 split}+E_{2 split} vs E_{gen} for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); 
+        fhMCGenEvsSplitENLocMaxN[i][j]   ->SetYTitle("E_{1 split}+E_{2 split} (GeV)");
+        fhMCGenEvsSplitENLocMaxN[i][j]   ->SetXTitle("E_{gen} (GeV)");
+        outputContainer->Add(fhMCGenEvsSplitENLocMaxN[i][j]) ; 
+        
+      }
       
+      if(fFillSSExtraHisto)
+      {
+        fhNCellNLocMax1[i][j]  = new TH2F(Form("hNCellNLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                          Form("#lambda_{0}^{2} vs E for N max  = 1 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                          nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
+        fhNCellNLocMax1[i][j] ->SetYTitle("N cells");
+        fhNCellNLocMax1[i][j] ->SetXTitle("E (GeV)");
+        outputContainer->Add(fhNCellNLocMax1[i][j]) ; 
+        
+        fhNCellNLocMax2[i][j]     = new TH2F(Form("hNCellNLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                             Form("#lambda_{0}^{2} vs E for N max  = 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                             nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
+        fhNCellNLocMax2[i][j]   ->SetYTitle("N cells");
+        fhNCellNLocMax2[i][j]   ->SetXTitle("E (GeV)");
+        outputContainer->Add(fhNCellNLocMax2[i][j]) ; 
+        
+        
+        fhNCellNLocMaxN[i][j]     = new TH2F(Form("hNCellNLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                             Form("#lambda_{0}^{2} vs E for N max  > 2 %s %s",ptype[i].Data(),sMatched[j].Data()),
+                                             nptbins,ptmin,ptmax,ncbins,ncmin,ncmax); 
+        fhNCellNLocMaxN[i][j]   ->SetYTitle("N cells");
+        fhNCellNLocMaxN[i][j]   ->SetXTitle("E (GeV)");
+        outputContainer->Add(fhNCellNLocMaxN[i][j]) ;
+      }
       
       fhM02Pi0LocMax1[i][j]     = new TH2F(Form("hM02Pi0LocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
-                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 1",fMassPi0Min,fMassPi0Max,ptype[i].Data()),
+                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 1",
+                                                GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
       fhM02Pi0LocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
       fhM02Pi0LocMax1[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhM02Pi0LocMax1[i][j]) ; 
       
       fhM02EtaLocMax1[i][j]     = new TH2F(Form("hM02EtaLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
-                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",fMassEtaMin,fMassEtaMax,ptype[i].Data()),
+                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
+                                                GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
       fhM02EtaLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
       fhM02EtaLocMax1[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhM02EtaLocMax1[i][j]) ; 
       
       fhM02ConLocMax1[i][j]    = new TH2F(Form("hM02ConLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
-                                          Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",fMassConMin,fMassConMax,ptype[i].Data()),
+                                          Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
+                                               GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
       fhM02ConLocMax1[i][j]   ->SetYTitle("#lambda_{0}^{2}");
       fhM02ConLocMax1[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhM02ConLocMax1[i][j]) ; 
       
       fhM02Pi0LocMax2[i][j]     = new TH2F(Form("hM02Pi0LocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
-                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 2",fMassPi0Min,fMassPi0Max,ptype[i].Data()),
+                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 2",
+                                                GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
       fhM02Pi0LocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
       fhM02Pi0LocMax2[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhM02Pi0LocMax2[i][j]) ; 
       
       fhM02EtaLocMax2[i][j]     = new TH2F(Form("hM02EtaLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
-                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",fMassEtaMin,fMassEtaMax,ptype[i].Data()),
+                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
+                                               GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
       fhM02EtaLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
       fhM02EtaLocMax2[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhM02EtaLocMax2[i][j]) ; 
       
       fhM02ConLocMax2[i][j]    = new TH2F(Form("hM02ConLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
-                                          Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",fMassConMin,fMassConMax,ptype[i].Data()),
+                                          Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
+                                               GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
       fhM02ConLocMax2[i][j]   ->SetYTitle("#lambda_{0}^{2}");
       fhM02ConLocMax2[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhM02ConLocMax2[i][j]) ; 
       
       fhM02Pi0LocMaxN[i][j]     = new TH2F(Form("hM02Pi0LocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
-                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max > 2",fMassPi0Min,fMassPi0Max,ptype[i].Data()),
+                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max > 2",
+                                                GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
       fhM02Pi0LocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
       fhM02Pi0LocMaxN[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhM02Pi0LocMaxN[i][j]) ; 
       
       fhM02EtaLocMaxN[i][j]     = new TH2F(Form("hM02EtaLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
-                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max > 2", fMassEtaMin,fMassEtaMax,ptype[i].Data()),
+                                           Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max > 2", 
+                                                GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
       fhM02EtaLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
       fhM02EtaLocMaxN[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhM02EtaLocMaxN[i][j]) ; 
       
       fhM02ConLocMaxN[i][j]    = new TH2F(Form("hM02ConLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
-                                          Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",fMassConMin,fMassConMax,ptype[i].Data()),
+                                          Form("#lambda_{0}^{2} vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
+                                               GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
       fhM02ConLocMaxN[i][j]   ->SetYTitle("#lambda_{0}^{2}");
       fhM02ConLocMaxN[i][j]   ->SetXTitle("E (GeV)");
       outputContainer->Add(fhM02ConLocMaxN[i][j]) ; 
+            
+      
+      fhMassPi0LocMax1[i][j]     = new TH2F(Form("hMassPi0LocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                           Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 1",
+                                                GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
+                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhMassPi0LocMax1[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
+      fhMassPi0LocMax1[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhMassPi0LocMax1[i][j]) ; 
+      
+      fhMassEtaLocMax1[i][j]     = new TH2F(Form("hMassEtaLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                           Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
+                                                GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
+                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhMassEtaLocMax1[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
+      fhMassEtaLocMax1[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhMassEtaLocMax1[i][j]) ; 
+      
+      fhMassConLocMax1[i][j]    = new TH2F(Form("hMassConLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                          Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
+                                               GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
+                                          nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhMassConLocMax1[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
+      fhMassConLocMax1[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhMassConLocMax1[i][j]) ; 
+      
+      fhMassPi0LocMax2[i][j]     = new TH2F(Form("hMassPi0LocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                           Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 2",
+                                                GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
+                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhMassPi0LocMax2[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
+      fhMassPi0LocMax2[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhMassPi0LocMax2[i][j]) ; 
+      
+      fhMassEtaLocMax2[i][j]     = new TH2F(Form("hMassEtaLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                           Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
+                                                GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
+                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhMassEtaLocMax2[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
+      fhMassEtaLocMax2[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhMassEtaLocMax2[i][j]) ; 
+      
+      fhMassConLocMax2[i][j]    = new TH2F(Form("hMassConLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                          Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
+                                               GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
+                                          nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhMassConLocMax2[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
+      fhMassConLocMax2[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhMassConLocMax2[i][j]) ; 
+      
+      fhMassPi0LocMaxN[i][j]     = new TH2F(Form("hMassPi0LocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                           Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max > 2",
+                                                GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
+                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhMassPi0LocMaxN[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
+      fhMassPi0LocMaxN[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhMassPi0LocMaxN[i][j]) ; 
+      
+      fhMassEtaLocMaxN[i][j]     = new TH2F(Form("hMassEtaLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                           Form("Mass vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max > 2", 
+                                                GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
+                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhMassEtaLocMaxN[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
+      fhMassEtaLocMaxN[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhMassEtaLocMaxN[i][j]) ; 
+      
+      fhMassConLocMaxN[i][j]    = new TH2F(Form("hMassConLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                          Form("Mass vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
+                                               GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
+                                          nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhMassConLocMaxN[i][j]   ->SetYTitle("Mass (MeV/c^{2})");
+      fhMassConLocMaxN[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhMassConLocMaxN[i][j]) ; 
+      
+      
+      fhAsyPi0LocMax1[i][j]     = new TH2F(Form("hAsyPi0LocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 1",
+                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
+                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhAsyPi0LocMax1[i][j]   ->SetYTitle("Asymmetry");
+      fhAsyPi0LocMax1[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsyPi0LocMax1[i][j]) ; 
+      
+      fhAsyEtaLocMax1[i][j]     = new TH2F(Form("hAsyEtaLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
+                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
+                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhAsyEtaLocMax1[i][j]   ->SetYTitle("Asymmetry");
+      fhAsyEtaLocMax1[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsyEtaLocMax1[i][j]) ; 
+      
+      fhAsyConLocMax1[i][j]    = new TH2F(Form("hAsyConLocMax1%s%s",pname[i].Data(),sMatched[j].Data()),
+                                           Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 1",
+                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
+                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhAsyConLocMax1[i][j]   ->SetYTitle("Asymmetry");
+      fhAsyConLocMax1[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsyConLocMax1[i][j]) ; 
+      
+      fhAsyPi0LocMax2[i][j]     = new TH2F(Form("hAsyPi0LocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max = 2",
+                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
+                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhAsyPi0LocMax2[i][j]   ->SetYTitle("Asymmetry");
+      fhAsyPi0LocMax2[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsyPi0LocMax2[i][j]) ; 
+      
+      fhAsyEtaLocMax2[i][j]     = new TH2F(Form("hAsyEtaLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
+                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
+                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhAsyEtaLocMax2[i][j]   ->SetYTitle("Asymmetry");
+      fhAsyEtaLocMax2[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsyEtaLocMax2[i][j]) ; 
+      
+      fhAsyConLocMax2[i][j]    = new TH2F(Form("hAsyConLocMax2%s%s",pname[i].Data(),sMatched[j].Data()),
+                                           Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max = 2",
+                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
+                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhAsyConLocMax2[i][j]   ->SetYTitle("Asymmetry");
+      fhAsyConLocMax2[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsyConLocMax2[i][j]) ; 
+      
+      fhAsyPi0LocMaxN[i][j]     = new TH2F(Form("hAsyPi0LocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2} %s, for N Local max > 2",
+                                                 GetCaloPID()->GetPi0MinMass(),GetCaloPID()->GetPi0MaxMass(),ptype[i].Data()),
+                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhAsyPi0LocMaxN[i][j]   ->SetYTitle("Asymmetry");
+      fhAsyPi0LocMaxN[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsyPi0LocMaxN[i][j]) ; 
+      
+      fhAsyEtaLocMaxN[i][j]     = new TH2F(Form("hAsyEtaLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                            Form("Asymmetry vs E for mass range [%2.2f-%2.2f] MeV/c^{2}, %s, for N Local max > 2", 
+                                                 GetCaloPID()->GetEtaMinMass(),GetCaloPID()->GetEtaMaxMass(),ptype[i].Data()),
+                                            nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhAsyEtaLocMaxN[i][j]   ->SetYTitle("Asymmetry");
+      fhAsyEtaLocMaxN[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsyEtaLocMaxN[i][j]) ; 
+      
+      fhAsyConLocMaxN[i][j]    = new TH2F(Form("hAsyConLocMaxN%s%s",pname[i].Data(),sMatched[j].Data()),
+                                           Form("Asymmetry vs E for mass range [%2.2f-%2.2f], %s, for N Local max > 2",
+                                                GetCaloPID()->GetPhotonMinMass(),GetCaloPID()->GetPhotonMaxMass(),ptype[i].Data()),
+                                           nptbins,ptmin,ptmax,mbins,mmin,mmax); 
+      fhAsyConLocMaxN[i][j]   ->SetYTitle("Asymmetry");
+      fhAsyConLocMaxN[i][j]   ->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAsyConLocMaxN[i][j]) ; 
       
     } // matched, not matched
     
-    
+      for(Int_t j = 0; j < 4; j++)
+      {  
+        
+        fhMassSplitEFractionNLocMax1Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax1%sEbin%d",pname[i].Data(),j),
+                                                           Form("Invariant mass of 2 highest energy cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
+                                                           120,0,1.2,mbins,mmin,mmax); 
+        fhMassSplitEFractionNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassSplitEFractionNLocMax1Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
+        outputContainer->Add(fhMassSplitEFractionNLocMax1Ebin[i][j]) ;   
+        
+        fhMassSplitEFractionNLocMax2Ebin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMax2%sEbin%d",pname[i].Data(),j),
+                                                           Form("Invariant mass of 2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
+                                                           120,0,1.2,mbins,mmin,mmax); 
+        fhMassSplitEFractionNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassSplitEFractionNLocMax2Ebin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
+        outputContainer->Add(fhMassSplitEFractionNLocMax2Ebin[i][j]) ;   
+        
+        fhMassSplitEFractionNLocMaxNEbin[i][j]  = new TH2F(Form("hMassSplitEFractionNLocMaxN%sEbin%d",pname[i].Data(),j),
+                                                           Form("Invariant mass of N>2 local maxima cells vs (E1+E2)/Ecluster, %s, E bin %d",ptype[i].Data(),j),
+                                                           120,0,1.2,mbins,mmin,mmax); 
+        fhMassSplitEFractionNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
+        fhMassSplitEFractionNLocMaxNEbin[i][j]->SetXTitle("(E_{split1}+E_{split2})/E_{cluster}");
+        outputContainer->Add(fhMassSplitEFractionNLocMaxNEbin[i][j]) ;   
+        
+        if(i>0 && fFillMCFractionHisto) // skip first entry in array, general case not filled
+        {
+          fhMCGenFracNLocMaxEbin[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%d",pname[i].Data(),j),
+                                                   Form("NLM vs E, %s, E bin %d",ptype[i].Data(),j),
+                                                   200,0,2,nMaxBins,0,nMaxBins); 
+          fhMCGenFracNLocMaxEbin[i][j]->SetYTitle("NLM");
+          fhMCGenFracNLocMaxEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
+          outputContainer->Add(fhMCGenFracNLocMaxEbin[i][j]) ;           
+          
+          fhMCGenFracNLocMaxEbinMatched[i][j]  = new TH2F(Form("hMCGenFracNLocMax%sEbin%dMatched",pname[i].Data(),j),
+                                                          Form("NLM vs E, %s, E bin %d, matched to a track",ptype[i].Data(),j),
+                                                          200,0,2,nMaxBins,0,nMaxBins); 
+          fhMCGenFracNLocMaxEbinMatched[i][j]->SetYTitle("NLM");
+          fhMCGenFracNLocMaxEbinMatched[i][j]->SetXTitle("E_{gen} / E_{reco}");
+          outputContainer->Add(fhMCGenFracNLocMaxEbinMatched[i][j]) ;   
+          
+          fhMassMCGenFracNLocMax1Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
+                                                        Form("Invariant mass of 2 highest energy cells vs E, %s, E bin %d",ptype[i].Data(),j),
+                                                        200,0,2,mbins,mmin,mmax); 
+          fhMassMCGenFracNLocMax1Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
+          fhMassMCGenFracNLocMax1Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
+          outputContainer->Add(fhMassMCGenFracNLocMax1Ebin[i][j]) ;   
+          
+          fhMassMCGenFracNLocMax2Ebin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
+                                                        Form("Invariant mass of 2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
+                                                        200,0,2,mbins,mmin,mmax); 
+          fhMassMCGenFracNLocMax2Ebin[i][j]->SetYTitle("M (GeV/c^{2})");
+          fhMassMCGenFracNLocMax2Ebin[i][j]->SetXTitle("E_{gen} / E_{reco}");
+          outputContainer->Add(fhMassMCGenFracNLocMax2Ebin[i][j]) ;   
+          
+          fhMassMCGenFracNLocMaxNEbin[i][j]  = new TH2F(Form("hMassMCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
+                                                        Form("Invariant mass of N>2 local maxima cells vs E, %s, E bin %d",ptype[i].Data(),j),
+                                                        200,0,2,mbins,mmin,mmax); 
+          fhMassMCGenFracNLocMaxNEbin[i][j]->SetYTitle("M (GeV/c^{2})");
+          fhMassMCGenFracNLocMaxNEbin[i][j]->SetXTitle("E_{gen} / E_{reco}");
+          outputContainer->Add(fhMassMCGenFracNLocMaxNEbin[i][j]) ;   
+          
+          fhM02MCGenFracNLocMax1Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax1%sEbin%d",pname[i].Data(),j),
+                                                          Form("#lambda_{0}^{2} vs E for N max  = 1 %s, E bin %d",ptype[i].Data(), j),
+                                                          200,0,2,ssbins,ssmin,ssmax); 
+          fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
+          fhM02MCGenFracNLocMax1Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
+          outputContainer->Add(fhM02MCGenFracNLocMax1Ebin[i][j]) ; 
+          
+          fhM02MCGenFracNLocMax2Ebin[i][j]     = new TH2F(Form("hM02MCGenFracNLocMax2%sEbin%d",pname[i].Data(),j),
+                                                          Form("#lambda_{0}^{2} vs E for N max  = 2 %s, E bin %d",ptype[i].Data(),j),
+                                                          200,0,2,ssbins,ssmin,ssmax); 
+          fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
+          fhM02MCGenFracNLocMax2Ebin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
+          outputContainer->Add(fhM02MCGenFracNLocMax2Ebin[i][j]) ; 
+          
+          fhM02MCGenFracNLocMaxNEbin[i][j]    = new TH2F(Form("hM02MCGenFracNLocMaxN%sEbin%d",pname[i].Data(),j),
+                                                         Form("#lambda_{0}^{2} vs E for N max  > 2 %s, E bin %d",ptype[i].Data(),j),
+                                                         200,0,2,ssbins,ssmin,ssmax); 
+          fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetYTitle("#lambda_{0}^{2}");
+          fhM02MCGenFracNLocMaxNEbin[i][j]   ->SetXTitle("E_{gen} / E_{reco}");
+          outputContainer->Add(fhM02MCGenFracNLocMaxNEbin[i][j]) ; 
+        }
+      }
   } // MC particle list
  
   for(Int_t i = 0; i < 4; i++)
   {  
+    if(IsDataMC())
+    {
+      fhMCAsymM02NLocMax1MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax1MCPi0Ebin%d",i),
+                                                 Form("Asymmetry of MC #pi^{0} of 2 highest energy cells #lambda_{0}^{2}, E bin %d",i),
+                                                 ssbins,ssmin,ssmax,100,0,1); 
+      fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
+      fhMCAsymM02NLocMax1MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
+      outputContainer->Add(fhMCAsymM02NLocMax1MCPi0Ebin[i]) ;   
+      
+      fhMCAsymM02NLocMax2MCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMax2MCPi0Ebin%d",i),
+                                                 Form("Asymmetry of MC #pi^{0} of 2 local maxima cells #lambda_{0}^{2}, E bin %d",i),
+                                                 ssbins,ssmin,ssmax,100,0,1); 
+      fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetYTitle("Decay asymmetry");
+      fhMCAsymM02NLocMax2MCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
+      outputContainer->Add(fhMCAsymM02NLocMax2MCPi0Ebin[i]) ;   
+      
+      fhMCAsymM02NLocMaxNMCPi0Ebin[i]  = new TH2F(Form("hMCAsymM02NLocMaxNMCPi0Ebin%d",i),
+                                                 Form("Asymmetry of MC #pi^{0} of N>2 local maxima cells vs #lambda_{0}^{2}, E bin %d",i),
+                                                 ssbins,ssmin,ssmax,100,0,1); 
+      fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetYTitle("Decay asymmetry");
+      fhMCAsymM02NLocMaxNMCPi0Ebin[i]->SetXTitle("#lambda_{0}^{2}");
+      outputContainer->Add(fhMCAsymM02NLocMaxNMCPi0Ebin[i]) ; 
+    }
+    
     fhMassM02NLocMax1Ebin[i]  = new TH2F(Form("hMassM02NLocMax1Ebin%d",i),
                                         Form("Invariant mass of 2 highest energy cells #lambda_{0}^{2}, E bin %d",i),
                                         ssbins,ssmin,ssmax,mbins,mmin,mmax); 
@@ -451,110 +994,207 @@ TList * AliAnaInsideClusterInvariantMass::GetCreateOutputObjects()
                                         ssbins,ssmin,ssmax,mbins,mmin,mmax); 
     fhMassM02NLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
     fhMassM02NLocMaxNEbin[i]->SetXTitle("#lambda_{0}^{2}");
-    outputContainer->Add(fhMassM02NLocMaxNEbin[i]) ;   
+    outputContainer->Add(fhMassM02NLocMaxNEbin[i]) ; 
+    
+    if(fFillSSExtraHisto)
+    {
+      fhMassDispEtaNLocMax1Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax1Ebin%d",i),
+                                               Form("Invariant mass of 2 highest energy cells #sigma_{#eta #eta}^{2}, E bin %d",i),
+                                               ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+      fhMassDispEtaNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
+      fhMassDispEtaNLocMax1Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
+      outputContainer->Add(fhMassDispEtaNLocMax1Ebin[i]) ;   
+      
+      fhMassDispEtaNLocMax2Ebin[i]  = new TH2F(Form("hMassDispEtaNLocMax2Ebin%d",i),
+                                               Form("Invariant mass of 2 local maxima cells #sigma_{#eta #eta}^{2}, E bin %d",i),
+                                               ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+      fhMassDispEtaNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
+      fhMassDispEtaNLocMax2Ebin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
+      outputContainer->Add(fhMassDispEtaNLocMax2Ebin[i]) ;   
+      
+      fhMassDispEtaNLocMaxNEbin[i]  = new TH2F(Form("hMassDispEtaNLocMaxNEbin%d",i),
+                                               Form("Invariant mass of N>2 local maxima cells vs #sigma_{#eta #eta}^{2}, E bin %d",i),
+                                               ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+      fhMassDispEtaNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
+      fhMassDispEtaNLocMaxNEbin[i]->SetXTitle("#sigma_{#eta #eta}^{2}");
+      outputContainer->Add(fhMassDispEtaNLocMaxNEbin[i]) ;   
+      
+      fhMassDispPhiNLocMax1Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax1Ebin%d",i),
+                                               Form("Invariant mass of 2 highest energy cells #sigma_{#phi #phi}^{2}, E bin %d",i),
+                                               ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+      fhMassDispPhiNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
+      fhMassDispPhiNLocMax1Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
+      outputContainer->Add(fhMassDispPhiNLocMax1Ebin[i]) ;   
+      
+      fhMassDispPhiNLocMax2Ebin[i]  = new TH2F(Form("hMassDispPhiNLocMax2Ebin%d",i),
+                                               Form("Invariant mass of 2 local maxima cells #sigma_{#phi #phi}^{2}, E bin %d",i),
+                                               ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+      fhMassDispPhiNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
+      fhMassDispPhiNLocMax2Ebin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
+      outputContainer->Add(fhMassDispPhiNLocMax2Ebin[i]) ;   
+      
+      fhMassDispPhiNLocMaxNEbin[i]  = new TH2F(Form("hMassDispPhiNLocMaxNEbin%d",i),
+                                               Form("Invariant mass of N>2 local maxima cells vs #sigma_{#phi #phi}^{2}, E bin %d",i),
+                                               ssbins,ssmin,ssmax,mbins,mmin,mmax); 
+      fhMassDispPhiNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
+      fhMassDispPhiNLocMaxNEbin[i]->SetXTitle("#sigma_{#phi #phi}^{2}");
+      outputContainer->Add(fhMassDispPhiNLocMaxNEbin[i]) ;   
+      
+      fhMassDispAsyNLocMax1Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax1Ebin%d",i),
+                                               Form("Invariant mass of 2 highest energy cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E bin %d",i),
+                                               200,-1,1,mbins,mmin,mmax); 
+      fhMassDispAsyNLocMax1Ebin[i]->SetYTitle("M (GeV/c^{2})");
+      fhMassDispAsyNLocMax1Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
+      outputContainer->Add(fhMassDispAsyNLocMax1Ebin[i]) ;   
+      
+      fhMassDispAsyNLocMax2Ebin[i]  = new TH2F(Form("hMassDispAsyNLocMax2Ebin%d",i),
+                                               Form("Invariant mass of 2 local maxima cells A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E bin %d",i),
+                                               200,-1,1,mbins,mmin,mmax); 
+      fhMassDispAsyNLocMax2Ebin[i]->SetYTitle("M (GeV/c^{2})");
+      fhMassDispAsyNLocMax2Ebin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
+      outputContainer->Add(fhMassDispAsyNLocMax2Ebin[i]) ;   
+      
+      fhMassDispAsyNLocMaxNEbin[i]  = new TH2F(Form("hMassDispAsyNLocMaxNEbin%d",i),
+                                               Form("Invariant mass of N>2 local maxima cells vs A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2}), E bin %d",i),
+                                               200,-1,1,mbins,mmin,mmax); 
+      fhMassDispAsyNLocMaxNEbin[i]->SetYTitle("M (GeV/c^{2})");
+      fhMassDispAsyNLocMaxNEbin[i]->SetXTitle("A = (#sigma_{#phi #phi}^{2} - #sigma_{#eta #eta}^{2}) / (#sigma_{#phi #phi}^{2} + #sigma_{#eta #eta}^{2})");
+      outputContainer->Add(fhMassDispAsyNLocMaxNEbin[i]) ;   
+    }
   }  
   
-  for(Int_t i = 0; i < n; i++)
-  {  
-    fhTrackMatchedDEtaLocMax1[i]  = new TH2F
-    (Form("hTrackMatchedDEtaLocMax1%s",pname[i].Data()),
-     Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
-     nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
-    fhTrackMatchedDEtaLocMax1[i]->SetYTitle("d#eta");
-    fhTrackMatchedDEtaLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
-    
-    fhTrackMatchedDPhiLocMax1[i]  = new TH2F
-    (Form("hTrackMatchedDPhiLocMax1%s",pname[i].Data()),
-     Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
-     nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
-    fhTrackMatchedDPhiLocMax1[i]->SetYTitle("d#phi (rad)");
-    fhTrackMatchedDPhiLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
-    
-    outputContainer->Add(fhTrackMatchedDEtaLocMax1[i]) ; 
-    outputContainer->Add(fhTrackMatchedDPhiLocMax1[i]) ;
-
-    fhTrackMatchedDEtaLocMax2[i]  = new TH2F
-    (Form("hTrackMatchedDEtaLocMax2%s",pname[i].Data()),
-     Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
-     nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
-    fhTrackMatchedDEtaLocMax2[i]->SetYTitle("d#eta");
-    fhTrackMatchedDEtaLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
-    
-    fhTrackMatchedDPhiLocMax2[i]  = new TH2F
-    (Form("hTrackMatchedDPhiLocMax2%s",pname[i].Data()),
-     Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
-     nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
-    fhTrackMatchedDPhiLocMax2[i]->SetYTitle("d#phi (rad)");
-    fhTrackMatchedDPhiLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
-    
-    outputContainer->Add(fhTrackMatchedDEtaLocMax2[i]) ; 
-    outputContainer->Add(fhTrackMatchedDPhiLocMax2[i]) ;
-
-    fhTrackMatchedDEtaLocMaxN[i]  = new TH2F
-    (Form("hTrackMatchedDEtaLocMaxN%s",pname[i].Data()),
-     Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
-     nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
-    fhTrackMatchedDEtaLocMaxN[i]->SetYTitle("d#eta");
-    fhTrackMatchedDEtaLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
-    
-    fhTrackMatchedDPhiLocMaxN[i]  = new TH2F
-    (Form("hTrackMatchedDPhiLocMaxN%s",pname[i].Data()),
-     Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
-     nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
-    fhTrackMatchedDPhiLocMaxN[i]->SetYTitle("d#phi (rad)");
-    fhTrackMatchedDPhiLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
-    
-    outputContainer->Add(fhTrackMatchedDEtaLocMaxN[i]) ; 
-    outputContainer->Add(fhTrackMatchedDPhiLocMaxN[i]) ;    
-    
+  if(fFillTMResidualHisto)
+  {
+    for(Int_t i = 0; i < n; i++)
+    {  
+      
+      fhTrackMatchedDEtaLocMax1[i]  = new TH2F
+      (Form("hTrackMatchedDEtaLocMax1%s",pname[i].Data()),
+       Form("d#eta of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
+       nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
+      fhTrackMatchedDEtaLocMax1[i]->SetYTitle("d#eta");
+      fhTrackMatchedDEtaLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
+      
+      fhTrackMatchedDPhiLocMax1[i]  = new TH2F
+      (Form("hTrackMatchedDPhiLocMax1%s",pname[i].Data()),
+       Form("d#phi of cluster-track vs cluster energy, 1 Local Maxima, %s",ptype[i].Data()),
+       nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
+      fhTrackMatchedDPhiLocMax1[i]->SetYTitle("d#phi (rad)");
+      fhTrackMatchedDPhiLocMax1[i]->SetXTitle("E_{cluster} (GeV)");
+      
+      outputContainer->Add(fhTrackMatchedDEtaLocMax1[i]) ; 
+      outputContainer->Add(fhTrackMatchedDPhiLocMax1[i]) ;
+      
+      fhTrackMatchedDEtaLocMax2[i]  = new TH2F
+      (Form("hTrackMatchedDEtaLocMax2%s",pname[i].Data()),
+       Form("d#eta of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
+       nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
+      fhTrackMatchedDEtaLocMax2[i]->SetYTitle("d#eta");
+      fhTrackMatchedDEtaLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
+      
+      fhTrackMatchedDPhiLocMax2[i]  = new TH2F
+      (Form("hTrackMatchedDPhiLocMax2%s",pname[i].Data()),
+       Form("d#phi of cluster-track vs cluster energy, 2 Local Maxima, %s",ptype[i].Data()),
+       nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
+      fhTrackMatchedDPhiLocMax2[i]->SetYTitle("d#phi (rad)");
+      fhTrackMatchedDPhiLocMax2[i]->SetXTitle("E_{cluster} (GeV)");
+      
+      outputContainer->Add(fhTrackMatchedDEtaLocMax2[i]) ; 
+      outputContainer->Add(fhTrackMatchedDPhiLocMax2[i]) ;
+      
+      fhTrackMatchedDEtaLocMaxN[i]  = new TH2F
+      (Form("hTrackMatchedDEtaLocMaxN%s",pname[i].Data()),
+       Form("d#eta of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
+       nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax); 
+      fhTrackMatchedDEtaLocMaxN[i]->SetYTitle("d#eta");
+      fhTrackMatchedDEtaLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
+      
+      fhTrackMatchedDPhiLocMaxN[i]  = new TH2F
+      (Form("hTrackMatchedDPhiLocMaxN%s",pname[i].Data()),
+       Form("d#phi of cluster-track vs cluster energy, N>2 Local Maxima, %s",ptype[i].Data()),
+       nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax); 
+      fhTrackMatchedDPhiLocMaxN[i]->SetYTitle("d#phi (rad)");
+      fhTrackMatchedDPhiLocMaxN[i]->SetXTitle("E_{cluster} (GeV)");
+      
+      outputContainer->Add(fhTrackMatchedDEtaLocMaxN[i]) ; 
+      outputContainer->Add(fhTrackMatchedDPhiLocMaxN[i]) ;    
+    }
+  }
+  
+  if(fFillAngleHisto)
+  {
+    for(Int_t j = 0; j < 2; j++)
+    {  
+      
+      fhAnglePairLocMax1[j]  = new TH2F(Form("hAnglePairLocMax1%s",sMatched[j].Data()),
+                                        Form("Opening angle of 2 highest energy cells vs pair Energy, %s",sMatched[j].Data()),
+                                        nptbins,ptmin,ptmax,200,0,0.2); 
+      fhAnglePairLocMax1[j]->SetYTitle("#alpha (rad)");
+      fhAnglePairLocMax1[j]->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAnglePairLocMax1[j]) ;   
+      
+      fhAnglePairLocMax2[j]  = new TH2F(Form("hAnglePairLocMax2%s",sMatched[j].Data()),
+                                        Form("Opening angle of 2 local maxima cells vs Energy, %s",sMatched[j].Data()),
+                                        nptbins,ptmin,ptmax,200,0,0.2); 
+      fhAnglePairLocMax2[j]->SetYTitle("#alpha (rad)");
+      fhAnglePairLocMax2[j]->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAnglePairLocMax2[j]) ;   
+      
+      fhAnglePairLocMaxN[j]  = new TH2F(Form("hAnglePairLocMaxN%s",sMatched[j].Data()),
+                                        Form("Opening angle of N>2 local maxima cells vs Energy, %s",sMatched[j].Data()),
+                                        nptbins,ptmin,ptmax,200,0,0.2); 
+      fhAnglePairLocMaxN[j]->SetYTitle("#alpha (rad)");
+      fhAnglePairLocMaxN[j]->SetXTitle("E (GeV)");
+      outputContainer->Add(fhAnglePairLocMaxN[j]) ;   
+      
+      fhAnglePairMassLocMax1[j]  = new TH2F(Form("hAnglePairMassLocMax1%s",sMatched[j].Data()),
+                                            Form("Opening angle of 2 highest energy cells vs Mass for E > 8 GeV, %s",sMatched[j].Data()),
+                                            mbins,mmin,mmax,200,0,0.2); 
+      fhAnglePairMassLocMax1[j]->SetXTitle("M (GeV/c^{2})");
+      fhAnglePairMassLocMax1[j]->SetYTitle("#alpha (rad)");
+      outputContainer->Add(fhAnglePairMassLocMax1[j]) ;   
+      
+      fhAnglePairMassLocMax2[j]  = new TH2F(Form("hAnglePairMassLocMax2%s",sMatched[j].Data()),
+                                            Form("Opening angle of 2 local maxima cells vs Mass for E > 8 GeV, %s",sMatched[j].Data()),
+                                            mbins,mmin,mmax,200,0,0.2); 
+      fhAnglePairMassLocMax2[j]->SetXTitle("M (GeV/c^{2})");
+      fhAnglePairMassLocMax2[j]->SetYTitle("#alpha (rad)");
+      outputContainer->Add(fhAnglePairMassLocMax2[j]) ;   
+      
+      fhAnglePairMassLocMaxN[j]  = new TH2F(Form("hAnglePairMassLocMaxN%s",sMatched[j].Data()),
+                                            Form("Opening angle of N>2 local maxima cells vs Mass for E > 8 GeV, %s",sMatched[j].Data()),
+                                            mbins,mmin,mmax,200,0,0.2); 
+      fhAnglePairMassLocMaxN[j]->SetXTitle("M (GeV/c^{2})");
+      fhAnglePairMassLocMaxN[j]->SetYTitle("#alpha (rad)");
+      outputContainer->Add(fhAnglePairMassLocMaxN[j]) ;  
+      
+    }
   }
   
   for(Int_t j = 0; j < 2; j++)
   {  
-    
-    fhAnglePairLocMax1[j]  = new TH2F(Form("hAnglePairLocMax1%s",sMatched[j].Data()),
-                                      Form("Opening angle of 2 highest energy cells vs pair Energy, %s",sMatched[j].Data()),
-                                      nptbins,ptmin,ptmax,200,0,0.2); 
-    fhAnglePairLocMax1[j]->SetYTitle("#alpha (rad)");
-    fhAnglePairLocMax1[j]->SetXTitle("E (GeV)");
-    outputContainer->Add(fhAnglePairLocMax1[j]) ;   
-    
-    fhAnglePairLocMax2[j]  = new TH2F(Form("hAnglePairLocMax2%s",sMatched[j].Data()),
-                                      Form("Opening angle of 2 local maxima cells vs Energy, %s",sMatched[j].Data()),
-                                      nptbins,ptmin,ptmax,200,0,0.2); 
-    fhAnglePairLocMax2[j]->SetYTitle("#alpha (rad)");
-    fhAnglePairLocMax2[j]->SetXTitle("E (GeV)");
-    outputContainer->Add(fhAnglePairLocMax2[j]) ;   
-    
-    fhAnglePairLocMaxN[j]  = new TH2F(Form("hAnglePairLocMaxN%s",sMatched[j].Data()),
-                                      Form("Opening angle of N>2 local maxima cells vs Energy, %s",sMatched[j].Data()),
-                                      nptbins,ptmin,ptmax,200,0,0.2); 
-    fhAnglePairLocMaxN[j]->SetYTitle("#alpha (rad)");
-    fhAnglePairLocMaxN[j]->SetXTitle("E (GeV)");
-    outputContainer->Add(fhAnglePairLocMaxN[j]) ;   
-    
-    fhAnglePairMassLocMax1[j]  = new TH2F(Form("hAnglePairMassLocMax1%s",sMatched[j].Data()),
-                                          Form("Opening angle of 2 highest energy cells vs Mass for E > 7 GeV, %s",sMatched[j].Data()),
-                                          mbins,mmin,mmax,200,0,0.2); 
-    fhAnglePairMassLocMax1[j]->SetXTitle("M (GeV/c^{2})");
-    fhAnglePairMassLocMax1[j]->SetYTitle("#alpha (rad)");
-    outputContainer->Add(fhAnglePairMassLocMax1[j]) ;   
-    
-    fhAnglePairMassLocMax2[j]  = new TH2F(Form("hAnglePairMassLocMax2%s",sMatched[j].Data()),
-                                          Form("Opening angle of 2 local maxima cells vs Mass for E > 7 GeV, %s",sMatched[j].Data()),
-                                          mbins,mmin,mmax,200,0,0.2); 
-    fhAnglePairMassLocMax2[j]->SetXTitle("M (GeV/c^{2})");
-    fhAnglePairMassLocMax2[j]->SetYTitle("#alpha (rad)");
-    outputContainer->Add(fhAnglePairMassLocMax2[j]) ;   
-    
-    fhAnglePairMassLocMaxN[j]  = new TH2F(Form("hAnglePairMassLocMaxN%s",sMatched[j].Data()),
-                                          Form("Opening angle of N>2 local maxima cells vs Mass for E > 7 GeV, %s",sMatched[j].Data()),
-                                          mbins,mmin,mmax,200,0,0.2); 
-    fhAnglePairMassLocMaxN[j]->SetXTitle("M (GeV/c^{2})");
-    fhAnglePairMassLocMaxN[j]->SetYTitle("#alpha (rad)");
-    outputContainer->Add(fhAnglePairMassLocMaxN[j]) ;  
-    
+  fhSplitEFractionvsAsyNLocMax1[j]     = new TH2F(Form("hSplitEFractionvsAsyNLocMax1%s",sMatched[j].Data()),
+                                                Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  = 1, E>8, %s",sMatched[j].Data()),
+                                                100,-1,1,120,0,1.2); 
+  fhSplitEFractionvsAsyNLocMax1[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
+  fhSplitEFractionvsAsyNLocMax1[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
+  outputContainer->Add(fhSplitEFractionvsAsyNLocMax1[j]) ; 
+  
+  fhSplitEFractionvsAsyNLocMax2[j]     = new TH2F(Form("hSplitEFractionvsAsyNLocMax2%s",sMatched[j].Data()),
+                                                Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  = 2,E>8, %s",sMatched[j].Data()),
+                                                100,-1,1,120,0,1.2); 
+  fhSplitEFractionvsAsyNLocMax2[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
+  fhSplitEFractionvsAsyNLocMax2[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
+  outputContainer->Add(fhSplitEFractionvsAsyNLocMax2[j]) ; 
+  
+  fhSplitEFractionvsAsyNLocMaxN[j]    = new TH2F(Form("hSplitEFractionvsAsyNLocMaxN%s",sMatched[j].Data()),
+                                               Form("(E1+E2)/E_{cluster} vs (E_{split1}-E_{split2})/(E_{split1}+E_{split2}) for N max  > 2, E>8, %s",sMatched[j].Data()),
+                                               100,-1,1,120,0,1.2); 
+  fhSplitEFractionvsAsyNLocMaxN[j]   ->SetXTitle("(E_{split1}-E_{split2})/(E_{split1}+E_{split2})");
+  fhSplitEFractionvsAsyNLocMaxN[j]   ->SetYTitle("(E_{split1}+E_{split2})/E_{cluster}");
+  outputContainer->Add(fhSplitEFractionvsAsyNLocMaxN[j]) ; 
   }
+   
   
   return outputContainer ;
   
@@ -565,16 +1205,19 @@ void AliAnaInsideClusterInvariantMass::Init()
 { 
   //Init
   //Do some checks
-  if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
+  if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD())
+  {
     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
     abort();
   }
-  else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
+  else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD())
+  {
     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
     abort();
   }
   
-  if( GetReader()->GetDataType() == AliCaloTrackReader::kMC ){
+  if( GetReader()->GetDataType() == AliCaloTrackReader::kMC )
+  {
     printf("AliAnaInsideClusterInvariantMass::Init() - !!STOP: You want to use pure MC data!!\n");
     abort();
     
@@ -594,16 +1237,8 @@ void AliAnaInsideClusterInvariantMass::InitParameters()
   fM02MaxCut   = 10 ;
   
   fMinNCells   = 4 ;
-  
-  fMassEtaMin  = 0.4;
-  fMassEtaMax  = 0.6;
-  
-  fMassPi0Min  = 0.08;
-  fMassPi0Max  = 0.20;
-  
-  fMassConMin  = 0.0;
-  fMassConMax  = 0.05;
-  
+  fMinBadDist  = 2 ;
+    
 }
 
 
@@ -616,69 +1251,99 @@ void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
   AliVCaloCells* cells = 0x0;
 
   //Select the Calorimeter of the photon
-  if(fCalorimeter == "PHOS"){
+  if(fCalorimeter == "PHOS")
+  {
     pl    = GetPHOSClusters();
     cells = GetPHOSCells();
   }
-  else if (fCalorimeter == "EMCAL"){
+  else if (fCalorimeter == "EMCAL")
+  {
     pl    = GetEMCALClusters();
     cells = GetEMCALCells();
   }
   
-  if(!pl || !cells) {
+  const Float_t ecut = 8.; // Fixed cut for some histograms
+  
+  if(!pl || !cells) 
+  {
     Info("MakeAnalysisFillHistograms","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
     return;
   }  
   
        if(fCalorimeter == "PHOS") return; // Not implemented for PHOS yet
 
-  for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++){
+  for(Int_t icluster = 0; icluster < pl->GetEntriesFast(); icluster++)
+  {
     AliVCluster * cluster = (AliVCluster*) (pl->At(icluster)); 
 
     // Study clusters with large shape parameter
     Float_t en = cluster->E();
     Float_t l0 = cluster->GetM02();
     Int_t   nc = cluster->GetNCells();
+    Float_t bd = cluster->GetDistanceToBadChannel() ; 
 
-    //If too small or big E or low number of cells, skip it
-    if( en < GetMinEnergy() || en > GetMaxEnergy() || nc < fMinNCells) continue ; 
-  
-    Int_t    absId1    = -1; Int_t absId2 = -1;
-    Int_t   *absIdList = new Int_t  [nc]; 
-    Float_t *maxEList  = new Float_t[nc]; 
-    Int_t    nMax      = GetCaloUtils()->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
-    Bool_t   matched   = IsTrackMatched(cluster,GetReader()->GetInputEvent());
     
+    //If too small or big E or low number of cells, or close to a bad channel skip it
+    if( en < GetMinEnergy() || en > GetMaxEnergy() || nc < fMinNCells || bd < fMinBadDist) continue ; 
+    
+    //printf("en %2.2f, GetMinEnergy() %2.2f, GetMaxEnergy() %2.2f, nc %d, fMinNCells %d,  bd %2.2f, fMinBadDist %2.2f\n",
+    //       en,GetMinEnergy(), GetMaxEnergy(), nc, fMinNCells, bd, fMinBadDist);
+    
+    // Get more Shower Shape parameters
+    Float_t ll0  = 0., ll1  = 0.;
+    Float_t disp= 0., dispEta = 0., dispPhi    = 0.; 
+    Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;  
+   
+    GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
+                                                                                 ll0, ll1, disp, dispEta, dispPhi, sEta, sPhi, sEtaPhi);
+    
+    Float_t dispAsy = -1;
+    if(dispEta+dispPhi >0 ) dispAsy = (dispPhi-dispEta) / (dispPhi+dispEta);
+    
+    Int_t    nMax = 0;
+    Double_t mass = 0., angle = 0.;
+    Double_t e1   = 0., e2    = 0.;
+    Int_t pidTag = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(cluster,cells,GetCaloUtils(),
+                                                                               GetVertex(0), nMax, mass, angle,e1,e2);    
     if (nMax <= 0) 
     {
-      printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found!\n");
+      if(GetDebug() > 0 )
+        printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - No local maximum found! It did not pass CaloPID selection criteria \n");
       
-      /*
-      for(Int_t iDigit  = 0; iDigit < cluster->GetNCells(); iDigit++ ) {
-        Float_t ec = cells->GetCellAmplitude(cluster->GetCellsAbsId()[iDigit]);
-        GetCaloUtils()->RecalibrateCellAmplitude(ec,fCalorimeter,cluster->GetCellsAbsId()[iDigit]);
-        printf("iDigit %d, absId %d, Ecell %f\n",iDigit,cluster->GetCellsAbsId()[iDigit], ec);
-      }
-      */
-      
-      delete [] absIdList ;
-      delete [] maxEList  ;
       return;
     }
     
-    fhNLocMax[0][matched]->Fill(en,nMax);
-    for(Int_t imax = 0; imax < nMax; imax++)
-    {
-      fhNLocMaxEMax [0][matched]->Fill(maxEList[imax]   ,nMax);
-      fhNLocMaxEFrac[0][matched]->Fill(maxEList[imax]/en,nMax);
-    }
+    Float_t splitFrac = (e1+e2)/en;
+    Float_t asym = -10;
+    if(e1+e2>0) asym = (e1-e2)/(e1+e2);
+    
+    Bool_t  matched   = IsTrackMatched(cluster,GetReader()->GetInputEvent());
     
+    fhNLocMax[0][matched]->Fill(en,nMax);
     
-    if     ( nMax == 1  ) { fhM02NLocMax1[0][matched]->Fill(en,l0) ; fhNCellNLocMax1[0][matched]->Fill(en,nc) ; }
-    else if( nMax == 2  ) { fhM02NLocMax2[0][matched]->Fill(en,l0) ; fhNCellNLocMax2[0][matched]->Fill(en,nc) ; }
-    else if( nMax >= 3  ) { fhM02NLocMaxN[0][matched]->Fill(en,l0) ; fhNCellNLocMaxN[0][matched]->Fill(en,nc) ; }
+    if     ( nMax == 1  ) 
+    { 
+      fhM02NLocMax1[0][matched]->Fill(en,l0) ; 
+      fhSplitEFractionNLocMax1[0][matched]->Fill(en,splitFrac) ; 
+      if(en > ecut) fhSplitEFractionvsAsyNLocMax1[matched]->Fill(asym,splitFrac) ; 
+      if(fFillSSExtraHisto) fhNCellNLocMax1[0][matched]->Fill(en,nc) ; 
+    }
+    else if( nMax == 2  ) 
+    { 
+      fhM02NLocMax2[0][matched]->Fill(en,l0) ; 
+      fhSplitEFractionNLocMax2[0][matched]->Fill(en,splitFrac) ; 
+      if(en > ecut) fhSplitEFractionvsAsyNLocMax2[matched]->Fill(asym,splitFrac) ; 
+      if(fFillSSExtraHisto) fhNCellNLocMax2[0][matched]->Fill(en,nc) ; }
+    else if( nMax >= 3  ) 
+    { 
+      fhM02NLocMaxN[0][matched]->Fill(en,l0) ; 
+      fhSplitEFractionNLocMaxN[0][matched]->Fill(en,splitFrac) ; 
+      if(en > ecut) fhSplitEFractionvsAsyNLocMaxN[matched]->Fill(asym,splitFrac) ; 
+      if(fFillSSExtraHisto) fhNCellNLocMaxN[0][matched]->Fill(en,nc) ; 
+    }
     else printf("N max smaller than 1 -> %d \n",nMax);
     
+    
     Float_t dZ  = cluster->GetTrackDz();
     Float_t dR  = cluster->GetTrackDx();
     
@@ -689,18 +1354,21 @@ void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
     }    
     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
     
-    if(TMath::Abs(dR) < 999)
+    if(TMath::Abs(dR) < 999 && fFillTMResidualHisto)
     {
       if     ( nMax == 1  ) { fhTrackMatchedDEtaLocMax1[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMax1[0]->Fill(en,dR); }
       else if( nMax == 2  ) { fhTrackMatchedDEtaLocMax2[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMax2[0]->Fill(en,dR); }
       else if( nMax >= 3  ) { fhTrackMatchedDEtaLocMaxN[0]->Fill(en,dZ); fhTrackMatchedDPhiLocMaxN[0]->Fill(en,dR); }
     }
-
+    
     // Play with the MC stack if available
     // Check origin of the candidates
-    Int_t mcindex = -1;
-    if(IsDataMC()){
-      
+    Int_t mcindex   = -1;
+    Float_t eprim   =  0;
+    Float_t asymGen = -2; 
+    Int_t mcLabel   = cluster->GetLabel();
+    if(IsDataMC())
+    {
       Int_t tag        = GetMCAnalysisUtils()->CheckOrigin(cluster->GetLabels(),cluster->GetNLabels(), GetReader(), 0);
             
       if      ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )      mcindex = kmcPi0;
@@ -711,225 +1379,347 @@ void  AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms()
                 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)) mcindex = kmcConversion;
       else if ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))   mcindex = kmcElectron;
       else                                                                                mcindex = kmcHadron;
-      
-      //GetMCAnalysisUtils()->PrintMCTag(tag);
-      //printf("\t MC index Assigned %d \n",mcindex);
-      fhNLocMaxNLabel[0]      [matched]->Fill(en,nMax);
-      fhNLocMaxNLabel[mcindex][matched]->Fill(cluster->GetNLabels(),nMax);
 
       fhNLocMax[mcindex][matched]->Fill(en,nMax);
-      for(Int_t imax = 0; imax < nMax; imax++)
-      {
-        fhNLocMaxEMax [mcindex][matched]->Fill(maxEList[imax]   ,nMax);
-        fhNLocMaxEFrac[mcindex][matched]->Fill(maxEList[imax]/en,nMax);
-      }
-      if     (nMax == 1 ) { fhM02NLocMax1[mcindex][matched]->Fill(en,l0) ; fhNCellNLocMax1[mcindex][matched]->Fill(en,nc) ; }
-      else if(nMax == 2 ) { fhM02NLocMax2[mcindex][matched]->Fill(en,l0) ; fhNCellNLocMax2[mcindex][matched]->Fill(en,nc) ; }
-      else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex][matched]->Fill(en,l0) ; fhNCellNLocMaxN[mcindex][matched]->Fill(en,nc) ; }
+            
+      if     (nMax == 1 ) { fhM02NLocMax1[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax1[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax1[mcindex][matched]->Fill(en,nc) ; }
+      else if(nMax == 2 ) { fhM02NLocMax2[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMax2[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMax2[mcindex][matched]->Fill(en,nc) ; }
+      else if(nMax >= 3 ) { fhM02NLocMaxN[mcindex][matched]->Fill(en,l0) ; fhSplitEFractionNLocMaxN[mcindex][matched]->Fill(en,splitFrac) ; if(fFillSSExtraHisto) fhNCellNLocMaxN[mcindex][matched]->Fill(en,nc) ; }
       
-      if(TMath::Abs(dR) < 999)
+      if(TMath::Abs(dR) < 999 && fFillTMResidualHisto)
       {
         if     ( nMax == 1  ) { fhTrackMatchedDEtaLocMax1[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMax1[mcindex]->Fill(en,dR); }
         else if( nMax == 2  ) { fhTrackMatchedDEtaLocMax2[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMax2[mcindex]->Fill(en,dR); }
         else if( nMax >= 3  ) { fhTrackMatchedDEtaLocMaxN[mcindex]->Fill(en,dZ); fhTrackMatchedDPhiLocMaxN[mcindex]->Fill(en,dR); }
       }
       
-    }  
+      Bool_t ok = kFALSE;
+      TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
+      eprim = primary.E();
+      
+      if(mcindex == kmcPi0 || mcindex == kmcEta)
+      {
+        if(mcindex == kmcPi0) 
+        {
+          asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,111,GetReader(),ok));
+          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok); 
+          if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
+        }
+        else 
+        {
+          asymGen = TMath::Abs(GetMCAnalysisUtils()->GetMCDecayAsymmetryForPDG(mcLabel,221,GetReader(),ok));
+          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok); 
+          if(grandmom.E() > 0 && ok) eprim =  grandmom.E();
+        }
+      }
+    } 
     
-    //---------------------------------------------------------------------
-    // From here only if M02 is large, fill histograms or split the cluster
-    //---------------------------------------------------------------------
+    Float_t efrac      = eprim/en;
+    Float_t efracSplit = 0;
+    if(e1+e2 > 0) efracSplit = eprim/(e1+e2);
 
-    if( l0 < fM02MinCut || l0 > fM02MaxCut ) 
+    //printf("e1 %2.2f, e2 %2.2f, eprim %2.2f, ereco %2.2f, esplit/ereco %2.2f, egen/ereco %2.2f, egen/esplit %2.2f\n",
+    //       e1,e2,eprim,en,splitFrac,efrac,efracSplit);
+    
+    Int_t ebin = -1;
+    if(en > 8  && en <= 12) ebin = 0; 
+    if(en > 12 && en <= 16) ebin = 1;
+    if(en > 16 && en <= 20) ebin = 2;
+    if(en > 20)             ebin = 3; 
+    
+    if(ebin >= 0 && IsDataMC() && fFillMCFractionHisto)
     {
-      delete [] absIdList ;
-      delete [] maxEList  ;
-      continue;    
+      if( !matched ) fhMCGenFracNLocMaxEbin       [mcindex][ebin]->Fill(efrac,nMax);
+      else           fhMCGenFracNLocMaxEbinMatched[mcindex][ebin]->Fill(efrac,nMax);
     }
-        
-    fhNLocMaxM02Cut[0][matched]->Fill(en,nMax);
-    if(IsDataMC()) fhNLocMaxM02Cut[mcindex][matched]->Fill(en,nMax);
     
-    //---------------------------------------------------------------------
-    // Get the 2 max indeces and do inv mass
-    //---------------------------------------------------------------------
-
-    if     ( nMax == 2 ) {
-      absId1 = absIdList[0];
-      absId2 = absIdList[1];
-    }
-    else if( nMax == 1 )
-    {
+    if     (nMax==1) 
+    { 
+      if( en > ecut ) 
+      {      
+        fhMassM02NLocMax1    [0][matched]->Fill(l0     ,  mass ); 
+        if(fFillSSExtraHisto)
+        {
+          fhMassDispEtaNLocMax1[0][matched]->Fill(dispEta,  mass ); 
+          fhMassDispPhiNLocMax1[0][matched]->Fill(dispPhi,  mass ); 
+          fhMassDispAsyNLocMax1[0][matched]->Fill(dispAsy,  mass );
+        }
+        
+        if(IsDataMC()) 
+        {
+          fhMassM02NLocMax1          [mcindex][matched]->Fill(l0     ,  mass  ); 
+          if(fFillMCFractionHisto)
+          {
+            fhMCGenFracNLocMax1        [mcindex][matched]->Fill(en     ,  efrac ); 
+            fhMCGenSplitEFracNLocMax1  [mcindex][matched]->Fill(en     ,  efracSplit ); 
+            fhMCGenEvsSplitENLocMax1   [mcindex][matched]->Fill(eprim  ,  e1+e2); 
+            fhMCGenEFracvsSplitEFracNLocMax1[mcindex][matched]->Fill(efrac,splitFrac ); 
+          }
+          
+          if(!matched && ebin >= 0)
+          {
+            if(fFillMCFractionHisto)
+            {
+              fhM02MCGenFracNLocMax1Ebin [mcindex][ebin]->Fill(efrac  ,  l0    ); 
+              fhMassMCGenFracNLocMax1Ebin[mcindex][ebin]->Fill(efrac  ,  mass  ); 
+            }
+            fhMCAsymM02NLocMax1MCPi0Ebin        [ebin]->Fill(l0     ,  asymGen );
+          }
+          
+          if(fFillSSExtraHisto)
+          {
+            fhMassDispEtaNLocMax1[mcindex][matched]->Fill(dispEta,  mass ); 
+            fhMassDispPhiNLocMax1[mcindex][matched]->Fill(dispPhi,  mass ); 
+            fhMassDispAsyNLocMax1[mcindex][matched]->Fill(dispAsy,  mass ); 
+          }
+        }
+      }
       
-      absId1 = absIdList[0];
-
-      //Find second highest energy cell
+      if(!matched && ebin >= 0)
+      {
+        fhMassSplitEFractionNLocMax1Ebin[0][ebin]->Fill(splitFrac,  mass);
+        if(IsDataMC())fhMassSplitEFractionNLocMax1Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
 
-      Float_t enmax = 0 ;
-      for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++){
-        Int_t absId = cluster->GetCellsAbsId()[iDigit];
-        if( absId == absId1 ) continue ; 
-        Float_t endig = cells->GetCellAmplitude(absId);
-        GetCaloUtils()->RecalibrateCellAmplitude(endig,fCalorimeter,absId); 
-        if(endig > enmax) {
-          enmax  = endig ;
-          absId2 = absId ;
+        fhMassM02NLocMax1Ebin    [ebin]->Fill(l0     ,  mass    );
+        if(fFillSSExtraHisto)
+        {
+          fhMassDispEtaNLocMax1Ebin[ebin]->Fill(dispEta,  mass );
+          fhMassDispPhiNLocMax1Ebin[ebin]->Fill(dispPhi,  mass );
+          fhMassDispAsyNLocMax1Ebin[ebin]->Fill(dispAsy,  mass );
         }
-      }// cell loop
-    }// 1 maxima 
-    else {  // n max > 2
-      // loop on maxima, find 2 highest
-      
-      // First max
-      Float_t enmax = 0 ;
-      for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
-        Float_t endig = maxEList[iDigit];
-        if(endig > enmax) {
-          enmax  = endig ;
-          absId1 = absIdList[iDigit];
+      }
+    }  
+    else if(nMax==2) 
+    {
+      if( en > ecut ) 
+      {      
+        fhMassM02NLocMax2    [0][matched]->Fill(l0     ,  mass ); 
+        if(fFillSSExtraHisto)
+        {
+          fhMassDispEtaNLocMax2[0][matched]->Fill(dispEta,  mass ); 
+          fhMassDispPhiNLocMax2[0][matched]->Fill(dispPhi,  mass ); 
+          fhMassDispAsyNLocMax2[0][matched]->Fill(dispAsy,  mass ); 
         }
-      }// first maxima loop
+        
+        if(IsDataMC()) 
+        {
+          fhMassM02NLocMax2        [mcindex][matched]->Fill(l0     ,  mass ); 
+          if(fFillMCFractionHisto)
+          {
+            fhMCGenFracNLocMax2      [mcindex][matched]->Fill(en     ,  efrac ); 
+            fhMCGenSplitEFracNLocMax2[mcindex][matched]->Fill(en     ,  efracSplit ); 
+            fhMCGenEvsSplitENLocMax2 [mcindex][matched]->Fill(eprim  ,  e1+e2); 
+            fhMCGenEFracvsSplitEFracNLocMax2[mcindex][matched]->Fill(efrac,splitFrac ); 
+          }
+          
+          if(!matched && ebin >= 0)
+          {
+            if(fFillMCFractionHisto)
+            {
+              fhM02MCGenFracNLocMax2Ebin [mcindex][ebin]->Fill(efrac  ,  l0    ); 
+              fhMassMCGenFracNLocMax2Ebin[mcindex][ebin]->Fill(efrac  ,  mass  ); 
+            }
+            fhMCAsymM02NLocMax2MCPi0Ebin        [ebin]->Fill(l0     ,  asymGen );
 
-      // Second max 
-      Float_t enmax2 = 0;
-      for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++){
-        if(absIdList[iDigit]==absId1) continue;
-        Float_t endig = maxEList[iDigit];
-        if(endig > enmax2) {
-          enmax2  = endig ;
-          absId2 = absIdList[iDigit];
+          }
+          if(fFillSSExtraHisto)
+          {
+            fhMassDispEtaNLocMax2[mcindex][matched]->Fill(dispEta,  mass ); 
+            fhMassDispPhiNLocMax2[mcindex][matched]->Fill(dispPhi,  mass ); 
+            fhMassDispAsyNLocMax2[mcindex][matched]->Fill(dispAsy,  mass ); 
+          }
         }
-      }// second maxima loop
-
-    } // n local maxima > 2
-    
-    //---------------------------------------------------------------------
-    // Split the cluster energy in 2, around the highest 2 local maxima
-    //---------------------------------------------------------------------
+      }
+      
+      if(!matched && ebin >= 0)
+      {
+        fhMassSplitEFractionNLocMax2Ebin[0][ebin]->Fill(splitFrac,  mass);
+        if(IsDataMC())fhMassSplitEFractionNLocMax2Ebin[mcindex][ebin]->Fill(splitFrac,  mass);
 
-    //Float_t en1 = 0, en2 = 0;
-    //SplitEnergy(absId1,absId2,cluster, cells, en1, en2, nMax /*absIdList, maxEList,*/);
+        fhMassM02NLocMax2Ebin    [ebin]->Fill(l0     ,  mass );
+        if(fFillSSExtraHisto)
+        {
+          fhMassDispEtaNLocMax2Ebin[ebin]->Fill(dispEta,  mass );
+          fhMassDispPhiNLocMax2Ebin[ebin]->Fill(dispPhi,  mass );
+          fhMassDispAsyNLocMax2Ebin[ebin]->Fill(dispAsy,  mass );
+        }
+      }   
+    }
+    else if(nMax > 2 ) 
+    {
+      if( en > ecut ) 
+      {      
+        fhMassM02NLocMaxN    [0][matched]->Fill(l0     ,  mass ); 
+        if(fFillSSExtraHisto)
+        {
+          fhMassDispEtaNLocMaxN[0][matched]->Fill(dispEta,  mass ); 
+          fhMassDispPhiNLocMaxN[0][matched]->Fill(dispPhi,  mass ); 
+          fhMassDispAsyNLocMaxN[0][matched]->Fill(dispAsy,  mass ); 
+        }
+        
+        if(IsDataMC()) 
+        {
+          fhMassM02NLocMaxN        [mcindex][matched]->Fill(l0     ,  mass ); 
+          if(fFillMCFractionHisto)
+          {
+            fhMCGenFracNLocMaxN      [mcindex][matched]->Fill(en     ,  efrac ); 
+            fhMCGenSplitEFracNLocMaxN[mcindex][matched]->Fill(en     ,  efracSplit ); 
+            fhMCGenEvsSplitENLocMaxN [mcindex][matched]->Fill(eprim  ,  e1+e2); 
+            fhMCGenEFracvsSplitEFracNLocMaxN[mcindex][matched]->Fill(efrac,  splitFrac ); 
+          }
+          
+          if(!matched && ebin >= 0)
+          {
+            if(fFillMCFractionHisto)
+            {
+              fhM02MCGenFracNLocMaxNEbin [mcindex][ebin]->Fill(efrac  ,  l0     ); 
+              fhMassMCGenFracNLocMaxNEbin[mcindex][ebin]->Fill(efrac  ,  mass   ); 
+            }
+            fhMCAsymM02NLocMaxNMCPi0Ebin        [ebin]->Fill(l0     ,  asymGen);
+          }
+          if(fFillSSExtraHisto)
+          {
+            fhMassDispEtaNLocMaxN[mcindex][matched]->Fill(dispEta,  mass ); 
+            fhMassDispPhiNLocMaxN[mcindex][matched]->Fill(dispPhi,  mass ); 
+            fhMassDispAsyNLocMaxN[mcindex][matched]->Fill(dispAsy,  mass ); 
+          }
+        }
+      }
+      
+      if(!matched && ebin >= 0)
+      {
+        fhMassSplitEFractionNLocMaxNEbin[0][ebin]->Fill(splitFrac,  mass);
+        if(IsDataMC())fhMassSplitEFractionNLocMaxNEbin[mcindex][ebin]->Fill(splitFrac,  mass);
 
-    AliAODCaloCluster *cluster1 = new AliAODCaloCluster(0, 0,NULL,0.,NULL,NULL,1,0);
-    AliAODCaloCluster *cluster2 = new AliAODCaloCluster(1, 0,NULL,0.,NULL,NULL,1,0);
+        fhMassM02NLocMaxNEbin    [ebin]->Fill(l0     ,  mass );
+        if(fFillSSExtraHisto)
+        {
+          fhMassDispEtaNLocMaxNEbin[ebin]->Fill(dispEta,  mass );
+          fhMassDispPhiNLocMaxNEbin[ebin]->Fill(dispPhi,  mass );
+          fhMassDispAsyNLocMaxNEbin[ebin]->Fill(dispAsy,  mass );
+        }
+      }
+    }
     
-    SplitEnergy(absId1,absId2,cluster, cells, cluster1, cluster2, nMax /*absIdList, maxEList,*/);
-
     //---------------------------------------------------------------------
-    // Get mass of pair of clusters
+    // From here only if M02 is large but not too large, fill histograms 
     //---------------------------------------------------------------------
-
-    // First set position of cluster as local maxima position, 
-    // assign splitted energy to calculate momentum
-    
-    //TLorentzVector cellMom1 = GetCellMomentum(absId1, en1, cells);
-    //TLorentzVector cellMom2 = GetCellMomentum(absId2, en2, cells);
-
-    TLorentzVector cellMom1; 
-    TLorentzVector cellMom2;  
     
-    cluster1->GetMomentum(cellMom1,GetVertex(0));
-    cluster2->GetMomentum(cellMom2,GetVertex(0));
+    if( l0 < fM02MinCut || l0 > fM02MaxCut ) continue ;
     
-    Float_t  mass  = (cellMom1+cellMom2).M();
-    Double_t angle = cellMom2.Angle(cellMom1.Vect());
-
+    fhNLocMaxM02Cut[0][matched]->Fill(en,nMax);
+    if(IsDataMC()) fhNLocMaxM02Cut[mcindex][matched]->Fill(en,nMax);
+        
     if     (nMax==1) 
     { 
-      fhAnglePairLocMax1[matched]->Fill(en,angle);
-      fhMassNLocMax1[0][matched] ->Fill(en,mass ); 
-
-      if( en > 7 ) 
-      {      
+      fhMassNLocMax1[0][matched]->Fill(en,mass ); 
+      fhAsymNLocMax1[0][matched]->Fill(en,asym ); 
+      // Effect of cuts in mass histograms 
+      if(splitFrac > 0.85 && !matched)
+      {
+        fhMassSplitECutNLocMax1->Fill(en,mass ); 
+        if(GetCaloPID()->IsInSplitM02Range(en,l0,nMax))
+        {
+          fhMassM02CutNLocMax1->Fill(en,mass);
+          if(TMath::Abs(asym) < 0.8) fhMassAsyCutNLocMax1->Fill(en,mass);
+        }
+      }
+      
+      if(fFillAngleHisto) 
+      {
+        fhAnglePairLocMax1[matched]->Fill(en,angle);
+      if( en > ecut ) 
         fhAnglePairMassLocMax1[matched]->Fill(mass,angle);
-        fhMassM02NLocMax1[0]  [matched]->Fill(l0,  mass ); 
       }
       
-      if(en > 6  && en <= 10) fhMassM02NLocMax1Ebin[0]->Fill(l0,  mass ); 
-      if(en > 10 && en <= 15) fhMassM02NLocMax1Ebin[1]->Fill(l0,  mass ); 
-      if(en > 15 && en <= 20) fhMassM02NLocMax1Ebin[2]->Fill(l0,  mass ); 
-      if(en > 20)             fhMassM02NLocMax1Ebin[3]->Fill(l0,  mass ); 
-
-      if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax1[0][matched]->Fill(en,l0);
-      else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax1[0][matched]->Fill(en,l0);
-      else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax1[0][matched]->Fill(en,l0);
+      if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMax1[0][matched]->Fill(en,l0); fhMassConLocMax1[0][matched]->Fill(en,mass);  fhAsyConLocMax1[0][matched]->Fill(en,asym); }
+      else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMax1[0][matched]->Fill(en,l0); fhMassPi0LocMax1[0][matched]->Fill(en,mass);  fhAsyPi0LocMax1[0][matched]->Fill(en,asym); }
+      else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMax1[0][matched]->Fill(en,l0); fhMassEtaLocMax1[0][matched]->Fill(en,mass);  fhAsyEtaLocMax1[0][matched]->Fill(en,asym); }
     }  
     else if(nMax==2) 
     {
-      fhAnglePairLocMax2[matched]->Fill(en,angle);
-      fhMassNLocMax2[0] [matched]->Fill(en,mass );
-
-      if( en > 7 )
+      fhMassNLocMax2[0][matched]->Fill(en,mass );
+      fhAsymNLocMax2[0][matched]->Fill(en,asym ); 
+      // Effect of cuts in mass histograms 
+      if(splitFrac > 0.85 && !matched)
       {
-        fhAnglePairMassLocMax2[matched]->Fill(mass,angle);
-        fhMassM02NLocMax2[0][matched]  ->Fill(l0,  mass ); 
-
+        fhMassSplitECutNLocMax2->Fill(en,mass); 
+        if(GetCaloPID()->IsInSplitM02Range(en,l0,nMax))
+        {
+          fhMassM02CutNLocMax2->Fill(en,mass);
+          if(TMath::Abs(asym) < 0.8) fhMassAsyCutNLocMax2->Fill(en,mass);
+        }
       }
       
-      if(en > 6  && en <= 10) fhMassM02NLocMax2Ebin[0]->Fill(l0,  mass ); 
-      if(en > 10 && en <= 15) fhMassM02NLocMax2Ebin[1]->Fill(l0,  mass ); 
-      if(en > 15 && en <= 20) fhMassM02NLocMax2Ebin[2]->Fill(l0,  mass ); 
-      if(en > 20)             fhMassM02NLocMax2Ebin[3]->Fill(l0,  mass );       
+      if(fFillAngleHisto) 
+      {
+        fhAnglePairLocMax2[matched]->Fill(en,angle);
+        if( en > ecut ) 
+          fhAnglePairMassLocMax2[matched]->Fill(mass,angle);
+      }
       
-      if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax2[0][matched]->Fill(en,l0);
-      else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax2[0][matched]->Fill(en,l0);
-      else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax2[0][matched]->Fill(en,l0);        
+      if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMax2[0][matched]->Fill(en,l0); fhMassConLocMax2[0][matched]->Fill(en,mass);  fhAsyConLocMax2[0][matched]->Fill(en,asym); }
+      else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMax2[0][matched]->Fill(en,l0); fhMassPi0LocMax2[0][matched]->Fill(en,mass);  fhAsyPi0LocMax2[0][matched]->Fill(en,asym); }
+      else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMax2[0][matched]->Fill(en,l0); fhMassEtaLocMax2[0][matched]->Fill(en,mass);  fhAsyEtaLocMax2[0][matched]->Fill(en,asym); }      
     }
     else if(nMax >2) 
     {
-      fhAnglePairLocMaxN[matched]->Fill(en,angle);
-      fhMassNLocMaxN[0] [matched]->Fill(en,mass );
-
-      if( en > 7 ) 
-      {      
-        fhAnglePairMassLocMaxN[matched]->Fill(mass,angle);
-        fhMassM02NLocMaxN[0]  [matched]->Fill(l0  ,mass );
+      fhMassNLocMaxN[0][matched]->Fill(en,mass);
+      fhAsymNLocMaxN[0][matched]->Fill(en,asym); 
+      // Effect of cuts in mass histograms
+      if(splitFrac > 0.85 && !matched)
+      {
+        fhMassSplitECutNLocMaxN->Fill(en,mass ); 
+        if(GetCaloPID()->IsInSplitM02Range(en,l0,nMax))
+        {
+          fhMassM02CutNLocMaxN->Fill(en,mass);
+          if(TMath::Abs(asym) < 0.8) fhMassAsyCutNLocMaxN->Fill(en,mass);
+        }
       }
       
-      if(en > 6  && en <= 10) fhMassM02NLocMaxNEbin[0]->Fill(l0,  mass ); 
-      if(en > 10 && en <= 15) fhMassM02NLocMaxNEbin[1]->Fill(l0,  mass ); 
-      if(en > 15 && en <= 20) fhMassM02NLocMaxNEbin[2]->Fill(l0,  mass ); 
-      if(en > 20)             fhMassM02NLocMaxNEbin[3]->Fill(l0,  mass );       
+      if(fFillAngleHisto) 
+      {
+        fhAnglePairLocMaxN[matched]->Fill(en,angle);
+        if( en > ecut ) 
+          fhAnglePairMassLocMaxN[matched]->Fill(mass,angle);
+      }
       
-      if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMaxN[0][matched]->Fill(en,l0);
-      else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMaxN[0][matched]->Fill(en,l0);
-      else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMaxN[0][matched]->Fill(en,l0);
+      if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMaxN[0][matched]->Fill(en,l0); fhMassConLocMaxN[0][matched]->Fill(en,mass);  fhAsyConLocMaxN[0][matched]->Fill(en,asym); }
+      else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMaxN[0][matched]->Fill(en,l0); fhMassPi0LocMaxN[0][matched]->Fill(en,mass);  fhAsyPi0LocMaxN[0][matched]->Fill(en,asym); }
+      else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMaxN[0][matched]->Fill(en,l0); fhMassEtaLocMaxN[0][matched]->Fill(en,mass);  fhAsyEtaLocMaxN[0][matched]->Fill(en,asym); } 
     }
     
-    if(IsDataMC()){
-            
+    
+    if(IsDataMC())
+    {
       if     (nMax==1) 
       { 
         fhMassNLocMax1[mcindex][matched]->Fill(en,mass); 
-        if( en > 7 )  fhMassM02NLocMax1[mcindex][matched]->Fill(l0,  mass ); 
-        if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax1[mcindex][matched]->Fill(en,l0);
-        else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax1[mcindex][matched]->Fill(en,l0);
-        else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax1[mcindex][matched]->Fill(en,l0);
+        fhAsymNLocMax1[mcindex][matched]->Fill(en,asym); 
+        if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMax1[mcindex][matched]->Fill(en,l0); fhMassConLocMax1[mcindex][matched]->Fill(en,mass); fhAsyConLocMax1[mcindex][matched]->Fill(en,asym); }
+        else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMax1[mcindex][matched]->Fill(en,l0); fhMassPi0LocMax1[mcindex][matched]->Fill(en,mass); fhAsyPi0LocMax1[mcindex][matched]->Fill(en,asym); }
+        else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMax1[mcindex][matched]->Fill(en,l0); fhMassEtaLocMax1[mcindex][matched]->Fill(en,mass); fhAsyEtaLocMax1[mcindex][matched]->Fill(en,asym); } 
       }  
       else if(nMax==2) 
       {
         fhMassNLocMax2[mcindex][matched]->Fill(en,mass);
-        if( en > 7 )  fhMassM02NLocMax2[mcindex][matched]->Fill(l0,  mass ); 
-        if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMax2[mcindex][matched]->Fill(en,l0);
-        else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMax2[mcindex][matched]->Fill(en,l0);
-        else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMax2[mcindex][matched]->Fill(en,l0);        
+        fhAsymNLocMax2[mcindex][matched]->Fill(en,asym); 
+        if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMax2[mcindex][matched]->Fill(en,l0); fhMassConLocMax2[mcindex][matched]->Fill(en,mass); fhAsyConLocMax2[mcindex][matched]->Fill(en,asym); }
+        else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMax2[mcindex][matched]->Fill(en,l0); fhMassPi0LocMax2[mcindex][matched]->Fill(en,mass); fhAsyPi0LocMax2[mcindex][matched]->Fill(en,asym); }
+        else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMax2[mcindex][matched]->Fill(en,l0); fhMassEtaLocMax2[mcindex][matched]->Fill(en,mass); fhAsyEtaLocMax2[mcindex][matched]->Fill(en,asym); } 
+        
       }
       else if(nMax >2) 
       {
         fhMassNLocMaxN[mcindex][matched]->Fill(en,mass);
-        if( en > 7 )  fhMassM02NLocMaxN[mcindex][matched]->Fill(l0,  mass ); 
-        if     (mass < fMassConMax && mass > fMassConMin) fhM02ConLocMaxN[mcindex][matched]->Fill(en,l0);
-        else if(mass < fMassPi0Max && mass > fMassPi0Min) fhM02Pi0LocMaxN[mcindex][matched]->Fill(en,l0);
-        else if(mass < fMassEtaMax && mass > fMassEtaMin) fhM02EtaLocMaxN[mcindex][matched]->Fill(en,l0);
+        fhAsymNLocMaxN[mcindex][matched]->Fill(en,asym); 
+        if     (pidTag==AliCaloPID::kPhoton) { fhM02ConLocMaxN[mcindex][matched]->Fill(en,l0); fhMassConLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyConLocMaxN[mcindex][matched]->Fill(en,asym); }
+        else if(pidTag==AliCaloPID::kPi0   ) { fhM02Pi0LocMaxN[mcindex][matched]->Fill(en,l0); fhMassPi0LocMaxN[mcindex][matched]->Fill(en,mass); fhAsyPi0LocMaxN[mcindex][matched]->Fill(en,asym); }
+        else if(pidTag==AliCaloPID::kEta)    { fhM02EtaLocMaxN[mcindex][matched]->Fill(en,l0); fhMassEtaLocMaxN[mcindex][matched]->Fill(en,mass); fhAsyEtaLocMaxN[mcindex][matched]->Fill(en,asym); } 
       }
       
     }//Work with MC truth first
     
-    delete cluster1 ;
-    delete cluster2 ;
-    delete [] absIdList ;
-    delete [] maxEList  ;
-
   }//loop
   
   if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::MakeAnalysisFillHistograms() - END \n");  
@@ -948,260 +1738,13 @@ void AliAnaInsideClusterInvariantMass::Print(const Option_t * opt) const
   printf("Calorimeter     =     %s\n",  fCalorimeter.Data()) ;
   printf("Loc. Max. E > %2.2f\n",       GetCaloUtils()->GetLocalMaximaCutE());
   printf("Loc. Max. E Diff > %2.2f\n",  GetCaloUtils()->GetLocalMaximaCutEDiff());
+  printf("Min. N Cells =%d \n",         fMinNCells) ;
+  printf("Min. Dist. to Bad =%1.1f \n", fMinBadDist) ;
   printf("%2.2f < lambda_0^2 <%2.2f \n",fM02MinCut,fM02MaxCut);
-  printf("pi0 : %2.2f<m<%2.2f \n",      fMassPi0Min,fMassPi0Max);
-  printf("eta : %2.2f<m<%2.2f \n",      fMassEtaMin,fMassEtaMax);
-  printf("conv: %2.2f<m<%2.2f \n",      fMassConMin,fMassConMax);
-
   printf("    \n") ;
   
 } 
 
 
 
-//________________________________________________________________________________________
-void AliAnaInsideClusterInvariantMass::SplitEnergy(const Int_t absId1, const Int_t absId2,
-                                                   AliVCluster* cluster, 
-                                                   AliVCaloCells* cells,
-                                                   //Float_t & e1, Float_t & e2,
-                                                   AliAODCaloCluster* cluster1,
-                                                   AliAODCaloCluster* cluster2,
-                                                   const Int_t nMax)
-{
-
-  // Split energy of cluster between the 2 local maxima, sum energy on 3x3, and if the 2 
-  // maxima are too close and have common cells, split the energy between the 2
-  
-  TH2F* hClusterMap    = 0 ;
-  TH2F* hClusterLocMax = 0 ;
-  TH2F* hCluster1      = 0 ;
-  TH2F* hCluster2      = 0 ;
-  
-  if(fPlotCluster)
-  {
-    hClusterMap    = new TH2F("hClusterMap","Cluster Map",48,0,48,24,0,24);
-    hClusterLocMax = new TH2F("hClusterLocMax","Cluster 2 highest local maxima",48,0,48,24,0,24);
-    hCluster1      = new TH2F("hCluster1","Cluster 1",48,0,48,24,0,24);
-    hCluster2      = new TH2F("hCluster2","Cluster 2",48,0,48,24,0,24);
-    hClusterMap    ->SetXTitle("column");
-    hClusterMap    ->SetYTitle("row");
-    hClusterLocMax ->SetXTitle("column");
-    hClusterLocMax ->SetYTitle("row");
-    hCluster1      ->SetXTitle("column");
-    hCluster1      ->SetYTitle("row");
-    hCluster2      ->SetXTitle("column");
-    hCluster2      ->SetYTitle("row");
-  }
-  
-  const Int_t ncells  = cluster->GetNCells();  
-  Int_t absIdList[ncells]; 
-  
-  Float_t e1 = 0,  e2   = 0 ;
-  Int_t icol = -1, irow = -1, iRCU = -1, sm = -1;  
-  Float_t eCluster = 0;
-  Float_t minCol = 100, minRow = 100, maxCol = -1, maxRow = -1; 
-  for(Int_t iDigit  = 0; iDigit < ncells; iDigit++ ) {
-    absIdList[iDigit] = cluster->GetCellsAbsId()[iDigit];
-    
-  
-    Float_t ec = cells->GetCellAmplitude(absIdList[iDigit]);
-    GetCaloUtils()->RecalibrateCellAmplitude(ec,fCalorimeter, absIdList[iDigit]);
-    eCluster+=ec;
-    
-    if(fPlotCluster) 
-    {
-      //printf("iDigit %d, absId %d, Ecell %f\n",iDigit,absIdList[iDigit], cells->GetCellAmplitude(absIdList[iDigit]));
-      sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList[iDigit], fCalorimeter, icol, irow, iRCU) ;
-      if(icol > maxCol) maxCol = icol;
-      if(icol < minCol) minCol = icol;
-      if(irow > maxRow) maxRow = irow;
-      if(irow < minRow) minRow = irow;
-      hClusterMap->Fill(icol,irow,ec);
-    }
-     
-  }
-    
-  // Init counters and variables
-  Int_t ncells1 = 1 ;
-  UShort_t absIdList1[9] ;  
-  Double_t fracList1 [9] ;  
-  absIdList1[0] = absId1 ;
-  fracList1 [0] = 1. ;
-  
-  Float_t ecell1 = cells->GetCellAmplitude(absId1);
-  GetCaloUtils()->RecalibrateCellAmplitude(ecell1, fCalorimeter, absId1);
-  e1 =  ecell1;  
-  
-  Int_t ncells2 = 1 ;
-  UShort_t absIdList2[9] ;  
-  Double_t fracList2 [9] ; 
-  absIdList2[0] = absId2 ;
-  fracList2 [0] = 1. ;
-
-  Float_t ecell2 = cells->GetCellAmplitude(absId2);
-  GetCaloUtils()->RecalibrateCellAmplitude(ecell2, fCalorimeter, absId2);
-  e2 =  ecell2;  
-  
-  if(fPlotCluster)
-  {
-    Int_t icol1 = -1, irow1 = -1, icol2 = -1, irow2 = -1;
-    sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId1, fCalorimeter, icol1, irow1, iRCU) ;
-    hClusterLocMax->Fill(icol1,irow1,ecell1);
-    sm = GetCaloUtils()->GetModuleNumberCellIndexes(absId2, fCalorimeter, icol2, irow2, iRCU) ;
-    hClusterLocMax->Fill(icol2,irow2,ecell2);
-  }
-  
-  // Very rough way to share the cluster energy
-  Float_t eRemain = (eCluster-ecell1-ecell2)/2;
-  Float_t shareFraction1 = ecell1/eCluster+eRemain/eCluster;
-  Float_t shareFraction2 = ecell2/eCluster+eRemain/eCluster;
-  for(Int_t iDigit = 0; iDigit < ncells; iDigit++){
-    Int_t absId = absIdList[iDigit];
-    
-    if(absId==absId1 || absId==absId2 || absId < 0) continue;
-    
-    Float_t ecell = cells->GetCellAmplitude(absId);
-    GetCaloUtils()->RecalibrateCellAmplitude(ecell, fCalorimeter, absId);
-    
-    if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absId ))
-    { 
-       absIdList1[ncells1]= absId;
-    
-      if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId2,absId ))
-      { 
-        fracList1[ncells1] = shareFraction1; 
-        e1 += ecell*shareFraction1;
-      }
-      else 
-      {
-        fracList1[ncells1] = 1.; 
-        e1 += ecell;
-      }
-      
-      ncells1++;
-      
-     } // neigbour to cell1
-    
-    if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId2,absId ))
-    { 
-      absIdList2[ncells2]= absId;
-     
-      if(GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absId ))
-      { 
-        fracList2[ncells2] = shareFraction2; 
-        e2 += ecell*shareFraction2;
-      }
-      else
-      { 
-        fracList2[ncells2] = 1.; 
-        e2 += ecell;
-      }
-
-      ncells2++;
-
-    } // neigbour to cell2
-
-  }
-  
-   if(GetDebug() > 1) printf("AliAnaInsideClusterInvariantMass::SplitEnergy() - n Local Max %d, Cluster energy  = %f, Ecell1 = %f, Ecell2 = %f, Enew1 = %f, Enew2 = %f, Remain %f, \n ncells %d, ncells1 %d, ncells2 %d, f1 %f, f2  %f, sum f12 = %f \n",
-         nMax, eCluster,ecell1,ecell2,e1,e2,eCluster-e1-e2,ncells,ncells1,ncells2,shareFraction1,shareFraction2,shareFraction1+shareFraction2);
-
-  cluster1->SetE(e1);
-  cluster2->SetE(e2);  
-  
-  cluster1->SetNCells(ncells1);
-  cluster2->SetNCells(ncells2);  
-  
-  cluster1->SetCellsAbsId(absIdList1);
-  cluster2->SetCellsAbsId(absIdList2);
-  
-  cluster1->SetCellsAmplitudeFraction(fracList1);
-  cluster2->SetCellsAmplitudeFraction(fracList2);
-  
-  GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster1);
-  GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterPosition(GetEMCALGeometry(), cells, cluster2);
-  
-  
-  if(fPlotCluster)
-  {
-    //printf("Cells of cluster1: ");
-    for(Int_t iDigit  = 0; iDigit < ncells1; iDigit++ ) 
-    {
-      //printf(" %d ",absIdList1[iDigit]);
-      
-      sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList1[iDigit], fCalorimeter, icol, irow, iRCU) ;
-      
-      if( GetCaloUtils()->AreNeighbours(fCalorimeter, absId2,absIdList1[iDigit]) )
-        hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit])*shareFraction1);
-      else 
-        hCluster1->Fill(icol,irow,cells->GetCellAmplitude(absIdList1[iDigit]));
-    }
-    
-    //printf(" \n ");
-    //printf("Cells of cluster2: ");
-    
-    for(Int_t iDigit  = 0; iDigit < ncells2; iDigit++ ) 
-    {
-      //printf(" %d ",absIdList2[iDigit]);
-      
-      sm = GetCaloUtils()->GetModuleNumberCellIndexes(absIdList2[iDigit], fCalorimeter, icol, irow, iRCU) ;
-      if( GetCaloUtils()->AreNeighbours(fCalorimeter, absId1,absIdList2[iDigit]) )
-        hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit])*shareFraction2);
-      else
-        hCluster2->Fill(icol,irow,cells->GetCellAmplitude(absIdList2[iDigit]));
-      
-    }
-    //printf(" \n ");
-    
-    gStyle->SetPadRightMargin(0.1);
-    gStyle->SetPadLeftMargin(0.1);
-    gStyle->SetOptStat(0);
-    gStyle->SetOptFit(000000);
-    
-    if(maxCol-minCol > maxRow-minRow)
-    {
-      maxRow+= maxCol-minCol;
-    }
-    else 
-    {
-      maxCol+= maxRow-minRow;
-    }
-
-    TCanvas  * c= new TCanvas("canvas", "canvas", 4000, 4000) ;
-    c->Divide(2,2);  
-    c->cd(1);
-    gPad->SetGridy();
-    gPad->SetGridx();
-    hClusterMap    ->SetAxisRange(minCol, maxCol,"X");
-    hClusterMap    ->SetAxisRange(minRow, maxRow,"Y");
-    hClusterMap    ->Draw("colz");
-    c->cd(2);
-    gPad->SetGridy();
-    gPad->SetGridx();
-    hClusterLocMax ->SetAxisRange(minCol, maxCol,"X");
-    hClusterLocMax ->SetAxisRange(minRow, maxRow,"Y");
-    hClusterLocMax ->Draw("colz");
-    c->cd(3);
-    gPad->SetGridy();
-    gPad->SetGridx();
-    hCluster1      ->SetAxisRange(minCol, maxCol,"X");
-    hCluster1      ->SetAxisRange(minRow, maxRow,"Y");
-    hCluster1      ->Draw("colz");
-    c->cd(4);
-    gPad->SetGridy();
-    gPad->SetGridx();
-    hCluster2      ->SetAxisRange(minCol, maxCol,"X");
-    hCluster2      ->SetAxisRange(minRow, maxRow,"Y");
-    hCluster2      ->Draw("colz");
-    
-    if(eCluster > 6 )c->Print(Form("clusterFigures/Event%d_E%1.0f_nMax%d_NCell1_%d_NCell2_%d.eps",GetEventNumber(),cluster->E(),nMax,ncells1,ncells2));
-    
-    delete c;
-    delete hClusterMap;
-    delete hClusterLocMax;
-    delete hCluster1;
-    delete hCluster2;
-  }
-}
-