]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx
change wrong histogram name setting
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaPi0EbE.cxx
index fad285cdfdae493b0466a6aa466020bd947715cc..f691f3b44edec04a899c0e22ef39b380384052d3 100755 (executable)
@@ -48,66 +48,80 @@ ClassImp(AliAnaPi0EbE)
 
 //____________________________
 AliAnaPi0EbE::AliAnaPi0EbE() :
-AliAnaCaloTrackCorrBaseClass(),fAnaType(kIMCalo),            fCalorimeter(""),
-fMinDist(0.),fMinDist2(0.),    fMinDist3(0.),
-fNLMCutMin(-1),                fNLMCutMax(10),
-fTimeCutMin(-10000),           fTimeCutMax(10000),
+AliAnaCaloTrackCorrBaseClass(),     fAnaType(kIMCalo),                  fCalorimeter(""),
+fMinDist(0.),fMinDist2(0.),         fMinDist3(0.),
+fNLMCutMin(-1),                     fNLMCutMax(10),
+fTimeCutMin(-10000),                fTimeCutMax(10000),
 fRejectTrackMatch(kTRUE),
 fFillPileUpHistograms(0),
-fFillWeightHistograms(kFALSE), fFillTMHisto(0),
-fFillSelectClHisto(0),         fFillOnlySimpleSSHisto(1),    fFillEMCALBCHistograms(0),
+fFillWeightHistograms(kFALSE),      fFillTMHisto(0),
+fFillSelectClHisto(0),              fFillOnlySimpleSSHisto(1),          fFillEMCALBCHistograms(0),
 fInputAODGammaConvName(""),
+fCheckSplitDistToBad(0),
 // Histograms
-fhPt(0),                       fhE(0),
-fhEEta(0),                     fhEPhi(0),
-fhPtEta(0),                    fhPtPhi(0),                   fhEtaPhi(0),
-fhEtaPhiEMCALBC0(0),           fhEtaPhiEMCALBC1(0),          fhEtaPhiEMCALBCN(0),
+fhPt(0),                            fhE(0),
+fhPtEta(0),                         fhPtPhi(0),                         fhEtaPhi(0),
+fhEtaPhiEMCALBC0(0),                fhEtaPhiEMCALBC1(0),                fhEtaPhiEMCALBCN(0),
 fhTimeTriggerEMCALBC0UMReMatchOpenTime(0),
 fhTimeTriggerEMCALBC0UMReMatchCheckNeigh(0),
 fhTimeTriggerEMCALBC0UMReMatchBoth(0),
-fhPtCentrality(),              fhPtEventPlane(0),
-fhPtReject(0),                 fhEReject(0),
-fhEEtaReject(0),               fhEPhiReject(0),              fhEtaPhiReject(0),
-fhMass(0),                     fhMassPt(0),                  fhMassSplitPt(0),
-fhSelectedMass(0),             fhSelectedMassPt(0),          fhSelectedMassSplitPt(0),
-fhAsymmetry(0),                fhSelectedAsymmetry(0),
-fhSplitE(0),                   fhSplitPt(0),
-fhSplitPtEta(0),               fhSplitPtPhi(0),
+fhPtCentrality(),                   fhPtEventPlane(0),
+fhPtReject(0),                      fhEReject(0),
+fhPtEtaReject(0),                   fhPtPhiReject(0),                   fhEtaPhiReject(0),
+fhMass(0),                          fhMassPt(0),                        fhMassSplitPt(0),
+fhSelectedMass(),                   fhSelectedMassPt(0),                fhSelectedMassSplitPt(0),
+fhMassNoOverlap(0),                 fhMassPtNoOverlap(0),               fhMassSplitPtNoOverlap(0),
+fhSelectedMassNoOverlap(0),         fhSelectedMassPtNoOverlap(0),       fhSelectedMassSplitPtNoOverlap(0),
+fhMCPi0PtRecoPtPrim(0),                       fhMCEtaPtRecoPtPrim(0),
+fhMCPi0PtRecoPtPrimNoOverlap(0),              fhMCEtaPtRecoPtPrimNoOverlap(0),
+fhMCPi0SplitPtRecoPtPrim(0),                  fhMCEtaSplitPtRecoPtPrim(0),
+fhMCPi0SplitPtRecoPtPrimNoOverlap(0),         fhMCEtaSplitPtRecoPtPrimNoOverlap(0),
+fhMCPi0SelectedPtRecoPtPrim(0),               fhMCEtaSelectedPtRecoPtPrim(0),
+fhMCPi0SelectedPtRecoPtPrimNoOverlap(0),      fhMCEtaSelectedPtRecoPtPrimNoOverlap(0),
+fhMCPi0SelectedSplitPtRecoPtPrim(0),          fhMCEtaSelectedSplitPtRecoPtPrim(0),
+fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap(0), fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap(0),
+fhAsymmetry(0),                     fhSelectedAsymmetry(0),
+fhSplitE(0),                        fhSplitPt(0),
+fhSplitPtEta(0),                    fhSplitPtPhi(0),
 fhNLocMaxSplitPt(0),
-fhPtDecay(0),                  fhEDecay(0),
+fhPtDecay(0),                       fhEDecay(0),
 // Shower shape histos
-fhEDispersion(0),              fhELambda0(0),                fhELambda1(0),
-fhELambda0NoTRD(0),            fhELambda0FracMaxCellCut(0),
-fhEFracMaxCell(0),             fhEFracMaxCellNoTRD(0),
-fhENCells(0),                  fhETime(0),                   fhEPairDiffTime(0),
-fhDispEtaE(0),                 fhDispPhiE(0),
-fhSumEtaE(0),                  fhSumPhiE(0),                 fhSumEtaPhiE(0),
-fhDispEtaPhiDiffE(0),          fhSphericityE(0),
+fhPtDispersion(0),                  fhPtLambda0(0),                     fhPtLambda1(0),
+fhPtLambda0NoTRD(0),                fhPtLambda0FracMaxCellCut(0),
+fhPtFracMaxCell(0),                 fhPtFracMaxCellNoTRD(0),
+fhPtNCells(0),                      fhPtTime(0),                        fhEPairDiffTime(0),
+fhPtDispEta(0),                     fhPtDispPhi(0),
+fhPtSumEta(0),                      fhPtSumPhi(0),                      fhPtSumEtaPhi(0),
+fhPtDispEtaPhiDiff(0),              fhPtSphericity(0),
 
 // MC histos
-fhMCE(),                       fhMCPt(),
-fhMCPhi(),                     fhMCEta(),
-fhMCEReject(),                 fhMCPtReject(),
+fhMCE(),                            fhMCPt(),
+fhMCPtPhi(),                        fhMCPtEta(),
+fhMCEReject(),                      fhMCPtReject(),
 fhMCPtCentrality(),
-fhMCPi0PtGenRecoFraction(0),   fhMCEtaPtGenRecoFraction(0),
-fhMCPi0DecayPt(0),             fhMCPi0DecayPtFraction(0),
-fhMCEtaDecayPt(0),             fhMCEtaDecayPtFraction(0),
+fhMCPi0PtGenRecoFraction(0),        fhMCEtaPtGenRecoFraction(0),
+fhMCPi0DecayPt(0),                  fhMCPi0DecayPtFraction(0),
+fhMCEtaDecayPt(0),                  fhMCEtaDecayPtFraction(0),
 fhMCOtherDecayPt(0),
-fhMassPairMCPi0(0),            fhMassPairMCEta(0),
-fhAnglePairMCPi0(0),           fhAnglePairMCEta(0),
+fhMassPairMCPi0(0),                 fhMassPairMCEta(0),
+fhAnglePairMCPi0(0),                fhAnglePairMCEta(0),
+fhMCPi0PtOrigin(0x0),               fhMCEtaPtOrigin(0x0),
+fhMCPi0ProdVertex(0),               fhMCEtaProdVertex(0),
+fhMCPi0ProdVertexInner(0),          fhMCEtaProdVertexInner(0),
+
 // Weight studies
-fhECellClusterRatio(0),        fhECellClusterLogRatio(0),
-fhEMaxCellClusterRatio(0),     fhEMaxCellClusterLogRatio(0),
-fhTrackMatchedDEta(0),         fhTrackMatchedDPhi(0),        fhTrackMatchedDEtaDPhi(0),
-fhTrackMatchedDEtaPos(0),      fhTrackMatchedDPhiPos(0),     fhTrackMatchedDEtaDPhiPos(0),
-fhTrackMatchedDEtaNeg(0),      fhTrackMatchedDPhiNeg(0),     fhTrackMatchedDEtaDPhiNeg(0),
-fhTrackMatchedMCParticleE(0),
-fhTrackMatchedMCParticleDEta(0), fhTrackMatchedMCParticleDPhi(0),
-fhdEdx(0),                     fhEOverP(0),                 fhEOverPNoTRD(0),
+fhECellClusterRatio(0),             fhECellClusterLogRatio(0),
+fhEMaxCellClusterRatio(0),          fhEMaxCellClusterLogRatio(0),
+fhTrackMatchedDEta(0),              fhTrackMatchedDPhi(0),              fhTrackMatchedDEtaDPhi(0),
+fhTrackMatchedDEtaPos(0),           fhTrackMatchedDPhiPos(0),           fhTrackMatchedDEtaDPhiPos(0),
+fhTrackMatchedDEtaNeg(0),           fhTrackMatchedDPhiNeg(0),           fhTrackMatchedDEtaDPhiNeg(0),
+fhTrackMatchedMCParticlePt(0),
+fhTrackMatchedMCParticleDEta(0),    fhTrackMatchedMCParticleDPhi(0),
+fhdEdx(0),                          fhEOverP(0),                        fhEOverPNoTRD(0),
 // Number of local maxima in cluster
-fhNLocMaxE(0),                 fhNLocMaxPt(0),
+fhNLocMaxPt(0),                     fhNLocMaxPtReject(0),
 // PileUp
-fhTimePtNoCut(0),                    fhTimePtSPD(0),           fhTimePtSPDMulti(0),
+fhTimePtNoCut(0),                   fhTimePtSPD(0),                     fhTimePtSPDMulti(0),
 fhTimeNPileUpVertSPD(0),            fhTimeNPileUpVertTrack(0),
 fhTimeNPileUpVertContributors(0),
 fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0),
@@ -121,36 +135,43 @@ fhPtNPileUpSPDVtxTimeCut2(0),       fhPtNPileUpTrkVtxTimeCut2(0)
   {
     fhMCE              [i] = 0;
     fhMCPt             [i] = 0;
-    fhMCNLocMaxPt      [i] = 0;
-    fhMCPhi            [i] = 0;
-    fhMCEta            [i] = 0;
+    fhMCPtPhi          [i] = 0;
+    fhMCPtEta          [i] = 0;
     fhMCPtCentrality   [i] = 0;
     
     fhMCSplitE         [i] = 0;
     fhMCSplitPt        [i] = 0;
     fhMCSplitPtPhi     [i] = 0;
     fhMCSplitPtEta     [i] = 0;
-    fhMCNLocMaxSplitPt [i] = 0;
-    
-    fhEMCLambda0       [i] = 0;
-    fhEMCLambda0NoTRD  [i] = 0;
-    fhEMCLambda0FracMaxCellCut[i]= 0;
-    fhEMCFracMaxCell   [i] = 0;
-    fhEMCLambda1       [i] = 0;
-    fhEMCDispersion    [i] = 0;
     
-    fhMCEDispEta       [i] = 0;
-    fhMCEDispPhi       [i] = 0;
-    fhMCESumEtaPhi     [i] = 0;
-    fhMCEDispEtaPhiDiff[i] = 0;
-    fhMCESphericity    [i] = 0;
-    fhMCEAsymmetry     [i] = 0;
+    fhMCNLocMaxPt      [i] = 0;
+    fhMCNLocMaxSplitPt [i] = 0;
+    fhMCNLocMaxPtReject[i] = 0;
+    
+    fhMCPtLambda0       [i] = 0;
+    fhMCPtLambda0NoTRD  [i] = 0;
+    fhMCPtLambda0FracMaxCellCut[i]= 0;
+    fhMCPtFracMaxCell   [i] = 0;
+    fhMCPtLambda1       [i] = 0;
+    fhMCPtDispersion    [i] = 0;
+    
+    fhMCPtDispEta       [i] = 0;
+    fhMCPtDispPhi       [i] = 0;
+    fhMCPtSumEtaPhi     [i] = 0;
+    fhMCPtDispEtaPhiDiff[i] = 0;
+    fhMCPtSphericity    [i] = 0;
+    fhMCPtAsymmetry     [i] = 0;
     
     fhMCMassPt             [i]=0;
     fhMCMassSplitPt        [i]=0;
     fhMCSelectedMassPt     [i]=0;
     fhMCSelectedMassSplitPt[i]=0;
     
+    fhMCMassPtNoOverlap             [i]=0;
+    fhMCMassSplitPtNoOverlap        [i]=0;
+    fhMCSelectedMassPtNoOverlap     [i]=0;
+    fhMCSelectedMassSplitPtNoOverlap[i]=0;
+    
     for(Int_t j = 0; j < 7; j++)
     {
       fhMCLambda0DispEta    [j][i] = 0;
@@ -176,15 +197,33 @@ fhPtNPileUpSPDVtxTimeCut2(0),       fhPtNPileUpTrkVtxTimeCut2(0)
   
   for(Int_t i = 0; i < 3; i++)
   {
-    fhELambda0LocMax       [i] = 0;
-    fhELambda1LocMax       [i] = 0;
-    fhEDispersionLocMax    [i] = 0;
-    fhEDispEtaLocMax       [i] = 0;
-    fhEDispPhiLocMax       [i] = 0;
-    fhESumEtaPhiLocMax     [i] = 0;
-    fhEDispEtaPhiDiffLocMax[i] = 0;
-    fhESphericityLocMax    [i] = 0;
-    fhEAsymmetryLocMax     [i] = 0;
+    fhPtLambda0LocMax       [i] = 0;
+    fhPtLambda1LocMax       [i] = 0;
+    fhPtDispersionLocMax    [i] = 0;
+    fhPtDispEtaLocMax       [i] = 0;
+    fhPtDispPhiLocMax       [i] = 0;
+    fhPtSumEtaPhiLocMax     [i] = 0;
+    fhPtDispEtaPhiDiffLocMax[i] = 0;
+    fhPtSphericityLocMax    [i] = 0;
+    fhPtAsymmetryLocMax     [i] = 0;
+    fhMassPtLocMax          [i] = 0;
+    fhSelectedMassPtLocMax  [i] = 0;
+    for(Int_t ipart = 0; ipart<6; ipart++)
+    {
+      fhMCPtLambda0LocMax     [ipart][i] = 0;
+      fhMCSelectedMassPtLocMax[ipart][i] = 0;
+    }
+    
+    fhMCPi0PtRecoPtPrimLocMax             [i] = 0;
+    fhMCEtaPtRecoPtPrimLocMax             [i] = 0;
+    fhMCPi0SplitPtRecoPtPrimLocMax        [i] = 0;
+    fhMCEtaSplitPtRecoPtPrimLocMax        [i] = 0;
+
+    fhMCPi0SelectedPtRecoPtPrimLocMax     [i] = 0;
+    fhMCEtaSelectedPtRecoPtPrimLocMax     [i] = 0;
+    fhMCPi0SelectedSplitPtRecoPtPrimLocMax[i] = 0;
+    fhMCEtaSelectedSplitPtRecoPtPrimLocMax[i] = 0;
+
   }
   
   //Weight studies
@@ -205,13 +244,22 @@ fhPtNPileUpSPDVtxTimeCut2(0),       fhPtNPileUpTrkVtxTimeCut2(0)
     
   }
   
+  for(Int_t iSM = 0; iSM < 22; iSM++)
+  {
+    fhNLocMaxPtSM[iSM] = 0;
+    for(Int_t inlm = 0; inlm < 3; inlm++)
+    {
+      fhSelectedMassPtLocMaxSM    [inlm][iSM] = 0;
+      fhSelectedLambda0PtLocMaxSM [inlm][iSM] = 0;
+    }
+  }
   //Initialize parameters
   InitParameters();
   
 }
 
