AliAnaPhoton: add histograms for conversion photons identification, MC and real data...
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 Feb 2011 15:58:55 +0000 (15:58 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 8 Feb 2011 15:58:55 +0000 (15:58 +0000)
PWG4/PartCorrDep/AliAnaPhoton.cxx
PWG4/PartCorrDep/AliAnaPhoton.h
PWG4/PartCorrDep/AliAnaPi0.cxx
PWG4/PartCorrDep/AliAnaPi0.h

index 9fa87fa..9f7810b 100755 (executable)
@@ -19,6 +19,8 @@
 // Class for the photon identification.
 // Clusters from calorimeters are identified as photons
 // and kept in the AOD. Few histograms produced.
+// Produces input for other analysis classes like AliAnaPi0, 
+// AliAnaParticleHadronCorrelation ... 
 //
 // -- Author: Gustavo Conesa (LNF-INFN) 
 //////////////////////////////////////////////////////////////////////////////
@@ -31,6 +33,7 @@
 #include <TObjString.h>
 //#include <Riostream.h>
 #include "TParticle.h"
+#include "TDatabasePDG.h"
 
 // --- Analysis system --- 
 #include "AliAnaPhoton.h" 
@@ -52,8 +55,12 @@ ClassImp(AliAnaPhoton)
     fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
     fCheckConversion(kFALSE),fAddConvertedPairsToAOD(kFALSE), fMassCut(0),
     fTimeCutMin(-1), fTimeCutMax(9999999), fNCellsCut(0),
-    fhVertex(0), fhNtraNclu(0),
-    fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
+    fConvDEtaCut(2.),fConvDPhiCut(6.),
+    //fhVertex(0), 
+    fhNtraNclu(0), fhNCellsPt(0),
+    fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0), fhEtaPhiPhoton(0), fhEtaPhi05Photon(0),
+    fhPtPhotonConv(0), fhEtaPhiPhotonConv(0),fhEtaPhi05PhotonConv(0),
+    fhConvDeltaEta(0),fhConvDeltaPhi(0), fhConvDeltaEtaPhi(0), fhConvAsym(0), fhConvPt(0),
     //MC
     fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),
     fhPtMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0), 
@@ -61,9 +68,15 @@ ClassImp(AliAnaPhoton)
     fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0), 
     fhPtISR(0),fhPhiISR(0),fhEtaISR(0), 
     fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0), 
-    fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0), 
-    fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0), 
-    fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
+    fhPtOtherDecay(0),  fhPhiOtherDecay(0),  fhEtaOtherDecay(0), 
+    fhPtConversion(0),  fhPhiConversion(0),  fhEtaConversion(0),fhPtConversionTagged(0),fhEtaPhiConversion(0),fhEtaPhi05Conversion(0),
+    fhPtAntiNeutron(0), fhPhiAntiNeutron(0), fhEtaAntiNeutron(0),
+    fhPtAntiProton(0),  fhPhiAntiProton(0),  fhEtaAntiProton(0), 
+    fhPtUnknown(0),     fhPhiUnknown(0),     fhEtaUnknown(0),
+    fhConvDeltaEtaMCConversion(0),  fhConvDeltaPhiMCConversion(0),  fhConvDeltaEtaPhiMCConversion(0),  fhConvAsymMCConversion(0),  fhConvPtMCConversion(0),  fhConvDispersionMCConversion(0), fhConvM02MCConversion(0),
+    fhConvDeltaEtaMCAntiNeutron(0), fhConvDeltaPhiMCAntiNeutron(0), fhConvDeltaEtaPhiMCAntiNeutron(0), fhConvAsymMCAntiNeutron(0), fhConvPtMCAntiNeutron(0), fhConvDispersionMCAntiNeutron(0),fhConvM02MCAntiNeutron(0),
+    fhConvDeltaEtaMCAntiProton(0),  fhConvDeltaPhiMCAntiProton(0),  fhConvDeltaEtaPhiMCAntiProton(0),  fhConvAsymMCAntiProton(0),  fhConvPtMCAntiProton(0),  fhConvDispersionMCAntiProton(0), fhConvM02MCAntiProton(0),
+    fhConvDeltaEtaMCString(0),      fhConvDeltaPhiMCString(0),      fhConvDeltaEtaPhiMCString(0),      fhConvAsymMCString(0),      fhConvPtMCString(0),      fhConvDispersionMCString(0),     fhConvM02MCString(0)
 {
   //default ctor
   
@@ -97,6 +110,8 @@ TObjString *  AliAnaPhoton::GetAnalysisCuts()
   parList+=onePar ;
   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
   parList+=onePar ;  
+  snprintf(onePar,buffersize,"Conversion Selection: fConvDEtaCut %f fConvDPhiCut %f\n",fConvDEtaCut, fConvDPhiCut) ;
+  parList+=onePar ; 
   
   //Get parameters set in base class.
   parList += GetBaseParametersList() ;
@@ -130,17 +145,22 @@ TList *  AliAnaPhoton::GetCreateOutputObjects()
   Float_t etamin = GetHistoEtaMin();   
   
   //Histograms of highest Photon identified in Event
-  fhVertex  = new TH3D ("Vertex","vertex position", 20,-10.,10., 20,-10.,10., 80,-40.,40.); 
-  fhVertex->SetXTitle("X");
-  fhVertex->SetYTitle("Y");
-  fhVertex->SetZTitle("Z");
-  outputContainer->Add(fhVertex);
+//  fhVertex  = new TH3D ("Vertex","vertex position", 20,-10.,10., 20,-10.,10., 80,-40.,40.); 
+//  fhVertex->SetXTitle("X");
+//  fhVertex->SetYTitle("Y");
+//  fhVertex->SetZTitle("Z");
+//  outputContainer->Add(fhVertex);
   
   fhNtraNclu  = new TH2F ("hNtracksNcluster","# of tracks vs # of clusters", 500,0,500, 500,0,500); 
   fhNtraNclu->SetXTitle("# of tracks");
   fhNtraNclu->SetYTitle("# of clusters");
   outputContainer->Add(fhNtraNclu);
   
+  fhNCellsPt  = new TH2F ("hNCellsPt","# of tracks vs # of clusters", nptbins,ptmin, ptmax, 10,0,10); 
+  fhNCellsPt->SetXTitle("p_{T}");
+  fhNCellsPt->SetYTitle("# of cells in cluster");
+  outputContainer->Add(fhNCellsPt);  
+  
   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
   fhPtPhoton->SetYTitle("N");
   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
@@ -148,7 +168,7 @@ TList *  AliAnaPhoton::GetCreateOutputObjects()
   
   fhPhiPhoton  = new TH2F
     ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
-  fhPhiPhoton->SetYTitle("#phi");
+  fhPhiPhoton->SetYTitle("#phi (rad)");
   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
   outputContainer->Add(fhPhiPhoton) ; 
   
@@ -158,6 +178,68 @@ TList *  AliAnaPhoton::GetCreateOutputObjects()
   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
   outputContainer->Add(fhEtaPhoton) ;
   
+  fhEtaPhiPhoton  = new TH2F
+  ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
+  fhEtaPhiPhoton->SetYTitle("#phi (rad)");
+  fhEtaPhiPhoton->SetXTitle("#eta");
+  outputContainer->Add(fhEtaPhiPhoton) ;
+  
+  fhEtaPhi05Photon  = new TH2F
+  ("hEtaPhi05Photon","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
+  fhEtaPhi05Photon->SetYTitle("#phi (rad)");
+  fhEtaPhi05Photon->SetXTitle("#eta");
+  outputContainer->Add(fhEtaPhi05Photon) ;
+  
+  
+  //Conversion
+
+  fhPtPhotonConv  = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax); 
+  fhPtPhotonConv->SetYTitle("N");
+  fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
+  outputContainer->Add(fhPtPhotonConv) ; 
+  
+  fhEtaPhiPhotonConv  = new TH2F
+  ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
+  fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
+  fhEtaPhiPhotonConv->SetXTitle("#eta");
+  outputContainer->Add(fhEtaPhiPhotonConv) ;
+  
+  fhEtaPhi05PhotonConv  = new TH2F
+  ("hEtaPhi05PhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
+  fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
+  fhEtaPhi05PhotonConv->SetXTitle("#eta");
+  outputContainer->Add(fhEtaPhi05PhotonConv) ;
+  
+  fhConvDeltaEta  = new TH2F
+  ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins,-0.5,0.5); 
+  fhConvDeltaEta->SetYTitle("#Delta #eta");
+  fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
+  outputContainer->Add(fhConvDeltaEta) ;
+  
+  fhConvDeltaPhi  = new TH2F
+  ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins,-0.5,0.5); 
+  fhConvDeltaPhi->SetYTitle("#Delta #phi");
+  fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
+  outputContainer->Add(fhConvDeltaPhi) ;
+  
+  fhConvDeltaEtaPhi  = new TH2F
+  ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
+  fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
+  fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
+  outputContainer->Add(fhConvDeltaEtaPhi) ;
+  
+  fhConvAsym  = new TH2F
+  ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1); 
+  fhConvAsym->SetYTitle("Asymmetry");
+  fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
+  outputContainer->Add(fhConvAsym) ;  
+  
+  fhConvPt  = new TH2F
+  ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.); 
+  fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
+  fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
+  outputContainer->Add(fhConvPt) ;
+  
   if(IsDataMC()){
     fhDeltaE  = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50); 
     fhDeltaE->SetXTitle("#Delta E (GeV)");
@@ -304,6 +386,57 @@ TList *  AliAnaPhoton::GetCreateOutputObjects()
     fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
     outputContainer->Add(fhEtaConversion) ;
     
+    fhPtConversionTagged  = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax); 
+    fhPtConversionTagged->SetYTitle("N");
+    fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtConversionTagged) ; 
+    
+    fhEtaPhiConversion  = new TH2F
+    ("hEtaPhiConversion","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
+    fhEtaPhiConversion->SetYTitle("#phi (rad)");
+    fhEtaPhiConversion->SetXTitle("#eta");
+    outputContainer->Add(fhEtaPhiConversion) ;
+    
+    fhEtaPhi05Conversion  = new TH2F
+    ("hEtaPhi05Conversion","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax); 
+    fhEtaPhi05Conversion->SetYTitle("#phi (rad)");
+    fhEtaPhi05Conversion->SetXTitle("#eta");
+    outputContainer->Add(fhEtaPhi05Conversion) ;
+    
+    fhPtAntiNeutron  = new TH1F("hPtMCAntiNeutron","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
+    fhPtAntiNeutron->SetYTitle("N");
+    fhPtAntiNeutron->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtAntiNeutron) ; 
+    
+    fhPhiAntiNeutron  = new TH2F
+    ("hPhiMCAntiNeutron","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+    fhPhiAntiNeutron->SetYTitle("#phi");
+    fhPhiAntiNeutron->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhPhiAntiNeutron) ; 
+    
+    fhEtaAntiNeutron  = new TH2F
+    ("hEtaMCAntiNeutron","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+    fhEtaAntiNeutron->SetYTitle("#eta");
+    fhEtaAntiNeutron->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhEtaAntiNeutron) ;
+        
+    fhPtAntiProton  = new TH1F("hPtMCAntiProtonConv","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
+    fhPtAntiProton->SetYTitle("N");
+    fhPtAntiProton->SetXTitle("p_{T #gamma}(GeV/c)");
+    outputContainer->Add(fhPtAntiProton) ; 
+    
+    fhPhiAntiProton  = new TH2F
+    ("hPhiMCAntiProton","#phi_{#gamma}, unknown origin",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
+    fhPhiAntiProton->SetYTitle("#phi");
+    fhPhiAntiProton->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhPhiAntiProton) ; 
+    
+    fhEtaAntiProton  = new TH2F
+    ("hEtaMCAntiProton","#eta_{#gamma}, unknown origin",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
+    fhEtaAntiProton->SetYTitle("#eta");
+    fhEtaAntiProton->SetXTitle("p_{T #gamma} (GeV/c)");
+    outputContainer->Add(fhEtaAntiProton) ;
+    
     fhPtUnknown  = new TH1F("hPtMCUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
     fhPtUnknown->SetYTitle("N");
     fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
@@ -321,6 +454,175 @@ TList *  AliAnaPhoton::GetCreateOutputObjects()
     fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
     outputContainer->Add(fhEtaUnknown) ;
        
+    fhConvDeltaEtaMCConversion  = new TH2F
+    ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5); 
+    fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
+    fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvDeltaEtaMCConversion) ;
+    
+    fhConvDeltaPhiMCConversion  = new TH2F
+    ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5); 
+    fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
+    fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvDeltaPhiMCConversion) ;
+    
+    fhConvDeltaEtaPhiMCConversion  = new TH2F
+    ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
+    fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
+    fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
+    outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
+    
+    fhConvAsymMCConversion  = new TH2F
+    ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1); 
+    fhConvAsymMCConversion->SetYTitle("Asymmetry");
+    fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvAsymMCConversion) ;
+    
+    fhConvPtMCConversion  = new TH2F
+    ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.); 
+    fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
+    fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvPtMCConversion) ;    
+    
+    fhConvDispersionMCConversion  = new TH2F
+    ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.); 
+    fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
+    fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
+    outputContainer->Add(fhConvDispersionMCConversion) ;   
+    
+    fhConvM02MCConversion  = new TH2F
+    ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
+    fhConvM02MCConversion->SetYTitle("M02 cluster 1");
+    fhConvM02MCConversion->SetXTitle("M02 cluster 2");
+    outputContainer->Add(fhConvM02MCConversion) ;           
+    
+    fhConvDeltaEtaMCAntiNeutron  = new TH2F
+    ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5); 
+    fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
+    fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
+    
+    fhConvDeltaPhiMCAntiNeutron  = new TH2F
+    ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5); 
+    fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
+    fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
+    
+    fhConvDeltaEtaPhiMCAntiNeutron  = new TH2F
+    ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
+    fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
+    fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
+    outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;    
+    
+    fhConvAsymMCAntiNeutron  = new TH2F
+    ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1); 
+    fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
+    fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvAsymMCAntiNeutron) ;
+    
+    fhConvPtMCAntiNeutron  = new TH2F
+    ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.); 
+    fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
+    fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvPtMCAntiNeutron) ;    
+    
+    fhConvDispersionMCAntiNeutron  = new TH2F
+    ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.); 
+    fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
+    fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
+    outputContainer->Add(fhConvDispersionMCAntiNeutron) ;       
+    
+    fhConvM02MCAntiNeutron  = new TH2F
+    ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
+    fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
+    fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
+    outputContainer->Add(fhConvM02MCAntiNeutron) ;  
+    
+    fhConvDeltaEtaMCAntiProton  = new TH2F
+    ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5); 
+    fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
+    fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
+    
+    fhConvDeltaPhiMCAntiProton  = new TH2F
+    ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5); 
+    fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
+    fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
+    
+    fhConvDeltaEtaPhiMCAntiProton  = new TH2F
+    ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
+    fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
+    fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
+    outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;    
+    
+    fhConvAsymMCAntiProton  = new TH2F
+    ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1); 
+    fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
+    fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvAsymMCAntiProton) ;
+    
+    fhConvPtMCAntiProton  = new TH2F
+    ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.); 
+    fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
+    fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvPtMCAntiProton) ;
+    
+    fhConvDispersionMCAntiProton  = new TH2F
+    ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.); 
+    fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
+    fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
+    outputContainer->Add(fhConvDispersionMCAntiProton) ;       
+    
+    fhConvM02MCAntiProton  = new TH2F
+    ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
+    fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
+    fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
+    outputContainer->Add(fhConvM02MCAntiProton) ;       
+    
+    fhConvDeltaEtaMCString  = new TH2F
+    ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5); 
+    fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
+    fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvDeltaEtaMCString) ;
+    
+    fhConvDeltaPhiMCString  = new TH2F
+    ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5); 
+    fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
+    fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvDeltaPhiMCString) ;
+    
+    fhConvDeltaEtaPhiMCString  = new TH2F
+    ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5); 
+    fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
+    fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
+    outputContainer->Add(fhConvDeltaEtaPhiMCString) ;    
+    
+    fhConvAsymMCString  = new TH2F
+    ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1); 
+    fhConvAsymMCString->SetYTitle("Asymmetry");
+    fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvAsymMCString) ;
+    
+    fhConvPtMCString  = new TH2F
+    ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.); 
+    fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
+    fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
+    outputContainer->Add(fhConvPtMCString) ;
+    
+    fhConvDispersionMCString  = new TH2F
+    ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
+    fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
+    fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
+    outputContainer->Add(fhConvDispersionMCString) ;       
+    
+    fhConvM02MCString  = new TH2F
+    ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.); 
+    fhConvM02MCString->SetYTitle("M02 cluster 1");
+    fhConvM02MCString->SetXTitle("M02 cluster 2");
+    outputContainer->Add(fhConvM02MCString) ;       
+    
+    
   }//Histos with MC
     
   return outputContainer ;
@@ -352,15 +654,15 @@ void AliAnaPhoton::InitParameters()
   //Initialize the parameters of the analysis.
   AddToHistogramsName("AnaPhoton_");
 
-  fCalorimeter = "PHOS" ;
-  fMinDist  = 2.;
-  fMinDist2 = 4.;
-  fMinDist3 = 5.;
-  fMassCut  = 0.03; //30 MeV
+  fCalorimeter = "EMCAL" ;
+  fMinDist     = 2.;
+  fMinDist2    = 4.;
+  fMinDist3    = 5.;
+  fMassCut     = 0.03; //30 MeV
        
   fTimeCutMin  = -1;
   fTimeCutMax  = 9999999;
-  fNCellsCut = 0;
+  fNCellsCut   = 0;
        
   fRejectTrackMatch       = kTRUE ;
   fCheckConversion        = kFALSE;
@@ -373,11 +675,9 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
 {
   //Do photon analysis and fill aods
   
-  //  Double_t vertex2[] = {0,0,0} ; //vertex from second input aod
-  //Get the vertex and check it is not too large in z, cut for SE
+  //Get the vertex 
   Double_t v[3] = {0,0,0}; //vertex ;
   GetReader()->GetVertex(v);
-  if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;  
   
   //Select the Calorimeter of the photon
   TObjArray * pl = 0x0; 
@@ -391,14 +691,20 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
     return;
   }
 
-  //Fill AODParticle with PHOS/EMCAL aods
+  //Init arrays, variables, get number of clusters
   TLorentzVector mom, mom2 ;
   Int_t nCaloClusters = pl->GetEntriesFast();
