PWGGA/CaloTrackCorrelations/AliAnaPi0.cxx

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