PWGGA/CaloTrackCorrelations/AliAnaPi0.cxx

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