-  if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
+  //List to be used in conversion analysis, to tag the cluster as candidate for conversion
   Bool_t * indexConverted = new Bool_t[nCaloClusters];
   for (Int_t i = 0; i < nCaloClusters; i++) 
     indexConverted[i] = kFALSE;
        
+  if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
+
+  //----------------------------------------------------
+  // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
+  //----------------------------------------------------
+  // Loop on clusters
   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){    
          
          AliVCluster * calo =  (AliVCluster*) (pl->At(icalo)); 
@@ -433,65 +739,104 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
     //    else if(input == 1) 
     //      calo->GetMomentum(mom,vertex2);//Assume that come from vertex in straight line  
     
+    //--------------------------------------
+    // Cluster selection
+    //--------------------------------------
+    if(GetDebug() > 2) 
+      printf("AliAnaPhoton::MakeAnalysisFillAOD() Current Event %d; Current File Name : %s, Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",GetReader()->GetEventNumber(),(GetReader()->GetCurrentFileName()).Data(), 
+         mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
+    //.......................................
     //If too small or big pt, skip it
     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; 
-    Double_t tof = calo->GetTOF()*1e9;
+    if(GetDebug() > 2) printf("\t Cluster %d Pass Pt Cut \n",icalo);
     
+    //.......................................
+    // TOF cut, BE CAREFUL WITH THIS CUT
+    Double_t tof = calo->GetTOF()*1e9;
     if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
-         
+         if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",icalo);
+    
+    //.......................................
     if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) continue;
-         
-    //printf("AliAnaPhoton::Current Event %d; Current File Name : %s, E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",GetReader()->GetEventNumber(),(GetReader()->GetCurrentFileName()).Data(), 
-      //     mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
+    if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",icalo);
     
+    //.......................................
     //Check acceptance selection
     if(IsFiducialCutOn()){
       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
       if(! in ) continue ;
     }
-    //printf("Fiducial cut passed \n");
-    
-    //Create AOD for analysis
-    AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
-    Int_t label = calo->GetLabel();
-    aodph.SetLabel(label);
-    //aodph.SetInputFileIndex(input);
+    if(GetDebug() > 2) printf("Fiducial cut passed \n");
     
-    //printf("Index %d, Id %d\n",icalo, calo->GetID());
-    //Set the indeces of the original caloclusters  
-    aodph.SetCaloLabel(calo->GetID(),-1);
-    aodph.SetDetector(fCalorimeter);
-    if(GetDebug() > 1) 
-      printf("AliAnaPhoton::MakeAnalysisFillAOD() - Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta()); 
+    //.......................................
+    //Skip matched clusters with tracks
+    if(fRejectTrackMatch){
+      if(IsTrackMatched(calo)) {
+        if(GetDebug() > 2) printf("\t Reject matched clusters\n");
+        continue ;
+      }
+      else  
+        if(GetDebug() > 2)  printf(" matching cut passed cut passed \n");
+    }// reject matched clusters
     
+    //.......................................
     //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)
+    if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
       continue ;
+    }
+    else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
     
-    if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
+    if(GetDebug() > 0) 
+      printf("AliAnaPhoton::MakeAnalysisFillAOD() Current Event %d; Current File Name : %s, After  selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",GetReader()->GetEventNumber(),(GetReader()->GetCurrentFileName()).Data(), 
+             mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
     
+    
+    //----------------------------
+    //Create AOD for analysis
+    //----------------------------
+    AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
+    
+    //...............................................
+    //Set the indeces of the original caloclusters (MC, ID), and calorimeter  
+    Int_t label = calo->GetLabel();
+    aodph.SetLabel(label);
+    //aodph.SetInputFileIndex(input);    
+    aodph.SetCaloLabel(calo->GetID(),-1);
+    aodph.SetDetector(fCalorimeter);
+    //printf("Index %d, Id %d\n",icalo, calo->GetID());
+
+    //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
+
+    //...............................................
+    //Set bad channel distance bit
     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; 
     else                         aodph.SetDistToBad(0) ;
     //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
-
-    //Skip matched clusters with tracks
-    if(fRejectTrackMatch && IsTrackMatched(calo)) continue ;
-    if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - TrackMatching cut passed \n");
     
-    //Check PID
+    //...............................................
+    //Set number of cells in this cluster
+    //Temporary patch FIXME
+    aodph.SetBtag(calo->GetNCells());
+    // MEFIX
+    
+    //-------------------------------------
     //PID selection or bit setting
+    //-------------------------------------
+    // MC
     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
       //Get most probable PID, check PID weights (in MC this option is mandatory)
       aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());        
       //If primary is not photon, skip it.
       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
-    }                                  
+    }  
+    //...............................................
+    // Data, PID check on
     else if(IsCaloPIDOn()){
-      
       //Get most probable PID, 2 options check PID weights 
       //or redo PID, recommended option for EMCal.             
       if(!IsCaloPIDRecalculationOn())
@@ -505,6 +850,8 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
       if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;                     
       
     }
+    //...............................................
+    // Data, PID check off
     else{
       //Set PID bits for later selection (AliAnaPi0 for example)
       //GetPDG already called in SetPIDBits.
@@ -514,36 +861,45 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
     
     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
     
+    //--------------------------------------------------------------------------------------
     //Play with the MC stack if available
+    //--------------------------------------------------------------------------------------
+
     //Check origin of the candidates
     if(IsDataMC()){
-      
       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
       if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
     }//Work with stack also   
     
-    
+    //--------------------------------------------------------------------------------------
+    // Conversions pairs analysis
     // Check if cluster comes from a conversion in the material in front of the calorimeter
     // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
-    
-    if(fCheckConversion && nCaloClusters > 1){
+    //--------------------------------------------------------------------------------------
+
+    // Do analysis only if there are more than one cluster
+    if( nCaloClusters > 1){
       Bool_t bConverted = kFALSE;
       Int_t id2 = -1;
                  
       //Check if set previously as converted couple, if so skip its use.
-      if (indexConverted[icalo]) continue;
+      if (fCheckConversion && indexConverted[icalo]) continue;
                  
+      // Second cluster loop
       for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
         //Check if set previously as converted couple, if so skip its use.
         if (indexConverted[jcalo]) continue;
         //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
         AliVCluster * calo2 =  (AliVCluster*) (pl->At(jcalo));              //Get cluster kinematics
+        
+        //Mixed event, get index of event
         Int_t evtIndex2 = 0 ; 
         if (GetMixedEvent()) {
           evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ; 
           
-        }        
-
+        }      
+        
+        //Get kinematics of second cluster
         if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
           calo2->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
         else{
@@ -556,23 +912,108 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
         if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
         if(!in2) continue;      
         
+        //................................................
         //Get mass of pair, if small, take this pair.
         //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut);
         if((mom+mom2).M() < fMassCut){  
-          bConverted = kTRUE;
+          
           id2 = calo2->GetID();
+          indexConverted[icalo]=kTRUE;
           indexConverted[jcalo]=kTRUE;
+          
           if(GetDebug() > 2)
             printf("AliAnaPhoton::MakeAnalysisFillAOD(): Pair with mass %f < %2.4f; \n    cluster1 id %d, e %2.3f  SM %d, eta %2.3f, phi %2.3f ; \n    cluster2 id %d, e %2.3f, SM %d,eta %2.3f, phi %2.3f\n",
                    (mom+mom2).M(),fMassCut,
                    calo->GetID(),calo->E(),GetCaloUtils()->GetModuleNumber(calo), mom.Eta(), mom.Phi(),
                    id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2),mom2.Eta(), mom2.Phi());
+          
+          //...............................................
+          //Fill few histograms with kinematics of the pair
+          //FIXME, move all this to MakeAnalysisFillHistograms ...
+          fhConvDeltaEta   ->Fill( (mom+mom2).M(), mom.Eta()-mom2.Eta() );
+          fhConvDeltaPhi   ->Fill( (mom+mom2).M(), mom.Phi()-mom2.Phi() );
+          fhConvAsym       ->Fill( (mom+mom2).M(), TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
+          fhConvDeltaEtaPhi->Fill( mom.Eta()-mom2.Eta(), mom.Phi()-mom2.Phi() );
+          fhConvPt         ->Fill( (mom+mom2).M(), (mom+mom2).Pt());
+          
+          //...............................................
+          //Select pairs in a eta-phi window
+          if(TMath::Abs(mom.Eta()-mom2.Eta()) < fConvDEtaCut && TMath::Abs(mom.Phi()-mom2.Phi()) < fConvDPhiCut )
+            bConverted = kTRUE;          
+          
+          //...........................................
+          //Fill more histograms, simulated data
+          //FIXME, move all this to MakeAnalysisFillHistograms ...
+          if(IsDataMC()){
+            
+            //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
+            Int_t ancPDG    = 0;
+            Int_t ancStatus = 0;
+            TLorentzVector momentum;
+            Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(calo->GetLabel(), calo2->GetLabel(), 
+                                                                        GetReader(), ancPDG, ancStatus, momentum);
+            
+            // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
+            //                          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
+            
+            Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(calo2->GetLabels(),calo2->GetNLabels(),GetReader(), 0);
+            if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCConversion)){
+              if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
+                fhConvDeltaEtaMCConversion   ->Fill( (mom+mom2).M(), mom.Eta()-mom2.Eta());
+                fhConvDeltaPhiMCConversion   ->Fill( (mom+mom2).M(), mom.Phi()-mom2.Phi());
+                fhConvAsymMCConversion       ->Fill( (mom+mom2).M(), TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
+                fhConvDeltaEtaPhiMCConversion->Fill( mom.Eta()-mom2.Eta(), mom.Phi()-mom2.Phi() );
+                fhConvPtMCConversion         ->Fill( (mom+mom2).M(), (mom+mom2).Pt());
+                fhConvDispersionMCConversion ->Fill( calo->GetDispersion(), calo2->GetDispersion());
+                fhConvM02MCConversion        ->Fill( calo->GetM02(), calo2->GetM02());
+
+              }              
+            }
+            else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiNeutron)){
+              if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
+                fhConvDeltaEtaMCAntiNeutron    ->Fill( (mom+mom2).M(), mom.Eta()-mom2.Eta());
+                fhConvDeltaPhiMCAntiNeutron    ->Fill( (mom+mom2).M(), mom.Phi()-mom2.Phi());
+                fhConvAsymMCAntiNeutron        ->Fill( (mom+mom2).M(), TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
+                fhConvDeltaEtaPhiMCAntiNeutron ->Fill( mom.Eta()-mom2.Eta(), mom.Phi()-mom2.Phi() );
+                fhConvPtMCAntiNeutron          ->Fill( (mom+mom2).M(), (mom+mom2).Pt());
+                fhConvDispersionMCAntiNeutron  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
+                fhConvM02MCAntiNeutron         ->Fill( calo->GetM02(), calo2->GetM02());
+              }
+            }
+            else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiProton)){
+              if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
+                fhConvDeltaEtaMCAntiProton    ->Fill( (mom+mom2).M(), mom.Eta()-mom2.Eta());
+                fhConvDeltaPhiMCAntiProton    ->Fill( (mom+mom2).M(), mom.Phi()-mom2.Phi());
+                fhConvAsymMCAntiProton        ->Fill( (mom+mom2).M(), TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
+                fhConvDeltaEtaPhiMCAntiProton ->Fill( mom.Eta()-mom2.Eta(), mom.Phi()-mom2.Phi() );
+                fhConvPtMCAntiProton          ->Fill( (mom+mom2).M(), (mom+mom2).Pt());
+                fhConvDispersionMCAntiProton  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
+                fhConvM02MCAntiProton         ->Fill( calo->GetM02(), calo2->GetM02());
+              }
+            }
+            
+            //Pairs coming from fragmenting pairs.
+            if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
+              fhConvDeltaEtaMCString    ->Fill( (mom+mom2).M(), mom.Eta()-mom2.Eta());
+              fhConvDeltaPhiMCString    ->Fill( (mom+mom2).M(), mom.Phi()-mom2.Phi());
+              fhConvAsymMCString        ->Fill( (mom+mom2).M(), TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
+              fhConvDeltaEtaPhiMCString ->Fill( mom.Eta()-mom2.Eta(), mom.Phi()-mom2.Phi() );
+              fhConvPtMCString          ->Fill( (mom+mom2).M(), (mom+mom2).Pt());
+              fhConvDispersionMCString  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
+              fhConvM02MCString         ->Fill( calo->GetM02(), calo2->GetM02());
+            }
+            
+          }// Data MC
+
           break;
         }
                          
       }//Mass loop
                  
+      //..........................................................................................................
+      //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
       if(bConverted){ 
+        //Add to AOD
         if(fAddConvertedPairsToAOD){
           //Create AOD of pair analysis
           TLorentzVector mpair = mom+mom2;
@@ -586,18 +1027,26 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
           aodpair.SetDetector(fCalorimeter);
           aodpair.SetPdg(aodph.GetPdg());
           aodpair.SetTag(aodph.GetTag());
-          
+          aodpair.SetTagged(kTRUE);
           //Add AOD with pair object to aod branch
           AddAODParticle(aodpair);
           //printf("\t \t both added pair\n");
         }
         
         //Do not add the current calocluster
-        continue;
+        if(fCheckConversion) continue;
+        else {
+          //printf("TAGGED\n");
+          //Tag this cluster as likely conversion
+          aodph.SetTagged(kTRUE);
+        }
       }//converted pair
     }//check conversion
     //printf("\t \t added single cluster %d\n",icalo);
          
+    //FIXME, this to MakeAnalysisFillHistograms ...
+    fhNCellsPt->Fill(aodph.Pt(),calo->GetNCells());
+    
     //Add AOD with photon object to aod branch
     AddAODParticle(aodph);
     
@@ -612,8 +1061,9 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
 //__________________________________________________________________
 void  AliAnaPhoton::MakeAnalysisFillHistograms() 
 {
-  //Do analysis and fill histograms
+  //Fill histograms
   
+  //-------------------------------------------------------------------
        // Access MC information in stack if requested, check that it exists.   
        AliStack * stack = 0x0;
        TParticle * primary = 0x0;   
@@ -647,11 +1097,15 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
                }
        }// is data and MC
        
-       //Loop on stored AOD photons
+  
+  // Get vertex
   Double_t v[3] = {0,0,0}; //vertex ;
   GetReader()->GetVertex(v);
-  fhVertex->Fill(v[0],v[1],v[2]);  
-  if(TMath::Abs(v[2]) > GetZvertexCut()) return ;  
+  //fhVertex->Fill(v[0],v[1],v[2]);  
+  if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
+  
+  //----------------------------------
+       //Loop on stored AOD photons
        Int_t naod = GetOutputAODBranch()->GetEntriesFast();
   fhNtraNclu->Fill(GetReader()->GetTrackMultiplicity(), naod);
        if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
@@ -670,6 +1124,7 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
          if(GetDebug() > 2) 
            printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
          
+    //................................
          //Fill photon histograms 
          Float_t ptcluster  = ph->Pt();
          Float_t phicluster = ph->Phi();
@@ -677,9 +1132,19 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
          Float_t ecluster   = ph->E();
          
          fhPtPhoton  ->Fill(ptcluster);
+    if(ph->IsTagged())fhPtPhotonConv->Fill(ptcluster);
          fhPhiPhoton ->Fill(ptcluster,phicluster);
          fhEtaPhoton ->Fill(ptcluster,etacluster);
-         
+    if(ptcluster > 0.5){
+      fhEtaPhiPhoton ->Fill(etacluster, phicluster);
+      if(ph->IsTagged())fhEtaPhiPhotonConv->Fill(etacluster, phicluster);
+    }
+    else {
+      fhEtaPhi05Photon ->Fill(etacluster, phicluster);
+      if(ph->IsTagged())fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
+    }
+
+    //.......................................
          //Play with the MC data if available
          if(IsDataMC()){
            
@@ -696,6 +1161,9 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
           fhPtConversion  ->Fill(ptcluster);
           fhPhiConversion ->Fill(ptcluster,phicluster);
           fhEtaConversion ->Fill(ptcluster,etacluster);
+          if(ph->IsTagged()) fhPtConversionTagged ->Fill(ptcluster);
+          if(ptcluster > 0.5)fhEtaPhiConversion   ->Fill(etacluster,phicluster);
+          else               fhEtaPhi05Conversion ->Fill(etacluster,phicluster);
         }                      
         
         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
@@ -728,6 +1196,18 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
           fhEtaOtherDecay ->Fill(ptcluster,etacluster);
         }
       }
+      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
+      {
+        fhPtAntiNeutron  ->Fill(ptcluster);
+        fhPhiAntiNeutron ->Fill(ptcluster,phicluster);
+        fhEtaAntiNeutron ->Fill(ptcluster,etacluster);
+      }
+      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
+      {
+        fhPtAntiProton  ->Fill(ptcluster);
+        fhPhiAntiProton ->Fill(ptcluster,phicluster);
+        fhEtaAntiProton ->Fill(ptcluster,etacluster);
+      }      
            else{
              fhPtUnknown  ->Fill(ptcluster);
              fhPhiUnknown ->Fill(ptcluster,phicluster);
@@ -742,7 +1222,7 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
         
            }
            
