]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx
fix compilation warning
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaPi0EbE.cxx
index e9e54e6bed8cef3ab4dc6c1fd83710fabcc82f3a..1828d95d3ed72318652f9347ec4262a7cc2de8b3 100755 (executable)
@@ -40,6 +40,7 @@
 #include "AliFiducialCut.h"
 #include "TParticle.h"
 #include "AliVCluster.h"
+#include "AliESDEvent.h"
 #include "AliAODEvent.h"
 #include "AliAODMCParticle.h"
 
@@ -49,12 +50,20 @@ ClassImp(AliAnaPi0EbE)
 AliAnaPi0EbE::AliAnaPi0EbE() : 
     AliAnaCaloTrackCorrBaseClass(),fAnaType(kIMCalo),            fCalorimeter(""),
     fMinDist(0.),fMinDist2(0.),    fMinDist3(0.),      
-    fTimeCutMin(-10000),           fTimeCutMax(10000),         
-    fFillWeightHistograms(kFALSE), fFillTMHisto(0),              fFillSelectClHisto(0),
+    fNLMCutMin(-1),                fNLMCutMax(10), 
+    fTimeCutMin(-10000),           fTimeCutMax(10000),
+    fFillPileUpHistograms(0),
+    fFillWeightHistograms(kFALSE), fFillTMHisto(0),              
+    fFillSelectClHisto(0),         fFillOnlySimpleSSHisto(1),
     fInputAODGammaConvName(""),
     // Histograms
     fhPt(0),                       fhE(0),                    
     fhEEta(0),                     fhEPhi(0),                    fhEtaPhi(0),
+    fhPtCentrality(),              fhPtEventPlane(0),
+    fhPtReject(0),                 fhEReject(0),
+    fhEEtaReject(0),               fhEPhiReject(0),              fhEtaPhiReject(0),
+    fhMass(0),                     fhAsymmetry(0), 
+    fhSelectedMass(0),             fhSelectedAsymmetry(0),
     fhPtDecay(0),                  fhEDecay(0),  
     // Shower shape histos
     fhEDispersion(0),              fhELambda0(0),                fhELambda1(0), 
@@ -63,11 +72,17 @@ AliAnaPi0EbE::AliAnaPi0EbE() :
     fhENCells(0),                  fhETime(0),                   fhEPairDiffTime(0),
     fhDispEtaE(0),                 fhDispPhiE(0),
     fhSumEtaE(0),                  fhSumPhiE(0),                 fhSumEtaPhiE(0),
-    fhDispEtaPhiDiffE(0),          fhSphericityE(0),
+    fhDispEtaPhiDiffE(0),          fhSphericityE(0),           
 
     // MC histos
-    fhPtMCNo(0),                   fhPhiMCNo(0),                 fhEtaMCNo(0), 
-    fhPtMC(0),                     fhPhiMC(0),                   fhEtaMC(0),
+    fhMCE(),                       fhMCPt(),        
+    fhMCPhi(),                     fhMCEta(),
+    fhMCEReject(),                 fhMCPtReject(),
+    fhMCPtCentrality(),            
+    fhMCPi0PtGenRecoFraction(0),   fhMCEtaPtGenRecoFraction(0),
+    fhMCPi0DecayPt(0),             fhMCPi0DecayPtFraction(0),      
+    fhMCEtaDecayPt(0),             fhMCEtaDecayPtFraction(0),
+    fhMCOtherDecayPt(0),           
     fhMassPairMCPi0(0),            fhMassPairMCEta(0),
     fhAnglePairMCPi0(0),           fhAnglePairMCEta(0),
     // Weight studies
@@ -77,12 +92,23 @@ AliAnaPi0EbE::AliAnaPi0EbE() :
     fhTrackMatchedMCParticle(0),   fhdEdx(0),                     
     fhEOverP(0),                   fhEOverPNoTRD(0),                
     // Number of local maxima in cluster
-    fhNLocMax(0)
+    fhNLocMax(0),
+    // PileUp
+    fhTimeENoCut(0),                    fhTimeESPD(0),           fhTimeESPDMulti(0),
+    fhTimeNPileUpVertSPD(0),            fhTimeNPileUpVertTrack(0),
+    fhTimeNPileUpVertContributors(0),
+    fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
 {
   //default ctor
   
   for(Int_t i = 0; i < 6; i++)
   {
+    fhMCE              [i] = 0;
+    fhMCPt             [i] = 0;
+    fhMCPhi            [i] = 0;                   
+    fhMCEta            [i] = 0;
+    fhMCPtCentrality   [i] = 0;
+    
     fhEMCLambda0       [i] = 0;
     fhEMCLambda0NoTRD  [i] = 0;
     fhEMCLambda0FracMaxCellCut[i]= 0;
@@ -90,12 +116,34 @@ AliAnaPi0EbE::AliAnaPi0EbE() :
     fhEMCLambda1       [i] = 0;
     fhEMCDispersion    [i] = 0;
     
-    fhMCEDispEta       [i]  = 0;
-    fhMCEDispPhi       [i]  = 0;
-    fhMCESumEtaPhi     [i]  = 0;
-    fhMCEDispEtaPhiDiff[i]  = 0;
-    fhMCESphericity    [i]  = 0;    
+    fhMCEDispEta       [i] = 0;
+    fhMCEDispPhi       [i] = 0;
+    fhMCESumEtaPhi     [i] = 0;
+    fhMCEDispEtaPhiDiff[i] = 0;
+    fhMCESphericity    [i] = 0;    
+    fhMCEAsymmetry     [i] = 0;          
+
+    for(Int_t j = 0; j < 7; j++)
+    {    
+      fhMCLambda0DispEta    [j][i] = 0;
+      fhMCLambda0DispPhi    [j][i] = 0;
+      fhMCDispEtaDispPhi    [j][i] = 0; 
+      fhMCAsymmetryLambda0  [j][i] = 0;    
+      fhMCAsymmetryDispEta  [j][i] = 0; 
+      fhMCAsymmetryDispPhi  [j][i] = 0;
+    }
+  }
   
+  for(Int_t j = 0; j < 7; j++)
+  {  
+    fhLambda0DispEta    [j] = 0;
+    fhLambda0DispPhi    [j] = 0;
+    fhDispEtaDispPhi    [j] = 0; 
+    fhAsymmetryLambda0  [j] = 0;    
+    fhAsymmetryDispEta  [j] = 0; 
+    fhAsymmetryDispPhi  [j] = 0;
+    
+    fhPtPi0PileUp       [j] = 0;
   }
   
   for(Int_t i = 0; i < 3; i++)
@@ -108,6 +156,7 @@ AliAnaPi0EbE::AliAnaPi0EbE() :
     fhESumEtaPhiLocMax     [i] = 0;
     fhEDispEtaPhiDiffLocMax[i] = 0;
     fhESphericityLocMax    [i] = 0;
+    fhEAsymmetryLocMax     [i] = 0;
   }
   
   //Weight studies