-//_______________________________________________________________________________________________
-void AliAnaPi0EbE::FillPileUpHistograms(const Float_t pt, const Float_t time, AliVCluster * calo)
+//___________________________________________________________________________________
+void AliAnaPi0EbE::FillPileUpHistograms(Float_t pt, Float_t time, AliVCluster * calo)
 {
   // Fill some histograms to understand pile-up
   if(!fFillPileUpHistograms) return;
@@ -384,8 +432,8 @@ void AliAnaPi0EbE::FillPileUpHistograms(const Float_t pt, const Float_t time, Al
 }
 
 
-//___________________________________________________________________________________________
-void AliAnaPi0EbE::FillRejectedClusterHistograms(const TLorentzVector mom, const Int_t mctag)
+//______________________________________________________________________________________________
+void AliAnaPi0EbE::FillRejectedClusterHistograms(TLorentzVector mom, Int_t mctag, Int_t nMaxima)
 {
   // Fill histograms that do not pass the identification (SS case only)
   
@@ -398,40 +446,41 @@ void AliAnaPi0EbE::FillRejectedClusterHistograms(const TLorentzVector mom, const
   fhPtReject     ->Fill(pt);
   fhEReject      ->Fill(ener);
   
-  fhEEtaReject   ->Fill(ener,eta);
-  fhEPhiReject   ->Fill(ener,phi);
+  fhPtEtaReject  ->Fill(ener,eta);
+  fhPtPhiReject  ->Fill(ener,phi);
   fhEtaPhiReject ->Fill(eta,phi);
   
+  fhNLocMaxPtReject->Fill(pt,nMaxima);
+
   if(IsDataMC())
   {
     Int_t mcIndex = GetMCIndex(mctag);
     fhMCEReject  [mcIndex] ->Fill(ener);
     fhMCPtReject [mcIndex] ->Fill(pt);
+    fhMCNLocMaxPtReject[mcIndex]->Fill(pt,nMaxima);
   }
 }
 
-//_____________________________________________________________________________________
-void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
-                                                 const Int_t nMaxima,
-                                                 const Int_t tag,
-                                                 const Float_t asy)
+//___________________________________________________________________________________
+void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, Float_t pt, Int_t nMaxima,
+                                                 Int_t tag, Float_t asy)
 {
   // Fill shower shape, timing and other histograms for selected clusters from decay
   
-  Float_t e    = cluster->E();
+  Float_t ener = cluster->E();
   Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
   Float_t l0   = cluster->GetM02();
   Float_t l1   = cluster->GetM20();
   Int_t   nSM  = GetModuleNumber(cluster);
   
-  Int_t ebin = -1;
-  if      (e < 2 ) ebin = 0;
-  else if (e < 4 ) ebin = 1;
-  else if (e < 6 ) ebin = 2;
-  else if (e < 10) ebin = 3;
-  else if (e < 15) ebin = 4;
-  else if (e < 20) ebin = 5;
-  else             ebin = 6;
+  Int_t ptbin = -1;
+  if      (pt < 2 ) ptbin = 0;
+  else if (pt < 4 ) ptbin = 1;
+  else if (pt < 6 ) ptbin = 2;
+  else if (pt < 10) ptbin = 3;
+  else if (pt < 15) ptbin = 4;
+  else if (pt < 20) ptbin = 5;
+  else              ptbin = 6;
   
   Int_t indexMax = -1;
   if     (nMaxima==1) indexMax = 0 ;
@@ -447,13 +496,13 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
   
   Float_t maxCellFraction = 0;
   GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
-  fhEFracMaxCell->Fill(e,maxCellFraction);
+  fhPtFracMaxCell->Fill(pt,maxCellFraction);
   
   FillWeightHistograms(cluster);
   
-  fhEDispersion->Fill(e, disp);
-  fhELambda0   ->Fill(e, l0  );
-  fhELambda1   ->Fill(e, l1  );
+  fhPtDispersion->Fill(pt, disp);
+  fhPtLambda0   ->Fill(pt, l0  );
+  fhPtLambda1   ->Fill(pt, l1  );
   
   Float_t ll0  = 0., ll1  = 0.;
   Float_t dispp= 0., dEta = 0., dPhi    = 0.;
@@ -463,55 +512,58 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
                                                                                  ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
     
-    fhDispEtaE        -> Fill(e,dEta);
-    fhDispPhiE        -> Fill(e,dPhi);
-    fhSumEtaE         -> Fill(e,sEta);
-    fhSumPhiE         -> Fill(e,sPhi);
-    fhSumEtaPhiE      -> Fill(e,sEtaPhi);
-    fhDispEtaPhiDiffE -> Fill(e,dPhi-dEta);
-    if(dEta+dPhi>0)fhSphericityE -> Fill(e,(dPhi-dEta)/(dEta+dPhi));
+    fhPtDispEta       -> Fill(pt,dEta);
+    fhPtDispPhi       -> Fill(pt,dPhi);
+    fhPtSumEta        -> Fill(pt,sEta);
+    fhPtSumPhi        -> Fill(pt,sPhi);
+    fhPtSumEtaPhi     -> Fill(pt,sEtaPhi);
+    fhPtDispEtaPhiDiff-> Fill(pt,dPhi-dEta);
+    if(dEta+dPhi>0)fhPtSphericity-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
     
-    fhDispEtaDispPhi[ebin]->Fill(dEta,dPhi);
-    fhLambda0DispEta[ebin]->Fill(l0  ,dEta);
-    fhLambda0DispPhi[ebin]->Fill(l0  ,dPhi);
+    fhDispEtaDispPhi[ptbin]->Fill(dEta,dPhi);
+    fhLambda0DispEta[ptbin]->Fill(l0  ,dEta);
+    fhLambda0DispPhi[ptbin]->Fill(l0  ,dPhi);
     
     if (fAnaType==kSSCalo)
     {
       // Asymmetry histograms
-      fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
-      fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
-      fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
+      fhAsymmetryLambda0[ptbin]->Fill(l0 ,asy);
+      fhAsymmetryDispEta[ptbin]->Fill(dEta,asy);
+      fhAsymmetryDispPhi[ptbin]->Fill(dPhi,asy);
     }
   }
   
-  fhNLocMaxE ->Fill(e ,nMaxima);
+  fhNLocMaxPt->Fill(pt,nMaxima);
   
-  fhELambda0LocMax   [indexMax]->Fill(e,l0);
-  fhELambda1LocMax   [indexMax]->Fill(e,l1);
-  fhEDispersionLocMax[indexMax]->Fill(e,disp);
+  if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
+    fhNLocMaxPtSM[nSM]->Fill(pt,nMaxima);
   
+  fhPtLambda0LocMax   [indexMax]->Fill(pt,l0);
+  fhPtLambda1LocMax   [indexMax]->Fill(pt,l1);
+  fhPtDispersionLocMax[indexMax]->Fill(pt,disp);
+
   if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
   {
-    fhEDispEtaLocMax       [indexMax]-> Fill(e,dEta);
-    fhEDispPhiLocMax       [indexMax]-> Fill(e,dPhi);
-    fhESumEtaPhiLocMax     [indexMax]-> Fill(e,sEtaPhi);
-    fhEDispEtaPhiDiffLocMax[indexMax]-> Fill(e,dPhi-dEta);
-    if(dEta+dPhi>0)       fhESphericityLocMax[indexMax]->Fill(e,(dPhi-dEta)/(dEta+dPhi));
-    if(fAnaType==kSSCalo) fhEAsymmetryLocMax [indexMax]->Fill(e  ,asy);
+    fhPtDispEtaLocMax       [indexMax]-> Fill(pt,dEta);
+    fhPtDispPhiLocMax       [indexMax]-> Fill(pt,dPhi);
+    fhPtSumEtaPhiLocMax     [indexMax]-> Fill(pt,sEtaPhi);
+    fhPtDispEtaPhiDiffLocMax[indexMax]-> Fill(pt,dPhi-dEta);
+    if(dEta+dPhi>0)       fhPtSphericityLocMax[indexMax]->Fill(pt,(dPhi-dEta)/(dEta+dPhi));
+    if(fAnaType==kSSCalo) fhPtAsymmetryLocMax [indexMax]->Fill(pt  ,asy);
     
   }
   
-  if(fCalorimeter=="EMCAL" && nSM < 6)
+  if(fCalorimeter=="EMCAL" && nSM < 6) // CAREFUL FOR 2012-13 runs change 6 to 4, -1 for 2015 ...
   {
-    fhELambda0NoTRD->Fill(e, l0  );
-    fhEFracMaxCellNoTRD->Fill(e,maxCellFraction);
+    fhPtLambda0NoTRD    ->Fill(pt, l0  );
+    fhPtFracMaxCellNoTRD->Fill(pt,maxCellFraction);
   }
   
   if(maxCellFraction < 0.5)
-    fhELambda0FracMaxCellCut->Fill(e, l0  );
+    fhPtLambda0FracMaxCellCut->Fill(pt, l0  );
   
-  fhETime  ->Fill(e, cluster->GetTOF()*1.e9);
-  fhENCells->Fill(e, cluster->GetNCells());
+  fhPtTime  ->Fill(pt, cluster->GetTOF()*1.e9);
+  fhPtNCells->Fill(pt, cluster->GetNCells());
   
   // Fill Track matching control histograms
   if(fFillTMHisto)
@@ -533,23 +585,23 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
     
     if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
     {
-      fhTrackMatchedDEta->Fill(e,dZ);
-      fhTrackMatchedDPhi->Fill(e,dR);
-      if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
+      fhTrackMatchedDEta->Fill(pt,dZ);
+      fhTrackMatchedDPhi->Fill(pt,dR);
+      if(ener > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
       
       if(track)
       {
         if(positive)
         {
-          fhTrackMatchedDEtaPos->Fill(cluster->E(),dZ);
-          fhTrackMatchedDPhiPos->Fill(cluster->E(),dR);
-          if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(dZ,dR);
+          fhTrackMatchedDEtaPos->Fill(pt,dZ);
+          fhTrackMatchedDPhiPos->Fill(pt,dR);
+          if(ener > 0.5) fhTrackMatchedDEtaDPhiPos->Fill(dZ,dR);
         }
         else
         {
-          fhTrackMatchedDEtaNeg->Fill(cluster->E(),dZ);
-          fhTrackMatchedDPhiNeg->Fill(cluster->E(),dR);
-          if(cluster->E() > 0.5) fhTrackMatchedDEtaDPhiNeg->Fill(dZ,dR);
+          fhTrackMatchedDEtaNeg->Fill(pt,dZ);
+          fhTrackMatchedDPhiNeg->Fill(pt,dR);
+          if(ener > 0.5) fhTrackMatchedDEtaDPhiNeg->Fill(dZ,dR);
         }
     }
     }
@@ -560,12 +612,13 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
       if(track)
       {
         Float_t dEdx = track->GetTPCsignal();
-        fhdEdx->Fill(e, dEdx);
+        fhdEdx->Fill(pt, dEdx);
         
-        Float_t eOverp = e/track->P();
-        fhEOverP->Fill(e,  eOverp);
+        Float_t eOverp = cluster->E()/track->P();
+        fhEOverP->Fill(pt,  eOverp);
         
-        if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(e,  eOverp);
+        // Change nSM for year > 2011 (< 4 in 2012-13, none after)
+        if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(pt,  eOverp);
         
       }
       //else
@@ -577,7 +630,7 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
         {
           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
-                    GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) mctag =  2.5 ;
+                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) mctag =  2.5 ;
           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) mctag =  0.5 ;
           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) mctag =  1.5 ;
           else                                                                                 mctag =  3.5 ;
@@ -586,13 +639,13 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
         else
         {
           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
-                    GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) mctag =  6.5 ;
+                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) mctag =  6.5 ;
           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) mctag =  4.5 ;
           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) mctag =  5.5 ;
           else                                                                                 mctag =  7.5 ;
         }
         
-        fhTrackMatchedMCParticleE   ->Fill(e , mctag);
+        fhTrackMatchedMCParticlePt   ->Fill(pt, mctag);
         fhTrackMatchedMCParticleDEta->Fill(dZ, mctag);
         fhTrackMatchedMCParticleDPhi->Fill(dR, mctag);
         
@@ -604,35 +657,38 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
   {
     Int_t mcIndex = GetMCIndex(tag);
     
-    fhEMCLambda0[mcIndex]    ->Fill(e, l0);
-    fhEMCLambda1[mcIndex]    ->Fill(e, l1);
-    fhEMCDispersion[mcIndex] ->Fill(e, disp);
-    fhEMCFracMaxCell[mcIndex]->Fill(e,maxCellFraction);
+    fhMCPtLambda0[mcIndex]    ->Fill(pt, l0);
+    fhMCPtLambda1[mcIndex]    ->Fill(pt, l1);
+    fhMCPtDispersion[mcIndex] ->Fill(pt, disp);
+    fhMCPtFracMaxCell[mcIndex]->Fill(pt,maxCellFraction);
     
+    fhMCPtLambda0LocMax     [mcIndex][indexMax]->Fill(pt,l0);
+
+    // Change nSM for year > 2011 (< 4 in 2012-13, none after)
     if(fCalorimeter=="EMCAL" && nSM < 6)
-      fhEMCLambda0NoTRD[mcIndex]->Fill(e, l0  );
+      fhMCPtLambda0NoTRD[mcIndex]->Fill(pt, l0  );
     
     if(maxCellFraction < 0.5)
-      fhEMCLambda0FracMaxCellCut[mcIndex]->Fill(e, l0  );
+      fhMCPtLambda0FracMaxCellCut[mcIndex]->Fill(pt, l0  );
     
     if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
     {
-      fhMCEDispEta        [mcIndex]-> Fill(e,dEta);
-      fhMCEDispPhi        [mcIndex]-> Fill(e,dPhi);
-      fhMCESumEtaPhi      [mcIndex]-> Fill(e,sEtaPhi);
-      fhMCEDispEtaPhiDiff [mcIndex]-> Fill(e,dPhi-dEta);
-      if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(e,(dPhi-dEta)/(dEta+dPhi));
+      fhMCPtDispEta        [mcIndex]-> Fill(pt,dEta);
+      fhMCPtDispPhi        [mcIndex]-> Fill(pt,dPhi);
+      fhMCPtSumEtaPhi      [mcIndex]-> Fill(pt,sEtaPhi);
+      fhMCPtDispEtaPhiDiff [mcIndex]-> Fill(pt,dPhi-dEta);
+      if(dEta+dPhi>0)fhMCPtSphericity[mcIndex]-> Fill(pt,(dPhi-dEta)/(dEta+dPhi));
       
       if (fAnaType==kSSCalo)
       {
-        fhMCAsymmetryLambda0[ebin][mcIndex]->Fill(l0 ,asy);
-        fhMCAsymmetryDispEta[ebin][mcIndex]->Fill(dEta,asy);
-        fhMCAsymmetryDispPhi[ebin][mcIndex]->Fill(dPhi,asy);
+        fhMCAsymmetryLambda0[ptbin][mcIndex]->Fill(l0 ,asy);
+        fhMCAsymmetryDispEta[ptbin][mcIndex]->Fill(dEta,asy);
+        fhMCAsymmetryDispPhi[ptbin][mcIndex]->Fill(dPhi,asy);
       }
       
-      fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta,dPhi);
-      fhMCLambda0DispEta[ebin][mcIndex]->Fill(l0  ,dEta);
-      fhMCLambda0DispPhi[ebin][mcIndex]->Fill(l0  ,dPhi);
+      fhMCDispEtaDispPhi[ptbin][mcIndex]->Fill(dEta,dPhi);
+      fhMCLambda0DispEta[ptbin][mcIndex]->Fill(l0  ,dEta);
+      fhMCLambda0DispPhi[ptbin][mcIndex]->Fill(l0  ,dPhi);
       
     }
     
@@ -785,47 +841,35 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
   
-  Int_t   ntimebins= GetHistogramRanges()->GetHistoTimeBins();
-  Float_t timemax  = GetHistogramRanges()->GetHistoTimeMax();
-  Float_t timemin  = GetHistogramRanges()->GetHistoTimeMin();
+  Int_t   ntimptbins  = GetHistogramRanges()->GetHistoTimeBins();
+  Float_t timemax     = GetHistogramRanges()->GetHistoTimeMax();
+  Float_t timemin     = GetHistogramRanges()->GetHistoTimeMin();
   