-           
+           //....................................................................
            // Access MC information in stack if requested, check that it exists.
            Int_t label =ph->GetLabel();
            if(label < 0) {
index a099d94..44ea8c0 100755 (executable)
@@ -9,6 +9,8 @@
 // Class for the photon identification.
 // Clusters from calorimeters are identified as photons
 // and kept in the AOD. Few histograms produced.
+// Produces input for other analysis classes like AliAnaPi0, 
+// AliAnaParticleHadronCorrelation ... 
 //
 
 //-- Author: Gustavo Conesa (INFN-LNF)
@@ -38,113 +40,191 @@ class AliAnaPhoton : public AliAnaPartCorrBaseClass {
 
  public:
        
+  //---------------------------------------
+  // General analysis frame methods
+  //---------------------------------------
+
   TObjString * GetAnalysisCuts();
+  
   TList      * GetCreateOutputObjects();
 
-  void Init();
+  void         Init();
 
-  void MakeAnalysisFillAOD()  ;
-  
-  void MakeAnalysisFillHistograms() ; 
+  void         InitParameters();
+
+  void         MakeAnalysisFillAOD()  ;
+
+  void         MakeAnalysisFillHistograms() ; 
   
-  void Print(const Option_t * opt)const;
+  void         Print(const Option_t * opt)const;
   
-  TString GetCalorimeter()   const {return fCalorimeter ; }
-  void SetCalorimeter(TString &det)    {fCalorimeter = det ; }
-
-  Bool_t IsTrackMatchRejectionOn()   const {return fRejectTrackMatch ; }
-  void SwitchOnTrackMatchRejection()  {fRejectTrackMatch = kTRUE ; }
-  void SwitchOffTrackMatchRejection() {fRejectTrackMatch = kFALSE ; }  
+  //---------------------------------------
+  // Analysis parameters setters getters
+  //---------------------------------------
 
-  Bool_t IsCheckConversionOn()   const {return fCheckConversion ; }
-  void SwitchOnConversionChecker()  {fCheckConversion = kTRUE ; }
-  void SwitchOffConversionChecker() {fCheckConversion = kFALSE ; }  
-       
-  Bool_t AreConvertedPairsInAOD()   const {return fAddConvertedPairsToAOD ; }
-  void SwitchOnAdditionConvertedPairsToAOD()  {fAddConvertedPairsToAOD = kTRUE ; }
-  void SwitchOffAdditionConvertedPairsToAOD() {fAddConvertedPairsToAOD = kFALSE ; }  
-       
-  void InitParameters();
-  void SetMinDistanceToBadChannel(Float_t m1, Float_t m2, Float_t m3) {
+  TString GetCalorimeter()                  const {return fCalorimeter ; }
+  void    SetCalorimeter(TString  & det)          {fCalorimeter = det  ; }
+  
+  
+  // ** Cluster selection methods **
+  
+  void    SetMinDistanceToBadChannel(Float_t m1, Float_t m2, Float_t m3) {
     fMinDist = m1;
     fMinDist2 = m2;
     fMinDist3 = m3;
   }
+
+  void     SetTimeCut(Double_t min, Double_t max) {fTimeCutMin = min; fTimeCutMax = max ; }
+  Double_t GetTimeCutMin()                  const {return fTimeCutMin ; }
+  Double_t GetTimeCutMax()                  const {return fTimeCutMax ; }      
        
-  Float_t GetMassCut()    const {return fMassCut ; }
-  void SetMassCut(Float_t m)    {fMassCut = m ; }
-       
-  void SetTimeCut(Double_t min, Double_t max) {fTimeCutMin = min; fTimeCutMax = max;}
-  Double_t GetTimeCutMin() const {return fTimeCutMin;}
-  Double_t GetTimeCutMax() const {return fTimeCutMax;} 
-       
-  void SetNCellCut(Int_t n) {fNCellsCut = n;}
-  Double_t GetNCellCut() const {return fNCellsCut;}
+  void     SetNCellCut(Int_t n)                   {fNCellsCut = n    ; }
+  Double_t GetNCellCut()                    const {return fNCellsCut ; }
+  
+  Bool_t   IsTrackMatchRejectionOn()        const {return fRejectTrackMatch   ; }
+  void     SwitchOnTrackMatchRejection()          {fRejectTrackMatch = kTRUE  ; }
+  void     SwitchOffTrackMatchRejection()         {fRejectTrackMatch = kFALSE ; }  
+
+  // ** Conversion pair analysis **
+  
+  Float_t  GetMassCut()                     const {return fMassCut ; }
+  void     SetMassCut(Float_t m)                  {fMassCut = m ; }
+  
+  Bool_t   IsCheckConversionOn()            const {return fCheckConversion   ; }
+  void     SwitchOnConversionChecker()            {fCheckConversion = kTRUE  ; }
+  void     SwitchOffConversionChecker()           {fCheckConversion = kFALSE ; }  
        
+  Bool_t   AreConvertedPairsInAOD()         const {return fAddConvertedPairsToAOD   ; }
+  void     SwitchOnAdditionConvertedPairsToAOD()  {fAddConvertedPairsToAOD = kTRUE  ; }
+  void     SwitchOffAdditionConvertedPairsToAOD() {fAddConvertedPairsToAOD = kFALSE ; }  
        
+  void     SetConvDEtaCut(Float_t c)              {fConvDEtaCut = c    ; }
+  Float_t  GetConvDEtaCut()                 const {return fConvDEtaCut ; }
+  void     SetConvDPhiCut(Float_t c)              {fConvDPhiCut = c    ; }
+  Float_t  GetConvDPhiCut()                 const {return fConvDPhiCut ; }
+  
   private:
  
-  TString fCalorimeter ; // Calorimeter where the gamma is searched;
-  Float_t fMinDist ;     // Minimal distance to bad channel to accept cluster
-  Float_t fMinDist2;     // Cuts on Minimal distance to study acceptance evaluation
-  Float_t fMinDist3;     // One more cut on distance used for acceptance-efficiency study
-  Bool_t  fRejectTrackMatch ;      //If PID on, reject clusters which have an associated TPC track
-  Bool_t  fCheckConversion;        // Combine pairs of clusters with mass close to 0
-  Bool_t  fAddConvertedPairsToAOD; // Put Converted pairs in AOD
-  Float_t fMassCut;                // Mass cut for the conversion pairs selection
-  Double_t fTimeCutMin  ;    // Remove clusters/cells with time smaller than this value, in ns
-  Double_t fTimeCutMax  ;    // Remove clusters/cells with time larger than this value, in ns
-  Int_t fNCellsCut ;     // Accept for the analysis clusters with more than fNCellsCut cells
-       
+  TString  fCalorimeter ;                // Calorimeter where the gamma is searched;
+  Float_t  fMinDist ;                    // Minimal distance to bad channel to accept cluster
+  Float_t  fMinDist2;                    // Cuts on Minimal distance to study acceptance evaluation
+  Float_t  fMinDist3;                    // One more cut on distance used for acceptance-efficiency study
+  Bool_t   fRejectTrackMatch ;           // If PID on, reject clusters which have an associated TPC track
+  Bool_t   fCheckConversion;             // Combine pairs of clusters with mass close to 0
+  Bool_t   fAddConvertedPairsToAOD;      // Put Converted pairs in AOD
+  Float_t  fMassCut;                     // Mass cut for the conversion pairs selection
+  Double_t fTimeCutMin  ;                // Remove clusters/cells with time smaller than this value, in ns
+  Double_t fTimeCutMax  ;                // Remove clusters/cells with time larger than this value, in ns
+  Int_t    fNCellsCut ;                  // Accept for the analysis clusters with more than fNCellsCut cells
+       Float_t  fConvDEtaCut;                 // Select conversion pairs when deta of pair smaller than cut
+  Float_t  fConvDPhiCut;                 // Select conversion pairs when dphi of pair smaller than cut
+  
   //Histograms 
-  TH3D * fhVertex; //! vertex distribution
-  TH2F * fhNtraNclu; //! track multiplicity distribution vs cluster multiplicity
-  TH1F * fhPtPhoton   ; //! Number of identified photon vs transerse momentum 
-  TH2F * fhPhiPhoton  ; //! Azimuthal angle of identified  photon vs transerse momentum 
-  TH2F * fhEtaPhoton  ; //! Pseudorapidity of identified  photon vs transerse momentum 
-               
-  //MC
-  TH1F * fhDeltaE  ; //! MC-Reco E distribution      
-  TH1F * fhDeltaPt ; //! MC-Reco pT distribution
-  TH1F * fhRatioE  ; //! Reco/MC E distribution      
-  TH1F * fhRatioPt ; //! Reco/MC pT distribution
-  TH2F * fh2E  ; //! E distribution, Reco vs MC
-  TH2F * fh2Pt ; //! pT distribution, Reco vs MC
+  TH2F * fhNtraNclu;                     //! track multiplicity distribution vs cluster multiplicity
+  TH2F * fhNCellsPt;                     //! number of cells in cluster vs pt 
+  TH1F * fhPtPhoton   ;                  //! Number of identified photon vs transerse momentum 
+  TH2F * fhPhiPhoton  ;                  //! Azimuthal angle of identified  photon vs transerse momentum 
+  TH2F * fhEtaPhoton  ;                  //! Pseudorapidity of identified  photon vs transerse momentum 
+  TH2F * fhEtaPhiPhoton  ;               //! Pseudorapidity vs Phi of identified  photon for transerse momentum > 0.5
+  TH2F * fhEtaPhi05Photon  ;             //! Pseudorapidity vs Phi of identified  photon for transerse momentum < 0.5
   
-  TH1F * fhPtMCPhoton;   //! Number of identified gamma 
-  TH2F * fhPhiMCPhoton;  //! Phi of identified gamma
-  TH2F * fhEtaMCPhoton;  //! eta of identified gamma   
+  //Conversion pairs
+  TH1F * fhPtPhotonConv   ;              //! Number of identified photon vs transerse momentum 
+  TH2F * fhEtaPhiPhotonConv  ;           //! Pseudorapidity vs Phi of identified  photon for transerse momentum > 0.5, for converted
+  TH2F * fhEtaPhi05PhotonConv  ;         //! Pseudorapidity vs Phi of identified  photon for transerse momentum < 0.5, for converted
+  TH2F * fhConvDeltaEta;                 //! Small mass photons, correlation in eta
+  TH2F * fhConvDeltaPhi;                 //! Small mass photons, correlation in phi
+  TH2F * fhConvDeltaEtaPhi;              //! Small mass photons, correlation in phi and eta
+  TH2F * fhConvAsym;                     //! Small mass photons, correlation in energy asymmetry
+  TH2F * fhConvPt;                       //! Small mass photons, pT of pair
+
+  //Fill MC dependent histograms
+  TH1F * fhDeltaE  ;                     //! MC-Reco E distribution      
+  TH1F * fhDeltaPt ;                     //! MC-Reco pT distribution
+  TH1F * fhRatioE  ;                     //! Reco/MC E distribution      
+  TH1F * fhRatioPt ;                     //! Reco/MC pT distribution
+  TH2F * fh2E  ;                         //! E distribution, Reco vs MC
+  TH2F * fh2Pt ;                         //! pT distribution, Reco vs MC
+  
+  //Origin of this cluster is ...
+  TH1F * fhPtMCPhoton;                   //! Number of identified gamma 
+  TH2F * fhPhiMCPhoton;                  //! Phi of identified gamma
+  TH2F * fhEtaMCPhoton;                  //! eta of identified gamma   
        
-  TH1F * fhPtPrompt;   //! Number of identified prompt gamma 
-  TH2F * fhPhiPrompt;  //! Phi of identified  prompt gamma
-  TH2F * fhEtaPrompt;  //! eta of identified  prompt gamma
-
-  TH1F * fhPtFragmentation;   //! Number of identified fragmentation gamma 
-  TH2F * fhPhiFragmentation;  //! Phi of identified  fragmentation gamma
-  TH2F * fhEtaFragmentation;  //! eta of identified  fragmentation gamma
-
-  TH1F * fhPtISR;   //! Number of identified initial state radiation gamma 
-  TH2F * fhPhiISR;  //! Phi of identified initial state radiation gamma
-  TH2F * fhEtaISR;  //! eta of identified initial state radiation gamma
-
-  TH1F * fhPtPi0Decay;   //! Number of identified Pi0Decay gamma 
-  TH2F * fhPhiPi0Decay;  //! Phi of identified  Pi0Decay gamma
-  TH2F * fhEtaPi0Decay;  //! eta of identified  Pi0Decay gamma
-
-  TH1F * fhPtOtherDecay;   //! Number of identified OtherDecay gamma 
-  TH2F * fhPhiOtherDecay;  //! Phi of identified  OtherDecay gamma
-  TH2F * fhEtaOtherDecay;  //! eta of identified  OtherDecay gamma
-
-  TH1F * fhPtConversion;   //! Number of identified Conversion gamma 
-  TH2F * fhPhiConversion;  //! Phi of identified  Conversion gamma
-  TH2F * fhEtaConversion;  //! eta of identified  Conversion gamma
-
-  TH1F * fhPtUnknown;   //! Number of identified Unknown gamma 
-  TH2F * fhPhiUnknown;  //! Phi of identified  Unknown gamma
-  TH2F * fhEtaUnknown;  //! eta of identified  Unknown gamma
-
-   ClassDef(AliAnaPhoton,8)
+  TH1F * fhPtPrompt;                     //! Number of identified prompt gamma 
+  TH2F * fhPhiPrompt;                    //! Phi of identified  prompt gamma
+  TH2F * fhEtaPrompt;                    //! eta of identified  prompt gamma
+
+  TH1F * fhPtFragmentation;              //! Number of identified fragmentation gamma 
+  TH2F * fhPhiFragmentation;             //! Phi of identified  fragmentation gamma
+  TH2F * fhEtaFragmentation;             //! eta of identified  fragmentation gamma
+
+  TH1F * fhPtISR;                        //! Number of identified initial state radiation gamma 
+  TH2F * fhPhiISR;                       //! Phi of identified initial state radiation gamma
+  TH2F * fhEtaISR;                       //! eta of identified initial state radiation gamma
+
+  TH1F * fhPtPi0Decay;                   //! Number of identified Pi0Decay gamma 
+  TH2F * fhPhiPi0Decay;                  //! Phi of identified  Pi0Decay gamma
+  TH2F * fhEtaPi0Decay;                  //! eta of identified  Pi0Decay gamma
+
+  TH1F * fhPtOtherDecay;                 //! Number of identified OtherDecay gamma 
+  TH2F * fhPhiOtherDecay;                //! Phi of identified  OtherDecay gamma
+  TH2F * fhEtaOtherDecay;                //! eta of identified  OtherDecay gamma
+
+  TH1F * fhPtConversion;                 //! Number of identified Conversion gamma 
+  TH2F * fhPhiConversion;                //! Phi of identified  Conversion gamma
+  TH2F * fhEtaConversion;                //! eta of identified  Conversion gamma
+  TH1F * fhPtConversionTagged;           //! Number of identified Conversion gamma, tagged as conversion 
+
+  TH2F * fhEtaPhiConversion  ;           //! Pseudorapidity vs Phi of identified  photon for transerse momentum > 0.5
+  TH2F * fhEtaPhi05Conversion  ;         //! Pseudorapidity vs Phi of identified  photon for transerse momentum < 0.5
+
+  TH1F * fhPtAntiNeutron;                //! Origin are AntiNeutrons or AntiProtons
+  TH2F * fhPhiAntiNeutron;               //! Origin are AntiNeutrons or AntiProtons
+  TH2F * fhEtaAntiNeutron;               //! Origin are AntiNeutrons or AntiProtons
+    
+  TH1F * fhPtAntiProton;                 //! Origin are AntiNeutrons or AntiProtons
+  TH2F * fhPhiAntiProton;                //! Origin are AntiNeutrons or AntiProtons
+  TH2F * fhEtaAntiProton;                //! Origin are AntiNeutrons or AntiProtons
+  
+  TH1F * fhPtUnknown;                    //! Number of identified Unknown gamma 
+  TH2F * fhPhiUnknown;                   //! Phi of identified  Unknown gamma
+  TH2F * fhEtaUnknown;                   //! eta of identified  Unknown gamma
+
+  //Conversion pairs analysis histograms
+  TH2F * fhConvDeltaEtaMCConversion;     //! Small mass cluster pairs, correlation in eta, origin of both clusters is conversion
+  TH2F * fhConvDeltaPhiMCConversion;     //! Small mass cluster pairs, correlation in phi, origin of both clusters is conversion
+  TH2F * fhConvDeltaEtaPhiMCConversion;  //! Small mass cluster pairs, correlation in eta-phi, origin of both clusters is conversion
+  TH2F * fhConvAsymMCConversion;         //! Small mass cluster pairs, correlation in energy asymmetry, origin of both clusters is conversion
+  TH2F * fhConvPtMCConversion;           //! Small mass cluster pairs, pt of pair, origin of both clusters is conversion
+  TH2F * fhConvDispersionMCConversion;   //! Small mass cluster pairs, dispersion of cluster 1 vs cluster 2
+  TH2F * fhConvM02MCConversion;          //! Small mass cluster pairs, m02 of cluster 1 vs cluster 2
+
+  TH2F * fhConvDeltaEtaMCAntiNeutron;    //! Small mass cluster pairs, correlation in eta, origin of both clusters is anti neutron
+  TH2F * fhConvDeltaPhiMCAntiNeutron;    //! Small mass cluster pairs, correlation in phi, origin of both clusters is anti neutron
+  TH2F * fhConvDeltaEtaPhiMCAntiNeutron; //! Small mass cluster pairs, correlation in eta-phi, origin of both clusters is anti neutron
+  TH2F * fhConvAsymMCAntiNeutron;        //! Small mass cluster pairs, correlation in energy asymmetry, origin of both clusters is anti neutron
+  TH2F * fhConvPtMCAntiNeutron;          //! Small mass cluster pairs, pt of pair, origin of both clusters is anti neutron
+  TH2F * fhConvDispersionMCAntiNeutron;  //! Small mass cluster pairs, dispersion of cluster 1 vs cluster 2, origin of both clusters is anti neutron
+  TH2F * fhConvM02MCAntiNeutron;         //! Small mass cluster pairs, m02 of cluster 1 vs cluster 2, origin of both clusters is anti neutron
+
+  TH2F * fhConvDeltaEtaMCAntiProton;     //! Small mass cluster pairs, correlation in eta, origin of both clusters is anti proton
+  TH2F * fhConvDeltaPhiMCAntiProton;     //! Small mass cluster pairs, correlation in phi, origin of both clusters is anti proton
+  TH2F * fhConvDeltaEtaPhiMCAntiProton;  //! Small mass cluster pairs, correlation in eta-phi, origin of both clusters is anti proton
+  TH2F * fhConvAsymMCAntiProton;         //! Small mass cluster pairs, correlation in energy asymmetry, origin of both clusters is anti proton
+  TH2F * fhConvPtMCAntiProton;           //! Small mass cluster pairs, pt of pairs, origin of both clusters is anti proton
+  TH2F * fhConvDispersionMCAntiProton;   //! Small mass cluster pairs, dispersion of cluster 1 vs cluster 2, origin of both clusters is anti proton
+  TH2F * fhConvM02MCAntiProton;          //! Small mass cluster pairs, m02 of cluster 1 vs cluster 2, origin of both clusters is anti proton
+
+  TH2F * fhConvDeltaEtaMCString;         //! Small mass cluster pairs, correlation in eta, origin of both clusters is string
+  TH2F * fhConvDeltaPhiMCString;         //! Small mass cluster pairs, correlation in phi, origin of both clusters is string
+  TH2F * fhConvDeltaEtaPhiMCString;      //! Small mass cluster pairs, correlation in eta-phi, origin of both clusters is string
+  TH2F * fhConvAsymMCString;             //! Small mass cluster pairs, correlation in energy asymmetry, origin of both clusters is string
+  TH2F * fhConvPtMCString;               //! Small mass cluster pairs, pt of pairs, origin of both clusters is string
+  TH2F * fhConvDispersionMCString;       //! Small mass cluster pairs, dispersion of cluster 1 vs cluster 2, origin of both clusters is string
+  TH2F * fhConvM02MCString;              //! Small mass cluster pairs, m02 of cluster 1 vs cluster 2, origin of both clusters is string
+
+   ClassDef(AliAnaPhoton,9)
 
 } ;
  
index 88313e9..336d8f9 100755 (executable)
@@ -16,7 +16,9 @@
 
 //_________________________________________________________________________
 // Class to collect two-photon invariant mass distributions for
-// extractin raw pi0 yield.
+// extracting raw pi0 yield.
+// Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles), 
+// it will do nothing if executed alone
 //
 //-- Author: Dmitri Peressounko (RRC "KI") 
 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
@@ -33,6 +35,7 @@
 #include "TROOT.h"
 #include "TClonesArray.h"
 #include "TObjString.h"
+#include "TDatabasePDG.h"
 
 //---- AliRoot system ----
 #include "AliAnaPi0.h"
@@ -47,7 +50,7 @@
 #include "AliAODEvent.h"
 #include "AliNeutralMesonSelection.h"
 #include "AliMixedEvent.h"
-
+#include "AliAODMCParticle.h"
 
 ClassImp(AliAnaPi0)
 
@@ -55,15 +58,22 @@ ClassImp(AliAnaPi0)
 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
 fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
 fNmaxMixEv(0), fCalorimeter(""),
-fNModules(12), fUseAngleCut(kFALSE), fEventsList(0x0), fMultiCutAna(kFALSE), 
-fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0), fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
-fhReMod(0x0),fhReDiffMod(0x0),
+fNModules(12), fUseAngleCut(kFALSE), fUseAngleEDepCut(kFALSE),fAngleCut(0), fAngleMaxCut(7.),fEventsList(0x0), fMultiCutAna(kFALSE), fMultiCutAnaSim(kFALSE),
+fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0),  fMakeInvPtPlots(kFALSE), fSameSM(kFALSE),
+fUseTrackMultBins(kFALSE),fUsePhotonMultBins(kFALSE),fUseAverClusterEBins(kFALSE),fUseAverCellEBins(kFALSE), fFillBadDistHisto(kFALSE),
+fhAverTotECluster(0),fhAverTotECell(0),
+fhReMod(0x0),   fhReDiffMod(0x0), fhMiMod(0x0),    fhMiDiffMod(0x0),
+fhReConv(0x0),   fhMiConv(0x0),   fhReConv2(0x0),  fhMiConv2(0x0),
 fhRe1(0x0),      fhMi1(0x0),      fhRe2(0x0),      fhMi2(0x0),      fhRe3(0x0),      fhMi3(0x0),
 fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
