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