]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx
new method to reject events triggered by clusters from other BC
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaPi0EbE.cxx
index 6703101baa9c083cf0fea3fddcf126a9577f74b5..6fd6309a92d3b913dd3c74cdb52c91c2e356cd66 100755 (executable)
@@ -50,27 +50,42 @@ ClassImp(AliAnaPi0EbE)
 AliAnaPi0EbE::AliAnaPi0EbE() : 
     AliAnaCaloTrackCorrBaseClass(),fAnaType(kIMCalo),            fCalorimeter(""),
     fMinDist(0.),fMinDist2(0.),    fMinDist3(0.),      
-    fNLMCutMin(-1),                fNLMCutMax(10),   
-    fTimeCutMin(-10000),           fTimeCutMax(10000),      
+    fNLMCutMin(-1),                fNLMCutMax(10), 
+    fTimeCutMin(-10000),           fTimeCutMax(10000),
+    fRejectTrackMatch(kTRUE),
     fFillPileUpHistograms(0),
     fFillWeightHistograms(kFALSE), fFillTMHisto(0),              
-    fFillSelectClHisto(0),         fFillOnlySimpleSSHisto(1),
+    fFillSelectClHisto(0),         fFillOnlySimpleSSHisto(1),    fFillEMCALBCHistograms(0),
     fInputAODGammaConvName(""),
     // Histograms
     fhPt(0),                       fhE(0),                    
-    fhEEta(0),                     fhEPhi(0),                    fhEtaPhi(0),
-    fhPtDecay(0),                  fhEDecay(0),  
+    fhEEta(0),                     fhEPhi(0),                    
+    fhPtEta(0),                    fhPtPhi(0),                   fhEtaPhi(0),
+    fhEtaPhiEMCALBC0(0),           fhEtaPhiEMCALBC1(0),          fhEtaPhiEMCALBCN(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),
+    fhNLocMaxSplitPt(0),
+    fhPtDecay(0),                  fhEDecay(0),
     // Shower shape histos
     fhEDispersion(0),              fhELambda0(0),                fhELambda1(0), 
     fhELambda0NoTRD(0),            fhELambda0FracMaxCellCut(0),  
-    fhEFracMaxCell(0),             fhEFracMaxCellNoTRD(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),             fhAsymmetryE(0), 
+    fhDispEtaPhiDiffE(0),          fhSphericityE(0),           
 
     // MC histos
-    fhMCPt(),                      fhMCPhi(),                    fhMCEta(),
+    fhMCE(),                       fhMCPt(),        
+    fhMCPhi(),                     fhMCEta(),
+    fhMCEReject(),                 fhMCPtReject(),
+    fhMCPtCentrality(),            
     fhMCPi0PtGenRecoFraction(0),   fhMCEtaPtGenRecoFraction(0),
     fhMCPi0DecayPt(0),             fhMCPi0DecayPtFraction(0),      
     fhMCEtaDecayPt(0),             fhMCEtaDecayPtFraction(0),
@@ -81,10 +96,11 @@ AliAnaPi0EbE::AliAnaPi0EbE() :
     fhECellClusterRatio(0),        fhECellClusterLogRatio(0),                 
     fhEMaxCellClusterRatio(0),     fhEMaxCellClusterLogRatio(0),
     fhTrackMatchedDEta(0),         fhTrackMatchedDPhi(0),        fhTrackMatchedDEtaDPhi(0),
-    fhTrackMatchedMCParticle(0),   fhdEdx(0),                     
-    fhEOverP(0),                   fhEOverPNoTRD(0),                
+    fhTrackMatchedMCParticleE(0),
+    fhTrackMatchedMCParticleDEta(0), fhTrackMatchedMCParticleDPhi(0),
+    fhdEdx(0),                     fhEOverP(0),                 fhEOverPNoTRD(0),
     // Number of local maxima in cluster
-    fhNLocMax(0),
+    fhNLocMaxE(0),                 fhNLocMaxPt(0),
     // PileUp
     fhTimeENoCut(0),                    fhTimeESPD(0),           fhTimeESPDMulti(0),
     fhTimeNPileUpVertSPD(0),            fhTimeNPileUpVertTrack(0),