-fhRePtNCellAsymCuts(0x0), fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),  
-fhEvents(0x0), fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0),
-fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0),
-fhPrimOpeningAngle(0x0),fhPrimCosOpeningAngle(0x0)
+fhRePtNCellAsymCuts(0x0), fhRePtNCellAsymCutsSM0(0x0), fhRePtNCellAsymCutsSM1(0x0), fhRePtNCellAsymCutsSM2(0x0), fhRePtNCellAsymCutsSM3(0x0), fhMiPtNCellAsymCuts(0x0),
+fhRePIDBits(0x0),fhRePtMult(0x0), fhRePtAsym(0x0), fhRePtAsymPi0(0x0),fhRePtAsymEta(0x0),  
+fhEvents(0x0), fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0), fhMixedOpeningAngle(0x0),fhMixedCosOpeningAngle(0x0),
+fhPrimPi0Pt(0x0), fhPrimPi0AccPt(0x0), fhPrimPi0Y(0x0), fhPrimPi0AccY(0x0), fhPrimPi0Phi(0x0), fhPrimPi0AccPhi(0x0),
+fhPrimPi0OpeningAngle(0x0),fhPrimPi0CosOpeningAngle(0x0),
+fhPrimEtaPt(0x0), fhPrimEtaAccPt(0x0), fhPrimEtaY(0x0), fhPrimEtaAccY(0x0), fhPrimEtaPhi(0x0), fhPrimEtaAccPhi(0x0),
+fhMCOrgMass(),fhMCOrgAsym(),  fhMCOrgDeltaEta(),fhMCOrgDeltaPhi(),
+fhMCPi0MassPtRec(), fhMCPi0MassPtTrue(), fhMCPi0PtTruePtRec(), fhMCEtaMassPtRec(), fhMCEtaMassPtTrue(), fhMCEtaPtTruePtRec()
 {
 //Default Ctor
  InitParameters();
@@ -105,7 +115,10 @@ void AliAnaPi0::InitParameters()
  
   fCalorimeter  = "PHOS";
   fUseAngleCut = kFALSE;
-  
+  fUseAngleEDepCut = kFALSE;
+  fAngleCut    = 0.; 
+  fAngleMaxCut = TMath::Pi(); 
+
   fMultiCutAna = kFALSE;
   
   fNPtCuts = 3;
@@ -144,7 +157,10 @@ TObjString * AliAnaPi0::GetAnalysisCuts()
   parList+=onePar ;
   snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
   parList+=onePar ;
-  snprintf(onePar,buffersize,"Pair in same Module: %d; Make InvPt plots %d \n",fSameSM, fMakeInvPtPlots) ;
+  snprintf(onePar,buffersize,"Pair in same Module: %d ; TrackMult as centrality: %d; PhotonMult as centrality: %d; cluster E as centrality: %d; cell as centrality: %d; Fill InvPt histos %d\n",
+           fSameSM, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
+  parList+=onePar ;
+  snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f,\n",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
   parList+=onePar ;
   snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
   for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
@@ -194,22 +210,30 @@ TList * AliAnaPi0::GetCreateOutputObjects()
   outputContainer->SetName(GetName()); 
        
   fhReMod     = new TH2D*[fNModules] ;
-  fhReDiffMod = new TH2D*[fNModules+1] ;
+  fhReDiffMod = new TH2D*[fNModules+3] ;
+  
+  fhMiMod     = new TH2D*[fNModules] ;
+  fhMiDiffMod = new TH2D*[fNModules+3] ;
   
   fhRe1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-  fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-  fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
   fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-  fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-  fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-  if(fMakeInvPtPlots ) {
+  if(fFillBadDistHisto){
+    fhRe2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+    fhRe3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+    fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+    fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+  }
+  if(fMakeInvPtPlots) {
     fhReInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-    fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-    fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
     fhMiInvPt1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-    fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-    fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+    if(fFillBadDistHisto){
+      fhReInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+      fhReInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+      fhMiInvPt2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+      fhMiInvPt3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+    }
   } 
+  
   const Int_t buffersize = 255;
   char key[buffersize] ;
   char title[buffersize] ;
@@ -233,130 +257,163 @@ TList * AliAnaPi0::GetCreateOutputObjects()
   Int_t ntrmbins  = GetHistoTrackMultiplicityBins();
   Int_t ntrmmax   = GetHistoTrackMultiplicityMax();
   Int_t ntrmmin   = GetHistoTrackMultiplicityMin(); 
-
+  
+  fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
+  fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
+  outputContainer->Add(fhAverTotECluster) ;
+  
+  fhAverTotECell    = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
+  fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
+  outputContainer->Add(fhAverTotECell) ;
+  
+  fhReConv = new TH2D("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+  fhReConv->SetXTitle("p_{T} (GeV/c)");
+  fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+  outputContainer->Add(fhReConv) ;
+  
+  fhMiConv = new TH2D("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+  fhMiConv->SetXTitle("p_{T} (GeV/c)");
+  fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+  outputContainer->Add(fhMiConv) ;
+  
+  fhReConv2 = new TH2D("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+  fhReConv2->SetXTitle("p_{T} (GeV/c)");
+  fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+  outputContainer->Add(fhReConv2) ;
+  
+  fhMiConv2 = new TH2D("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+  fhMiConv2->SetXTitle("p_{T} (GeV/c)");
+  fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+  outputContainer->Add(fhMiConv2) ;
+  
   for(Int_t ic=0; ic<fNCentrBin; ic++){
-    for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
-      for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
-        Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
-        //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
-        //Distance to bad module 1
-        snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
-        snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
-                 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-        fhRe1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-        fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
-        fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-        //printf("name: %s\n ",fhRe1[index]->GetName());
-        outputContainer->Add(fhRe1[index]) ;
-        
-        //Distance to bad module 2
-        snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
-        snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
-                 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-        fhRe2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-        fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
-        fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-        outputContainer->Add(fhRe2[index]) ;
-        
-        //Distance to bad module 3
-        snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
-        snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
-                 ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-        fhRe3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-        fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
-        fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-        outputContainer->Add(fhRe3[index]) ;
-        
-        if(fMakeInvPtPlots ) {
-          //Inverse pT 
+      for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
+        for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
+          Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
+          //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
           //Distance to bad module 1
-          snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
+          snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-          fhReInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-          fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
-          fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-          outputContainer->Add(fhReInvPt1[index]) ;
-          
-          //Distance to bad module 2
-          snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
-          snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
-                   ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-          fhReInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-          fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
-          fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-          outputContainer->Add(fhReInvPt2[index]) ;
-          
-          //Distance to bad module 3
-          snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
-          snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
-                   ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-          fhReInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-          fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
-          fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-          outputContainer->Add(fhReInvPt3[index]) ;
-        }
-        
-        if(fDoOwnMix){
-          //Distance to bad module 1
-          snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
-          snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
-                   ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-          fhMi1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-          fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
-          fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-          outputContainer->Add(fhMi1[index]) ;
+          fhRe1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+          fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
+          fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+          //printf("name: %s\n ",fhRe1[index]->GetName());
+          outputContainer->Add(fhRe1[index]) ;
           
-          //Distance to bad module 2
-          snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
-          snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
-                   ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-          fhMi2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-          fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
-          fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-          outputContainer->Add(fhMi2[index]) ;
-          
-          //Distance to bad module 3
-          snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
-          snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
-                   ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-          fhMi3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-          fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
-          fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-          outputContainer->Add(fhMi3[index]) ;
-          if(fMakeInvPtPlots ) {
-            //Inverse pT
-            //Distance to bad module 1
-            snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
-            snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
-                     ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-            fhMiInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-            fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
-            fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-            outputContainer->Add(fhMiInvPt1[index]) ;
-            
+          if(fFillBadDistHisto){
             //Distance to bad module 2
-            snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
-            snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
+            snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
+            snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-            fhMiInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-            fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
-            fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-            outputContainer->Add(fhMiInvPt2[index]) ;
+            fhRe2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+            fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
+            fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+            outputContainer->Add(fhRe2[index]) ;
             
             //Distance to bad module 3
-            snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
-            snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
+            snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
+            snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
-            fhMiInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-            fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
-            fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
-            outputContainer->Add(fhMiInvPt3[index]) ;
+            fhRe3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+            fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
+            fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+            outputContainer->Add(fhRe3[index]) ;
           }
-        } 
+          
+          //Inverse pT 
+          if(fMakeInvPtPlots){
+            //Distance to bad module 1
+            snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
+            snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
+                     ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
+            fhReInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+            fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
+            fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+            outputContainer->Add(fhReInvPt1[index]) ;
+            
+            if(fFillBadDistHisto){
+              //Distance to bad module 2
+              snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
+              snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
+                       ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
+              fhReInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+              fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
+              fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+              outputContainer->Add(fhReInvPt2[index]) ;
+              
+              //Distance to bad module 3
+              snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
+              snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
+                       ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
+              fhReInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+              fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
+              fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+              outputContainer->Add(fhReInvPt3[index]) ;
+            }
+          }
+          if(fDoOwnMix){
+            //Distance to bad module 1
+            snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
+            snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
+                     ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
+            fhMi1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+            fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
+            fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+            outputContainer->Add(fhMi1[index]) ;
+            if(fFillBadDistHisto){
+              //Distance to bad module 2
+              snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
+              snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
+                       ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
+              fhMi2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+              fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
+              fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+              outputContainer->Add(fhMi2[index]) ;
+              
+              //Distance to bad module 3
+              snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
+              snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
+                       ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
+              fhMi3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+              fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
+              fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+              outputContainer->Add(fhMi3[index]) ;
+            }
+            //Inverse pT
+            if(fMakeInvPtPlots){
+              //Distance to bad module 1
+              snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
+              snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
+                       ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
+              fhMiInvPt1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+              fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
+              fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+              outputContainer->Add(fhMiInvPt1[index]) ;
+              if(fFillBadDistHisto){
+                //Distance to bad module 2
+                snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
+                snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
+                         ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
+                fhMiInvPt2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+                fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
+                fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+                outputContainer->Add(fhMiInvPt2[index]) ;
+                
+                //Distance to bad module 3
+                snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
+                snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
+                         ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
+                fhMiInvPt3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+                fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
+                fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+                outputContainer->Add(fhMiInvPt3[index]) ;
+              }
+            }
+          } 
+        }
       }
     }
-  }
   
   fhRePtAsym = new TH2D("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
   fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
@@ -385,7 +442,12 @@ TList * AliAnaPi0::GetCreateOutputObjects()
       outputContainer->Add(fhRePIDBits[ipid]) ;
     }// pid bit loop
     
-    fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+    fhRePtNCellAsymCuts    = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+    fhRePtNCellAsymCutsSM0 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+    fhRePtNCellAsymCutsSM1 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+    fhRePtNCellAsymCutsSM2 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+    fhRePtNCellAsymCutsSM3 = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+    fhMiPtNCellAsymCuts    = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
     for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
       for(Int_t icell=0; icell<fNCellNCuts; icell++){
         for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
@@ -397,6 +459,42 @@ TList * AliAnaPi0::GetCreateOutputObjects()
           fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
           fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
           outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
+                    
+          snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d_SM0",ipt,icell,iasym) ;
+          snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 0 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
+          fhRePtNCellAsymCutsSM0[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+          fhRePtNCellAsymCutsSM0[index]->SetXTitle("p_{T} (GeV/c)");
+          fhRePtNCellAsymCutsSM0[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+          outputContainer->Add(fhRePtNCellAsymCutsSM0[index]) ;
+          
+          snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d_SM1",ipt,icell,iasym) ;
+          snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 1 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
+          fhRePtNCellAsymCutsSM1[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+          fhRePtNCellAsymCutsSM1[index]->SetXTitle("p_{T} (GeV/c)");
+          fhRePtNCellAsymCutsSM1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+          outputContainer->Add(fhRePtNCellAsymCutsSM1[index]) ;
+          
+          snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d_SM2",ipt,icell,iasym) ;
+          snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 2 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
+          fhRePtNCellAsymCutsSM2[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+          fhRePtNCellAsymCutsSM2[index]->SetXTitle("p_{T} (GeV/c)");
+          fhRePtNCellAsymCutsSM2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+          outputContainer->Add(fhRePtNCellAsymCutsSM2[index]) ;
+          
+          snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d_SM3",ipt,icell,iasym) ;
+          snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM 3 ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
+          fhRePtNCellAsymCutsSM3[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+          fhRePtNCellAsymCutsSM3[index]->SetXTitle("p_{T} (GeV/c)");
+          fhRePtNCellAsymCutsSM3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+          outputContainer->Add(fhRePtNCellAsymCutsSM3[index]) ;
+          
+          snprintf(key,   buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
+          snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
+          fhMiPtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+          fhMiPtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
+          fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+          outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
+          
         }
       }
     }
@@ -418,59 +516,202 @@ TList * AliAnaPi0::GetCreateOutputObjects()
   outputContainer->Add(fhEvents) ;
        
   fhRealOpeningAngle  = new TH2D
-  ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5); 
+  ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi()); 
   fhRealOpeningAngle->SetYTitle("#theta(rad)");
   fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
   outputContainer->Add(fhRealOpeningAngle) ;
   
   fhRealCosOpeningAngle  = new TH2D
-  ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1); 
+  ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1); 
   fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
   fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
   outputContainer->Add(fhRealCosOpeningAngle) ;
        
+  if(fDoOwnMix){
+    
+    fhMixedOpeningAngle  = new TH2D
+    ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi()); 
+    fhMixedOpeningAngle->SetYTitle("#theta(rad)");
+    fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
+    outputContainer->Add(fhMixedOpeningAngle) ;
+    
+    fhMixedCosOpeningAngle  = new TH2D
+    ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1); 
+    fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
+    fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
+    outputContainer->Add(fhMixedCosOpeningAngle) ;
+    
+  }
+  
   //Histograms filled only if MC data is requested     
   if(IsDataMC()){
-
-    fhPrimPt     = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
-    fhPrimAccPt  = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
-    outputContainer->Add(fhPrimPt) ;
-    outputContainer->Add(fhPrimAccPt) ;
+    //Pi0
+    fhPrimPi0Pt     = new TH1D("hPrimPi0Pt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
+    fhPrimPi0AccPt  = new TH1D("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
+    outputContainer->Add(fhPrimPi0Pt) ;
+    outputContainer->Add(fhPrimPi0AccPt) ;
+    
+    fhPrimPi0Y      = new TH1D("hPrimPi0Rapidity","Rapidity of primary pi0",netabins,etamin,etamax) ; 
+    outputContainer->Add(fhPrimPi0Y) ;
+    
+    fhPrimPi0AccY   = new TH1D("hPrimPi0AccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ; 
+    outputContainer->Add(fhPrimPi0AccY) ;
     
-    fhPrimY      = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ; 
-    outputContainer->Add(fhPrimY) ;
+    fhPrimPi0Phi    = new TH1D("hPrimPi0Phi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; 
+    outputContainer->Add(fhPrimPi0Phi) ;
     
-    fhPrimAccY   = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ; 
-    outputContainer->Add(fhPrimAccY) ;
+    fhPrimPi0AccPhi = new TH1D("hPrimPi0AccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; 
+    outputContainer->Add(fhPrimPi0AccPhi) ;
     
-    fhPrimPhi    = new TH1D("hPrimaryPhi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; 
-    outputContainer->Add(fhPrimPhi) ;
+    //Eta
+    fhPrimEtaPt     = new TH1D("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
+    fhPrimEtaAccPt  = new TH1D("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
+    outputContainer->Add(fhPrimEtaPt) ;
+    outputContainer->Add(fhPrimEtaAccPt) ;
     
-    fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; 
-    outputContainer->Add(fhPrimAccPhi) ;
+    fhPrimEtaY      = new TH1D("hPrimEtaRapidity","Rapidity of primary eta",netabins,etamin,etamax) ; 
+    outputContainer->Add(fhPrimEtaY) ;
     
+    fhPrimEtaAccY   = new TH1D("hPrimEtaAccRapidity","Rapidity of primary eta",netabins,etamin,etamax) ; 
+    outputContainer->Add(fhPrimEtaAccY) ;
     
-    fhPrimOpeningAngle  = new TH2D
-    ("hPrimOpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5); 
-    fhPrimOpeningAngle->SetYTitle("#theta(rad)");
-    fhPrimOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
-    outputContainer->Add(fhPrimOpeningAngle) ;
+    fhPrimEtaPhi    = new TH1D("hPrimEtaPhi","Azimithal of primary eta",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; 
+    outputContainer->Add(fhPrimEtaPhi) ;
+    
+    fhPrimEtaAccPhi = new TH1D("hPrimEtaAccPhi","Azimithal of primary eta with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ; 
+    outputContainer->Add(fhPrimEtaAccPhi) ;
+        
     
-    fhPrimCosOpeningAngle  = new TH2D
-    ("hPrimCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1); 
-    fhPrimCosOpeningAngle->SetYTitle("cos (#theta) ");
-    fhPrimCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
-    outputContainer->Add(fhPrimCosOpeningAngle) ;
+    fhPrimPi0OpeningAngle  = new TH2D
+    ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5); 
+    fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
+    fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
+    outputContainer->Add(fhPrimPi0OpeningAngle) ;
+    
+    fhPrimPi0CosOpeningAngle  = new TH2D
+    ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1); 
+    fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
+    fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
+    outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
+    
+    for(Int_t i = 0; i<13; i++){
+      fhMCOrgMass[i] = new TH2D(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+      fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
+      fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+      outputContainer->Add(fhMCOrgMass[i]) ;
+      
+      fhMCOrgAsym[i]= new TH2D(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
+      fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
+      fhMCOrgAsym[i]->SetYTitle("A");
+      outputContainer->Add(fhMCOrgAsym[i]) ;
+      
+      fhMCOrgDeltaEta[i] = new TH2D(Form("hMCOrgDeltaEta_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
+      fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
+      fhMCOrgDeltaEta[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+      outputContainer->Add(fhMCOrgDeltaEta[i]) ;
+      
+      fhMCOrgDeltaPhi[i]= new TH2D(Form("hMCOrgDeltaPhi_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
+      fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
+      fhMCOrgDeltaPhi[i]->SetYTitle("A");
+      outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
+      
+    }
     
+    if(fMultiCutAnaSim){
+      fhMCPi0MassPtTrue  = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      fhMCPi0MassPtRec   = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      fhMCEtaMassPtRec   = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      fhMCEtaMassPtTrue  = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
+        for(Int_t icell=0; icell<fNCellNCuts; icell++){
+          for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
+            Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
+            
+            fhMCPi0MassPtRec[index] = new TH2D(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
+                                               Form("Reconstructed Mass vs reconstructed p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+                                               nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+            fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
+            fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+            outputContainer->Add(fhMCPi0MassPtRec[index]) ;    
+            
+            fhMCPi0MassPtTrue[index] = new TH2D(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
+                                                Form("Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+                                                nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+            fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
+            fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+            outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
+            
+            fhMCPi0PtTruePtRec[index] = new TH2D(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
+                                                 Form("Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+                                                 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
+            fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
+            fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
+            outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
+            
+            fhMCEtaMassPtRec[index] = new TH2D(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
+                                               Form("Reconstructed Mass vs reconstructed p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+                                               nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+            fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
+            fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+            outputContainer->Add(fhMCEtaMassPtRec[index]) ;
+            
+            fhMCEtaMassPtTrue[index] = new TH2D(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
+                                                Form("Reconstructed Mass vs generated p_T of true #eta cluster pairs for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+                                                nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+            fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
+            fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+            outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
+            
+            fhMCEtaPtTruePtRec[index] = new TH2D(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
+                                                 Form("Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2} for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
+                                                 nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
+            fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
+            fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
+            outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
+          }
+        }
+      }  
+    }//multi cut ana
+    else {
+      fhMCPi0MassPtTrue  = new TH2D*[1];
+      fhMCPi0PtTruePtRec = new TH2D*[1];
+      fhMCEtaMassPtTrue  = new TH2D*[1];
+      fhMCEtaPtTruePtRec = new TH2D*[1];
+      
+      fhMCPi0MassPtTrue[0] = new TH2D("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+      fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
+      fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+      outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
+      
+      fhMCPi0PtTruePtRec[0]= new TH2D("hMCPi0PtTruePtRec","Generated vs reconstructed p_T of true #pi^{0} cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
+      fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
+      fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
+      outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
+      
+      fhMCEtaMassPtTrue[0] = new TH2D("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+      fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
+      fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+      outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
+      
+      fhMCEtaPtTruePtRec[0]= new TH2D("hMCEtaPtTruePtRec","Generated vs reconstructed p_T of true #eta cluster pairs, 0.01 < rec. mass < 0.17 MeV/c^{2}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
+      fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
+      fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
+      outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
+    }
   }
   
-  TString * pairname = new TString[fNModules];
+  TString * pairname = new TString[fNModules+3];
   if(fCalorimeter=="EMCAL"){
     pairname[0]="A side (0-2)"; 
     pairname[1]="C side (1-3)";
     pairname[2]="Sector 0 (0-1)"; 
     pairname[3]="Sector 1 (2-3)";
-    for(Int_t i = 4 ; i < fNModules ; i++) pairname[i]="";}
+    pairname[4]="Cluster in different SM";
+    pairname[5]="SM 0 and SM3"; 
+    pairname[6]="SM 1 and SM2";    
+    for(Int_t i = 7 ; i < fNModules ; i++) pairname[i]="";}
   if(fCalorimeter=="PHOS") {
     pairname[0]="(0-1)"; 
     pairname[1]="(0-2)";
@@ -492,15 +733,45 @@ TList * AliAnaPi0::GetCreateOutputObjects()
     fhReDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
     fhReDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
     outputContainer->Add(fhReDiffMod[imod]) ;
+    
+    if(fDoOwnMix){ 
+      snprintf(key, buffersize,"hMiMod_%d",imod) ;
+      snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
+      fhMiMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+      fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
+      fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+      outputContainer->Add(fhMiMod[imod]) ;
+      
+      snprintf(key, buffersize,"hMiDiffMod_%d",imod) ;
+      snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
+      fhMiDiffMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+      fhMiDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
+      fhMiDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+      outputContainer->Add(fhMiDiffMod[imod]) ;
+    }
+    
   }
   
-  delete [] pairname;
+  for (Int_t imod=4; imod<7; imod++) {
+    
+    snprintf(key, buffersize,"hReDiffMod_%d",imod) ;
+    snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
+    fhReDiffMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+    fhReDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
+    fhReDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+    outputContainer->Add(fhReDiffMod[imod]) ;
   
-  snprintf(key, buffersize,"hReDiffMod_%d",fNModules) ;
-  snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for all Modules Combination") ;
-  fhReDiffMod[fNModules]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
-  outputContainer->Add(fhReDiffMod[fNModules]) ;
+    if(fDoOwnMix){ 
+      snprintf(key, buffersize,"hMiDiffMod_%d",imod) ;
+      snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
+      fhMiDiffMod[imod]  = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+      fhMiDiffMod[imod]->SetXTitle("p_{T} (GeV/c)");
+      fhMiDiffMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
+      outputContainer->Add(fhMiDiffMod[imod]) ;
+    }
+  }
   
+  delete [] pairname;
   
 //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
 //  
@@ -526,7 +797,7 @@ void AliAnaPi0::Print(const Option_t * /*opt*/) const
   printf("Cuts: \n") ;
   printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ;
   printf("Number of modules:             %d \n",fNModules) ;
-  printf("Select pairs with their angle: %d \n",fUseAngleCut) ;
+  printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
   printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
   printf("\tasymmetry < ");
   for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
@@ -556,24 +827,34 @@ void AliAnaPi0::Print(const Option_t * /*opt*/) const
 void AliAnaPi0::FillAcceptanceHistograms(){
   //Fill acceptance histograms if MC data is available
   
-  if(IsDataMC() && GetReader()->ReadStack()){  
+  if(GetReader()->ReadStack()){        
     AliStack * stack = GetMCStack();
-    if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
+    if(stack){
       for(Int_t i=0 ; i<stack->GetNprimary(); i++){
         TParticle * prim = stack->Particle(i) ;
-        if(prim->GetPdgCode() == 111){
+        Int_t pdg = prim->GetPdgCode();
+        if( pdg == 111 || pdg == 221){
           Double_t pi0Pt = prim->Pt() ;
           //printf("pi0, pt %2.2f\n",pi0Pt);
           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception         
           Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
           Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
-          if(TMath::Abs(pi0Y) < 0.5){
-            fhPrimPt->Fill(pi0Pt) ;
+          if(pdg == 111){
+            if(TMath::Abs(pi0Y) < 0.5){
+              fhPrimPi0Pt->Fill(pi0Pt) ;
+            }
+            fhPrimPi0Y  ->Fill(pi0Y) ;
+            fhPrimPi0Phi->Fill(phi) ;
+          }
+          else if(pdg == 221){
+            if(TMath::Abs(pi0Y) < 0.5){
+              fhPrimEtaPt->Fill(pi0Pt) ;
+            }
+            fhPrimEtaY  ->Fill(pi0Y) ;
+            fhPrimEtaPhi->Fill(phi) ;
           }
-          fhPrimY  ->Fill(pi0Y) ;
-          fhPrimPhi->Fill(phi) ;
-          
           //Check if both photons hit Calorimeter
+          if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
           Int_t iphot1=prim->GetFirstDaughter() ;
           Int_t iphot2=prim->GetLastDaughter() ;
           if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
@@ -618,14 +899,19 @@ void AliAnaPi0::FillAcceptanceHistograms(){
               }          
               
               if(inacceptance){
-                
-                fhPrimAccPt->Fill(pi0Pt) ;
-                fhPrimAccPhi->Fill(phi) ;
-                fhPrimAccY->Fill(pi0Y) ;
-                Double_t angle  = lv1.Angle(lv2.Vect());
-                fhPrimOpeningAngle   ->Fill(pi0Pt,angle);
-                fhPrimCosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
-                
+                if(pdg==111){
+                  fhPrimPi0AccPt->Fill(pi0Pt) ;
+                  fhPrimPi0AccPhi->Fill(phi) ;
+                  fhPrimPi0AccY->Fill(pi0Y) ;
+                  Double_t angle  = lv1.Angle(lv2.Vect());
+                  fhPrimPi0OpeningAngle   ->Fill(pi0Pt,angle);
+                  fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
+                }
+                else if(pdg==221){
+                  fhPrimEtaAccPt->Fill(pi0Pt) ;
+                  fhPrimEtaAccPhi->Fill(phi) ;
+                  fhPrimEtaAccY->Fill(pi0Y) ;
+                }
               }//Accepted
             }// 2 photons      
           }//Check daughters exist
@@ -634,39 +920,348 @@ void AliAnaPi0::FillAcceptanceHistograms(){
     }//stack exists and data is MC
   }//read stack
   else if(GetReader()->ReadAODMCParticles()){
-    if(GetDebug() >= 0)  printf("AliAnaPi0::FillAcceptanceHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
-  }    
+    
+    TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
+    if(mcparticles){
+      Int_t nprim = mcparticles->GetEntriesFast();
+      for(Int_t i=0 ; i < nprim; i++){
+        AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);        
+        Int_t pdg = prim->GetPdgCode();
+        if( pdg == 111 || pdg == 221){
+          Double_t pi0Pt = prim->Pt() ;
+          //printf("pi0, pt %2.2f\n",pi0Pt);
+          if(prim->E() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception      
+          Double_t pi0Y  = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
+          Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
+          if(pdg == 111){
+            if(TMath::Abs(pi0Y) < 0.5){
+              fhPrimPi0Pt->Fill(pi0Pt) ;
+            }
+            fhPrimPi0Y  ->Fill(pi0Y) ;
+            fhPrimPi0Phi->Fill(phi) ;
+          }
+          else if(pdg == 221){
+            if(TMath::Abs(pi0Y) < 0.5){
+              fhPrimEtaPt->Fill(pi0Pt) ;
+            }
+            fhPrimEtaY  ->Fill(pi0Y) ;
+            fhPrimEtaPhi->Fill(phi) ;
+          }
+          //Check if both photons hit Calorimeter
+          if(prim->GetNDaughters()!=2) return; //Only interested in 2 gamma decay
+          Int_t iphot1=prim->GetDaughter(0) ;
+          Int_t iphot2=prim->GetDaughter(1) ;
+          if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim){
+            AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);   
+            AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);   
+            if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
+              //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
+              //       phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
+              
+              TLorentzVector lv1, lv2;
+              lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
+              lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
+
+              Bool_t inacceptance = kFALSE;
+              if(fCalorimeter == "PHOS"){
+                if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
+                  Int_t mod ;
+                  Double_t x,z ;
+                  Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
+                  Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
+                  if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) && 
+                     GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x)) 
+                    inacceptance = kTRUE;
+                  if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+                }
+                else{
+                  
+                  if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter)) 
+                    inacceptance = kTRUE ;
+                  if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+                }
+                
+              }           
+              else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
+                if(GetEMCALGeometry()){
+                  TVector3 vtx(phot1->Xv(),phot1->Yv(),phot1->Zv());
+                  TVector3 vimpact(0,0,0);
+                  Int_t absID1=0;
+                  
+                  GetEMCALGeometry()->ImpactOnEmcal(vtx,phot1->Theta(),phot1->Phi(),absID1,vimpact);
+                  TVector3 vtx2(phot2->Xv(),phot2->Yv(),phot2->Zv());
+                  TVector3 vimpact2(0,0,0);
+                  Int_t absID2=0;
+                  GetEMCALGeometry()->ImpactOnEmcal(vtx2,phot2->Theta(),phot2->Phi(),absID2,vimpact2);
+//                  if(TMath::Abs(phot1->Eta()) < 0.7 && phot1->Phi() > 80*TMath::DegToRad() && phot1->Phi() < 120*TMath::DegToRad() ) 
+//                    printf("Phot1 accepted? %d\n",absID1);
+//                  if(TMath::Abs(phot2->Eta()) < 0.7 && phot2->Phi() > 80*TMath::DegToRad() && phot2->Phi() < 120*TMath::DegToRad() ) 
+//                    printf("Phot2 accepted? %d\n",absID2);
+
+                  if( absID1 >= 0 && absID2 >= 0) 
+                    inacceptance = kTRUE;
+                  if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+                }
+                else{
+                  if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter)) 
+                    inacceptance = kTRUE ;
+                  if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
+                }
+              }          
+              
+              if(inacceptance){
+                if(pdg==111){
+                  fhPrimPi0AccPt->Fill(pi0Pt) ;
+                  fhPrimPi0AccPhi->Fill(phi) ;
+                  fhPrimPi0AccY->Fill(pi0Y) ;
+                  Double_t angle  = lv1.Angle(lv2.Vect());
+                  fhPrimPi0OpeningAngle   ->Fill(pi0Pt,angle);
+                  fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
+                }
+                else if(pdg==221){
+                  fhPrimEtaAccPt->Fill(pi0Pt) ;
+                  fhPrimEtaAccPhi->Fill(phi) ;
+                  fhPrimEtaAccY->Fill(pi0Y) ;
+                }
+              }//Accepted
+            }// 2 photons      
+          }//Check daughters exist
+        }// Primary pi0
+      }//loop on primaries     
+    }//stack exists and data is MC
+    
+    
+  }    // read AOD MC
 }
 
+//_____________________________________________________________
+void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1,  const Int_t index2, 
+                                              const Float_t pt1,   const Float_t pt2, 
+                                              const Int_t ncell1,  const Int_t ncell2,
+                                              const Double_t mass, const Double_t pt, const Double_t asym, 
+                                              const Double_t deta, const Double_t dphi){
+  //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
+  //Adjusted for Pythia, need to see what to do for other generators.
+  //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles, 
+  // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
+  
+  Int_t ancPDG    = 0;
+  Int_t ancStatus = 0;
+  TLorentzVector ancMomentum;
+  Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2, 
+                                                              GetReader(), ancPDG, ancStatus,ancMomentum);
+  
+  if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
+                            ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
+    
+    if(ancLabel > -1){
+      if(ancPDG==22){//gamma
+        fhMCOrgMass[0]->Fill(pt,mass);
+        fhMCOrgAsym[0]->Fill(pt,asym);
+        fhMCOrgDeltaEta[0]->Fill(pt,deta);
+        fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
+      }              
+      else if(TMath::Abs(ancPDG)==11){//e
+        fhMCOrgMass[1]->Fill(pt,mass);
+        fhMCOrgAsym[1]->Fill(pt,asym);
+        fhMCOrgDeltaEta[1]->Fill(pt,deta);
+        fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
+      }          
+      else if(ancPDG==111){//Pi0
+        fhMCOrgMass[2]->Fill(pt,mass);
+        fhMCOrgAsym[2]->Fill(pt,asym);
+        fhMCOrgDeltaEta[2]->Fill(pt,deta);
+        fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
+        if(fMultiCutAnaSim){
+          for(Int_t ipt=0; ipt<fNPtCuts; ipt++){          
+            for(Int_t icell=0; icell<fNCellNCuts; icell++){
+              for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
+                Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
+                if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        && 
+                   asym   <  fAsymCuts[iasym]                                   && 
+                   ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){ 
+                  fhMCPi0MassPtRec [index]->Fill(pt,mass);
+                  fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
+                  if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
+                }//pass the different cuts
+              }// pid bit cut loop
+            }// icell loop
+          }// pt cut loop
+        }//Multi cut ana sim
+        else {
+          fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
+          if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt); 
+        }
+      }
+      else if(ancPDG==221){//Eta
+        fhMCOrgMass[3]->Fill(pt,mass);
+        fhMCOrgAsym[3]->Fill(pt,asym);
+        fhMCOrgDeltaEta[3]->Fill(pt,deta);
+        fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
+        if(fMultiCutAnaSim){
+          for(Int_t ipt=0; ipt<fNPtCuts; ipt++){          
+            for(Int_t icell=0; icell<fNCellNCuts; icell++){
+              for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
+                Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
+                if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        && 
+                   asym   <  fAsymCuts[iasym]                                   && 
+                   ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){ 
+                  fhMCEtaMassPtRec [index]->Fill(pt,mass);
+                  fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
+                  if(mass < 0.17 && mass > 0.1) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
+                }//pass the different cuts
+              }// pid bit cut loop
+            }// icell loop
+          }// pt cut loop
+        } //Multi cut ana sim
+        else {
+          fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
+          if(mass < 0.17 && mass > 0.1) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt); 
+        }
+      }
+      else if(ancPDG==-2212){//AProton
+        fhMCOrgMass[4]->Fill(pt,mass);
+        fhMCOrgAsym[4]->Fill(pt,asym);
+        fhMCOrgDeltaEta[4]->Fill(pt,deta);
+        fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
+      }   
+      else if(ancPDG==-2112){//ANeutron
+        fhMCOrgMass[5]->Fill(pt,mass);
+        fhMCOrgAsym[5]->Fill(pt,asym);
+        fhMCOrgDeltaEta[5]->Fill(pt,deta);
+        fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
+      }       
+      else if(TMath::Abs(ancPDG)==13){//muons
+        fhMCOrgMass[6]->Fill(pt,mass);
+        fhMCOrgAsym[6]->Fill(pt,asym);
+        fhMCOrgDeltaEta[6]->Fill(pt,deta);
+        fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
+      }                   
+      else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
+        if(ancStatus==1){//Stable particles, converted? not decayed resonances
+          fhMCOrgMass[6]->Fill(pt,mass);
+          fhMCOrgAsym[6]->Fill(pt,asym);
+          fhMCOrgDeltaEta[6]->Fill(pt,deta);
+          fhMCOrgDeltaPhi[6]->Fill(pt,dphi);  
+        }
+        else{//resonances and other decays, more hadron conversions?
+          fhMCOrgMass[7]->Fill(pt,mass);
+          fhMCOrgAsym[7]->Fill(pt,asym);
+          fhMCOrgDeltaEta[7]->Fill(pt,deta);
+          fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
+        }
+      }
+      else {//Partons, colliding protons, strings, intermediate corrections
+        if(ancStatus==11 || ancStatus==12){//String fragmentation
+          fhMCOrgMass[8]->Fill(pt,mass);
+          fhMCOrgAsym[8]->Fill(pt,asym);
+          fhMCOrgDeltaEta[8]->Fill(pt,deta);
+          fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
+        }
+        else if (ancStatus==21){
+          if(ancLabel < 2) {//Colliding protons
+            fhMCOrgMass[11]->Fill(pt,mass);
+            fhMCOrgAsym[11]->Fill(pt,asym);
+            fhMCOrgDeltaEta[11]->Fill(pt,deta);
+            fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
+          }//colliding protons  
+          else if(ancLabel < 6){//partonic initial states interactions
+            fhMCOrgMass[9]->Fill(pt,mass);
+            fhMCOrgAsym[9]->Fill(pt,asym);
+            fhMCOrgDeltaEta[9]->Fill(pt,deta);
+            fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
+          }
+          else if(ancLabel < 8){//Final state partons radiations?
+            fhMCOrgMass[10]->Fill(pt,mass);
+            fhMCOrgAsym[10]->Fill(pt,asym);
+            fhMCOrgDeltaEta[10]->Fill(pt,deta);
+            fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
+          }
+          else {
+            printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
+                   ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
+          }
+        }//status 21
+        else {
+          printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
+                 ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
+        }
+      }////Partons, colliding protons, strings, intermediate corrections
+    }//ancLabel > -1 
+    else { //ancLabel <= -1
+      //printf("Not related at all label = %d\n",ancLabel);
+      fhMCOrgMass[12]->Fill(pt,mass);
+      fhMCOrgAsym[12]->Fill(pt,asym);
+      fhMCOrgDeltaEta[12]->Fill(pt,deta);
+      fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
+    }
+}  
+
 //____________________________________________________________________________________________________________________________________________________
 void AliAnaPi0::MakeAnalysisFillHistograms() 
 {
   //Process one event and extract photons from AOD branch 
   // filled with AliAnaPhoton and fill histos with invariant mass
   
-  //In case of MC data, fill acceptance histograms
-  FillAcceptanceHistograms();
-  
-  //Apply some cuts on event: vertex position and centrality range  
-  Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
-  if(IsBadRun(iRun)) return ;  
+  //In case of simulated data, fill acceptance histograms
+  if(IsDataMC())FillAcceptanceHistograms();
+
+  //Init some variables
+//Int_t   iRun     = (GetReader()->GetInputEvent())->GetRunNumber() ;
+  Int_t   nPhot    = GetInputAODBranch()->GetEntriesFast() ;
+  Int_t   nClus    = 0;
+  Int_t   nCell    = 0;
+  Float_t eClusTot = 0;
+  Float_t eCellTot = 0;
   
-  Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
+  // Count the number of clusters and cells, in case multiplicity bins dependent on such numbers
+  // are requested
+  if(fCalorimeter=="EMCAL"){ 
+    nClus = GetAODEMCAL()  ->GetEntriesFast();
+    nCell = GetEMCALCells()->GetNumberOfCells();
+    for(Int_t icl=0; icl < nClus; icl++) eClusTot += ((AliVCluster*) GetAODEMCAL()->At(icl))->E();
+    for(Int_t jce=0; jce < nCell; jce++) eCellTot +=  GetEMCALCells()->GetAmplitude(jce);
+
+  }
+  else {                     
+    nClus = GetAODPHOS()  ->GetEntriesFast();
+    nCell = GetPHOSCells()->GetNumberOfCells();
+    for(Int_t icl=0; icl < nClus; icl++) eClusTot +=  ((AliVCluster*)GetAODPHOS()->At(icl))->E();
+    for(Int_t jce=0; jce < nCell; jce++) eCellTot +=  GetPHOSCells()->GetAmplitude(jce);
+  }
+
+  //Fill the average number of cells or clusters per SM
+  eClusTot /=fNModules;
+  eCellTot /=fNModules;
+  fhAverTotECluster->Fill(eClusTot);
+  fhAverTotECell   ->Fill(eCellTot);
+
   if(GetDebug() > 1) 
     printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
-  if(nPhot < 2 )
-    return ; 
-  Int_t module1 = -1;
-  Int_t module2 = -1;
-  Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex 
-  Int_t evtIndex1 = 0 ; 
-  Int_t currentEvtIndex = -1 ; 
-  Int_t curCentrBin = 0 ; 
-  Int_t curRPBin    = 0 ; 
-  Int_t curZvertBin = 0 ;
+
+  //If less than photon 2 entries in the list, skip this event
+  if(nPhot < 2 ) return ; 
   
+  if(GetDebug() > 1) 
+    printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
+  
+  //Init variables
+  Int_t module1         = -1;
+  Int_t module2         = -1;
+  Double_t vert[]       = {0.0, 0.0, 0.0} ; //vertex 
+  Int_t evtIndex1       = 0 ; 
+  Int_t currentEvtIndex = -1; 
+  Int_t curCentrBin     = 0 ; 
+  Int_t curRPBin        = 0 ; 
+  Int_t curZvertBin     = 0 ;
+  
+  //---------------------------------
+  //First loop on photons/clusters
+  //---------------------------------
   for(Int_t i1=0; i1<nPhot-1; i1++){
     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
+    //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
+
     // get the event index in the mixed buffer where the photon comes from 
     // in case of mixing with analysis frame, not own mixing
     evtIndex1 = GetEventIndex(p1, vert) ; 
@@ -676,25 +1271,87 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
     if ( evtIndex1 == -2 )
       continue ; 
     if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ;   //vertex cut
+    
+    //----------------------------------------------------------------------------
+    // Get the multiplicity bin. Different cases: centrality (PbPb), 
+    // average cluster multiplicity, average cell multiplicity, track multiplicity 
+    // default is centrality bins
+    //----------------------------------------------------------------------------
     if (evtIndex1 != currentEvtIndex) {
-      curCentrBin = GetEventCentrality();
-      //In case too few centrality bins defined with respect the selected centrality class, put in the last bin the higher centrality classes
-      if(curCentrBin >= fNCentrBin) curCentrBin=fNCentrBin-1;
+      if(fUseTrackMultBins){ // Track multiplicity bins
+        //printf("track  mult %d\n",GetTrackMultiplicity());
+        curCentrBin = (GetTrackMultiplicity()-1)/5; 
+        if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+        //printf("track mult bin %d\n",curCentrBin);
+      }
+      else if(fUsePhotonMultBins){ // Photon multiplicity bins
+        //printf("photon  mult %d cluster mult %d\n",nPhot, nClus);
+        curCentrBin = (nClus-1)/3; 
+        if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+        //printf("photon mult bin %d\n",curCentrBin);        
+      }
+      else if(fUseAverClusterEBins){ // Cluster multiplicity bins
+        //Bins for pp, if needed can be done in a more general way
+        if     (eClusTot < 0.5 )curCentrBin = 0;
+        else if(eClusTot < 1.0) curCentrBin = 1;
+        else if(eClusTot < 1.5) curCentrBin = 2;
+        else if(eClusTot < 2.0) curCentrBin = 3;
+        else if(eClusTot < 3.0) curCentrBin = 4;
+        else if(eClusTot < 4.0) curCentrBin = 5;
+        else if(eClusTot < 5.0) curCentrBin = 6;
+        else if(eClusTot < 7.5) curCentrBin = 7;
+        else if(eClusTot < 10.) curCentrBin = 8;
+        else if(eClusTot < 15.) curCentrBin = 9;
+        else if(eClusTot < 20.) curCentrBin = 10;
+        if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+        //printf("cluster E average %f, bin %d \n",eClusTot,curCentrBin);
+      }
+      else if(fUseAverCellEBins){ // Cell multiplicity bins
+        //Bins for pp, if needed can be done in a more general way
+        if     (eCellTot < 0.5) curCentrBin = 0;
+        else if(eCellTot < 1.0) curCentrBin = 1;
+        else if(eCellTot < 1.5) curCentrBin = 2;
+        else if(eCellTot < 2.0) curCentrBin = 3;
+        else if(eCellTot < 3.0) curCentrBin = 4;
+        else if(eCellTot < 4.0) curCentrBin = 5;
+        else if(eCellTot < 5.0) curCentrBin = 6;
+        else if(eCellTot < 7.5) curCentrBin = 7;
+        else if(eCellTot < 10.) curCentrBin = 8;
+        else if(eCellTot < 15.) curCentrBin = 9;
+        else if(eCellTot < 20.) curCentrBin = 10; 
+        if(curCentrBin > fNCentrBin-1) curCentrBin=fNCentrBin-1;
+        //printf("cell E average %f, bin %d \n",eCellTot,curCentrBin);
+      }
+      else { //Event centrality
+        curCentrBin = GetEventCentrality();
+      }
+
+      //Get vertex z bin
       curRPBin    = 0 ;
       curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
+      
+      //Fill event bin info
       fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
       currentEvtIndex = evtIndex1 ; 
       if(GetDebug() > 1) 
-        printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d\n",curCentrBin,curRPBin,curZvertBin);
+        printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
     }
     
     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
     
+    //Get the momentum of this cluster
     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
-    //Get Module number
+    
+    //Get (Super)Module number of this cluster
     module1 = GetModuleNumber(p1);
+    
+    //---------------------------------
+    //Second loop on photons/clusters
+    //---------------------------------
     for(Int_t i2=i1+1; i2<nPhot; i2++){
       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
+      
+      //In case of mixing frame, check we are not in the same event as the first cluster
       Int_t evtIndex2 = GetEventIndex(p2, vert) ; 
       if ( evtIndex2 == -1 )
         return ; 
@@ -702,35 +1359,59 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
         continue ;    
       if (GetMixedEvent() && (evtIndex1 == evtIndex2))
         continue ;
+      
       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
+      //Get the momentum of this cluster
       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
       //Get module number
-      module2 = GetModuleNumber(p2);
-      Double_t m  = (photon1 + photon2).M() ;
-      Double_t pt = (photon1 + photon2).Pt();
-      Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
+      module2       = GetModuleNumber(p2);
+      
+      //---------------------------------
+      // Get pair kinematics
+      //---------------------------------
+      Double_t m    = (photon1 + photon2).M() ;
+      Double_t pt   = (photon1 + photon2).Pt();
+      Double_t deta = photon1.Eta() - photon2.Eta();
+      Double_t dphi = photon1.Phi() - photon2.Phi();
+      Double_t a    = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
+      
       if(GetDebug() > 2)
-        printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
-               p1->Pt(), p2->Pt(), pt,m,a);
+        printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
+      
+      //--------------------------------
+      // Opening angle selection
+      //--------------------------------
       //Check if opening angle is too large or too small compared to what is expected  
       Double_t angle   = photon1.Angle(photon2.Vect());
-      //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
-      //printf("angle %f\n",angle);
-      if(fUseAngleCut && angle < 0.1) 
+      if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
+        if(GetDebug() > 2)
+          printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
         continue;
+      }
 
+      if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
+        if(GetDebug() > 2)
+          printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
+        continue;
+      }
+
+      //-------------------------------------------------------------------------------------------------
       //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
+      //-------------------------------------------------------------------------------------------------
       if(a < fAsymCuts[0]){
         if(module1==module2 && module1 >=0 && module1<fNModules)
           fhReMod[module1]->Fill(pt,m) ;
         else  
-          fhReDiffMod[fNModules]->Fill(pt,m) ;
+          fhReDiffMod[fNModules+2]->Fill(pt,m) ;
         
         if(fCalorimeter=="EMCAL"){
           if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[0]->Fill(pt,m) ; 
           if((module1==1 && module2==3) || (module1==3 && module2==1)) fhReDiffMod[1]->Fill(pt,m) ; 
           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[2]->Fill(pt,m) ;
-          if((module1==2 && module2==3) || (module1==3 && module2==2)) fhReDiffMod[3]->Fill(pt,m) ; 
+          if((module1==2 && module2==3) || (module1==3 && module2==2)) fhReDiffMod[3]->Fill(pt,m) ;
+          if((module1==0 && module2==3) || (module1==3 && module2==0)) fhReDiffMod[4]->Fill(pt,m) ;
+          if((module1==2 && module2==1) || (module1==1 && module2==2)) fhReDiffMod[5]->Fill(pt,m) ;
         }
         else {
           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[0]->Fill(pt,m) ; 
@@ -743,25 +1424,32 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
       Bool_t ok = kTRUE;
       if(fSameSM && module1!=module2) ok=kFALSE;
       if(ok){
+        
+        //Check if one of the clusters comes from a conversion 
+        if     (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
+        else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
+        
         //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
         for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
           if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){ 
             for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
               if(a < fAsymCuts[iasym]){
                 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
-                //printf("cen %d, pid %d, asy %d, Index %d\n",curCentrBin,ipid,iasym,index);
+                //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym);
                 fhRe1     [index]->Fill(pt,m);
                 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
-                if(p1->DistToBad()>0 && p2->DistToBad()>0){
-                  fhRe2     [index]->Fill(pt,m) ;
-                  if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
-                  if(p1->DistToBad()>1 && p2->DistToBad()>1){
-                    fhRe3     [index]->Fill(pt,m) ;
-                    if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
-                  }//assymetry cut
-                }// asymmetry cut loop
-              }// bad 3
-            }// bad2
+                if(fFillBadDistHisto){
+                  if(p1->DistToBad()>0 && p2->DistToBad()>0){
+                    fhRe2     [index]->Fill(pt,m) ;
+                    if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
+                    if(p1->DistToBad()>1 && p2->DistToBad()>1){
+                      fhRe3     [index]->Fill(pt,m) ;
+                      if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
+                    }// bad 3
+                  }// bad2
+                }// Fill bad dist histos
+              }//assymetry cut
+            }// asymmetry cut loop
           }// bad 1
         }// pid bit loop
         
@@ -771,24 +1459,16 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
         
         //Fill histograms with pair assymmetry
         fhRePtAsym->Fill(pt,a);
-        if(m > 0.10 && m < 0.16) fhRePtAsymPi0->Fill(pt,a);
+        if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
         if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
         
-        //Multi cuts analysis 
-        if(fMultiCutAna){
-          //Histograms for different PID bits selection
-          for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
-            
-            if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    && 
-               p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
-            
-            //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
-          } // pid bit cut loop
-          
-          //Several pt,ncell and asymmetry cuts
-          //Get the number of cells
-          Int_t ncell1 = 0;
-          Int_t ncell2 = 0;
+        //-------------------------------------------------------
+        //Get the number of cells needed for multi cut analysis.
+        //-------------------------------------------------------        
+        Int_t ncell1 = 0;
+        Int_t ncell2 = 0;
+        if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim)){
+
           AliVEvent * event = GetReader()->GetInputEvent();
           if(event){
             for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
@@ -808,62 +1488,148 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
             }
             //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
           }
+        }
+        
+        //---------
+        // MC data
+        //---------
+        //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
+        if(IsDataMC()) FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);    
+        
+        //-----------------------
+        //Multi cuts analysis 
+        //-----------------------
+        if(fMultiCutAna){
+          //Histograms for different PID bits selection
+          for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
+            
+            if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    && 
+               p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
+            
+            //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
+          } // pid bit cut loop
+          
+          //Several pt,ncell and asymmetry cuts
           for(Int_t ipt=0; ipt<fNPtCuts; ipt++){          
             for(Int_t icell=0; icell<fNCellNCuts; icell++){
               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
                 if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        && 
                    a        <   fAsymCuts[iasym]                                    && 
-                   ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell]) fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
-                
-                //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
+                   ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell]){
+                    fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
+                  //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
+                  if(module1==module2){
+                    if     (module1==0)  fhRePtNCellAsymCutsSM0[index]->Fill(pt,m) ;
+                    else if(module1==1)  fhRePtNCellAsymCutsSM1[index]->Fill(pt,m) ;
+                    else if(module1==2)  fhRePtNCellAsymCutsSM2[index]->Fill(pt,m) ;
+                    else if(module1==3)  fhRePtNCellAsymCutsSM3[index]->Fill(pt,m) ;
+                    else printf("AliAnaPi0::FillHistograms() - WRONG SM NUMBER\n");
+                  }
+                }
               }// pid bit cut loop
             }// icell loop
           }// pt cut loop
           for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
             if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
           }
-          
         }// multiple cuts analysis
       }// ok if same sm
     }// second same event particle
   }// first cluster