@@ -122,10 +171,111 @@ AliAnaPi0EbE::AliAnaPi0EbE() :
   
 }
 
+//_______________________________________________________________________________
+void AliAnaPi0EbE::FillPileUpHistograms(const Float_t energy, const Float_t time) 
+{
+  // Fill some histograms to understand pile-up
+  if(!fFillPileUpHistograms) return;
+  
+  //printf("E %f, time %f\n",energy,time);
+  AliVEvent * event = GetReader()->GetInputEvent();
+  
+  fhTimeENoCut->Fill(energy,time);
+  if(GetReader()->IsPileUpFromSPD())     fhTimeESPD     ->Fill(energy,time);
+  if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
+  
+  if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
+  
+  AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
+  AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
+  
+  // N pile up vertices
+  Int_t nVerticesSPD    = -1;
+  Int_t nVerticesTracks = -1;
+  
+  if      (esdEv)
+  {
+    nVerticesSPD    = esdEv->GetNumberOfPileupVerticesSPD();
+    nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
+    
+  }//ESD
+  else if (aodEv)
+  {
+    nVerticesSPD    = aodEv->GetNumberOfPileupVerticesSPD();
+    nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
+  }//AOD
+  
+  fhTimeNPileUpVertSPD  ->Fill(time,nVerticesSPD);
+  fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
+  
+  //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n", 
+  //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
+  
+  Int_t ncont = -1;
+  Float_t z1 = -1, z2 = -1;
+  Float_t diamZ = -1;
+  for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
+  {
+    if      (esdEv)
+    {
+      const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
+      ncont=pv->GetNContributors();
+      z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
+      z2 = pv->GetZ();
+      diamZ = esdEv->GetDiamondZ();
+    }//ESD
+    else if (aodEv)
+    {
+      AliAODVertex *pv=aodEv->GetVertex(iVert);
+      if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
+      ncont=pv->GetNContributors();
+      z1=aodEv->GetPrimaryVertexSPD()->GetZ();
+      z2=pv->GetZ();
+      diamZ = aodEv->GetDiamondZ();
+    }// AOD
+    
+    Double_t distZ  = TMath::Abs(z2-z1);
+    diamZ  = TMath::Abs(z2-diamZ);
+    
+    fhTimeNPileUpVertContributors  ->Fill(time,ncont);
+    fhTimePileUpMainVertexZDistance->Fill(time,distZ);
+    fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
+    
+  }// loop
+}
+
+
+//___________________________________________________________________________________________
+void AliAnaPi0EbE::FillRejectedClusterHistograms(const TLorentzVector mom, const Int_t mctag)
+{
+ // Fill histograms that do not pass the identification (SS case only)
+  
+  Float_t ener  = mom.E();
+  Float_t pt    = mom.Pt();
+  Float_t phi   = mom.Phi();
+  if(phi < 0) phi+=TMath::TwoPi();
+  Float_t eta = mom.Eta();
+  
+  fhPtReject     ->Fill(pt);
+  fhEReject      ->Fill(ener);
+  
+  fhEEtaReject   ->Fill(ener,eta);
+  fhEPhiReject   ->Fill(ener,phi);
+  fhEtaPhiReject ->Fill(eta,phi);
+  
+  if(IsDataMC())
+  {
+    Int_t mcIndex = GetMCIndex(mctag);
+    fhMCEReject  [mcIndex] ->Fill(ener);
+    fhMCPtReject [mcIndex] ->Fill(pt);
+  }    
+}
+
 //_____________________________________________________________________________________
 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster, 
                                                  const Int_t nMaxima,
-                                                 const Int_t tag)
+                                                 const Int_t tag, 
+                                                 const Float_t asy)
 {
   // Fill shower shape, timing and other histograms for selected clusters from decay
   
@@ -135,6 +285,21 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
   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 indexMax = -1;
+  if     (nMaxima==1) indexMax = 0 ;
+  else if(nMaxima==2) indexMax = 1 ; 
+  else                indexMax = 2 ; 
+  
+  
   AliVCaloCells * cell = 0x0; 
   if(fCalorimeter == "PHOS") 
     cell = GetPHOSCells();
@@ -154,44 +319,47 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
   Float_t ll0  = 0., ll1  = 0.;
   Float_t dispp= 0., dEta = 0., dPhi    = 0.; 
   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;  
-  if(fCalorimeter == "EMCAL")
+  if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
   {
     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));
+    if(dEta+dPhi>0)fhSphericityE -> Fill(e,(dPhi-dEta)/(dEta+dPhi));
     
-    if      (e < 2 ) fhDispEtaDispPhiEBin[0]->Fill(dEta,dPhi);
-    else if (e < 4 ) fhDispEtaDispPhiEBin[1]->Fill(dEta,dPhi);
-    else if (e < 6 ) fhDispEtaDispPhiEBin[2]->Fill(dEta,dPhi);
-    else if (e < 10) fhDispEtaDispPhiEBin[3]->Fill(dEta,dPhi);
-    else             fhDispEtaDispPhiEBin[4]->Fill(dEta,dPhi);
+    fhDispEtaDispPhi[ebin]->Fill(dEta,dPhi);
+    fhLambda0DispEta[ebin]->Fill(l0  ,dEta);
+    fhLambda0DispPhi[ebin]->Fill(l0  ,dPhi);
     
+    if (fAnaType==kSSCalo)
+    {
+      // Asymmetry histograms
+      fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
+      fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
+      fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
+    }
   }  
   
   fhNLocMax->Fill(e,nMaxima);
-  Int_t indexMax = -1;
-  
-  if     (nMaxima==1) indexMax = 0 ;
-  else if(nMaxima==2) indexMax = 1 ; 
-  else                indexMax = 2 ; 
 
   fhELambda0LocMax   [indexMax]->Fill(e,l0); 
   fhELambda1LocMax   [indexMax]->Fill(e,l1);
   fhEDispersionLocMax[indexMax]->Fill(e,disp);
-  if(fCalorimeter=="EMCAL") 
+  
+  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(dEta+dPhi>0)       fhESphericityLocMax[indexMax]->Fill(e,(dPhi-dEta)/(dEta+dPhi));
+    if(fAnaType==kSSCalo) fhEAsymmetryLocMax [indexMax]->Fill(e  ,asy);
+    
   }
   
   if(fCalorimeter=="EMCAL" && nSM < 6) 
