PWGGA/CaloTrackCorrelations/AliAnaPi0.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //_________________________________________________________________________
  17 // Class to collect two-photon invariant mass distributions for
  18 // extracting raw pi0 yield.
  19 // Input is produced by AliAnaPhoton (or any other analysis producing output AliAODPWG4Particles),
  20 // it will do nothing if executed alone
  21 //
  22 //-- Author: Dmitri Peressounko (RRC "KI")
  23 //-- Adapted to CaloTrackCorr frame by Lamia Benhabib (SUBATECH)
  24 //-- and Gustavo Conesa (INFN-Frascati)
  25 //_________________________________________________________________________
  26
  27
  28 // --- ROOT system ---
  29 #include "TH3.h"
  30 #include "TH2F.h"
  31 //#include "Riostream.h"
  32 #include "TCanvas.h"
  33 #include "TPad.h"
  34 #include "TROOT.h"
  35 #include "TClonesArray.h"
  36 #include "TObjString.h"
  37 #include "TDatabasePDG.h"
  38
  39 //---- AliRoot system ----
  40 #include "AliAnaPi0.h"
  41 #include "AliCaloTrackReader.h"
  42 #include "AliCaloPID.h"
  43 #include "AliStack.h"
  44 #include "AliFiducialCut.h"
  45 #include "TParticle.h"
  46 #include "AliVEvent.h"
  47 #include "AliESDCaloCluster.h"
  48 #include "AliESDEvent.h"
  49 #include "AliAODEvent.h"
  50 #include "AliNeutralMesonSelection.h"
  51 #include "AliMixedEvent.h"
  52 #include "AliAODMCParticle.h"
  53
  54 // --- Detectors ---
  55 #include "AliPHOSGeoUtils.h"
  56 #include "AliEMCALGeometry.h"
  57
  58 ClassImp(AliAnaPi0)
  59
  60 //______________________________________________________
  61 AliAnaPi0::AliAnaPi0() : AliAnaCaloTrackCorrBaseClass(),
  62 fEventsList(0x0),
  63 fCalorimeter(""),            fNModules(22),
  64 fUseAngleCut(kFALSE),        fUseAngleEDepCut(kFALSE),     fAngleCut(0),                 fAngleMaxCut(7.),
  65 fMultiCutAna(kFALSE),        fMultiCutAnaSim(kFALSE),
  66 fNPtCuts(0),                 fNAsymCuts(0),                fNCellNCuts(0),               fNPIDBits(0),
  67 fMakeInvPtPlots(kFALSE),     fSameSM(kFALSE),
  68 fFillSMCombinations(kFALSE), fCheckConversion(kFALSE),
  69 fFillBadDistHisto(kFALSE),   fFillSSCombinations(kFALSE),
  70 fFillAngleHisto(kFALSE),     fFillAsymmetryHisto(kFALSE),  fFillOriginHisto(0),          fFillArmenterosThetaStar(0),
  71 fCheckAccInSector(kFALSE),
  72 //Histograms
  73 fhAverTotECluster(0),        fhAverTotECell(0),            fhAverTotECellvsCluster(0),
  74 fhEDensityCluster(0),        fhEDensityCell(0),            fhEDensityCellvsCluster(0),
  75 fhReMod(0x0),                fhReSameSideEMCALMod(0x0),    fhReSameSectorEMCALMod(0x0),  fhReDiffPHOSMod(0x0),
  76 fhMiMod(0x0),                fhMiSameSideEMCALMod(0x0),    fhMiSameSectorEMCALMod(0x0),  fhMiDiffPHOSMod(0x0),
  77 fhReConv(0x0),               fhMiConv(0x0),                fhReConv2(0x0),  fhMiConv2(0x0),
  78 fhRe1(0x0),                  fhMi1(0x0),                   fhRe2(0x0),                   fhMi2(0x0),
  79 fhRe3(0x0),                  fhMi3(0x0),                   fhReInvPt1(0x0),              fhMiInvPt1(0x0),
  80 fhReInvPt2(0x0),             fhMiInvPt2(0x0),              fhReInvPt3(0x0),              fhMiInvPt3(0x0),
  81 fhRePtNCellAsymCuts(0x0),    fhMiPtNCellAsymCuts(0x0),     fhRePtNCellAsymCutsSM(),
  82 fhRePIDBits(0x0),            fhRePtMult(0x0),              fhReSS(),
  83 fhRePtAsym(0x0),             fhRePtAsymPi0(0x0),           fhRePtAsymEta(0x0),
  84 fhEventBin(0),               fhEventMixBin(0),
  85 fhCentrality(0x0),           fhCentralityNoPair(0x0),
  86 fhEventPlaneResolution(0x0),
  87 fhRealOpeningAngle(0x0),     fhRealCosOpeningAngle(0x0),   fhMixedOpeningAngle(0x0),     fhMixedCosOpeningAngle(0x0),
  88 // MC histograms
  89 fhPrimPi0E(0x0),             fhPrimPi0Pt(0x0),             fhPrimPi0PtRejected(0x0),
  90 fhPrimPi0AccE(0x0),          fhPrimPi0AccPt(0x0),
  91 fhPrimPi0Y(0x0),             fhPrimPi0AccY(0x0),
  92 fhPrimPi0Yeta(0x0),          fhPrimPi0YetaYcut(0x0),       fhPrimPi0AccYeta(0x0),
  93 fhPrimPi0Phi(0x0),           fhPrimPi0AccPhi(0x0),
  94 fhPrimPi0OpeningAngle(0x0),  fhPrimPi0OpeningAngleAsym(0x0),fhPrimPi0CosOpeningAngle(0x0),
  95 fhPrimPi0PtCentrality(0),    fhPrimPi0PtEventPlane(0),
  96 fhPrimPi0AccPtCentrality(0), fhPrimPi0AccPtEventPlane(0),
  97 fhPrimEtaE(0x0),             fhPrimEtaPt(0x0),             fhPrimEtaPtRejected(0x0),
  98 fhPrimEtaAccE(0x0),          fhPrimEtaAccPt(0x0),
  99 fhPrimEtaY(0x0),             fhPrimEtaAccY(0x0),
 100 fhPrimEtaYeta(0x0),          fhPrimEtaYetaYcut(0x0),       fhPrimEtaAccYeta(0x0),
 101 fhPrimEtaPhi(0x0),           fhPrimEtaAccPhi(0x0),
 102 fhPrimEtaOpeningAngle(0x0),  fhPrimEtaOpeningAngleAsym(0x0),fhPrimEtaCosOpeningAngle(0x0),
 103 fhPrimEtaPtCentrality(0),    fhPrimEtaPtEventPlane(0),
 104 fhPrimEtaAccPtCentrality(0), fhPrimEtaAccPtEventPlane(0),
 105 fhPrimPi0PtOrigin(0x0),      fhPrimEtaPtOrigin(0x0),
 106 fhMCOrgMass(),               fhMCOrgAsym(),                fhMCOrgDeltaEta(),            fhMCOrgDeltaPhi(),
 107 fhMCPi0MassPtRec(),          fhMCPi0MassPtTrue(),          fhMCPi0PtTruePtRec(),
 108 fhMCEtaMassPtRec(),          fhMCEtaMassPtTrue(),          fhMCEtaPtTruePtRec(),
 109 fhMCPi0PtOrigin(0x0),        fhMCEtaPtOrigin(0x0),
 110 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     fhPrimPi0PtRejected = new TH1F("hPrimPi0PtRejected","Primary #pi^{0} pt",nptbins,ptmin,ptmax) ;
 699     fhPrimPi0Pt     = new TH1F("hPrimPi0Pt","Primary #pi^{0} #it{p}_{T} , |#it{Y}|<1",nptbins,ptmin,ptmax) ;
 700     fhPrimPi0AccPt  = new TH1F("hPrimPi0AccPt","Primary #pi^{0} #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
 701     fhPrimPi0Pt   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 702     fhPrimPi0AccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 703     outputContainer->Add(fhPrimPi0PtRejected) ;
 704     outputContainer->Add(fhPrimPi0Pt) ;
 705     outputContainer->Add(fhPrimPi0AccPt) ;
 706
 707     Int_t netabinsopen =  TMath::Nint(netabins*4/(etamax-etamin));
 708     fhPrimPi0Y      = new TH2F("hPrimPi0Rapidity","Rapidity of primary #pi^{0}",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
 709     fhPrimPi0Y   ->SetYTitle("#it{Rapidity}");
 710     fhPrimPi0Y   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 711     outputContainer->Add(fhPrimPi0Y) ;
 712
 713     fhPrimPi0Yeta      = new TH2F("hPrimPi0PseudoRapidity","PseudoRapidity of primary #pi^{0}",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
 714     fhPrimPi0Yeta   ->SetYTitle("#eta");
 715     fhPrimPi0Yeta   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 716     outputContainer->Add(fhPrimPi0Yeta) ;
 717
 718     fhPrimPi0YetaYcut      = new TH2F("hPrimPi0PseudoRapidityYcut","PseudoRapidity of primary #pi^{0}, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
 719     fhPrimPi0YetaYcut   ->SetYTitle("#eta");
 720     fhPrimPi0YetaYcut   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 721     outputContainer->Add(fhPrimPi0YetaYcut) ;
 722
 723     fhPrimPi0AccY   = new TH2F("hPrimPi0AccRapidity","Rapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
 724     fhPrimPi0AccY->SetYTitle("Rapidity");
 725     fhPrimPi0AccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 726     outputContainer->Add(fhPrimPi0AccY) ;
 727
 728     fhPrimPi0AccYeta      = new TH2F("hPrimPi0AccPseudoRapidity","PseudoRapidity of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,netabins,etamin,etamax) ;
 729     fhPrimPi0AccYeta   ->SetYTitle("#eta");
 730     fhPrimPi0AccYeta   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 731     outputContainer->Add(fhPrimPi0AccYeta) ;
 732
 733     Int_t nphibinsopen = TMath::Nint(nphibins*TMath::TwoPi()/(phimax-phimin));
 734     fhPrimPi0Phi    = new TH2F("hPrimPi0Phi","#phi of primary #pi^{0}, |#it{Y}|<1",nptbins,ptmin,ptmax,nphibinsopen,0,360) ;
 735     fhPrimPi0Phi->SetYTitle("#phi (deg)");
 736     fhPrimPi0Phi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 737     outputContainer->Add(fhPrimPi0Phi) ;
 738
 739     fhPrimPi0AccPhi = new TH2F("hPrimPi0AccPhi","#phi of primary #pi^{0} with accepted daughters",nptbins,ptmin,ptmax,
 740                                nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 741     fhPrimPi0AccPhi->SetYTitle("#phi (deg)");
 742     fhPrimPi0AccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 743     outputContainer->Add(fhPrimPi0AccPhi) ;
 744
 745     fhPrimPi0PtCentrality     = new TH2F("hPrimPi0PtCentrality","Primary #pi^{0} #it{p}_{T} vs reco centrality, |#it{Y}|<1",
 746                                          nptbins,ptmin,ptmax, 100, 0, 100) ;
 747     fhPrimPi0AccPtCentrality  = new TH2F("hPrimPi0AccPtCentrality","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco centrality",
 748                                          nptbins,ptmin,ptmax, 100, 0, 100) ;
 749     fhPrimPi0PtCentrality   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 750     fhPrimPi0AccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 751     fhPrimPi0PtCentrality   ->SetYTitle("Centrality");
 752     fhPrimPi0AccPtCentrality->SetYTitle("Centrality");
 753     outputContainer->Add(fhPrimPi0PtCentrality) ;
 754     outputContainer->Add(fhPrimPi0AccPtCentrality) ;
 755
 756     fhPrimPi0PtEventPlane     = new TH2F("hPrimPi0PtEventPlane","Primary #pi^{0} #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",
 757                                          nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
 758     fhPrimPi0AccPtEventPlane  = new TH2F("hPrimPi0AccPtEventPlane","Primary #pi^{0} with both photons in acceptance #it{p}_{T} vs reco event plane angle",
 759                                          nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
 760     fhPrimPi0PtEventPlane   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 761     fhPrimPi0AccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 762     fhPrimPi0PtEventPlane   ->SetYTitle("Event Plane Angle (rad)");
 763     fhPrimPi0AccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
 764     outputContainer->Add(fhPrimPi0PtEventPlane) ;
 765     outputContainer->Add(fhPrimPi0AccPtEventPlane) ;
 766
 767     //Eta
 768
 769     fhPrimEtaE     = new TH1F("hPrimEtaE","Primary eta E",nptbins,ptmin,ptmax) ;
 770     fhPrimEtaAccE  = new TH1F("hPrimEtaAccE","Primary #eta #it{E} with both photons in acceptance",nptbins,ptmin,ptmax) ;
 771     fhPrimEtaE   ->SetXTitle("#it{E} (GeV)");
 772     fhPrimEtaAccE->SetXTitle("#it{E} (GeV)");
 773     outputContainer->Add(fhPrimEtaE) ;
 774     outputContainer->Add(fhPrimEtaAccE) ;
 775
 776     fhPrimEtaPtRejected = new TH1F("hPrimEtaPtRejected","Primary #eta #it{p}_{T}",nptbins,ptmin,ptmax) ;
 777     fhPrimEtaPt     = new TH1F("hPrimEtaPt","Primary #eta #it{p}_{T}",nptbins,ptmin,ptmax) ;
 778     fhPrimEtaAccPt  = new TH1F("hPrimEtaAccPt","Primary eta #it{p}_{T} with both photons in acceptance",nptbins,ptmin,ptmax) ;
 779     fhPrimEtaPt   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 780     fhPrimEtaAccPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 781     outputContainer->Add(fhPrimEtaPtRejected) ;
 782     outputContainer->Add(fhPrimEtaPt) ;
 783     outputContainer->Add(fhPrimEtaAccPt) ;
 784
 785     fhPrimEtaY      = new TH2F("hPrimEtaRapidity","Rapidity of primary #eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
 786     fhPrimEtaY->SetYTitle("#it{Rapidity}");
 787     fhPrimEtaY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 788     outputContainer->Add(fhPrimEtaY) ;
 789
 790     fhPrimEtaYeta      = new TH2F("hPrimEtaPseudoRapidityEta","PseudoRapidity of primary #eta",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
 791     fhPrimEtaYeta->SetYTitle("#it{Rapidity}");
 792     fhPrimEtaYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 793     outputContainer->Add(fhPrimEtaYeta) ;
 794
 795     fhPrimEtaYetaYcut      = new TH2F("hPrimEtaPseudoRapidityEtaYcut","PseudoRapidity of primary #eta, |#it{Y}|<1",nptbins,ptmin,ptmax,netabinsopen,-2, 2) ;
 796     fhPrimEtaYetaYcut->SetYTitle("#it{Pseudorapidity}");
 797     fhPrimEtaYetaYcut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 798     outputContainer->Add(fhPrimEtaYetaYcut) ;
 799
 800     fhPrimEtaAccY   = new TH2F("hPrimEtaAccRapidity","Rapidity of primary #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
 801     fhPrimEtaAccY->SetYTitle("#it{Rapidity}");
 802     fhPrimEtaAccY->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 803     outputContainer->Add(fhPrimEtaAccY) ;
 804
 805     fhPrimEtaAccYeta  = new TH2F("hPrimEtaAccPseudoRapidity","PseudoRapidity of primary #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax) ;
 806     fhPrimEtaAccYeta->SetYTitle("#it{Pseudorapidity}");
 807     fhPrimEtaAccYeta->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 808     outputContainer->Add(fhPrimEtaAccYeta) ;
 809
 810     fhPrimEtaPhi    = new TH2F("hPrimEtaPhi","Azimuthal of primary #eta",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 811     fhPrimEtaPhi->SetYTitle("#phi (deg)");
 812     fhPrimEtaPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 813     outputContainer->Add(fhPrimEtaPhi) ;
 814
 815     fhPrimEtaAccPhi = new TH2F("hPrimEtaAccPhi","Azimuthal of primary #eta with accepted daughters",nptbins,ptmin,ptmax, nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 816     fhPrimEtaAccPhi->SetYTitle("#phi (deg)");
 817     fhPrimEtaAccPhi->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 818     outputContainer->Add(fhPrimEtaAccPhi) ;
 819
 820     fhPrimEtaPtCentrality     = new TH2F("hPrimEtaPtCentrality","Primary #eta #it{p}_{T} vs reco centrality, |#it{Y}|<1",nptbins,ptmin,ptmax, 100, 0, 100) ;
 821     fhPrimEtaAccPtCentrality  = new TH2F("hPrimEtaAccPtCentrality","Primary #eta with both photons in acceptance #it{p}_{T} vs reco centrality",nptbins,ptmin,ptmax, 100, 0, 100) ;
 822     fhPrimEtaPtCentrality   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 823     fhPrimEtaAccPtCentrality->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 824     fhPrimEtaPtCentrality   ->SetYTitle("Centrality");
 825     fhPrimEtaAccPtCentrality->SetYTitle("Centrality");
 826     outputContainer->Add(fhPrimEtaPtCentrality) ;
 827     outputContainer->Add(fhPrimEtaAccPtCentrality) ;
 828
 829     fhPrimEtaPtEventPlane     = new TH2F("hPrimEtaPtEventPlane","Primary #eta #it{p}_{T} vs reco event plane angle, |#it{Y}|<1",nptbins,ptmin,ptmax, 100, 0, TMath::Pi()) ;
 830     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()) ;
 831     fhPrimEtaPtEventPlane   ->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 832     fhPrimEtaAccPtEventPlane->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 833     fhPrimEtaPtEventPlane   ->SetYTitle("Event Plane Angle (rad)");
 834     fhPrimEtaAccPtEventPlane->SetYTitle("Event Plane Angle (rad)");
 835     outputContainer->Add(fhPrimEtaPtEventPlane) ;
 836     outputContainer->Add(fhPrimEtaAccPtEventPlane) ;
 837
 838      if(fFillAngleHisto)
 839     {
 840       fhPrimPi0OpeningAngle  = new TH2F
 841       ("hPrimPi0OpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
 842       fhPrimPi0OpeningAngle->SetYTitle("#theta(rad)");
 843       fhPrimPi0OpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 844       outputContainer->Add(fhPrimPi0OpeningAngle) ;
 845
 846       fhPrimPi0OpeningAngleAsym  = new TH2F
 847       ("hPrimPi0OpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{c}",100,0,1,100,0,0.5);
 848       fhPrimPi0OpeningAngleAsym->SetXTitle("|A|=| (E_{1}-E_{2}) / (E_{1}+E_{2}) |");
 849       fhPrimPi0OpeningAngleAsym->SetYTitle("#theta(rad)");
 850       outputContainer->Add(fhPrimPi0OpeningAngleAsym) ;
 851
 852       fhPrimPi0CosOpeningAngle  = new TH2F
 853       ("hPrimPi0CosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
 854       fhPrimPi0CosOpeningAngle->SetYTitle("cos (#theta) ");
 855       fhPrimPi0CosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 856       outputContainer->Add(fhPrimPi0CosOpeningAngle) ;
 857
 858       fhPrimEtaOpeningAngle  = new TH2F
 859       ("hPrimEtaOpeningAngle","Angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,0,0.5);
 860       fhPrimEtaOpeningAngle->SetYTitle("#theta(rad)");
 861       fhPrimEtaOpeningAngle->SetXTitle("E_{#eta} (GeV)");
 862       outputContainer->Add(fhPrimEtaOpeningAngle) ;
 863
 864       fhPrimEtaOpeningAngleAsym  = new TH2F
 865       ("hPrimEtaOpeningAngleAsym","Angle between all primary #gamma pair vs #it{Asymmetry}, #it{p}_{T}>5 GeV/#it{c}",100,0,1,100,0,0.5);
 866       fhPrimEtaOpeningAngleAsym->SetXTitle("|#it{A}|=| (#it{E}_{1}-#it{E}_{2}) / (#it{E}_{1}+#it{E}_{2}) |");
 867       fhPrimEtaOpeningAngleAsym->SetYTitle("#theta(rad)");
 868       outputContainer->Add(fhPrimEtaOpeningAngleAsym) ;
 869
 870
 871       fhPrimEtaCosOpeningAngle  = new TH2F
 872       ("hPrimEtaCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#eta}",nptbins,ptmin,ptmax,100,-1,1);
 873       fhPrimEtaCosOpeningAngle->SetYTitle("cos (#theta) ");
 874       fhPrimEtaCosOpeningAngle->SetXTitle("#it{E}_{ #eta} (GeV)");
 875       outputContainer->Add(fhPrimEtaCosOpeningAngle) ;
 876
 877
 878     }
 879
 880     if(fFillOriginHisto)
 881     {
 882       //Prim origin
 883       //Pi0
 884       fhPrimPi0PtOrigin     = new TH2F("hPrimPi0PtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
 885       fhPrimPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 886       fhPrimPi0PtOrigin->SetYTitle("Origin");
 887       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 888       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 889       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances ");
 890       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 891       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
 892       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
 893       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
 894       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
 895       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
 896       fhPrimPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
 897       outputContainer->Add(fhPrimPi0PtOrigin) ;
 898
 899       fhMCPi0PtOrigin     = new TH2F("hMCPi0PtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,11,0,11) ;
 900       fhMCPi0PtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 901       fhMCPi0PtOrigin->SetYTitle("Origin");
 902       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 903       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 904       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
 905       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 906       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(5 ,"#rho");
 907       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(6 ,"#omega");
 908       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(7 ,"K");
 909       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(8 ,"Other");
 910       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(9 ,"#eta");
 911       fhMCPi0PtOrigin->GetYaxis()->SetBinLabel(10 ,"#eta prime");
 912       outputContainer->Add(fhMCPi0PtOrigin) ;
 913
 914       //Eta
 915       fhPrimEtaPtOrigin     = new TH2F("hPrimEtaPtOrigin","Primary #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
 916       fhPrimEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 917       fhPrimEtaPtOrigin->SetYTitle("Origin");
 918       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 919       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 920       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
 921       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 922       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
 923       fhPrimEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime ");
 924
 925       outputContainer->Add(fhPrimEtaPtOrigin) ;
 926
 927       fhMCEtaPtOrigin     = new TH2F("hMCEtaPtOrigin","Reconstructed pair from generated #pi^{0} #it{p}_{T} vs origin",nptbins,ptmin,ptmax,7,0,7) ;
 928       fhMCEtaPtOrigin->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 929       fhMCEtaPtOrigin->SetYTitle("Origin");
 930       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(1 ,"Status 21");
 931       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(2 ,"Quark");
 932       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(3 ,"qq Resonances");
 933       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(4 ,"Resonances");
 934       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(5 ,"Other");
 935       fhMCEtaPtOrigin->GetYaxis()->SetBinLabel(6 ,"#eta prime");
 936
 937       outputContainer->Add(fhMCEtaPtOrigin) ;
 938
 939       fhMCPi0ProdVertex = new TH2F("hMCPi0ProdVertex","Selected reco pair from generated #pi^{0} #it{p}_{T} vs production vertex",
 940                                    200,0.,20.,5000,0,500) ;
 941       fhMCPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 942       fhMCPi0ProdVertex->SetYTitle("#it{R} (cm)");
 943       outputContainer->Add(fhMCPi0ProdVertex) ;
 944
 945       fhMCEtaProdVertex = new TH2F("hMCEtaProdVertex","Selected reco pair from generated #eta #it{p}_{T} vs production vertex",
 946                                    200,0.,20.,5000,0,500) ;
 947       fhMCEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 948       fhMCEtaProdVertex->SetYTitle("#it{R} (cm)");
 949       outputContainer->Add(fhMCEtaProdVertex) ;
 950
 951       fhPrimPi0ProdVertex = new TH2F("hPrimPi0ProdVertex","generated #pi^{0} #it{p}_{T} vs production vertex",
 952                                    200,0.,20.,5000,0,500) ;
 953       fhPrimPi0ProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 954       fhPrimPi0ProdVertex->SetYTitle("#it{R} (cm)");
 955       outputContainer->Add(fhPrimPi0ProdVertex) ;
 956
 957       fhPrimEtaProdVertex = new TH2F("hPrimEtaProdVertex","generated #eta #it{p}_{T} vs production vertex",
 958                                    200,0.,20.,5000,0,500) ;
 959       fhPrimEtaProdVertex->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 960       fhPrimEtaProdVertex->SetYTitle("#it{R} (cm)");
 961       outputContainer->Add(fhPrimEtaProdVertex) ;
 962
 963       for(Int_t i = 0; i<13; i++)
 964       {
 965         fhMCOrgMass[i] = new TH2F(Form("hMCOrgMass_%d",i),Form("#it{M} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 966         fhMCOrgMass[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 967         fhMCOrgMass[i]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
 968         outputContainer->Add(fhMCOrgMass[i]) ;
 969
 970         fhMCOrgAsym[i]= new TH2F(Form("hMCOrgAsym_%d",i),Form("#it{Asymmetry} vs #it{p}_{T}, origin %d",i),nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax) ;
 971         fhMCOrgAsym[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 972         fhMCOrgAsym[i]->SetYTitle("A");
 973         outputContainer->Add(fhMCOrgAsym[i]) ;
 974
 975         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) ;
 976         fhMCOrgDeltaEta[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 977         fhMCOrgDeltaEta[i]->SetYTitle("#Delta #eta");
 978         outputContainer->Add(fhMCOrgDeltaEta[i]) ;
 979
 980         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) ;
 981         fhMCOrgDeltaPhi[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 982         fhMCOrgDeltaPhi[i]->SetYTitle("#Delta #phi (rad)");
 983         outputContainer->Add(fhMCOrgDeltaPhi[i]) ;
 984
 985       }
 986
 987       if(fMultiCutAnaSim)
 988       {
 989         fhMCPi0MassPtTrue  = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 990         fhMCPi0MassPtRec   = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 991         fhMCPi0PtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 992         fhMCEtaMassPtRec   = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 993         fhMCEtaMassPtTrue  = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 994         fhMCEtaPtTruePtRec = new TH2F*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 995         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
 996           for(Int_t icell=0; icell<fNCellNCuts; icell++){
 997             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 998               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
 999
1000               fhMCPi0MassPtRec[index] = new TH2F(Form("hMCPi0MassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1001                                                  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]),
1002                                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1003               fhMCPi0MassPtRec[index]->SetXTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1004               fhMCPi0MassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1005               outputContainer->Add(fhMCPi0MassPtRec[index]) ;
1006
1007               fhMCPi0MassPtTrue[index] = new TH2F(Form("hMCPi0MassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1008                                                   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]),
1009                                                   nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1010               fhMCPi0MassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1011               fhMCPi0MassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1012               outputContainer->Add(fhMCPi0MassPtTrue[index]) ;
1013
1014               fhMCPi0PtTruePtRec[index] = new TH2F(Form("hMCPi0PtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1015                                                    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]),
1016                                                    nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1017               fhMCPi0PtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1018               fhMCPi0PtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1019               outputContainer->Add(fhMCPi0PtTruePtRec[index]) ;
1020
1021               fhMCEtaMassPtRec[index] = new TH2F(Form("hMCEtaMassPtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1022                                                  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]),
1023                                                  nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1024               fhMCEtaMassPtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1025               fhMCEtaMassPtRec[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1026               outputContainer->Add(fhMCEtaMassPtRec[index]) ;
1027
1028               fhMCEtaMassPtTrue[index] = new TH2F(Form("hMCEtaMassPtTrue_pt%d_cell%d_asym%d",ipt,icell,iasym),
1029                                                   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]),
1030                                                   nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1031               fhMCEtaMassPtTrue[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1032               fhMCEtaMassPtTrue[index]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1033               outputContainer->Add(fhMCEtaMassPtTrue[index]) ;
1034
1035               fhMCEtaPtTruePtRec[index] = new TH2F(Form("hMCEtaPtTruePtRec_pt%d_cell%d_asym%d",ipt,icell,iasym),
1036                                                    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]),
1037                                                    nptbins,ptmin,ptmax,nptbins,ptmin,ptmax) ;
1038               fhMCEtaPtTruePtRec[index]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1039               fhMCEtaPtTruePtRec[index]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1040               outputContainer->Add(fhMCEtaPtTruePtRec[index]) ;
1041             }
1042           }
1043         }
1044       }//multi cut ana
1045       else
1046       {
1047         fhMCPi0MassPtTrue  = new TH2F*[1];
1048         fhMCPi0PtTruePtRec = new TH2F*[1];
1049         fhMCEtaMassPtTrue  = new TH2F*[1];
1050         fhMCEtaPtTruePtRec = new TH2F*[1];
1051
1052         fhMCPi0MassPtTrue[0] = new TH2F("hMCPi0MassPtTrue","Reconstructed Mass vs generated p_T of true #pi^{0} cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1053         fhMCPi0MassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1054         fhMCPi0MassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1055         outputContainer->Add(fhMCPi0MassPtTrue[0]) ;
1056
1057         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) ;
1058         fhMCPi0PtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1059         fhMCPi0PtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1060         outputContainer->Add(fhMCPi0PtTruePtRec[0]) ;
1061
1062         fhMCEtaMassPtTrue[0] = new TH2F("hMCEtaMassPtTrue","Reconstructed Mass vs generated p_T of true #eta cluster pairs",nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1063         fhMCEtaMassPtTrue[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1064         fhMCEtaMassPtTrue[0]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1065         outputContainer->Add(fhMCEtaMassPtTrue[0]) ;
1066
1067         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) ;
1068         fhMCEtaPtTruePtRec[0]->SetXTitle("#it{p}_{T, generated} (GeV/#it{c})");
1069         fhMCEtaPtTruePtRec[0]->SetYTitle("#it{p}_{T, reconstructed} (GeV/#it{c})");
1070         outputContainer->Add(fhMCEtaPtTruePtRec[0]) ;
1071       }
1072     }
1073   }
1074
1075   if(fFillSMCombinations)
1076   {
1077     TString pairnamePHOS[] = {"(0-1)","(0-2)","(1-2)","(0-3)","(0-4)","(1-3)","(1-4)","(2-3)","(2-4)","(3-4)"};
1078     for(Int_t imod=0; imod<fNModules; imod++)
1079     {
1080       //Module dependent invariant mass
1081       snprintf(key, buffersize,"hReMod_%d",imod) ;
1082       snprintf(title, buffersize,"Real #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1083       fhReMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1084       fhReMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1085       fhReMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1086       outputContainer->Add(fhReMod[imod]) ;
1087       if(fCalorimeter=="PHOS")
1088       {
1089         snprintf(key, buffersize,"hReDiffPHOSMod_%d",imod) ;
1090         snprintf(title, buffersize,"Real pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1091         fhReDiffPHOSMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1092         fhReDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1093         fhReDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1094         outputContainer->Add(fhReDiffPHOSMod[imod]) ;
1095       }
1096       else{//EMCAL
1097         if(imod<fNModules/2)
1098         {
1099           snprintf(key, buffersize,"hReSameSectorEMCAL_%d",imod) ;
1100           snprintf(title, buffersize,"Real pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1101           fhReSameSectorEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1102           fhReSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1103           fhReSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1104           outputContainer->Add(fhReSameSectorEMCALMod[imod]) ;
1105         }
1106         if(imod<fNModules-2)
1107         {
1108           snprintf(key, buffersize,"hReSameSideEMCAL_%d",imod) ;
1109           snprintf(title, buffersize,"Real pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1110           fhReSameSideEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1111           fhReSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1112           fhReSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1113           outputContainer->Add(fhReSameSideEMCALMod[imod]) ;
1114         }
1115       }//EMCAL
1116
1117       if(DoOwnMix())
1118       {
1119         snprintf(key, buffersize,"hMiMod_%d",imod) ;
1120         snprintf(title, buffersize,"Mixed #it{M}_{#gamma#gamma} distr. for Module %d",imod) ;
1121         fhMiMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1122         fhMiMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1123         fhMiMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1124         outputContainer->Add(fhMiMod[imod]) ;
1125
1126         if(fCalorimeter=="PHOS"){
1127           snprintf(key, buffersize,"hMiDiffPHOSMod_%d",imod) ;
1128           snprintf(title, buffersize,"Mixed pairs PHOS, clusters in different Modules: %s",(pairnamePHOS[imod]).Data()) ;
1129           fhMiDiffPHOSMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1130           fhMiDiffPHOSMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1131           fhMiDiffPHOSMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1132           outputContainer->Add(fhMiDiffPHOSMod[imod]) ;
1133         }//PHOS
1134         else{//EMCAL
1135           if(imod<fNModules/2)
1136           {
1137             snprintf(key, buffersize,"hMiSameSectorEMCALMod_%d",imod) ;
1138             snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same sector, SM(%d,%d)",imod*2,imod*2+1) ;
1139             fhMiSameSectorEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1140             fhMiSameSectorEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1141             fhMiSameSectorEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1142             outputContainer->Add(fhMiSameSectorEMCALMod[imod]) ;
1143           }
1144           if(imod<fNModules-2){
1145
1146             snprintf(key, buffersize,"hMiSameSideEMCALMod_%d",imod) ;
1147             snprintf(title, buffersize,"Mixed pairs EMCAL, clusters in same side SM(%d,%d)",imod, imod+2) ;
1148             fhMiSameSideEMCALMod[imod]  = new TH2F(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
1149             fhMiSameSideEMCALMod[imod]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
1150             fhMiSameSideEMCALMod[imod]->SetYTitle("#it{M}_{#gamma,#gamma} (GeV/#it{c}^{2})");
1151             outputContainer->Add(fhMiSameSideEMCALMod[imod]) ;
1152           }
1153         }//EMCAL
1154       }// own mix
1155     }//loop combinations
1156   } // SM combinations
1157
1158   if(fFillArmenterosThetaStar && IsDataMC())
1159   {
1160     TString ebin[] = {"8 < E < 12 GeV","12 < E < 16 GeV", "16 < E < 20 GeV", "E > 20 GeV" };
1161     Int_t narmbins = 400;
1162     Float_t armmin = 0;
1163     Float_t armmax = 0.4;
1164
1165     for(Int_t i = 0; i < 4; i++)
1166     {
1167       fhArmPrimPi0[i] =  new TH2F(Form("hArmenterosPrimPi0EBin%d",i),
1168                                   Form("Armenteros of primary #pi^{0}, %s",ebin[i].Data()),
1169                                   200, -1, 1, narmbins,armmin,armmax);
1170       fhArmPrimPi0[i]->SetYTitle("#it{p}_{T}^{Arm}");
1171       fhArmPrimPi0[i]->SetXTitle("#alpha^{Arm}");
1172       outputContainer->Add(fhArmPrimPi0[i]) ;
1173
1174       fhArmPrimEta[i] =  new TH2F(Form("hArmenterosPrimEtaEBin%d",i),
1175                                       Form("Armenteros of primary #eta, %s",ebin[i].Data()),
1176                                       200, -1, 1, narmbins,armmin,armmax);
1177       fhArmPrimEta[i]->SetYTitle("#it{p}_{T}^{Arm}");
1178       fhArmPrimEta[i]->SetXTitle("#alpha^{Arm}");
1179       outputContainer->Add(fhArmPrimEta[i]) ;
1180
1181     }
1182
1183     // Same as asymmetry ...
1184     fhCosThStarPrimPi0  = new TH2F
1185     ("hCosThStarPrimPi0","cos(#theta *) for primary #pi^{0}",nptbins,ptmin,ptmax,200,-1,1);
1186     fhCosThStarPrimPi0->SetYTitle("cos(#theta *)");
1187     fhCosThStarPrimPi0->SetXTitle("E_{ #pi^{0}} (GeV)");
1188     outputContainer->Add(fhCosThStarPrimPi0) ;
1189
1190     fhCosThStarPrimEta  = new TH2F
1191     ("hCosThStarPrimEta","cos(#theta *) for primary #eta",nptbins,ptmin,ptmax,200,-1,1);
1192     fhCosThStarPrimEta->SetYTitle("cos(#theta *)");
1193     fhCosThStarPrimEta->SetXTitle("E_{ #eta} (GeV)");
1194     outputContainer->Add(fhCosThStarPrimEta) ;
1195
1196   }
1197
1198   //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
1199   //
1200   //    printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
1201   //
1202   //  }
1203
1204   return outputContainer;
1205 }
1206
1207 //___________________________________________________
1208 void AliAnaPi0::Print(const Option_t * /*opt*/) const
1209 {
1210   //Print some relevant parameters set for the analysis
1211   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1212   AliAnaCaloTrackCorrBaseClass::Print(" ");
1213
1214   printf("Number of bins in Centrality:  %d \n",GetNCentrBin()) ;
1215   printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
1216   printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
1217   printf("Depth of event buffer: %d \n",GetNMaxEvMix()) ;
1218   printf("Pair in same Module: %d \n",fSameSM) ;
1219   printf("Cuts: \n") ;
1220   // printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ; //It crashes here, why?
1221   printf("Number of modules:             %d \n",fNModules) ;
1222   printf("Select pairs with their angle: %d, edep %d, min angle %2.3f, max angle %2.3f \n",fUseAngleCut, fUseAngleEDepCut, fAngleCut, fAngleMaxCut) ;
1223   printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
1224   printf("\tasymmetry < ");
1225   for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
1226   printf("\n");
1227
1228   printf("PID selection bits: n = %d, \n",fNPIDBits) ;
1229   printf("\tPID bit = ");
1230   for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
1231   printf("\n");
1232
1233   if(fMultiCutAna){
1234     printf("pT cuts: n = %d, \n",fNPtCuts) ;
1235     printf("\tpT > ");
1236     for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
1237     printf("GeV/c\n");
1238
1239     printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
1240     printf("\tnCell > ");
1241     for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
1242     printf("\n");
1243
1244   }
1245   printf("------------------------------------------------------\n") ;
1246 }
1247
1248 //________________________________________
1249 void AliAnaPi0::FillAcceptanceHistograms()
1250 {
1251   //Fill acceptance histograms if MC data is available
1252
1253   Float_t cen = GetEventCentrality();
1254   Float_t ep  = GetEventPlaneAngle();
1255
1256   if(GetReader()->ReadStack())
1257   {
1258     AliStack * stack = GetMCStack();
1259     if(stack)
1260     {
1261       for(Int_t i=0 ; i<stack->GetNtrack(); i++)
1262       {
1263         if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1264
1265         TParticle * prim = stack->Particle(i) ;
1266         Int_t pdg = prim->GetPdgCode();
1267         //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
1268         //                             prim->GetName(), prim->GetPdgCode());
1269
1270         if( pdg == 111 || pdg == 221)
1271         {
1272           Double_t mesonPt = prim->Pt() ;
1273           Double_t mesonE  = prim->Energy() ;
1274           if(mesonE == TMath::Abs(prim->Pz()))
1275           {
1276             if( pdg == 111 ) fhPrimPi0PtRejected->Fill(mesonPt);
1277             else             fhPrimEtaPtRejected->Fill(mesonPt);
1278             continue ; //Protection against floating point exception
1279           }
1280
1281           Double_t mesonY    = 0.5*TMath::Log((mesonE-prim->Pz())/(mesonE+prim->Pz())) ;
1282           Double_t mesonYeta = prim->Eta() ;
1283           Double_t phi     = TMath::RadToDeg()*prim->Phi() ;
1284
1285           if(pdg == 111)
1286           {
1287             if(TMath::Abs(mesonY) < 1.0)
1288             {
1289               fhPrimPi0E  ->Fill(mesonE ) ;
1290               fhPrimPi0Pt ->Fill(mesonPt) ;
1291               fhPrimPi0Phi->Fill(mesonPt, phi) ;
1292               fhPrimPi0YetaYcut    ->Fill(mesonPt,mesonYeta) ;
1293               fhPrimPi0PtCentrality->Fill(mesonPt,cen) ;
1294               fhPrimPi0PtEventPlane->Fill(mesonPt,ep ) ;
1295             }
1296             fhPrimPi0Y   ->Fill(mesonPt, mesonY) ;
1297             fhPrimPi0Yeta->Fill(mesonPt, mesonYeta) ;
1298           }
1299           else if(pdg == 221)
1300           {
1301             if(TMath::Abs(mesonY) < 1.0)
1302             {
1303               fhPrimEtaE  ->Fill(mesonE ) ;
1304               fhPrimEtaPt ->Fill(mesonPt) ;
1305               fhPrimEtaPhi->Fill(mesonPt, phi) ;
1306               fhPrimEtaYetaYcut    ->Fill(mesonPt,mesonYeta) ;
1307               fhPrimEtaPtCentrality->Fill(mesonPt,cen) ;
1308               fhPrimEtaPtEventPlane->Fill(mesonPt,ep ) ;
1309             }
1310             fhPrimEtaY   ->Fill(mesonPt, mesonY) ;
1311             fhPrimEtaYeta->Fill(mesonPt, mesonYeta) ;
1312           }
1313
1314           //Origin of meson
1315           if(fFillOriginHisto && TMath::Abs(mesonY) < 0.7)
1316           {
1317             Int_t momindex  = prim->GetFirstMother();
1318             if(momindex >= 0)
1319             {
1320               TParticle* mother = stack->Particle(momindex);
1321               Int_t mompdg    = TMath::Abs(mother->GetPdgCode());
1322               Int_t momstatus = mother->GetStatusCode();
1323               if(pdg == 111)
1324               {
1325                 if     (momstatus  == 21)fhPrimPi0PtOrigin->Fill(mesonPt,0.5);//parton
1326                 else if(mompdg     < 22 ) fhPrimPi0PtOrigin->Fill(mesonPt,1.5);//quark
1327                 else if(mompdg     > 2100  && mompdg   < 2210) fhPrimPi0PtOrigin->Fill(mesonPt,2.5);// resonances
1328                 else if(mompdg    == 221) fhPrimPi0PtOrigin->Fill(mesonPt,8.5);//eta
1329                 else if(mompdg    == 331) fhPrimPi0PtOrigin->Fill(mesonPt,9.5);//eta prime
1330                 else if(mompdg    == 213) fhPrimPi0PtOrigin->Fill(mesonPt,4.5);//rho
1331                 else if(mompdg    == 223) fhPrimPi0PtOrigin->Fill(mesonPt,5.5);//omega
1332                 else if(mompdg    >= 310   && mompdg    <= 323) fhPrimPi0PtOrigin->Fill(mesonPt,6.5);//k0S, k+-,k*
1333                 else if(mompdg    == 130) fhPrimPi0PtOrigin->Fill(mesonPt,6.5);//k0L
1334                 else if(momstatus == 11 || momstatus  == 12 ) fhPrimPi0PtOrigin->Fill(mesonPt,3.5);//resonances
1335                 else                      fhPrimPi0PtOrigin->Fill(mesonPt,7.5);//other?
1336
1337                 //printf("Prim Pi0: index %d, pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1338                 //                   momindex, mesonPt,mother->R(),mompdg,momstatus,mother->GetUniqueID());
1339
1340                 fhPrimPi0ProdVertex->Fill(mesonPt,mother->R());
1341
1342               }//pi0
1343               else
1344               {
1345                 if     (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(mesonPt,0.5);//parton
1346                 else if(mompdg    < 22  ) fhPrimEtaPtOrigin->Fill(mesonPt,1.5);//quark
1347                 else if(mompdg    > 2100  && mompdg   < 2210) fhPrimEtaPtOrigin->Fill(mesonPt,2.5);//qq resonances
1348                 else if(mompdg    == 331) fhPrimEtaPtOrigin->Fill(mesonPt,5.5);//eta prime
1349                 else if(momstatus == 11 || momstatus  == 12 ) fhPrimEtaPtOrigin->Fill(mesonPt,3.5);//resonances
1350                 else fhPrimEtaPtOrigin->Fill(mesonPt,4.5);//stable, conversions?
1351                 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1352
1353                 fhPrimEtaProdVertex->Fill(mesonPt,mother->R());
1354
1355               }
1356             } // pi0 has mother
1357           }
1358
1359           //Check if both photons hit Calorimeter
1360           if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
1361           Int_t iphot1=prim->GetFirstDaughter() ;
1362           Int_t iphot2=prim->GetLastDaughter() ;
1363           if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack())
1364           {
1365             TParticle * phot1 = stack->Particle(iphot1) ;
1366             TParticle * phot2 = stack->Particle(iphot2) ;
1367             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
1368             {
1369               //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
1370               //        phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
1371
1372               TLorentzVector lv1, lv2,lvmeson;
1373               phot1->Momentum(lv1);
1374               phot2->Momentum(lv2);
1375               prim ->Momentum(lvmeson);
1376
1377               if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
1378
1379               Bool_t inacceptance = kFALSE;
1380               if(fCalorimeter == "PHOS")
1381               {
1382                 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
1383                 {
1384                   Int_t mod ;
1385                   Double_t x,z ;
1386                   if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
1387                     inacceptance = kTRUE;
1388                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1389                 }
1390                 else
1391                 {
1392                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1393                     inacceptance = kTRUE ;
1394                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1395                 }
1396               }
1397               else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
1398               {
1399                 if(GetEMCALGeometry())
1400                 {
1401                   Int_t absID1=0;
1402                   Int_t absID2=0;
1403
1404                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1405                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1406
1407                   if( absID1 >= 0 && absID2 >= 0)
1408                     inacceptance = kTRUE;
1409
1410                   if(inacceptance && fCheckAccInSector)
1411                   {
1412                     Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
1413                     Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
1414
1415                     Int_t j=0;
1416                     Bool_t sameSector = kFALSE;
1417                     for(Int_t isector = 0; isector < fNModules/2; isector++)
1418                     {
1419                       j=2*isector;
1420                       if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
1421                     }
1422
1423                     if(sm1!=sm2 && !sameSector)  inacceptance = kFALSE;
1424
1425                     //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
1426                   }
1427
1428                   //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
1429                   //                    inacceptance = kTRUE;
1430                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1431                 }
1432                 else
1433                 {
1434                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1435                     inacceptance = kTRUE ;
1436                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1437                 }
1438               }
1439
1440               if(inacceptance)
1441               {
1442                 Float_t  asym  = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
1443                 Double_t angle = lv1.Angle(lv2.Vect());
1444
1445                 if(pdg==111)
1446                 {
1447                   fhPrimPi0AccE   ->Fill(mesonE) ;
1448                   fhPrimPi0AccPt  ->Fill(mesonPt) ;
1449                   fhPrimPi0AccPhi ->Fill(mesonPt, phi) ;
1450                   fhPrimPi0AccY   ->Fill(mesonPt, mesonY) ;
1451                   fhPrimPi0AccYeta->Fill(mesonPt,mesonYeta) ;
1452                   fhPrimPi0AccPtCentrality->Fill(mesonPt,cen) ;
1453                   fhPrimPi0AccPtEventPlane->Fill(mesonPt,ep ) ;
1454
1455                   if(fFillAngleHisto)
1456                   {
1457                     fhPrimPi0OpeningAngle    ->Fill(mesonPt,angle);
1458                     if(mesonPt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
1459                     fhPrimPi0CosOpeningAngle ->Fill(mesonPt,TMath::Cos(angle));
1460                   }
1461                 }
1462                 else if(pdg==221)
1463                 {
1464                   fhPrimEtaAccE   ->Fill(mesonE ) ;
1465                   fhPrimEtaAccPt  ->Fill(mesonPt) ;
1466                   fhPrimEtaAccPhi ->Fill(mesonPt, phi) ;
1467                   fhPrimEtaAccY   ->Fill(mesonPt, mesonY) ;
1468                   fhPrimEtaAccYeta->Fill(mesonPt, mesonYeta) ;
1469                   fhPrimEtaAccPtCentrality->Fill(mesonPt,cen) ;
1470                   fhPrimEtaAccPtEventPlane->Fill(mesonPt,ep ) ;
1471
1472                   if(fFillAngleHisto)
1473                   {
1474                     fhPrimEtaOpeningAngle    ->Fill(mesonPt,angle);
1475                     if(mesonPt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
1476                     fhPrimEtaCosOpeningAngle ->Fill(mesonPt,TMath::Cos(angle));
1477                   }
1478                 }
1479               }//Accepted
1480             }// 2 photons
1481           }//Check daughters exist
1482         }// Primary pi0 or eta
1483       }//loop on primaries
1484     }//stack exists and data is MC
1485   }//read stack
1486   else if(GetReader()->ReadAODMCParticles())
1487   {
1488     TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1489     if(mcparticles)
1490     {
1491       Int_t nprim = mcparticles->GetEntriesFast();
1492
1493       for(Int_t i=0; i < nprim; i++)
1494       {
1495         if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
1496
1497         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
1498
1499         // Only generator particles, when they come from PYTHIA, PHOJET, HERWIG ...
1500         //if( prim->GetStatus() == 0 && (GetMCAnalysisUtils()->GetMCGenerator()).Length()!=0) break;
1501
1502         Int_t pdg = prim->GetPdgCode();
1503         if( pdg == 111 || pdg == 221)
1504         {
1505           Double_t mesonPt = prim->Pt() ;
1506           Double_t mesonE  = prim->E() ;
1507           //printf("pi0, pt %2.2f, eta %f, phi %f\n",mesonPt, prim->Eta(), prim->Phi());
1508
1509           if(mesonE == TMath::Abs(prim->Pz()))
1510           {
1511             if( pdg == 111 ) fhPrimPi0PtRejected->Fill(mesonPt);
1512             else             fhPrimEtaPtRejected->Fill(mesonPt);
1513             continue ; //Protection against floating point exception
1514           }
1515
1516           Double_t mesonY    = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
1517           Double_t mesonYeta = prim->Eta() ;
1518           Double_t phi     = TMath::RadToDeg()*prim->Phi() ;
1519
1520           if(pdg == 111)
1521           {
1522             if(TMath::Abs(mesonY) < 1)
1523             {
1524               fhPrimPi0E  ->Fill(mesonE ) ;
1525               fhPrimPi0Pt ->Fill(mesonPt) ;
1526               fhPrimPi0Phi->Fill(mesonPt, phi) ;
1527               fhPrimPi0YetaYcut    ->Fill(mesonPt,mesonYeta) ;
1528               fhPrimPi0PtCentrality->Fill(mesonPt,cen) ;
1529               fhPrimPi0PtEventPlane->Fill(mesonPt,ep ) ;
1530             }
1531             fhPrimPi0Y   ->Fill(mesonPt, mesonY   ) ;
1532             fhPrimPi0Yeta->Fill(mesonPt, mesonYeta) ;
1533           }
1534           else if(pdg == 221)
1535           {
1536             if(TMath::Abs(mesonY) < 1)
1537             {
1538               fhPrimEtaE  ->Fill(mesonE ) ;
1539               fhPrimEtaPt ->Fill(mesonPt) ;
1540               fhPrimEtaPhi->Fill(mesonPt, phi) ;
1541               fhPrimEtaYetaYcut    ->Fill(mesonPt,mesonYeta) ;
1542               fhPrimEtaPtCentrality->Fill(mesonPt,cen) ;
1543               fhPrimEtaPtEventPlane->Fill(mesonPt,ep ) ;
1544             }
1545             fhPrimEtaY   ->Fill(mesonPt, mesonY   ) ;
1546             fhPrimEtaYeta->Fill(mesonPt, mesonYeta) ;
1547           }
1548
1549           //Origin of meson
1550           Int_t momindex  = prim->GetMother();
1551           if(momindex >= 0 && TMath::Abs(mesonY) < 0.7)
1552           {
1553             AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1554             Int_t mompdg    = TMath::Abs(mother->GetPdgCode());
1555             Int_t momstatus = mother->GetStatus();
1556             if(fFillOriginHisto)
1557             {
1558               if(pdg == 111)
1559               {
1560                 if     (momstatus  == 21) fhPrimPi0PtOrigin->Fill(mesonPt,0.5);//parton
1561                 else if(mompdg     < 22 ) fhPrimPi0PtOrigin->Fill(mesonPt,1.5);//quark
1562                 else if(mompdg     > 2100  && mompdg   < 2210) fhPrimPi0PtOrigin->Fill(mesonPt,2.5);// resonances
1563                 else if(mompdg    == 221) fhPrimPi0PtOrigin->Fill(mesonPt,8.5);//eta
1564                 else if(mompdg    == 331) fhPrimPi0PtOrigin->Fill(mesonPt,9.5);//eta prime
1565                 else if(mompdg    == 213) fhPrimPi0PtOrigin->Fill(mesonPt,4.5);//rho
1566                 else if(mompdg    == 223) fhPrimPi0PtOrigin->Fill(mesonPt,5.5);//omega
1567                 else if(mompdg    >= 310   && mompdg    <= 323) fhPrimPi0PtOrigin->Fill(mesonPt,6.5);//k0S, k+-,k*
1568                 else if(mompdg    == 130) fhPrimPi0PtOrigin->Fill(mesonPt,6.5);//k0L
1569                 else if(momstatus == 11 || momstatus  == 12 ) fhPrimPi0PtOrigin->Fill(mesonPt,3.5);//resonances
1570                 else                      fhPrimPi0PtOrigin->Fill(mesonPt,7.5);//other?
1571
1572                 Float_t prodR =TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1573                //printf("Prim Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1574                //        mesonPt,prodR,mompdg,momstatus,mother->GetUniqueID());
1575
1576                 fhPrimPi0ProdVertex->Fill(mesonPt,prodR);
1577
1578               }//pi0
1579               else
1580               {
1581                 if     (momstatus == 21 ) fhPrimEtaPtOrigin->Fill(mesonPt,0.5);//parton
1582                 else if(mompdg    < 22  ) fhPrimEtaPtOrigin->Fill(mesonPt,1.5);//quark
1583                 else if(mompdg    > 2100  && mompdg   < 2210) fhPrimEtaPtOrigin->Fill(mesonPt,2.5);//qq resonances
1584                 else if(mompdg    == 331) fhPrimEtaPtOrigin->Fill(mesonPt,5.5);//eta prime
1585                 else if(momstatus == 11 || momstatus  == 12 ) fhPrimEtaPtOrigin->Fill(mesonPt,3.5);//resonances
1586                 else fhPrimEtaPtOrigin->Fill(mesonPt,4.5);//stable, conversions?
1587                 //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1588
1589                 Float_t prodR =TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1590                 //printf("Prim Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1591                 //        mesonPt,prodR,mompdg,momstatus,mother->GetUniqueID());
1592
1593                 fhPrimEtaProdVertex->Fill(mesonPt,prodR);
1594
1595               }
1596             }
1597           }//pi0 has mother
1598
1599           //Check if both photons hit Calorimeter
1600           if(prim->GetNDaughters()!=2) continue; //Only interested in 2 gamma decay
1601           Int_t iphot1=prim->GetDaughter(0) ;
1602           Int_t iphot2=prim->GetDaughter(1) ;
1603           if(iphot1>-1 && iphot1<nprim && iphot2>-1 && iphot2<nprim)
1604           {
1605             AliAODMCParticle * phot1 = (AliAODMCParticle *) mcparticles->At(iphot1);
1606             AliAODMCParticle * phot2 = (AliAODMCParticle *) mcparticles->At(iphot2);
1607             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22)
1608             {
1609               TLorentzVector lv1, lv2,lvmeson;
1610               lv1.SetPxPyPzE(phot1->Px(),phot1->Py(),phot1->Pz(),phot1->E());
1611               lv2.SetPxPyPzE(phot2->Px(),phot2->Py(),phot2->Pz(),phot2->E());
1612               lvmeson.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
1613
1614                if(fFillArmenterosThetaStar) FillArmenterosThetaStar(pdg,lvmeson,lv1,lv2);
1615
1616               Bool_t inacceptance = kFALSE;
1617               if(fCalorimeter == "PHOS")
1618               {
1619                 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet())
1620                 {
1621                   Int_t mod ;
1622                   Double_t x,z ;
1623                   Double_t vtx []={phot1->Xv(),phot1->Yv(),phot1->Zv()};
1624                   Double_t vtx2[]={phot2->Xv(),phot2->Yv(),phot2->Zv()};
1625                   if(GetPHOSGeometry()->ImpactOnEmc(vtx, phot1->Theta(),phot1->Phi(),mod,z,x) &&
1626                      GetPHOSGeometry()->ImpactOnEmc(vtx2,phot2->Theta(),phot2->Phi(),mod,z,x))
1627                     inacceptance = kTRUE;
1628                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1629                 }
1630                 else
1631                 {
1632                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1633                     inacceptance = kTRUE ;
1634                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1635                 }
1636               }
1637               else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet())
1638               {
1639                 if(GetEMCALGeometry())
1640                 {
1641                   Int_t absID1=0;
1642                   Int_t absID2=0;
1643
1644                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot1->Eta(),phot1->Phi(),absID1);
1645                   GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(phot2->Eta(),phot2->Phi(),absID2);
1646
1647                   if( absID1 >= 0 && absID2 >= 0)
1648                     inacceptance = kTRUE;
1649
1650                   if(inacceptance && fCheckAccInSector)
1651                   {
1652                     Int_t sm1 = GetEMCALGeometry()->GetSuperModuleNumber(absID1);
1653                     Int_t sm2 = GetEMCALGeometry()->GetSuperModuleNumber(absID2);
1654
1655                     Int_t j=0;
1656                     Bool_t sameSector = kFALSE;
1657                     for(Int_t isector = 0; isector < fNModules/2; isector++)
1658                     {
1659                       j=2*isector;
1660                       if((sm1==j && sm2==j+1) || (sm1==j+1 && sm2==j)) sameSector = kTRUE;
1661                     }
1662
1663                     if(sm1!=sm2 && !sameSector)  inacceptance = kFALSE;
1664
1665                     //if(sm1!=sm2)printf("sm1 %d, sm2 %d, same sector %d, in acceptance %d\n",sm1,sm2,sameSector,inacceptance);
1666                   }
1667
1668                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1669                 }
1670                 else
1671                 {
1672                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
1673                     inacceptance = kTRUE ;
1674                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
1675                 }
1676               }
1677
1678               if(inacceptance)
1679               {
1680                 Float_t  asym  = TMath::Abs((lv1.E()-lv2.E()) / (lv1.E()+lv2.E()));
1681                 Double_t angle = lv1.Angle(lv2.Vect());
1682
1683                 if(pdg==111)
1684                 {
1685                   //                printf("ACCEPTED pi0: pt %2.2f, phi %3.2f, eta %1.2f\n",mesonPt,phi,mesonY);
1686                   fhPrimPi0AccE   ->Fill(mesonE ) ;
1687                   fhPrimPi0AccPt  ->Fill(mesonPt) ;
1688                   fhPrimPi0AccPhi ->Fill(mesonPt, phi) ;
1689                   fhPrimPi0AccY   ->Fill(mesonPt, mesonY) ;
1690                   fhPrimPi0AccYeta->Fill(mesonPt, mesonYeta) ;
1691                   fhPrimPi0AccPtCentrality->Fill(mesonPt,cen) ;
1692                   fhPrimPi0AccPtEventPlane->Fill(mesonPt,ep ) ;
1693
1694                   if(fFillAngleHisto)
1695                   {
1696                     fhPrimPi0OpeningAngle    ->Fill(mesonPt,angle);
1697                     if(mesonPt > 5)fhPrimPi0OpeningAngleAsym->Fill(asym,angle);
1698                     fhPrimPi0CosOpeningAngle ->Fill(mesonPt,TMath::Cos(angle));
1699                   }
1700                 }
1701                 else if(pdg==221)
1702                 {
1703                   fhPrimEtaAccE   ->Fill(mesonE ) ;
1704                   fhPrimEtaAccPt  ->Fill(mesonPt) ;
1705                   fhPrimEtaAccPhi ->Fill(mesonPt, phi) ;
1706                   fhPrimEtaAccY   ->Fill(mesonPt, mesonY) ;
1707                   fhPrimEtaAccYeta->Fill(mesonPt, mesonYeta) ;
1708                   fhPrimEtaAccPtCentrality->Fill(mesonPt,cen) ;
1709                   fhPrimEtaAccPtEventPlane->Fill(mesonPt,ep ) ;
1710
1711                   if(fFillAngleHisto)
1712                   {
1713                     fhPrimEtaOpeningAngle    ->Fill(mesonPt,angle);
1714                     if(mesonPt > 5)fhPrimEtaOpeningAngleAsym->Fill(asym,angle);
1715                     fhPrimEtaCosOpeningAngle ->Fill(mesonPt,TMath::Cos(angle));
1716                   }
1717                 }
1718               }//Accepted
1719             }// 2 photons
1720           }//Check daughters exist
1721         }// Primary pi0 or eta
1722       }//loop on primaries
1723     }//stack exists and data is MC
1724
1725
1726   }     // read AOD MC
1727 }
1728
1729 //__________________________________________________________________________________
1730 void AliAnaPi0::FillArmenterosThetaStar(Int_t pdg,             TLorentzVector meson,
1731                                         TLorentzVector daugh1, TLorentzVector daugh2)
1732 {
1733   // Fill armenteros plots
1734
1735   // Get pTArm and AlphaArm
1736   Float_t momentumSquaredMother = meson.P()*meson.P();
1737   Float_t momentumDaughter1AlongMother = 0.;
1738   Float_t momentumDaughter2AlongMother = 0.;
1739
1740   if (momentumSquaredMother > 0.)
1741   {
1742     momentumDaughter1AlongMother = (daugh1.Px()*meson.Px() + daugh1.Py()*meson.Py()+ daugh1.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
1743     momentumDaughter2AlongMother = (daugh2.Px()*meson.Px() + daugh2.Py()*meson.Py()+ daugh2.Pz()*meson.Pz()) / sqrt(momentumSquaredMother);
1744   }
1745
1746   Float_t momentumSquaredDaughter1 = daugh1.P()*daugh1.P();
1747   Float_t ptArmSquared = momentumSquaredDaughter1 - momentumDaughter1AlongMother*momentumDaughter1AlongMother;
1748
1749   Float_t pTArm = 0.;
1750   if (ptArmSquared > 0.)
1751     pTArm = sqrt(ptArmSquared);
1752
1753   Float_t alphaArm = 0.;
1754   if(momentumDaughter1AlongMother +momentumDaughter2AlongMother > 0)
1755     alphaArm = (momentumDaughter1AlongMother -momentumDaughter2AlongMother) / (momentumDaughter1AlongMother + momentumDaughter2AlongMother);
1756
1757   TLorentzVector daugh1Boost = daugh1;
1758   daugh1Boost.Boost(-meson.BoostVector());
1759   Float_t  cosThStar=TMath::Cos(daugh1Boost.Vect().Angle(meson.Vect()));
1760
1761   Float_t en   = meson.Energy();
1762   Int_t   ebin = -1;
1763   if(en > 8  && en <= 12) ebin = 0;
1764   if(en > 12 && en <= 16) ebin = 1;
1765   if(en > 16 && en <= 20) ebin = 2;
1766   if(en > 20)             ebin = 3;
1767   if(ebin < 0 || ebin > 3) return ;
1768
1769
1770   if(pdg==111)
1771   {
1772     fhCosThStarPrimPi0->Fill(en,cosThStar);
1773     fhArmPrimPi0[ebin]->Fill(alphaArm,pTArm);
1774   }
1775   else
1776   {
1777     fhCosThStarPrimEta->Fill(en,cosThStar);
1778     fhArmPrimEta[ebin]->Fill(alphaArm,pTArm);
1779   }
1780
1781   if(GetDebug() > 2 )
1782   {
1783     Float_t asym = 0;
1784     if(daugh1.E()+daugh2.E() > 0) asym = TMath::Abs(daugh1.E()-daugh2.E())/(daugh1.E()+daugh2.E());
1785
1786     printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f, asymmetry %1.3f\n",
1787          en,alphaArm,pTArm,cosThStar,asym);
1788   }
1789 }
1790
1791 //_______________________________________________________________________________________
1792 void AliAnaPi0::FillMCVersusRecDataHistograms(Int_t index1,  Int_t index2,
1793                                               Float_t pt1,   Float_t pt2,
1794                                               Int_t ncell1,  Int_t ncell2,
1795                                               Double_t mass, Double_t pt,  Double_t asym,
1796                                               Double_t deta, Double_t dphi)
1797 {
1798   //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1799   //Adjusted for Pythia, need to see what to do for other generators.
1800   //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1801   // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1802
1803   if(!fFillOriginHisto) return;
1804
1805   Int_t ancPDG    = 0;
1806   Int_t ancStatus = 0;
1807   TLorentzVector ancMomentum;
1808   TVector3 prodVertex;
1809   Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1810                                                               GetReader(), ancPDG, ancStatus,ancMomentum, prodVertex);
1811
1812   Int_t momindex  = -1;
1813   Int_t mompdg    = -1;
1814   Int_t momstatus = -1;
1815   if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1816                             ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1817
1818   Float_t prodR = -1;
1819
1820   if(ancLabel > -1)
1821   {
1822     if(ancPDG==22){//gamma
1823       fhMCOrgMass[0]->Fill(pt,mass);
1824       fhMCOrgAsym[0]->Fill(pt,asym);
1825       fhMCOrgDeltaEta[0]->Fill(pt,deta);
1826       fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1827     }
1828     else if(TMath::Abs(ancPDG)==11){//e
1829       fhMCOrgMass[1]->Fill(pt,mass);
1830       fhMCOrgAsym[1]->Fill(pt,asym);
1831       fhMCOrgDeltaEta[1]->Fill(pt,deta);
1832       fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1833     }
1834     else if(ancPDG==111){//Pi0
1835       fhMCOrgMass[2]->Fill(pt,mass);
1836       fhMCOrgAsym[2]->Fill(pt,asym);
1837       fhMCOrgDeltaEta[2]->Fill(pt,deta);
1838       fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1839       if(fMultiCutAnaSim)
1840       {
1841         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1842           for(Int_t icell=0; icell<fNCellNCuts; icell++){
1843             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1844               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1845               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1846                  asym   <  fAsymCuts[iasym]                                   &&
1847                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1848                 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1849                 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1850                 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1851               }//pass the different cuts
1852             }// pid bit cut loop
1853           }// icell loop
1854         }// pt cut loop
1855       }//Multi cut ana sim
1856       else
1857       {
1858         fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1859
1860         if(mass < 0.17 && mass > 0.1)
1861         {
1862           fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1863
1864           if(fFillOriginHisto)
1865           {
1866             //Int_t uniqueId = -1;
1867             if(GetReader()->ReadStack())
1868             {
1869               TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1870               momindex  = ancestor->GetFirstMother();
1871               if(momindex < 0) return;
1872               TParticle* mother = GetMCStack()->Particle(momindex);
1873               mompdg    = TMath::Abs(mother->GetPdgCode());
1874               momstatus = mother->GetStatusCode();
1875               prodR = mother->R();
1876               //uniqueId = mother->GetUniqueID();
1877             }
1878             else
1879             {
1880               TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1881               AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1882               momindex  = ancestor->GetMother();
1883               if(momindex < 0) return;
1884               AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1885               mompdg    = TMath::Abs(mother->GetPdgCode());
1886               momstatus = mother->GetStatus();
1887               prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1888               //uniqueId = mother->GetUniqueID();
1889             }
1890
1891 //            printf("Reco Pi0: pt %2.2f Prod vertex %3.3f, origin pdg %d, origin status %d, origin UI %d\n",
1892 //                   pt,prodR,mompdg,momstatus,uniqueId);
1893
1894             fhMCPi0ProdVertex->Fill(pt,prodR);
1895
1896             if     (momstatus  == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1897             else if(mompdg     < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1898             else if(mompdg     > 2100  && mompdg   < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1899             else if(mompdg    == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1900             else if(mompdg    == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1901             else if(mompdg    == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1902             else if(mompdg    == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1903             else if(mompdg    >= 310   && mompdg    <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1904             else if(mompdg    == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1905             else if(momstatus == 11 || momstatus  == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1906             else                      fhMCPi0PtOrigin->Fill(pt,7.5);//other?
1907
1908           }
1909         }//pi0 mass region
1910       }
1911     }
1912     else if(ancPDG==221){//Eta
1913       fhMCOrgMass[3]->Fill(pt,mass);
1914       fhMCOrgAsym[3]->Fill(pt,asym);
1915       fhMCOrgDeltaEta[3]->Fill(pt,deta);
1916       fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1917       if(fMultiCutAnaSim){
1918         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1919           for(Int_t icell=0; icell<fNCellNCuts; icell++){
1920             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1921               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1922               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1923                  asym   <  fAsymCuts[iasym]                                   &&
1924                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1925                 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1926                 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1927                 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1928               }//pass the different cuts
1929             }// pid bit cut loop
1930           }// icell loop
1931         }// pt cut loop
1932       } //Multi cut ana sim
1933       else
1934       {
1935         fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1936         if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1937
1938         if(fFillOriginHisto)
1939         {
1940           if(GetReader()->ReadStack())
1941           {
1942             TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1943             momindex  = ancestor->GetFirstMother();
1944             if(momindex < 0) return;
1945             TParticle* mother = GetMCStack()->Particle(momindex);
1946             mompdg    = TMath::Abs(mother->GetPdgCode());
1947             momstatus = mother->GetStatusCode();
1948             prodR = mother->R();
1949           }
1950           else
1951           {
1952             TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1953             AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1954             momindex  = ancestor->GetMother();
1955             if(momindex < 0) return;
1956             AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1957             mompdg    = TMath::Abs(mother->GetPdgCode());
1958             momstatus = mother->GetStatus();
1959             prodR = TMath::Sqrt(mother->Xv()*mother->Xv()+mother->Yv()*mother->Yv());
1960           }
1961
1962           fhMCEtaProdVertex->Fill(pt,prodR);
1963
1964           if     (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1965           else if(mompdg    < 22  ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1966           else if(mompdg    > 2100  && mompdg   < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1967           else if(mompdg    == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1968           else if(momstatus == 11 || momstatus  == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1969           else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1970           //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1971         }
1972       }// eta mass region
1973     }
1974     else if(ancPDG==-2212){//AProton
1975       fhMCOrgMass[4]->Fill(pt,mass);
1976       fhMCOrgAsym[4]->Fill(pt,asym);
1977       fhMCOrgDeltaEta[4]->Fill(pt,deta);
1978       fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1979     }
1980     else if(ancPDG==-2112){//ANeutron
1981       fhMCOrgMass[5]->Fill(pt,mass);
1982       fhMCOrgAsym[5]->Fill(pt,asym);
1983       fhMCOrgDeltaEta[5]->Fill(pt,deta);
1984       fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1985     }
1986     else if(TMath::Abs(ancPDG)==13){//muons
1987       fhMCOrgMass[6]->Fill(pt,mass);
1988       fhMCOrgAsym[6]->Fill(pt,asym);
1989       fhMCOrgDeltaEta[6]->Fill(pt,deta);
1990       fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1991     }
1992     else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1993       if(ancStatus==1){//Stable particles, converted? not decayed resonances
1994         fhMCOrgMass[6]->Fill(pt,mass);
1995         fhMCOrgAsym[6]->Fill(pt,asym);
1996         fhMCOrgDeltaEta[6]->Fill(pt,deta);
1997         fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1998       }
1999       else{//resonances and other decays, more hadron conversions?
2000         fhMCOrgMass[7]->Fill(pt,mass);
2001         fhMCOrgAsym[7]->Fill(pt,asym);
2002         fhMCOrgDeltaEta[7]->Fill(pt,deta);
2003         fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
2004       }
2005     }
2006     else {//Partons, colliding protons, strings, intermediate corrections
2007       if(ancStatus==11 || ancStatus==12){//String fragmentation
2008         fhMCOrgMass[8]->Fill(pt,mass);
2009         fhMCOrgAsym[8]->Fill(pt,asym);
2010         fhMCOrgDeltaEta[8]->Fill(pt,deta);
2011         fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
2012       }
2013       else if (ancStatus==21){
2014         if(ancLabel < 2) {//Colliding protons
2015           fhMCOrgMass[11]->Fill(pt,mass);
2016           fhMCOrgAsym[11]->Fill(pt,asym);
2017           fhMCOrgDeltaEta[11]->Fill(pt,deta);
2018           fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
2019         }//colliding protons
2020         else if(ancLabel < 6){//partonic initial states interactions
2021           fhMCOrgMass[9]->Fill(pt,mass);
2022           fhMCOrgAsym[9]->Fill(pt,asym);
2023           fhMCOrgDeltaEta[9]->Fill(pt,deta);
2024           fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
2025         }
2026         else if(ancLabel < 8){//Final state partons radiations?
2027           fhMCOrgMass[10]->Fill(pt,mass);
2028           fhMCOrgAsym[10]->Fill(pt,asym);
2029           fhMCOrgDeltaEta[10]->Fill(pt,deta);
2030           fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
2031         }
2032         // else {
2033         //   printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2034         //          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2035         // }
2036       }//status 21
2037       //else {
2038       //  printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
2039       //         ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2040       // }
2041     }////Partons, colliding protons, strings, intermediate corrections
2042   }//ancLabel > -1
2043   else { //ancLabel <= -1
2044     //printf("Not related at all label = %d\n",ancLabel);
2045     fhMCOrgMass[12]->Fill(pt,mass);
2046     fhMCOrgAsym[12]->Fill(pt,asym);
2047     fhMCOrgDeltaEta[12]->Fill(pt,deta);
2048     fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
2049   }
2050 }
2051
2052 //__________________________________________
2053 void AliAnaPi0::MakeAnalysisFillHistograms()
2054 {
2055   //Process one event and extract photons from AOD branch
2056   // filled with AliAnaPhoton and fill histos with invariant mass
2057
2058   //In case of simulated data, fill acceptance histograms
2059   if(IsDataMC())FillAcceptanceHistograms();
2060
2061   //if (GetReader()->GetEventNumber()%10000 == 0)
2062   // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
2063
2064   if(!GetInputAODBranch())
2065   {
2066     printf("AliAnaPi0::MakeAnalysisFillHistograms() - No input aod photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
2067     abort();
2068   }
2069
2070   //Init some variables
2071   Int_t   nPhot    = GetInputAODBranch()->GetEntriesFast() ;
2072
2073   if(GetDebug() > 1)
2074     printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
2075
2076   //If less than photon 2 entries in the list, skip this event
2077   if(nPhot < 2 ) {
2078
2079     if(GetDebug() > 2)
2080       printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
2081
2082     if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
2083
2084     return ;
2085   }
2086
2087   Int_t ncentr = GetNCentrBin();
2088   if(GetNCentrBin()==0) ncentr = 1;
2089
2090   //Init variables
2091   Int_t module1         = -1;
2092   Int_t module2         = -1;
2093   Double_t vert[]       = {0.0, 0.0, 0.0} ; //vertex
2094   Int_t evtIndex1       = 0 ;
2095   Int_t currentEvtIndex = -1;
2096   Int_t curCentrBin     = GetEventCentralityBin();
2097   //Int_t curVzBin        = GetEventVzBin();
2098   //Int_t curRPBin        = GetEventRPBin();
2099   Int_t eventbin        = GetEventMixBin();
2100
2101   if(eventbin > GetNCentrBin()*GetNZvertBin()*GetNRPBin())
2102   {
2103      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());
2104     return;
2105   }
2106
2107   //Get shower shape information of clusters
2108   TObjArray *clusters = 0;
2109   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2110   else if(fCalorimeter=="PHOS" ) clusters = GetPHOSClusters() ;
2111
2112   //---------------------------------
2113   //First loop on photons/clusters
2114   //---------------------------------
2115   for(Int_t i1=0; i1<nPhot-1; i1++)
2116   {
2117     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2118     //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
2119
2120     // get the event index in the mixed buffer where the photon comes from
2121     // in case of mixing with analysis frame, not own mixing
2122     evtIndex1 = GetEventIndex(p1, vert) ;
2123     //printf("charge = %d\n", track->Charge());
2124     if ( evtIndex1 == -1 )
2125       return ;
2126     if ( evtIndex1 == -2 )
2127       continue ;
2128
2129     //printf("z vertex %f < %f\n",vert[2],GetZvertexCut());
2130     if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ;   //vertex cut
2131
2132
2133     if (evtIndex1 != currentEvtIndex)
2134     {
2135       //Fill event bin info
2136       if(DoOwnMix()) fhEventBin->Fill(eventbin) ;
2137       if(GetNCentrBin() > 1)
2138       {
2139         fhCentrality->Fill(curCentrBin);
2140         if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
2141       }
2142       currentEvtIndex = evtIndex1 ;
2143     }
2144
2145     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
2146
2147     //Get the momentum of this cluster
2148     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2149
2150     //Get (Super)Module number of this cluster
2151     module1 = GetModuleNumber(p1);
2152
2153     //------------------------------------------
2154     //Get index in VCaloCluster array
2155     AliVCluster *cluster1 = 0;
2156     Bool_t bFound1        = kFALSE;
2157     Int_t  caloLabel1     = p1->GetCaloLabel(0);
2158     Bool_t iclus1         =-1;
2159     if(clusters)
2160     {
2161       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
2162         AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
2163         if(cluster)
2164         {
2165           if     (cluster->GetID()==caloLabel1)
2166           {
2167             bFound1  = kTRUE  ;
2168             cluster1 = cluster;
2169             iclus1   = iclus;
2170           }
2171         }
2172         if(bFound1) break;
2173       }
2174     }// calorimeter clusters loop
2175
2176     //---------------------------------
2177     //Second loop on photons/clusters
2178     //---------------------------------
2179     for(Int_t i2=i1+1; i2<nPhot; i2++)
2180     {
2181       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
2182
2183       //In case of mixing frame, check we are not in the same event as the first cluster
2184       Int_t evtIndex2 = GetEventIndex(p2, vert) ;
2185       if ( evtIndex2 == -1 )
2186         return ;
2187       if ( evtIndex2 == -2 )
2188         continue ;
2189       if (GetMixedEvent() && (evtIndex1 == evtIndex2))
2190         continue ;
2191
2192       //------------------------------------------
2193       //Get index in VCaloCluster array
2194       AliVCluster *cluster2 = 0;
2195       Bool_t bFound2        = kFALSE;
2196       Int_t caloLabel2      = p2->GetCaloLabel(0);
2197       if(clusters){
2198         for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
2199           AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
2200           if(cluster){
2201             if(cluster->GetID()==caloLabel2) {
2202               bFound2  = kTRUE  ;
2203               cluster2 = cluster;
2204             }
2205           }
2206           if(bFound2) break;
2207         }// calorimeter clusters loop
2208       }
2209
2210       Float_t tof1  = -1;
2211       Float_t l01   = -1;
2212       if(cluster1 && bFound1){
2213         tof1  = cluster1->GetTOF()*1e9;
2214         l01   = cluster1->GetM02();
2215       }
2216       //      else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
2217       //                   p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
2218
2219       Float_t tof2  = -1;
2220       Float_t l02   = -1;
2221       if(cluster2 && bFound2)
2222       {
2223         tof2  = cluster2->GetTOF()*1e9;
2224         l02   = cluster2->GetM02();
2225
2226       }
2227       //      else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
2228       //                  p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
2229
2230       if(clusters)
2231       {
2232         Double_t t12diff = tof1-tof2;
2233         if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2234       }
2235       //------------------------------------------
2236
2237       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
2238
2239       //Get the momentum of this cluster
2240       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2241       //Get module number
2242       module2       = GetModuleNumber(p2);
2243
2244       //---------------------------------
2245       // Get pair kinematics
2246       //---------------------------------
2247       Double_t m    = (photon1 + photon2).M() ;
2248       Double_t pt   = (photon1 + photon2).Pt();
2249       Double_t deta = photon1.Eta() - photon2.Eta();
2250       Double_t dphi = photon1.Phi() - photon2.Phi();
2251       Double_t a    = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2252
2253       if(GetDebug() > 2)
2254         printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
2255
2256       //--------------------------------
2257       // Opening angle selection
2258       //--------------------------------
2259       //Check if opening angle is too large or too small compared to what is expected
2260       Double_t angle   = photon1.Angle(photon2.Vect());
2261       if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
2262         if(GetDebug() > 2)
2263           printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2264         continue;
2265       }
2266
2267       if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2268         if(GetDebug() > 2)
2269           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
2270         continue;
2271       }
2272
2273       //-------------------------------------------------------------------------------------------------
2274       //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2275       //-------------------------------------------------------------------------------------------------
2276       if(a < fAsymCuts[0] && fFillSMCombinations)
2277       {
2278         if(module1==module2 && module1 >=0 && module1<fNModules)
2279           fhReMod[module1]->Fill(pt,m) ;
2280
2281         if(fCalorimeter=="EMCAL")
2282         {
2283           // Same sector
2284           Int_t j=0;
2285           for(Int_t i = 0; i < fNModules/2; i++)
2286           {
2287             j=2*i;
2288             if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
2289           }
2290
2291           // Same side
2292           for(Int_t i = 0; i < fNModules-2; i++){
2293             if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
2294           }
2295         }//EMCAL
2296         else {//PHOS
2297           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
2298           if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
2299           if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
2300         }//PHOS
2301       }
2302
2303       //In case we want only pairs in same (super) module, check their origin.
2304       Bool_t ok = kTRUE;
2305       if(fSameSM && module1!=module2) ok=kFALSE;
2306       if(ok)
2307       {
2308         //Check if one of the clusters comes from a conversion
2309         if(fCheckConversion)
2310         {
2311           if     (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
2312           else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
2313         }
2314
2315         // Fill shower shape cut histograms
2316         if(fFillSSCombinations)
2317         {
2318           if     ( l01 > 0.01 && l01 < 0.4  &&
2319                    l02 > 0.01 && l02 < 0.4 )               fhReSS[0]->Fill(pt,m); // Tight
2320           else if( l01 > 0.4  && l02 > 0.4 )               fhReSS[1]->Fill(pt,m); // Loose
2321           else if( l01 > 0.01 && l01 < 0.4  && l02 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2322           else if( l02 > 0.01 && l02 < 0.4  && l01 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2323         }
2324
2325         //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2326         for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2327           if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
2328             for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
2329               if(a < fAsymCuts[iasym])
2330               {
2331                 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2332                 //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);
2333
2334                 if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2335
2336                 fhRe1     [index]->Fill(pt,m);
2337                 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
2338                 if(fFillBadDistHisto){
2339                   if(p1->DistToBad()>0 && p2->DistToBad()>0){
2340                     fhRe2     [index]->Fill(pt,m) ;
2341                     if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
2342                     if(p1->DistToBad()>1 && p2->DistToBad()>1){
2343                       fhRe3     [index]->Fill(pt,m) ;
2344                       if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
2345                     }// bad 3
2346                   }// bad2
2347                 }// Fill bad dist histos
2348               }//assymetry cut
2349             }// asymmetry cut loop
2350           }// bad 1
2351         }// pid bit loop
2352
2353         //Fill histograms with opening angle
2354         if(fFillAngleHisto)
2355         {
2356           fhRealOpeningAngle   ->Fill(pt,angle);
2357           fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2358         }
2359
2360         //Fill histograms with pair assymmetry
2361         if(fFillAsymmetryHisto)
2362         {
2363           fhRePtAsym->Fill(pt,a);
2364           if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
2365           if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
2366         }
2367
2368         //-------------------------------------------------------
2369         //Get the number of cells needed for multi cut analysis.
2370         //-------------------------------------------------------
2371         Int_t ncell1 = 0;
2372         Int_t ncell2 = 0;
2373         if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim))
2374         {
2375           AliVEvent * event = GetReader()->GetInputEvent();
2376           if(event){
2377             for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++)
2378             {
2379               AliVCluster *cluster = event->GetCaloCluster(iclus);
2380
2381               Bool_t is = kFALSE;
2382               if     (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
2383               else if(fCalorimeter == "PHOS"  && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
2384
2385               if(is){
2386                 if      (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
2387                 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
2388               } // PHOS or EMCAL cluster as requested in analysis
2389
2390               if(ncell2 > 0 &&  ncell1 > 0) break; // No need to continue the iteration
2391
2392             }
2393             //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
2394           }
2395         }
2396
2397         //---------
2398         // MC data
2399         //---------
2400         //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2401         if(IsDataMC())
2402         {
2403           if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2404              GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2405           {
2406             fhReMCFromConversion->Fill(pt,m);
2407           }
2408           else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2409                   !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2410           {
2411             fhReMCFromNotConversion->Fill(pt,m);
2412           }
2413           else
2414           {
2415             fhReMCFromMixConversion->Fill(pt,m);
2416           }
2417
2418           FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2419         }
2420
2421         //-----------------------
2422         //Multi cuts analysis
2423         //-----------------------
2424         if(fMultiCutAna)
2425         {
2426           //Histograms for different PID bits selection
2427           for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2428
2429             if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    &&
2430                p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
2431
2432             //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2433           } // pid bit cut loop
2434
2435           //Several pt,ncell and asymmetry cuts
2436           for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2437             for(Int_t icell=0; icell<fNCellNCuts; icell++){
2438               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2439                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2440                 if(p1->E() >   fPtCuts[ipt]      && p2->E() > fPtCuts[ipt]        &&
2441                    a        <   fAsymCuts[iasym]                                    &&
2442                    ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell]){
2443                   fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
2444                   //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2445                   if(fFillSMCombinations && module1==module2){
2446                     fhRePtNCellAsymCutsSM[module1][index]->Fill(pt,m) ;
2447                   }
2448                 }
2449               }// pid bit cut loop
2450             }// icell loop
2451           }// pt cut loop
2452           if(GetHistogramRanges()->GetHistoTrackMultiplicityBins()){
2453             for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
2454               if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
2455             }
2456           }
2457         }// multiple cuts analysis
2458       }// ok if same sm
2459     }// second same event particle
2460   }// first cluster
2461
2462   //-------------------------------------------------------------
2463   // Mixing
2464   //-------------------------------------------------------------
2465   if(DoOwnMix())
2466   {
2467     //Recover events in with same characteristics as the current event
2468
2469     //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2470     if(eventbin < 0) return ;
2471
2472     TList * evMixList=fEventsList[eventbin] ;
2473
2474     if(!evMixList)
2475     {
2476       printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mix event list not available, bin %d \n",eventbin);
2477       return;
2478     }
2479
2480     Int_t nMixed = evMixList->GetSize() ;
2481     for(Int_t ii=0; ii<nMixed; ii++)
2482     {
2483       TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2484       Int_t nPhot2=ev2->GetEntriesFast() ;
2485       Double_t m = -999;
2486       if(GetDebug() > 1)
2487         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, GetEventCentralityBin());
2488
2489       fhEventMixBin->Fill(eventbin) ;
2490
2491       //---------------------------------
2492       //First loop on photons/clusters
2493       //---------------------------------
2494       for(Int_t i1=0; i1<nPhot; i1++){
2495         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2496         if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2497
2498         //Get kinematics of cluster and (super) module of this cluster
2499         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2500         module1 = GetModuleNumber(p1);
2501
2502         //---------------------------------
2503         //First loop on photons/clusters
2504         //---------------------------------
2505         for(Int_t i2=0; i2<nPhot2; i2++){
2506           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2507
2508           //Get kinematics of second cluster and calculate those of the pair
2509           TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2510           m           = (photon1+photon2).M() ;
2511           Double_t pt = (photon1 + photon2).Pt();
2512           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2513
2514           //Check if opening angle is too large or too small compared to what is expected
2515           Double_t angle   = photon1.Angle(photon2.Vect());
2516           if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
2517             if(GetDebug() > 2)
2518               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2519             continue;
2520           }
2521           if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2522             if(GetDebug() > 2)
2523               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2524             continue;
2525
2526           }
2527
2528           if(GetDebug() > 2)
2529             printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
2530                    p1->Pt(), p2->Pt(), pt,m,a);
2531
2532           //In case we want only pairs in same (super) module, check their origin.
2533           module2 = GetModuleNumber(p2);
2534
2535           //-------------------------------------------------------------------------------------------------
2536           //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2537           //-------------------------------------------------------------------------------------------------
2538           if(a < fAsymCuts[0] && fFillSMCombinations){
2539             if(module1==module2 && module1 >=0 && module1<fNModules)
2540               fhMiMod[module1]->Fill(pt,m) ;
2541
2542             if(fCalorimeter=="EMCAL"){
2543
2544               // Same sector
2545               Int_t j=0;
2546               for(Int_t i = 0; i < fNModules/2; i++){
2547                 j=2*i;
2548                 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
2549               }
2550
2551               // Same side
2552               for(Int_t i = 0; i < fNModules-2; i++){
2553                 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
2554               }
2555             }//EMCAL
2556             else {//PHOS
2557               if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2558               if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2559               if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
2560             }//PHOS
2561
2562
2563           }
2564
2565           Bool_t ok = kTRUE;
2566           if(fSameSM && module1!=module2) ok=kFALSE;
2567           if(ok){
2568
2569             //Check if one of the clusters comes from a conversion
2570             if(fCheckConversion){
2571               if     (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2572               else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2573             }
2574             //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2575             for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2576               if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
2577               {
2578                 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2579                 {
2580                   if(a < fAsymCuts[iasym])
2581                   {
2582                     Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2583
2584                     if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2585
2586                     fhMi1     [index]->Fill(pt,m) ;
2587                     if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
2588                     if(fFillBadDistHisto)
2589                     {
2590                       if(p1->DistToBad()>0 && p2->DistToBad()>0)
2591                       {
2592                         fhMi2     [index]->Fill(pt,m) ;
2593                         if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
2594                         if(p1->DistToBad()>1 && p2->DistToBad()>1)
2595                         {
2596                           fhMi3     [index]->Fill(pt,m) ;
2597                           if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
2598                         }
2599                       }
2600                     }// Fill bad dist histo
2601                   }//Asymmetry cut
2602                 }// Asymmetry loop
2603               }//PID cut
2604             }//loop for histograms
2605
2606             //-----------------------
2607             //Multi cuts analysis
2608             //-----------------------
2609             if(fMultiCutAna){
2610               //Several pt,ncell and asymmetry cuts
2611
2612               for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2613                 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2614                   for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2615                     Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2616                     if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        &&
2617                        a        <   fAsymCuts[iasym]                                    //&&
2618                        //p1->GetBtag() >=  fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2619                        ){
2620                       fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
2621                       //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2622                     }
2623                   }// pid bit cut loop
2624                 }// icell loop
2625               }// pt cut loop
2626             } // Multi cut ana
2627
2628             //Fill histograms with opening angle
2629             if(fFillAngleHisto){
2630               fhMixedOpeningAngle   ->Fill(pt,angle);
2631               fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2632             }
2633
2634           }//ok
2635         }// second cluster loop
2636       }//first cluster loop
2637     }//loop on mixed events
2638
2639     //--------------------------------------------------------
2640     //Add the current event to the list of events for mixing
2641     //--------------------------------------------------------
2642     TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2643     //Add current event to buffer and Remove redundant events
2644     if(currentEvent->GetEntriesFast()>0){
2645       evMixList->AddFirst(currentEvent) ;
2646       currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2647       if(evMixList->GetSize() >= GetNMaxEvMix())
2648       {
2649         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2650         evMixList->RemoveLast() ;
2651         delete tmp ;
2652       }
2653     }
2654     else{ //empty event
2655       delete currentEvent ;
2656       currentEvent=0 ;
2657     }
2658   }// DoOwnMix
2659
2660 }
2661
2662 //________________________________________________________________________
2663 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2664 {
2665   // retieves the event index and checks the vertex
2666   //    in the mixed buffer returns -2 if vertex NOK
2667   //    for normal events   returns 0 if vertex OK and -1 if vertex NOK
2668
2669   Int_t evtIndex = -1 ;
2670   if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2671
2672     if (GetMixedEvent()){
2673
2674       evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2675       GetVertex(vert,evtIndex);
2676
2677       if(TMath::Abs(vert[2])> GetZvertexCut())
2678         evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2679     } else {// Single event
2680
2681       GetVertex(vert);
2682
2683       if(TMath::Abs(vert[2])> GetZvertexCut())
2684         evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2685       else
2686         evtIndex = 0 ;
2687     }
2688   }//No MC reader
2689   else {
2690     evtIndex = 0;
2691     vert[0] = 0. ;
2692     vert[1] = 0. ;
2693     vert[2] = 0. ;
2694   }
2695
2696   return evtIndex ;
2697 }
2698