-  
+        
+  //-------------------------------------------------------------
+  // Mixing
+  //-------------------------------------------------------------
   if(fDoOwnMix){
-    //Fill mixed
+    //Recover events in with same characteristics as the current event
     TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
     Int_t nMixed = evMixList->GetSize() ;
     for(Int_t ii=0; ii<nMixed; ii++){  
       TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
       Int_t nPhot2=ev2->GetEntriesFast() ;
       Double_t m = -999;
-      if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
+      if(GetDebug() > 1) 
+        printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, curCentrBin);
       
+      //---------------------------------
+      //First loop on photons/clusters
+      //---------------------------------      
       for(Int_t i1=0; i1<nPhot; i1++){
         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
+        if(fSameSM && GetModuleNumber(p1)!=module1) continue;
+        
+        //Get kinematics of cluster and (super) module of this cluster
         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
         module1 = GetModuleNumber(p1);
+        
+        //---------------------------------
+        //First loop on photons/clusters
+        //---------------------------------        
         for(Int_t i2=0; i2<nPhot2; i2++){
           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
           
+          //Get kinematics of second cluster and calculate those of the pair
           TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
-          m =           (photon1+photon2).M() ; 
+          m           = (photon1+photon2).M() ; 
           Double_t pt = (photon1 + photon2).Pt();
           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
           
           //Check if opening angle is too large or too small compared to what is expected
           Double_t angle   = photon1.Angle(photon2.Vect());
-          //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
-          if(fUseAngleCut && angle < 0.1) continue;  
+          if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){ 
+            if(GetDebug() > 2)
+              printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
+            continue;
+          }
+          if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
+            if(GetDebug() > 2)
+              printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
+            continue; 
+          
+          } 
           
           if(GetDebug() > 2)
             printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
                    p1->Pt(), p2->Pt(), pt,m,a);        
