PWGGA/CaloTrackCorrelations/AliAnaPi0.cxx

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