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),             fhPrimPi0PtRejected(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),             fhPrimEtaPtRejected(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 fhMCPi0ProdVertexInner(0),   fhMCEtaProdVertexInner(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, Y<1",nptbins,ptmin,ptmax) ;
 692     fhPrimPi0AccE  = new TH1F("hPrimPi0AccE","Primary #pi^{0} 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     fhPrimPi0PtRejected = new TH1F("hPrimPi0PtRejected","Primary #pi^{0} pt",nptbins,ptmin,ptmax) ;
 699     fhPrimPi0Pt     = new TH1F("hPrimPi0Pt","Primary #pi^{0} #it{p}_{T} , #it{Y}<1",nptbins,ptmin,ptmax) ;
 700     fhPrimPi0AccPt  = new TH1F("hPrimPi0AccPt","Primary #pi^{0} pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
 701     fhPrimPi0Pt   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 702     fhPrimPi0AccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 703     outputContainer->Add(fhPrimPi0PtRejected) ;
 704     outputContainer->Add(fhPrimPi0Pt) ;
 705     outputContainer->Add(fhPrimPi0AccPt) ;
 706
 707     Int_t netabinsopen =  TMath::Nint(netabins*4/(etamax-etamin));
 708     fhPrimPi0Y      = new TH2F("hPrimPi0Rapidity","Rapidity of primary pi0",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
 709     fhPrimPi0Y   ->SetYTitle("#it{Rapidity}");
 710     fhPrimPi0Y   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 711     outputContainer->Add(fhPrimPi0Y) ;
 712
 713     fhPrimPi0Yeta      = new TH2F("hPrimPi0PseudoRapidity","PseudoRapidity of primary pi0",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
 714     fhPrimPi0Yeta   ->SetYTitle("#eta");
 715     fhPrimPi0Yeta   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 716     outputContainer->Add(fhPrimPi0Yeta) ;
 717
 718     fhPrimPi0YetaYcut      = new TH2F("hPrimPi0PseudoRapidityYcut","PseudoRapidity of primary pi0, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
 719     fhPrimPi0YetaYcut   ->SetYTitle("#eta");
 720     fhPrimPi0YetaYcut   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 721     outputContainer->Add(fhPrimPi0YetaYcut) ;
 722
 723     fhPrimPi0AccY   = new TH2F("hPrimPi0AccRapidity","Rapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
 724     fhPrimPi0AccY->SetYTitle("Rapidity");
 725     fhPrimPi0AccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 726     outputContainer->Add(fhPrimPi0AccY) ;
 727
 728     fhPrimPi0AccYeta      = new TH2F("hPrimPi0AccPseudoRapidity","PseudoRapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
 729     fhPrimPi0AccYeta   ->SetYTitle("#eta");
 730     fhPrimPi0AccYeta   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 731     outputContainer->Add(fhPrimPi0AccYeta) ;
 732
 733     Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
 734     fhPrimPi0Phi    = new TH2F("hPrimPi0Phi","Azimuthal of primary pi0, Y<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
 735     fhPrimPi0Phi->SetYTitle("#phi (deg)");
 736     fhPrimPi0Phi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 737     outputContainer->Add(fhPrimPi0Phi) ;
 738
 739     fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","Azimuthal of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,
 740                                nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 741     fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
 742     fhPrimPi0AccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 743     outputContainer->Add(fhPrimPi0AccPhi) ;
 744
 745     fhPrimPi0PtCentrality     = new TH2F("hPrimPi0PtCentrality","Primary #pi^{0} #it{p}_{T} vs reco centrality, Y<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
 746     fhPrimPi0AccPtCentrality  = new TH2F("hPrimPi0AccPtCentrality","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco centrality",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, Y<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
 755     fhPrimPi0AccPtEventPlane  = new TH2F("hPrimPi0AccPtEventPlane","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco event plane angle",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
 756     fhPrimPi0PtEventPlane   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 757     fhPrimPi0AccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 758     fhPrimPi0PtEventPlane   ->SetYTitle("Event Plane Angle (rad)");
 759     fhPrimPi0AccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
 760     outputContainer->Add(fhPrimPi0PtEventPlane) ;
 761     outputContainer->Add(fhPrimPi0AccPtEventPlane) ;
 762
 763     //Eta
 764
 765     fhPrimEtaE     = new TH1F("hPrimEtaE","Primary eta E",nptbins,ptmin,ptmax) ;
 766     fhPrimEtaAccE  = new TH1F("hPrimEtaAccE","Primary eta E with both photons in acceptance",nptbins,ptmin,ptmax) ;
 767     fhPrimEtaE   ->SetXTitle("#it{E} (GeV)");
 768     fhPrimEtaAccE->SetXTitle("#it{E} (GeV)");
 769     outputContainer->Add(fhPrimEtaE) ;
 770     outputContainer->Add(fhPrimEtaAccE) ;
 771
 772     fhPrimEtaPtRejected = new TH1F("hPrimEtaPtRejected","Primary eta pt",nptbins,ptmin,ptmax) ;
 773     fhPrimEtaPt     = new TH1F("hPrimEtaPt","Primary eta pt",nptbins,ptmin,ptmax) ;
 774     fhPrimEtaAccPt  = new TH1F("hPrimEtaAccPt","Primary eta #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
 775     fhPrimEtaPt   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 776     fhPrimEtaAccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 777     outputContainer->Add(fhPrimEtaPtRejected) ;
 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, |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, 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, Y<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
 826     fhPrimEtaAccPtEventPlane  = new TH2F("hPrimEtaAccPtEventPlane","Primary #eta with both photons 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",nptbins,ptmin,ptmax,2000,0,500) ;
 936       fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 937       fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
 938       outputContainer->Add(fhMCPi0ProdVertex) ;
 939
 940       fhMCPi0ProdVertexInner = new TH2F("hMCPi0ProdVertexInner","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",nptbins,ptmin,ptmax,500,0,50) ;
 941       fhMCPi0ProdVertexInner->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 942       fhMCPi0ProdVertexInner->SetYTitle("#it{R} (cm)");
 943       outputContainer->Add(fhMCPi0ProdVertexInner) ;
 944
 945       fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",nptbins,ptmin,ptmax,2000,0,500) ;
 946       fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 947       fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
 948       outputContainer->Add(fhMCEtaProdVertex) ;
 949
 950       fhMCEtaProdVertexInner = new TH2F("hMCEtaProdVertexInner","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",nptbins,ptmin,ptmax,500,0,50) ;
 951       fhMCEtaProdVertexInner->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 952       fhMCEtaProdVertexInner->SetYTitle("#it{R} (cm)");
 953       outputContainer->Add(fhMCEtaProdVertexInner) ;
 954
 955       for(Int_t i = 0; i<13; i++)
 956       {
 957         fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("#it{M} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 958         fhMCOrgMass[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 959         fhMCOrgMass[i]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
 960         outputContainer->Add(fhMCOrgMass[i]) ;
 961
 962         fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("#it{Asymmetry} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
 963         fhMCOrgAsym[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 964         fhMCOrgAsym[i]->SetYTitle("A");
 965         outputContainer->Add(fhMCOrgAsym[i]) ;
 966
 967         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) ;
 968         fhMCOrgDeltaEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 969         fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
 970         outputContainer->Add(fhMCOrgDeltaEta[i]) ;
 971
 972         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) ;
 973         fhMCOrgDeltaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 974         fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
 975         outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
 976
 977       }
 978
 979       if(fMultiCutAnaSim)
 980       {
 981         fhMCPi0MassPtTrue  = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 982         fhMCPi0MassPtRec   = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 983         fhMCPi0PtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 984         fhMCEtaMassPtRec   = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 985         fhMCEtaMassPtTrue  = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 986         fhMCEtaPtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 987         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
 988           for(Int_t icell=0; icell<fNCellNCuts; icell++){
 989             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 990               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
 991
 992               fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
 993                                                  Form("Reconstructed Mass vs reconstructed #it{p}_{T} of true #pi^{0} cluster pairs for #it{p}_{T} >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
 994                                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 995               fhMCPi0MassPtRec[index]->SetXTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
 996               fhMCPi0MassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
 997               outputContainer->Add(fhMCPi0MassPtRec[index]) ;
 998
 999               fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1000                                                   Form("Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs for #it{p}_{T} >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1001                                                   nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1002               fhMCPi0MassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1003               fhMCPi0MassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1004               outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
1005
1006               fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1007                                                    Form("Generated vs reconstructed p_T of true #pi^{0} 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]),
1008                                                    nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1009               fhMCPi0PtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1010               fhMCPi0PtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1011               outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
1012
1013               fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1014                                                  Form("Reconstructed Mass vs reconstructed p_T of true #eta cluster pairs for #it{p}_{T} >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1015                                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1016               fhMCEtaMassPtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1017               fhMCEtaMassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1018               outputContainer->Add(fhMCEtaMassPtRec[index]) ;
1019
1020               fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1021                                                   Form("Reconstructed Mass vs generated p_T of true #eta cluster pairs for #it{p}_{T} >%2.2f, ncell>%d and asym >%1.2f",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]),
1022                                                   nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1023               fhMCEtaMassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1024               fhMCEtaMassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1025               outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
1026
1027               fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1028                                                    Form("Generated vs reconstructed 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]),
1029                                                    nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1030               fhMCEtaPtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1031               fhMCEtaPtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1032               outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
1033             }
1034           }
1035         }
1036       }//multi cut ana
1037       else
1038       {
1039         fhMCPi0MassPtTrue  = new TH2F*[1];
1040         fhMCPi0PtTruePtRec = new TH2F*[1];
1041         fhMCEtaMassPtTrue  = new TH2F*[1];
1042         fhMCEtaPtTruePtRec = new TH2F*[1];
1043
1044         fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1045         fhMCPi0MassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1046         fhMCPi0MassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1047         outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
1048
1049         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) ;
1050         fhMCPi0PtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1051         fhMCPi0PtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1052         outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
1053
1054         fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1055         fhMCEtaMassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1056         fhMCEtaMassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1057         outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
1058
1059         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) ;
1060         fhMCEtaPtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1061         fhMCEtaPtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1062         outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
1063       }
1064     }
1065   }
1066
1067   if(fFillSMCombinations)
1068   {
1069     TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
1070     for(Int_t imod=0; imod<fNModules; imod++)
1071     {
1072       //Module dependent invariant mass
1073       snprintf(key, buffersize,"hReMod_%d",imod) ;
1074       snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1075       fhReMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1076       fhReMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1077       fhReMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1078       outputContainer->Add(fhReMod[imod]) ;
1079       if(fCalorimeter=="PHOS")
1080       {
1081         snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
1082         snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1083         fhReDiffPHOSMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1084         fhReDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1085         fhReDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1086         outputContainer->Add(fhReDiffPHOSMod[imod]) ;
1087       }
1088       else{//EMCAL
1089         if(imod<fNModules/2)
1090         {
1091           snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
1092           snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1093           fhReSameSectorEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1094           fhReSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1095           fhReSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1096           outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
1097         }
1098         if(imod<fNModules-2)
1099         {
1100           snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
1101           snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1102           fhReSameSideEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1103           fhReSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1104           fhReSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1105           outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
1106         }
1107       }//EMCAL
1108
1109       if(DoOwnMix())
1110       {
1111         snprintf(key, buffersize,"hMiMod_%d",imod) ;
1112         snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1113         fhMiMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1114         fhMiMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1115         fhMiMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1116         outputContainer->Add(fhMiMod[imod]) ;
1117
1118         if(fCalorimeter=="PHOS"){
1119           snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
1120           snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1121           fhMiDiffPHOSMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1122           fhMiDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1123           fhMiDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1124           outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
1125         }//PHOS
1126         else{//EMCAL
1127           if(imod<fNModules/2)
1128           {
1129             snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
1130             snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1131             fhMiSameSectorEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1132             fhMiSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1133             fhMiSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1134             outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
1135           }
1136           if(imod<fNModules-2){
1137
1138             snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
1139             snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1140             fhMiSameSideEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1141             fhMiSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1142             fhMiSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1143             outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
1144           }
1145         }//EMCAL
1146       }// own mix
1147     }//loop combinations
1148   } // SM combinations
1149
1150   if(fFillArmenterosThetaStar && IsDataMC())
1151   {
1152     TString ebin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
1153     Int_t narmbins = 400;
1154     Float_t armmin = 0;
1155     Float_t armmax = 0.4;
1156
1157     for(Int_t i = 0; i < 4; i++)
1158     {
1159       fhArmPrimPi0[i] =  new TH2F(Form("hArmenterosPrimPi0EBin%d",i),
1160                                   Form("Armenteros of primary #pi^{0}, %s",ebin[i].Data()),
1161                                   200, -1, 1, narmbins,armmin,armmax);
1162       fhArmPrimPi0[i]->SetYTitle("#it{p}_{T}^{Arm}");
1163       fhArmPrimPi0[i]->SetXTitle("#alpha^{Arm}");
1164       outputContainer->Add(fhArmPrimPi0[i]) ;
1165
1166       fhArmPrimEta[i] =  new TH2F(Form("hArmenterosPrimEtaEBin%d",i),
1167                                       Form("Armenteros of primary #eta, %s",ebin[i].Data()),
1168                                       200, -1, 1, narmbins,armmin,armmax);
1169       fhArmPrimEta[i]->SetYTitle("#it{p}_{T}^{Arm}");
1170       fhArmPrimEta[i]->SetXTitle("#alpha^{Arm}");
1171       outputContainer->Add(fhArmPrimEta[i]) ;
1172
1173     }
1174
1175     // Same as asymmetry ...
1176     fhCosThStarPrimPi0  = new TH2F
1177     ("hCosThStarPrimPi0","cos(#theta *) for primary #pi^{0}",nptbins,ptmin,ptmax,200,-1,1);
1178     fhCosThStarPrimPi0->SetYTitle("cos(#theta *)");
1179     fhCosThStarPrimPi0->SetXTitle("E_{ #pi^{0}} (GeV)");
1180     outputContainer->Add(fhCosThStarPrimPi0) ;
1181
1182     fhCosThStarPrimEta  = new TH2F
1183     ("hCosThStarPrimEta","cos(#theta *) for primary #eta",nptbins,ptmin,ptmax,200,-1,1);
1184     fhCosThStarPrimEta->SetYTitle("cos(#theta *)");
1185     fhCosThStarPrimEta->SetXTitle("E_{ #eta} (GeV)");
1186     outputContainer->Add(fhCosThStarPrimEta) ;
1187
1188   }
1189
1190   //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
1191   //
1192   //    printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
1193   //
1194   //  }
1195
1196   return outputContainer;
1197 }
1198
1199 //___________________________________________________
1200 void AliAnaPi0::Print(const Option_t * /*opt*/) const
1201 {
1202   //Print some relevant parameters set for the analysis
1203   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1204   AliAnaCaloTrackCorrBaseClass::Print(" ");
1205
1206   printf("Number of bins in Centrality:  %d \n",GetNCentrBin()) ;
1207   printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
1208   printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
1209   printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
1210   printf("Pair in same Module: %d \n",fSameSM) ;
1211   printf("Cuts: \n") ;
1212   // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
1213   printf("Number of modules:             %d \n",fNModules) ;
1214   printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
1215   printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
1216   printf("\tasymmetry < ");
1217   for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
1218   printf("\n");
1219
1220   printf("PID selection bits: n = %d, \n",fNPIDBits) ;
1221   printf("\tPID bit = ");
1222   for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
1223   printf("\n");
1224
1225   if(fMultiCutAna){
1226     printf("pT cuts: n = %d, \n",fNPtCuts) ;
1227     printf("\tpT > ");
1228     for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
1229     printf("GeV/c\n");
1230
1231     printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
1232     printf("\tnCell > ");
1233     for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
1234     printf("\n");
1235
1236   }
1237   printf("------------------------------------------------------\n") ;
1238 }
1239
1240 //________________________________________
1241 void AliAnaPi0::FillAcceptanceHistograms()
1242 {
1243   //Fill acceptance histograms if MC data is available
1244
1245   Float_t cen = GetEventCentrality();
1246   Float_t ep  = GetEventPlaneAngle();
1247
1248   if(GetReader()->ReadStack())
1249   {
1250     AliStack * stack = GetMCStack();
1251     if(stack)
1252     {
1253       for(Int_t i=0 ; i<stack->GetNtrack(); i++)
1254       {
1255         if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1256
1257         TParticle * prim = stack->Particle(i) ;
1258         Int_t pdg = prim->GetPdgCode();
1259         //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
1260         //                             prim->GetName(), prim->GetPdgCode());
1261
1262         if( pdg == 111 || pdg == 221)
1263         {
1264           Double_t pi0Pt = prim->Pt() ;
1265           Double_t pi0E  = prim->Energy() ;
1266           if(pi0E == TMath::Abs(prim->Pz()))
1267           {
1268             if( pdg == 111 ) fhPrimPi0PtRejected->Fill(pi0Pt);
1269             else             fhPrimEtaPtRejected->Fill(pi0Pt);
1270             continue ; //Protection against floating point exception
1271           }
1272
1273           Double_t pi0Y    = 0.5*TMath::Log((pi0E-prim->Pz())/(pi0E+prim->Pz())) ;
1274           Double_t pi0Yeta = prim->Eta() ;
1275           Double_t phi     = TMath::RadToDeg()*prim->Phi() ;
1276
1277           if(pdg == 111)
1278           {
1279             if(TMath::Abs(pi0Y) < 1.0)
1280             {
1281               fhPrimPi0E  ->Fill(pi0E ) ;
1282               fhPrimPi0Pt ->Fill(pi0Pt) ;
1283               fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1284               fhPrimPi0YetaYcut    ->Fill(pi0Pt,pi0Yeta) ;
1285               fhPrimPi0PtCentrality->Fill(pi0Pt,cen) ;
1286               fhPrimPi0PtEventPlane->Fill(pi0Pt,ep ) ;
1287             }
1288             fhPrimPi0Y   ->Fill(pi0Pt, pi0Y) ;
1289             fhPrimPi0Yeta->Fill(pi0Pt, pi0Yeta) ;
1290           }
1291           else if(pdg == 221)
1292           {
1293             if(TMath::Abs(pi0Y) < 1.0)
1294             {
1295               fhPrimEtaE  ->Fill(pi0E ) ;
1296               fhPrimEtaPt ->Fill(pi0Pt) ;
1297               fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1298               fhPrimEtaYetaYcut    ->Fill(pi0Pt,pi0Yeta) ;
1299               fhPrimEtaPtCentrality->Fill(pi0Pt,cen) ;
1300               fhPrimEtaPtEventPlane->Fill(pi0Pt,ep ) ;
1301             }
1302             fhPrimEtaY   ->Fill(pi0Pt, pi0Y) ;
1303             fhPrimEtaYeta->Fill(pi0Pt, pi0Yeta) ;
1304           }
1305
1306           //Origin of meson
1307           if(fFillOriginHisto)
1308           {
1309             Int_t momindex  = prim->GetFirstMother();
1310             if(momindex >= 0)
1311             {
1312               TParticle* mother = stack->Particle(momindex);
1313               Int_t mompdg    = TMath::Abs(mother->GetPdgCode());
1314               Int_t momstatus = mother->GetStatusCode();
1315               if(pdg == 111)
1316               {
1317                 if     (momstatus  == 21)fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1318                 else if(mompdg     < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1319                 else if(mompdg     > 2100  && mompdg   < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1320                 else if(mompdg    == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1321                 else if(mompdg    == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1322                 else if(mompdg    == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1323                 else if(mompdg    == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1324                 else if(mompdg    >= 310   && mompdg    <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1325                 else if(mompdg    == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1326                 else if(momstatus == 11 || momstatus  == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1327                 else                      fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1328               }//pi0
1329               else
1330               {
1331                 if     (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1332                 else if(mompdg    < 22  ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1333                 else if(mompdg    > 2100  && mompdg   < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1334                 else if(mompdg    == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1335                 else if(momstatus == 11 || momstatus  == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1336                 else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1337                 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1338               }
1339             } // pi0 has mother
1340           }
1341
1342           //Check if both photons hit Calorimeter
1343           if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
1344           Int_t iphot1=prim->GetFirstDaughter() ;
1345           Int_t iphot2=prim->GetLastDaughter() ;
1346           if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack())
1347           {
1348             TParticle * phot1 = stack->Particle(iphot1) ;
1349             TParticle * phot2 = stack->Particle(iphot2) ;
1350             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
1351             {
1352               //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
1353               //        phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
1354
1355               TLorentzVector lv1, lv2,lvmeson;
1356               phot1->Momentum(lv1);
1357               phot2->Momentum(lv2);
1358               prim ->Momentum(lvmeson);
1359
1360               if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
1361
1362               Bool_t inacceptance = kFALSE;
1363               if(fCalorimeter == "PHOS")
1364               {
1365                 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
1366                 {
1367                   Int_t mod ;
1368                   Double_t x,z ;
1369                   if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
1370                     inacceptance = kTRUE;
1371                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1372                 }
1373                 else
1374                 {
1375                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1376                     inacceptance = kTRUE ;
1377                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1378                 }
1379               }
1380               else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
1381               {
1382                 if(GetEMCALGeometry())
1383                 {
1384                   Int_t absID1=0;
1385                   Int_t absID2=0;
1386
1387                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1388                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1389
1390                   if( absID1 >= 0 && absID2 >= 0)
1391                     inacceptance = kTRUE;
1392
1393                   if(inacceptance && fCheckAccInSector)
1394                   {
1395                     Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
1396                     Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
1397
1398                     Int_t j=0;
1399                     Bool_t sameSector = kFALSE;
1400                     for(Int_t isector = 0; isector < fNModules/2; isector++)
1401                     {
1402                       j=2*isector;
1403                       if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
1404                     }
1405
1406                     if(sm1!=sm2 && !sameSector)  inacceptance = kFALSE;
1407
1408                     //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
1409                   }
1410
1411                   //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1412                   //                    inacceptance = kTRUE;
1413                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1414                 }
1415                 else
1416                 {
1417                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1418                     inacceptance = kTRUE ;
1419                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1420                 }
1421               }
1422
1423               if(inacceptance)
1424               {
1425                 Float_t  asym  = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
1426                 Double_t angle = lv1.Angle(lv2.Vect());
1427
1428                 if(pdg==111)
1429                 {
1430                   fhPrimPi0AccE   ->Fill(pi0E) ;
1431                   fhPrimPi0AccPt  ->Fill(pi0Pt) ;
1432                   fhPrimPi0AccPhi ->Fill(pi0Pt, phi) ;
1433                   fhPrimPi0AccY   ->Fill(pi0Pt, pi0Y) ;
1434                   fhPrimPi0AccYeta->Fill(pi0Pt,pi0Yeta) ;
1435                   fhPrimPi0AccPtCentrality->Fill(pi0Pt,cen) ;
1436                   fhPrimPi0AccPtEventPlane->Fill(pi0Pt,ep ) ;
1437
1438                   if(fFillAngleHisto)
1439                   {
1440                     fhPrimPi0OpeningAngle    ->Fill(pi0Pt,angle);
1441                     if(pi0Pt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
1442                     fhPrimPi0CosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1443                   }
1444                 }
1445                 else if(pdg==221)
1446                 {
1447                   fhPrimEtaAccE   ->Fill(pi0E ) ;
1448                   fhPrimEtaAccPt  ->Fill(pi0Pt) ;
1449                   fhPrimEtaAccPhi ->Fill(pi0Pt, phi) ;
1450                   fhPrimEtaAccY   ->Fill(pi0Pt, pi0Y) ;
1451                   fhPrimEtaAccYeta->Fill(pi0Pt, pi0Yeta) ;
1452                   fhPrimEtaAccPtCentrality->Fill(pi0Pt,cen) ;
1453                   fhPrimEtaAccPtEventPlane->Fill(pi0Pt,ep ) ;
1454
1455                   if(fFillAngleHisto)
1456                   {
1457                     fhPrimEtaOpeningAngle    ->Fill(pi0Pt,angle);
1458                     if(pi0Pt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
1459                     fhPrimEtaCosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1460                   }
1461                 }
1462               }//Accepted
1463             }// 2 photons
1464           }//Check daughters exist
1465         }// Primary pi0 or eta
1466       }//loop on primaries
1467     }//stack exists and data is MC
1468   }//read stack
1469   else if(GetReader()->ReadAODMCParticles())
1470   {
1471     TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1472     if(mcparticles)
1473     {
1474       Int_t nprim = mcparticles->GetEntriesFast();
1475
1476       for(Int_t i=0; i < nprim; i++)
1477       {
1478         if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1479
1480         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
1481
1482         // Only generator particles, when they come from PYTHIA, PHOJET, HERWIG ...
1483         //if( prim->GetStatus() == 0 && (GetMCAnalysisUtils()->GetMCGenerator()).Length()!=0) break;
1484
1485         Int_t pdg = prim->GetPdgCode();
1486         if( pdg == 111 || pdg == 221)
1487         {
1488           Double_t pi0Pt = prim->Pt() ;
1489           Double_t pi0E  = prim->E() ;
1490           //printf("pi0, pt %2.2f, eta %f, phi %f\n",pi0Pt, prim->Eta(), prim->Phi());
1491
1492           if(pi0E == TMath::Abs(prim->Pz()))
1493           {
1494             if( pdg == 111 ) fhPrimPi0PtRejected->Fill(pi0Pt);
1495             else             fhPrimEtaPtRejected->Fill(pi0Pt);
1496             continue ; //Protection against floating point exception
1497           }
1498
1499           Double_t pi0Y    = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
1500           Double_t pi0Yeta = prim->Eta() ;
1501           Double_t phi     = TMath::RadToDeg()*prim->Phi() ;
1502
1503           if(pdg == 111)
1504           {
1505             if(TMath::Abs(pi0Y) < 1)
1506             {
1507               fhPrimPi0E  ->Fill(pi0E ) ;
1508               fhPrimPi0Pt ->Fill(pi0Pt) ;
1509               fhPrimPi0Phi->Fill(pi0Pt, phi) ;
1510               fhPrimPi0YetaYcut    ->Fill(pi0Pt,pi0Yeta) ;
1511               fhPrimPi0PtCentrality->Fill(pi0Pt,cen) ;
1512               fhPrimPi0PtEventPlane->Fill(pi0Pt,ep ) ;
1513             }
1514             fhPrimPi0Y   ->Fill(pi0Pt, pi0Y   ) ;
1515             fhPrimPi0Yeta->Fill(pi0Pt, pi0Yeta) ;
1516           }
1517           else if(pdg == 221)
1518           {
1519             if(TMath::Abs(pi0Y) < 1)
1520             {
1521               fhPrimEtaE  ->Fill(pi0E ) ;
1522               fhPrimEtaPt ->Fill(pi0Pt) ;
1523               fhPrimEtaPhi->Fill(pi0Pt, phi) ;
1524               fhPrimEtaYetaYcut    ->Fill(pi0Pt,pi0Yeta) ;
1525               fhPrimEtaPtCentrality->Fill(pi0Pt,cen) ;
1526               fhPrimEtaPtEventPlane->Fill(pi0Pt,ep ) ;
1527             }
1528             fhPrimEtaY   ->Fill(pi0Pt, pi0Y   ) ;
1529             fhPrimEtaYeta->Fill(pi0Pt, pi0Yeta) ;
1530           }
1531
1532           //Origin of meson
1533           Int_t momindex  = prim->GetMother();
1534           if(momindex >= 0)
1535           {
1536             AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1537             Int_t mompdg    = TMath::Abs(mother->GetPdgCode());
1538             Int_t momstatus = mother->GetStatus();
1539             if(fFillOriginHisto)
1540             {
1541               if(pdg == 111)
1542               {
1543                 if     (momstatus  == 21) fhPrimPi0PtOrigin->Fill(pi0Pt,0.5);//parton
1544                 else if(mompdg     < 22 ) fhPrimPi0PtOrigin->Fill(pi0Pt,1.5);//quark
1545                 else if(mompdg     > 2100  && mompdg   < 2210) fhPrimPi0PtOrigin->Fill(pi0Pt,2.5);// resonances
1546                 else if(mompdg    == 221) fhPrimPi0PtOrigin->Fill(pi0Pt,8.5);//eta
1547                 else if(mompdg    == 331) fhPrimPi0PtOrigin->Fill(pi0Pt,9.5);//eta prime
1548                 else if(mompdg    == 213) fhPrimPi0PtOrigin->Fill(pi0Pt,4.5);//rho
1549                 else if(mompdg    == 223) fhPrimPi0PtOrigin->Fill(pi0Pt,5.5);//omega
1550                 else if(mompdg    >= 310   && mompdg    <= 323) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0S, k+-,k*
1551                 else if(mompdg    == 130) fhPrimPi0PtOrigin->Fill(pi0Pt,6.5);//k0L
1552                 else if(momstatus == 11 || momstatus  == 12 ) fhPrimPi0PtOrigin->Fill(pi0Pt,3.5);//resonances
1553                 else                      fhPrimPi0PtOrigin->Fill(pi0Pt,7.5);//other?
1554               }//pi0
1555               else
1556               {
1557                 if     (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(pi0Pt,0.5);//parton
1558                 else if(mompdg    < 22  ) fhPrimEtaPtOrigin->Fill(pi0Pt,1.5);//quark
1559                 else if(mompdg    > 2100  && mompdg   < 2210) fhPrimEtaPtOrigin->Fill(pi0Pt,2.5);//qq resonances
1560                 else if(mompdg    == 331) fhPrimEtaPtOrigin->Fill(pi0Pt,5.5);//eta prime
1561                 else if(momstatus == 11 || momstatus  == 12 ) fhPrimEtaPtOrigin->Fill(pi0Pt,3.5);//resonances
1562                 else fhPrimEtaPtOrigin->Fill(pi0Pt,4.5);//stable, conversions?
1563                 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1564               }
1565             }
1566           }//pi0 has mother
1567
1568           //Check if both photons hit Calorimeter
1569           if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
1570           Int_t iphot1=prim->GetDaughter(0) ;
1571           Int_t iphot2=prim->GetDaughter(1) ;
1572           if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim)
1573           {
1574             AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
1575             AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
1576             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
1577             {
1578               TLorentzVector lv1, lv2,lvmeson;
1579               lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1580               lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1581               lvmeson.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
1582
1583                if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
1584
1585               Bool_t inacceptance = kFALSE;
1586               if(fCalorimeter == "PHOS")
1587               {
1588                 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
1589                 {
1590                   Int_t mod ;
1591                   Double_t x,z ;
1592                   Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
1593                   Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
1594                   if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) &&
1595                      GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x))
1596                     inacceptance = kTRUE;
1597                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1598                 }
1599                 else
1600                 {
1601                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1602                     inacceptance = kTRUE ;
1603                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1604                 }
1605               }
1606               else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
1607               {
1608                 if(GetEMCALGeometry())
1609                 {
1610                   Int_t absID1=0;
1611                   Int_t absID2=0;
1612
1613                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1614                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1615
1616                   if( absID1 >= 0 && absID2 >= 0)
1617                     inacceptance = kTRUE;
1618
1619                   if(inacceptance && fCheckAccInSector)
1620                   {
1621                     Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
1622                     Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
1623
1624                     Int_t j=0;
1625                     Bool_t sameSector = kFALSE;
1626                     for(Int_t isector = 0; isector < fNModules/2; isector++)
1627                     {
1628                       j=2*isector;
1629                       if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
1630                     }
1631
1632                     if(sm1!=sm2 && !sameSector)  inacceptance = kFALSE;
1633
1634                     //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
1635                   }
1636
1637                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1638                 }
1639                 else
1640                 {
1641                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1642                     inacceptance = kTRUE ;
1643                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1644                 }
1645               }
1646
1647               if(inacceptance)
1648               {
1649                 Float_t  asym  = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
1650                 Double_t angle = lv1.Angle(lv2.Vect());
1651
1652                 if(pdg==111)
1653                 {
1654                   //                printf("ACCEPTED pi0: pt %2.2f, phi %3.2f, eta %1.2f\n",pi0Pt,phi,pi0Y);
1655                   fhPrimPi0AccE   ->Fill(pi0E ) ;
1656                   fhPrimPi0AccPt  ->Fill(pi0Pt) ;
1657                   fhPrimPi0AccPhi ->Fill(pi0Pt, phi) ;
1658                   fhPrimPi0AccY   ->Fill(pi0Pt, pi0Y) ;
1659                   fhPrimPi0AccYeta->Fill(pi0Pt, pi0Yeta) ;
1660                   fhPrimPi0AccPtCentrality->Fill(pi0Pt,cen) ;
1661                   fhPrimPi0AccPtEventPlane->Fill(pi0Pt,ep ) ;
1662
1663                   if(fFillAngleHisto)
1664                   {
1665                     fhPrimPi0OpeningAngle    ->Fill(pi0Pt,angle);
1666                     if(pi0Pt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
1667                     fhPrimPi0CosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1668                   }
1669                 }
1670                 else if(pdg==221)
1671                 {
1672                   fhPrimEtaAccE   ->Fill(pi0E ) ;
1673                   fhPrimEtaAccPt  ->Fill(pi0Pt) ;
1674                   fhPrimEtaAccPhi ->Fill(pi0Pt, phi) ;
1675                   fhPrimEtaAccY   ->Fill(pi0Pt, pi0Y) ;
1676                   fhPrimEtaAccYeta->Fill(pi0Pt, pi0Yeta) ;
1677                   fhPrimEtaAccPtCentrality->Fill(pi0Pt,cen) ;
1678                   fhPrimEtaAccPtEventPlane->Fill(pi0Pt,ep ) ;
1679
1680                   if(fFillAngleHisto)
1681                   {
1682                     fhPrimEtaOpeningAngle    ->Fill(pi0Pt,angle);
1683                     if(pi0Pt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
1684                     fhPrimEtaCosOpeningAngle ->Fill(pi0Pt,TMath::Cos(angle));
1685                   }
1686                 }
1687               }//Accepted
1688             }// 2 photons
1689           }//Check daughters exist
1690         }// Primary pi0 or eta
1691       }//loop on primaries
1692     }//stack exists and data is MC
1693
1694
1695   }     // read AOD MC
1696 }
1697
1698 //__________________________________________________________________________________
1699 void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg,             TLorentzVector meson,
1700                                         TLorentzVector daugh1, TLorentzVector daugh2)
1701 {
1702   // Fill armenteros plots
1703
1704   // Get pTArm and AlphaArm
1705   Float_t momentumSquaredMother = meson.P()*meson.P();
1706   Float_t momentumDaughter1AlongMother = 0.;
1707   Float_t momentumDaughter2AlongMother = 0.;
1708
1709   if (momentumSquaredMother > 0.)
1710   {
1711     momentumDaughter1AlongMother = (daugh1.Px()*meson.Px() + daugh1.Py()*meson.Py()+ daugh1.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
1712     momentumDaughter2AlongMother = (daugh2.Px()*meson.Px() + daugh2.Py()*meson.Py()+ daugh2.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
1713   }
1714
1715   Float_t momentumSquaredDaughter1 = daugh1.P()*daugh1.P();
1716   Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
1717
1718   Float_t pTArm = 0.;
1719   if (ptArmSquared > 0.)
1720     pTArm = sqrt(ptArmSquared);
1721
1722   Float_t alphaArm = 0.;
1723   if(momentumDaughter1AlongMother +momentumDaughter2AlongMother > 0)
1724     alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
1725
1726   TLorentzVector daugh1Boost = daugh1;
1727   daugh1Boost.Boost(-meson.BoostVector());
1728   Float_t  cosThStar=TMath::Cos(daugh1Boost.Vect().Angle(meson.Vect()));
1729
1730   Float_t en   = meson.Energy();
1731   Int_t   ebin = -1;
1732   if(en > 8  && en <= 12) ebin = 0;
1733   if(en > 12 && en <= 16) ebin = 1;
1734   if(en > 16 && en <= 20) ebin = 2;
1735   if(en > 20)             ebin = 3;
1736   if(ebin < 0 || ebin > 3) return ;
1737
1738
1739   if(pdg==111)
1740   {
1741     fhCosThStarPrimPi0->Fill(en,cosThStar);
1742     fhArmPrimPi0[ebin]->Fill(alphaArm,pTArm);
1743   }
1744   else
1745   {
1746     fhCosThStarPrimEta->Fill(en,cosThStar);
1747     fhArmPrimEta[ebin]->Fill(alphaArm,pTArm);
1748   }
1749
1750   if(GetDebug() > 2 )
1751   {
1752     Float_t asym = 0;
1753     if(daugh1.E()+daugh2.E() > 0) asym = TMath::Abs(daugh1.E()-daugh2.E())/(daugh1.E()+daugh2.E());
1754
1755     printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
1756          en,alphaArm,pTArm,cosThStar,asym);
1757   }
1758 }
1759
1760 //_______________________________________________________________________________________
1761 void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1,  Int_t index2,
1762                                               Float_t pt1,   Float_t pt2,
1763                                               Int_t ncell1,  Int_t ncell2,
1764                                               Double_t mass, Double_t pt,  Double_t asym,
1765                                               Double_t deta, Double_t dphi)
1766 {
1767   //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1768   //Adjusted for Pythia, need to see what to do for other generators.
1769   //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1770   // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1771
1772   if(!fFillOriginHisto) return;
1773
1774   Int_t ancPDG    = 0;
1775   Int_t ancStatus = 0;
1776   TLorentzVector ancMomentum;
1777   TVector3 prodVertex;
1778   Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1779                                                               GetReader(), ancPDG, ancStatus,ancMomentum, prodVertex);
1780
1781   Int_t momindex  = -1;
1782   Int_t mompdg    = -1;
1783   Int_t momstatus = -1;
1784   if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1785                             ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1786
1787   Float_t prodR = -1;
1788
1789   if(ancLabel > -1)
1790   {
1791     if(ancPDG==22){//gamma
1792       fhMCOrgMass[0]->Fill(pt,mass);
1793       fhMCOrgAsym[0]->Fill(pt,asym);
1794       fhMCOrgDeltaEta[0]->Fill(pt,deta);
1795       fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1796     }
1797     else if(TMath::Abs(ancPDG)==11){//e
1798       fhMCOrgMass[1]->Fill(pt,mass);
1799       fhMCOrgAsym[1]->Fill(pt,asym);
1800       fhMCOrgDeltaEta[1]->Fill(pt,deta);
1801       fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1802     }
1803     else if(ancPDG==111){//Pi0
1804       fhMCOrgMass[2]->Fill(pt,mass);
1805       fhMCOrgAsym[2]->Fill(pt,asym);
1806       fhMCOrgDeltaEta[2]->Fill(pt,deta);
1807       fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1808       if(fMultiCutAnaSim)
1809       {
1810         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1811           for(Int_t icell=0; icell<fNCellNCuts; icell++){
1812             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1813               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1814               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1815                  asym   <  fAsymCuts[iasym]                                   &&
1816                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1817                 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1818                 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1819                 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1820               }//pass the different cuts
1821             }// pid bit cut loop
1822           }// icell loop
1823         }// pt cut loop
1824       }//Multi cut ana sim
1825       else
1826       {
1827         fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1828
1829         if(mass < 0.17 && mass > 0.1)
1830         {
1831           fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1832
1833           if(fFillOriginHisto)
1834           {
1835             if(GetReader()->ReadStack())
1836             {
1837               TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1838               momindex  = ancestor->GetFirstMother();
1839               if(momindex < 0) return;
1840               TParticle* mother = GetMCStack()->Particle(momindex);
1841               mompdg    = TMath::Abs(mother->GetPdgCode());
1842               momstatus = mother->GetStatusCode();
1843               prodR = mother->R();
1844             }
1845             else
1846             {
1847               TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1848               AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1849               momindex  = ancestor->GetMother();
1850               if(momindex < 0) return;
1851               AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1852               mompdg    = TMath::Abs(mother->GetPdgCode());
1853               momstatus = mother->GetStatus();
1854               prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1855             }
1856
1857             fhMCPi0ProdVertex->Fill(pt,prodR);
1858             if(prodR < 50)fhMCPi0ProdVertexInner->Fill(pt,prodR);
1859
1860             if     (momstatus  == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1861             else if(mompdg     < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1862             else if(mompdg     > 2100  && mompdg   < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1863             else if(mompdg    == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1864             else if(mompdg    == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1865             else if(mompdg    == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1866             else if(mompdg    == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1867             else if(mompdg    >= 310   && mompdg    <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1868             else if(mompdg    == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1869             else if(momstatus == 11 || momstatus  == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1870             else                      fhMCPi0PtOrigin->Fill(pt,7.5);//other?
1871
1872           }
1873         }//pi0 mass region
1874       }
1875     }
1876     else if(ancPDG==221){//Eta
1877       fhMCOrgMass[3]->Fill(pt,mass);
1878       fhMCOrgAsym[3]->Fill(pt,asym);
1879       fhMCOrgDeltaEta[3]->Fill(pt,deta);
1880       fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1881       if(fMultiCutAnaSim){
1882         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1883           for(Int_t icell=0; icell<fNCellNCuts; icell++){
1884             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1885               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1886               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1887                  asym   <  fAsymCuts[iasym]                                   &&
1888                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1889                 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1890                 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1891                 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1892               }//pass the different cuts
1893             }// pid bit cut loop
1894           }// icell loop
1895         }// pt cut loop
1896       } //Multi cut ana sim
1897       else
1898       {
1899         fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1900         if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1901
1902         if(fFillOriginHisto)
1903         {
1904           if(GetReader()->ReadStack())
1905           {
1906             TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1907             momindex  = ancestor->GetFirstMother();
1908             if(momindex < 0) return;
1909             TParticle* mother = GetMCStack()->Particle(momindex);
1910             mompdg    = TMath::Abs(mother->GetPdgCode());
1911             momstatus = mother->GetStatusCode();
1912             prodR = mother->R();
1913           }
1914           else
1915           {
1916             TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1917             AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1918             momindex  = ancestor->GetMother();
1919             if(momindex < 0) return;
1920             AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1921             mompdg    = TMath::Abs(mother->GetPdgCode());
1922             momstatus = mother->GetStatus();
1923             prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1924           }
1925
1926           fhMCEtaProdVertex->Fill(pt,prodR);
1927           if(prodR < 50)fhMCEtaProdVertexInner->Fill(pt,prodR);
1928
1929           if     (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1930           else if(mompdg    < 22  ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1931           else if(mompdg    > 2100  && mompdg   < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1932           else if(mompdg    == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1933           else if(momstatus == 11 || momstatus  == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1934           else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1935           //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1936         }
1937       }// eta mass region
1938     }
1939     else if(ancPDG==-2212){//AProton
1940       fhMCOrgMass[4]->Fill(pt,mass);
1941       fhMCOrgAsym[4]->Fill(pt,asym);
1942       fhMCOrgDeltaEta[4]->Fill(pt,deta);
1943       fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1944     }
1945     else if(ancPDG==-2112){//ANeutron
1946       fhMCOrgMass[5]->Fill(pt,mass);
1947       fhMCOrgAsym[5]->Fill(pt,asym);
1948       fhMCOrgDeltaEta[5]->Fill(pt,deta);
1949       fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1950     }
1951     else if(TMath::Abs(ancPDG)==13){//muons
1952       fhMCOrgMass[6]->Fill(pt,mass);
1953       fhMCOrgAsym[6]->Fill(pt,asym);
1954       fhMCOrgDeltaEta[6]->Fill(pt,deta);
1955       fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1956     }
1957     else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1958       if(ancStatus==1){//Stable particles, converted? not decayed resonances
1959         fhMCOrgMass[6]->Fill(pt,mass);
1960         fhMCOrgAsym[6]->Fill(pt,asym);
1961         fhMCOrgDeltaEta[6]->Fill(pt,deta);
1962         fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1963       }
1964       else{//resonances and other decays, more hadron conversions?
1965         fhMCOrgMass[7]->Fill(pt,mass);
1966         fhMCOrgAsym[7]->Fill(pt,asym);
1967         fhMCOrgDeltaEta[7]->Fill(pt,deta);
1968         fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
1969       }
1970     }
1971     else {//Partons, colliding protons, strings, intermediate corrections
1972       if(ancStatus==11 || ancStatus==12){//String fragmentation
1973         fhMCOrgMass[8]->Fill(pt,mass);
1974         fhMCOrgAsym[8]->Fill(pt,asym);
1975         fhMCOrgDeltaEta[8]->Fill(pt,deta);
1976         fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1977       }
1978       else if (ancStatus==21){
1979         if(ancLabel < 2) {//Colliding protons
1980           fhMCOrgMass[11]->Fill(pt,mass);
1981           fhMCOrgAsym[11]->Fill(pt,asym);
1982           fhMCOrgDeltaEta[11]->Fill(pt,deta);
1983           fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1984         }//colliding protons
1985         else if(ancLabel < 6){//partonic initial states interactions
1986           fhMCOrgMass[9]->Fill(pt,mass);
1987           fhMCOrgAsym[9]->Fill(pt,asym);
1988           fhMCOrgDeltaEta[9]->Fill(pt,deta);
1989           fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1990         }
1991         else if(ancLabel < 8){//Final state partons radiations?
1992           fhMCOrgMass[10]->Fill(pt,mass);
1993           fhMCOrgAsym[10]->Fill(pt,asym);
1994           fhMCOrgDeltaEta[10]->Fill(pt,deta);
1995           fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1996         }
1997         // else {
1998         //   printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1999         //          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2000         // }
2001       }//status 21
2002       //else {
2003       //  printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2004       //         ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2005       // }
2006     }////Partons, colliding protons, strings, intermediate corrections
2007   }//ancLabel > -1
2008   else { //ancLabel <= -1
2009     //printf("Not related at all label = %d\n",ancLabel);
2010     fhMCOrgMass[12]->Fill(pt,mass);
2011     fhMCOrgAsym[12]->Fill(pt,asym);
2012     fhMCOrgDeltaEta[12]->Fill(pt,deta);
2013     fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
2014   }
2015 }
2016
2017 //__________________________________________
2018 void AliAnaPi0::MakeAnalysisFillHistograms()
2019 {
2020   //Process one event and extract photons from AOD branch
2021   // filled with AliAnaPhoton and fill histos with invariant mass
2022
2023   //In case of simulated data, fill acceptance histograms
2024   if(IsDataMC())FillAcceptanceHistograms();
2025
2026   //if (GetReader()->GetEventNumber()%10000 == 0)
2027   // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
2028
2029   if(!GetInputAODBranch())
2030   {
2031     printf("AliAnaPi0::MakeAnalysisFillHistograms() - No input aod photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
2032     abort();
2033   }
2034
2035   //Init some variables
2036   Int_t   nPhot    = GetInputAODBranch()->GetEntriesFast() ;
2037
2038   if(GetDebug() > 1)
2039     printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
2040
2041   //If less than photon 2 entries in the list, skip this event
2042   if(nPhot < 2 ) {
2043
2044     if(GetDebug() > 2)
2045       printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
2046
2047     if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
2048
2049     return ;
2050   }
2051
2052   Int_t ncentr = GetNCentrBin();
2053   if(GetNCentrBin()==0) ncentr = 1;
2054
2055   //Init variables
2056   Int_t module1         = -1;
2057   Int_t module2         = -1;
2058   Double_t vert[]       = {0.0, 0.0, 0.0} ; //vertex
2059   Int_t evtIndex1       = 0 ;
2060   Int_t currentEvtIndex = -1;
2061   Int_t curCentrBin     = GetEventCentralityBin();
2062   //Int_t curVzBin        = GetEventVzBin();
2063   //Int_t curRPBin        = GetEventRPBin();
2064   Int_t eventbin        = GetEventMixBin();
2065
2066   if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
2067   {
2068      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());
2069     return;
2070   }
2071
2072   //Get shower shape information of clusters
2073   TObjArray *clusters = 0;
2074   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2075   else if(fCalorimeter=="PHOS" ) clusters = GetPHOSClusters() ;
2076
2077   //---------------------------------
2078   //First loop on photons/clusters
2079   //---------------------------------
2080   for(Int_t i1=0; i1<nPhot-1; i1++)
2081   {
2082     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2083     //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
2084
2085     // get the event index in the mixed buffer where the photon comes from
2086     // in case of mixing with analysis frame, not own mixing
2087     evtIndex1 = GetEventIndex(p1, vert) ;
2088     //printf("charge = %d\n", track->Charge());
2089     if ( evtIndex1 == -1 )
2090       return ;
2091     if ( evtIndex1 == -2 )
2092       continue ;
2093
2094     //printf("z vertex %f < %f\n",vert[2],GetZvertexCut());
2095     if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ;   //vertex cut
2096
2097
2098     if (evtIndex1 != currentEvtIndex)
2099     {
2100       //Fill event bin info
2101       if(DoOwnMix()) fhEventBin->Fill(eventbin) ;
2102       if(GetNCentrBin() > 1)
2103       {
2104         fhCentrality->Fill(curCentrBin);
2105         if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
2106       }
2107       currentEvtIndex = evtIndex1 ;
2108     }
2109
2110     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
2111
2112     //Get the momentum of this cluster
2113     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2114
2115     //Get (Super)Module number of this cluster
2116     module1 = GetModuleNumber(p1);
2117
2118     //------------------------------------------
2119     //Get index in VCaloCluster array
2120     AliVCluster *cluster1 = 0;
2121     Bool_t bFound1        = kFALSE;
2122     Int_t  caloLabel1     = p1->GetCaloLabel(0);
2123     Bool_t iclus1         =-1;
2124     if(clusters)
2125     {
2126       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
2127         AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
2128         if(cluster)
2129         {
2130           if     (cluster->GetID()==caloLabel1)
2131           {
2132             bFound1  = kTRUE  ;
2133             cluster1 = cluster;
2134             iclus1   = iclus;
2135           }
2136         }
2137         if(bFound1) break;
2138       }
2139     }// calorimeter clusters loop
2140
2141     //---------------------------------
2142     //Second loop on photons/clusters
2143     //---------------------------------
2144     for(Int_t i2=i1+1; i2<nPhot; i2++)
2145     {
2146       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
2147
2148       //In case of mixing frame, check we are not in the same event as the first cluster
2149       Int_t evtIndex2 = GetEventIndex(p2, vert) ;
2150       if ( evtIndex2 == -1 )
2151         return ;
2152       if ( evtIndex2 == -2 )
2153         continue ;
2154       if (GetMixedEvent() && (evtIndex1 == evtIndex2))
2155         continue ;
2156
2157       //------------------------------------------
2158       //Get index in VCaloCluster array
2159       AliVCluster *cluster2 = 0;
2160       Bool_t bFound2        = kFALSE;
2161       Int_t caloLabel2      = p2->GetCaloLabel(0);
2162       if(clusters){
2163         for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
2164           AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
2165           if(cluster){
2166             if(cluster->GetID()==caloLabel2) {
2167               bFound2  = kTRUE  ;
2168               cluster2 = cluster;
2169             }
2170           }
2171           if(bFound2) break;
2172         }// calorimeter clusters loop
2173       }
2174
2175       Float_t tof1  = -1;
2176       Float_t l01   = -1;
2177       if(cluster1 && bFound1){
2178         tof1  = cluster1->GetTOF()*1e9;
2179         l01   = cluster1->GetM02();
2180       }
2181       //      else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
2182       //                   p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
2183
2184       Float_t tof2  = -1;
2185       Float_t l02   = -1;
2186       if(cluster2 && bFound2)
2187       {
2188         tof2  = cluster2->GetTOF()*1e9;
2189         l02   = cluster2->GetM02();
2190
2191       }
2192       //      else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
2193       //                  p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
2194
2195       if(clusters)
2196       {
2197         Double_t t12diff = tof1-tof2;
2198         if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2199       }
2200       //------------------------------------------
2201
2202       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
2203
2204       //Get the momentum of this cluster
2205       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2206       //Get module number
2207       module2       = GetModuleNumber(p2);
2208
2209       //---------------------------------
2210       // Get pair kinematics
2211       //---------------------------------
2212       Double_t m    = (photon1 + photon2).M() ;
2213       Double_t pt   = (photon1 + photon2).Pt();
2214       Double_t deta = photon1.Eta() - photon2.Eta();
2215       Double_t dphi = photon1.Phi() - photon2.Phi();
2216       Double_t a    = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2217
2218       if(GetDebug() > 2)
2219         printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
2220
2221       //--------------------------------
2222       // Opening angle selection
2223       //--------------------------------
2224       //Check if opening angle is too large or too small compared to what is expected
2225       Double_t angle   = photon1.Angle(photon2.Vect());
2226       if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
2227         if(GetDebug() > 2)
2228           printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2229         continue;
2230       }
2231
2232       if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2233         if(GetDebug() > 2)
2234           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
2235         continue;
2236       }
2237
2238       //-------------------------------------------------------------------------------------------------
2239       //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2240       //-------------------------------------------------------------------------------------------------
2241       if(a < fAsymCuts[0] && fFillSMCombinations)
2242       {
2243         if(module1==module2 && module1 >=0 && module1<fNModules)
2244           fhReMod[module1]->Fill(pt,m) ;
2245
2246         if(fCalorimeter=="EMCAL")
2247         {
2248           // Same sector
2249           Int_t j=0;
2250           for(Int_t i = 0; i < fNModules/2; i++)
2251           {
2252             j=2*i;
2253             if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
2254           }
2255
2256           // Same side
2257           for(Int_t i = 0; i < fNModules-2; i++){
2258             if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
2259           }
2260         }//EMCAL
2261         else {//PHOS
2262           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
2263           if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
2264           if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
2265         }//PHOS
2266       }
2267
2268       //In case we want only pairs in same (super) module, check their origin.
2269       Bool_t ok = kTRUE;
2270       if(fSameSM && module1!=module2) ok=kFALSE;
2271       if(ok)
2272       {
2273         //Check if one of the clusters comes from a conversion
2274         if(fCheckConversion)
2275         {
2276           if     (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
2277           else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
2278         }
2279
2280         // Fill shower shape cut histograms
2281         if(fFillSSCombinations)
2282         {
2283           if     ( l01 > 0.01 && l01 < 0.4  &&
2284                    l02 > 0.01 && l02 < 0.4 )               fhReSS[0]->Fill(pt,m); // Tight
2285           else if( l01 > 0.4  && l02 > 0.4 )               fhReSS[1]->Fill(pt,m); // Loose
2286           else if( l01 > 0.01 && l01 < 0.4  && l02 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2287           else if( l02 > 0.01 && l02 < 0.4  && l01 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2288         }
2289
2290         //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2291         for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2292           if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
2293             for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
2294               if(a < fAsymCuts[iasym])
2295               {
2296                 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2297                 //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);
2298
2299                 if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2300
2301                 fhRe1     [index]->Fill(pt,m);
2302                 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
2303                 if(fFillBadDistHisto){
2304                   if(p1->DistToBad()>0 && p2->DistToBad()>0){
2305                     fhRe2     [index]->Fill(pt,m) ;
2306                     if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
2307                     if(p1->DistToBad()>1 && p2->DistToBad()>1){
2308                       fhRe3     [index]->Fill(pt,m) ;
2309                       if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
2310                     }// bad 3
2311                   }// bad2
2312                 }// Fill bad dist histos
2313               }//assymetry cut
2314             }// asymmetry cut loop
2315           }// bad 1
2316         }// pid bit loop
2317
2318         //Fill histograms with opening angle
2319         if(fFillAngleHisto)
2320         {
2321           fhRealOpeningAngle   ->Fill(pt,angle);
2322           fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2323         }
2324
2325         //Fill histograms with pair assymmetry
2326         if(fFillAsymmetryHisto)
2327         {
2328           fhRePtAsym->Fill(pt,a);
2329           if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
2330           if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
2331         }
2332
2333         //-------------------------------------------------------
2334         //Get the number of cells needed for multi cut analysis.
2335         //-------------------------------------------------------
2336         Int_t ncell1 = 0;
2337         Int_t ncell2 = 0;
2338         if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim))
2339         {
2340           AliVEvent * event = GetReader()->GetInputEvent();
2341           if(event){
2342             for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++)
2343             {
2344               AliVCluster *cluster = event->GetCaloCluster(iclus);
2345
2346               Bool_t is = kFALSE;
2347               if     (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
2348               else if(fCalorimeter == "PHOS"  && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
2349
2350               if(is){
2351                 if      (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
2352                 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
2353               } // PHOS or EMCAL cluster as requested in analysis
2354
2355               if(ncell2 > 0 &&  ncell1 > 0) break; // No need to continue the iteration
2356
2357             }
2358             //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
2359           }
2360         }
2361
2362         //---------
2363         // MC data
2364         //---------
2365         //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2366         if(IsDataMC())
2367         {
2368           if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2369              GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2370           {
2371             fhReMCFromConversion->Fill(pt,m);
2372           }
2373           else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2374                   !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2375           {
2376             fhReMCFromNotConversion->Fill(pt,m);
2377           }
2378           else
2379           {
2380             fhReMCFromMixConversion->Fill(pt,m);
2381           }
2382
2383           FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2384         }
2385
2386         //-----------------------
2387         //Multi cuts analysis
2388         //-----------------------
2389         if(fMultiCutAna)
2390         {
2391           //Histograms for different PID bits selection
2392           for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2393
2394             if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    &&
2395                p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
2396
2397             //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2398           } // pid bit cut loop
2399
2400           //Several pt,ncell and asymmetry cuts
2401           for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2402             for(Int_t icell=0; icell<fNCellNCuts; icell++){
2403               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2404                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2405                 if(p1->E() >   fPtCuts[ipt]      && p2->E() > fPtCuts[ipt]        &&
2406                    a        <   fAsymCuts[iasym]                                    &&
2407                    ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell]){
2408                   fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
2409                   //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2410                   if(fFillSMCombinations && module1==module2){
2411                     fhRePtNCellAsymCutsSM[module1][index]->Fill(pt,m) ;
2412                   }
2413                 }
2414               }// pid bit cut loop
2415             }// icell loop
2416           }// pt cut loop
2417           if(GetHistogramRanges()->GetHistoTrackMultiplicityBins()){
2418             for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
2419               if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
2420             }
2421           }
2422         }// multiple cuts analysis
2423       }// ok if same sm
2424     }// second same event particle
2425   }// first cluster
2426
2427   //-------------------------------------------------------------
2428   // Mixing
2429   //-------------------------------------------------------------
2430   if(DoOwnMix())
2431   {
2432     //Recover events in with same characteristics as the current event
2433
2434     //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2435     if(eventbin < 0) return ;
2436
2437     TList * evMixList=fEventsList[eventbin] ;
2438
2439     if(!evMixList)
2440     {
2441       printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mix event list not available, bin %d \n",eventbin);
2442       return;
2443     }
2444
2445     Int_t nMixed = evMixList->GetSize() ;
2446     for(Int_t ii=0; ii<nMixed; ii++)
2447     {
2448       TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2449       Int_t nPhot2=ev2->GetEntriesFast() ;
2450       Double_t m = -999;
2451       if(GetDebug() > 1)
2452         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, GetEventCentralityBin());
2453
2454       fhEventMixBin->Fill(eventbin) ;
2455
2456       //---------------------------------
2457       //First loop on photons/clusters
2458       //---------------------------------
2459       for(Int_t i1=0; i1<nPhot; i1++){
2460         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2461         if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2462
2463         //Get kinematics of cluster and (super) module of this cluster
2464         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2465         module1 = GetModuleNumber(p1);
2466
2467         //---------------------------------
2468         //First loop on photons/clusters
2469         //---------------------------------
2470         for(Int_t i2=0; i2<nPhot2; i2++){
2471           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2472
2473           //Get kinematics of second cluster and calculate those of the pair
2474           TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2475           m           = (photon1+photon2).M() ;
2476           Double_t pt = (photon1 + photon2).Pt();
2477           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2478
2479           //Check if opening angle is too large or too small compared to what is expected
2480           Double_t angle   = photon1.Angle(photon2.Vect());
2481           if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
2482             if(GetDebug() > 2)
2483               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2484             continue;
2485           }
2486           if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2487             if(GetDebug() > 2)
2488               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2489             continue;
2490
2491           }
2492
2493           if(GetDebug() > 2)
2494             printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
2495                    p1->Pt(), p2->Pt(), pt,m,a);
2496
2497           //In case we want only pairs in same (super) module, check their origin.
2498           module2 = GetModuleNumber(p2);
2499
2500           //-------------------------------------------------------------------------------------------------
2501           //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2502           //-------------------------------------------------------------------------------------------------
2503           if(a < fAsymCuts[0] && fFillSMCombinations){
2504             if(module1==module2 && module1 >=0 && module1<fNModules)
2505               fhMiMod[module1]->Fill(pt,m) ;
2506
2507             if(fCalorimeter=="EMCAL"){
2508
2509               // Same sector
2510               Int_t j=0;
2511               for(Int_t i = 0; i < fNModules/2; i++){
2512                 j=2*i;
2513                 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
2514               }
2515
2516               // Same side
2517               for(Int_t i = 0; i < fNModules-2; i++){
2518                 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
2519               }
2520             }//EMCAL
2521             else {//PHOS
2522               if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2523               if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2524               if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
2525             }//PHOS
2526
2527
2528           }
2529
2530           Bool_t ok = kTRUE;
2531           if(fSameSM && module1!=module2) ok=kFALSE;
2532           if(ok){
2533
2534             //Check if one of the clusters comes from a conversion
2535             if(fCheckConversion){
2536               if     (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2537               else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2538             }
2539             //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2540             for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2541               if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
2542               {
2543                 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2544                 {
2545                   if(a < fAsymCuts[iasym])
2546                   {
2547                     Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2548
2549                     if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2550
2551                     fhMi1     [index]->Fill(pt,m) ;
2552                     if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
2553                     if(fFillBadDistHisto)
2554                     {
2555                       if(p1->DistToBad()>0 && p2->DistToBad()>0)
2556                       {
2557                         fhMi2     [index]->Fill(pt,m) ;
2558                         if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
2559                         if(p1->DistToBad()>1 && p2->DistToBad()>1)
2560                         {
2561                           fhMi3     [index]->Fill(pt,m) ;
2562                           if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
2563                         }
2564                       }
2565                     }// Fill bad dist histo
2566                   }//Asymmetry cut
2567                 }// Asymmetry loop
2568               }//PID cut
2569             }//loop for histograms
2570
2571             //-----------------------
2572             //Multi cuts analysis
2573             //-----------------------
2574             if(fMultiCutAna){
2575               //Several pt,ncell and asymmetry cuts
2576
2577               for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2578                 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2579                   for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2580                     Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2581                     if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        &&
2582                        a        <   fAsymCuts[iasym]                                    //&&
2583                        //p1->GetBtag() >=  fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2584                        ){
2585                       fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
2586                       //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2587                     }
2588                   }// pid bit cut loop
2589                 }// icell loop
2590               }// pt cut loop
2591             } // Multi cut ana
2592
2593             //Fill histograms with opening angle
2594             if(fFillAngleHisto){
2595               fhMixedOpeningAngle   ->Fill(pt,angle);
2596               fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2597             }
2598
2599           }//ok
2600         }// second cluster loop
2601       }//first cluster loop
2602     }//loop on mixed events
2603
2604     //--------------------------------------------------------
2605     //Add the current event to the list of events for mixing
2606     //--------------------------------------------------------
2607     TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2608     //Add current event to buffer and Remove redundant events
2609     if(currentEvent->GetEntriesFast()>0){
2610       evMixList->AddFirst(currentEvent) ;
2611       currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2612       if(evMixList->GetSize() >= GetNMaxEvMix())
2613       {
2614         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2615         evMixList->RemoveLast() ;
2616         delete tmp ;
2617       }
2618     }
2619     else{ //empty event
2620       delete currentEvent ;
2621       currentEvent=0 ;
2622     }
2623   }// DoOwnMix
2624
2625 }
2626
2627 //________________________________________________________________________
2628 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2629 {
2630   // retieves the event index and checks the vertex
2631   //    in the mixed buffer returns -2 if vertex NOK
2632   //    for normal events   returns 0 if vertex OK and -1 if vertex NOK
2633
2634   Int_t evtIndex = -1 ;
2635   if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2636
2637     if (GetMixedEvent()){
2638
2639       evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2640       GetVertex(vert,evtIndex);
2641
2642       if(TMath::Abs(vert[2])> GetZvertexCut())
2643         evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2644     } else {// Single event
2645
2646       GetVertex(vert);
2647
2648       if(TMath::Abs(vert[2])> GetZvertexCut())
2649         evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2650       else
2651         evtIndex = 0 ;
2652     }
2653   }//No MC reader
2654   else {
2655     evtIndex = 0;
2656     vert[0] = 0. ;
2657     vert[1] = 0. ;
2658     vert[2] = 0. ;
2659   }
2660
2661   return evtIndex ;
2662 }
2663