+          
           //In case we want only pairs in same (super) module, check their origin.
           module2 = GetModuleNumber(p2);
+          
+          //-------------------------------------------------------------------------------------------------
+          //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
+          //-------------------------------------------------------------------------------------------------          
+          if(a < fAsymCuts[0]){
+            if(module1==module2 && module1 >=0 && module1<fNModules)
+              fhMiMod[module1]->Fill(pt,m) ;
+            else  
+              fhMiDiffMod[fNModules+2]->Fill(pt,m) ;
+            
+            if(fCalorimeter=="EMCAL"){
+              if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffMod[0]->Fill(pt,m) ; 
+              if((module1==1 && module2==3) || (module1==3 && module2==1)) fhMiDiffMod[1]->Fill(pt,m) ; 
+              if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffMod[2]->Fill(pt,m) ;
+              if((module1==2 && module2==3) || (module1==3 && module2==2)) fhMiDiffMod[3]->Fill(pt,m) ; 
+              if((module1==0 && module2==3) || (module1==3 && module2==0)) fhMiDiffMod[4]->Fill(pt,m) ;
+              if((module1==2 && module2==1) || (module1==1 && module2==2)) fhMiDiffMod[5]->Fill(pt,m) ;
+
+            }
+            else {
+              if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffMod[0]->Fill(pt,m) ; 
+              if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffMod[1]->Fill(pt,m) ; 
+              if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffMod[2]->Fill(pt,m) ;
+            }
+          }
+          
           Bool_t ok = kTRUE;
           if(fSameSM && module1!=module2) ok=kFALSE;
           if(ok){
+            
+            //Check if one of the clusters comes from a conversion 
+            if     (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
+            else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
+
+            //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
             for(Int_t ipid=0; ipid<fNPIDBits; ipid++){ 
               if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ 
                 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
@@ -871,23 +1637,52 @@ void AliAnaPi0::MakeAnalysisFillHistograms()
                     Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
                     fhMi1     [index]->Fill(pt,m) ;
                     if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
-                    if(p1->DistToBad()>0 && p2->DistToBad()>0){
-                      fhMi2     [index]->Fill(pt,m) ;
-                      if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
-                      if(p1->DistToBad()>1 && p2->DistToBad()>1){
-                        fhMi3     [index]->Fill(pt,m) ;
-                        if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
+                    if(fFillBadDistHisto){
+                      if(p1->DistToBad()>0 && p2->DistToBad()>0){
+                        fhMi2     [index]->Fill(pt,m) ;
+                        if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
+                        if(p1->DistToBad()>1 && p2->DistToBad()>1){
+                          fhMi3     [index]->Fill(pt,m) ;
+                          if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
+                        }
                       }
-                    }
+                    }// Fill bad dist histo
                   }//Asymmetry cut
                 }// Asymmetry loop
               }//PID cut
             }//loop for histograms