@@ -95,9 +111,19 @@ AliAnaPi0EbE::AliAnaPi0EbE() :
   
   for(Int_t i = 0; i < 6; i++)
   {
+    fhMCE              [i] = 0;
     fhMCPt             [i] = 0;
-    fhMCPhi            [i] = 0;                   
+    fhMCNLocMaxPt      [i] = 0;
+    fhMCPhi            [i] = 0;
     fhMCEta            [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;
@@ -112,6 +138,11 @@ AliAnaPi0EbE::AliAnaPi0EbE() :
     fhMCESphericity    [i] = 0;    
     fhMCEAsymmetry     [i] = 0;          
 
+    fhMCMassPt             [i]=0;
+    fhMCMassSplitPt        [i]=0;
+    fhMCSelectedMassPt     [i]=0;
+    fhMCSelectedMassSplitPt[i]=0;
+    
     for(Int_t j = 0; j < 7; j++)
     {    
       fhMCLambda0DispEta    [j][i] = 0;
@@ -131,7 +162,9 @@ AliAnaPi0EbE::AliAnaPi0EbE() :
     fhAsymmetryLambda0  [j] = 0;    
     fhAsymmetryDispEta  [j] = 0; 
     fhAsymmetryDispPhi  [j] = 0;
-  }  
+    
+    fhPtPi0PileUp       [j] = 0;
+  }
   
   for(Int_t i = 0; i < 3; i++)
   {
@@ -153,13 +186,19 @@ AliAnaPi0EbE::AliAnaPi0EbE() :
     if(i<8)fhMassPairLocMax[i] = 0;
   }
   
+  for(Int_t i = 0; i < 12; i++)
+  {
+    fhEtaPhiTriggerEMCALBC[i] = 0 ;
+    fhTimeTriggerEMCALBC  [i] = 0 ;
+  }
+  
   //Initialize parameters
   InitParameters();
   
 }
 
-//___________________________________________________________________
-void AliAnaPi0EbE::FillPileUpHistograms(Float_t energy, Float_t time) 
+//_______________________________________________________________________________
+void AliAnaPi0EbE::FillPileUpHistograms(const Float_t energy, const Float_t time) 
 {
   // Fill some histograms to understand pile-up
   if(!fFillPileUpHistograms) return;
@@ -231,6 +270,33 @@ void AliAnaPi0EbE::FillPileUpHistograms(Float_t energy, Float_t time)
   }// 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,
@@ -299,14 +365,13 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
     if (fAnaType==kSSCalo)
     {
       // Asymmetry histograms
-      fhAsymmetryE            ->Fill(e  ,asy);
       fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
       fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
       fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
     }
   }  
   
-  fhNLocMax->Fill(e,nMaxima);
+  fhNLocMaxE ->Fill(e ,nMaxima);
 
   fhELambda0LocMax   [indexMax]->Fill(e,l0); 
   fhELambda1LocMax   [indexMax]->Fill(e,l1);
@@ -379,24 +444,30 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
       
       if(IsDataMC())
       {
+        Float_t mctag = -1;
         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
         {
           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
-                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle->Fill(e, 2.5 );
-          else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle->Fill(e, 0.5 );
-          else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle->Fill(e, 1.5 );
-          else                                                                                 fhTrackMatchedMCParticle->Fill(e, 3.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 ;
           
         }
         else
         {
           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
-                     GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle->Fill(e, 6.5 );
-          else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle->Fill(e, 4.5 );
-          else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle->Fill(e, 5.5 );
-          else                                                                                 fhTrackMatchedMCParticle->Fill(e, 7.5 );
-        }        
-      }  // MC              
+                     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);
+        fhTrackMatchedMCParticleDEta->Fill(dZ, mctag);
+        fhTrackMatchedMCParticleDPhi->Fill(dR, mctag);
+        
+      }  // MC
     }
   }// Track matching histograms   
   