@@ -273,33 +441,7 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
   
   if(IsDataMC()) 
   {
-    Int_t mcIndex = 0;
-    
-    if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )
-    {
-      mcIndex = kmcPi0 ;      
-    }//pi0
-    else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )
-    {
-      mcIndex = kmcEta ; 
-    }//eta          
-    else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
-               GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
-    {
-      mcIndex = kmcConversion ; 
-    }//conversion photon
-    else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
-    {
-      mcIndex = kmcPhoton ; 
-    }//photon   no conversion
-    else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
-    {
-      mcIndex = kmcElectron ; 
-    }//electron
-    else 
-    {
-      mcIndex = kmcHadron ; 
-    }//other particles 
+    Int_t mcIndex = GetMCIndex(tag);
     
     fhEMCLambda0[mcIndex]    ->Fill(e, l0);
     fhEMCLambda1[mcIndex]    ->Fill(e, l1);
@@ -308,27 +450,33 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
     
     if(fCalorimeter=="EMCAL" && nSM < 6) 
       fhEMCLambda0NoTRD[mcIndex]->Fill(e, l0  );
+    
     if(maxCellFraction < 0.5) 
       fhEMCLambda0FracMaxCellCut[mcIndex]->Fill(e, l0  );  
     
-    if(fCalorimeter == "EMCAL")
+    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));  
+      if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(e,(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);
+      }
       
-      if      (e < 2 ) fhDispEtaDispPhiEBin[0]->Fill(dEta,dPhi);
-      else if (e < 4 ) fhDispEtaDispPhiEBin[1]->Fill(dEta,dPhi);
-      else if (e < 6 ) fhDispEtaDispPhiEBin[2]->Fill(dEta,dPhi);
-      else if (e < 10) fhDispEtaDispPhiEBin[3]->Fill(dEta,dPhi);
-      else             fhDispEtaDispPhiEBin[4]->Fill(dEta,dPhi);
+      fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta,dPhi);
+      fhMCLambda0DispEta[ebin][mcIndex]->Fill(l0  ,dEta);
+      fhMCLambda0DispPhi[ebin][mcIndex]->Fill(l0  ,dPhi);
       
     }
     
-    
   }//MC
+  
 }
 
 //________________________________________________________
@@ -442,95 +590,6 @@ TObjString * AliAnaPi0EbE::GetAnalysisCuts()
   return new TObjString(parList) ;
 }
 
-//__________________________________________________________________
-void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1, 
-                                       AliAODPWG4Particle * photon2, 
-                                       Int_t & label, Int_t & tag)
-{
-  // Check the labels of pare in case mother was same pi0 or eta
-  // Set the new AOD accordingly
-  
-  Int_t  label1 = photon1->GetLabel();
-  Int_t  label2 = photon2->GetLabel();
-  
-  if(label1 < 0 || label2 < 0 ) return ;
-  
-  //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
-  //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
-  Int_t tag1 = photon1->GetTag();
-  Int_t tag2 = photon2->GetTag();
-
-  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)    )
-     )
-  {
-    
-    //Check if pi0/eta mother is the same
-    if(GetReader()->ReadStack())
-    { 
-      if(label1>=0)
-      {
-        TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
-        label1 = mother1->GetFirstMother();
-        //mother1 = GetMCStack()->Particle(label1);//pi0
-      }
-      if(label2>=0)
-      {
-        TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
-        label2 = mother2->GetFirstMother();
-        //mother2 = GetMCStack()->Particle(label2);//pi0
-      }
-    } // STACK
-    else if(GetReader()->ReadAODMCParticles())
-    {//&& (input > -1)){
-      if(label1>=0)
-      {
-        AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
-        label1 = mother1->GetMother();
-        //mother1 = GetMCStack()->Particle(label1);//pi0
-      }
-      if(label2>=0)
-      {
-        AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
-        label2 = mother2->GetMother();
-        //mother2 = GetMCStack()->Particle(label2);//pi0
-      }
-    }// AOD
-    
-    //printf("mother1 %d, mother2 %d\n",label1,label2);
-    if( label1 == label2 && label1>=0 )
-    {
-      
-      label = label1;
-
-      TLorentzVector mom1 = *(photon1->Momentum());
-      TLorentzVector mom2 = *(photon2->Momentum());
-      
-      Double_t angle = mom2.Angle(mom1.Vect());
-      Double_t mass  = (mom1+mom2).M();
-      Double_t epair = (mom1+mom2).E();
-
-      if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
-      {
-        fhMassPairMCPi0 ->Fill(epair,mass);
-        fhAnglePairMCPi0->Fill(epair,angle);
-        GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
-      }
-      else 
-      {
-        fhMassPairMCEta ->Fill(epair,mass);
-        fhAnglePairMCEta->Fill(epair,angle);
-        GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
-      }
-      
-    } // same label
-  } // both from eta or pi0 decay
-  
-}   
-
 //_____________________________________________
 TList *  AliAnaPi0EbE::GetCreateOutputObjects()
 {  
@@ -565,7 +624,15 @@ 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();      
   
+  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)");
@@ -594,6 +661,59 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   fhEtaPhi->SetXTitle("#eta");
   outputContainer->Add(fhEtaPhi) ; 
   
