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