@@ -425,7 +496,6 @@ void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
 
       if (fAnaType==kSSCalo)
       {
-        fhMCEAsymmetry            [mcIndex]->Fill(e  ,asy);
         fhMCAsymmetryLambda0[ebin][mcIndex]->Fill(l0 ,asy);
         fhMCAsymmetryDispEta[ebin][mcIndex]->Fill(dEta,asy);
         fhMCAsymmetryDispPhi[ebin][mcIndex]->Fill(dPhi,asy);
@@ -616,6 +686,18 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   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);
+  fhPtPhi->SetYTitle("#phi (rad)");
+  fhPtPhi->SetXTitle("p_{T} (GeV/c)");
+  outputContainer->Add(fhPtPhi) ;
+  
+  fhPtEta  = new TH2F
+  ("hPtEta","Selected #pi^{0} (#eta) pairs: p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
+  fhPtEta->SetYTitle("#eta");
+  fhPtEta->SetXTitle("p_{T} (GeV/c)");
+  outputContainer->Add(fhPtEta) ;
   
   fhEtaPhi  = new TH2F
   ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax); 
@@ -623,6 +705,119 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   fhEtaPhi->SetXTitle("#eta");
   outputContainer->Add(fhEtaPhi) ; 
   
+  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);
+    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);
+    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);
+    fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
+    fhEtaPhiEMCALBCN->SetXTitle("#eta");
+    outputContainer->Add(fhEtaPhiEMCALBCN) ;
+    
+    for(Int_t i = 0; i < 12; i++)
+    {
+      fhEtaPhiTriggerEMCALBC[i] = new TH2F
+      (Form("hEtaPhiTriggerEMCALBC%d",i-5),
+       Form("cluster,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");
+      outputContainer->Add(fhEtaPhiTriggerEMCALBC[i]) ;
+      
+      fhTimeTriggerEMCALBC[i] = new TH2F
+      (Form("hTimeTriggerEMCALBC%d",i-5),
+       Form("time of cluster vs E of clusters, Trigger EMCAL-BC=%d",i-5),
+       nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
+      fhTimeTriggerEMCALBC[i]->SetXTitle("E (GeV)");
+      fhTimeTriggerEMCALBC[i]->SetYTitle("time (ns)");
+      outputContainer->Add(fhTimeTriggerEMCALBC[i]);
+      
+      fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
+      (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-5),
+       Form("time of cluster vs E of clusters, Trigger EMCAL-BC=%d",i-5),
+       nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
+      fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("E (GeV)");
+      fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("time (ns)");
+      outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
+    }
+  }
+  
+  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) ; 
+
+  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)
   {
     fhPtDecay  = new TH1F("hPtDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax); 
@@ -750,11 +945,20 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
       }
     }    
     
-    fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
-                         nptbins,ptmin,ptmax,10,0,10); 
-    fhNLocMax ->SetYTitle("N maxima");
-    fhNLocMax ->SetXTitle("E (GeV)");
-    outputContainer->Add(fhNLocMax) ;  
+    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) ;
+    
+    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) ;
+    }
     
     for (Int_t i = 0; i < 3; i++) 
     {
@@ -831,13 +1035,14 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
     
   }
   