-  TString nlm[]   ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
-  TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
-  TString pname[] ={"Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
+  TString nlm[]   = {"1 Local Maxima","2 Local Maxima", "NLM > 2"};
+  TString ptype[] = {"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
+  TString pname[] = {"Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
   Int_t   bin[]   = {0,2,4,6,10,15,20,100}; // energy bins
   
   fhPt  = new TH1F("hPt","Number of identified  #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
-  fhPt->SetYTitle("N");
-  fhPt->SetXTitle("p_{T} (GeV/c)");
+  fhPt->SetYTitle("#it{N}");
+  fhPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
   outputContainer->Add(fhPt) ;
   
   fhE  = new TH1F("hE","Number of identified  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
-  fhE->SetYTitle("N");
-  fhE->SetXTitle("E (GeV)");
+  fhE->SetYTitle("#it{N}");
+  fhE->SetXTitle("#it{E} (GeV)");
   outputContainer->Add(fhE) ;
   
-  fhEPhi  = new TH2F
-  ("hEPhi","Selected #pi^{0} (#eta) pairs: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
-  fhEPhi->SetYTitle("#phi (rad)");
-  fhEPhi->SetXTitle("E (GeV)");
-  outputContainer->Add(fhEPhi) ;
-  
-  fhEEta  = new TH2F
-  ("hEEta","Selected #pi^{0} (#eta) pairs: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
-  fhEEta->SetYTitle("#eta");
-  fhEEta->SetXTitle("E (GeV)");
-  outputContainer->Add(fhEEta) ;
-  
   fhPtPhi  = new TH2F
-  ("hPtPhi","Selected #pi^{0} (#eta) pairs: p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
+  ("hPtPhi","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
   fhPtPhi->SetYTitle("#phi (rad)");
-  fhPtPhi->SetXTitle("p_{T} (GeV/c)");
+  fhPtPhi->SetXTitle("#it{E} (GeV)");
   outputContainer->Add(fhPtPhi) ;
   
   fhPtEta  = new TH2F
-  ("hPtEta","Selected #pi^{0} (#eta) pairs: p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
+  ("hPtEta","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
   fhPtEta->SetYTitle("#eta");
-  fhPtEta->SetXTitle("p_{T} (GeV/c)");
+  fhPtEta->SetXTitle("#it{E} (GeV)");
   outputContainer->Add(fhPtEta) ;
   
   fhEtaPhi  = new TH2F
@@ -837,19 +881,19 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
   {
     fhEtaPhiEMCALBC0  = new TH2F
-    ("hEtaPhiEMCALBC0","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
+    ("hEtaPhiEMCALBC0","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with |#it{t}| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
     fhEtaPhiEMCALBC0->SetYTitle("#phi (rad)");
     fhEtaPhiEMCALBC0->SetXTitle("#eta");
     outputContainer->Add(fhEtaPhiEMCALBC0) ;
     
     fhEtaPhiEMCALBC1  = new TH2F
-    ("hEtaPhiEMCALBC1","cluster,E > 2 GeV, #eta vs #phi, for clusters with 25 < |time| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
+    ("hEtaPhiEMCALBC1","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with 25 < |#it{t}| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
     fhEtaPhiEMCALBC1->SetYTitle("#phi (rad)");
     fhEtaPhiEMCALBC1->SetXTitle("#eta");
     outputContainer->Add(fhEtaPhiEMCALBC1) ;
     
     fhEtaPhiEMCALBCN  = new TH2F
-    ("hEtaPhiEMCALBCN","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
+    ("hEtaPhiEMCALBCN","cluster, #it{E} > 2 GeV, #eta vs #phi, for clusters with |#it{t}| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
     fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
     fhEtaPhiEMCALBCN->SetXTitle("#eta");
     outputContainer->Add(fhEtaPhiEMCALBCN) ;
@@ -858,7 +902,7 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
     {
       fhEtaPhiTriggerEMCALBC[i] = new TH2F
       (Form("hEtaPhiTriggerEMCALBC%d",i-5),
-       Form("meson E > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-5),
+       Form("meson #it{E} > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-5),
        netabins,etamin,etamax,nphibins,phimin,phimax);
       fhEtaPhiTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
       fhEtaPhiTriggerEMCALBC[i]->SetXTitle("#eta");
@@ -866,23 +910,23 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
       
       fhTimeTriggerEMCALBC[i] = new TH2F
       (Form("hTimeTriggerEMCALBC%d",i-5),
-       Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
-       nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
-      fhTimeTriggerEMCALBC[i]->SetXTitle("E (GeV)");
-      fhTimeTriggerEMCALBC[i]->SetYTitle("time (ns)");
+       Form("meson #it{t} vs #it{E}, Trigger EMCAL-BC=%d",i-5),
+       nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
+      fhTimeTriggerEMCALBC[i]->SetXTitle("#it{E} (GeV)");
+      fhTimeTriggerEMCALBC[i]->SetYTitle("#it{t} (ns)");
       outputContainer->Add(fhTimeTriggerEMCALBC[i]);
       
       fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
       (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-5),
-       Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
-       nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
-      fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("E (GeV)");
-      fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("time (ns)");
+       Form("meson #it{t} vs #it{E}, Trigger EMCAL-BC=%d",i-5),
+       nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
+      fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("#it{E} (GeV)");
+      fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("#it{t} (ns)");
       outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
       
       fhEtaPhiTriggerEMCALBCUM[i] = new TH2F
       (Form("hEtaPhiTriggerEMCALBC%d_UnMatch",i-5),
-       Form("meson E > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-5),
+       Form("meson #it{E} > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-5),
        netabins,etamin,etamax,nphibins,phimin,phimax);
       fhEtaPhiTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
       fhEtaPhiTriggerEMCALBCUM[i]->SetXTitle("#eta");
@@ -890,71 +934,71 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
       
       fhTimeTriggerEMCALBCUM[i] = new TH2F
       (Form("hTimeTriggerEMCALBC%d_UnMatch",i-5),
-       Form("meson time vs E, unmatched trigger EMCAL-BC=%d",i-5),
-       nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
-      fhTimeTriggerEMCALBCUM[i]->SetXTitle("E (GeV)");
-      fhTimeTriggerEMCALBCUM[i]->SetYTitle("time (ns)");
+       Form("meson #it{t} vs #it{E}, unmatched trigger EMCAL-BC=%d",i-5),
+       nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
+      fhTimeTriggerEMCALBCUM[i]->SetXTitle("#it{E} (GeV)");
+      fhTimeTriggerEMCALBCUM[i]->SetYTitle("#it{t} (ns)");
       outputContainer->Add(fhTimeTriggerEMCALBCUM[i]);
       
     }
     
     fhTimeTriggerEMCALBC0UMReMatchOpenTime = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_OpenTime",
-                                                      "cluster time vs E of clusters, no match, rematch open time",
-                                                      nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
-    fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("E (GeV)");
-    fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("time (ns)");
+                                                      "cluster #it{t} vs #it{E} of clusters, no match, rematch open time",
+                                                      nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
+    fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetXTitle("#it{E} (GeV)");
+    fhTimeTriggerEMCALBC0UMReMatchOpenTime->SetYTitle("#it{t} (ns)");
     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchOpenTime);
     
     
     fhTimeTriggerEMCALBC0UMReMatchCheckNeigh = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_CheckNeighbours",
-                                                        "cluster time vs E of clusters, no match, rematch with neigbour parches",
-                                                        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
-    fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("E (GeV)");
-    fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("time (ns)");
+                                                        "cluster #it{t} vs #it{E} of clusters, no match, rematch with neigbour parches",
+                                                        nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
+    fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetXTitle("#it{E} (GeV)");
+    fhTimeTriggerEMCALBC0UMReMatchCheckNeigh->SetYTitle("#it{t} (ns)");
     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchCheckNeigh);
     
     fhTimeTriggerEMCALBC0UMReMatchBoth = new TH2F("hTimeTriggerBC0_UnMatch_ReMatch_Both",
-                                                  "cluster time vs E of clusters, no match, rematch open time and neigbour",
-                                                  nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
-    fhTimeTriggerEMCALBC0UMReMatchBoth->SetXTitle("E (GeV)");
-    fhTimeTriggerEMCALBC0UMReMatchBoth->SetYTitle("time (ns)");
+                                                  "cluster #it{t} vs #it{E} of clusters, no match, rematch open time and neigbour",
+                                                  nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
+    fhTimeTriggerEMCALBC0UMReMatchBoth->SetXTitle("#it{E} (GeV)");
+    fhTimeTriggerEMCALBC0UMReMatchBoth->SetYTitle("#it{t} (ns)");
     outputContainer->Add(fhTimeTriggerEMCALBC0UMReMatchBoth);
     
   }
   
-  fhPtCentrality  = new TH2F("hPtCentrality","centrality vs p_{T}",nptbins,ptmin,ptmax, 100,0,100);
+  fhPtCentrality  = new TH2F("hPtCentrality","centrality vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,100);
   fhPtCentrality->SetYTitle("centrality");
-  fhPtCentrality->SetXTitle("p_{T}(GeV/c)");
+  fhPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
   outputContainer->Add(fhPtCentrality) ;
   
-  fhPtEventPlane  = new TH2F("hPtEventPlane","event plane angle vs p_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
+  fhPtEventPlane  = new TH2F("hPtEventPlane","event plane angle vs #it{p}_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
   fhPtEventPlane->SetYTitle("Event plane angle (rad)");
-  fhPtEventPlane->SetXTitle("p_{T} (GeV/c)");
+  fhPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
   outputContainer->Add(fhPtEventPlane) ;
   
   if(fAnaType == kSSCalo)
   {
     fhPtReject  = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
-    fhPtReject->SetYTitle("N");
-    fhPtReject->SetXTitle("p_{T} (GeV/c)");
+    fhPtReject->SetYTitle("#it{N}");
+    fhPtReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhPtReject) ;
     
     fhEReject  = new TH1F("hEReject","Number of rejected as  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
-    fhEReject->SetYTitle("N");
-    fhEReject->SetXTitle("E (GeV)");
+    fhEReject->SetYTitle("#it{N}");
+    fhEReject->SetXTitle("#it{E} (GeV)");
     outputContainer->Add(fhEReject) ;
     
-    fhEPhiReject  = new TH2F
-    ("hEPhiReject","Rejected #pi^{0} (#eta) cluster: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
-    fhEPhiReject->SetYTitle("#phi (rad)");
-    fhEPhiReject->SetXTitle("E (GeV)");
-    outputContainer->Add(fhEPhiReject) ;
+    fhPtPhiReject  = new TH2F
+    ("hPtPhiReject","Rejected #pi^{0} (#eta) cluster: #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
+    fhPtPhiReject->SetYTitle("#phi (rad)");
+    fhPtPhiReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtPhiReject) ;
     
-    fhEEtaReject  = new TH2F
-    ("hEEtaReject","Rejected #pi^{0} (#eta) cluster: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
-    fhEEtaReject->SetYTitle("#eta");
-    fhEEtaReject->SetXTitle("E (GeV)");
-    outputContainer->Add(fhEEtaReject) ;
+    fhPtEtaReject  = new TH2F
+    ("hPtEtaReject","Rejected #pi^{0} (#eta) cluster: #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
+    fhPtEtaReject->SetYTitle("#eta");
+    fhPtEtaReject->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtEtaReject) ;
     
     fhEtaPhiReject  = new TH2F
     ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
@@ -964,39 +1008,114 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   }
   
   fhMass  = new TH2F
-  ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
-  fhMass->SetYTitle("mass (GeV/c^{2})");
-  fhMass->SetXTitle("E (GeV)");
+  ("hMass","all pairs #it{M}: #it{E} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+  fhMass->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+  fhMass->SetXTitle("#it{E} (GeV)");
   outputContainer->Add(fhMass) ;
   
   fhSelectedMass  = new TH2F
-  ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
-  fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
-  fhSelectedMass->SetXTitle("E (GeV)");
+  ("hSelectedMass","Selected #pi^{0} (#eta) pairs #it{M}: E vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+  fhSelectedMass->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+  fhSelectedMass->SetXTitle("#it{E} (GeV)");
   outputContainer->Add(fhSelectedMass) ;
   
-  fhMassPt  = new TH2F
-  ("hMassPt","all pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
-  fhMassPt->SetYTitle("mass (GeV/c^{2})");
-  fhMassPt->SetXTitle("p_{T} (GeV/c)");
-  outputContainer->Add(fhMassPt) ;
-  
-  fhSelectedMassPt  = new TH2F
-  ("hSelectedMassPt","Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
-  fhSelectedMassPt->SetYTitle("mass (GeV/c^{2})");
-  fhSelectedMassPt->SetXTitle("p_{T} (GeV/c)");
-  outputContainer->Add(fhSelectedMassPt) ;
+  if(fAnaType == kSSCalo)
+  {
+    
+    fhMassPt  = new TH2F
+    ("hMassPt","all pairs #it{M}: #it{p}_{T} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+    fhMassPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+    fhMassPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhMassPt) ;
+    
+    fhSelectedMassPt  = new TH2F
+    ("hSelectedMassPt","Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+    fhSelectedMassPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+    fhSelectedMassPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhSelectedMassPt) ;
+    
+    for(Int_t inlm = 0; inlm < 3; inlm++)
+    {
+      fhMassPtLocMax[inlm]  = new TH2F
+      (Form("hMassPtNLocMax%d",inlm+1),Form("all pairs #it{M}: #it{p}_{T} vs #it{M} and NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+      fhMassPtLocMax[inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+      fhMassPtLocMax[inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhMassPtLocMax[inlm]) ;
+      
+      fhSelectedMassPtLocMax[inlm]  = new TH2F
+      (Form("hSelectedMassPtLocMax%d",inlm+1),Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s",nlm[inlm].Data()),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+      fhSelectedMassPtLocMax[inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+      fhSelectedMassPtLocMax[inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhSelectedMassPtLocMax[inlm]) ;
+      
+      for(Int_t iSM = 0; iSM < GetCaloUtils()->GetNumberOfSuperModulesUsed(); iSM++)
+      {
+        fhSelectedMassPtLocMaxSM[inlm][iSM]  = new TH2F
+        (Form("hSelectedMassPtLocMax%d_SM%d",inlm+1,iSM),Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s for SM=%d",nlm[inlm].Data(),iSM),nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+        fhSelectedMassPtLocMaxSM[inlm][iSM]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+        fhSelectedMassPtLocMaxSM[inlm][iSM]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhSelectedMassPtLocMaxSM[inlm][iSM]) ;
+
+        fhSelectedLambda0PtLocMaxSM[inlm][iSM]  = new TH2F
+        (Form("hSelectedLambda0PtLocMax%d_SM%d",inlm+1,iSM),Form("Selected #pi^{0} (#eta) pairs #lambda_{0}^{2}: #it{p}_{T} vs #it{M}, NLM=%s for SM=%d",nlm[inlm].Data(),iSM),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+        fhSelectedLambda0PtLocMaxSM[inlm][iSM]->SetYTitle("#lambda_{0}^{2}");
+        fhSelectedLambda0PtLocMaxSM[inlm][iSM]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhSelectedLambda0PtLocMaxSM[inlm][iSM]) ;
+      }
+      
+      if(IsDataMC())
+      {
+        for(Int_t ipart = 0; ipart < 6; ipart++)
+        {
+          fhMCSelectedMassPtLocMax[ipart][inlm]  = new TH2F
+          (Form("hSelectedMassPtLocMax%d_MC%s",inlm+1,pname[ipart].Data()),
+           Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, NLM=%s, %s",nlm[inlm].Data(),pname[ipart].Data()),
+           nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+          fhMCSelectedMassPtLocMax[ipart][inlm]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+          fhMCSelectedMassPtLocMax[ipart][inlm]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+          outputContainer->Add(fhMCSelectedMassPtLocMax[ipart][inlm]) ;
+        }
+      }
+    }
+    
+    if(IsDataMC())
+    {
+      fhMassNoOverlap  = new TH2F
+      ("hMassNoOverlap","all pairs #it{M}: #it{E} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+      fhMassNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+      fhMassNoOverlap->SetXTitle("#it{E} (GeV)");
+      outputContainer->Add(fhMassNoOverlap) ;
+      
+      fhSelectedMassNoOverlap  = new TH2F
+      ("hSelectedMassNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: #it{E} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+      fhSelectedMassNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+      fhSelectedMassNoOverlap->SetXTitle("#it{E} (GeV)");
+      outputContainer->Add(fhSelectedMassNoOverlap) ;
+      
+      fhMassPtNoOverlap  = new TH2F
+      ("hMassPtNoOverlap","all pairs #it{M}: #it{p}_{T} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+      fhMassPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+      fhMassPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhMassPtNoOverlap) ;
+      
+      fhSelectedMassPtNoOverlap  = new TH2F
+      ("hSelectedMassPtNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M}, no overlap",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+      fhSelectedMassPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+      fhSelectedMassPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhSelectedMassPtNoOverlap) ;
+    }
+  }
   
   if(fAnaType != kSSCalo)
   {
     fhPtDecay  = new TH1F("hPtDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
-    fhPtDecay->SetYTitle("N");
-    fhPtDecay->SetXTitle("p_{T} (GeV/c)");
+    fhPtDecay->SetYTitle("#it{N}");
+    fhPtDecay->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhPtDecay) ;
     
     fhEDecay  = new TH1F("hEDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
-    fhEDecay->SetYTitle("N");
-    fhEDecay->SetXTitle("E (GeV)");
+    fhEDecay->SetYTitle("#it{N}");
+    fhEDecay->SetXTitle("#it{E} (GeV)");
     outputContainer->Add(fhEDecay) ;
   }
   
@@ -1004,91 +1123,90 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   
   if( fFillSelectClHisto )
   {
-    
-    fhEDispersion  = new TH2F
-    ("hEDispersion","Selected #pi^{0} (#eta) pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-    fhEDispersion->SetYTitle("D^{2}");
-    fhEDispersion->SetXTitle("E (GeV)");
-    outputContainer->Add(fhEDispersion) ;
-    
-    fhELambda0  = new TH2F
-    ("hELambda0","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-    fhELambda0->SetYTitle("#lambda_{0}^{2}");
-    fhELambda0->SetXTitle("E (GeV)");
-    outputContainer->Add(fhELambda0) ;
-    
-    fhELambda1  = new TH2F
-    ("hELambda1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-    fhELambda1->SetYTitle("#lambda_{1}^{2}");
-    fhELambda1->SetXTitle("E (GeV)");
-    outputContainer->Add(fhELambda1) ;
-    
-    fhELambda0FracMaxCellCut  = new TH2F
-    ("hELambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-    fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
-    fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
-    outputContainer->Add(fhELambda0FracMaxCellCut) ;
-    
-    fhEFracMaxCell  = new TH2F
-    ("hEFracMaxCell","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
-    fhEFracMaxCell->SetYTitle("Fraction");
-    fhEFracMaxCell->SetXTitle("E (GeV)");
-    outputContainer->Add(fhEFracMaxCell) ;
+    fhPtDispersion  = new TH2F
+    ("hPtDispersion","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+    fhPtDispersion->SetYTitle("D^{2}");
+    fhPtDispersion->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtDispersion) ;
+    
+    fhPtLambda0  = new TH2F
+    ("hPtLambda0","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+    fhPtLambda0->SetYTitle("#lambda_{0}^{2}");
+    fhPtLambda0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtLambda0) ;
+    
+    fhPtLambda1  = new TH2F
+    ("hPtLambda1","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+    fhPtLambda1->SetYTitle("#lambda_{1}^{2}");
+    fhPtLambda1->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtLambda1) ;
+    
+    fhPtLambda0FracMaxCellCut  = new TH2F
+    ("hPtLambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+    fhPtLambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
+    fhPtLambda0FracMaxCellCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtLambda0FracMaxCellCut) ;
+    
+    fhPtFracMaxCell  = new TH2F
+    ("hPtFracMaxCell","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
+    fhPtFracMaxCell->SetYTitle("Fraction");
+    fhPtFracMaxCell->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhPtFracMaxCell) ;
     
     if(fCalorimeter=="EMCAL")
     {
-      fhELambda0NoTRD  = new TH2F
-      ("hELambda0NoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-      fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
-      fhELambda0NoTRD->SetXTitle("E (GeV)");
-      outputContainer->Add(fhELambda0NoTRD) ;
-      
-      fhEFracMaxCellNoTRD  = new TH2F
-      ("hEFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
-      fhEFracMaxCellNoTRD->SetYTitle("Fraction");
-      fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
-      outputContainer->Add(fhEFracMaxCellNoTRD) ;
+      fhPtLambda0NoTRD  = new TH2F
+      ("hPtLambda0NoTRD","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+      fhPtLambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
+      fhPtLambda0NoTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhPtLambda0NoTRD) ;
+      
+      fhPtFracMaxCellNoTRD  = new TH2F
+      ("hPtFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
+      fhPtFracMaxCellNoTRD->SetYTitle("Fraction");
+      fhPtFracMaxCellNoTRD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhPtFracMaxCellNoTRD) ;
       
       if(!fFillOnlySimpleSSHisto)
       {
-        fhDispEtaE  = new TH2F ("hDispEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
-        fhDispEtaE->SetXTitle("E (GeV)");
-        fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
-        outputContainer->Add(fhDispEtaE);
-        
-        fhDispPhiE  = new TH2F ("hDispPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
-        fhDispPhiE->SetXTitle("E (GeV)");
-        fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
-        outputContainer->Add(fhDispPhiE);
-        
-        fhSumEtaE  = new TH2F ("hSumEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
-        fhSumEtaE->SetXTitle("E (GeV)");
-        fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
-        outputContainer->Add(fhSumEtaE);
-        
-        fhSumPhiE  = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
+        fhPtDispEta  = new TH2F ("hPtDispEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs #it{p}_{T}",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
+        fhPtDispEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        fhPtDispEta->SetYTitle("#sigma^{2}_{#eta #eta}");
+        outputContainer->Add(fhPtDispEta);
+        
+        fhPtDispPhi  = new TH2F ("hPtDispPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs #it{p}_{T}",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
+        fhPtDispPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        fhPtDispPhi->SetYTitle("#sigma^{2}_{#phi #phi}");
+        outputContainer->Add(fhPtDispPhi);
+        
+        fhPtSumEta  = new TH2F ("hPtSumEta","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs #it{p}_{T}",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
+        fhPtSumEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        fhPtSumEta->SetYTitle("#delta^{2}_{#eta #eta}");
+        outputContainer->Add(fhPtSumEta);
+        
+        fhPtSumPhi  = new TH2F ("hPtSumPhi","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs #it{p}_{T}",
                                nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
-        fhSumPhiE->SetXTitle("E (GeV)");
-        fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
-        outputContainer->Add(fhSumPhiE);
+        fhPtSumPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        fhPtSumPhi->SetYTitle("#delta^{2}_{#phi #phi}");
+        outputContainer->Add(fhPtSumPhi);
         
-        fhSumEtaPhiE  = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
+        fhPtSumEtaPhi  = new TH2F ("hPtSumEtaPhi","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs #it{p}_{T}",
                                   nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
-        fhSumEtaPhiE->SetXTitle("E (GeV)");
-        fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
-        outputContainer->Add(fhSumEtaPhiE);
+        fhPtSumEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        fhPtSumEtaPhi->SetYTitle("#delta^{2}_{#eta #phi}");
+        outputContainer->Add(fhPtSumEtaPhi);
         
-        fhDispEtaPhiDiffE  = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
+        fhPtDispEtaPhiDiff  = new TH2F ("hPtDispEtaPhiDiff","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs #it{p}_{T}",
                                        nptbins,ptmin,ptmax,200, -10,10);
-        fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
-        fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
-        outputContainer->Add(fhDispEtaPhiDiffE);
+        fhPtDispEtaPhiDiff->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        fhPtDispEtaPhiDiff->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
+        outputContainer->Add(fhPtDispEtaPhiDiff);
         
-        fhSphericityE  = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
+        fhPtSphericity  = new TH2F ("hPtSphericity","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs #it{p}_{T} (GeV/#it{c})",
                                    nptbins,ptmin,ptmax, 200, -1,1);
-        fhSphericityE->SetXTitle("E (GeV)");
-        fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
-        outputContainer->Add(fhSphericityE);
+        fhPtSphericity->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        fhPtSphericity->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
+        outputContainer->Add(fhPtSphericity);
         
         for(Int_t i = 0; i < 7; i++)
         {
@@ -1113,100 +1231,124 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
         }
       }
     }
-    
-    fhNLocMaxE = new TH2F("hNLocMaxE","Number of local maxima in cluster",
-                          nptbins,ptmin,ptmax,10,0,10);
-    fhNLocMaxE ->SetYTitle("N maxima");
-    fhNLocMaxE ->SetXTitle("E (GeV)");
-    outputContainer->Add(fhNLocMaxE) ;
+
+    fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster, selected clusters",
+                           nptbins,ptmin,ptmax,20,0,20);
+    fhNLocMaxPt ->SetYTitle("N maxima");
+    fhNLocMaxPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    outputContainer->Add(fhNLocMaxPt) ;
+
+    for(Int_t iSM = 0; iSM < GetCaloUtils()->GetNumberOfSuperModulesUsed(); iSM++)
+    {
+      fhNLocMaxPtSM[iSM] = new TH2F(Form("hNLocMaxPt_SM%d",iSM),Form("Number of local maxima in cluster, selected clusters in SM %d",iSM),
+                               nptbins,ptmin,ptmax,20,0,20);
+      fhNLocMaxPtSM[iSM] ->SetYTitle("N maxima");
+      fhNLocMaxPtSM[iSM] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhNLocMaxPtSM[iSM]) ;
+    }
     
     if(fAnaType == kSSCalo)
     {
-      fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster",
-                             nptbins,ptmin,ptmax,10,0,10);
-      fhNLocMaxPt ->SetYTitle("N maxima");
-      fhNLocMaxPt ->SetXTitle("p_{T} (GeV/c)");
-      outputContainer->Add(fhNLocMaxPt) ;
+
+      fhNLocMaxPtReject = new TH2F("hNLocMaxPtReject","Number of local maxima in cluster, rejected clusters",
+                             nptbins,ptmin,ptmax,20,0,20);
+      fhNLocMaxPtReject ->SetYTitle("N maxima");
+      fhNLocMaxPtReject ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhNLocMaxPtReject) ;
     }
     
     for (Int_t i = 0; i < 3; i++)
     {
-      fhELambda0LocMax[i]  = new TH2F(Form("hELambda0LocMax%d",i+1),
-                                      Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, %s",nlm[i].Data()),
+      fhPtLambda0LocMax[i]  = new TH2F(Form("hPtLambda0LocMax%d",i+1),
+                                      Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, NLM=%s",nlm[i].Data()),
                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-      fhELambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
-      fhELambda0LocMax[i]->SetXTitle("E (GeV)");
-      outputContainer->Add(fhELambda0LocMax[i]) ;
+      fhPtLambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
+      fhPtLambda0LocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhPtLambda0LocMax[i]) ;
+
+      if(IsDataMC())
+      {
+        for(Int_t ipart = 0; ipart < 6; ipart++)
+        {
+          fhMCPtLambda0LocMax[ipart][i]  = new TH2F
+          (Form("hPtLambda0LocMax%d_MC%s",i+1,pname[ipart].Data()),
+           Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{0}, NLM=%s, MC %s",nlm[i].Data(),pname[ipart].Data()),
+           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+          fhMCPtLambda0LocMax[ipart][i]->SetYTitle("#lambda_{0}^{2}");
+          fhMCPtLambda0LocMax[ipart][i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+          outputContainer->Add(fhMCPtLambda0LocMax[ipart][i]) ;
+        }
+      }
       
-      fhELambda1LocMax[i]  = new TH2F(Form("hELambda1LocMax%d",i+1),
-                                      Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, %s",nlm[i].Data()),
+      fhPtLambda1LocMax[i]  = new TH2F(Form("hPtLambda1LocMax%d",i+1),
+                                      Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #lambda_{1}, %s",nlm[i].Data()),
                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-      fhELambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
-      fhELambda1LocMax[i]->SetXTitle("E (GeV)");
-      outputContainer->Add(fhELambda1LocMax[i]) ;
+      fhPtLambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
+      fhPtLambda1LocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhPtLambda1LocMax[i]) ;
       
-      fhEDispersionLocMax[i]  = new TH2F(Form("hEDispersionLocMax%d",i+1),
-                                         Form("Selected #pi^{0} (#eta) pairs: E vs dispersion^{2}, %s",nlm[i].Data()),
+      fhPtDispersionLocMax[i]  = new TH2F(Form("hPtDispersionLocMax%d",i+1),
+                                         Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs dispersion^{2}, %s",nlm[i].Data()),
                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-      fhEDispersionLocMax[i]->SetYTitle("dispersion^{2}");
-      fhEDispersionLocMax[i]->SetXTitle("E (GeV)");
-      outputContainer->Add(fhEDispersionLocMax[i]) ;
+      fhPtDispersionLocMax[i]->SetYTitle("dispersion^{2}");
+      fhPtDispersionLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhPtDispersionLocMax[i]) ;
       
       if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
       {
-        fhEDispEtaLocMax[i]  = new TH2F(Form("hEDispEtaLocMax%d",i+1),
-                                        Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
+        fhPtDispEtaLocMax[i]  = new TH2F(Form("hPtDispEtaLocMax%d",i+1),
+                                        Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-        fhEDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
-        fhEDispEtaLocMax[i]->SetXTitle("E (GeV)");
-        outputContainer->Add(fhEDispEtaLocMax[i]) ;
+        fhPtDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
+        fhPtDispEtaLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhPtDispEtaLocMax[i]) ;
         
-        fhEDispPhiLocMax[i]  = new TH2F(Form("hEDispPhiLocMax%d",i+1),
-                                        Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
+        fhPtDispPhiLocMax[i]  = new TH2F(Form("hPtDispPhiLocMax%d",i+1),
+                                        Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-        fhEDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
-        fhEDispPhiLocMax[i]->SetXTitle("E (GeV)");
-        outputContainer->Add(fhEDispPhiLocMax[i]) ;
+        fhPtDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
+        fhPtDispPhiLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhPtDispPhiLocMax[i]) ;
         
-        fhESumEtaPhiLocMax[i]  = new TH2F(Form("hESumEtaPhiLocMax%d",i+1),
-                                          Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
+        fhPtSumEtaPhiLocMax[i]  = new TH2F(Form("hPtSumEtaPhiLocMax%d",i+1),
+                                          Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
                                           nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
-        fhESumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
-        fhESumEtaPhiLocMax[i]->SetXTitle("E (GeV)");
-        outputContainer->Add(fhESumEtaPhiLocMax[i]) ;
+        fhPtSumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
+        fhPtSumEtaPhiLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhPtSumEtaPhiLocMax[i]) ;
         
-        fhEDispEtaPhiDiffLocMax[i]  = new TH2F(Form("hEDispEtaPhiDiffLocMax%d",i+1),
-                                               Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
+        fhPtDispEtaPhiDiffLocMax[i]  = new TH2F(Form("hPtDispEtaPhiDiffLocMax%d",i+1),
+                                               Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
                                                nptbins,ptmin,ptmax,200, -10,10);
-        fhEDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
-        fhEDispEtaPhiDiffLocMax[i]->SetXTitle("E (GeV)");
-        outputContainer->Add(fhEDispEtaPhiDiffLocMax[i]) ;
+        fhPtDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
+        fhPtDispEtaPhiDiffLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhPtDispEtaPhiDiffLocMax[i]) ;
         
-        fhESphericityLocMax[i]  = new TH2F(Form("hESphericityLocMax%d",i+1),
-                                           Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
+        fhPtSphericityLocMax[i]  = new TH2F(Form("hPtSphericityLocMax%d",i+1),
+                                           Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
                                            nptbins,ptmin,ptmax,200, -1,1);
-        fhESphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
-        fhESphericityLocMax[i]->SetXTitle("E (GeV)");
-        outputContainer->Add(fhESphericityLocMax[i]) ;
+        fhPtSphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
+        fhPtSphericityLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhPtSphericityLocMax[i]) ;
       }
       
     }
     
-    fhENCells  = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
-    fhENCells->SetXTitle("E (GeV)");
-    fhENCells->SetYTitle("# of cells in cluster");
-    outputContainer->Add(fhENCells);
+    fhPtNCells  = new TH2F ("hPtNCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
+    fhPtNCells->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhPtNCells->SetYTitle("# of cells in cluster");
+    outputContainer->Add(fhPtNCells);
     
-    fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
-    fhETime->SetXTitle("E (GeV)");
-    fhETime->SetYTitle("t (ns)");
-    outputContainer->Add(fhETime);
+    fhPtTime = new TH2F("hPtTime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
+    fhPtTime->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhPtTime->SetYTitle("t (ns)");
+    outputContainer->Add(fhPtTime);
     
   }
   
   
   fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
-  fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
+  fhEPairDiffTime->SetXTitle("#it{E}_{pair} (GeV)");
   fhEPairDiffTime->SetYTitle("#Delta t (ns)");
   outputContainer->Add(fhEPairDiffTime);
   
@@ -1225,10 +1367,10 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
       
       fhMassPairLocMax[i]  = new TH2F
       (Form("MassPairLocMax%s",combiName[i].Data()),
-       Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
+       Form("#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
        nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
-      fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
-      fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
+      fhMassPairLocMax[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+      fhMassPairLocMax[i]->SetXTitle("#it{E}_{pair} (GeV)");
       outputContainer->Add(fhMassPairLocMax[i]) ;
     }
   }
@@ -1237,21 +1379,21 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   {
     fhTrackMatchedDEta  = new TH2F
     ("hTrackMatchedDEta",
-     "d#eta of cluster-track vs cluster energy",
+     "d#eta of cluster-track vs cluster #it{p}_{T}",
      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
     fhTrackMatchedDEta->SetYTitle("d#eta");
-    fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
+    fhTrackMatchedDEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     
     fhTrackMatchedDPhi  = new TH2F
     ("hTrackMatchedDPhi",
-     "d#phi of cluster-track vs cluster energy",
+     "d#phi of cluster-track vs cluster #it{p}_{T}",
      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
-    fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
+    fhTrackMatchedDPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     
     fhTrackMatchedDEtaDPhi  = new TH2F
     ("hTrackMatchedDEtaDPhi",
-     "d#eta vs d#phi of cluster-track vs cluster energy",
+     "d#eta vs d#phi of cluster-track",
      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
@@ -1262,21 +1404,21 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
 
     fhTrackMatchedDEtaPos  = new TH2F
     ("hTrackMatchedDEtaPos",
-     "d#eta of cluster-track vs cluster energy",
+     "d#eta of cluster-track vs cluster #it{p}_{T}",
      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
     fhTrackMatchedDEtaPos->SetYTitle("d#eta");
-    fhTrackMatchedDEtaPos->SetXTitle("E_{cluster} (GeV)");
+    fhTrackMatchedDEtaPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     
     fhTrackMatchedDPhiPos  = new TH2F
     ("hTrackMatchedDPhiPos",
-     "d#phi of cluster-track vs cluster energy",
+     "d#phi of cluster-track vs cluster #it{p}_{T}",
      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
     fhTrackMatchedDPhiPos->SetYTitle("d#phi (rad)");
-    fhTrackMatchedDPhiPos->SetXTitle("E_{cluster} (GeV)");
+    fhTrackMatchedDPhiPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     
     fhTrackMatchedDEtaDPhiPos  = new TH2F
     ("hTrackMatchedDEtaDPhiPos",
-     "d#eta vs d#phi of cluster-track vs cluster energy",
+     "d#eta vs d#phi of cluster-track",
      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
     fhTrackMatchedDEtaDPhiPos->SetYTitle("d#phi (rad)");
     fhTrackMatchedDEtaDPhiPos->SetXTitle("d#eta");
@@ -1287,21 +1429,21 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
 
     fhTrackMatchedDEtaNeg  = new TH2F
     ("hTrackMatchedDEtaNeg",
-     "d#eta of cluster-track vs cluster energy",
+     "d#eta of cluster-track vs cluster #it{p}_{T}",
      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
     fhTrackMatchedDEtaNeg->SetYTitle("d#eta");
-    fhTrackMatchedDEtaNeg->SetXTitle("E_{cluster} (GeV)");
+    fhTrackMatchedDEtaNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     
     fhTrackMatchedDPhiNeg  = new TH2F
     ("hTrackMatchedDPhiNeg",
-     "d#phi of cluster-track vs cluster energy",
+     "d#phi of cluster-track vs cluster #it{p}_{T}",
      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
     fhTrackMatchedDPhiNeg->SetYTitle("d#phi (rad)");
-    fhTrackMatchedDPhiNeg->SetXTitle("E_{cluster} (GeV)");
+    fhTrackMatchedDPhiNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     
     fhTrackMatchedDEtaDPhiNeg  = new TH2F
     ("hTrackMatchedDEtaDPhiNeg",
-     "d#eta vs d#phi of cluster-track vs cluster energy",
+     "d#eta vs d#phi of cluster-track",
      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
     fhTrackMatchedDEtaDPhiNeg->SetYTitle("d#phi (rad)");
     fhTrackMatchedDEtaDPhiNeg->SetXTitle("d#eta");
@@ -1310,43 +1452,43 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
     outputContainer->Add(fhTrackMatchedDPhiNeg) ;
     outputContainer->Add(fhTrackMatchedDEtaDPhiNeg) ;
     
-    fhdEdx  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
-    fhdEdx->SetXTitle("E (GeV)");
-    fhdEdx->SetYTitle("<dE/dx>");
+    fhdEdx  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster #it{p}_{T}", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
+    fhdEdx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhdEdx->SetYTitle("<#it{dE}/#it{dx}>");
     outputContainer->Add(fhdEdx);
     
-    fhEOverP  = new TH2F ("hEOverP","matched track E/p vs cluster ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
-    fhEOverP->SetXTitle("E (GeV)");
-    fhEOverP->SetYTitle("E/p");
+    fhEOverP  = new TH2F ("hEOverP","matched track E/p vs cluster #it{p}_{T}", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
+    fhEOverP->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhEOverP->SetYTitle("#it{E}/#it{p}");
     outputContainer->Add(fhEOverP);
     
     if(fCalorimeter=="EMCAL")
     {
       fhEOverPNoTRD  = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
-      fhEOverPNoTRD->SetXTitle("E (GeV)");
-      fhEOverPNoTRD->SetYTitle("E/p");
+      fhEOverPNoTRD->SetXTitle("#it{E} (GeV)");
+      fhEOverPNoTRD->SetYTitle("#it{E}/#it{p}");
       outputContainer->Add(fhEOverPNoTRD);
     }
     
     if(IsDataMC() && fFillTMHisto)
     {
-      fhTrackMatchedMCParticleE  = new TH2F
-      ("hTrackMatchedMCParticleE",
+      fhTrackMatchedMCParticlePt  = new TH2F
+      ("hTrackMatchedMCParticlePt",
        "Origin of particle vs energy",
        nptbins,ptmin,ptmax,8,0,8);
-      fhTrackMatchedMCParticleE->SetXTitle("E (GeV)");
-      //fhTrackMatchedMCParticleE->SetYTitle("Particle type");
+      fhTrackMatchedMCParticlePt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      //fhTrackMatchedMCParticlePt->SetYTitle("Particle type");
       
-      fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(1 ,"Photon");
-      fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(2 ,"Electron");
-      fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
-      fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(4 ,"Rest");
-      fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
-      fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
-      fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
-      fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
+      fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(1 ,"Photon");
+      fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(2 ,"Electron");
+      fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
+      fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(4 ,"Rest");
+      fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
+      fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
+      fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
+      fhTrackMatchedMCParticlePt->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
       
-      outputContainer->Add(fhTrackMatchedMCParticleE);
+      outputContainer->Add(fhTrackMatchedMCParticlePt);
       
       fhTrackMatchedMCParticleDEta  = new TH2F
       ("hTrackMatchedMCParticleDEta",
@@ -1392,39 +1534,39 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   {
     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
                                      nptbins,ptmin,ptmax, 100,0,1.);
-    fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
-    fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
+    fhECellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
+    fhECellClusterRatio->SetYTitle("#it{E}_{cell i}/#it{E}_{cluster}");
     outputContainer->Add(fhECellClusterRatio);
     
     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
                                         nptbins,ptmin,ptmax, 100,-10,0);
-    fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
-    fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
+    fhECellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
+    fhECellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
     outputContainer->Add(fhECellClusterLogRatio);
     
     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
                                         nptbins,ptmin,ptmax, 100,0,1.);
-    fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
-    fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
+    fhEMaxCellClusterRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
+    fhEMaxCellClusterRatio->SetYTitle("#it{E}_{max cell}/#it{E}_{cluster}");
     outputContainer->Add(fhEMaxCellClusterRatio);
     
     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
                                            nptbins,ptmin,ptmax, 100,-10,0);
-    fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
-    fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
+    fhEMaxCellClusterLogRatio->SetXTitle("#it{E}_{cluster} (GeV) ");
+    fhEMaxCellClusterLogRatio->SetYTitle("Log (#it{E}_{max cell}/#it{E}_{cluster})");
     outputContainer->Add(fhEMaxCellClusterLogRatio);
     
     for(Int_t iw = 0; iw < 14; iw++)
     {
       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected decay photons from neutral meson",1+0.5*iw),
                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-      fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
+      fhLambda0ForW0[iw]->SetXTitle("#it{E}_{cluster}");
       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
       outputContainer->Add(fhLambda0ForW0[iw]);
       
       //      fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected decay photons from neutral meson",0.5+0.5*iw),
       //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-      //      fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
+      //      fhLambda1ForW0[iw]->SetXTitle("#it{E}_{cluster}");
       //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
       //      outputContainer->Add(fhLambda1ForW0[iw]);
       
@@ -1433,45 +1575,93 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   
   if(IsDataMC())
   {
+    // Origin
+    
+    fhMCPi0PtOrigin     = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
+    fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhMCPi0PtOrigin->SetYTitle("Origin");
+    fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
+    fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
+    fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
+    fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
+    fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
+    fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
+    fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
+    fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
+    fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
+    fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
+    outputContainer->Add(fhMCPi0PtOrigin) ;
+    
+    fhMCEtaPtOrigin     = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
+    fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhMCEtaPtOrigin->SetYTitle("Origin");
+    fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
+    fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
+    fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
+    fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
+    fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
+    fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
+    outputContainer->Add(fhMCEtaPtOrigin) ;
+    
+    fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",nptbins,ptmin,ptmax,2000,0,500) ;
+    fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
+    outputContainer->Add(fhMCPi0ProdVertex) ;
+    
+    fhMCPi0ProdVertexInner = new TH2F("hMCPi0ProdVertexInner","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",nptbins,ptmin,ptmax,500,0,50) ;
+    fhMCPi0ProdVertexInner->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhMCPi0ProdVertexInner->SetYTitle("#it{R} (cm)");
+    outputContainer->Add(fhMCPi0ProdVertexInner) ;
+    
+    fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",nptbins,ptmin,ptmax,2000,0,500) ;
+    fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
+    outputContainer->Add(fhMCEtaProdVertex) ;
+    
+    fhMCEtaProdVertexInner = new TH2F("hMCEtaProdVertexInner","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",nptbins,ptmin,ptmax,500,0,50) ;
+    fhMCEtaProdVertexInner->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhMCEtaProdVertexInner->SetYTitle("#it{R} (cm)");
+    outputContainer->Add(fhMCEtaProdVertexInner) ;
+    
     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
     {
-      fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), pT versus E primary #pi^{0} / E reco",
+      fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), #it{p}_{T} versus E primary #pi^{0} / E reco",
                                           nptbins,ptmin,ptmax,200,0,2);
-      fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
-      fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
+      fhMCPi0PtGenRecoFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
+      fhMCPi0PtGenRecoFraction->SetYTitle("#it{E}^{#pi^{0} mother} / #it{E}^{rec}");
       outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
       
-      fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
+      fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),#it{p}_{T} versus E primary #eta / E reco",
                                           nptbins,ptmin,ptmax,200,0,2);
-      fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
-      fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
+      fhMCEtaPtGenRecoFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
+      fhMCEtaPtGenRecoFraction->SetYTitle("#it{E}^{ #eta mother} / #it{E}^{rec}");
       outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
       
       fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
-      fhMCPi0DecayPt->SetYTitle("N");
-      fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
+      fhMCPi0DecayPt->SetYTitle("#it{N}");
+      fhMCPi0DecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
       outputContainer->Add(fhMCPi0DecayPt) ;
       
-      fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta), pT versus E primary  #gamma / E primary #pi^{0}",
+      fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta), #it{p}_{T} versus E primary  #gamma / #it{E} primary #pi^{0}",
                                         nptbins,ptmin,ptmax,100,0,1);
-      fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
+      fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/#it{c})");
       fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
       outputContainer->Add(fhMCPi0DecayPtFraction) ;
       
       fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
-      fhMCEtaDecayPt->SetYTitle("N");
-      fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
+      fhMCEtaDecayPt->SetYTitle("#it{N}");
+      fhMCEtaDecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
       outputContainer->Add(fhMCEtaDecayPt) ;
       
-      fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay  identified as #pi^{0} (#eta), pT versus E primary  #gamma / E primary #eta",
+      fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay  identified as #pi^{0} (#eta), #it{p}_{T} versus E primary  #gamma / E primary #eta",
                                         nptbins,ptmin,ptmax,100,0,1);
-      fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
-      fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
+      fhMCEtaDecayPtFraction->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
+      fhMCEtaDecayPtFraction->SetYTitle("#it{E}^{gen} / #it{E}^{gen-mother}");
       outputContainer->Add(fhMCEtaDecayPtFraction) ;
       
       fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0})  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
-      fhMCOtherDecayPt->SetYTitle("N");
-      fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
+      fhMCOtherDecayPt->SetYTitle("#it{N}");
+      fhMCOtherDecayPt->SetXTitle("#it{p}^{rec}_{T} (GeV/#it{c})");
       outputContainer->Add(fhMCOtherDecayPt) ;
       
     }
@@ -1482,32 +1672,32 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
       
       fhAnglePairMCPi0  = new TH2F
       ("AnglePairMCPi0",
-       "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
+       "Angle between decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
       fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
-      fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
+      fhAnglePairMCPi0->SetXTitle("#it{E}_{pair} (GeV)");
       outputContainer->Add(fhAnglePairMCPi0) ;
       
       if (fAnaType!= kSSCalo)
       {
         fhAnglePairMCEta  = new TH2F
         ("AnglePairMCEta",
-         "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
+         "Angle between decay #gamma pair vs #it{E}_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
         fhAnglePairMCEta->SetYTitle("#alpha (rad)");
-        fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
+        fhAnglePairMCEta->SetXTitle("#it{E}_{pair} (GeV)");
         outputContainer->Add(fhAnglePairMCEta) ;
         
         fhMassPairMCPi0  = new TH2F
         ("MassPairMCPi0",
-         "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
-        fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
-        fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
+         "#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
+        fhMassPairMCPi0->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+        fhMassPairMCPi0->SetXTitle("#it{E}_{pair} (GeV)");
         outputContainer->Add(fhMassPairMCPi0) ;
         
         fhMassPairMCEta  = new TH2F
         ("MassPairMCEta",
-         "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
-        fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
-        fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
+         "#it{M} for decay #gamma pair vs #it{E}_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
+        fhMassPairMCEta->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+        fhMassPairMCEta->SetXTitle("#it{E}_{pair} (GeV)");
         outputContainer->Add(fhMassPairMCEta) ;
       }
       
@@ -1519,8 +1709,8 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
          Form("Identified as #pi^{0} (#eta), cluster from %s",
               ptype[i].Data()),
          nptbins,ptmin,ptmax);
-        fhMCE[i]->SetYTitle("N");
-        fhMCE[i]->SetXTitle("E (GeV)");
+        fhMCE[i]->SetYTitle("#it{N}");
+        fhMCE[i]->SetXTitle("#it{E} (GeV)");
         outputContainer->Add(fhMCE[i]) ;
         
         fhMCPt[i]  = new TH1F
@@ -1528,8 +1718,8 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
          Form("Identified as #pi^{0} (#eta), cluster from %s",
               ptype[i].Data()),
          nptbins,ptmin,ptmax);
-        fhMCPt[i]->SetYTitle("N");
-        fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
+        fhMCPt[i]->SetYTitle("#it{N}");
+        fhMCPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhMCPt[i]) ;
         
         fhMCPtCentrality[i]  = new TH2F
@@ -1538,27 +1728,34 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
               ptype[i].Data()),
          nptbins,ptmin,ptmax, 100,0,100);
         fhMCPtCentrality[i]->SetYTitle("centrality");
-        fhMCPtCentrality[i]->SetXTitle("p_{T} (GeV/c)");
+        fhMCPtCentrality[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhMCPtCentrality[i]) ;
         
         if(fAnaType == kSSCalo)
         {
-          
           fhMCNLocMaxPt[i] = new TH2F
           (Form("hNLocMaxPt_MC%s",pname[i].Data()),
-           Form("cluster from %s, pT of cluster, for NLM",ptype[i].Data()),
-           nptbins,ptmin,ptmax,10,0,10);
-          fhMCNLocMaxPt[i] ->SetYTitle("N maxima");
-          fhMCNLocMaxPt[i] ->SetXTitle("p_{T} (GeV/c)");
+           Form("cluster from %s, #it{p}_{T} of cluster vs NLM, accepted",ptype[i].Data()),
+           nptbins,ptmin,ptmax,20,0,20);
+          fhMCNLocMaxPt[i] ->SetYTitle("#it{NLM}");
+          fhMCNLocMaxPt[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
           outputContainer->Add(fhMCNLocMaxPt[i]) ;
+          fhMCNLocMaxPtReject[i] = new TH2F
+          (Form("hNLocMaxPtReject_MC%s",pname[i].Data()),
+           Form("cluster from %s, #it{p}_{T} of cluster vs NLM, rejected",ptype[i].Data()),
+           nptbins,ptmin,ptmax,20,0,20);
+          fhMCNLocMaxPtReject[i] ->SetYTitle("#it{NLM}");
+          fhMCNLocMaxPtReject[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+          outputContainer->Add(fhMCNLocMaxPtReject[i]) ;
           
           fhMCEReject[i]  = new TH1F
           (Form("hEReject_MC%s",pname[i].Data()),
            Form("Rejected as #pi^{0} (#eta), cluster from %s",
                 ptype[i].Data()),
            nptbins,ptmin,ptmax);
-          fhMCEReject[i]->SetYTitle("N");
-          fhMCEReject[i]->SetXTitle("E (GeV)");
+          fhMCEReject[i]->SetYTitle("#it{N}");
+          fhMCEReject[i]->SetXTitle("#it{E} (GeV)");
           outputContainer->Add(fhMCEReject[i]) ;
           
           fhMCPtReject[i]  = new TH1F
@@ -1566,112 +1763,130 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
            Form("Rejected as #pi^{0} (#eta), cluster from %s",
                 ptype[i].Data()),
            nptbins,ptmin,ptmax);
-          fhMCPtReject[i]->SetYTitle("N");
-          fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
+          fhMCPtReject[i]->SetYTitle("#it{N}");
+          fhMCPtReject[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
           outputContainer->Add(fhMCPtReject[i]) ;
         }
         
-        fhMCPhi[i]  = new TH2F
-        (Form("hPhi_MC%s",pname[i].Data()),
+        fhMCPtPhi[i]  = new TH2F
+        (Form("hPtPhi_MC%s",pname[i].Data()),
          Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
-        fhMCPhi[i]->SetYTitle("#phi");
-        fhMCPhi[i]->SetXTitle("p_{T} (GeV/c)");
-        outputContainer->Add(fhMCPhi[i]) ;
+        fhMCPtPhi[i]->SetYTitle("#phi");
+        fhMCPtPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhMCPtPhi[i]) ;
         
-        fhMCEta[i]  = new TH2F
-        (Form("hEta_MC%s",pname[i].Data()),
+        fhMCPtEta[i]  = new TH2F
+        (Form("hPtEta_MC%s",pname[i].Data()),
          Form("Identified as #pi^{0} (#eta), cluster from %s",
               ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
-        fhMCEta[i]->SetYTitle("#eta");
-        fhMCEta[i]->SetXTitle("p_{T} (GeV/c)");
-        outputContainer->Add(fhMCEta[i]) ;
+        fhMCPtEta[i]->SetYTitle("#eta");
+        fhMCPtEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhMCPtEta[i]) ;
         
         fhMCMassPt[i]  = new TH2F
         (Form("hMassPt_MC%s",pname[i].Data()),
-         Form("all pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
+         Form("all pairs #it{M}: #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
-        fhMCMassPt[i]->SetYTitle("mass (GeV/c^{2})");
-        fhMCMassPt[i]->SetXTitle("p_{T} (GeV/c)");
+        fhMCMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+        fhMCMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhMCMassPt[i]) ;
         
         fhMCSelectedMassPt[i]  = new TH2F
         (Form("hSelectedMassPt_MC%s",pname[i].Data()),
-         Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
+         Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
-        fhMCSelectedMassPt[i]->SetYTitle("mass (GeV/c^{2})");
-        fhMCSelectedMassPt[i]->SetXTitle("p_{T} (GeV/c)");
+        fhMCSelectedMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+        fhMCSelectedMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhMCSelectedMassPt[i]) ;
         
+        if(fAnaType == kSSCalo)
+        {
+          fhMCMassPtNoOverlap[i]  = new TH2F
+          (Form("hMassPtNoOverlap_MC%s",pname[i].Data()),
+           Form("all pairs #it{M}: #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
+           nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+          fhMCMassPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+          fhMCMassPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+          outputContainer->Add(fhMCMassPtNoOverlap[i]) ;
+          
+          fhMCSelectedMassPtNoOverlap[i]  = new TH2F
+          (Form("hSelectedMassPtNoOverlap_MC%s",pname[i].Data()),
+           Form("Selected #pi^{0} (#eta) pairs #it{M}: #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
+           nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+          fhMCSelectedMassPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+          fhMCSelectedMassPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+          outputContainer->Add(fhMCSelectedMassPtNoOverlap[i]) ;
+        }
         
         if( fFillSelectClHisto )
         {
-          fhEMCLambda0[i]  = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
-                                      Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
+          fhMCPtLambda0[i]  = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
+                                      Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}",ptype[i].Data()),
                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-          fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
-          fhEMCLambda0[i]->SetXTitle("E (GeV)");
-          outputContainer->Add(fhEMCLambda0[i]) ;
+          fhMCPtLambda0[i]->SetYTitle("#lambda_{0}^{2}");
+          fhMCPtLambda0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+          outputContainer->Add(fhMCPtLambda0[i]) ;
           
-          fhEMCLambda1[i]  = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
-                                      Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
+          fhMCPtLambda1[i]  = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
+                                      Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{1}^{2}",ptype[i].Data()),
                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-          fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
-          fhEMCLambda1[i]->SetXTitle("E (GeV)");
-          outputContainer->Add(fhEMCLambda1[i]) ;
+          fhMCPtLambda1[i]->SetYTitle("#lambda_{1}^{2}");
+          fhMCPtLambda1[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+          outputContainer->Add(fhMCPtLambda1[i]) ;
           
-          fhEMCDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
-                                         Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
+          fhMCPtDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
+                                         Form("Selected pair, cluster from %s : #it{p}_{T} vs dispersion^{2}",ptype[i].Data()),
                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-          fhEMCDispersion[i]->SetYTitle("D^{2}");
-          fhEMCDispersion[i]->SetXTitle("E (GeV)");
-          outputContainer->Add(fhEMCDispersion[i]) ;
+          fhMCPtDispersion[i]->SetYTitle("#it{D}^{2}");
+          fhMCPtDispersion[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+          outputContainer->Add(fhMCPtDispersion[i]) ;
           
           if(fCalorimeter=="EMCAL")
           {
-            fhEMCLambda0NoTRD[i]  = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
-                                             Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
+            fhMCPtLambda0NoTRD[i]  = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
+                                             Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
                                              nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-            fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
-            fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
-            outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
+            fhMCPtLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
+            fhMCPtLambda0NoTRD[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+            outputContainer->Add(fhMCPtLambda0NoTRD[i]) ;
             
             if(!fFillOnlySimpleSSHisto)
             {
-              fhMCEDispEta[i]  = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
-                                           Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
+              fhMCPtDispEta[i]  = new TH2F (Form("hPtDispEta_MC%s",pname[i].Data()),
+                                           Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs #it{p}_{T}",ptype[i].Data()),
                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
-              fhMCEDispEta[i]->SetXTitle("E (GeV)");
-              fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
-              outputContainer->Add(fhMCEDispEta[i]);
+              fhMCPtDispEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+              fhMCPtDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
+              outputContainer->Add(fhMCPtDispEta[i]);
               
-              fhMCEDispPhi[i]  = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
-                                           Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
+              fhMCPtDispPhi[i]  = new TH2F (Form("hPtDispPhi_MC%s",pname[i].Data()),
+                                           Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs #it{p}_{T}",ptype[i].Data()),
                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
-              fhMCEDispPhi[i]->SetXTitle("E (GeV)");
-              fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
-              outputContainer->Add(fhMCEDispPhi[i]);
+              fhMCPtDispPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+              fhMCPtDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+              outputContainer->Add(fhMCPtDispPhi[i]);
               
-              fhMCESumEtaPhi[i]  = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
-                                             Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptype[i].Data()),
+              fhMCPtSumEtaPhi[i]  = new TH2F (Form("hPtSumEtaPhi_MC%s",pname[i].Data()),
+                                             Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs #it{p}_{T}",ptype[i].Data()),
                                              nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
-              fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
-              fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
-              outputContainer->Add(fhMCESumEtaPhi[i]);
+              fhMCPtSumEtaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+              fhMCPtSumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
+              outputContainer->Add(fhMCPtSumEtaPhi[i]);
               
-              fhMCEDispEtaPhiDiff[i]  = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
-                                                  Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
+              fhMCPtDispEtaPhiDiff[i]  = new TH2F (Form("hPtDispEtaPhiDiff_MC%s",pname[i].Data()),
+                                                  Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs #it{p}_{T}",ptype[i].Data()),
                                                   nptbins,ptmin,ptmax,200,-10,10);
-              fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
-              fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
-              outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
+              fhMCPtDispEtaPhiDiff[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+              fhMCPtDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
+              outputContainer->Add(fhMCPtDispEtaPhiDiff[i]);
               
-              fhMCESphericity[i]  = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
+              fhMCPtSphericity[i]  = new TH2F (Form("hPtSphericity_MC%s",pname[i].Data()),
                                               Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptype[i].Data()),
                                               nptbins,ptmin,ptmax, 200,-1,1);
-              fhMCESphericity[i]->SetXTitle("E (GeV)");
-              fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
-              outputContainer->Add(fhMCESphericity[i]);
+              fhMCPtSphericity[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+              fhMCPtSphericity[i]->SetYTitle("#it{s} = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
+              outputContainer->Add(fhMCPtSphericity[i]);
               
               for(Int_t ie = 0; ie < 7; ie++)
               {
@@ -1700,19 +1915,19 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
             }
           }
           
-          fhEMCLambda0FracMaxCellCut[i]  = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
-                                                    Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
+          fhMCPtLambda0FracMaxCellCut[i]  = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
+                                                    Form("Selected pair, cluster from %s : #it{p}_{T} vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
-          fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
-          fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
-          outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
+          fhMCPtLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
+          fhMCPtLambda0FracMaxCellCut[i]->SetXTitle("#it{E} (GeV)");
+          outputContainer->Add(fhMCPtLambda0FracMaxCellCut[i]) ;
           
-          fhEMCFracMaxCell[i]  = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
-                                          Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
+          fhMCPtFracMaxCell[i]  = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
+                                          Form("Selected pair, cluster from %s : #it{p}_{T} vs Max cell fraction of energy",ptype[i].Data()),
                                           nptbins,ptmin,ptmax,100,0,1);
-          fhEMCFracMaxCell[i]->SetYTitle("Fraction");
-          fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
-          outputContainer->Add(fhEMCFracMaxCell[i]) ;
+          fhMCPtFracMaxCell[i]->SetYTitle("#it{Fraction}");
+          fhMCPtFracMaxCell[i]->SetXTitle("#it{E} (GeV)");
+          outputContainer->Add(fhMCPtFracMaxCell[i]) ;
           
         }//
       } // shower shape histo
@@ -1722,90 +1937,280 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   
   if(fAnaType==kSSCalo)
   {
-    fhAsymmetry  = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
+    fhAsymmetry  = new TH2F ("hAsymmetry","#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",
                              nptbins,ptmin,ptmax, 200, -1,1);
-    fhAsymmetry->SetXTitle("E (GeV)");
-    fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+    fhAsymmetry->SetXTitle("#it{E} (GeV)");
+    fhAsymmetry->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
     outputContainer->Add(fhAsymmetry);
     
-    fhSelectedAsymmetry  = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
+    fhSelectedAsymmetry  = new TH2F ("hSelectedAsymmetry","#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",
                                      nptbins,ptmin,ptmax, 200, -1,1);
-    fhSelectedAsymmetry->SetXTitle("E (GeV)");
-    fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+    fhSelectedAsymmetry->SetXTitle("#it{E} (GeV)");
+    fhSelectedAsymmetry->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
     outputContainer->Add(fhSelectedAsymmetry);
     
     fhSplitE  = new TH1F
     ("hSplitE","Selected #pi^{0} (#eta) pairs energy sum of split sub-clusters",nptbins,ptmin,ptmax);
     fhSplitE->SetYTitle("counts");
-    fhSplitE->SetXTitle("E (GeV)");
+    fhSplitE->SetXTitle("#it{E} (GeV)");
     outputContainer->Add(fhSplitE) ;
     
     fhSplitPt  = new TH1F
-    ("hSplitPt","Selected #pi^{0} (#eta) pairs pT sum of split sub-clusters",nptbins,ptmin,ptmax);
+    ("hSplitPt","Selected #pi^{0} (#eta) pairs #it{p}_{T} sum of split sub-clusters",nptbins,ptmin,ptmax);
     fhSplitPt->SetYTitle("counts");
-    fhSplitPt->SetXTitle("p_{T} (GeV/c)");
+    fhSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhSplitPt) ;
     
     
     fhSplitPtPhi  = new TH2F
-    ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
+    ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster #it{p}_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
     fhSplitPtPhi->SetYTitle("#phi (rad)");
-    fhSplitPtPhi->SetXTitle("p_{T} (GeV/c)");
+    fhSplitPtPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhSplitPtPhi) ;
     
     fhSplitPtEta  = new TH2F
-    ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
+    ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster #it{p}_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
     fhSplitPtEta->SetYTitle("#eta");
-    fhSplitPtEta->SetXTitle("p_{T} (GeV/c)");
+    fhSplitPtEta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhSplitPtEta) ;
     
     
     fhNLocMaxSplitPt = new TH2F("hNLocMaxSplitPt","Number of local maxima in cluster",
-                                nptbins,ptmin,ptmax,10,0,10);
-    fhNLocMaxSplitPt ->SetYTitle("N maxima");
-    fhNLocMaxSplitPt ->SetXTitle("p_{T} (GeV/c)");
+                                nptbins,ptmin,ptmax,20,0,20);
+    fhNLocMaxSplitPt ->SetYTitle("#it{NLM}");
+    fhNLocMaxSplitPt ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhNLocMaxSplitPt) ;
     
     
     fhMassSplitPt  = new TH2F
-    ("hMassSplitPt","all pairs mass: sum split sub-cluster p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
-    fhMassSplitPt->SetYTitle("mass (GeV/c^{2})");
-    fhMassSplitPt->SetXTitle("p_{T} (GeV/c)");
+    ("hMassSplitPt","all pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}",
+     nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+    fhMassSplitPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+    fhMassSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhMassSplitPt) ;
     
     fhSelectedMassSplitPt  = new TH2F
-    ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs mass: sum split sub-cluster p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
-    fhSelectedMassSplitPt->SetYTitle("mass (GeV/c^{2})");
-    fhSelectedMassSplitPt->SetXTitle("p_{T} (GeV/c)");
+    ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}",
+     nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+    fhSelectedMassSplitPt->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+    fhSelectedMassSplitPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhSelectedMassSplitPt) ;
     
-    
-    
     if(IsDataMC())
     {
+      fhMassSplitPtNoOverlap  = new TH2F
+      ("hMassSplitPtNoOverlap","all pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}, no overlap",
+       nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+      fhMassSplitPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+      fhMassSplitPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhMassSplitPtNoOverlap) ;
+      
+      fhSelectedMassSplitPtNoOverlap  = new TH2F
+      ("hSelectedMassSplitPtNoOverlap","Selected #pi^{0} (#eta) pairs #it{M}: sum split sub-cluster #it{p}_{T} vs #it{M}, no overlap",
+       nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+      fhSelectedMassSplitPtNoOverlap->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+      fhSelectedMassSplitPtNoOverlap->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhSelectedMassSplitPtNoOverlap) ;
+
+      
+      fhMCPi0PtRecoPtPrim  = new TH2F
+      ("hMCPi0PtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCPi0PtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCPi0PtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCPi0PtRecoPtPrim ) ;
+      
+      fhMCPi0PtRecoPtPrimNoOverlap  = new TH2F
+      ("hMCPi0PtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCPi0PtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCPi0PtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCPi0PtRecoPtPrimNoOverlap ) ;
+      
+      fhMCPi0SelectedPtRecoPtPrim  = new TH2F
+      ("hMCPi0SelectedPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCPi0SelectedPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCPi0SelectedPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCPi0SelectedPtRecoPtPrim ) ;
+      
+      fhMCPi0SelectedPtRecoPtPrimNoOverlap  = new TH2F
+      ("hMCPi0SelectedPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCPi0SelectedPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimNoOverlap ) ;
+
+      
+      fhMCPi0SplitPtRecoPtPrim  = new TH2F
+      ("hMCPi0SplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCPi0SplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCPi0SplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCPi0SplitPtRecoPtPrim ) ;
+      
+      fhMCPi0SplitPtRecoPtPrimNoOverlap  = new TH2F
+      ("hMCPi0SplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCPi0SplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCPi0SplitPtRecoPtPrimNoOverlap ) ;
+      
+      fhMCPi0SelectedSplitPtRecoPtPrim  = new TH2F
+      ("hMCPi0SelectedSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCPi0SelectedSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCPi0SelectedSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrim ) ;
+      
+      fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap  = new TH2F
+      ("hMCPi0SelectedSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap ) ;
+
+      fhMCEtaPtRecoPtPrim  = new TH2F
+      ("hMCEtaPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCEtaPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCEtaPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCEtaPtRecoPtPrim ) ;
+      
+      fhMCEtaPtRecoPtPrimNoOverlap  = new TH2F
+      ("hMCEtaPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCEtaPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCEtaPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCEtaPtRecoPtPrimNoOverlap ) ;
+      
+      fhMCEtaSelectedPtRecoPtPrim  = new TH2F
+      ("hMCEtaSelectedPtRecoPtPrim","#it{p}_{T,reco} vs #it{p}_{T,gen}",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCEtaSelectedPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCEtaSelectedPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCEtaSelectedPtRecoPtPrim ) ;
+      
+      fhMCEtaSelectedPtRecoPtPrimNoOverlap  = new TH2F
+      ("hMCEtaSelectedPtRecoPtPrimNoOverlap","#it{p}_{T,reco} vs #it{p}_{T,gen}, no overlap",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCEtaSelectedPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimNoOverlap ) ;
+      
+      
+      fhMCEtaSplitPtRecoPtPrim  = new TH2F
+      ("hMCEtaSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCEtaSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCEtaSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCEtaSplitPtRecoPtPrim ) ;
+      
+      fhMCEtaSplitPtRecoPtPrimNoOverlap  = new TH2F
+      ("hMCEtaSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCEtaSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCEtaSplitPtRecoPtPrimNoOverlap ) ;
+      
+      fhMCEtaSelectedSplitPtRecoPtPrim  = new TH2F
+      ("hMCEtaSelectedSplitPtRecoPtPrim","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCEtaSelectedSplitPtRecoPtPrim ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCEtaSelectedSplitPtRecoPtPrim ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrim ) ;
+      
+      fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap  = new TH2F
+      ("hMCEtaSelectedSplitPtRecoPtPrimNoOverlap","#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, no overlap",
+       nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+      fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+      outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap ) ;
+      
+      
+      for(Int_t inlm = 0; inlm < 3; inlm++)
+      {
+        fhMCPi0PtRecoPtPrimLocMax[inlm]  = new TH2F
+        (Form("hMCPi0PtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
+         nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+        fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+        fhMCPi0PtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+        outputContainer->Add(fhMCPi0PtRecoPtPrimLocMax[inlm] ) ;
+        
+        fhMCPi0SelectedPtRecoPtPrimLocMax[inlm]  = new TH2F
+        (Form("hMCPi0SelectedPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
+         nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+        fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+        fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+        outputContainer->Add(fhMCPi0SelectedPtRecoPtPrimLocMax[inlm] ) ;
+        
+        fhMCPi0SplitPtRecoPtPrimLocMax[inlm]  = new TH2F
+        (Form("hMCPi0SplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
+         nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+        fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+        fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+        outputContainer->Add(fhMCPi0SplitPtRecoPtPrimLocMax[inlm] ) ;
+        
+        fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm]  = new TH2F
+        (Form("hMCPi0SelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
+         nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+        fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+        fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+        outputContainer->Add(fhMCPi0SelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
+        
+        fhMCEtaPtRecoPtPrimLocMax[inlm]  = new TH2F
+        (Form("hMCEtaPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
+         nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+        fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+        fhMCEtaPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+        outputContainer->Add(fhMCEtaPtRecoPtPrimLocMax[inlm] ) ;
+        
+        fhMCEtaSelectedPtRecoPtPrimLocMax[inlm]  = new TH2F
+        (Form("hMCEtaSelectedPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
+         nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+        fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+        fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+        outputContainer->Add(fhMCEtaSelectedPtRecoPtPrimLocMax[inlm] ) ;
+        
+        fhMCEtaSplitPtRecoPtPrimLocMax[inlm]  = new TH2F
+        (Form("hMCEtaSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
+         nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+        fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+        fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+        outputContainer->Add(fhMCEtaSplitPtRecoPtPrimLocMax[inlm] ) ;
+        
+        fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm]  = new TH2F
+        (Form("hMCEtaSelectedSplitPtRecoPtPrimLocMax%d",inlm+1),Form("#it{p}_{T,reco} (split sum) vs #it{p}_{T,gen}, %s",nlm[inlm].Data()),
+         nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+        fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetYTitle("#it{p}_{T,gen} (GeV/#it{c})");
+        fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ->SetXTitle("#it{p}_{T,reco} (GeV/#it{c})");
+        outputContainer->Add(fhMCEtaSelectedSplitPtRecoPtPrimLocMax[inlm] ) ;
+        
+      }
+      
       for(Int_t i = 0; i< 6; i++)
       {
-        fhMCEAsymmetry[i]  = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
-                                       Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
+        fhMCPtAsymmetry[i]  = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
+                                       Form("cluster from %s : #it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ) vs #it{E}",ptype[i].Data()),
                                        nptbins,ptmin,ptmax, 200,-1,1);
-        fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
-        fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
-        outputContainer->Add(fhMCEAsymmetry[i]);
+        fhMCPtAsymmetry[i]->SetXTitle("#it{E} (GeV)");
+        fhMCPtAsymmetry[i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
+        outputContainer->Add(fhMCPtAsymmetry[i]);
         
         fhMCSplitE[i]  = new TH1F
         (Form("hSplitE_MC%s",pname[i].Data()),
          Form("cluster from %s, energy sum of split sub-clusters",ptype[i].Data()),
          nptbins,ptmin,ptmax);
         fhMCSplitE[i]->SetYTitle("counts");
-        fhMCSplitE[i]->SetXTitle("E (GeV)");
+        fhMCSplitE[i]->SetXTitle("#it{E} (GeV)");
         outputContainer->Add(fhMCSplitE[i]) ;
         
         fhMCSplitPt[i]  = new TH1F
         (Form("hSplitPt_MC%s",pname[i].Data()),
-         Form("cluster from %s, pT sum of split sub-clusters",ptype[i].Data()),
+         Form("cluster from %s, #it{p}_{T} sum of split sub-clusters",ptype[i].Data()),
          nptbins,ptmin,ptmax);
         fhMCSplitPt[i]->SetYTitle("counts");
-        fhMCSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
+        fhMCSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhMCSplitPt[i]) ;
         
         
@@ -1814,7 +2219,7 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
          Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
         fhMCSplitPtPhi[i]->SetYTitle("#phi");
-        fhMCSplitPtPhi[i]->SetXTitle("p_{T} (GeV/c)");
+        fhMCSplitPtPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhMCSplitPtPhi[i]) ;
         
         fhMCSplitPtEta[i]  = new TH2F
@@ -1822,34 +2227,49 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
          Form("Identified as #pi^{0} (#eta), cluster from %s",
               ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
         fhMCSplitPtEta[i]->SetYTitle("#eta");
-        fhMCSplitPtEta[i]->SetXTitle("p_{T} (GeV/c)");
+        fhMCSplitPtEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhMCSplitPtEta[i]) ;
         
         
         fhMCNLocMaxSplitPt[i] = new TH2F
         (Form("hNLocMaxSplitPt_MC%s",pname[i].Data()),
-         Form("cluster from %s, pT sum of split sub-clusters, for NLM",ptype[i].Data()),
-         nptbins,ptmin,ptmax,10,0,10);
-        fhMCNLocMaxSplitPt[i] ->SetYTitle("N maxima");
-        fhMCNLocMaxSplitPt[i] ->SetXTitle("p_{T} (GeV/c)");
+         Form("cluster from %s, #it{p}_{T} sum of split sub-clusters, for NLM",ptype[i].Data()),
+         nptbins,ptmin,ptmax,20,0,20);
+        fhMCNLocMaxSplitPt[i] ->SetYTitle("#it{NLM}");
+        fhMCNLocMaxSplitPt[i] ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhMCNLocMaxSplitPt[i]) ;
         
         fhMCMassSplitPt[i]  = new TH2F
         (Form("hMassSplitPt_MC%s",pname[i].Data()),
-         Form("all pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
+         Form("all pairs #it{M}: split #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
-        fhMCMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
-        fhMCMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
+        fhMCMassSplitPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+        fhMCMassSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhMCMassSplitPt[i]) ;
         
         fhMCSelectedMassSplitPt[i]  = new TH2F
         (Form("hSelectedMassSplitPt_MC%s",pname[i].Data()),
-         Form("Selected #pi^{0} (#eta) pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
+         Form("Selected #pi^{0} (#eta) pairs #it{M}: split #it{p}_{T} vs #it{M} from %s",ptype[i].Data()),
          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
-        fhMCSelectedMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
-        fhMCSelectedMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
+        fhMCSelectedMassSplitPt[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+        fhMCSelectedMassSplitPt[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
+
+        fhMCMassSplitPtNoOverlap[i]  = new TH2F
+        (Form("hMassSplitPtNoOverlap_MC%s",pname[i].Data()),
+         Form("all pairs #it{M}: split #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
+         nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+        fhMCMassSplitPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+        fhMCMassSplitPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhMCMassSplitPtNoOverlap[i]) ;
         
+        fhMCSelectedMassSplitPtNoOverlap[i]  = new TH2F
+        (Form("hSelectedMassSplitPtNoOverlap_MC%s",pname[i].Data()),
+         Form("Selected #pi^{0} (#eta) pairs #it{M}: split #it{p}_{T} vs #it{M} from %s, no overlap",ptype[i].Data()),
+         nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+        fhMCSelectedMassSplitPtNoOverlap[i]->SetYTitle("#it{M} (GeV/#it{c}^{2})");
+        fhMCSelectedMassSplitPtNoOverlap[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhMCSelectedMassSplitPtNoOverlap[i]) ;
       }
     }
   }
@@ -1860,36 +2280,36 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
     
     for(Int_t i = 0; i< 3; i++)
     {
-      fhEAsymmetryLocMax[i]  = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
-                                        Form("Selected #pi^{0} (#eta) pairs: E vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
+      fhPtAsymmetryLocMax[i]  = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
+                                        Form("Selected #pi^{0} (#eta) pairs: #it{p}_{T} vs #it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} ), %s",nlm[i].Data()),
                                         nptbins,ptmin,ptmax,200, -1,1);
-      fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
-      fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
-      outputContainer->Add(fhEAsymmetryLocMax[i]) ;
+      fhPtAsymmetryLocMax[i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
+      fhPtAsymmetryLocMax[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhPtAsymmetryLocMax[i]) ;
     }
     
     for(Int_t ie = 0; ie< 7; ie++)
     {
       
       fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
-                                         Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
+                                         Form("#lambda_{0}^{2} vs A for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
                                          ssbins,ssmin,ssmax , 200,-1,1);
       fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
       fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
       outputContainer->Add(fhAsymmetryLambda0[ie]);
       
       fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
-                                         Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
+                                         Form("#sigma^{2}_{#eta #eta} vs #it{A} for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
                                          ssbins,ssmin,ssmax , 200,-1,1);
       fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
-      fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+      fhAsymmetryDispEta[ie]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
       outputContainer->Add(fhAsymmetryDispEta[ie]);
       
       fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
-                                         Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
+                                         Form("#sigma^{2}_{#phi #phi} vs #it{A} for %d < #it{E} < %d GeV",bin[ie],bin[ie+1]),
                                          ssbins,ssmin,ssmax , 200,-1,1);
       fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
-      fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+      fhAsymmetryDispPhi[ie]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
       outputContainer->Add(fhAsymmetryDispPhi[ie]);
     }
     
@@ -1901,24 +2321,24 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
         for(Int_t ie = 0; ie < 7; ie++)
         {
           fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
-                                                  Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
+                                                  Form("cluster from %s : #lambda_{0}^{2} vs A for %d < #it{E} < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
                                                   ssbins,ssmin,ssmax , 200,-1,1);
           fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
-          fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+          fhMCAsymmetryLambda0[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
           outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
           
           fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
-                                                  Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
+                                                  Form("cluster from %s : #sigma^{2}_{#eta #eta} vs #it{A} for %d < #it{E} < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
                                                   ssbins,ssmin,ssmax , 200,-1,1);
           fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
-          fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+          fhMCAsymmetryDispEta[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
           outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
           
           fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
-                                                  Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
+                                                  Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #it{A} for %d < #it{E} < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
                                                   ssbins,ssmin,ssmax , 200,-1,1);
           fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
-          fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+          fhMCAsymmetryDispPhi[ie][i]->SetYTitle("#it{A} = ( #it{E}_{1} - #it{E}_{2} ) / ( #it{E}_{1} + #it{E}_{2} )");
           outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
         }
       }
@@ -1933,100 +2353,100 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
     for(Int_t i = 0 ; i < 7 ; i++)
     {
       fhPtPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
-                                   Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
-      fhPtPileUp[i]->SetXTitle("p_{T} (GeV/c)");
+                                   Form("Selected #pi^{0} (#eta) #it{p}_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
+      fhPtPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
       outputContainer->Add(fhPtPileUp[i]);
       
       fhPtCellTimePileUp[i]  = new TH2F(Form("hPtCellTimePileUp%s",pileUpName[i].Data()),
                                              Form("Pt vs cell time in cluster, %s Pile-Up event",pileUpName[i].Data()),
-                                             nptbins,ptmin,ptmax,ntimebins,timemin,timemax);
-      fhPtCellTimePileUp[i]->SetXTitle("p_{T} (GeV/c)");
-      fhPtCellTimePileUp[i]->SetYTitle("t_{cell} (ns)");
+                                             nptbins,ptmin,ptmax,ntimptbins,timemin,timemax);
+      fhPtCellTimePileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      fhPtCellTimePileUp[i]->SetYTitle("#it{t}_{cell} (ns)");
       outputContainer->Add(fhPtCellTimePileUp[i]);
       
       fhPtTimeDiffPileUp[i]  = new TH2F(Form("hPtTimeDiffPileUp%s",pileUpName[i].Data()),
                                              Form("Pt vs t_{max}-t_{cell} in cluster, %s Pile-Up event",pileUpName[i].Data()),
                                              nptbins,ptmin,ptmax,400,-200,200);
-      fhPtTimeDiffPileUp[i]->SetXTitle("p_{T} (GeV/c");
-      fhPtTimeDiffPileUp[i]->SetYTitle("t_{max}-t_{cell} (ns)");
+      fhPtTimeDiffPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      fhPtTimeDiffPileUp[i]->SetYTitle("#it{t}_{max}-#it{t}_{cell} (ns)");
       outputContainer->Add(fhPtTimeDiffPileUp[i]);
 
     }
     
-    fhTimePtNoCut  = new TH2F ("hTimePt_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
-    fhTimePtNoCut->SetXTitle("p_{T} (GeV/c)");
-    fhTimePtNoCut->SetYTitle("time (ns)");
+    fhTimePtNoCut  = new TH2F ("hTimePt_NoCut","#it{t} of cluster vs #it{E} of clusters, no cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
+    fhTimePtNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhTimePtNoCut->SetYTitle("#it{t} (ns)");
     outputContainer->Add(fhTimePtNoCut);
     
-    fhTimePtSPD  = new TH2F ("hTimePt_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
-    fhTimePtSPD->SetXTitle("p_{T} (GeV/c)");
-    fhTimePtSPD->SetYTitle("time (ns)");
+    fhTimePtSPD  = new TH2F ("hTimePt_SPD","#it{t} of cluster vs #it{E} of clusters, SPD cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
+    fhTimePtSPD->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhTimePtSPD->SetYTitle("#it{t} (ns)");
     outputContainer->Add(fhTimePtSPD);
     
-    fhTimePtSPDMulti  = new TH2F ("hTimePt_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
-    fhTimePtSPDMulti->SetXTitle("p_{T} (GeV/c)");
-    fhTimePtSPDMulti->SetYTitle("time (ns)");
+    fhTimePtSPDMulti  = new TH2F ("hTimePt_SPDMulti","time of cluster vs #it{E} of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimptbins,timemin,timemax);
+    fhTimePtSPDMulti->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+    fhTimePtSPDMulti->SetYTitle("#it{t} (ns)");
     outputContainer->Add(fhTimePtSPDMulti);
     
-    fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
+    fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","#it{t} of cluster vs #it{N} pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
-    fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
+    fhTimeNPileUpVertSPD->SetXTitle("#it{t} (ns)");
     outputContainer->Add(fhTimeNPileUpVertSPD);
     
-    fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
+    fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","#it{t} of cluster vs #it{N} pile-up Tracks vertex", ntimptbins,timemin,timemax, 50,0,50 );
     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
-    fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
+    fhTimeNPileUpVertTrack->SetXTitle("#it{t} (ns)");
     outputContainer->Add(fhTimeNPileUpVertTrack);
     
-    fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
+    fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","#it{t} of cluster vs #it{N} constributors to pile-up SPD vertex", ntimptbins,timemin,timemax,50,0,50);
     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
-    fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
+    fhTimeNPileUpVertContributors->SetXTitle("#it{t} (ns)");
     outputContainer->Add(fhTimeNPileUpVertContributors);
     
-    fhTimePileUpMainVertexZDistance  = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimebins,timemin,timemax,100,0,50);
-    fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
-    fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
+    fhTimePileUpMainVertexZDistance  = new TH2F ("hTime_PileUpMainVertexZDistance","#it{t} of cluster vs distance in #it{Z} pile-up SPD vertex - main SPD vertex",ntimptbins,timemin,timemax,100,0,50);
+    fhTimePileUpMainVertexZDistance->SetYTitle("distance #it{Z} (cm) ");
+    fhTimePileUpMainVertexZDistance->SetXTitle("#it{t} (ns)");
     outputContainer->Add(fhTimePileUpMainVertexZDistance);
     
-    fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
-    fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
-    fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
+    fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","#it{t} of cluster vs distance in #it{Z} pile-up SPD vertex - z diamond",ntimptbins,timemin,timemax,100,0,50);
+    fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance #it{Z} (cm) ");
+    fhTimePileUpMainVertexZDiamond->SetXTitle("#it{t} (ns)");
     outputContainer->Add(fhTimePileUpMainVertexZDiamond);
                
-               fhPtNPileUpSPDVtx  = new TH2F ("hPt_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
+               fhPtNPileUpSPDVtx  = new TH2F ("hPt_NPileUpVertSPD","#it{p}_{T} of cluster vs #it{N} pile-up SPD vertex",
                                                                                                                                         nptbins,ptmin,ptmax,20,0,20);
                fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
-               fhPtNPileUpSPDVtx->SetXTitle("p_{T} (GeV/c)");
+               fhPtNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
                outputContainer->Add(fhPtNPileUpSPDVtx);
          
-               fhPtNPileUpTrkVtx  = new TH2F ("hPt_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
+               fhPtNPileUpTrkVtx  = new TH2F ("hPt_NPileUpVertTracks","#it{p}_{T} of cluster vs #it{N} pile-up Tracks vertex",
                                                                                                                                         nptbins,ptmin,ptmax, 20,0,20 );
                fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
-               fhPtNPileUpTrkVtx->SetXTitle("p_{T} (GeV/c)");
+               fhPtNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
                outputContainer->Add(fhPtNPileUpTrkVtx);
                
-               fhPtNPileUpSPDVtxTimeCut  = new TH2F ("hPt_NPileUpVertSPD_TimeCut","pT of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
+               fhPtNPileUpSPDVtxTimeCut  = new TH2F ("hPt_NPileUpVertSPD_TimeCut","#it{p}_{T} of cluster vs N pile-up SPD vertex, |tof| < 25 ns",
                                           nptbins,ptmin,ptmax,20,0,20);
                fhPtNPileUpSPDVtxTimeCut->SetYTitle("# vertex ");
-               fhPtNPileUpSPDVtxTimeCut->SetXTitle("p_{T} (GeV/c)");
+               fhPtNPileUpSPDVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
                outputContainer->Add(fhPtNPileUpSPDVtxTimeCut);
          
-               fhPtNPileUpTrkVtxTimeCut  = new TH2F ("hPt_NPileUpVertTracks_TimeCut","pT of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
+               fhPtNPileUpTrkVtxTimeCut  = new TH2F ("hPt_NPileUpVertTracks_TimeCut","#it{p}_{T} of cluster vs N pile-up Tracks vertex, |tof| < 25 ns",
                                                                                                                                                                        nptbins,ptmin,ptmax, 20,0,20 );
                fhPtNPileUpTrkVtxTimeCut->SetYTitle("# vertex ");
-               fhPtNPileUpTrkVtxTimeCut->SetXTitle("p_{T} (GeV/c)");
+               fhPtNPileUpTrkVtxTimeCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
                outputContainer->Add(fhPtNPileUpTrkVtxTimeCut);
     
-    fhPtNPileUpSPDVtxTimeCut2  = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","pT of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
+    fhPtNPileUpSPDVtxTimeCut2  = new TH2F ("hPt_NPileUpVertSPD_TimeCut2","#it{p}_{T} of cluster vs N pile-up SPD vertex, -25 < tof < 75 ns",
                                            nptbins,ptmin,ptmax,20,0,20);
                fhPtNPileUpSPDVtxTimeCut2->SetYTitle("# vertex ");
-               fhPtNPileUpSPDVtxTimeCut2->SetXTitle("p_{T} (GeV/c)");
+               fhPtNPileUpSPDVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
                outputContainer->Add(fhPtNPileUpSPDVtxTimeCut2);
          
-               fhPtNPileUpTrkVtxTimeCut2  = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","pT of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
+               fhPtNPileUpTrkVtxTimeCut2  = new TH2F ("hPt_NPileUpVertTracks_TimeCut2","#it{p}_{T} of cluster vs N pile-up Tracks vertex, -25 < tof < 75 ns",
                                            nptbins,ptmin,ptmax, 20,0,20 );
                fhPtNPileUpTrkVtxTimeCut2->SetYTitle("# vertex ");
-               fhPtNPileUpTrkVtxTimeCut2->SetXTitle("p_{T} (GeV/c)");
+               fhPtNPileUpTrkVtxTimeCut2->SetXTitle("#it{p}_{T} (GeV/#it{c})");
                outputContainer->Add(fhPtNPileUpTrkVtxTimeCut2);
     
   }
@@ -2063,7 +2483,7 @@ Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
     return kmcEta ;
   }//eta
   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
-            GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
+             GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
   {
     return kmcConversion ;
   }//conversion photon
@@ -2103,8 +2523,8 @@ void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
   if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)    ) ||
-     (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
-      GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay)    )
+      (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
+       GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay)    )
      )
   {
     
@@ -2191,7 +2611,7 @@ void AliAnaPi0EbE::Init()
 void AliAnaPi0EbE::InitParameters()
 {
   //Initialize the parameters of the analysis.
-  AddToHistogramsName("AnaPi0EbE_");
+  AddToHistogramsName("AnaPi0Eb#it{E}_");
   
   fInputAODGammaConvName = "PhotonsCTS" ;
   fAnaType = kIMCalo ;
@@ -2203,7 +2623,6 @@ void AliAnaPi0EbE::InitParameters()
   fNLMECutMin[0] = 10.;
   fNLMECutMin[1] = 6. ;
   fNLMECutMin[2] = 6. ;
-  
 }
 
 //__________________________________________________________________
@@ -2355,8 +2774,8 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
         //Fill some histograms about shower shape
         if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
         {
-          FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
-          FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
+          FillSelectedClusterHistograms(cluster1, mom1.Pt(), nMaxima1, photon1->GetTag());
+          FillSelectedClusterHistograms(cluster2, mom2.Pt(), nMaxima2, photon2->GetTag());
         }
         
         // Tag both photons as decay
@@ -2502,7 +2921,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
         //Fill some histograms about shower shape
         if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
         {
-          FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
+          FillSelectedClusterHistograms(cluster, mom1.Pt(), nMaxima, photon1->GetTag());
         }
         
         // Tag both photons as decay
@@ -2611,19 +3030,6 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     if(GetDebug() > 1)
       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",mom.Pt(),mom.Phi(),mom.Eta());
     
-    //Check Distance to Bad channel, set bit.
-    Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
-    if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
-    if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
-      continue ;
-    
-    if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
-    
-    //.......................................
-    // TOF cut, BE CAREFUL WITH THIS CUT
-    Double_t tof = calo->GetTOF()*1e9;
-    if(tof < fTimeCutMin || tof > fTimeCutMax) continue ;
-    
     //Play with the MC stack if available
     //Check origin of the candidates
     Int_t tag  = 0 ;
@@ -2634,57 +3040,175 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
     }
     
-    //Skip matched clusters with tracks
-    if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
+    //Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
+    
+    //Check Distance to Bad channel, set bit.
+    Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
+    if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
+    if(distBad < fMinDist){ //In bad channel (PHOS cristal size 2.2x2.2 cm)
+      //FillRejectedClusterHistograms(mom,tag,nMaxima);
+      continue ;
+    }
+    if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
+    
+    //If too low number of cells, skip it
+    if ( calo->GetNCells() < GetCaloPID()->GetClusterSplittingMinNCells())
     {
-      FillRejectedClusterHistograms(mom,tag);
+      //FillRejectedClusterHistograms(mom,tag,nMaxima);
       continue ;
     }
     
+    if(GetDebug() > 1)
+      printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: N cells cut passed %d > %d\n",
+             calo->GetNCells(), GetCaloPID()->GetClusterSplittingMinNCells());
+    
+    //.......................................
+    // TOF cut, BE CAREFUL WITH THIS CUT
+    Double_t tof = calo->GetTOF()*1e9;
+    if(tof < fTimeCutMin || tof > fTimeCutMax)
+    {
+      //FillRejectedClusterHistograms(mom,tag,nMaxima);
+      continue ;
+    }
+
     //Check PID
     //PID selection or bit setting
-    Int_t    nMaxima = 0 ;
-    Double_t mass    = 0 , angle = 0;
+    Int_t    nMaxima  = 0;
+    Double_t mass     = 0, angle    = 0;
+    Int_t    absId1   =-1, absId2   =-1;
+    Float_t  distbad1 =-1, distbad2 =-1;
+    Bool_t   fidcut1  = 0, fidcut2  = 0;
     TLorentzVector    l1, l2;
-    Int_t    absId1 = -1; Int_t absId2 = -1;
-    
+
     Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
                                                                                    GetVertex(evtIndex),nMaxima,
-                                                                                   mass,angle,l1,l2,absId1,absId2) ;
+                                                                                   mass,angle,l1,l2,absId1,absId2,
+                                                                                   distbad1,distbad2,fidcut1,fidcut2) ;
+    
     
     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
     
     
+    // Skip events where one of the new clusters (lowest energy) is close to an EMCal border or a bad channel
+    if( (fCheckSplitDistToBad) &&
+       (!fidcut2 || !fidcut1 || distbad1 < fMinDist || distbad2 < fMinDist))
+    {
+      if(GetDebug() > 1)
+        Info("MakeShowerShapeIdentification", "Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cl1 %d, cl2 %d \n",
+               calo->GetDistanceToBadChannel(),distbad1,distbad2, fidcut1,fidcut2);
+      
+      //FillRejectedClusterHistograms(mom,tag,nMaxima);
+      continue ;
+    }
+    
     //Skip events with too few or too many  NLM
     if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
     {
-      FillRejectedClusterHistograms(mom,tag);
+      //FillRejectedClusterHistograms(mom,tag,nMaxima);
       continue ;
     }
     
-    if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
-    if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
-    if(nMaxima >  2 && fNLMECutMin[2] > mom.E()) continue;
-    
     if(GetDebug() > 1)
       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
     
+    //Skip matched clusters with tracks
+    if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
+    {
+      FillRejectedClusterHistograms(mom,tag,nMaxima);
+      continue ;
+    }
+
     Float_t e1 = l1.Energy();
     Float_t e2 = l2.Energy();
     TLorentzVector l12 = l1+l2;
     Float_t ptSplit = l12.Pt();
     Float_t  eSplit = e1+e2;
-    Int_t mcIndex = GetMCIndex(tag);
+    
+    Int_t   mcIndex   =-1;
+    Int_t   noverlaps = 0;
+    Float_t ptprim    = 0;
+    if(IsDataMC())
+    {
+      mcIndex = GetMCIndex(tag);
+      
+      Bool_t ok      = kFALSE;
+      Int_t  mcLabel = calo->GetLabel();
+      
+      TLorentzVector primary = GetMCAnalysisUtils()->GetMother(mcLabel,GetReader(),ok);
+      
+      Int_t mesonLabel = -1;
+      
+      if(mcIndex == kmcPi0 || mcIndex == kmcEta)
+      {
+        if(mcIndex == kmcPi0)
+        {
+          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,111,GetReader(),ok,mesonLabel);
+          if(grandmom.E() > 0 && ok) ptprim =  grandmom.Pt();
+        }
+        else
+        {
+          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(mcLabel,221,GetReader(),ok,mesonLabel);
+          if(grandmom.E() > 0 && ok) ptprim =  grandmom.Pt();
+        }
+      }
+            
+      const UInt_t nlabels = calo->GetNLabels();
+      Int_t overpdg[nlabels];
+      noverlaps = GetMCAnalysisUtils()->GetNOverlaps(calo->GetLabels(), nlabels,tag,mesonLabel,GetReader(),overpdg);
+    }
     
     //mass of all clusters
-    fhMass       ->Fill(mom.E(),mass);
+    fhMass       ->Fill(mom.E() ,mass);
     fhMassPt     ->Fill(mom.Pt(),mass);
-    fhMassSplitPt->Fill(ptSplit,mass);
+    fhMassSplitPt->Fill(ptSplit ,mass);
     
+    Int_t indexMax = -1;
+    if     (nMaxima==1) indexMax = 0 ;
+    else if(nMaxima==2) indexMax = 1 ;
+    else                indexMax = 2 ;
+    fhMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
+
     if(IsDataMC())
     {
       fhMCMassPt[mcIndex]     ->Fill(mom.Pt(),mass);
-      fhMCMassSplitPt[mcIndex]->Fill(ptSplit,mass);
+      fhMCMassSplitPt[mcIndex]->Fill(ptSplit ,mass);
+      if(mcIndex==kmcPi0)
+      {
+        fhMCPi0PtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
+        fhMCPi0SplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
+        fhMCPi0PtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
+        fhMCPi0SplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
+
+      }
+      else if(mcIndex==kmcEta)
+      {
+        fhMCEtaPtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
+        fhMCEtaSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
+        fhMCEtaPtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
+        fhMCEtaSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
+      }
+
+      if(noverlaps==0)
+      {
+        if(mcIndex==kmcPi0)
+        {
+          fhMCPi0PtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
+          fhMCPi0SplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
+        }
+        else if(mcIndex==kmcEta)
+        {
+          fhMCEtaPtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
+          fhMCEtaSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
+        }
+        
+        fhMassNoOverlap       ->Fill(mom.E() ,mass);
+        fhMassPtNoOverlap     ->Fill(mom.Pt(),mass);
+        fhMassSplitPtNoOverlap->Fill(ptSplit ,mass);
+        
+        fhMCMassPtNoOverlap[mcIndex]     ->Fill(mom.Pt(),mass);
+        fhMCMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit ,mass);
+      }
     }
     
     // Asymmetry of all clusters
@@ -2693,29 +3217,28 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
     fhAsymmetry->Fill(mom.E(),asy);
     
-    
     if(IsDataMC())
     {
-      fhMCEAsymmetry[mcIndex]->Fill(mom.E(),asy);
+      fhMCPtAsymmetry[mcIndex]->Fill(mom.Pt(),asy);
     }
     
     // If cluster does not pass pid, not pi0/eta, skip it.
     if     (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
     {
-      if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Pi0\n");
-      FillRejectedClusterHistograms(mom,tag);
+      if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Pi0\n");
+      FillRejectedClusterHistograms(mom,tag,nMaxima);
       continue ;
     }
     
     else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
     {
-      if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Eta\n");
-      FillRejectedClusterHistograms(mom,tag);
+      if(GetDebug() > 1) Info("MakeShowerShapeIdentification","Cluster is not Eta\n");
+      FillRejectedClusterHistograms(mom,tag,nMaxima);
       continue ;
     }
     
     if(GetDebug() > 1)
-      printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
+      Info("MakeShowerShapeIdentification","Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
              mom.Pt(), idPartType);
     
     //Mass and asymmetry of selected pairs
@@ -2723,6 +3246,50 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     fhSelectedMass       ->Fill(mom.E() ,mass);
     fhSelectedMassPt     ->Fill(mom.Pt(),mass);
     fhSelectedMassSplitPt->Fill(ptSplit ,mass);
+    fhSelectedMassPtLocMax[indexMax]->Fill(mom.Pt(),mass);
+    
+    Int_t   nSM  = GetModuleNumber(calo);
+    if(nSM < GetCaloUtils()->GetNumberOfSuperModulesUsed() && nSM >=0)
+    {
+      fhSelectedMassPtLocMaxSM   [indexMax][nSM]->Fill(mom.Pt(),mass);
+      fhSelectedLambda0PtLocMaxSM[indexMax][nSM]->Fill(mom.Pt(),calo->GetM02());
+    }
+    
+    if(IsDataMC())
+    {
+      if(mcIndex==kmcPi0)
+      {
+        fhMCPi0SelectedPtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
+        fhMCPi0SelectedSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
+        fhMCPi0SelectedPtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
+        fhMCPi0SelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
+      }
+      else if(mcIndex==kmcEta)
+      {
+        fhMCEtaSelectedPtRecoPtPrim                     ->Fill(mom.Pt(),ptprim);
+        fhMCEtaSelectedSplitPtRecoPtPrim                ->Fill(ptSplit ,ptprim);
+        fhMCEtaSelectedPtRecoPtPrimLocMax     [indexMax]->Fill(mom.Pt(),ptprim);
+        fhMCEtaSelectedSplitPtRecoPtPrimLocMax[indexMax]->Fill(ptSplit ,ptprim);
+      }
+      
+      if(noverlaps==0)
+      {
+        fhSelectedMassNoOverlap       ->Fill(mom.E() ,mass);
+        fhSelectedMassPtNoOverlap     ->Fill(mom.Pt(),mass);
+        fhSelectedMassSplitPtNoOverlap->Fill(ptSplit ,mass);
+        
+        if(mcIndex==kmcPi0)
+        {
+          fhMCPi0SelectedPtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
+          fhMCPi0SelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
+        }
+        else if(mcIndex==kmcEta)
+        {
+          fhMCEtaSelectedPtRecoPtPrimNoOverlap     ->Fill(mom.Pt(),ptprim);
+          fhMCEtaSelectedSplitPtRecoPtPrimNoOverlap->Fill(ptSplit ,ptprim);
+        }
+      }
+    }
     
     fhSplitE        ->Fill( eSplit);
     fhSplitPt       ->Fill(ptSplit);
@@ -2731,7 +3298,6 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     fhSplitPtPhi    ->Fill(ptSplit,phi);
     fhSplitPtEta    ->Fill(ptSplit,mom.Eta());
     fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
-    fhNLocMaxPt     ->Fill(mom.Pt(),nMaxima);
     
     //Check split-clusters with good time window difference
     Double_t tof1  = cells->GetCellTime(absId1);
@@ -2756,11 +3322,22 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       
       fhMCSelectedMassPt     [mcIndex]->Fill(mom.Pt(),mass);
       fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
+      fhMCSelectedMassPtLocMax[mcIndex][indexMax]->Fill(mom.Pt(),mass);
+
+      if(noverlaps==0)
+      {
+        fhMCSelectedMassPtNoOverlap     [mcIndex]->Fill(mom.Pt(),mass);
+        fhMCSelectedMassSplitPtNoOverlap[mcIndex]->Fill(ptSplit,mass);
+      }
     }
     
     //-----------------------
     //Create AOD for analysis
     
+    if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
+    if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
+    if(nMaxima >  2 && fNLMECutMin[2] > mom.E()) continue;
+    
     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
     aodpi0.SetLabel(calo->GetLabel());
     
@@ -2783,7 +3360,7 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     //Fill some histograms about shower shape
     if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
     {
-      FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
+      FillSelectedClusterHistograms(calo, aodpi0.Pt(), nMaxima, tag, asy);
     }
     
     // Fill histograms to undertand pile-up before other cuts applied
@@ -2829,7 +3406,7 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
          }
       }
       else if(TMath::Abs(bc) >= 6)
-        printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Trigger BC not expected = %d\n",bc);
+        Info("MakeShowerShapeIdentification","Trigger BC not expected = %d\n",bc);
     }
     
     //Add AOD with pi0 object to aod branch
@@ -2837,7 +3414,7 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     
   }//loop
   
-  if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
+  if(GetDebug() > 1) Info("MakeShowerShapeIdentification","End fill AODs \n");
   
 }
 //______________________________________________
@@ -2847,19 +3424,17 @@ void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
   
   if(!GetOutputAODBranch())
   {
-    printf("AliAnaPi0EbE::MakeAnalysisFillHistograms()  - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
-    abort();
+    AliFatal(Form("No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data()));
   }
   //Loop on stored AOD pi0
   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
-  if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
+  if(GetDebug() > 0) Info("MakeAnalysisFillHistograms","aod branch entries %d\n", naod);
   
   Float_t cen = GetEventCentrality();
   Float_t ep  = GetEventPlaneAngle();
   
   for(Int_t iaod = 0; iaod < naod ; iaod++)
   {
-    
     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
     Int_t pdg = pi0->GetIdentifiedParticleType();
          
@@ -2875,8 +3450,6 @@ void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
     fhPt     ->Fill(pt  );
     fhE      ->Fill(ener);
     
-    fhEEta   ->Fill(ener,eta);
-    fhEPhi   ->Fill(ener,phi);
     fhPtEta  ->Fill(pt  ,eta);
     fhPtPhi  ->Fill(pt  ,phi);
     fhEtaPhi ->Fill(eta ,phi);
@@ -2887,19 +3460,19 @@ void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
     if(IsDataMC())
     {
       Int_t tag     = pi0->GetTag();
+      Int_t label   = pi0->GetLabel();
       Int_t mcIndex = GetMCIndex(tag);
       
-      fhMCE  [mcIndex] ->Fill(ener);
-      fhMCPt [mcIndex] ->Fill(pt);
-      fhMCPhi[mcIndex] ->Fill(pt,phi);
-      fhMCEta[mcIndex] ->Fill(pt,eta);
+      fhMCE    [mcIndex] ->Fill(ener);
+      fhMCPt   [mcIndex] ->Fill(pt);
+      fhMCPtPhi[mcIndex] ->Fill(pt,phi);
+      fhMCPtEta[mcIndex] ->Fill(pt,eta);
       
       fhMCPtCentrality[mcIndex]->Fill(pt,cen);
       
       if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
       {
         Float_t efracMC   = 0;
-        Int_t   label     = pi0->GetLabel();
         Int_t   momlabel  = -1;
         Bool_t  ok        = kFALSE;
         
@@ -2951,6 +3524,71 @@ void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
         
       }
       
+      if( mcIndex==kmcPi0 || mcIndex==kmcEta )
+      {
+        Float_t prodR     = -1;
+        Int_t   momindex  = -1;
+        Int_t   mompdg    = -1;
+        Int_t   momstatus = -1;
+
+        if(GetReader()->ReadStack())
+        {
+          TParticle* ancestor = GetMCStack()->Particle(label);
+          momindex  = ancestor->GetFirstMother();
+          if(momindex < 0) return;
+          TParticle* mother = GetMCStack()->Particle(momindex);
+          mompdg    = TMath::Abs(mother->GetPdgCode());
+          momstatus = mother->GetStatusCode();
+          prodR = mother->R();
+        }
+        else
+        {
+          TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
+          AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(label);
+          momindex  = ancestor->GetMother();
+          if(momindex < 0) return;
+          AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
+          mompdg    = TMath::Abs(mother->GetPdgCode());
+          momstatus = mother->GetStatus();
+          prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
+        }
+        
+        if( mcIndex==kmcPi0 )
+        {
+          fhMCPi0ProdVertex->Fill(pt,prodR);
+          if(prodR < 50)fhMCPi0ProdVertexInner->Fill(pt,prodR);
+          
+          if     (momstatus  == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
+          else if(mompdg     < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
+          else if(mompdg     > 2100  && mompdg   < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
+          else if(mompdg    == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
+          else if(mompdg    == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
+          else if(mompdg    == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
+          else if(mompdg    == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
+          else if(mompdg    >= 310   && mompdg    <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
+          else if(mompdg    == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
+          else if(momstatus == 11 || momstatus  == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
+          else                      fhMCPi0PtOrigin->Fill(pt,7.5);//other?
+        }
+        else if (mcIndex==kmcEta )
+        {
+          fhMCEtaProdVertex->Fill(pt,prodR);
+          if(prodR < 50)fhMCEtaProdVertexInner->Fill(pt,prodR);
+          
+          if     (momstatus  == 21) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
+          else if(mompdg     < 22 ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
+          else if(mompdg     > 2100  && mompdg   < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);// resonances
+          else if(mompdg    == 221) fhMCEtaPtOrigin->Fill(pt,8.5);//eta
+          else if(mompdg    == 331) fhMCEtaPtOrigin->Fill(pt,9.5);//eta prime
+          else if(mompdg    == 213) fhMCEtaPtOrigin->Fill(pt,4.5);//rho
+          else if(mompdg    == 223) fhMCEtaPtOrigin->Fill(pt,5.5);//omega
+          else if(mompdg    >= 310   && mompdg    <= 323) fhMCEtaPtOrigin->Fill(pt,6.5);//k0S, k+-,k*
+          else if(mompdg    == 130) fhMCEtaPtOrigin->Fill(pt,6.5);//k0L
+          else if(momstatus == 11 || momstatus  == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
+          else                      fhMCEtaPtOrigin->Fill(pt,7.5);//other?
+        }
+      }
+
     }//Histograms with MC
     
   }// aod loop