PWGGA/CaloTrackCorrelations/AliAnaPi0.cxx

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