+  fhPtCentrality  = new TH2F("hPtCentrality","centrality vs p_{T}",nptbins,ptmin,ptmax, 100,0,100);
+  fhPtCentrality->SetYTitle("centrality");
+  fhPtCentrality->SetXTitle("p_{T}(GeV/c)");
+  outputContainer->Add(fhPtCentrality) ;
+  
+  fhPtEventPlane  = new TH2F("hPtEventPlane","event plane angle vs p_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
+  fhPtEventPlane->SetYTitle("Event plane angle (rad)");
+  fhPtEventPlane->SetXTitle("p_{T} (GeV/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)");
+    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)");
+    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) ; 
+    
+    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) ; 
+    
+    fhEtaPhiReject  = new TH2F
+    ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax); 
+    fhEtaPhiReject->SetYTitle("#phi (rad)");
+    fhEtaPhiReject->SetXTitle("#eta");
+    outputContainer->Add(fhEtaPhiReject) ; 
+  }
+  
+  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)");
+  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)");
+  outputContainer->Add(fhSelectedMass) ; 
+  
   if(fAnaType != kSSCalo)
   {
     fhPtDecay  = new TH1F("hPtDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax); 
@@ -656,54 +776,68 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
       fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
       outputContainer->Add(fhEFracMaxCellNoTRD) ; 
       
-      
-      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("#sigma'^{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",  
-                             nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
-      fhSumPhiE->SetXTitle("E (GeV)");
-      fhSumPhiE->SetYTitle("#sigma'^{2}_{#phi #phi}");
-      outputContainer->Add(fhSumPhiE);  
-      
-      fhSumEtaPhiE  = new TH2F ("hSumEtaPhiE","#sigma'^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",  
-                                nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax); 
-      fhSumEtaPhiE->SetXTitle("E (GeV)");
-      fhSumEtaPhiE->SetYTitle("#sigma'^{2}_{#eta #phi}");
-      outputContainer->Add(fhSumEtaPhiE);
-      
-      fhDispEtaPhiDiffE  = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E", 
-                                     nptbins,ptmin,ptmax,200, -10,10); 
-      fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
-      fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
-      outputContainer->Add(fhDispEtaPhiDiffE);    
-      
-      fhSphericityE  = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",  
-                                 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);
-      
-      Int_t bin[] = {0,2,4,6,10,1000};
-      for(Int_t i = 0; i < 5; i++)
+      if(!fFillOnlySimpleSSHisto)
       {
-        fhDispEtaDispPhiEBin[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]), 
-                                            ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
-        fhDispEtaDispPhiEBin[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
-        fhDispEtaDispPhiEBin[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
-        outputContainer->Add(fhDispEtaDispPhiEBin[i]); 
+        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",  
+                               nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
+        fhSumPhiE->SetXTitle("E (GeV)");
+        fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
+        outputContainer->Add(fhSumPhiE);  
+        
+        fhSumEtaPhiE  = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",  
+                                  nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax); 
+        fhSumEtaPhiE->SetXTitle("E (GeV)");
+        fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
+        outputContainer->Add(fhSumEtaPhiE);
+        
+        fhDispEtaPhiDiffE  = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E", 
+                                       nptbins,ptmin,ptmax,200, -10,10); 
+        fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
+        fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
+        outputContainer->Add(fhDispEtaPhiDiffE);    
+        
+        fhSphericityE  = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",  
+                                   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);
+        
+        for(Int_t i = 0; i < 7; i++)
+        {
+          fhDispEtaDispPhi[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]), 
+                                          ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
+          fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
+          fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+          outputContainer->Add(fhDispEtaDispPhi[i]); 
+          
+          fhLambda0DispEta[i] = new TH2F (Form("hLambda0DispEta_EBin%d",i),Form("#lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]), 
+                                          ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
+          fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
+          fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
+          outputContainer->Add(fhLambda0DispEta[i]);       
+          
+          fhLambda0DispPhi[i] = new TH2F (Form("hLambda0DispPhi_EBin%d",i),Form("#lambda^{2}_{0}} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",bin[i],bin[i+1]), 
+                                          ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
+          fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
+          fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+          outputContainer->Add(fhLambda0DispPhi[i]); 
+          
+        }
       }
     }    
     
@@ -713,7 +847,6 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
     fhNLocMax ->SetXTitle("E (GeV)");
     outputContainer->Add(fhNLocMax) ;  
     
-    TString nlm[] ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
     for (Int_t i = 0; i < 3; i++) 
     {
       fhELambda0LocMax[i]  = new TH2F(Form("hELambda0LocMax%d",i+1),
@@ -737,7 +870,7 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
       fhEDispersionLocMax[i]->SetXTitle("E (GeV)");
       outputContainer->Add(fhEDispersionLocMax[i]) ; 
       
-      if(fCalorimeter == "EMCAL")
+      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()),
@@ -787,7 +920,7 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
     fhETime->SetYTitle("t (ns)");
     outputContainer->Add(fhETime);    
     
-  }// Invariant mass analysis in calorimeters and calorimeter + conversion photons
+  }
   
   if(fAnaType == kIMCalo)
   {
@@ -862,7 +995,7 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
       outputContainer->Add(fhEOverPNoTRD);   
     }   
     
-    if(IsDataMC())
+    if(IsDataMC() && fFillTMHisto)
     {
       fhTrackMatchedMCParticle  = new TH2F
       ("hTrackMatchedMCParticle",
@@ -886,7 +1019,6 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   
   if(fFillWeightHistograms)
   {
-    
     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) ");
@@ -930,43 +1062,52 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   
   if(IsDataMC()) 
   {
-    if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
-       GetReader()->GetDataType() != AliCaloTrackReader::kMC)
+    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",
+                                   nptbins,ptmin,ptmax,200,0,2); 
+      fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
+      fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / 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",
+                                   nptbins,ptmin,ptmax,200,0,2); 
+      fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
+      fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
+      outputContainer->Add(fhMCEtaPtGenRecoFraction) ; 
       
-      fhPtMC  = new TH1F("hPtMC","Identified #pi^{0} (#eta) from #pi^{0} (#eta)",nptbins,ptmin,ptmax); 
-      fhPtMC->SetYTitle("N");
-      fhPtMC->SetXTitle("p_{T} (GeV/c)");
-      outputContainer->Add(fhPtMC) ; 
+      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)");
+      outputContainer->Add(fhMCPi0DecayPt) ; 
       
-      fhPhiMC  = new TH2F
-      ("hPhiMC","Identified #pi^{0} (#eta) from #pi^{0} (#eta)",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-      fhPhiMC->SetYTitle("#phi");
-      fhPhiMC->SetXTitle("p_{T} (GeV/c)");
-      outputContainer->Add(fhPhiMC) ; 
+      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}",
+                                        nptbins,ptmin,ptmax,100,0,1); 
+      fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
+      fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
+      outputContainer->Add(fhMCPi0DecayPtFraction) ; 
       
-      fhEtaMC  = new TH2F
-      ("hEtaMC","Identified #pi^{0} (#eta) from #pi^{0} (#eta)",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-      fhEtaMC->SetYTitle("#eta");
-      fhEtaMC->SetXTitle("p_{T} (GeV/c)");
-      outputContainer->Add(fhEtaMC) ;
+      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)");
+      outputContainer->Add(fhMCEtaDecayPt) ; 
       
-      fhPtMCNo  = new TH1F("hPtMCNo","Identified #pi^{0} (#eta) not from #pi^{0} (#eta)",nptbins,ptmin,ptmax); 
-      fhPtMCNo->SetYTitle("N");
-      fhPtMCNo->SetXTitle("p_{T} (GeV/c)");
-      outputContainer->Add(fhPtMCNo) ; 
+      fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay  identified as #pi^{0} (#eta), pT 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}");
+      outputContainer->Add(fhMCEtaDecayPtFraction) ; 
       