-  if(fAnaType == kIMCalo)
-  {
+
     fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
     fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
     fhEPairDiffTime->SetYTitle("#Delta t (ns)");
     outputContainer->Add(fhEPairDiffTime);
-    
+  
+  if(fAnaType == kIMCalo)
+  {
     TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
     TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
       "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
@@ -906,23 +1111,61 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
     
     if(IsDataMC() && fFillTMHisto)
     {
-      fhTrackMatchedMCParticle  = new TH2F
-      ("hTrackMatchedMCParticle",
+      fhTrackMatchedMCParticleE  = new TH2F
+      ("hTrackMatchedMCParticleE",
        "Origin of particle vs energy",
        nptbins,ptmin,ptmax,8,0,8); 
-      fhTrackMatchedMCParticle->SetXTitle("E (GeV)");   
-      //fhTrackMatchedMCParticle->SetYTitle("Particle type");
-      
-      fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(1 ,"Photon");
-      fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(2 ,"Electron");
-      fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
-      fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(4 ,"Rest");
-      fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
-      fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
-      fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
-      fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
-      
-      outputContainer->Add(fhTrackMatchedMCParticle);   
+      fhTrackMatchedMCParticleE->SetXTitle("E (GeV)");   
+      //fhTrackMatchedMCParticleE->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");
+      
+      outputContainer->Add(fhTrackMatchedMCParticleE);
+      
+      fhTrackMatchedMCParticleDEta  = new TH2F
+      ("hTrackMatchedMCParticleDEta",
+       "Origin of particle vs #eta residual",
+       nresetabins,resetamin,resetamax,8,0,8);
+      fhTrackMatchedMCParticleDEta->SetXTitle("#Delta #eta");
+      //fhTrackMatchedMCParticleDEta->SetYTitle("Particle type");
+      
+      fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(1 ,"Photon");
+      fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(2 ,"Electron");
+      fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
+      fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(4 ,"Rest");
+      fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
+      fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
+      fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
+      fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
+      
+      outputContainer->Add(fhTrackMatchedMCParticleDEta);
+
+      fhTrackMatchedMCParticleDPhi  = new TH2F
+      ("hTrackMatchedMCParticleDPhi",
+       "Origin of particle vs #phi residual",
+       nresphibins,resphimin,resphimax,8,0,8);
+      fhTrackMatchedMCParticleDPhi->SetXTitle("#Delta #phi");
+      //fhTrackMatchedMCParticleDPhi->SetYTitle("Particle type");
+      
+      fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(1 ,"Photon");
+      fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(2 ,"Electron");
+      fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
+      fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(4 ,"Rest");
+      fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
+      fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
+      fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
+      fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
+      
+      outputContainer->Add(fhTrackMatchedMCParticleDPhi);
+
+      
     }
   }  
   
@@ -1052,6 +1295,15 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
       for(Int_t i = 0; i < 6; 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]) ; 
+        
         fhMCPt[i]  = new TH1F
         (Form("hPt_MC%s",pname[i].Data()),
          Form("Identified as #pi^{0} (#eta), cluster from %s",
@@ -1061,6 +1313,45 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
         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)
+        {
+          
+          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)");
+          outputContainer->Add(fhMCNLocMaxPt[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)");
+          outputContainer->Add(fhMCEReject[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()),
@@ -1076,7 +1367,23 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
         fhMCEta[i]->SetYTitle("#eta");
         fhMCEta[i]->SetXTitle("p_{T} (GeV/c)");
         outputContainer->Add(fhMCEta[i]) ;
+
+        fhMCMassPt[i]  = new TH2F
+        (Form("hMassPt_MC%s",pname[i].Data()),
+         Form("all pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
+         nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+        fhMCMassPt[i]->SetYTitle("mass (GeV/c^{2})");
+        fhMCMassPt[i]->SetXTitle("p_{T} (GeV/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()),
+         nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+        fhMCSelectedMassPt[i]->SetYTitle("mass (GeV/c^{2})");
+        fhMCSelectedMassPt[i]->SetXTitle("p_{T} (GeV/c)");
+        outputContainer->Add(fhMCSelectedMassPt[i]) ;
+
         
         if( fFillSelectClHisto )
         {
@@ -1194,15 +1501,143 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
     } //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);
+    
+    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)");
+    outputContainer->Add(fhSplitE) ;
+    
+    fhSplitPt  = new TH1F
+    ("hSplitPt","Selected #pi^{0} (#eta) pairs pT sum of split sub-clusters",nptbins,ptmin,ptmax);
+    fhSplitPt->SetYTitle("counts");
+    fhSplitPt->SetXTitle("p_{T} (GeV/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);
+    fhSplitPtPhi->SetYTitle("#phi (rad)");
+    fhSplitPtPhi->SetXTitle("p_{T} (GeV/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);
+    fhSplitPtEta->SetYTitle("#eta");
+    fhSplitPtEta->SetXTitle("p_{T} (GeV/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)");
+    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)");
+    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)");
+    outputContainer->Add(fhSelectedMassSplitPt) ;
+    
+
+    
+    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]);
+        
+        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)");
+        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()),
+         nptbins,ptmin,ptmax);
+        fhMCSplitPt[i]->SetYTitle("counts");
+        fhMCSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
+        outputContainer->Add(fhMCSplitPt[i]) ;
+        
+        
+        fhMCSplitPtPhi[i]  = new TH2F
+        (Form("hSplitPtPhi_MC%s",pname[i].Data()),
+         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)");
+        outputContainer->Add(fhMCSplitPtPhi[i]) ;
+        
+        fhMCSplitPtEta[i]  = new TH2F
+        (Form("hSplitPtEta_MC%s",pname[i].Data()),
+         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)");
+        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)");
+        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()),
+         nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+        fhMCMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
+        fhMCMassSplitPt[i]->SetXTitle("p_{T} (GeV/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()),
+         nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
+        fhMCSelectedMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
+        fhMCSelectedMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
+        outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
+        
+      } 
+    }
+  }
   
   if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
   {
     
-    fhAsymmetryE  = new TH2F ("hAsymmetryE","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",  
-                               nptbins,ptmin,ptmax, 200, -1,1); 
-    fhAsymmetryE->SetXTitle("E (GeV)");
-    fhAsymmetryE->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
-    outputContainer->Add(fhAsymmetryE);
     
     for(Int_t i = 0; i< 3; i++)
     {
@@ -1244,13 +1679,6 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
     {
       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]);
-        
         for(Int_t ie = 0; ie < 7; ie++)
         {
           fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
@@ -1280,6 +1708,17 @@ TList *  AliAnaPi0EbE::GetCreateOutputObjects()
   
   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)");
@@ -1386,8 +1825,8 @@ void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
   
   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 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
+  //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
   Int_t tag1 = photon1->GetTag();
   Int_t tag2 = photon2->GetTag();
   
@@ -1419,13 +1858,13 @@ void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
     {//&& (input > -1)){
       if(label1>=0)
       {
-        AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
+        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(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
+        AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
         label2 = mother2->GetMother();
         //mother2 = GetMCStack()->Particle(label2);//pi0
       }
@@ -1491,6 +1930,10 @@ void AliAnaPi0EbE::InitParameters()
   fMinDist2 = 4.;
   fMinDist3 = 5.;
   
+  fNLMECutMin[0] = 10.;
+  fNLMECutMin[1] = 6. ;
+  fNLMECutMin[2] = 6. ;
+
 }
 
 //__________________________________________________________________
@@ -1630,6 +2073,8 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
       
       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))
