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 )
1689   {
1690     if(ancLabel >= 0) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1691                              ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1692     else              printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor not found \n");
1693   }
1694
1695   Float_t prodR = -1;
1696   Int_t mcIndex = -1;
1697
1698   if(ancLabel > -1)
1699   {
1700     if(ancPDG==22)
1701     {//gamma
1702       mcIndex = 0;
1703     }
1704     else if(TMath::Abs(ancPDG)==11)
1705     {//e
1706       mcIndex = 1;
1707     }
1708     else if(ancPDG==111)
1709     {//Pi0
1710       mcIndex = 2;
1711       if(fMultiCutAnaSim)
1712       {
1713         for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1714         {
1715           for(Int_t icell=0; icell<fNCellNCuts; icell++)
1716           {
1717             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1718             {
1719               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1720               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1721                  asym   <  fAsymCuts[iasym]                                   &&
1722                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
1723               {
1724                 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1725                 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1726                 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1727               }//pass the different cuts
1728             }// pid bit cut loop
1729           }// icell loop
1730         }// pt cut loop
1731       }//Multi cut ana sim
1732       else
1733       {
1734         fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1735
1736         if(mass < 0.17 && mass > 0.1)
1737         {
1738           fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1739
1740           //Int_t uniqueId = -1;
1741           if(GetReader()->ReadStack())
1742           {
1743             TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1744             momindex  = ancestor->GetFirstMother();
1745             if(momindex < 0) return;
1746             TParticle* mother = GetMCStack()->Particle(momindex);
1747             mompdg    = TMath::Abs(mother->GetPdgCode());
1748             momstatus = mother->GetStatusCode();
1749             prodR = mother->R();
1750             //uniqueId = mother->GetUniqueID();
1751           }
1752           else
1753           {
1754             TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1755             AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1756             momindex  = ancestor->GetMother();
1757             if(momindex < 0) return;
1758             AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1759             mompdg    = TMath::Abs(mother->GetPdgCode());
1760             momstatus = mother->GetStatus();
1761             prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1762             //uniqueId = mother->GetUniqueID();
1763           }
1764
1765           //            printf("Reco Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1766           //                   pt,prodR,mompdg,momstatus,uniqueId);
1767
1768           fhMCPi0ProdVertex->Fill(pt,prodR);
1769
1770           if     (momstatus  == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1771           else if(mompdg     < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1772           else if(mompdg     > 2100  && mompdg   < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1773           else if(mompdg    == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1774           else if(mompdg    == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1775           else if(mompdg    == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1776           else if(mompdg    == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1777           else if(mompdg    >= 310   && mompdg    <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1778           else if(mompdg    == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1779           else if(momstatus == 11 || momstatus  == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1780           else                      fhMCPi0PtOrigin->Fill(pt,7.5);//other?
1781
1782
1783         }//pi0 mass region
1784       }
1785     }
1786     else if(ancPDG==221)
1787     {//Eta
1788       mcIndex = 3;
1789       if(fMultiCutAnaSim)
1790       {
1791         for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
1792         {
1793           for(Int_t icell=0; icell<fNCellNCuts; icell++)
1794           {
1795             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
1796             {
1797               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1798               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1799                  asym   <  fAsymCuts[iasym]                                   &&
1800                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell])
1801               {
1802                 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1803                 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1804                 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1805               }//pass the different cuts
1806             }// pid bit cut loop
1807           }// icell loop
1808         }// pt cut loop
1809       } //Multi cut ana sim
1810       else
1811       {
1812         fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1813         if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1814
1815         if(GetReader()->ReadStack())
1816         {
1817           TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1818           momindex  = ancestor->GetFirstMother();
1819           if(momindex < 0) return;
1820           TParticle* mother = GetMCStack()->Particle(momindex);
1821           mompdg    = TMath::Abs(mother->GetPdgCode());
1822           momstatus = mother->GetStatusCode();
1823           prodR = mother->R();
1824         }
1825         else
1826         {
1827           TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1828           AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1829           momindex  = ancestor->GetMother();
1830           if(momindex < 0) return;
1831           AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1832           mompdg    = TMath::Abs(mother->GetPdgCode());
1833           momstatus = mother->GetStatus();
1834           prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1835         }
1836
1837         fhMCEtaProdVertex->Fill(pt,prodR);
1838
1839         if     (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1840         else if(mompdg    < 22  ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1841         else if(mompdg    > 2100  && mompdg   < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1842         else if(mompdg    == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1843         else if(momstatus == 11 || momstatus  == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1844         else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1845         //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1846
1847       }// eta mass region
1848     }
1849     else if(ancPDG==-2212){//AProton
1850       mcIndex = 4;
1851     }
1852     else if(ancPDG==-2112){//ANeutron
1853       mcIndex = 5;
1854     }
1855     else if(TMath::Abs(ancPDG)==13){//muons
1856       mcIndex = 6;
1857     }
1858     else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7)
1859     {
1860       if(ancStatus==1)
1861       {//Stable particles, converted? not decayed resonances
1862         mcIndex = 6;
1863       }
1864       else
1865       {//resonances and other decays, more hadron conversions?
1866         mcIndex = 7;
1867       }
1868     }
1869     else
1870     {//Partons, colliding protons, strings, intermediate corrections
1871       if(ancStatus==11 || ancStatus==12)
1872       {//String fragmentation
1873         mcIndex = 8;
1874       }
1875       else if (ancStatus==21){
1876         if(ancLabel < 2)
1877         {//Colliding protons
1878           mcIndex = 11;
1879         }//colliding protons
1880         else if(ancLabel < 6)
1881         {//partonic initial states interactions
1882           mcIndex = 9;
1883         }
1884         else if(ancLabel < 8)
1885         {//Final state partons radiations?
1886           mcIndex = 10;
1887         }
1888         // else {
1889         //   printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1890         //          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1891         // }
1892       }//status 21
1893       //else {
1894       //  printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1895       //         ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1896       // }
1897     }////Partons, colliding protons, strings, intermediate corrections
1898   }//ancLabel > -1
1899   else { //ancLabel <= -1
1900     //printf("Not related at all label = %d\n",ancLabel);
1901     mcIndex = 12;
1902   }
1903
1904   if(mcIndex >=0 && mcIndex < 13)
1905   {
1906     fhMCOrgMass[mcIndex]->Fill(pt,mass);
1907     fhMCOrgAsym[mcIndex]->Fill(pt,asym);
1908     fhMCOrgDeltaEta[mcIndex]->Fill(pt,deta);
1909     fhMCOrgDeltaPhi[mcIndex]->Fill(pt,dphi);
1910   }
1911
1912 }
1913
1914 //__________________________________________
1915 void AliAnaPi0::MakeAnalysisFillHistograms()
1916 {
1917   //Process one event and extract photons from AOD branch
1918   // filled with AliAnaPhoton and fill histos with invariant mass
1919
1920   //In case of simulated data, fill acceptance histograms
1921   if(IsDataMC())FillAcceptanceHistograms();
1922
1923   //if (GetReader()->GetEventNumber()%10000 == 0)
1924   // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1925
1926   if(!GetInputAODBranch())
1927   {
1928     printf("AliAnaPi0::MakeAnalysisFillHistograms() - No input aod photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1929     abort();
1930   }
1931
1932   //Init some variables
1933   Int_t   nPhot    = GetInputAODBranch()->GetEntriesFast() ;
1934
1935   if(GetDebug() > 1)
1936     printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
1937
1938   //If less than photon 2 entries in the list, skip this event
1939   if(nPhot < 2 )
1940   {
1941     if(GetDebug() > 2)
1942       printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
1943
1944     if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
1945
1946     return ;
1947   }
1948
1949   Int_t ncentr = GetNCentrBin();
1950   if(GetNCentrBin()==0) ncentr = 1;
1951
1952   //Init variables
1953   Int_t module1         = -1;
1954   Int_t module2         = -1;
1955   Double_t vert[]       = {0.0, 0.0, 0.0} ; //vertex
1956   Int_t evtIndex1       = 0 ;
1957   Int_t currentEvtIndex = -1;
1958   Int_t curCentrBin     = GetEventCentralityBin();
1959   //Int_t curVzBin        = GetEventVzBin();
1960   //Int_t curRPBin        = GetEventRPBin();
1961   Int_t eventbin        = GetEventMixBin();
1962
1963   if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
1964   {
1965      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());
1966     return;
1967   }
1968
1969   //Get shower shape information of clusters
1970   TObjArray *clusters = 0;
1971   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1972   else if(fCalorimeter=="PHOS" ) clusters = GetPHOSClusters() ;
1973
1974   //---------------------------------
1975   //First loop on photons/clusters
1976   //---------------------------------
1977   for(Int_t i1=0; i1<nPhot-1; i1++)
1978   {
1979     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
1980     //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d/%d\n",i1,nPhot-1);
1981
1982     // get the event index in the mixed buffer where the photon comes from
1983     // in case of mixing with analysis frame, not own mixing
1984     evtIndex1 = GetEventIndex(p1, vert) ;
1985     if ( evtIndex1 == -1 )
1986       return ;
1987     if ( evtIndex1 == -2 )
1988       continue ;
1989
1990     // Only effective in case of mixed event frame
1991     if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ;   //vertex cut
1992
1993     if (evtIndex1 != currentEvtIndex)
1994     {
1995       //Fill event bin info
1996       if(DoOwnMix()) fhEventBin->Fill(eventbin) ;
1997       if(GetNCentrBin() > 1)
1998       {
1999         fhCentrality->Fill(curCentrBin);
2000         if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
2001       }
2002       currentEvtIndex = evtIndex1 ;
2003     }
2004
2005     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
2006
2007     //Get the momentum of this cluster
2008     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2009
2010     //Get (Super)Module number of this cluster
2011     module1 = GetModuleNumber(p1);
2012
2013     //------------------------------------------
2014     // Recover original cluster
2015     Int_t iclus1 = -1 ;
2016     AliVCluster * cluster1 = FindCluster(clusters,p1->GetCaloLabel(0),iclus1);
2017     if(!cluster1) printf("AliAnaPi0 - Cluster1 not found!\n");
2018
2019     //---------------------------------
2020     //Second loop on photons/clusters
2021     //---------------------------------
2022     for(Int_t i2=i1+1; i2<nPhot; i2++)
2023     {
2024       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
2025       //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster2 i %d/%d\n",i2,nPhot);
2026
2027       //In case of mixing frame, check we are not in the same event as the first cluster
2028       Int_t evtIndex2 = GetEventIndex(p2, vert) ;
2029       if ( evtIndex2 == -1 )
2030         return ;
2031       if ( evtIndex2 == -2 )
2032         continue ;
2033       if (GetMixedEvent() && (evtIndex1 == evtIndex2))
2034         continue ;
2035
2036       //------------------------------------------
2037       // Recover original cluster
2038       Int_t iclus2 = -1;
2039       AliVCluster * cluster2 = FindCluster(clusters,p2->GetCaloLabel(0),iclus2,iclus1+1);
2040       // start new loop from iclus1+1 to gain some time
2041       if(!cluster2) printf("AliAnaPi0 - Cluster2 not found!\n");
2042
2043       // Get the TOF,l0 and ncells from the clusters
2044       Float_t tof1  = -1;
2045       Float_t l01   = -1;
2046       Int_t ncell1  = 0;
2047       if(cluster1)
2048       {
2049         tof1   = cluster1->GetTOF()*1e9;
2050         l01    = cluster1->GetM02();
2051         ncell1 = cluster1->GetNCells();
2052         //printf("cluster1: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster1->E(),p1->E(),l01,tof1);
2053       }
2054       //else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
2055       //            p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
2056
2057       Float_t tof2  = -1;
2058       Float_t l02   = -1;
2059       Int_t ncell2  = 0;
2060       if(cluster2)
2061       {
2062         tof2   = cluster2->GetTOF()*1e9;
2063         l02    = cluster2->GetM02();
2064         ncell2 = cluster2->GetNCells();
2065         //printf("cluster2: E %2.2f (%2.2f), l0 %2.2f, tof %2.2f\n",cluster2->E(),p2->E(),l02,tof2);
2066       }
2067       //else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
2068       //            p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
2069
2070       if(cluster1 && cluster2)
2071       {
2072         Double_t t12diff = tof1-tof2;
2073         if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2074       }
2075       //------------------------------------------
2076
2077       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
2078
2079       //Get the momentum of this cluster
2080       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2081       //Get module number
2082       module2       = GetModuleNumber(p2);
2083
2084       //---------------------------------
2085       // Get pair kinematics
2086       //---------------------------------
2087       Double_t m    = (photon1 + photon2).M() ;
2088       Double_t pt   = (photon1 + photon2).Pt();
2089       Double_t deta = photon1.Eta() - photon2.Eta();
2090       Double_t dphi = photon1.Phi() - photon2.Phi();
2091       Double_t a    = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2092
2093       if(GetDebug() > 2)
2094         printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
2095
2096       //--------------------------------
2097       // Opening angle selection
2098       //--------------------------------
2099       //Check if opening angle is too large or too small compared to what is expected
2100       Double_t angle   = photon1.Angle(photon2.Vect());
2101       if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05))
2102       {
2103         if(GetDebug() > 2)
2104           printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2105         continue;
2106       }
2107
2108       if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
2109       {
2110         if(GetDebug() > 2)
2111           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
2112         continue;
2113       }
2114
2115       //-------------------------------------------------------------------------------------------------
2116       //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2117       //-------------------------------------------------------------------------------------------------
2118       if(a < fAsymCuts[0] && fFillSMCombinations)
2119       {
2120         if(module1==module2 && module1 >=0 && module1<fNModules)
2121           fhReMod[module1]->Fill(pt,m) ;
2122
2123         if(fCalorimeter=="EMCAL")
2124         {
2125           // Same sector
2126           Int_t j=0;
2127           for(Int_t i = 0; i < fNModules/2; i++)
2128           {
2129             j=2*i;
2130             if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
2131           }
2132
2133           // Same side
2134           for(Int_t i = 0; i < fNModules-2; i++){
2135             if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
2136           }
2137         }//EMCAL
2138         else {//PHOS
2139           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
2140           if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
2141           if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
2142         }//PHOS
2143       }
2144
2145       //In case we want only pairs in same (super) module, check their origin.
2146       Bool_t ok = kTRUE;
2147       if(fSameSM && module1!=module2) ok=kFALSE;
2148       if(ok)
2149       {
2150         //Check if one of the clusters comes from a conversion
2151         if(fCheckConversion)
2152         {
2153           if     (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
2154           else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
2155         }
2156
2157         // Fill shower shape cut histograms
2158         if(fFillSSCombinations)
2159         {
2160           if     ( l01 > 0.01 && l01 < 0.4  &&
2161                    l02 > 0.01 && l02 < 0.4 )               fhReSS[0]->Fill(pt,m); // Tight
2162           else if( l01 > 0.4  && l02 > 0.4 )               fhReSS[1]->Fill(pt,m); // Loose
2163           else if( l01 > 0.01 && l01 < 0.4  && l02 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2164           else if( l02 > 0.01 && l02 < 0.4  && l01 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2165         }
2166
2167         //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2168         for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2169         {
2170           if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)))
2171           {
2172             for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2173             {
2174               if(a < fAsymCuts[iasym])
2175               {
2176                 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2177                 //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);
2178
2179                 if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2180
2181                 fhRe1     [index]->Fill(pt,m);
2182                 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
2183                 if(fFillBadDistHisto)
2184                 {
2185                   if(p1->DistToBad()>0 && p2->DistToBad()>0)
2186                   {
2187                     fhRe2     [index]->Fill(pt,m) ;
2188                     if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
2189                     if(p1->DistToBad()>1 && p2->DistToBad()>1)
2190                     {
2191                       fhRe3     [index]->Fill(pt,m) ;
2192                       if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
2193                     }// bad 3
2194                   }// bad2
2195                 }// Fill bad dist histos
2196               }//assymetry cut
2197             }// asymmetry cut loop
2198           }// bad 1
2199         }// pid bit loop
2200
2201         //Fill histograms with opening angle
2202         if(fFillAngleHisto)
2203         {
2204           fhRealOpeningAngle   ->Fill(pt,angle);
2205           fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2206         }
2207
2208         //Fill histograms with pair assymmetry
2209         if(fFillAsymmetryHisto)
2210         {
2211           fhRePtAsym->Fill(pt,a);
2212           if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
2213           if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
2214         }
2215
2216         //---------
2217         // MC data
2218         //---------
2219         //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2220         if(IsDataMC())
2221         {
2222           if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2223              GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2224           {
2225             fhReMCFromConversion->Fill(pt,m);
2226           }
2227           else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2228                   !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2229           {
2230             fhReMCFromNotConversion->Fill(pt,m);
2231           }
2232           else
2233           {
2234             fhReMCFromMixConversion->Fill(pt,m);
2235           }
2236
2237           if(fFillOriginHisto)
2238             FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2239         }
2240
2241         //-----------------------
2242         //Multi cuts analysis
2243         //-----------------------
2244         if(fMultiCutAna)
2245         {
2246           //Histograms for different PID bits selection
2247           for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2248           {
2249             if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    &&
2250                p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
2251
2252             //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2253           } // pid bit cut loop
2254
2255           //Several pt,ncell and asymmetry cuts
2256           for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2257           {
2258             for(Int_t icell=0; icell<fNCellNCuts; icell++)
2259             {
2260               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2261               {
2262                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2263                 if(p1->E() >   fPtCuts[ipt]      && p2->E() > fPtCuts[ipt]        &&
2264                    a        <   fAsymCuts[iasym]                                  &&
2265                    ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell])
2266                 {
2267                   fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
2268                   //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2269                   if(fFillSMCombinations && module1==module2)
2270                   {
2271                     fhRePtNCellAsymCutsSM[module1][index]->Fill(pt,m) ;
2272                   }
2273                 }
2274               }// pid bit cut loop
2275             }// icell loop
2276           }// pt cut loop
2277
2278           if(GetHistogramRanges()->GetHistoTrackMultiplicityBins())
2279           {
2280             for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++)
2281             {
2282               if(a < fAsymCuts[iasym]) fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
2283             }
2284           }
2285         }// multiple cuts analysis
2286
2287       }// ok if same sm
2288
2289     }// second same event particle
2290   }// first cluster
2291
2292   //-------------------------------------------------------------
2293   // Mixing
2294   //-------------------------------------------------------------
2295   if(DoOwnMix())
2296   {
2297     //Recover events in with same characteristics as the current event
2298
2299     //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2300     if(eventbin < 0) return ;
2301
2302     TList * evMixList=fEventsList[eventbin] ;
2303
2304     if(!evMixList)
2305     {
2306       printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mix event list not available, bin %d \n",eventbin);
2307       return;
2308     }
2309
2310     Int_t nMixed = evMixList->GetSize() ;
2311     for(Int_t ii=0; ii<nMixed; ii++)
2312     {
2313       TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2314       Int_t nPhot2=ev2->GetEntriesFast() ;
2315       Double_t m = -999;
2316       if(GetDebug() > 1)
2317         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n",
2318                ii, nPhot2, GetEventCentralityBin());
2319
2320       fhEventMixBin->Fill(eventbin) ;
2321
2322       //---------------------------------
2323       //First loop on photons/clusters
2324       //---------------------------------
2325       for(Int_t i1=0; i1<nPhot; i1++)
2326       {
2327         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2328
2329         if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2330
2331         //Get kinematics of cluster and (super) module of this cluster
2332         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2333         module1 = GetModuleNumber(p1);
2334
2335         //---------------------------------
2336         //First loop on photons/clusters
2337         //---------------------------------
2338         for(Int_t i2=0; i2<nPhot2; i2++)
2339         {
2340           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2341
2342           //Get kinematics of second cluster and calculate those of the pair
2343           TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2344           m           = (photon1+photon2).M() ;
2345           Double_t pt = (photon1 + photon2).Pt();
2346           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2347
2348           //Check if opening angle is too large or too small compared to what is expected
2349           Double_t angle   = photon1.Angle(photon2.Vect());
2350           if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05))
2351           {
2352             if(GetDebug() > 2)
2353               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2354             continue;
2355           }
2356
2357           if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut))
2358           {
2359             if(GetDebug() > 2)
2360               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2361             continue;
2362           }
2363
2364           if(GetDebug() > 2)
2365             printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
2366                    p1->Pt(), p2->Pt(), pt,m,a);
2367
2368           //In case we want only pairs in same (super) module, check their origin.
2369           module2 = GetModuleNumber(p2);
2370
2371           //-------------------------------------------------------------------------------------------------
2372           //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2373           //-------------------------------------------------------------------------------------------------
2374           if(a < fAsymCuts[0] && fFillSMCombinations){
2375             if(module1==module2 && module1 >=0 && module1<fNModules)
2376               fhMiMod[module1]->Fill(pt,m) ;
2377
2378             if(fCalorimeter=="EMCAL")
2379             {
2380               // Same sector
2381               Int_t j=0;
2382               for(Int_t i = 0; i < fNModules/2; i++)
2383               {
2384                 j=2*i;
2385                 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
2386               }
2387
2388               // Same side
2389               for(Int_t i = 0; i < fNModules-2; i++)
2390               {
2391                 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
2392               }
2393             }//EMCAL
2394             else
2395             {//PHOS
2396               if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2397               if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2398               if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
2399             }//PHOS
2400
2401
2402           }
2403
2404           Bool_t ok = kTRUE;
2405           if(fSameSM && module1!=module2) ok=kFALSE;
2406           if(ok){
2407
2408             //Check if one of the clusters comes from a conversion
2409             if(fCheckConversion)
2410             {
2411               if     (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2412               else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2413             }
2414             //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2415             for(Int_t ipid=0; ipid<fNPIDBits; ipid++)
2416             {
2417               if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
2418               {
2419                 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2420                 {
2421                   if(a < fAsymCuts[iasym])
2422                   {
2423                     Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2424
2425                     if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2426
2427                     fhMi1     [index]->Fill(pt,m) ;
2428
2429                     if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
2430
2431                     if(fFillBadDistHisto)
2432                     {
2433                       if(p1->DistToBad()>0 && p2->DistToBad()>0)
2434                       {
2435                         fhMi2     [index]->Fill(pt,m) ;
2436                         if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
2437                         if(p1->DistToBad()>1 && p2->DistToBad()>1)
2438                         {
2439                           fhMi3     [index]->Fill(pt,m) ;
2440                           if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
2441                         }
2442                       }
2443                     }// Fill bad dist histo
2444                   }//Asymmetry cut
2445                 }// Asymmetry loop
2446               }//PID cut
2447             }//loop for histograms
2448
2449             //-----------------------
2450             //Multi cuts analysis
2451             //-----------------------
2452             if(fMultiCutAna){
2453               //Several pt,ncell and asymmetry cuts
2454
2455               for(Int_t ipt=0; ipt<fNPtCuts; ipt++)
2456               {
2457                 for(Int_t icell=0; icell<fNCellNCuts; icell++)
2458                 {
2459                   for(Int_t iasym=0; iasym<fNAsymCuts; iasym++)
2460                   {
2461                     Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2462                     if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        &&
2463                        a        <   fAsymCuts[iasym]                                    //&&
2464                        //p1->GetBtag() >=  fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2465                        )
2466                     {
2467                       fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
2468                       //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2469                     }
2470                   }// pid bit cut loop
2471                 }// icell loop
2472               }// pt cut loop
2473             } // Multi cut ana
2474
2475             //Fill histograms with opening angle
2476             if(fFillAngleHisto)
2477             {
2478               fhMixedOpeningAngle   ->Fill(pt,angle);
2479               fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2480             }
2481
2482           }//ok
2483         }// second cluster loop
2484       }//first cluster loop
2485     }//loop on mixed events
2486
2487     //--------------------------------------------------------
2488     // Add the current event to the list of events for mixing
2489     //--------------------------------------------------------
2490
2491     TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2492     //Add current event to buffer and Remove redundant events
2493     if( currentEvent->GetEntriesFast() > 0 )
2494     {
2495       evMixList->AddFirst(currentEvent) ;
2496       currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2497       if( evMixList->GetSize() >= GetNMaxEvMix() )
2498       {
2499         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2500         evMixList->RemoveLast() ;
2501         delete tmp ;
2502       }
2503     }
2504     else
2505     { //empty event
2506       delete currentEvent ;
2507       currentEvent=0 ;
2508     }
2509   }// DoOwnMix
2510
2511   if(GetDebug() > 0) printf("AliAnaPi0::MakeAnalysisFillHistograms() - End fill histograms\n");
2512 }
2513
2514 //________________________________________________________________________
2515 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2516 {
2517   // retieves the event index and checks the vertex
2518   //    in the mixed buffer returns -2 if vertex NOK
2519   //    for normal events   returns 0 if vertex OK and -1 if vertex NOK
2520
2521   Int_t evtIndex = -1 ;
2522   if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2523   {
2524     if (GetMixedEvent())
2525     {
2526       evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2527       GetVertex(vert,evtIndex);
2528
2529       if(TMath::Abs(vert[2])> GetZvertexCut())
2530         evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2531     }
2532     else
2533     {
2534       // Single event
2535       GetVertex(vert);
2536
2537       if(TMath::Abs(vert[2])> GetZvertexCut())
2538         evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2539       else
2540         evtIndex = 0 ;
2541     }
2542   }//No MC reader
2543   else
2544   {
2545     evtIndex = 0;
2546     vert[0] = 0. ;
2547     vert[1] = 0. ;
2548     vert[2] = 0. ;
2549   }
2550
2551   return evtIndex ;
2552 }
2553