-      fhPhiMCNo  = new TH2F
-      ("hPhiMCNo","Identified #pi^{0} (#eta) not from #pi^{0} (#eta)",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-      fhPhiMCNo->SetYTitle("#phi");
-      fhPhiMCNo->SetXTitle("p_{T} (GeV/c)");
-      outputContainer->Add(fhPhiMCNo) ; 
+      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)");
+      outputContainer->Add(fhMCOtherDecayPt) ; 
       
-      fhEtaMCNo  = new TH2F
-      ("hEtaMCNo","Identified #pi^{0} (#eta) not from #pi^{0} (#eta)",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
-      fhEtaMCNo->SetYTitle("#eta");
-      fhEtaMCNo->SetXTitle("p_{T} (GeV/c)");
-      outputContainer->Add(fhEtaMCNo) ;
+    }
+       
+    if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
+       GetReader()->GetDataType() != AliCaloTrackReader::kMC)
+    {
       
       fhAnglePairMCPi0  = new TH2F
       ("AnglePairMCPi0",
@@ -975,144 +1116,508 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
       fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
       outputContainer->Add(fhAnglePairMCPi0) ; 
 
-      fhAnglePairMCEta  = new TH2F
-      ("AnglePairMCEta",
-       "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5); 
-      fhAnglePairMCEta->SetYTitle("#alpha (rad)");
-      fhAnglePairMCEta->SetXTitle("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)");
-      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)");
-      outputContainer->Add(fhMassPairMCEta) ; 
-
-      TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"}; 
-      TString pname[] ={"Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
+      if (fAnaType!= kSSCalo)
+      {
+        fhAnglePairMCEta  = new TH2F
+        ("AnglePairMCEta",
+         "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5); 
+        fhAnglePairMCEta->SetYTitle("#alpha (rad)");
+        fhAnglePairMCEta->SetXTitle("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)");
+        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)");
+        outputContainer->Add(fhMassPairMCEta) ; 
+      }
+      
       for(Int_t i = 0; i < 6; i++)
       { 
-        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()),
-                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
-        fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
-        fhEMCLambda0[i]->SetXTitle("E (GeV)");
-        outputContainer->Add(fhEMCLambda0[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()),
-                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
-        fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
-        fhEMCLambda1[i]->SetXTitle("E (GeV)");
-        outputContainer->Add(fhEMCLambda1[i]) ; 
+        fhMCE[i]  = new TH1F
+        (Form("hE_MC%s",pname[i].Data()),
+         Form("Identified as #pi^{0} (#eta), cluster from %s",
+              ptype[i].Data()),
+         nptbins,ptmin,ptmax); 
+        fhMCE[i]->SetYTitle("N");
+        fhMCE[i]->SetXTitle("E (GeV)");
+        outputContainer->Add(fhMCE[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()),
-                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
-        fhEMCDispersion[i]->SetYTitle("D^{2}");
-        fhEMCDispersion[i]->SetXTitle("E (GeV)");
-        outputContainer->Add(fhEMCDispersion[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()),
-                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
-          fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
-          fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
-          outputContainer->Add(fhEMCLambda0NoTRD[i]) ; 
-          
-          
-          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()),
-                                       nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
-          fhMCEDispEta[i]->SetXTitle("E (GeV)");
-          fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
-          outputContainer->Add(fhMCEDispEta[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()),
-                                       nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
-          fhMCEDispPhi[i]->SetXTitle("E (GeV)");
-          fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
-          outputContainer->Add(fhMCEDispPhi[i]);  
+        fhMCPt[i]  = new TH1F
+        (Form("hPt_MC%s",pname[i].Data()),
+         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)");
+        outputContainer->Add(fhMCPt[i]) ; 
+        
+        fhMCPtCentrality[i]  = new TH2F
+        (Form("hPtCentrality_MC%s",pname[i].Data()),
+         Form("Identified as #pi^{0} (#eta), cluster from %s",
+              ptype[i].Data()),
+         nptbins,ptmin,ptmax, 100,0,100);
+        fhMCPtCentrality[i]->SetYTitle("centrality");
+        fhMCPtCentrality[i]->SetXTitle("p_{T} (GeV/c)");
+        outputContainer->Add(fhMCPtCentrality[i]) ;
+        
+        if(fAnaType == kSSCalo)
+        {
+          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)");
+          outputContainer->Add(fhMCEReject[i]) ; 
           
-          fhMCESumEtaPhi[i]  = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
-                                         Form("cluster from %s : #sigma'^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptype[i].Data()),  
-                                         nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax); 
-          fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
-          fhMCESumEtaPhi[i]->SetYTitle("#sigma'^{2}_{#eta #phi}");
-          outputContainer->Add(fhMCESumEtaPhi[i]);
+          fhMCPtReject[i]  = new TH1F
+          (Form("hPtReject_MC%s",pname[i].Data()),
+           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)");
+          outputContainer->Add(fhMCPtReject[i]) ;           
+        }
+        
+        fhMCPhi[i]  = new TH2F
+        (Form("hPhi_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]) ; 
+        
+        fhMCEta[i]  = new TH2F
+        (Form("hEta_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]) ;
+        
+        
+        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()),
+                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
+          fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
+          fhEMCLambda0[i]->SetXTitle("E (GeV)");
+          outputContainer->Add(fhEMCLambda0[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()),  
-                                              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]);    
+          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()),
+                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
+          fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
+          fhEMCLambda1[i]->SetXTitle("E (GeV)");
+          outputContainer->Add(fhEMCLambda1[i]) ; 
           
-          fhMCESphericity[i]  = new TH2F (Form("hESphericity_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]);
+          fhEMCDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
+                                         Form("Selected pair, cluster from %s : E 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]) ; 
           
-          Int_t bin[] = {0,2,4,6,10,1000};
-          for(Int_t ie = 0; ie < 5; ie++)
+          if(fCalorimeter=="EMCAL")
           {
-            fhMCDispEtaDispPhiEBin[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
+            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()),
+                                             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
+            fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
+            fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
+            outputContainer->Add(fhEMCLambda0NoTRD[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()),
+                                           nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
+              fhMCEDispEta[i]->SetXTitle("E (GeV)");
+              fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
+              outputContainer->Add(fhMCEDispEta[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()),
+                                           nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
+              fhMCEDispPhi[i]->SetXTitle("E (GeV)");
+              fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+              outputContainer->Add(fhMCEDispPhi[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()),  
+                                             nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax); 
+              fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
+              fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
+              outputContainer->Add(fhMCESumEtaPhi[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()),  
+                                                  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]);    
+              
+              fhMCESphericity[i]  = new TH2F (Form("hESphericity_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]);
+              
+              for(Int_t ie = 0; ie < 7; ie++)
+              {
+                fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
                                                       Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]), 
                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
-            fhMCDispEtaDispPhiEBin[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
-            fhMCDispEtaDispPhiEBin[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
-            outputContainer->Add(fhMCDispEtaDispPhiEBin[ie][i]); 
-          }            
-        }
-        
-        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()),
-                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
-        fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
-        fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
-        outputContainer->Add(fhEMCLambda0FracMaxCellCut[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()),
-                                        nptbins,ptmin,ptmax,100,0,1); 
-        fhEMCFracMaxCell[i]->SetYTitle("Fraction");
-        fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
-        outputContainer->Add(fhEMCFracMaxCell[i]) ;           
+                fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
+                fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+                outputContainer->Add(fhMCDispEtaDispPhi[ie][i]); 
                 
-      }//
+                fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
+                                                      Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]), 
+                                                      ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
+                fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
+                fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+                outputContainer->Add(fhMCLambda0DispEta[ie][i]);       
+                
+                fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
+                                                      Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]), 
+                                                      ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
+                fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
+                fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+                outputContainer->Add(fhMCLambda0DispPhi[ie][i]); 
+                
+              }            
+            }
+          }
+          
+          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()),
+                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
+          fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
+          fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
+          outputContainer->Add(fhEMCLambda0FracMaxCellCut[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()),
+                                          nptbins,ptmin,ptmax,100,0,1); 
+          fhEMCFracMaxCell[i]->SetYTitle("Fraction");
+          fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
+          outputContainer->Add(fhEMCFracMaxCell[i]) ;           
+          
+        }//
+      } // shower shape histo
       
     } //Not MC reader
   }//Histos with MC
   
