PWG4/PartCorrDep/AliAnaPi0.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id: $ */
  16
  17 //_________________________________________________________________________
  18 // Class to collect two-photon invariant mass distributions for
  19 // extracting raw pi0 yield.
  20 // Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles),
  21 // it will do nothing if executed alone
  22 //
  23 //-- Author: Dmitri Peressounko (RRC "KI")
  24 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
  25 //-- and Gustavo Conesa (INFN-Frascati)
  26 //_________________________________________________________________________
  27
  28
  29 // --- ROOT system ---
  30 #include "TH3.h"
  31 #include "TH2F.h"
  32 //#include "Riostream.h"
  33 #include "TCanvas.h"
  34 #include "TPad.h"
  35 #include "TROOT.h"
  36 #include "TClonesArray.h"
  37 #include "TObjString.h"
  38 #include "TDatabasePDG.h"
  39
  40 //---- AliRoot system ----
  41 #include "AliAnaPi0.h"
  42 #include "AliCaloTrackReader.h"
  43 #include "AliCaloPID.h"
  44 #include "AliStack.h"
  45 #include "AliFiducialCut.h"
  46 #include "TParticle.h"
  47 #include "AliVEvent.h"
  48 #include "AliESDCaloCluster.h"
  49 #include "AliESDEvent.h"
  50 #include "AliAODEvent.h"
  51 #include "AliNeutralMesonSelection.h"
  52 #include "AliMixedEvent.h"
  53 #include "AliAODMCParticle.h"
  54
  55 // --- Detectors ---
  56 #include "AliPHOSGeoUtils.h"
  57 #include "AliEMCALGeometry.h"
  58
  59 ClassImp(AliAnaPi0)
  60
  61 //________________________________________________________________________________________________________________________________________________
  62 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
  63 fDoOwnMix(kFALSE),           fEventsList(0x0),
  64 fCalorimeter(""),            fNModules(12),
  65 fUseAngleCut(kFALSE),        fUseAngleEDepCut(kFALSE),     fAngleCut(0),                 fAngleMaxCut(7.),
  66 fMultiCutAna(kFALSE),        fMultiCutAnaSim(kFALSE),
  67 fNPtCuts(0),                 fNAsymCuts(0),                fNCellNCuts(0),               fNPIDBits(0),
  68 fMakeInvPtPlots(kFALSE),     fSameSM(kFALSE),              fFillSMCombinations(kFALSE),  fCheckConversion(kFALSE),
  69 fUseTrackMultBins(kFALSE),   fUsePhotonMultBins(kFALSE),   fUseAverCellEBins(kFALSE),    fUseAverClusterEBins(kFALSE),
  70 fUseAverClusterEDenBins(0),  fFillBadDistHisto(kFALSE),
  71 fhAverTotECluster(0),        fhAverTotECell(0),            fhAverTotECellvsCluster(0),
  72 fhEDensityCluster(0),        fhEDensityCell(0),            fhEDensityCellvsCluster(0),
  73 //Histograms
  74 fhReMod(0x0),                fhReSameSideEMCALMod(0x0),    fhReSameSectorEMCALMod(0x0),  fhReDiffPHOSMod(0x0),
  75 fhMiMod(0x0),                fhMiSameSideEMCALMod(0x0),    fhMiSameSectorEMCALMod(0x0),  fhMiDiffPHOSMod(0x0),
  76 fhReConv(0x0),               fhMiConv(0x0),                fhReConv2(0x0),  fhMiConv2(0x0),
  77 fhRe1(0x0),                  fhMi1(0x0),                   fhRe2(0x0),                   fhMi2(0x0),
  78 fhRe3(0x0),                  fhMi3(0x0),                   fhReInvPt1(0x0),              fhMiInvPt1(0x0),
  79 fhReInvPt2(0x0),             fhMiInvPt2(0x0),              fhReInvPt3(0x0),              fhMiInvPt3(0x0),
  80 fhRePtNCellAsymCuts(0x0),    fhMiPtNCellAsymCuts(0x0),     fhRePtNCellAsymCutsSM(),
  81 fhRePIDBits(0x0),            fhRePtMult(0x0),              fhReSS(),
  82 fhRePtAsym(0x0),             fhRePtAsymPi0(0x0),           fhRePtAsymEta(0x0),
  83 fhEvents(0x0),               fhCentrality(0x0),            fhCentralityNoPair(0x0),
  84 fhEventPlaneAngle(0x0),      fhEventPlaneResolution(0x0),
  85 fhRealOpeningAngle(0x0),     fhRealCosOpeningAngle(0x0),   fhMixedOpeningAngle(0x0),     fhMixedCosOpeningAngle(0x0),
  86 // MC histograms
  87 fhPrimPi0Pt(0x0),            fhPrimPi0AccPt(0x0),          fhPrimPi0Y(0x0),              fhPrimPi0AccY(0x0),
  88 fhPrimPi0Phi(0x0),           fhPrimPi0AccPhi(0x0),         fhPrimPi0OpeningAngle(0x0),   fhPrimPi0CosOpeningAngle(0x0),
  89 fhPrimEtaPt(0x0),            fhPrimEtaAccPt(0x0),          fhPrimEtaY(0x0),              fhPrimEtaAccY(0x0),
  90 fhPrimEtaPhi(0x0),           fhPrimEtaAccPhi(0x0),         fhPrimPi0PtOrigin(0x0),       fhPrimEtaPtOrigin(0x0),
  91 fhMCOrgMass(),               fhMCOrgAsym(),                fhMCOrgDeltaEta(),            fhMCOrgDeltaPhi(),
  92 fhMCPi0MassPtRec(),          fhMCPi0MassPtTrue(),          fhMCPi0PtTruePtRec(),
  93 fhMCEtaMassPtRec(),          fhMCEtaMassPtTrue(),          fhMCEtaPtTruePtRec(),
  94 fhMCPi0PtOrigin(0x0),        fhMCEtaPtOrigin(0x0),
  95 fhReMCFromConversion(0),     fhReMCFromNotConversion(0),   fhReMCFromMixConversion(0)
  96 {
  97 //Default Ctor
  98  InitParameters();
  99
 100 }
 101
 102 //________________________________________________________________________________________________________________________________________________
 103 AliAnaPi0::~AliAnaPi0() {
 104   // Remove event containers
 105
 106   if(fDoOwnMix && fEventsList){
 107     for(Int_t ic=0; ic<GetNCentrBin(); ic++){
 108       for(Int_t iz=0; iz<GetNZvertBin(); iz++){
 109         for(Int_t irp=0; irp<GetNRPBin(); irp++){
 110           fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
 111           delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
 112         }
 113       }
 114     }
 115     delete[] fEventsList;
 116   }
 117
 118 }
 119
 120 //________________________________________________________________________________________________________________________________________________
 121 void AliAnaPi0::InitParameters()
 122 {
 123   //Init parameters when first called the analysis
 124   //Set default parameters
 125   SetInputAODName("PWG4Particle");
 126
 127   AddToHistogramsName("AnaPi0_");
 128   fNModules = 12; // set maximum to maximum number of EMCAL modules
 129
 130   fCalorimeter  = "PHOS";
 131   fUseAngleCut = kFALSE;
 132   fUseAngleEDepCut = kFALSE;
 133   fAngleCut    = 0.;
 134   fAngleMaxCut = TMath::Pi();
 135
 136   fMultiCutAna = kFALSE;
 137
 138   fNPtCuts = 1;
 139   fPtCuts[0] = 0.; fPtCuts[1] = 0.3;   fPtCuts[2] = 0.5;
 140   for(Int_t i = fNPtCuts; i < 10; i++)fPtCuts[i] = 0.;
 141
 142   fNAsymCuts = 2;
 143   fAsymCuts[0] = 1.;  fAsymCuts[1] = 0.7; //fAsymCuts[2] = 0.6; //  fAsymCuts[3] = 0.1;
 144   for(Int_t i = fNAsymCuts; i < 10; i++)fAsymCuts[i] = 0.;
 145
 146   fNCellNCuts = 1;
 147   fCellNCuts[0] = 0; fCellNCuts[1] = 1;   fCellNCuts[2] = 2;
 148   for(Int_t i = fNCellNCuts; i < 10; i++)fCellNCuts[i]  = 0;
 149
 150   fNPIDBits = 1;
 151   fPIDBits[0] = 0;   fPIDBits[1] = 2; //  fPIDBits[2] = 4; fPIDBits[3] = 6;// check, no cut,  dispersion, neutral, dispersion&&neutral
 152   for(Int_t i = fNPIDBits; i < 10; i++)fPIDBits[i] = 0;
 153
 154 }
 155
 156
 157 //________________________________________________________________________________________________________________________________________________
 158 TObjString * AliAnaPi0::GetAnalysisCuts()
 159 {
 160   //Save parameters used for analysis
 161   TString parList ; //this will be list of parameters used for this analysis.
 162   const Int_t buffersize = 255;
 163   char onePar[buffersize] ;
 164   snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
 165   parList+=onePar ;
 166   snprintf(onePar,buffersize,"Number of bins in Centrality:  %d \n",GetNCentrBin()) ;
 167   parList+=onePar ;
 168   snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
 169   parList+=onePar ;
 170   snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
 171   parList+=onePar ;
 172   snprintf(onePar,buffersize,"Depth of event buffer: %d \n",GetNMaxEvMix()) ;
 173   parList+=onePar ;
 174   snprintf(onePar,buffersize,"Pair in same Module: %d ; Fill Different SM histos %d; CheckConversions %d; TrackMult as centrality: %d; PhotonMult as centrality: %d; cluster E as centrality: %d; cell as centrality: %d; Fill InvPt histos %d\n",
 175            fSameSM, fFillSMCombinations, fCheckConversion, fUseTrackMultBins, fUsePhotonMultBins, fUseAverClusterEBins, fUseAverCellEBins, fMakeInvPtPlots) ;
 176   parList+=onePar ;
 177   snprintf(onePar,buffersize,"Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f,\n",fUseAngleCut, fUseAngleEDepCut,fAngleCut,fAngleMaxCut) ;
 178   parList+=onePar ;
 179   snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
 180   for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
 181   parList+=onePar ;
 182   snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
 183   for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
 184   parList+=onePar ;
 185   snprintf(onePar,buffersize,"Cuts: \n") ;
 186   parList+=onePar ;
 187   snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
 188   parList+=onePar ;
 189   snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
 190   parList+=onePar ;
 191   snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
 192   parList+=onePar ;
 193   if(fMultiCutAna){
 194     snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
 195     for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
 196     parList+=onePar ;
 197     snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
 198     for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
 199     parList+=onePar ;
 200   }
 201
 202   return new TObjString(parList) ;
 203 }
 204
 205 //________________________________________________________________________________________________________________________________________________
 206 TList * AliAnaPi0::GetCreateOutputObjects()
 207 {
 208   // Create histograms to be saved in output file and
 209   // store them in fOutputContainer
 210
 211   //create event containers
 212   fEventsList = new TList*[GetNCentrBin()*GetNZvertBin()*GetNRPBin()] ;
 213
 214   for(Int_t ic=0; ic<GetNCentrBin(); ic++){
 215     for(Int_t iz=0; iz<GetNZvertBin(); iz++){
 216       for(Int_t irp=0; irp<GetNRPBin(); irp++){
 217         fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
 218         fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->SetOwner(kFALSE);
 219       }
 220     }
 221   }
 222
 223   TList * outputContainer = new TList() ;
 224   outputContainer->SetName(GetName());
 225
 226   fhReMod                = new TH2F*[fNModules]   ;
 227   fhMiMod                = new TH2F*[fNModules]   ;
 228
 229   if(fCalorimeter == "PHOS"){
 230     fhReDiffPHOSMod        = new TH2F*[fNModules]   ;
 231     fhMiDiffPHOSMod        = new TH2F*[fNModules]   ;
 232   }
 233   else{
 234     fhReSameSectorEMCALMod = new TH2F*[fNModules/2] ;
 235     fhReSameSideEMCALMod   = new TH2F*[fNModules-2] ;
 236     fhMiSameSectorEMCALMod = new TH2F*[fNModules/2] ;
 237     fhMiSameSideEMCALMod   = new TH2F*[fNModules-2] ;
 238   }
 239
 240
 241   fhRe1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 242   fhMi1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 243   if(fFillBadDistHisto){
 244     fhRe2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 245     fhRe3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 246     fhMi2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 247     fhMi3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 248   }
 249   if(fMakeInvPtPlots) {
 250     fhReInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 251     fhMiInvPt1 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 252     if(fFillBadDistHisto){
 253       fhReInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 254       fhReInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 255       fhMiInvPt2 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 256       fhMiInvPt3 = new TH2F*[GetNCentrBin()*fNPIDBits*fNAsymCuts] ;
 257     }
 258   }
 259
 260   const Int_t buffersize = 255;
 261   char key[buffersize] ;
 262   char title[buffersize] ;
 263
 264   Int_t nptbins   = GetHistoPtBins();
 265   Int_t nphibins  = GetHistoPhiBins();
 266   Int_t netabins  = GetHistoEtaBins();
 267   Float_t ptmax   = GetHistoPtMax();
 268   Float_t phimax  = GetHistoPhiMax();
 269   Float_t etamax  = GetHistoEtaMax();
 270   Float_t ptmin   = GetHistoPtMin();
 271   Float_t phimin  = GetHistoPhiMin();
 272   Float_t etamin  = GetHistoEtaMin();
 273
 274   Int_t nmassbins = GetHistoMassBins();
 275   Int_t nasymbins = GetHistoAsymmetryBins();
 276   Float_t massmax = GetHistoMassMax();
 277   Float_t asymmax = GetHistoAsymmetryMax();
 278   Float_t massmin = GetHistoMassMin();
 279   Float_t asymmin = GetHistoAsymmetryMin();
 280   Int_t ntrmbins  = GetHistoTrackMultiplicityBins();
 281   Int_t ntrmmax   = GetHistoTrackMultiplicityMax();
 282   Int_t ntrmmin   = GetHistoTrackMultiplicityMin();
 283
 284   if(GetNCentrBin() > 1 && (fUseAverCellEBins||fUseAverClusterEBins||fUseAverClusterEDenBins)){
 285
 286     fhAverTotECluster = new TH1F("hAverTotECluster","hAverTotECluster",200,0,50) ;
 287     fhAverTotECluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
 288     outputContainer->Add(fhAverTotECluster) ;
 289
 290     fhAverTotECell    = new TH1F("hAverTotECell","hAverTotECell",200,0,50) ;
 291     fhAverTotECell->SetXTitle("E_{cell, aver. SM} (GeV)");
 292     outputContainer->Add(fhAverTotECell) ;
 293
 294     fhAverTotECellvsCluster    = new TH2F("hAverTotECellvsCluster","hAverTotECellvsCluster",200,0,50,200,0,50) ;
 295     fhAverTotECellvsCluster->SetYTitle("E_{cell, aver. SM} (GeV)");
 296     fhAverTotECellvsCluster->SetXTitle("E_{cluster, aver. SM} (GeV)");
 297     outputContainer->Add(fhAverTotECellvsCluster) ;
 298
 299     fhEDensityCluster = new TH1F("hEDensityCluster","hEDensityCluster",200,0,50) ;
 300     fhEDensityCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
 301     outputContainer->Add(fhEDensityCluster) ;
 302
 303     fhEDensityCell    = new TH1F("hEDensityCell","hEDensityCell",200,0,50) ;
 304     fhEDensityCell->SetXTitle("#Sigma E_{cell} / N_{cell} (GeV)");
 305     outputContainer->Add(fhEDensityCell) ;
 306
 307     fhEDensityCellvsCluster    = new TH2F("hEDensityCellvsCluster","hEDensityCellvsCluster",200,0,50,200,0,50) ;
 308     fhEDensityCellvsCluster->SetYTitle("#Sigma E_{cell} / N_{cell} (GeV)");
 309     fhEDensityCellvsCluster->SetXTitle("#Sigma E_{cluster} / N_{cluster} (GeV)");
 310     outputContainer->Add(fhEDensityCellvsCluster) ;
 311
 312   }//counting and average histograms
 313
 314   if(fCheckConversion){
 315     fhReConv = new TH2F("hReConv","Real Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 316     fhReConv->SetXTitle("p_{T} (GeV/c)");
 317     fhReConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 318     outputContainer->Add(fhReConv) ;
 319
 320     fhReConv2 = new TH2F("hReConv2","Real Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 321     fhReConv2->SetXTitle("p_{T} (GeV/c)");
 322     fhReConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 323     outputContainer->Add(fhReConv2) ;
 324
 325     if(fDoOwnMix){
 326       fhMiConv = new TH2F("hMiConv","Mixed Pair with one recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 327       fhMiConv->SetXTitle("p_{T} (GeV/c)");
 328       fhMiConv->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 329       outputContainer->Add(fhMiConv) ;
 330
 331       fhMiConv2 = new TH2F("hMiConv2","Mixed Pair with 2 recombined conversion ",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 332       fhMiConv2->SetXTitle("p_{T} (GeV/c)");
 333       fhMiConv2->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 334       outputContainer->Add(fhMiConv2) ;
 335     }
 336   }
 337
 338   for(Int_t ic=0; ic<GetNCentrBin(); ic++){
 339     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
 340       for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 341         Int_t index = ((ic*fNPIDBits)+ipid)*fNAsymCuts + iasym;
 342         //printf("cen %d, pid %d, asy %d, Index %d\n",ic,ipid,iasym,index);
 343         //Distance to bad module 1
 344         snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
 345         snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
 346                  ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 347         fhRe1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 348         fhRe1[index]->SetXTitle("p_{T} (GeV/c)");
 349         fhRe1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 350         //printf("name: %s\n ",fhRe1[index]->GetName());
 351         outputContainer->Add(fhRe1[index]) ;
 352
 353         if(fFillBadDistHisto){
 354           //Distance to bad module 2
 355           snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
 356           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
 357                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 358           fhRe2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 359           fhRe2[index]->SetXTitle("p_{T} (GeV/c)");
 360           fhRe2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 361           outputContainer->Add(fhRe2[index]) ;
 362
 363           //Distance to bad module 3
 364           snprintf(key, buffersize,"hRe_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
 365           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
 366                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 367           fhRe3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 368           fhRe3[index]->SetXTitle("p_{T} (GeV/c)");
 369           fhRe3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 370           outputContainer->Add(fhRe3[index]) ;
 371         }
 372
 373         //Inverse pT
 374         if(fMakeInvPtPlots){
 375           //Distance to bad module 1
 376           snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
 377           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
 378                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 379           fhReInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 380           fhReInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
 381           fhReInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 382           outputContainer->Add(fhReInvPt1[index]) ;
 383
 384           if(fFillBadDistHisto){
 385             //Distance to bad module 2
 386             snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
 387             snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
 388                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 389             fhReInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 390             fhReInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
 391             fhReInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 392             outputContainer->Add(fhReInvPt2[index]) ;
 393
 394             //Distance to bad module 3
 395             snprintf(key, buffersize,"hReInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
 396             snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
 397                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 398             fhReInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 399             fhReInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
 400             fhReInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 401             outputContainer->Add(fhReInvPt3[index]) ;
 402           }
 403         }
 404         if(fDoOwnMix){
 405           //Distance to bad module 1
 406           snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
 407           snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
 408                    ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 409           fhMi1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 410           fhMi1[index]->SetXTitle("p_{T} (GeV/c)");
 411           fhMi1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 412           outputContainer->Add(fhMi1[index]) ;
 413           if(fFillBadDistHisto){
 414             //Distance to bad module 2
 415             snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
 416             snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
 417                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 418             fhMi2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 419             fhMi2[index]->SetXTitle("p_{T} (GeV/c)");
 420             fhMi2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 421             outputContainer->Add(fhMi2[index]) ;
 422
 423             //Distance to bad module 3
 424             snprintf(key, buffersize,"hMi_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
 425             snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 3",
 426                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 427             fhMi3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 428             fhMi3[index]->SetXTitle("p_{T} (GeV/c)");
 429             fhMi3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 430             outputContainer->Add(fhMi3[index]) ;
 431           }
 432           //Inverse pT
 433           if(fMakeInvPtPlots){
 434             //Distance to bad module 1
 435             snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist1",ic,ipid,iasym) ;
 436             snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 1",
 437                      ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 438             fhMiInvPt1[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 439             fhMiInvPt1[index]->SetXTitle("p_{T} (GeV/c)");
 440             fhMiInvPt1[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 441             outputContainer->Add(fhMiInvPt1[index]) ;
 442             if(fFillBadDistHisto){
 443               //Distance to bad module 2
 444               snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist2",ic,ipid,iasym) ;
 445               snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f, dist bad 2",
 446                        ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 447               fhMiInvPt2[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 448               fhMiInvPt2[index]->SetXTitle("p_{T} (GeV/c)");
 449               fhMiInvPt2[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 450               outputContainer->Add(fhMiInvPt2[index]) ;
 451
 452               //Distance to bad module 3
 453               snprintf(key, buffersize,"hMiInvPt_cen%d_pidbit%d_asy%d_dist3",ic,ipid,iasym) ;
 454               snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d, PID bit=%d and asymmetry %1.2f,dist bad 3",
 455                        ic,fPIDBits[ipid], fAsymCuts[iasym]) ;
 456               fhMiInvPt3[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 457               fhMiInvPt3[index]->SetXTitle("p_{T} (GeV/c)");
 458               fhMiInvPt3[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 459               outputContainer->Add(fhMiInvPt3[index]) ;
 460             }
 461           }
 462         }
 463       }
 464     }
 465   }
 466
 467   fhRePtAsym = new TH2F("hRePtAsym","Asymmetry vs pt, for pairs",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
 468   fhRePtAsym->SetXTitle("p_{T} (GeV/c)");
 469   fhRePtAsym->SetYTitle("Asymmetry");
 470   outputContainer->Add(fhRePtAsym);
 471
 472   fhRePtAsymPi0 = new TH2F("hRePtAsymPi0","Asymmetry vs pt, for pairs close to #pi^{0} mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
 473   fhRePtAsymPi0->SetXTitle("p_{T} (GeV/c)");
 474   fhRePtAsymPi0->SetYTitle("Asymmetry");
 475   outputContainer->Add(fhRePtAsymPi0);
 476
 477   fhRePtAsymEta = new TH2F("hRePtAsymEta","Asymmetry vs pt, for pairs close to #eta mass",nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
 478   fhRePtAsymEta->SetXTitle("p_{T} (GeV/c)");
 479   fhRePtAsymEta->SetYTitle("Asymmetry");
 480   outputContainer->Add(fhRePtAsymEta);
 481
 482   if(fMultiCutAna){
 483
 484     fhRePIDBits         = new TH2F*[fNPIDBits];
 485     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
 486       snprintf(key,   buffersize,"hRe_pidbit%d",ipid) ;
 487       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
 488       fhRePIDBits[ipid] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 489       fhRePIDBits[ipid]->SetXTitle("p_{T} (GeV/c)");
 490       fhRePIDBits[ipid]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 491       outputContainer->Add(fhRePIDBits[ipid]) ;
 492     }// pid bit loop
 493
 494     fhRePtNCellAsymCuts    = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 495     fhMiPtNCellAsymCuts    = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 496
 497     if(fFillSMCombinations){
 498       for(Int_t iSM = 0; iSM < fNModules; iSM++) fhRePtNCellAsymCutsSM[iSM] = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 499
 500     }
 501
 502     for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
 503       for(Int_t icell=0; icell<fNCellNCuts; icell++){
 504         for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 505           snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
 506           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
 507           Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
 508           //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
 509           fhRePtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 510           fhRePtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
 511           fhRePtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 512           outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
 513
 514           snprintf(key,   buffersize,"hMi_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
 515           snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
 516           fhMiPtNCellAsymCuts[index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 517           fhMiPtNCellAsymCuts[index]->SetXTitle("p_{T} (GeV/c)");
 518           fhMiPtNCellAsymCuts[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 519           outputContainer->Add(fhMiPtNCellAsymCuts[index]) ;
 520
 521           if(fFillSMCombinations){
 522             for(Int_t iSM = 0; iSM < fNModules; iSM++){
 523               snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d_SM%d",ipt,icell,iasym,iSM) ;
 524               snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%1.2f, SM %d ",
 525                        fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym],iSM) ;
 526               fhRePtNCellAsymCutsSM[iSM][index] = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 527               fhRePtNCellAsymCutsSM[iSM][index]->SetXTitle("p_{T} (GeV/c)");
 528               fhRePtNCellAsymCutsSM[iSM][index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 529               outputContainer->Add(fhRePtNCellAsymCutsSM[iSM][index]) ;
 530             }
 531
 532           }
 533         }
 534       }
 535     }
 536
 537     if(ntrmbins!=0){
 538       fhRePtMult = new TH3F*[fNAsymCuts] ;
 539       for(Int_t iasym = 0; iasym<fNAsymCuts; iasym++){
 540         fhRePtMult[iasym] = new TH3F(Form("hRePtMult_asym%d",iasym),Form("(p_{T},C,M)_{#gamma#gamma}, A<%1.2f",fAsymCuts[iasym]),
 541                                      nptbins,ptmin,ptmax,ntrmbins,ntrmmin,ntrmmax,nmassbins,massmin,massmax);
 542         fhRePtMult[iasym]->SetXTitle("p_{T} (GeV/c)");
 543         fhRePtMult[iasym]->SetYTitle("Track multiplicity");
 544         fhRePtMult[iasym]->SetZTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 545         outputContainer->Add(fhRePtMult[iasym]) ;
 546       }
 547     }
 548   }// multi cuts analysis
 549
 550   fhReSS[0] = new TH2F("hRe_SS_Tight"," 0.01 < #lambda_{0}^{2} < 0.4",
 551                        nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 552   fhReSS[0]->SetXTitle("p_{T} (GeV/c)");
 553   fhReSS[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 554   outputContainer->Add(fhReSS[0]) ;
 555
 556
 557   fhReSS[1] = new TH2F("hRe_SS_Loose"," #lambda_{0}^{2} > 0.4",
 558                        nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 559   fhReSS[1]->SetXTitle("p_{T} (GeV/c)");
 560   fhReSS[1]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 561   outputContainer->Add(fhReSS[1]) ;
 562
 563
 564   fhReSS[2] = new TH2F("hRe_SS_Both"," cluster_{1} #lambda_{0}^{2} > 0.4; cluster_{2} 0.01 < #lambda_{0}^{2} < 0.4",
 565                        nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 566   fhReSS[2]->SetXTitle("p_{T} (GeV/c)");
 567   fhReSS[2]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 568   outputContainer->Add(fhReSS[2]) ;
 569
 570   fhEvents=new TH3F("hEvents","Number of events",GetNCentrBin(),0.,1.*GetNCentrBin(),
 571                     GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
 572
 573   fhEvents->SetXTitle("Centrality bin");
 574   fhEvents->SetYTitle("Z vertex bin bin");
 575   fhEvents->SetZTitle("RP bin");
 576   outputContainer->Add(fhEvents) ;
 577
 578   if(GetNCentrBin()>1){
 579     fhCentrality=new TH1F("hCentralityBin","Number of events in centrality bin",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
 580     fhCentrality->SetXTitle("Centrality bin");
 581     outputContainer->Add(fhCentrality) ;
 582
 583     fhCentralityNoPair=new TH1F("hCentralityBinNoPair","Number of events in centrality bin, with no cluster pairs",GetNCentrBin(),0.,1.*GetNCentrBin()) ;
 584     fhCentralityNoPair->SetXTitle("Centrality bin");
 585     outputContainer->Add(fhCentralityNoPair) ;
 586   }
 587
 588   if(GetNRPBin() > 1 ){
 589
 590     fhEventPlaneAngle=new TH1F("hEventPlaneAngleBin","Number of events in centrality bin",100,0.,TMath::TwoPi()) ;
 591     fhEventPlaneAngle->SetXTitle("EP angle (rad)");
 592     outputContainer->Add(fhEventPlaneAngle) ;
 593
 594     if(GetNCentrBin()>1){
 595       fhEventPlaneResolution=new TH2F("hEventPlaneResolution","Event plane resolution",GetNCentrBin(),0,GetNCentrBin(),100,0.,TMath::TwoPi()) ;
 596       fhEventPlaneResolution->SetYTitle("Resolution");
 597       fhEventPlaneResolution->SetXTitle("Centrality Bin");
 598       outputContainer->Add(fhEventPlaneResolution) ;
 599     }
 600   }
 601
 602   fhRealOpeningAngle  = new TH2F
 603   ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,300,0,TMath::Pi());
 604   fhRealOpeningAngle->SetYTitle("#theta(rad)");
 605   fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 606   outputContainer->Add(fhRealOpeningAngle) ;
 607
 608   fhRealCosOpeningAngle  = new TH2F
 609   ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,1);
 610   fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
 611   fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 612   outputContainer->Add(fhRealCosOpeningAngle) ;
 613
 614   if(fDoOwnMix){
 615
 616     fhMixedOpeningAngle  = new TH2F
 617     ("hMixedOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,300,0,TMath::Pi());
 618     fhMixedOpeningAngle->SetYTitle("#theta(rad)");
 619     fhMixedOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 620     outputContainer->Add(fhMixedOpeningAngle) ;
 621
 622     fhMixedCosOpeningAngle  = new TH2F
 623     ("hMixedCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}, Mixed pairs",nptbins,ptmin,ptmax,100,0,1);
 624     fhMixedCosOpeningAngle->SetYTitle("cos (#theta) ");
 625     fhMixedCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 626     outputContainer->Add(fhMixedCosOpeningAngle) ;
 627
 628   }
 629
 630   //Histograms filled only if MC data is requested
 631   if(IsDataMC()){
 632
 633     fhReMCFromConversion = new TH2F("hReMCFromConversion","Invariant mass of 2 clusters originated in conversions",
 634                          nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 635     fhReMCFromConversion->SetXTitle("p_{T} (GeV/c)");
 636     fhReMCFromConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 637     outputContainer->Add(fhReMCFromConversion) ;
 638
 639     fhReMCFromNotConversion = new TH2F("hReMCNotFromConversion","Invariant mass of 2 clusters not originated in conversions",
 640                                     nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 641     fhReMCFromNotConversion->SetXTitle("p_{T} (GeV/c)");
 642     fhReMCFromNotConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 643     outputContainer->Add(fhReMCFromNotConversion) ;
 644
 645     fhReMCFromMixConversion = new TH2F("hReMCFromMixConversion","Invariant mass of 2 clusters one from conversion and the other not",
 646                                     nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 647     fhReMCFromMixConversion->SetXTitle("p_{T} (GeV/c)");
 648     fhReMCFromMixConversion->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 649     outputContainer->Add(fhReMCFromMixConversion) ;
 650
 651     //Pi0
 652     fhPrimPi0Pt     = new TH1F("hPrimPi0Pt","Primary pi0 pt, Y<1",nptbins,ptmin,ptmax) ;
 653     fhPrimPi0AccPt  = new TH1F("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
 654     fhPrimPi0Pt   ->SetXTitle("p_{T} (GeV/c)");
 655     fhPrimPi0AccPt->SetXTitle("p_{T} (GeV/c)");
 656     outputContainer->Add(fhPrimPi0Pt) ;
 657     outputContainer->Add(fhPrimPi0AccPt) ;
 658
 659     Int_t netabinsopen =  TMath::Nint(netabins*4/(etamax-etamin));
 660     fhPrimPi0Y      = new TH2F("hPrimPi0Rapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
 661     fhPrimPi0Y   ->SetYTitle("Rapidity");
 662     fhPrimPi0Y   ->SetXTitle("p_{T} (GeV/c)");
 663     outputContainer->Add(fhPrimPi0Y) ;
 664
 665     fhPrimPi0AccY   = new TH2F("hPrimPi0AccRapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
 666     fhPrimPi0AccY->SetYTitle("Rapidity");
 667     fhPrimPi0AccY->SetXTitle("p_{T} (GeV/c)");
 668     outputContainer->Add(fhPrimPi0AccY) ;
 669
 670     Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
 671     fhPrimPi0Phi    = new TH2F("hPrimPi0Phi","Azimuthal of primary pi0, Y<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
 672     fhPrimPi0Phi->SetYTitle("#phi (deg)");
 673     fhPrimPi0Phi->SetXTitle("p_{T} (GeV/c)");
 674     outputContainer->Add(fhPrimPi0Phi) ;
 675
 676     fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","Azimuthal of primary pi0 with accepted daughters",nptbins,ptmin,ptmax,
 677                                nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 678     fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
 679     fhPrimPi0AccPhi->SetXTitle("p_{T} (GeV/c)");
 680     outputContainer->Add(fhPrimPi0AccPhi) ;
 681
 682     //Eta
 683     fhPrimEtaPt     = new TH1F("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
 684     fhPrimEtaAccPt  = new TH1F("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
 685     fhPrimEtaPt   ->SetXTitle("p_{T} (GeV/c)");
 686     fhPrimEtaAccPt->SetXTitle("p_{T} (GeV/c)");
 687     outputContainer->Add(fhPrimEtaPt) ;
 688     outputContainer->Add(fhPrimEtaAccPt) ;
 689
 690     fhPrimEtaY      = new TH2F("hPrimEtaRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
 691     fhPrimEtaY->SetYTitle("Rapidity");
 692     fhPrimEtaY->SetXTitle("p_{T} (GeV/c)");
 693     outputContainer->Add(fhPrimEtaY) ;
 694
 695     fhPrimEtaAccY   = new TH2F("hPrimEtaAccRapidity","Rapidity of primary eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
 696     fhPrimEtaAccY->SetYTitle("Rapidity");
 697     fhPrimEtaAccY->SetXTitle("p_{T} (GeV/c)");
 698     outputContainer->Add(fhPrimEtaAccY) ;
 699
 700     fhPrimEtaPhi    = new TH2F("hPrimEtaPhi","Azimuthal of primary eta",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 701     fhPrimEtaPhi->SetYTitle("#phi (deg)");
 702     fhPrimEtaPhi->SetXTitle("p_{T} (GeV/c)");
 703     outputContainer->Add(fhPrimEtaPhi) ;
 704
 705     fhPrimEtaAccPhi = new TH2F("hPrimEtaAccPhi","Azimuthal of primary eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 706     fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
 707     fhPrimEtaAccPhi->SetXTitle("p_{T} (GeV/c)");
 708     outputContainer->Add(fhPrimEtaAccPhi) ;
 709
 710
 711     //Prim origin
 712     //Pi0
 713     fhPrimPi0PtOrigin     = new TH2F("hPrimPi0PtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
 714     fhPrimPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
 715     fhPrimPi0PtOrigin->SetYTitle("Origin");
 716     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 717     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 718     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
 719     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 720     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
 721     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
 722     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
 723     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
 724     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
 725     fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
 726     outputContainer->Add(fhPrimPi0PtOrigin) ;
 727
 728     fhMCPi0PtOrigin     = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,11,0,11) ;
 729     fhMCPi0PtOrigin->SetXTitle("p_{T} (GeV/c)");
 730     fhMCPi0PtOrigin->SetYTitle("Origin");
 731     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 732     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 733     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
 734     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 735     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
 736     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
 737     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
 738     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
 739     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
 740     fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
 741     outputContainer->Add(fhMCPi0PtOrigin) ;
 742
 743     //Eta
 744     fhPrimEtaPtOrigin     = new TH2F("hPrimEtaPtOrigin","Primary pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
 745     fhPrimEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
 746     fhPrimEtaPtOrigin->SetYTitle("Origin");
 747     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 748     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 749     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
 750     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 751     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
 752     fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
 753
 754     outputContainer->Add(fhPrimEtaPtOrigin) ;
 755
 756     fhMCEtaPtOrigin     = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated pi0 pt vs origin",nptbins,ptmin,ptmax,7,0,7) ;
 757     fhMCEtaPtOrigin->SetXTitle("p_{T} (GeV/c)");
 758     fhMCEtaPtOrigin->SetYTitle("Origin");
 759     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 760     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 761     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
 762     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 763     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
 764     fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
 765
 766     outputContainer->Add(fhMCEtaPtOrigin) ;
 767
 768
 769     fhPrimPi0OpeningAngle  = new TH2F
 770     ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
 771     fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
 772     fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 773     outputContainer->Add(fhPrimPi0OpeningAngle) ;
 774
 775     fhPrimPi0CosOpeningAngle  = new TH2F
 776     ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
 777     fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
 778     fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 779     outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
 780
 781     for(Int_t i = 0; i<13; i++){
 782       fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("mass vs pt, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 783       fhMCOrgMass[i]->SetXTitle("p_{T} (GeV/c)");
 784       fhMCOrgMass[i]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 785       outputContainer->Add(fhMCOrgMass[i]) ;
 786
 787       fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("asymmetry vs pt, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
 788       fhMCOrgAsym[i]->SetXTitle("p_{T} (GeV/c)");
 789       fhMCOrgAsym[i]->SetYTitle("A");
 790       outputContainer->Add(fhMCOrgAsym[i]) ;
 791
 792       fhMCOrgDeltaEta[i] = new TH2F(Form("hMCOrgDeltaEta_%d",i),Form("#Delta #eta of pair vs pt, origin %d",i),nptbins,ptmin,ptmax,netabins,-1.4,1.4) ;
 793       fhMCOrgDeltaEta[i]->SetXTitle("p_{T} (GeV/c)");
 794       fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
 795       outputContainer->Add(fhMCOrgDeltaEta[i]) ;
 796
 797       fhMCOrgDeltaPhi[i]= new TH2F(Form("hMCOrgDeltaPhi_%d",i),Form("#Delta #phi of pair vs p_{T}, origin %d",i),nptbins,ptmin,ptmax,nphibins,-0.7,0.7) ;
 798       fhMCOrgDeltaPhi[i]->SetXTitle("p_{T} (GeV/c)");
 799       fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
 800       outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
 801
 802     }
 803
 804     if(fMultiCutAnaSim){
 805       fhMCPi0MassPtTrue  = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 806       fhMCPi0MassPtRec   = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 807       fhMCPi0PtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 808       fhMCEtaMassPtRec   = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 809       fhMCEtaMassPtTrue  = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 810       fhMCEtaPtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 811       for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
 812         for(Int_t icell=0; icell<fNCellNCuts; icell++){
 813           for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 814             Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
 815
 816             fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
 817                                                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]),
 818                                                nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 819             fhMCPi0MassPtRec[index]->SetXTitle("p_{T, reconstructed} (GeV/c)");
 820             fhMCPi0MassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 821             outputContainer->Add(fhMCPi0MassPtRec[index]) ;
 822
 823             fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
 824                                                 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]),
 825                                                 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 826             fhMCPi0MassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
 827             fhMCPi0MassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 828             outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
 829
 830             fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
 831                                                  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]),
 832                                                  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
 833             fhMCPi0PtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
 834             fhMCPi0PtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
 835             outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
 836
 837             fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
 838                                                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]),
 839                                                nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 840             fhMCEtaMassPtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
 841             fhMCEtaMassPtRec[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 842             outputContainer->Add(fhMCEtaMassPtRec[index]) ;
 843
 844             fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
 845                                                 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]),
 846                                                 nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 847             fhMCEtaMassPtTrue[index]->SetXTitle("p_{T, generated} (GeV/c)");
 848             fhMCEtaMassPtTrue[index]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 849             outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
 850
 851             fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
 852                                                  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]),
 853                                                  nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
 854             fhMCEtaPtTruePtRec[index]->SetXTitle("p_{T, generated} (GeV/c)");
 855             fhMCEtaPtTruePtRec[index]->SetYTitle("p_{T, reconstructed} (GeV/c)");
 856             outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
 857           }
 858         }
 859       }
 860     }//multi cut ana
 861     else {
 862       fhMCPi0MassPtTrue  = new TH2F*[1];
 863       fhMCPi0PtTruePtRec = new TH2F*[1];
 864       fhMCEtaMassPtTrue  = new TH2F*[1];
 865       fhMCEtaPtTruePtRec = new TH2F*[1];
 866
 867       fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 868       fhMCPi0MassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
 869       fhMCPi0MassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 870       outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
 871
 872       fhMCPi0PtTruePtRec[0]= new TH2F("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) ;
 873       fhMCPi0PtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
 874       fhMCPi0PtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
 875       outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
 876
 877       fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 878       fhMCEtaMassPtTrue[0]->SetXTitle("p_{T, generated} (GeV/c)");
 879       fhMCEtaMassPtTrue[0]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 880       outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
 881
 882       fhMCEtaPtTruePtRec[0]= new TH2F("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) ;
 883       fhMCEtaPtTruePtRec[0]->SetXTitle("p_{T, generated} (GeV/c)");
 884       fhMCEtaPtTruePtRec[0]->SetYTitle("p_{T, reconstructed} (GeV/c)");
 885       outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
 886     }
 887   }
 888
 889   if(fFillSMCombinations){
 890     TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
 891     for(Int_t imod=0; imod<fNModules; imod++){
 892       //Module dependent invariant mass
 893       snprintf(key, buffersize,"hReMod_%d",imod) ;
 894       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
 895       fhReMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 896       fhReMod[imod]->SetXTitle("p_{T} (GeV/c)");
 897       fhReMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 898       outputContainer->Add(fhReMod[imod]) ;
 899       if(fCalorimeter=="PHOS"){
 900         snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
 901         snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
 902         fhReDiffPHOSMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 903         fhReDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
 904         fhReDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 905         outputContainer->Add(fhReDiffPHOSMod[imod]) ;
 906       }
 907       else{//EMCAL
 908         if(imod<fNModules/2){
 909           snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
 910           snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
 911           fhReSameSectorEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 912           fhReSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
 913           fhReSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 914           outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
 915         }
 916         if(imod<fNModules-2){
 917           snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
 918           snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
 919           fhReSameSideEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 920           fhReSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
 921           fhReSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 922           outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
 923         }
 924       }//EMCAL
 925
 926       if(fDoOwnMix){
 927         snprintf(key, buffersize,"hMiMod_%d",imod) ;
 928         snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for Module %d",imod) ;
 929         fhMiMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 930         fhMiMod[imod]->SetXTitle("p_{T} (GeV/c)");
 931         fhMiMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 932         outputContainer->Add(fhMiMod[imod]) ;
 933
 934         if(fCalorimeter=="PHOS"){
 935           snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
 936           snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
 937           fhMiDiffPHOSMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 938           fhMiDiffPHOSMod[imod]->SetXTitle("p_{T} (GeV/c)");
 939           fhMiDiffPHOSMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 940           outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
 941         }//PHOS
 942         else{//EMCAL
 943           if(imod<fNModules/2){
 944             snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
 945             snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
 946             fhMiSameSectorEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 947             fhMiSameSectorEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
 948             fhMiSameSectorEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 949             outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
 950           }
 951           if(imod<fNModules-2){
 952             snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
 953             snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
 954             fhMiSameSideEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 955             fhMiSameSideEMCALMod[imod]->SetXTitle("p_{T} (GeV/c)");
 956             fhMiSameSideEMCALMod[imod]->SetYTitle("m_{#gamma,#gamma} (GeV/c^{2})");
 957             outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
 958           }
 959         }//EMCAL
 960       }// own mix
 961     }//loop combinations
 962   } // SM combinations
 963
 964   //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
 965   //
 966   //    printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
 967   //
 968   //  }
 969
 970   return outputContainer;
 971 }
 972
 973 //_________________________________________________________________________________________________________________________________________________
 974 void AliAnaPi0::Print(const Option_t * /*opt*/) const
 975 {
 976   //Print some relevant parameters set for the analysis
 977   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 978   AliAnaPartCorrBaseClass::Print(" ");
 979
 980   printf("Number of bins in Centrality:  %d \n",GetNCentrBin()) ;
 981   printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
 982   printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
 983   printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
 984   printf("Pair in same Module: %d \n",fSameSM) ;
 985   printf("Cuts: \n") ;
 986   // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
 987   printf("Number of modules:             %d \n",fNModules) ;
 988   printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
 989   printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
 990   printf("\tasymmetry < ");
 991   for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
 992   printf("\n");
 993
 994   printf("PID selection bits: n = %d, \n",fNPIDBits) ;
 995   printf("\tPID bit = ");
 996   for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
 997   printf("\n");
 998
 999   if(fMultiCutAna){
1000     printf("pT cuts: n = %d, \n",fNPtCuts) ;
1001     printf("\tpT > ");
1002     for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
1003     printf("GeV/c\n");
1004
1005     printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
1006     printf("\tnCell > ");
1007     for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
1008     printf("\n");
1009
1010   }
1011   printf("------------------------------------------------------\n") ;
1012 }
1013
1014 //_____________________________________________________________
1015 void AliAnaPi0::FillAcceptanceHistograms(){
1016   //Fill acceptance histograms if MC data is available
1017
1018   if(GetReader()->ReadStack()){
1019     AliStack * stack = GetMCStack();
1020     if(stack){
1021       for(Int_t i=0 ; i<stack->GetNtrack(); i++){
1022         TParticle * prim = stack->Particle(i) ;
1023         Int_t pdg = prim->GetPdgCode();
1024         //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
1025         //                             prim->GetName(), prim->GetPdgCode());
1026
1027         if( pdg == 111 || pdg == 221){
1028           Double_t pi0Pt = prim->Pt() ;
1029           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
1030           Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
1031           Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
1032           if(pdg == 111){
1033             if(TMath::Abs(pi0Y) < 1.0){
1034               fhPrimPi0Pt ->Fill(pi0Pt) ;
1035               fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1036             }
1037             fhPrimPi0Y  ->Fill(pi0Pt, pi0Y) ;
1038           }
1039           else if(pdg == 221){
1040             if(TMath::Abs(pi0Y) < 1.0){
1041               fhPrimEtaPt ->Fill(pi0Pt) ;
1042               fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1043             }
1044             fhPrimEtaY  ->Fill(pi0Pt, pi0Y) ;
1045           }
1046
1047           //Origin of meson
1048           Int_t momindex  = prim->GetFirstMother();
1049           if(momindex >= 0) {
1050             TParticle* mother = stack->Particle(momindex);
1051             Int_t mompdg    = TMath::Abs(mother->GetPdgCode());
1052             Int_t momstatus = mother->GetStatusCode();
1053             if(pdg == 111){
1054               if     (momstatus  == 21)fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1055               else if(mompdg     < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1056               else if(mompdg     > 2100  && mompdg   < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1057               else if(mompdg    == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1058               else if(mompdg    == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1059               else if(mompdg    == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1060               else if(mompdg    == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1061               else if(mompdg    >= 310   && mompdg    <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1062               else if(mompdg    == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1063               else if(momstatus == 11 || momstatus  == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1064               else                      fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1065             }//pi0
1066             else {
1067               if     (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1068               else if(mompdg    < 22  ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1069               else if(mompdg    > 2100  && mompdg   < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1070               else if(mompdg    == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1071               else if(momstatus == 11 || momstatus  == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1072               else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1073               //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1074             }
1075           } // pi0 has mother
1076
1077           //Check if both photons hit Calorimeter
1078           if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
1079           Int_t iphot1=prim->GetFirstDaughter() ;
1080           Int_t iphot2=prim->GetLastDaughter() ;
1081           if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
1082             TParticle * phot1 = stack->Particle(iphot1) ;
1083             TParticle * phot2 = stack->Particle(iphot2) ;
1084             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
1085               //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",
1086               //        phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
1087
1088               TLorentzVector lv1, lv2;
1089               phot1->Momentum(lv1);
1090               phot2->Momentum(lv2);
1091               Bool_t inacceptance = kFALSE;
1092               if(fCalorimeter == "PHOS"){
1093                 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
1094                   Int_t mod ;
1095                   Double_t x,z ;
1096                   if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
1097                     inacceptance = kTRUE;
1098                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1099                 }
1100                 else{
1101
1102                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1103                     inacceptance = kTRUE ;
1104                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1105                 }
1106
1107               }
1108               else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
1109                 if(GetEMCALGeometry()){
1110
1111                   Int_t absID1=0;
1112                   Int_t absID2=0;
1113
1114                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1115                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1116
1117                   if( absID1 >= 0 && absID2 >= 0)
1118                     inacceptance = kTRUE;
1119
1120                   //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1121                   //                    inacceptance = kTRUE;
1122                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1123                 }
1124                 else{
1125                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1126                     inacceptance = kTRUE ;
1127                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1128                 }
1129               }
1130
1131               if(inacceptance){
1132                 if(pdg==111){
1133                   fhPrimPi0AccPt ->Fill(pi0Pt) ;
1134                   fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1135                   fhPrimPi0AccY  ->Fill(pi0Pt, pi0Y) ;
1136                   Double_t angle  = lv1.Angle(lv2.Vect());
1137                   fhPrimPi0OpeningAngle   ->Fill(pi0Pt,angle);
1138                   fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1139                 }
1140                 else if(pdg==221){
1141                   fhPrimEtaAccPt ->Fill(pi0Pt) ;
1142                   fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1143                   fhPrimEtaAccY  ->Fill(pi0Pt, pi0Y) ;
1144                 }
1145               }//Accepted
1146             }// 2 photons
1147           }//Check daughters exist
1148         }// Primary pi0 or eta
1149       }//loop on primaries
1150     }//stack exists and data is MC
1151   }//read stack
1152   else if(GetReader()->ReadAODMCParticles()){
1153     TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1154     if(mcparticles){
1155       Int_t nprim = mcparticles->GetEntriesFast();
1156
1157       for(Int_t i=0; i < nprim; i++)
1158       {
1159         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
1160
1161         // Only generator particles, when they come from PYTHIA, PHOJET, HERWIG ...
1162         //if( prim->GetStatus() == 0 && (GetMCAnalysisUtils()->GetMCGenerator()).Length()!=0) break;
1163
1164         Int_t pdg = prim->GetPdgCode();
1165         if( pdg == 111 || pdg == 221){
1166           Double_t pi0Pt = prim->Pt() ;
1167           //printf("pi0, pt %2.2f, eta %f, phi %f\n",pi0Pt, prim->Eta(), prim->Phi());
1168           if(prim->E() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
1169
1170           Double_t pi0Y  = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
1171           Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
1172           if(pdg == 111){
1173             if(TMath::Abs(pi0Y) < 1){
1174               fhPrimPi0Pt->Fill(pi0Pt) ;
1175               fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1176             }
1177             fhPrimPi0Y  ->Fill(pi0Pt, pi0Y) ;
1178           }
1179           else if(pdg == 221){
1180             if(TMath::Abs(pi0Y) < 1){
1181               fhPrimEtaPt->Fill(pi0Pt) ;
1182               fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1183             }
1184             fhPrimEtaY  ->Fill(pi0Pt, pi0Y) ;
1185           }
1186
1187           //Origin of meson
1188           Int_t momindex  = prim->GetMother();
1189           if(momindex >= 0) {
1190             AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1191             Int_t mompdg    = TMath::Abs(mother->GetPdgCode());
1192             Int_t momstatus = mother->GetStatus();
1193             if(pdg == 111){
1194               if     (momstatus  == 21) fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1195               else if(mompdg     < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1196               else if(mompdg     > 2100  && mompdg   < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1197               else if(mompdg    == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1198               else if(mompdg    == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1199               else if(mompdg    == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1200               else if(mompdg    == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1201               else if(mompdg    >= 310   && mompdg    <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1202               else if(mompdg    == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1203               else if(momstatus == 11 || momstatus  == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1204               else                      fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1205             }//pi0
1206             else {
1207               if     (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1208               else if(mompdg    < 22  ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1209               else if(mompdg    > 2100  && mompdg   < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1210               else if(mompdg    == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1211               else if(momstatus == 11 || momstatus  == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1212               else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1213               //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1214             }
1215           }//pi0 has mother
1216
1217           //Check if both photons hit Calorimeter
1218           if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
1219           Int_t iphot1=prim->GetDaughter(0) ;
1220           Int_t iphot2=prim->GetDaughter(1) ;
1221           if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim){
1222             AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
1223             AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
1224             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
1225               TLorentzVector lv1, lv2;
1226               lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1227               lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1228
1229               Bool_t inacceptance = kFALSE;
1230               if(fCalorimeter == "PHOS"){
1231                 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
1232                   Int_t mod ;
1233                   Double_t x,z ;
1234                   Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
1235                   Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
1236                   if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) &&
1237                      GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x))
1238                     inacceptance = kTRUE;
1239                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1240                 }
1241                 else{
1242
1243                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1244                     inacceptance = kTRUE ;
1245                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1246                 }
1247
1248               }
1249               else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
1250                 if(GetEMCALGeometry()){
1251
1252                   Int_t absID1=0;
1253                   Int_t absID2=0;
1254
1255                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1256                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1257
1258                   if( absID1 >= 0 && absID2 >= 0)
1259                     inacceptance = kTRUE;
1260
1261
1262                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1263                 }
1264                 else{
1265                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1266                     inacceptance = kTRUE ;
1267                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1268                 }
1269               }
1270
1271               if(inacceptance){
1272                 if(pdg==111){
1273                   //                printf("ACCEPTED pi0: pt %2.2f, phi %3.2f, eta %1.2f\n",pi0Pt,phi,pi0Y);
1274                   fhPrimPi0AccPt ->Fill(pi0Pt) ;
1275                   fhPrimPi0AccPhi->Fill(pi0Pt, phi) ;
1276                   fhPrimPi0AccY  ->Fill(pi0Pt, pi0Y) ;
1277                   Double_t angle  = lv1.Angle(lv2.Vect());
1278                   fhPrimPi0OpeningAngle   ->Fill(pi0Pt,angle);
1279                   fhPrimPi0CosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
1280                 }
1281                 else if(pdg==221){
1282                   fhPrimEtaAccPt ->Fill(pi0Pt) ;
1283                   fhPrimEtaAccPhi->Fill(pi0Pt, phi) ;
1284                   fhPrimEtaAccY  ->Fill(pi0Pt, pi0Y) ;
1285                 }
1286               }//Accepted
1287             }// 2 photons
1288           }//Check daughters exist
1289         }// Primary pi0 or eta
1290       }//loop on primaries
1291     }//stack exists and data is MC
1292
1293
1294   }     // read AOD MC
1295 }
1296
1297 //_____________________________________________________________
1298 void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1,  const Int_t index2,
1299                                               const Float_t pt1,   const Float_t pt2,
1300                                               const Int_t ncell1,  const Int_t ncell2,
1301                                               const Double_t mass, const Double_t pt, const Double_t asym,
1302                                               const Double_t deta, const Double_t dphi){
1303   //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1304   //Adjusted for Pythia, need to see what to do for other generators.
1305   //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1306   // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1307
1308   Int_t ancPDG    = 0;
1309   Int_t ancStatus = 0;
1310   TLorentzVector ancMomentum;
1311   TVector3 prodVertex;
1312   Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1313                                                               GetReader(), ancPDG, ancStatus,ancMomentum, prodVertex);
1314
1315   Int_t momindex  = -1;
1316   Int_t mompdg    = -1;
1317   Int_t momstatus = -1;
1318   if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1319                             ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1320
1321   if(ancLabel > -1){
1322     if(ancPDG==22){//gamma
1323       fhMCOrgMass[0]->Fill(pt,mass);
1324       fhMCOrgAsym[0]->Fill(pt,asym);
1325       fhMCOrgDeltaEta[0]->Fill(pt,deta);
1326       fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1327     }
1328     else if(TMath::Abs(ancPDG)==11){//e
1329       fhMCOrgMass[1]->Fill(pt,mass);
1330       fhMCOrgAsym[1]->Fill(pt,asym);
1331       fhMCOrgDeltaEta[1]->Fill(pt,deta);
1332       fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1333     }
1334     else if(ancPDG==111){//Pi0
1335       fhMCOrgMass[2]->Fill(pt,mass);
1336       fhMCOrgAsym[2]->Fill(pt,asym);
1337       fhMCOrgDeltaEta[2]->Fill(pt,deta);
1338       fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1339       if(fMultiCutAnaSim){
1340         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1341           for(Int_t icell=0; icell<fNCellNCuts; icell++){
1342             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1343               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1344               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1345                  asym   <  fAsymCuts[iasym]                                   &&
1346                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1347                 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1348                 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1349                 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1350               }//pass the different cuts
1351             }// pid bit cut loop
1352           }// icell loop
1353         }// pt cut loop
1354       }//Multi cut ana sim
1355       else {
1356         fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1357         if(mass < 0.17 && mass > 0.1) {
1358           fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1359
1360           if(GetReader()->ReadStack()){
1361             TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1362             momindex  = ancestor->GetFirstMother();
1363             if(momindex < 0) return;
1364             TParticle* mother = GetMCStack()->Particle(momindex);
1365             mompdg    = TMath::Abs(mother->GetPdgCode());
1366             momstatus = mother->GetStatusCode();
1367           }
1368           else {
1369             TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1370             AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1371             momindex  = ancestor->GetMother();
1372             if(momindex < 0) return;
1373             AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1374             mompdg    = TMath::Abs(mother->GetPdgCode());
1375             momstatus = mother->GetStatus();
1376           }
1377
1378           if     (momstatus  == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1379           else if(mompdg     < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1380           else if(mompdg     > 2100  && mompdg   < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1381           else if(mompdg    == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1382           else if(mompdg    == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1383           else if(mompdg    == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1384           else if(mompdg    == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1385           else if(mompdg    >= 310   && mompdg    <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1386           else if(mompdg    == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1387           else if(momstatus == 11 || momstatus  == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1388           else                      fhMCPi0PtOrigin->Fill(pt,7.5);//other?
1389
1390         }//pi0 mass region
1391
1392       }
1393     }
1394     else if(ancPDG==221){//Eta
1395       fhMCOrgMass[3]->Fill(pt,mass);
1396       fhMCOrgAsym[3]->Fill(pt,asym);
1397       fhMCOrgDeltaEta[3]->Fill(pt,deta);
1398       fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1399       if(fMultiCutAnaSim){
1400         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1401           for(Int_t icell=0; icell<fNCellNCuts; icell++){
1402             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1403               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1404               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1405                  asym   <  fAsymCuts[iasym]                                   &&
1406                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1407                 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1408                 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1409                 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1410               }//pass the different cuts
1411             }// pid bit cut loop
1412           }// icell loop
1413         }// pt cut loop
1414       } //Multi cut ana sim
1415       else {
1416         fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1417         if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1418
1419         if(GetReader()->ReadStack()){
1420           TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1421           momindex  = ancestor->GetFirstMother();
1422           if(momindex < 0) return;
1423           TParticle* mother = GetMCStack()->Particle(momindex);
1424           mompdg    = TMath::Abs(mother->GetPdgCode());
1425           momstatus = mother->GetStatusCode();
1426         }
1427         else {
1428           TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
1429           AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1430           momindex  = ancestor->GetMother();
1431           if(momindex < 0) return;
1432           AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1433           mompdg    = TMath::Abs(mother->GetPdgCode());
1434           momstatus = mother->GetStatus();
1435         }
1436
1437         if     (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1438         else if(mompdg    < 22  ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1439         else if(mompdg    > 2100  && mompdg   < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1440         else if(mompdg    == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1441         else if(momstatus == 11 || momstatus  == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1442         else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1443         //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1444       }// eta mass region
1445     }
1446     else if(ancPDG==-2212){//AProton
1447       fhMCOrgMass[4]->Fill(pt,mass);
1448       fhMCOrgAsym[4]->Fill(pt,asym);
1449       fhMCOrgDeltaEta[4]->Fill(pt,deta);
1450       fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1451     }
1452     else if(ancPDG==-2112){//ANeutron
1453       fhMCOrgMass[5]->Fill(pt,mass);
1454       fhMCOrgAsym[5]->Fill(pt,asym);
1455       fhMCOrgDeltaEta[5]->Fill(pt,deta);
1456       fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1457     }
1458     else if(TMath::Abs(ancPDG)==13){//muons
1459       fhMCOrgMass[6]->Fill(pt,mass);
1460       fhMCOrgAsym[6]->Fill(pt,asym);
1461       fhMCOrgDeltaEta[6]->Fill(pt,deta);
1462       fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1463     }
1464     else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1465       if(ancStatus==1){//Stable particles, converted? not decayed resonances
1466         fhMCOrgMass[6]->Fill(pt,mass);
1467         fhMCOrgAsym[6]->Fill(pt,asym);
1468         fhMCOrgDeltaEta[6]->Fill(pt,deta);
1469         fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1470       }
1471       else{//resonances and other decays, more hadron conversions?
1472         fhMCOrgMass[7]->Fill(pt,mass);
1473         fhMCOrgAsym[7]->Fill(pt,asym);
1474         fhMCOrgDeltaEta[7]->Fill(pt,deta);
1475         fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
1476       }
1477     }
1478     else {//Partons, colliding protons, strings, intermediate corrections
1479       if(ancStatus==11 || ancStatus==12){//String fragmentation
1480         fhMCOrgMass[8]->Fill(pt,mass);
1481         fhMCOrgAsym[8]->Fill(pt,asym);
1482         fhMCOrgDeltaEta[8]->Fill(pt,deta);
1483         fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1484       }
1485       else if (ancStatus==21){
1486         if(ancLabel < 2) {//Colliding protons
1487           fhMCOrgMass[11]->Fill(pt,mass);
1488           fhMCOrgAsym[11]->Fill(pt,asym);
1489           fhMCOrgDeltaEta[11]->Fill(pt,deta);
1490           fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1491         }//colliding protons
1492         else if(ancLabel < 6){//partonic initial states interactions
1493           fhMCOrgMass[9]->Fill(pt,mass);
1494           fhMCOrgAsym[9]->Fill(pt,asym);
1495           fhMCOrgDeltaEta[9]->Fill(pt,deta);
1496           fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1497         }
1498         else if(ancLabel < 8){//Final state partons radiations?
1499           fhMCOrgMass[10]->Fill(pt,mass);
1500           fhMCOrgAsym[10]->Fill(pt,asym);
1501           fhMCOrgDeltaEta[10]->Fill(pt,deta);
1502           fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1503         }
1504         // else {
1505         //   printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1506         //          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1507         // }
1508       }//status 21
1509       //else {
1510       //  printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1511       //         ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1512       // }
1513     }////Partons, colliding protons, strings, intermediate corrections
1514   }//ancLabel > -1
1515   else { //ancLabel <= -1
1516     //printf("Not related at all label = %d\n",ancLabel);
1517     fhMCOrgMass[12]->Fill(pt,mass);
1518     fhMCOrgAsym[12]->Fill(pt,asym);
1519     fhMCOrgDeltaEta[12]->Fill(pt,deta);
1520     fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
1521   }
1522 }
1523
1524 //____________________________________________________________________________________________________________________________________________________
1525 void AliAnaPi0::CountAndGetAverages(Int_t &nClus,Int_t &nCell, Float_t &eClusTot,Float_t &eCellTot, Float_t &eDenClus,Float_t &eDenCell) {
1526   // Count the number of clusters and cells, deposited energy, and do some averages in case multiplicity bins dependent on such numbers
1527   // are requested
1528   if(fCalorimeter=="EMCAL"){
1529     nClus = GetEMCALClusters()  ->GetEntriesFast();
1530     nCell = GetEMCALCells()->GetNumberOfCells();
1531     for(Int_t icl=0; icl < nClus; icl++) {
1532       Float_t e1 = ((AliVCluster*)GetEMCALClusters()->At(icl))->E();
1533       eClusTot +=  e1;
1534     }// first cluster
1535
1536     for(Int_t jce=0; jce < nCell; jce++) eCellTot +=  GetEMCALCells()->GetAmplitude(jce);
1537   }
1538   else {
1539     nClus = GetPHOSClusters()->GetEntriesFast();
1540     nCell = GetPHOSCells()   ->GetNumberOfCells();
1541     for(Int_t icl=0; icl < nClus; icl++) {
1542       Float_t e1 = ((AliVCluster*)GetPHOSClusters()->At(icl))->E();
1543       eClusTot +=  e1;
1544     }// first cluster
1545     for(Int_t jce=0; jce < nCell; jce++) eCellTot +=  GetPHOSCells()->GetAmplitude(jce);
1546   }
1547   if(GetDebug() > 1)
1548     printf("AliAnaPi0::MakeAnalysisFillHistograms() - # Clusters %d, sum cluster E per SM %f,# Cells %d, sum cell E per SM %f\n", nClus,eClusTot,nCell,eCellTot);
1549
1550   //Fill histograms with "energy density", ncell and nclust will be > 0 since there are at least 2 "photons"
1551   eDenClus = eClusTot/nClus;
1552   eDenCell = eCellTot/nCell;
1553   fhEDensityCluster      ->Fill(eDenClus);
1554   fhEDensityCell         ->Fill(eDenCell);
1555   fhEDensityCellvsCluster->Fill(eDenClus, eDenCell);
1556   //Fill the average number of cells or clusters per SM
1557   eClusTot /=fNModules;
1558   eCellTot /=fNModules;
1559   fhAverTotECluster      ->Fill(eClusTot);
1560   fhAverTotECell         ->Fill(eCellTot);
1561   fhAverTotECellvsCluster->Fill(eClusTot, eCellTot);
1562   //printf("Average Cluster: E %f, density %f;  Average Cell E %f, density  %f\n ",eClusTot,eDenClus,eCellTot,eDenCell);
1563 }
1564
1565 //____________________________________________________________________________________________________________________________________________________
1566 void AliAnaPi0::MakeAnalysisFillHistograms()
1567 {
1568   //Process one event and extract photons from AOD branch
1569   // filled with AliAnaPhoton and fill histos with invariant mass
1570
1571   //In case of simulated data, fill acceptance histograms
1572   if(IsDataMC())FillAcceptanceHistograms();
1573
1574   //if (GetReader()->GetEventNumber()%10000 == 0)
1575   // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1576
1577   //Init some variables
1578   Int_t   nPhot    = GetInputAODBranch()->GetEntriesFast() ;
1579   Int_t   nClus    = 0;
1580   Int_t   nCell    = 0;
1581   Float_t eClusTot = 0;
1582   Float_t eCellTot = 0;
1583   Float_t eDenClus = 0;
1584   Float_t eDenCell = 0;
1585
1586   if(GetNCentrBin() > 1 && (fUseAverCellEBins||fUseAverClusterEBins||fUseAverClusterEDenBins))
1587     CountAndGetAverages(nClus,nCell,eClusTot,eCellTot,eDenClus,eDenCell);
1588
1589
1590   if(GetDebug() > 1)
1591     printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
1592
1593   //If less than photon 2 entries in the list, skip this event
1594   if(nPhot < 2 ) {
1595
1596     if(GetDebug() > 2)
1597       printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
1598
1599     if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
1600
1601     return ;
1602   }
1603
1604   //Init variables
1605   Int_t module1         = -1;
1606   Int_t module2         = -1;
1607   Double_t vert[]       = {0.0, 0.0, 0.0} ; //vertex
1608   Int_t evtIndex1       = 0 ;
1609   Int_t currentEvtIndex = -1;
1610   Int_t curCentrBin     = 0 ;
1611   Int_t curRPBin        = 0 ;
1612   Int_t curZvertBin     = 0 ;
1613
1614   //Get shower shape information of clusters
1615   TObjArray *clusters = 0;
1616   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1617   else if(fCalorimeter=="PHOS" ) clusters = GetPHOSClusters() ;
1618
1619   //---------------------------------
1620   //First loop on photons/clusters
1621   //---------------------------------
1622   for(Int_t i1=0; i1<nPhot-1; i1++){
1623     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
1624     //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
1625
1626     // get the event index in the mixed buffer where the photon comes from
1627     // in case of mixing with analysis frame, not own mixing
1628     evtIndex1 = GetEventIndex(p1, vert) ;
1629     //printf("charge = %d\n", track->Charge());
1630     if ( evtIndex1 == -1 )
1631       return ;
1632     if ( evtIndex1 == -2 )
1633       continue ;
1634
1635     //printf("z vertex %f < %f\n",vert[2],GetZvertexCut());
1636     if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ;   //vertex cut
1637
1638
1639     //----------------------------------------------------------------------------
1640     // Get the multiplicity bin. Different cases: centrality (PbPb),
1641     // average cluster multiplicity, average cell multiplicity, track multiplicity
1642     // default is centrality bins
1643     //----------------------------------------------------------------------------
1644     if (evtIndex1 != currentEvtIndex) {
1645       if(fUseTrackMultBins){ // Track multiplicity bins
1646         //printf("track  mult %d\n",GetTrackMultiplicity());
1647         curCentrBin = (GetTrackMultiplicity()-1)/5;
1648         if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
1649         //printf("track mult bin %d\n",curCentrBin);
1650       }
1651       else if(fUsePhotonMultBins){ // Photon multiplicity bins
1652         //printf("photon  mult %d cluster mult %d\n",nPhot, nClus);
1653         curCentrBin = nPhot-2;
1654         if(curCentrBin > GetNCentrBin() -1) curCentrBin=GetNCentrBin()-1;
1655         //printf("photon mult bin %d\n",curRPBin);
1656       }
1657       else if(fUseAverClusterEBins){ // Cluster average energy bins
1658         //Bins for pp, if needed can be done in a more general way
1659         curCentrBin = (Int_t) eClusTot/10 * GetNCentrBin();
1660         if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
1661         //printf("cluster E average %f, bin %d \n",eClusTot,curCentrBin);
1662       }
1663       else if(fUseAverCellEBins){ // Cell average energy bins
1664         //Bins for pp, if needed can be done in a more general way
1665         curCentrBin = (Int_t) eCellTot/10*GetNCentrBin();
1666         if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
1667         //printf("cell E average %f, bin %d \n",eCellTot,curCentrBin);
1668       }
1669       else if(fUseAverClusterEDenBins){ // Energy density bins
1670         //Bins for pp, if needed can be done in a more general way
1671         curCentrBin = (Int_t) eDenClus/10*GetNCentrBin();
1672         if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
1673         //printf("cluster Eden average %f, bin %d \n",eDenClus,curCentrBin);
1674       }
1675       else { //Event centrality
1676         // Centrality task returns at maximum 10, 20 or 100, depending on option chosen and
1677         // number of bins, the bin has to be corrected
1678         curCentrBin = GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt();
1679         if(GetDebug() > 0 )printf("AliAnaPi0::MakeAnalysisFillHistograms() - curCentrBin %d, centrality %d, n bins %d, max bin from centrality %d\n",
1680                                   curCentrBin, GetEventCentrality(), GetNCentrBin(), GetReader()->GetCentralityOpt());
1681       }
1682
1683       if (curCentrBin < 0 || curCentrBin >= GetNCentrBin()){
1684         if(GetDebug() > 0)
1685           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality bin <%d> not expected, n bins <%d> , return\n",curCentrBin,GetNCentrBin());
1686         return;
1687       }
1688
1689       //Reaction plane bin
1690       curRPBin    = 0 ;
1691       if(GetNRPBin()>1 && GetEventPlane()){
1692         Float_t epAngle = GetEventPlane()->GetEventplane(GetEventPlaneMethod());
1693         fhEventPlaneAngle->Fill(epAngle);
1694         curRPBin = TMath::Nint(epAngle*(GetNRPBin()-1)/TMath::Pi());
1695         if(curRPBin >= GetNRPBin()) printf("RP Bin %d out of range %d\n",curRPBin,GetNRPBin());
1696         //printf("RP: %d, %f, angle %f, n bin %d\n", curRPBin,epAngle*(GetNRPBin()-1)/TMath::Pi(),epAngle,GetNRPBin());
1697       }
1698
1699       //Get vertex z bin
1700       curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
1701
1702       //Fill event bin info
1703       fhEvents    ->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
1704       if(GetNCentrBin() > 1) {
1705         fhCentrality->Fill(curCentrBin);
1706         if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
1707       }
1708       currentEvtIndex = evtIndex1 ;
1709       if(GetDebug() > 1)
1710         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Centrality %d, Vertex Bin %d, RP bin %d \n",curCentrBin,curRPBin,curZvertBin);
1711     }
1712
1713     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
1714
1715     //Get the momentum of this cluster
1716     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
1717
1718     //Get (Super)Module number of this cluster
1719     module1 = GetModuleNumber(p1);
1720
1721     //------------------------------------------
1722     //Get index in VCaloCluster array
1723     AliVCluster *cluster1 = 0;
1724     Bool_t bFound1        = kFALSE;
1725     Int_t  caloLabel1     = p1->GetCaloLabel(0);
1726     Bool_t iclus1         =-1;
1727     if(clusters){
1728       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
1729         AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
1730         if(cluster){
1731           if     (cluster->GetID()==caloLabel1) {
1732             bFound1  = kTRUE  ;
1733             cluster1 = cluster;
1734             iclus1   = iclus;
1735           }
1736         }
1737         if(bFound1) break;
1738       }
1739     }// calorimeter clusters loop
1740
1741     //---------------------------------
1742     //Second loop on photons/clusters
1743     //---------------------------------
1744     for(Int_t i2=i1+1; i2<nPhot; i2++){
1745       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
1746
1747       //In case of mixing frame, check we are not in the same event as the first cluster
1748       Int_t evtIndex2 = GetEventIndex(p2, vert) ;
1749       if ( evtIndex2 == -1 )
1750         return ;
1751       if ( evtIndex2 == -2 )
1752         continue ;
1753       if (GetMixedEvent() && (evtIndex1 == evtIndex2))
1754         continue ;
1755
1756       //------------------------------------------
1757       //Get index in VCaloCluster array
1758       AliVCluster *cluster2 = 0;
1759       Bool_t bFound2        = kFALSE;
1760       Int_t caloLabel2      = p2->GetCaloLabel(0);
1761       if(clusters){
1762         for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
1763           AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
1764           if(cluster){
1765             if(cluster->GetID()==caloLabel2) {
1766               bFound2  = kTRUE  ;
1767               cluster2 = cluster;
1768             }
1769           }
1770           if(bFound2) break;
1771         }// calorimeter clusters loop
1772       }
1773
1774       Float_t tof1  = -1;
1775       Float_t l01   = -1;
1776       if(cluster1 && bFound1){
1777         tof1  = cluster1->GetTOF()*1e9;
1778         l01   = cluster1->GetM02();
1779       }
1780       //      else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
1781       //                   p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
1782
1783       Float_t tof2  = -1;
1784       Float_t l02   = -1;
1785       if(cluster2 && bFound2){
1786         tof2  = cluster2->GetTOF()*1e9;
1787         l02   = cluster2->GetM02();
1788
1789       }
1790       //      else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
1791       //                  p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
1792
1793       if(clusters){
1794         Double_t t12diff = tof1-tof2;
1795         if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
1796       }
1797       //------------------------------------------
1798
1799       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
1800
1801       //Get the momentum of this cluster
1802       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
1803       //Get module number
1804       module2       = GetModuleNumber(p2);
1805
1806       //---------------------------------
1807       // Get pair kinematics
1808       //---------------------------------
1809       Double_t m    = (photon1 + photon2).M() ;
1810       Double_t pt   = (photon1 + photon2).Pt();
1811       Double_t deta = photon1.Eta() - photon2.Eta();
1812       Double_t dphi = photon1.Phi() - photon2.Phi();
1813       Double_t a    = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
1814
1815       if(GetDebug() > 2)
1816         printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
1817
1818       //--------------------------------
1819       // Opening angle selection
1820       //--------------------------------
1821       //Check if opening angle is too large or too small compared to what is expected
1822       Double_t angle   = photon1.Angle(photon2.Vect());
1823       if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
1824         if(GetDebug() > 2)
1825           printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
1826         continue;
1827       }
1828
1829       if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
1830         if(GetDebug() > 2)
1831           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
1832         continue;
1833       }
1834
1835       //-------------------------------------------------------------------------------------------------
1836       //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
1837       //-------------------------------------------------------------------------------------------------
1838       if(a < fAsymCuts[0] && fFillSMCombinations){
1839         if(module1==module2 && module1 >=0 && module1<fNModules)
1840           fhReMod[module1]->Fill(pt,m) ;
1841
1842         if(fCalorimeter=="EMCAL"){
1843
1844           // Same sector
1845           Int_t j=0;
1846           for(Int_t i = 0; i < fNModules/2; i++){
1847             j=2*i;
1848             if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
1849           }
1850
1851           // Same side
1852           for(Int_t i = 0; i < fNModules-2; i++){
1853             if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
1854           }
1855         }//EMCAL
1856         else {//PHOS
1857           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
1858           if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
1859           if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
1860         }//PHOS
1861       }
1862
1863       //In case we want only pairs in same (super) module, check their origin.
1864       Bool_t ok = kTRUE;
1865       if(fSameSM && module1!=module2) ok=kFALSE;
1866       if(ok){
1867
1868         //Check if one of the clusters comes from a conversion
1869         if(fCheckConversion){
1870           if     (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
1871           else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
1872         }
1873
1874         // Fill shower shape cut histograms
1875         if     ( l01 > 0.01 && l01 < 0.4  &&
1876                  l02 > 0.01 && l02 < 0.4 )               fhReSS[0]->Fill(pt,m); // Tight
1877         else if( l01 > 0.4  && l02 > 0.4 )               fhReSS[1]->Fill(pt,m); // Loose
1878         else if( l01 > 0.01 && l01 < 0.4  && l02 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
1879         else if( l02 > 0.01 && l02 < 0.4  && l01 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
1880
1881         //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
1882         for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1883           if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
1884             for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
1885               if(a < fAsymCuts[iasym]){
1886                 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
1887                 //printf("index %d :(cen %d * nPID %d + ipid %d)*nasym %d + iasym %d\n",index,curCentrBin,fNPIDBits,ipid,fNAsymCuts,iasym);
1888                 fhRe1     [index]->Fill(pt,m);
1889                 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
1890                 if(fFillBadDistHisto){
1891                   if(p1->DistToBad()>0 && p2->DistToBad()>0){
1892                     fhRe2     [index]->Fill(pt,m) ;
1893                     if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
1894                     if(p1->DistToBad()>1 && p2->DistToBad()>1){
1895                       fhRe3     [index]->Fill(pt,m) ;
1896                       if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
1897                     }// bad 3
1898                   }// bad2
1899                 }// Fill bad dist histos
1900               }//assymetry cut
1901             }// asymmetry cut loop
1902           }// bad 1
1903         }// pid bit loop
1904
1905         //Fill histograms with opening angle
1906         fhRealOpeningAngle   ->Fill(pt,angle);
1907         fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
1908
1909         //Fill histograms with pair assymmetry
1910         fhRePtAsym->Fill(pt,a);
1911         if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
1912         if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
1913
1914         //-------------------------------------------------------
1915         //Get the number of cells needed for multi cut analysis.
1916         //-------------------------------------------------------
1917         Int_t ncell1 = 0;
1918         Int_t ncell2 = 0;
1919         if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim)){
1920
1921           AliVEvent * event = GetReader()->GetInputEvent();
1922           if(event){
1923             for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
1924               AliVCluster *cluster = event->GetCaloCluster(iclus);
1925
1926               Bool_t is = kFALSE;
1927               if     (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
1928               else if(fCalorimeter == "PHOS"  && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
1929
1930               if(is){
1931                 if      (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
1932                 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
1933               } // PHOS or EMCAL cluster as requested in analysis
1934
1935               if(ncell2 > 0 &&  ncell1 > 0) break; // No need to continue the iteration
1936
1937             }
1938             //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
1939           }
1940         }
1941
1942         //---------
1943         // MC data
1944         //---------
1945         //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1946         if(IsDataMC()) {
1947
1948           if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
1949              GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
1950           {
1951             fhReMCFromConversion->Fill(pt,m);
1952           }
1953           else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
1954                   !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
1955           {
1956             fhReMCFromNotConversion->Fill(pt,m);
1957           }
1958           else
1959           {
1960             fhReMCFromMixConversion->Fill(pt,m);
1961           }
1962
1963           FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
1964         }
1965
1966         //-----------------------
1967         //Multi cuts analysis
1968         //-----------------------
1969         if(fMultiCutAna){
1970           //Histograms for different PID bits selection
1971           for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
1972
1973             if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    &&
1974                p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
1975
1976             //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
1977           } // pid bit cut loop
1978
1979           //Several pt,ncell and asymmetry cuts
1980           for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1981             for(Int_t icell=0; icell<fNCellNCuts; icell++){
1982               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1983                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1984                 if(p1->E() >   fPtCuts[ipt]      && p2->E() > fPtCuts[ipt]        &&
1985                    a        <   fAsymCuts[iasym]                                    &&
1986                    ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell]){
1987                   fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
1988                   //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
1989                   if(fFillSMCombinations && module1==module2){
1990                     fhRePtNCellAsymCutsSM[module1][index]->Fill(pt,m) ;
1991                   }
1992                 }
1993               }// pid bit cut loop
1994             }// icell loop
1995           }// pt cut loop
1996           if(GetHistoTrackMultiplicityBins()){
1997             for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
1998               if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
1999             }
2000           }
2001         }// multiple cuts analysis
2002       }// ok if same sm
2003     }// second same event particle
2004   }// first cluster
2005
2006   //-------------------------------------------------------------
2007   // Mixing
2008   //-------------------------------------------------------------
2009   if(fDoOwnMix){
2010     //printf("Cen bin %d, RP bin %d, e aver %f, mult %d\n",curCentrBin,curRPBin, eClusTot, nPhot);
2011     //Recover events in with same characteristics as the current event
2012     TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
2013     Int_t nMixed = evMixList->GetSize() ;
2014     for(Int_t ii=0; ii<nMixed; ii++){
2015       TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2016       Int_t nPhot2=ev2->GetEntriesFast() ;
2017       Double_t m = -999;
2018       if(GetDebug() > 1)
2019         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, curCentrBin);
2020
2021       //---------------------------------
2022       //First loop on photons/clusters
2023       //---------------------------------
2024       for(Int_t i1=0; i1<nPhot; i1++){
2025         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2026         if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2027
2028         //Get kinematics of cluster and (super) module of this cluster
2029         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2030         module1 = GetModuleNumber(p1);
2031
2032         //---------------------------------
2033         //First loop on photons/clusters
2034         //---------------------------------
2035         for(Int_t i2=0; i2<nPhot2; i2++){
2036           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2037
2038           //Get kinematics of second cluster and calculate those of the pair
2039           TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2040           m           = (photon1+photon2).M() ;
2041           Double_t pt = (photon1 + photon2).Pt();
2042           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2043
2044           //Check if opening angle is too large or too small compared to what is expected
2045           Double_t angle   = photon1.Angle(photon2.Vect());
2046           if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
2047             if(GetDebug() > 2)
2048               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2049             continue;
2050           }
2051           if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2052             if(GetDebug() > 2)
2053               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2054             continue;
2055
2056           }
2057
2058           if(GetDebug() > 2)
2059             printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
2060                    p1->Pt(), p2->Pt(), pt,m,a);
2061
2062           //In case we want only pairs in same (super) module, check their origin.
2063           module2 = GetModuleNumber(p2);
2064
2065           //-------------------------------------------------------------------------------------------------
2066           //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2067           //-------------------------------------------------------------------------------------------------
2068           if(a < fAsymCuts[0] && fFillSMCombinations){
2069             if(module1==module2 && module1 >=0 && module1<fNModules)
2070               fhMiMod[module1]->Fill(pt,m) ;
2071
2072             if(fCalorimeter=="EMCAL"){
2073
2074               // Same sector
2075               Int_t j=0;
2076               for(Int_t i = 0; i < fNModules/2; i++){
2077                 j=2*i;
2078                 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
2079               }
2080
2081               // Same side
2082               for(Int_t i = 0; i < fNModules-2; i++){
2083                 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
2084               }
2085             }//EMCAL
2086             else {//PHOS
2087               if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2088               if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2089               if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
2090             }//PHOS
2091
2092
2093           }
2094
2095           Bool_t ok = kTRUE;
2096           if(fSameSM && module1!=module2) ok=kFALSE;
2097           if(ok){
2098
2099             //Check if one of the clusters comes from a conversion
2100             if(fCheckConversion){
2101               if     (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2102               else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2103             }
2104             //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2105             for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2106               if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
2107                 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
2108                   if(a < fAsymCuts[iasym]){
2109                     Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2110                     fhMi1     [index]->Fill(pt,m) ;
2111                     if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
2112                     if(fFillBadDistHisto){
2113                       if(p1->DistToBad()>0 && p2->DistToBad()>0){
2114                         fhMi2     [index]->Fill(pt,m) ;
2115                         if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
2116                         if(p1->DistToBad()>1 && p2->DistToBad()>1){
2117                           fhMi3     [index]->Fill(pt,m) ;
2118                           if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
2119                         }
2120                       }
2121                     }// Fill bad dist histo
2122                   }//Asymmetry cut
2123                 }// Asymmetry loop
2124               }//PID cut
2125             }//loop for histograms
2126
2127             //-----------------------
2128             //Multi cuts analysis
2129             //-----------------------
2130             if(fMultiCutAna){
2131               //Several pt,ncell and asymmetry cuts
2132
2133               for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2134                 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2135                   for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2136                     Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2137                     if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        &&
2138                        a        <   fAsymCuts[iasym]                                    //&&
2139                        //p1->GetBtag() >=  fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2140                        ){
2141                       fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
2142                       //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2143                     }
2144                   }// pid bit cut loop
2145                 }// icell loop
2146               }// pt cut loop
2147             } // Multi cut ana
2148
2149             //Fill histograms with opening angle
2150             fhMixedOpeningAngle   ->Fill(pt,angle);
2151             fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2152           }//ok
2153         }// second cluster loop
2154       }//first cluster loop
2155     }//loop on mixed events
2156
2157     //--------------------------------------------------------
2158     //Add the current event to the list of events for mixing
2159     //--------------------------------------------------------
2160     TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2161     //Add current event to buffer and Remove redundant events
2162     if(currentEvent->GetEntriesFast()>0){
2163       evMixList->AddFirst(currentEvent) ;
2164       currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2165       if(evMixList->GetSize() >= GetNMaxEvMix())
2166       {
2167         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2168         evMixList->RemoveLast() ;
2169         delete tmp ;
2170       }
2171     }
2172     else{ //empty event
2173       delete currentEvent ;
2174       currentEvent=0 ;
2175     }
2176   }// DoOwnMix
2177
2178 }
2179
2180 //____________________________________________________________________________________________________________________________________________________
2181 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2182 {
2183   // retieves the event index and checks the vertex
2184   //    in the mixed buffer returns -2 if vertex NOK
2185   //    for normal events   returns 0 if vertex OK and -1 if vertex NOK
2186
2187   Int_t evtIndex = -1 ;
2188   if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2189
2190     if (GetMixedEvent()){
2191
2192       evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2193       GetVertex(vert,evtIndex);
2194
2195       if(TMath::Abs(vert[2])> GetZvertexCut())
2196         evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2197     } else {// Single event
2198
2199       GetVertex(vert);
2200
2201       if(TMath::Abs(vert[2])> GetZvertexCut())
2202         evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2203       else
2204         evtIndex = 0 ;
2205     }
2206   }//No MC reader
2207   else {
2208     evtIndex = 0;
2209     vert[0] = 0. ;
2210     vert[1] = 0. ;
2211     vert[2] = 0. ;
2212   }
2213
2214   return evtIndex ;
2215 }
2216