@@ -1657,6 +2102,9 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
         //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);        
@@ -1768,11 +2216,14 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
       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))
       {
@@ -1795,9 +2246,12 @@ 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
-        FillPileUpHistograms(mom.E(),cluster->GetTOF()*1e9);     
+        if(cluster)FillPileUpHistograms(mom.E(),cluster->GetTOF()*1e9);
         
         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
         
@@ -1900,37 +2354,93 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     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);
+    }      
+    
     //Skip matched clusters with tracks
-    if(IsTrackMatched(calo, GetReader()->GetInputEvent())) continue ;
-
+    if(fRejectTrackMatch && 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;
-    Double_t e1      = 0 , e2    = 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,e1,e2) ;   
+                                                                                   mass,angle,l1,l2,absId1,absId2) ;
     
     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) continue ;
-
+    if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) 
+    {
+      FillRejectedClusterHistograms(mom,tag);
+      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);
     
+    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);
+
+    //mass of all clusters
+    fhMass       ->Fill(mom.E(),mass);
+    fhMassPt     ->Fill(mom.Pt(),mass);
+    fhMassSplitPt->Fill(ptSplit,mass);
+
+    if(IsDataMC())
+    {
+      fhMCMassPt[mcIndex]     ->Fill(mom.Pt(),mass);
+      fhMCMassSplitPt[mcIndex]->Fill(ptSplit,mass);
+    }
+
+    // Asymmetry of all clusters
+    Float_t asy =-10;
+
+    if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
+    fhAsymmetry->Fill(mom.E(),asy);
+
+
+    if(IsDataMC())
+    {
+      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 ;
     }  
     