+  if(fAnaType==kSSCalo)
+  {
+    fhAsymmetry  = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",  
+                              nptbins,ptmin,ptmax, 200, -1,1); 
+    fhAsymmetry->SetXTitle("E (GeV)");
+    fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+    outputContainer->Add(fhAsymmetry);
+    
+    fhSelectedAsymmetry  = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",  
+                             nptbins,ptmin,ptmax, 200, -1,1); 
+    fhSelectedAsymmetry->SetXTitle("E (GeV)");
+    fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+    outputContainer->Add(fhSelectedAsymmetry);
+    
+    if(IsDataMC())
+    {
+      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()),  
+                                       nptbins,ptmin,ptmax, 200,-1,1); 
+        fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
+        fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+        outputContainer->Add(fhMCEAsymmetry[i]);
+      } 
+    }
+  }
+  
+  if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
+  {
+    
+    
+    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()),
+                                        nptbins,ptmin,ptmax,200, -1,1); 
+      fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+      fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
+      outputContainer->Add(fhEAsymmetryLocMax[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]), 
+                                         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]), 
+                                         ssbins,ssmin,ssmax , 200,-1,1); 
+      fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
+      fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+      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]), 
+                                         ssbins,ssmin,ssmax , 200,-1,1); 
+      fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
+      fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+      outputContainer->Add(fhAsymmetryDispPhi[ie]);           
+    }        
+    
+    
+    if(IsDataMC()) 
+    {
+      for(Int_t i = 0; i< 6; i++)
+      {
+        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]), 
+                                                  ssbins,ssmin,ssmax , 200,-1,1); 
+          fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
+          fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+          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]), 
+                                                  ssbins,ssmin,ssmax , 200,-1,1); 
+          fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
+          fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+          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]), 
+                                                  ssbins,ssmin,ssmax , 200,-1,1); 
+          fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
+          fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
+          outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);     
+        }        
+      }
+    }
+  }
+  
+  if(fFillPileUpHistograms)
+  {
+    
+    TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
+
+    for(Int_t i = 0 ; i < 7 ; i++)
+    {
+      fhPtPi0PileUp[i]  = new TH1F(Form("hPtPi0PileUp%s",pileUpName[i].Data()),
+                                      Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
+      fhPtPi0PileUp[i]->SetXTitle("p_{T} (GeV/c)");
+      outputContainer->Add(fhPtPi0PileUp[i]);
+    }
+    
+    fhTimeENoCut  = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
+    fhTimeENoCut->SetXTitle("E (GeV)");
+    fhTimeENoCut->SetYTitle("time (ns)");
+    outputContainer->Add(fhTimeENoCut);  
+    
+    fhTimeESPD  = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
+    fhTimeESPD->SetXTitle("E (GeV)");
+    fhTimeESPD->SetYTitle("time (ns)");
+    outputContainer->Add(fhTimeESPD);  
+    
+    fhTimeESPDMulti  = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax); 
+    fhTimeESPDMulti->SetXTitle("E (GeV)");
+    fhTimeESPDMulti->SetYTitle("time (ns)");
+    outputContainer->Add(fhTimeESPDMulti);  
+    
+    fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50); 
+    fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
+    fhTimeNPileUpVertSPD->SetXTitle("time (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->SetYTitle("# vertex ");
+    fhTimeNPileUpVertTrack->SetXTitle("time (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->SetYTitle("# vertex ");
+    fhTimeNPileUpVertContributors->SetXTitle("time (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)");
+    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)");
+    outputContainer->Add(fhTimePileUpMainVertexZDiamond);  
+    
+  }
   
   //Keep neutral meson selection histograms if requiered
   //Setting done in AliNeutralMesonSelection
   
-  if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){
-    
+  if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
+  {
     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
+    
     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
+    
     delete nmsHistos;
-         
   }
   
   return outputContainer ;
   
 }
 
