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