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