PWGGA/CaloTrackCorrelations/AliAnaPi0.cxx

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