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