+//_____________________________________________
+Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
+{ 
+  
+  // Assign mc index depending on MC bit set
+  
+  if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )
+  {
+    return kmcPi0 ;      
+  }//pi0
+  else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )
+  {
+    return kmcEta ; 
+  }//eta          
+  else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
+             GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
+  {
+    return kmcConversion ; 
+  }//conversion photon
+  else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
+  {
+    return kmcPhoton ; 
+  }//photon   no conversion
+  else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
+  {
+    return kmcElectron ; 
+  }//electron
+  else 
+  {
+    return kmcHadron ; 
+  }//other particles 
+  
+}
+
+//__________________________________________________________________
+void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1, 
+                                       AliAODPWG4Particle * photon2, 
+                                       Int_t & label, Int_t & tag)
+{
+  // Check the labels of pare in case mother was same pi0 or eta
+  // Set the new AOD accordingly
+  
+  Int_t  label1 = photon1->GetLabel();
+  Int_t  label2 = photon2->GetLabel();
+  
+  if(label1 < 0 || label2 < 0 ) return ;
+  
+  //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
+  //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
+  Int_t tag1 = photon1->GetTag();
+  Int_t tag2 = photon2->GetTag();
+  
+  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)    )
+     )
+  {
+    
+    //Check if pi0/eta mother is the same
+    if(GetReader()->ReadStack())
+    { 
+      if(label1>=0)
+      {
+        TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
+        label1 = mother1->GetFirstMother();
+        //mother1 = GetMCStack()->Particle(label1);//pi0
+      }
+      if(label2>=0)
+      {
+        TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
+        label2 = mother2->GetFirstMother();
+        //mother2 = GetMCStack()->Particle(label2);//pi0
+      }
+    } // STACK
+    else if(GetReader()->ReadAODMCParticles())
+    {//&& (input > -1)){
+      if(label1>=0)
+      {
+        AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
+        label1 = mother1->GetMother();
+        //mother1 = GetMCStack()->Particle(label1);//pi0
+      }
+      if(label2>=0)
+      {
+        AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
+        label2 = mother2->GetMother();
+        //mother2 = GetMCStack()->Particle(label2);//pi0
+      }
+    }// AOD
+    
+    //printf("mother1 %d, mother2 %d\n",label1,label2);
+    if( label1 == label2 && label1>=0 )
+    {
+      
+      label = label1;
+      
+      TLorentzVector mom1 = *(photon1->Momentum());
+      TLorentzVector mom2 = *(photon2->Momentum());
+      
+      Double_t angle = mom2.Angle(mom1.Vect());
+      Double_t mass  = (mom1+mom2).M();
+      Double_t epair = (mom1+mom2).E();
+      
+      if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
+      {
+        fhMassPairMCPi0 ->Fill(epair,mass);
+        fhAnglePairMCPi0->Fill(epair,angle);
+        GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
+      }
+      else 
+      {
+        fhMassPairMCEta ->Fill(epair,mass);
+        fhAnglePairMCEta->Fill(epair,angle);
+        GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
+      }
+      
+    } // same label
+  } // both from eta or pi0 decay
+  
+}   
+
 //____________________________________________________________________________
 void AliAnaPi0EbE::Init()
 { 
@@ -1276,6 +1781,14 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
       if(nMaxima1 >  1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass); 
       if(nMaxima2 >  1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass); 
       
+      //Skip events with too few or too many  NLM
+      if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
+      
+      if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
+      
+      //Mass of all pairs
+      fhMass->Fill(epair,(mom1+mom2).M());
+      
       //Select good pair (good phi, pt cuts, aperture and invariant mass)
       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
       {
@@ -1298,9 +1811,16 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
         
         fhPtDecay->Fill(photon2->Pt());
         fhEDecay ->Fill(photon2->E() );
-
+        
         //Create AOD for analysis
         mom = mom1+mom2;
+                
+        //Mass of selected pairs
+        fhSelectedMass->Fill(epair,mom.M());
+        
+        // Fill histograms to undertand pile-up before other cuts applied
+        // Remember to relax time cuts in the reader
+        FillPileUpHistograms(mom.E(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9) /2);        
         
         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
         
@@ -1368,7 +1888,8 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
   
   Int_t nCTS  = inputAODGammaConv->GetEntriesFast();
   Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
-  if(nCTS<=0 || nCalo <=0) {
+  if(nCTS<=0 || nCalo <=0) 
+  {
     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
     return;
   }
@@ -1401,15 +1922,21 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
       else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
       else                fhMassPairLocMax[2]->Fill(epair,mass);
       
+      if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
+      if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
+
       //Play with the MC stack if available
       if(IsDataMC())
       {
         Int_t  label2 = photon2->GetLabel();
-        if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex()));
+        if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
         
         HasPairSameMCMother(photon1, photon2, label, tag) ;
       }
       
+      //Mass of selected pairs
+      fhMass->Fill(epair,(mom1+mom2).M());
+      
       //Select good pair (good phi, pt cuts, aperture and invariant mass)
       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
       {
@@ -1432,6 +1959,13 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
         
         mom = mom1+mom2;
         
+        //Mass of selected pairs
+        fhSelectedMass->Fill(epair,mom.M());
+        
+        // Fill histograms to undertand pile-up before other cuts applied
+        // Remember to relax time cuts in the reader
+        if(cluster)FillPileUpHistograms(mom.E(),cluster->GetTOF()*1e9);
+        
         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
         
         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
@@ -1497,18 +2031,18 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
     
     //Get Momentum vector, 
+    Double_t vertex[]={0,0,0};
     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
     {
       calo->GetMomentum(mom,GetVertex(evtIndex)) ;
     }//Assume that come from vertex in straight line
     else
     {
-      Double_t vertex[]={0,0,0};
       calo->GetMomentum(mom,vertex) ;
     }
          
     //If too small or big pt, skip it
-    if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; 
+    if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ; 
     
     //Check acceptance selection
     if(IsFiducialCutOn())
@@ -1517,76 +2051,130 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       if(! in ) continue ;
     }
     
-    //Create AOD for analysis
-    AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
-    aodpi0.SetLabel(calo->GetLabel());
-    
-    //Set the indeces of the original caloclusters  
-    aodpi0.SetCaloLabel(calo->GetID(),-1);
-    aodpi0.SetDetector(fCalorimeter);
     if(GetDebug() > 1) 
-      printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodpi0.Pt(),aodpi0.Phi(),aodpi0.Eta());   
-    
+      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 ;
+    if(IsDataMC())
+    {
+      tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
+      //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
+      if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
+    }      
     
-    if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
-    
-    if     (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
-    else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ; 
-    else                         aodpi0.SetDistToBad(0) ;
+    //Skip matched clusters with tracks
+    if(IsTrackMatched(calo, GetReader()->GetInputEvent())) 
+    {
+      FillRejectedClusterHistograms(mom,tag);
+      continue ;
+    }
+        
     
     //Check PID
     //PID selection or bit setting
     Int_t    nMaxima = 0 ; 
     Double_t mass    = 0 , angle = 0;
-
-    //Skip matched clusters with tracks
-    if(IsTrackMatched(calo, GetReader()->GetInputEvent())) continue ;
-      
-    // Check if cluster is pi0 via cluster splitting
-    aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
-                                                                                                 GetVertex(evtIndex),
-                                                                                                 nMaxima,mass,angle)); 
+    Double_t e1      = 0 , e2    = 0;
+    Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
+                                                                                   GetVertex(evtIndex),nMaxima,
+                                                                                   mass,angle,e1,e2) ;   
     