+            
+            //-----------------------
+            //Multi cuts analysis 
+            //-----------------------            
+            if(fMultiCutAna){
+              //Several pt,ncell and asymmetry cuts
+              for(Int_t ipt=0; ipt<fNPtCuts; ipt++){          
+                for(Int_t icell=0; icell<fNCellNCuts; icell++){
+                  for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
+                    Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
+                    if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        && 
+                       a        <   fAsymCuts[iasym]                                    && 
+                       p1->GetBtag() >=  fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell]){
+                      fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
+                      //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
+                    }
+                  }// pid bit cut loop
+                }// icell loop
+              }// pt cut loop
+            } // Multi cut ana
+            
+            //Fill histograms with opening angle
+            fhMixedOpeningAngle   ->Fill(pt,angle);
+            fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));          
           }//ok
         }// second cluster loop
       }//first cluster loop
     }//loop on mixed events
     
+    //--------------------------------------------------------
+    //Add the current event to the list of events for mixing
+    //--------------------------------------------------------
     TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
     //Add current event to buffer and Remove redundant events 
     if(currentEvent->GetEntriesFast()>0){
@@ -924,16 +1719,21 @@ void AliAnaPi0::ReadHistograms(TList* outputList)
   if(!fhMi1) fhMi1 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
   if(!fhMi2) fhMi2 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
   if(!fhMi3) fhMi3 = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;      
-  if(fMakeInvPtPlots){
-    if(!fhReInvPt1) fhReInvPt1  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-    if(!fhReInvPt2) fhReInvPt2  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-    if(!fhReInvPt3) fhReInvPt3  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-    if(!fhMiInvPt1) fhMiInvPt1  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-    if(!fhMiInvPt2) fhMiInvPt2  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
-    if(!fhMiInvPt3) fhMiInvPt3  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ; 
-  }
+  if(!fhReInvPt1) fhReInvPt1  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+  if(!fhReInvPt2) fhReInvPt2  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+  if(!fhReInvPt3) fhReInvPt3  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+  if(!fhMiInvPt1) fhMiInvPt1  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+  if(!fhMiInvPt2) fhMiInvPt2  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;
+  if(!fhMiInvPt3) fhMiInvPt3  = new TH2D*[fNCentrBin*fNPIDBits*fNAsymCuts] ;   
   if(!fhReMod)    fhReMod     = new TH2D*[fNModules]   ;       
   if(!fhReDiffMod)fhReDiffMod = new TH2D*[fNModules+1] ;       
+  if(!fhMiMod)    fhReMod     = new TH2D*[fNModules]   ;       
+  if(!fhMiDiffMod)fhReDiffMod = new TH2D*[fNModules+1] ;       
+  
+  fhReConv  = (TH2D*) outputList->At(index++);
+  fhMiConv  = (TH2D*) outputList->At(index++);
+  fhReConv2 = (TH2D*) outputList->At(index++);
+  fhMiConv2 = (TH2D*) outputList->At(index++);
 
   for(Int_t ic=0; ic<fNCentrBin; ic++){
     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
@@ -943,21 +1743,19 @@ void AliAnaPi0::ReadHistograms(TList* outputList)
         fhRe1[ihisto] = (TH2D*) outputList->At(index++);
         fhRe2[ihisto] = (TH2D*) outputList->At(index++);
         fhRe3[ihisto] = (TH2D*) outputList->At(index++);
-        if(fMakeInvPtPlots){
-          fhReInvPt1[ihisto] = (TH2D*) outputList->At(index++);
-          fhReInvPt2[ihisto] = (TH2D*) outputList->At(index++);
-          fhReInvPt3[ihisto] = (TH2D*) outputList->At(index++);
-        }
+      
+        fhReInvPt1[ihisto] = (TH2D*) outputList->At(index++);
+        fhReInvPt2[ihisto] = (TH2D*) outputList->At(index++);
+        fhReInvPt3[ihisto] = (TH2D*) outputList->At(index++);
       
         if(fDoOwnMix){
           fhMi1[ihisto] = (TH2D*) outputList->At(index++);
           fhMi2[ihisto] = (TH2D*) outputList->At(index++);
           fhMi3[ihisto] = (TH2D*) outputList->At(index++);
-          if(fMakeInvPtPlots){
-            fhMiInvPt1[ihisto] = (TH2D*) outputList->At(index++);
-            fhMiInvPt2[ihisto] = (TH2D*) outputList->At(index++);
-            fhMiInvPt3[ihisto] = (TH2D*) outputList->At(index++); 
-          }
+      
+          fhMiInvPt1[ihisto] = (TH2D*) outputList->At(index++);
+          fhMiInvPt2[ihisto] = (TH2D*) outputList->At(index++);
+          fhMiInvPt3[ihisto] = (TH2D*) outputList->At(index++); 
         }//Own mix
       }//asymmetry loop
     }// pid loop
@@ -993,20 +1791,66 @@ void AliAnaPi0::ReadHistograms(TList* outputList)
   
   fhRealOpeningAngle     = (TH2D*)  outputList->At(index++);
   fhRealCosOpeningAngle  = (TH2D*)  outputList->At(index++);
+  if(fDoOwnMix){
+    fhMixedOpeningAngle     = (TH2D*)  outputList->At(index++);
+    fhMixedCosOpeningAngle  = (TH2D*)  outputList->At(index++);
+  }
   
   //Histograms filled only if MC data is requested     
   if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
-    fhPrimPt     = (TH1D*)  outputList->At(index++);
-    fhPrimAccPt  = (TH1D*)  outputList->At(index++);
-    fhPrimY      = (TH1D*)  outputList->At(index++);
-    fhPrimAccY   = (TH1D*)  outputList->At(index++);
-    fhPrimPhi    = (TH1D*)  outputList->At(index++);
-    fhPrimAccPhi = (TH1D*)  outputList->At(index++);
+    fhPrimPi0Pt     = (TH1D*)  outputList->At(index++);
+    fhPrimPi0AccPt  = (TH1D*)  outputList->At(index++);
+    fhPrimPi0Y      = (TH1D*)  outputList->At(index++);
+    fhPrimPi0AccY   = (TH1D*)  outputList->At(index++);
+    fhPrimPi0Phi    = (TH1D*)  outputList->At(index++);
+    fhPrimPi0AccPhi = (TH1D*)  outputList->At(index++);
+    for(Int_t i = 0; i<13; i++){
+      fhMCOrgMass[i]     = (TH2D*)  outputList->At(index++);
+      fhMCOrgAsym[i]     = (TH2D*)  outputList->At(index++);
+      fhMCOrgDeltaEta[i] = (TH2D*)  outputList->At(index++);
+      fhMCOrgDeltaPhi[i] = (TH2D*)  outputList->At(index++);
+    }
+    
+    if(fMultiCutAnaSim){
+      fhMCPi0MassPtTrue  = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      fhMCPi0MassPtRec   = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      fhMCPi0PtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      fhMCEtaMassPtTrue  = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      fhMCEtaMassPtRec   = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      fhMCEtaPtTruePtRec = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
+      for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
+        for(Int_t icell=0; icell<fNCellNCuts; icell++){
+          for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
+            Int_t in = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
+            fhMCPi0MassPtTrue[in]  = (TH2D*)  outputList->At(index++);
+            fhMCPi0PtTruePtRec[in] = (TH2D*)  outputList->At(index++);
+            fhMCEtaMassPtTrue[in]  = (TH2D*)  outputList->At(index++);
+            fhMCEtaPtTruePtRec[in] = (TH2D*)  outputList->At(index++);
+          }
+        }
+      }
+    }
+    else{
+      fhMCPi0MassPtTrue  = new TH2D*[1];
+      fhMCPi0PtTruePtRec = new TH2D*[1];
+      fhMCEtaMassPtTrue  = new TH2D*[1];
+      fhMCEtaPtTruePtRec = new TH2D*[1];
+      
+      fhMCPi0MassPtTrue[0]  = (TH2D*)  outputList->At(index++);
+      fhMCPi0PtTruePtRec[0] = (TH2D*)  outputList->At(index++);
+      fhMCEtaMassPtTrue[0]  = (TH2D*)  outputList->At(index++);
+      fhMCEtaPtTruePtRec[0] = (TH2D*)  outputList->At(index++);
+    }
   }
   
-  for(Int_t imod=0; imod < fNModules; imod++)
-    fhReMod[imod] = (TH2D*) outputList->At(index++);
-  
+  for(Int_t imod=0; imod < fNModules; imod++){
+    fhReMod[imod]     = (TH2D*) outputList->At(index++);
+    fhReDiffMod[imod] = (TH2D*) outputList->At(index++);
+    if(fDoOwnMix){
+      fhMiMod[imod]     = (TH2D*) outputList->At(index++);
+      fhMiDiffMod[imod] = (TH2D*) outputList->At(index++);
+    }
+  }
   
 }
 
index 982c533..97a8a53 100755 (executable)
@@ -5,8 +5,10 @@
 /* $Id: $ */
 
 //_________________________________________________________________________
-// Class to fill two-photon invariant mass hisograms 
+// Class to fill two-photon invariant mass histograms 
 // to be used to extract pi0 raw yield.
+// Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles), 
+// it will do nothing if executed alone
 //
 //-- Author: Dmitri Peressounko (RRC "KI")
 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
@@ -35,134 +37,233 @@ class AliAnaPi0 : public AliAnaPartCorrBaseClass {
   
  public:
 
+  //-------------------------------
+  // General analysis frame methods
+  //-------------------------------
+
   TObjString * GetAnalysisCuts();
   TList      * GetCreateOutputObjects(); 
-  
-  void Print(const Option_t * opt) const;
-  
-  //void Init();
-  void InitParameters();
-  
-  void FillAcceptanceHistograms();
+  void Terminate(TList* outputList);
+  void ReadHistograms(TList * outputList); //Fill histograms with histograms in ouput list, needed in Terminate.
+  void Print(const Option_t * opt) const;  
   //void MakeAnalysisFillAOD() {;} //Not needed
   void MakeAnalysisFillHistograms();
-  
-  //    void SetBadRunsList(){;} ;     //Set list of runs which can be used for this analysis
-  //To be defined in future.
-    
-  //Setters for parameters of event buffers
-  void SetNCentrBin(Int_t n=5) {fNCentrBin=n ;} //number of bins in centrality 
-//  void SetNZvertBin(Int_t n=5) {fNZvertBin=n ;} //number of bins for vertex position
-//  void SetNRPBin(Int_t n=6)    {fNrpBin=n ;}    //number of bins in reaction plain
-  void SetNMaxEvMix(Int_t n=20){fNmaxMixEv=n ;} //Maximal number of events for mixing
-  
-  //Setters for event selection
- // void SetZvertexCut(Float_t zcut=40.){fZvtxCut=zcut ;} //cut on vertex position
-  
+  //void Init();
+  void InitParameters();
+
+  //Calorimeter options
   TString GetCalorimeter()   const {return fCalorimeter ; }
   void SetCalorimeter(TString & det)    {fCalorimeter = det ; }
-       
-  void Terminate(TList* outputList);
-  void ReadHistograms(TList * outputList); //Fill histograms with histograms in ouput list, needed in Terminate.
-       
   void SetNumberOfModules(Int_t nmod) {fNModules = nmod;}
-       
-  void SwitchOnAngleSelection()    {fUseAngleCut = kTRUE ; }
-  void SwitchOffAngleSelection()   {fUseAngleCut = kFALSE ; }
   
+  //-------------------------------
+  // EVENT Bin Methods
+  //-------------------------------
+
   virtual Int_t GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)  ;
 
-  void SwitchOnInvPtWeight()    {fMakeInvPtPlots = kTRUE  ; }
-  void SwitchOffInvPtWeight()   {fMakeInvPtPlots = kFALSE ; }
+  //Setters for parameters of event buffers
+  void SetNCentrBin(Int_t n=5)    {fNCentrBin=n ;} //number of bins in centrality 
+//void SetNZvertBin(Int_t n=5)    {fNZvertBin=n ;} //number of bins for vertex position
+//void SetNRPBin(Int_t n=6)       {fNrpBin=n ;}    //number of bins in reaction plain
+  void SetNMaxEvMix(Int_t n=20)   {fNmaxMixEv=n ;} //Maximal number of events for mixing
+  
+  //Switchs for event multiplicity bin option, by default, centrality
+  void SwitchOnTrackMultBins()    {fUseTrackMultBins = kTRUE  ; }
+  void SwitchOffTrackMultBins()   {fUseTrackMultBins = kFALSE ; }
+  
+  void SwitchOnPhotonMultBins()   {fUsePhotonMultBins = kTRUE  ; }
+  void SwitchOffPhotonMultBins()  {fUsePhotonMultBins = kFALSE ; }
+  
+  void SwitchOnClusterEBins()     {fUseAverClusterEBins = kTRUE  ; }
+  void SwitchOffClusterEBins()    {fUseAverClusterEBins = kFALSE ; }
   
-  void SwitchOnOwnMix()    {fDoOwnMix = kTRUE  ; }
-  void SwitchOffOwnMix()   {fDoOwnMix = kFALSE ; }
+  void SwitchOnCellEBins()        {fUseAverCellEBins = kTRUE  ; }
+  void SwitchOffCellEBins()       {fUseAverCellEBins = kFALSE ; }
 
-  void SwitchOnSameSM()    {fSameSM = kTRUE  ; }
-  void SwitchOffSameSM()   {fSameSM = kFALSE ; }
+  //-------------------------------
+       //Opening angle pair selection
+  //-------------------------------
+  void SwitchOnAngleSelection()      {fUseAngleCut = kTRUE      ; }
+  void SwitchOffAngleSelection()     {fUseAngleCut = kFALSE     ; }
+  void SwitchOnAngleEDepSelection()  {fUseAngleEDepCut = kTRUE  ; }
+  void SwitchOffAngleEDepSelection() {fUseAngleEDepCut = kFALSE ; }
+  void SetAngleCut(Float_t a)        {fAngleCut    = a          ; }
+  void SetAngleMaxCut(Float_t a)     {fAngleMaxCut = a          ; }
+
+  //-------------------------------
+  // Use mixing code of this class
+  //-------------------------------
+  void SwitchOnOwnMix()              {fDoOwnMix = kTRUE  ; }
+  void SwitchOffOwnMix()             {fDoOwnMix = kFALSE ; }
+
+  //------------------------------------------
+  //Do analysis only with clusters in same SM
+  //------------------------------------------
+  void SwitchOnSameSM()              {fSameSM = kTRUE    ; }
+  void SwitchOffSameSM()             {fSameSM = kFALSE   ; }
+  
+  //-------------------------------
+  //Histogram filling options off by default
+  //-------------------------------
+  void SwitchOnInvPtWeight()         {fMakeInvPtPlots = kTRUE  ; }
+  void SwitchOffInvPtWeight()        {fMakeInvPtPlots = kFALSE ; }
   
