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