-    if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
+    if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
+  
+        
+    //Skip events with too few or too many  NLM
+    if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) 
+    {
+      FillRejectedClusterHistograms(mom,tag);
+      continue ;
+    }
     
-    // If cluster does not pass pid, not pi0, skip it.
-    // TO DO, add option for Eta ... or conversions
-    if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;              
+    if(GetDebug() > 1)
+      printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
     
-    if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",
-                              aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
+    //mass of all clusters
+    fhMass->Fill(mom.E(),mass);
+
+    // Asymmetry of all clusters
+    Float_t asy =-10;      
+    if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
+    fhAsymmetry->Fill(mom.E(),asy);
     
-    //Play with the MC stack if available
-    //Check origin of the candidates
-    Int_t tag  = 0 ;
     if(IsDataMC())
     {
-      if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
-          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
-        //aodpi0.SetInputFileIndex(input);
-        tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(), aodpi0.GetInputFileIndex());
-        //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
-        aodpi0.SetTag(tag);
-        if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
-      }
-    }//Work with stack also   
+      Int_t mcIndex = GetMCIndex(tag);
+      fhMCEAsymmetry[mcIndex]->Fill(mom.E(),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);
+      continue ;
+    }  
+    
+    else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)     
+    { 
+      if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Eta\n");
+      FillRejectedClusterHistograms(mom,tag);
+      continue ;
+    }  
+    
+    if(GetDebug() > 1) 
+      printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
+                              mom.Pt(), idPartType);
+    
+    //Mass and asymmetry of selected pairs
+    fhSelectedAsymmetry->Fill(mom.E(),asy);
+    fhSelectedMass     ->Fill(mom.E(),mass);
+
+    //-----------------------
+    //Create AOD for analysis
+    
+    AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
+    aodpi0.SetLabel(calo->GetLabel());
+    
+    //Set the indeces of the original caloclusters  
+    aodpi0.SetCaloLabel(calo->GetID(),-1);
+    aodpi0.SetDetector(fCalorimeter);
+
+    if     (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
+    else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ; 
+    else                         aodpi0.SetDistToBad(0) ;
+    
+    // Check if cluster is pi0 via cluster splitting
+    aodpi0.SetIdentifiedParticleType(idPartType); 
+        
+    // Add number of local maxima to AOD, method name in AOD to be FIXED
+    aodpi0.SetFiducialArea(nMaxima);
+    
+    aodpi0.SetTag(tag);
     
     //Fill some histograms about shower shape
     if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
     {
-      FillSelectedClusterHistograms(calo, nMaxima, tag);
-    }         
+      FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
+    }  
+    
+    // Fill histograms to undertand pile-up before other cuts applied
+    // Remember to relax time cuts in the reader
+    FillPileUpHistograms(calo->E(),calo->GetTOF()*1e9);
     
     //Add AOD with pi0 object to aod branch
     AddAODParticle(aodpi0);
@@ -1610,6 +2198,9 @@ void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
   
+  Float_t cen = GetEventCentrality();
+  Float_t ep  = GetEventPlaneAngle();
+  
   for(Int_t iaod = 0; iaod < naod ; iaod++)
   {
     
@@ -1625,30 +2216,94 @@ void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
     if(phi < 0) phi+=TMath::TwoPi();
     Float_t eta = pi0->Eta();
     
-    fhPt     ->Fill(pt);
+    fhPt     ->Fill(pt  );
     fhE      ->Fill(ener);
     
     fhEEta   ->Fill(ener,eta);
     fhEPhi   ->Fill(ener,phi);
-    fhEtaPhi ->Fill(eta,phi);
+    fhEtaPhi ->Fill(eta ,phi);
 
+    fhPtCentrality ->Fill(pt,cen) ;
+    fhPtEventPlane ->Fill(pt,ep ) ;
+    
+    if(fFillPileUpHistograms)
+    {
+      if(GetReader()->IsPileUpFromSPD())               fhPtPi0PileUp[0]->Fill(pt);
+      if(GetReader()->IsPileUpFromEMCal())             fhPtPi0PileUp[1]->Fill(pt);
+      if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtPi0PileUp[2]->Fill(pt);
+      if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtPi0PileUp[3]->Fill(pt);
+      if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtPi0PileUp[4]->Fill(pt);
+      if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtPi0PileUp[5]->Fill(pt);
+      if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPi0PileUp[6]->Fill(pt);
+    }
+
+    
     if(IsDataMC())
     {
-      if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) || 
-         GetReader()->GetDataType() != AliCaloTrackReader::kMC){
-        if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0))
+      Int_t tag     = pi0->GetTag();
+      Int_t mcIndex = GetMCIndex(tag);
+
+      fhMCE  [mcIndex] ->Fill(ener);
+      fhMCPt [mcIndex] ->Fill(pt);
+      fhMCPhi[mcIndex] ->Fill(pt,phi);
+      fhMCEta[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();
+        
+        Bool_t ok = kFALSE;
+        TLorentzVector mom   = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok); 
+        if(!ok) continue;
+        
+        if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
+        {
+          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok); 
+          if(grandmom.E() > 0 && ok) 
+          {
+            efracMC =  grandmom.E()/ener;
+            fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
+          }
+        }        
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
         {
-          fhPtMC  ->Fill(pt);
-          fhPhiMC ->Fill(pt,phi);
-          fhEtaMC ->Fill(pt,eta);
+          fhMCPi0DecayPt->Fill(pt);
+          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
+          if(grandmom.E() > 0 && ok) 
+          {
+            efracMC =  mom.E()/grandmom.E();
+            fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
+          }
         }
-        else
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
+        {
+          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok); 
+          if(grandmom.E() > 0 && ok) 
+          {
+            efracMC =  grandmom.E()/ener;
+            fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
+          }
+        }        
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
         {
-          fhPtMCNo  ->Fill(pt);
-          fhPhiMCNo ->Fill(pt,phi);
-          fhEtaMCNo ->Fill(pt,eta);
+          fhMCEtaDecayPt->Fill(pt);
+          TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok); 
+          if(grandmom.E() > 0 && ok) 
+          {
+            efracMC =  mom.E()/grandmom.E();
+            fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
+          }
         }
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
+        {
+          fhMCOtherDecayPt->Fill(pt);
+        }
+        
       }
+      
     }//Histograms with MC
     
   }// aod loop