-  //Cuts for multiple analysis
-  void SwitchOnMultipleCutAnalysis()   {fMultiCutAna = kTRUE;}
-  void SwitchOffMultipleCutAnalysis()  {fMultiCutAna = kFALSE;}
+  void SwitchOnFillBadDistHisto()    {fFillBadDistHisto    = kTRUE;}
+  void SwitchOffFillBadDistHisto()   {fFillBadDistHisto    = kFALSE;}
+  
+  //-------------------------------------------
+  //Cuts for multiple analysis, off by default
+  //-------------------------------------------
+  void SwitchOnMultipleCutAnalysis()          {fMultiCutAna    = kTRUE;}
+  void SwitchOffMultipleCutAnalysis()         {fMultiCutAna    = kFALSE;}
 
-  void SetNPtCuts   (Int_t size)  {if(size <= 10)fNPtCuts    = size; }
-  void SetNAsymCuts (Int_t size)  {if(size <= 10)fNAsymCuts  = size; }
-  void SetNNCellCuts(Int_t size)  {if(size <= 10)fNCellNCuts = size; }
-  void SetNPIDBits  (Int_t size)  {if(size <= 10)fNPIDBits   = size; }
+  void SetNPtCuts   (Int_t size)              {if(size <= 10)fNPtCuts    = size; }
+  void SetNAsymCuts (Int_t size)              {if(size <= 10)fNAsymCuts  = size; }
+  void SetNNCellCuts(Int_t size)              {if(size <= 10)fNCellNCuts = size; }
+  void SetNPIDBits  (Int_t size)              {if(size <= 10)fNPIDBits   = size; }
   
   void SetPtCutsAt   (Int_t pos,Float_t val)  {if(pos < 10)fPtCuts[pos]    = val;}
   void SetAsymCutsAt (Int_t pos,Float_t val)  {if(pos < 10)fAsymCuts[pos]  = val;}
   void SetNCellCutsAt(Int_t pos,Int_t val)    {if(pos < 10)fCellNCuts[pos] = val;}
   void SetPIDBitsAt  (Int_t pos,Int_t val)    {if(pos < 10)fPIDBits[pos]   = val;}
   
-  private:
-  Bool_t IsBadRun(Int_t /*iRun*/) const {return kFALSE;} //Tests if this run bad according to private list
+  //MC analysis related methods
+  void FillAcceptanceHistograms();
+  void FillMCVersusRecDataHistograms(const Int_t    index1,  const Int_t    index2,
+                                     const Float_t  pt1,     const Float_t  pt2, 
+                                     const Int_t    ncells1, const Int_t    ncells2, 
+                                     const Double_t mass,    const Double_t pt,  const Double_t asym,    
+                                     const Double_t deta,    const Double_t dphi);
+  
+  void SwitchOnMultipleCutAnalysisInSimulation()   {fMultiCutAnaSim = kTRUE;}
+  void SwitchOffMultipleCutAnalysisInSimulation()  {fMultiCutAnaSim = kFALSE;}
+
   
   private:
-  Bool_t   fDoOwnMix;      // Do combinatorial background not the one provided by the frame
-  Int_t    fNCentrBin ;           // Number of bins in event container for centrality
-  // Int_t    fNrpBin ;            // Number of bins in event container for reaction plain
-  Int_t    fNmaxMixEv ;           // Maximal number of events stored in buffer for mixing
-  TString  fCalorimeter ;  // Select Calorimeter for IM
-  Int_t    fNModules ;     // Number of EMCAL/PHOS modules, set as many histogras as modules 
-  Bool_t   fUseAngleCut ;  // Select pairs depending on their opening angle
-  TList ** fEventsList ;   //[fNCentrBin*GetNZvertBin()*GetNRPBin()]! Containers for photons in stored events
-  Bool_t   fMultiCutAna;   // Do analysis with several or fixed cut
-  Int_t    fNPtCuts;       // number of pt cuts
-  Float_t  fPtCuts[10];    // array with different pt cuts
-  Int_t    fNAsymCuts;     // number of assymmetry cuts
-  Float_t  fAsymCuts[10];  // array with different assymetry cuts
-  Int_t    fNCellNCuts;    // number of cuts with number of cells in cluster
-  Int_t    fCellNCuts[10]; // array with different cell number cluster cuts
-  Int_t    fNPIDBits ;          // number of possible PID bit combinations
-  Int_t    fPIDBits[10];   // array with different PID bits
-
-  Bool_t   fMakeInvPtPlots;// do plots with inverse pt weight
-  Bool_t   fSameSM;        // select only pairs in same SM;
+  Bool_t   fDoOwnMix;            // Do combinatorial background not the one provided by the frame
+  Int_t    fNCentrBin ;                 // Number of bins in event container for centrality
+//Int_t    fNZvertBin ;                 // Number of bins in event container for vertex position
+//Int_t    fNrpBin ;              // Number of bins in event container for reaction plain
+  Int_t    fNmaxMixEv ;                 // Maximal number of events stored in buffer for mixing
+  TString  fCalorimeter ;        // Select Calorimeter for IM
+  Int_t    fNModules ;           // Number of EMCAL/PHOS modules, set as many histogras as modules 
+  Bool_t   fUseAngleCut ;        // Select pairs depending on their opening angle
+  Bool_t   fUseAngleEDepCut ;    // Select pairs depending on their opening angle
+  Float_t  fAngleCut ;           // Select pairs with opening angle larger than a threshold
+  Float_t  fAngleMaxCut ;        // Select pairs with opening angle smaller than a threshold
+  TList ** fEventsList ;         //![fNCentrBin*GetNZvertBin()*GetNRPBin()] Containers for photons in stored events
+  
+  //Multiple cuts analysis
+  Bool_t   fMultiCutAna;         // Do analysis with several or fixed cut
+  Bool_t   fMultiCutAnaSim;      // Do analysis with several or fixed cut, in the simulation related part
+  Int_t    fNPtCuts;             // Number of pt cuts
+  Float_t  fPtCuts[10];          // Array with different pt cuts
+  Int_t    fNAsymCuts;           // Number of assymmetry cuts
+  Float_t  fAsymCuts[10];        // Array with different assymetry cuts
+  Int_t    fNCellNCuts;          // Number of cuts with number of cells in cluster
+  Int_t    fCellNCuts[10];       // Array with different cell number cluster cuts
+  Int_t    fNPIDBits ;                // Number of possible PID bit combinations
+  Int_t    fPIDBits[10];         // Array with different PID bits
+  
+  //Switchs of different analysis options
+  Bool_t   fMakeInvPtPlots;      // D plots with inverse pt weight
+  Bool_t   fSameSM;              // Select only pairs in same SM;
+  Bool_t   fUseTrackMultBins;    // Use track multiplicity and not centrality bins
+  Bool_t   fUsePhotonMultBins;   // Use photon multiplicity and not centrality bins
+  Bool_t   fUseAverClusterEBins; // Use cluster average energy and not centrality bins
+  Bool_t   fUseAverCellEBins;    // Use cell average energy and not centrality bins
+  Bool_t   fFillBadDistHisto;    // Do plots for different distances to bad channels
   
   //Histograms
-  TH2D ** fhReMod ;     //[fNModules] REAL two-photon invariant mass distribution for different calorimeter modules.
-  TH2D ** fhReDiffMod ; //[fNModules+1] REAL two-photon invariant mass distribution for different clusters in different calorimeter modules.
-
-  TH2D ** fhRe1 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry 
-  TH2D ** fhMi1 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry
-  TH2D ** fhRe2 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry 
-  TH2D ** fhMi2 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry
-  TH2D ** fhRe3 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry 
-  TH2D ** fhMi3 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry
-  
-  TH2D ** fhReInvPt1 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT
-  TH2D ** fhMiInvPt1 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT
-  TH2D ** fhReInvPt2 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT 
-  TH2D ** fhMiInvPt2 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT
-  TH2D ** fhReInvPt3 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT
-  TH2D ** fhMiInvPt3 ;  //[fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT
-  
-  //Multiple cuts
-  TH2D ** fhRePtNCellAsymCuts ; //[fNPtCuts*fNAsymCuts*fNCellNCuts] REAL two-photon invariant mass distribution for different pt cut, n cell cuts and assymetry
-  TH2D ** fhRePIDBits ;         //[fNPIDBits] REAL two-photon invariant mass distribution for different PID bits
-  TH3D ** fhRePtMult ;          //[fNAsymCuts] REAL two-photon invariant mass distribution for different track multiplicity and assymetry cuts
-  
-  TH2D *  fhRePtAsym    ;       //!REAL two-photon pt vs asymmetry
-  TH2D *  fhRePtAsymPi0 ;       //!REAL two-photon pt vs asymmetry, close to pi0 mass
-  TH2D *  fhRePtAsymEta ;       //!REAL two-photon pt vs asymmetry, close to eta mass
-
-  TH3D * fhEvents; //!Number of events per centrality, RP, zbin
-
-  TH2D * fhRealOpeningAngle ;    //! Opening angle of pair versus pair energy
-  TH2D * fhRealCosOpeningAngle ; //! Cosinus of opening angle of pair version pair energy
-
-  //Acceptance
-  TH1D * fhPrimPt ;    //! Spectrum of Primary 
-  TH1D * fhPrimAccPt ; //! Spectrum of primary with accepted daughters 
-  TH1D * fhPrimY ;     //! Rapidity distribution of primary particles
-  TH1D * fhPrimAccY ;  //! Rapidity distribution of primary with accepted daughters
-  TH1D * fhPrimPhi ;   //! Azimutal distribution of primary particles
-  TH1D * fhPrimAccPhi; //! Azimutal distribution of primary with accepted daughters    
-  TH2D * fhPrimOpeningAngle ;    //! Opening angle of pair versus pair energy, primaries
-  TH2D * fhPrimCosOpeningAngle ; //! Cosinus of opening angle of pair version pair energy, primaries
-       
-  ClassDef(AliAnaPi0,12)
+  
+  TH1F* fhAverTotECluster;          //! Average number of clusters in SM
+  TH1F* fhAverTotECell;             //! Average number of cells    in SM
+  
+  TH2D ** fhReMod ;                 //![fNModules]   REAL  two-photon invariant mass distribution for different calorimeter modules.
+  TH2D ** fhReDiffMod ;             //![fNModules+1] REAL  two-photon invariant mass distribution for different clusters in different calorimeter modules.
+  TH2D ** fhMiMod ;                 //![fNModules]   MIXED two-photon invariant mass distribution for different calorimeter modules.
+  TH2D ** fhMiDiffMod ;             //![fNModules+1] MIXED two-photon invariant mass distribution for different clusters in different calorimeter modules.
+  
+  // Pairs with at least one cluster tagged as conversion
+  TH2D * fhReConv ;                 //! REAL  two-photon invariant mass distribution one of the pair was 2 clusters with small mass 
+  TH2D * fhMiConv ;                 //! MIXED two-photon invariant mass distribution one of the pair was 2 clusters with small mass
+  TH2D * fhReConv2 ;                //! REAL  two-photon invariant mass distribution both pair photons recombined from 2 clusters with small mass 
+  TH2D * fhMiConv2 ;                //! MIXED two-photon invariant mass distribution both pair photons recombined from 2 clusters with small mass
+
+  TH2D ** fhRe1 ;                   //![fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry 
+  TH2D ** fhMi1 ;                   //![fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry
+  TH2D ** fhRe2 ;                   //![fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry 
+  TH2D ** fhMi2 ;                   //![fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry
+  TH2D ** fhRe3 ;                   //![fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry 
+  TH2D ** fhMi3 ;                   //![fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry
+  
+  //Histograms weighted by inverse pT
+  TH2D ** fhReInvPt1 ;              //![fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT
+  TH2D ** fhMiInvPt1 ;              //![fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT
+  TH2D ** fhReInvPt2 ;              //![fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT 
+  TH2D ** fhMiInvPt2 ;              //![fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT
+  TH2D ** fhReInvPt3 ;              //![fNCentrBin*fNPIDBits*fNAsymCuts] REAL  two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT
+  TH2D ** fhMiInvPt3 ;              //![fNCentrBin*fNPIDBits*fNAsymCuts] MIXED two-photon invariant mass distribution for different centralities and Asymmetry, inverse pT
+  
+  //Multiple cuts: Assymmetry, pt, n cells, PID
+  TH2D ** fhRePtNCellAsymCuts ;     //![fNPtCuts*fNAsymCuts*fNCellNCuts] REAL two-photon invariant mass distribution for different pt cut, n cell cuts and assymetry
+  TH2D ** fhRePtNCellAsymCutsSM0 ;  //![fNPtCuts*fNAsymCuts*fNCellNCuts] REAL two-photon invariant mass distribution for different pt cut, n cell cuts and assymetry
+  TH2D ** fhRePtNCellAsymCutsSM1 ;  //![fNPtCuts*fNAsymCuts*fNCellNCuts] REAL two-photon invariant mass distribution for different pt cut, n cell cuts and assymetry
+  TH2D ** fhRePtNCellAsymCutsSM2 ;  //![fNPtCuts*fNAsymCuts*fNCellNCuts] REAL two-photon invariant mass distribution for different pt cut, n cell cuts and assymetry
+  TH2D ** fhRePtNCellAsymCutsSM3 ;  //![fNPtCuts*fNAsymCuts*fNCellNCuts] REAL two-photon invariant mass distribution for different pt cut, n cell cuts and assymetry
+  TH2D ** fhMiPtNCellAsymCuts ;     //![fNPtCuts*fNAsymCuts*fNCellNCuts] Mixed two-photon invariant mass distribution for different pt cut, n cell cuts and assymetry
+  TH2D ** fhRePIDBits ;             //![fNPIDBits]  REAL two-photon invariant mass distribution for different PID bits
+  TH3D ** fhRePtMult ;              //![fNAsymCuts] REAL two-photon invariant mass distribution for different track multiplicity and assymetry cuts
+  
+  // Asymmetry vs pt, in pi0/eta regions
+  TH2D *  fhRePtAsym    ;           //! REAL two-photon pt vs asymmetry
+  TH2D *  fhRePtAsymPi0 ;           //! REAL two-photon pt vs asymmetry, close to pi0 mass
+  TH2D *  fhRePtAsymEta ;           //! REAL two-photon pt vs asymmetry, close to eta mass
+
+  TH3D * fhEvents;                  //! Number of events per centrality, RP, zbin
+  
+  // Pair opening angle
+  TH2D * fhRealOpeningAngle ;       //! Opening angle of pair versus pair energy
+  TH2D * fhRealCosOpeningAngle ;    //! Cosinus of opening angle of pair version pair energy
+  TH2D * fhMixedOpeningAngle ;      //! Opening angle of pair versus pair energy
+  TH2D * fhMixedCosOpeningAngle ;   //! Cosinus of opening angle of pair version pair energy
+  
+  //MC analysis histograms
+  //Pi0 Acceptance
+  TH1D * fhPrimPi0Pt ;              //! Spectrum of Primary 
+  TH1D * fhPrimPi0AccPt ;           //! Spectrum of primary with accepted daughters 
+  TH1D * fhPrimPi0Y ;               //! Rapidity distribution of primary particles
+  TH1D * fhPrimPi0AccY ;            //! Rapidity distribution of primary with accepted daughters
+  TH1D * fhPrimPi0Phi ;             //! Azimutal distribution of primary particles
+  TH1D * fhPrimPi0AccPhi;           //! Azimutal distribution of primary with accepted daughters       
+  TH2D * fhPrimPi0OpeningAngle ;    //! Opening angle of pair versus pair energy, primaries
+  TH2D * fhPrimPi0CosOpeningAngle ; //! Cosinus of opening angle of pair version pair energy, primaries
+  //Eta acceptance
+  TH1D * fhPrimEtaPt ;              //! Spectrum of Primary 
+  TH1D * fhPrimEtaAccPt ;           //! Spectrum of primary with accepted daughters 
+  TH1D * fhPrimEtaY ;               //! Rapidity distribution of primary particles
+  TH1D * fhPrimEtaAccY ;            //! Rapidity distribution of primary with accepted daughters
+  TH1D * fhPrimEtaPhi ;             //! Azimutal distribution of primary particles
+  TH1D * fhPrimEtaAccPhi;           //! Azimutal distribution of primary with accepted daughters       
+  
+  //Pair origin
+  //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles, 
+  // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
+  TH2D * fhMCOrgMass[13];           //! Mass vs pt of real pairs, check common origin of pair
+  TH2D * fhMCOrgAsym[13];           //! Asymmetry vs pt of real pairs, check common origin of pair
+  TH2D * fhMCOrgDeltaEta[13];       //! Delta Eta vs pt of real pairs, check common origin of pair
+  TH2D * fhMCOrgDeltaPhi[13];       //! Delta Phi vs pt of real pairs, check common origin of pair
+  
+  //Multiple cuts in simulation, origin pi0 or eta
+  TH2D ** fhMCPi0MassPtRec;         //![fNPtCuts*fNAsymCuts*fNCellNCuts] Real pi0 pairs, reconstructed mass vs reconstructed pt of original pair  
+  TH2D ** fhMCPi0MassPtTrue;        //![fNPtCuts*fNAsymCuts*fNCellNCuts] Real pi0 pairs, reconstructed mass vs generated pt of original pair  
+  TH2D ** fhMCPi0PtTruePtRec;       //![fNPtCuts*fNAsymCuts*fNCellNCuts] Real pi0 pairs, reconstructed pt vs generated pt of pair
+  TH2D ** fhMCEtaMassPtRec;         //![fNPtCuts*fNAsymCuts*fNCellNCuts] Real eta pairs, reconstructed mass vs reconstructed pt of original pair  
+  TH2D ** fhMCEtaMassPtTrue;        //![fNPtCuts*fNAsymCuts*fNCellNCuts] Real eta pairs, reconstructed mass vs generated pt of original pair  
+  TH2D ** fhMCEtaPtTruePtRec;       //![fNPtCuts*fNAsymCuts*fNCellNCuts] Real eta pairs, reconstructed pt vs generated pt of pair
+
+  ClassDef(AliAnaPi0,13)
 } ;