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