@@ -1938,6 +2448,47 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       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);
+    fhSelectedMassPt     ->Fill(mom.Pt(),mass);
+    fhSelectedMassSplitPt->Fill(ptSplit ,mass);
+
+    fhSplitE        ->Fill( eSplit);
+    fhSplitPt       ->Fill(ptSplit);
+    Float_t phi = mom.Phi();
+    if(phi<0) phi+=TMath::TwoPi();
+    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);
+    GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
+    tof1*=1.e9;
+    
+    Double_t tof2  = cells->GetCellTime(absId2);
+    GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
+    tof2*=1.e9;
+    
+    Double_t t12diff = tof1-tof2;
+    fhEPairDiffTime->Fill(e1+e2,    t12diff);
+    
+    if(IsDataMC())
+    {
+      fhMCSplitE        [mcIndex]->Fill( eSplit);
+      fhMCSplitPt       [mcIndex]->Fill(ptSplit);
+      fhMCSplitPtPhi    [mcIndex]->Fill(ptSplit,phi);
+      fhMCSplitPtEta    [mcIndex]->Fill(ptSplit,mom.Eta());
+      fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
+      fhMCNLocMaxPt     [mcIndex]->Fill(mom.Pt(),nMaxima);
+      
+      fhMCSelectedMassPt     [mcIndex]->Fill(mom.Pt(),mass);
+      fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
+
+    }
+    
     //-----------------------
     //Create AOD for analysis
     
@@ -1958,28 +2509,42 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     // Add number of local maxima to AOD, method name in AOD to be FIXED
     aodpi0.SetFiducialArea(nMaxima);
     
-    //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(), 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   
+    aodpi0.SetTag(tag);
     
     //Fill some histograms about shower shape
     if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
     {
-      Float_t asy =-10;      
-      if(e1+e2 > 0 ) asy = (e1-e2) / (e1+e2);
       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);
+    Double_t tofcluster   = calo->GetTOF()*1e9;
+    Double_t tofclusterUS = TMath::Abs(tofcluster);
+
+    FillPileUpHistograms(calo->E(),tofcluster);
+    if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL")
+    {
+      Float_t phicluster = aodpi0.Phi();
+      if(phicluster < 0) phicluster+=TMath::TwoPi();
+      
+      if(calo->E() > 2)
+      {
+        if      (tofclusterUS < 25) fhEtaPhiEMCALBC0->Fill(aodpi0.Eta(), phicluster);
+        else if (tofclusterUS < 75) fhEtaPhiEMCALBC1->Fill(aodpi0.Eta(), phicluster);
+        else                        fhEtaPhiEMCALBCN->Fill(aodpi0.Eta(), phicluster);
+      }
+      
+      Int_t bc = GetReader()->IsPileUpClusterTriggeredEvent();
+      if(bc > -7 && bc < 8)
+      {
+        if(calo->E() > 2) fhEtaPhiTriggerEMCALBC[bc+5]->Fill(aodpi0.Eta(), phicluster);
+        fhTimeTriggerEMCALBC[bc+5]->Fill(calo->E(), tofcluster);
+        if(GetReader()->IsPileUpFromSPD()) fhTimeTriggerEMCALBCPileUpSPD[bc+5]->Fill(calo->E(), tofcluster);
+      }
+      else printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Trigger BC not expected = %d\n",bc);
+    }
     
     //Add AOD with pi0 object to aod branch
     AddAODParticle(aodpi0);
@@ -2003,6 +2568,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++)
   {
     
@@ -2018,22 +2586,42 @@ 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);
+    fhPtEta  ->Fill(pt  ,eta);
+    fhPtPhi  ->Fill(pt  ,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())
     {
       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;