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(const Int_t pdg, const TLorentzVector meson,
1570                                const TLorentzVector daugh1, const 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   Float_t asym = 0;
1609   if(daugh1.E()+daugh2.E() > 0) asym = TMath::Abs(daugh1.E()-daugh2.E())/(daugh1.E()+daugh2.E());
1610
1611   if(pdg==111)
1612   {
1613     fhCosThStarPrimPi0->Fill(en,cosThStar);
1614     fhArmPrimPi0[ebin]->Fill(alphaArm,pTArm);
1615   }
1616   else
1617   {
1618     fhCosThStarPrimEta->Fill(en,cosThStar);
1619     fhArmPrimEta[ebin]->Fill(alphaArm,pTArm);
1620   }
1621
1622   if(GetDebug() > 2 ) printf("AliAnaPi0::FillArmenterosThetaStar() - E %f, alphaArm %f, pTArm %f, cos(theta*) %f\n",en,alphaArm,pTArm,cosThStar);
1623
1624 }
1625
1626
1627 //_____________________________________________________________
1628 void AliAnaPi0::FillMCVersusRecDataHistograms(const Int_t index1,  const Int_t index2,
1629                                               const Float_t pt1,   const Float_t pt2,
1630                                               const Int_t ncell1,  const Int_t ncell2,
1631                                               const Double_t mass, const Double_t pt, const Double_t asym,
1632                                               const Double_t deta, const Double_t dphi){
1633   //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
1634   //Adjusted for Pythia, need to see what to do for other generators.
1635   //Array of histograms ordered as follows: 0-Photon, 1-electron, 2-pi0, 3-eta, 4-a-proton, 5-a-neutron, 6-stable particles,
1636   // 7-other decays, 8-string, 9-final parton, 10-initial parton, intermediate, 11-colliding proton, 12-unrelated
1637
1638   if(!fFillOriginHisto) return;
1639
1640   Int_t ancPDG    = 0;
1641   Int_t ancStatus = 0;
1642   TLorentzVector ancMomentum;
1643   TVector3 prodVertex;
1644   Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(index1, index2,
1645                                                               GetReader(), ancPDG, ancStatus,ancMomentum, prodVertex);
1646
1647   Int_t momindex  = -1;
1648   Int_t mompdg    = -1;
1649   Int_t momstatus = -1;
1650   if(GetDebug() > 1) printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
1651                             ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1652
1653   if(ancLabel > -1){
1654     if(ancPDG==22){//gamma
1655       fhMCOrgMass[0]->Fill(pt,mass);
1656       fhMCOrgAsym[0]->Fill(pt,asym);
1657       fhMCOrgDeltaEta[0]->Fill(pt,deta);
1658       fhMCOrgDeltaPhi[0]->Fill(pt,dphi);
1659     }
1660     else if(TMath::Abs(ancPDG)==11){//e
1661       fhMCOrgMass[1]->Fill(pt,mass);
1662       fhMCOrgAsym[1]->Fill(pt,asym);
1663       fhMCOrgDeltaEta[1]->Fill(pt,deta);
1664       fhMCOrgDeltaPhi[1]->Fill(pt,dphi);
1665     }
1666     else if(ancPDG==111){//Pi0
1667       fhMCOrgMass[2]->Fill(pt,mass);
1668       fhMCOrgAsym[2]->Fill(pt,asym);
1669       fhMCOrgDeltaEta[2]->Fill(pt,deta);
1670       fhMCOrgDeltaPhi[2]->Fill(pt,dphi);
1671       if(fMultiCutAnaSim){
1672         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1673           for(Int_t icell=0; icell<fNCellNCuts; icell++){
1674             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1675               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1676               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1677                  asym   <  fAsymCuts[iasym]                                   &&
1678                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1679                 fhMCPi0MassPtRec [index]->Fill(pt,mass);
1680                 fhMCPi0MassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1681                 if(mass < 0.17 && mass > 0.1) fhMCPi0PtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1682               }//pass the different cuts
1683             }// pid bit cut loop
1684           }// icell loop
1685         }// pt cut loop
1686       }//Multi cut ana sim
1687       else {
1688         fhMCPi0MassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1689         if(mass < 0.17 && mass > 0.1) {
1690           fhMCPi0PtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1691           if(fFillOriginHisto)
1692           {
1693             if(GetReader()->ReadStack())
1694             {
1695               TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1696               momindex  = ancestor->GetFirstMother();
1697               if(momindex < 0) return;
1698               TParticle* mother = GetMCStack()->Particle(momindex);
1699               mompdg    = TMath::Abs(mother->GetPdgCode());
1700               momstatus = mother->GetStatusCode();
1701             }
1702             else
1703             {
1704               TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1705               AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1706               momindex  = ancestor->GetMother();
1707               if(momindex < 0) return;
1708               AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1709               mompdg    = TMath::Abs(mother->GetPdgCode());
1710               momstatus = mother->GetStatus();
1711             }
1712
1713             if     (momstatus  == 21) fhMCPi0PtOrigin->Fill(pt,0.5);//parton
1714             else if(mompdg     < 22 ) fhMCPi0PtOrigin->Fill(pt,1.5);//quark
1715             else if(mompdg     > 2100  && mompdg   < 2210) fhMCPi0PtOrigin->Fill(pt,2.5);// resonances
1716             else if(mompdg    == 221) fhMCPi0PtOrigin->Fill(pt,8.5);//eta
1717             else if(mompdg    == 331) fhMCPi0PtOrigin->Fill(pt,9.5);//eta prime
1718             else if(mompdg    == 213) fhMCPi0PtOrigin->Fill(pt,4.5);//rho
1719             else if(mompdg    == 223) fhMCPi0PtOrigin->Fill(pt,5.5);//omega
1720             else if(mompdg    >= 310   && mompdg    <= 323) fhMCPi0PtOrigin->Fill(pt,6.5);//k0S, k+-,k*
1721             else if(mompdg    == 130) fhMCPi0PtOrigin->Fill(pt,6.5);//k0L
1722             else if(momstatus == 11 || momstatus  == 12 ) fhMCPi0PtOrigin->Fill(pt,3.5);//resonances
1723             else                      fhMCPi0PtOrigin->Fill(pt,7.5);//other?
1724           }
1725         }//pi0 mass region
1726
1727       }
1728     }
1729     else if(ancPDG==221){//Eta
1730       fhMCOrgMass[3]->Fill(pt,mass);
1731       fhMCOrgAsym[3]->Fill(pt,asym);
1732       fhMCOrgDeltaEta[3]->Fill(pt,deta);
1733       fhMCOrgDeltaPhi[3]->Fill(pt,dphi);
1734       if(fMultiCutAnaSim){
1735         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
1736           for(Int_t icell=0; icell<fNCellNCuts; icell++){
1737             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
1738               Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
1739               if(pt1    >  fPtCuts[ipt]      && pt2    >  fPtCuts[ipt]        &&
1740                  asym   <  fAsymCuts[iasym]                                   &&
1741                  ncell1 >= fCellNCuts[icell] && ncell2 >= fCellNCuts[icell]){
1742                 fhMCEtaMassPtRec [index]->Fill(pt,mass);
1743                 fhMCEtaMassPtTrue[index]->Fill(ancMomentum.Pt(),mass);
1744                 if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[index]->Fill(ancMomentum.Pt(),pt);
1745               }//pass the different cuts
1746             }// pid bit cut loop
1747           }// icell loop
1748         }// pt cut loop
1749       } //Multi cut ana sim
1750       else {
1751         fhMCEtaMassPtTrue[0]->Fill(ancMomentum.Pt(),mass);
1752         if(mass < 0.65 && mass > 0.45) fhMCEtaPtTruePtRec[0]->Fill(ancMomentum.Pt(),pt);
1753
1754         if(fFillOriginHisto)
1755         {
1756           if(GetReader()->ReadStack())
1757           {
1758             TParticle* ancestor = GetMCStack()->Particle(ancLabel);
1759             momindex  = ancestor->GetFirstMother();
1760             if(momindex < 0) return;
1761             TParticle* mother = GetMCStack()->Particle(momindex);
1762             mompdg    = TMath::Abs(mother->GetPdgCode());
1763             momstatus = mother->GetStatusCode();
1764           }
1765           else
1766           {
1767             TClonesArray * mcparticles = GetReader()->GetAODMCParticles();
1768             AliAODMCParticle* ancestor = (AliAODMCParticle *) mcparticles->At(ancLabel);
1769             momindex  = ancestor->GetMother();
1770             if(momindex < 0) return;
1771             AliAODMCParticle* mother = (AliAODMCParticle *) mcparticles->At(momindex);
1772             mompdg    = TMath::Abs(mother->GetPdgCode());
1773             momstatus = mother->GetStatus();
1774           }
1775
1776           if     (momstatus == 21 ) fhMCEtaPtOrigin->Fill(pt,0.5);//parton
1777           else if(mompdg    < 22  ) fhMCEtaPtOrigin->Fill(pt,1.5);//quark
1778           else if(mompdg    > 2100  && mompdg   < 2210) fhMCEtaPtOrigin->Fill(pt,2.5);//qq resonances
1779           else if(mompdg    == 331) fhMCEtaPtOrigin->Fill(pt,5.5);//eta prime
1780           else if(momstatus == 11 || momstatus  == 12 ) fhMCEtaPtOrigin->Fill(pt,3.5);//resonances
1781           else fhMCEtaPtOrigin->Fill(pt,4.5);//stable, conversions?
1782           //printf("Other Meson pdg %d, Mother %s, pdg %d, status %d\n",pdg, TDatabasePDG::Instance()->GetParticle(mompdg)->GetName(),mompdg, momstatus );
1783         }
1784       }// eta mass region
1785     }
1786     else if(ancPDG==-2212){//AProton
1787       fhMCOrgMass[4]->Fill(pt,mass);
1788       fhMCOrgAsym[4]->Fill(pt,asym);
1789       fhMCOrgDeltaEta[4]->Fill(pt,deta);
1790       fhMCOrgDeltaPhi[4]->Fill(pt,dphi);
1791     }
1792     else if(ancPDG==-2112){//ANeutron
1793       fhMCOrgMass[5]->Fill(pt,mass);
1794       fhMCOrgAsym[5]->Fill(pt,asym);
1795       fhMCOrgDeltaEta[5]->Fill(pt,deta);
1796       fhMCOrgDeltaPhi[5]->Fill(pt,dphi);
1797     }
1798     else if(TMath::Abs(ancPDG)==13){//muons
1799       fhMCOrgMass[6]->Fill(pt,mass);
1800       fhMCOrgAsym[6]->Fill(pt,asym);
1801       fhMCOrgDeltaEta[6]->Fill(pt,deta);
1802       fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1803     }
1804     else if (TMath::Abs(ancPDG) > 100 && ancLabel > 7) {
1805       if(ancStatus==1){//Stable particles, converted? not decayed resonances
1806         fhMCOrgMass[6]->Fill(pt,mass);
1807         fhMCOrgAsym[6]->Fill(pt,asym);
1808         fhMCOrgDeltaEta[6]->Fill(pt,deta);
1809         fhMCOrgDeltaPhi[6]->Fill(pt,dphi);
1810       }
1811       else{//resonances and other decays, more hadron conversions?
1812         fhMCOrgMass[7]->Fill(pt,mass);
1813         fhMCOrgAsym[7]->Fill(pt,asym);
1814         fhMCOrgDeltaEta[7]->Fill(pt,deta);
1815         fhMCOrgDeltaPhi[7]->Fill(pt,dphi);
1816       }
1817     }
1818     else {//Partons, colliding protons, strings, intermediate corrections
1819       if(ancStatus==11 || ancStatus==12){//String fragmentation
1820         fhMCOrgMass[8]->Fill(pt,mass);
1821         fhMCOrgAsym[8]->Fill(pt,asym);
1822         fhMCOrgDeltaEta[8]->Fill(pt,deta);
1823         fhMCOrgDeltaPhi[8]->Fill(pt,dphi);
1824       }
1825       else if (ancStatus==21){
1826         if(ancLabel < 2) {//Colliding protons
1827           fhMCOrgMass[11]->Fill(pt,mass);
1828           fhMCOrgAsym[11]->Fill(pt,asym);
1829           fhMCOrgDeltaEta[11]->Fill(pt,deta);
1830           fhMCOrgDeltaPhi[11]->Fill(pt,dphi);
1831         }//colliding protons
1832         else if(ancLabel < 6){//partonic initial states interactions
1833           fhMCOrgMass[9]->Fill(pt,mass);
1834           fhMCOrgAsym[9]->Fill(pt,asym);
1835           fhMCOrgDeltaEta[9]->Fill(pt,deta);
1836           fhMCOrgDeltaPhi[9]->Fill(pt,dphi);
1837         }
1838         else if(ancLabel < 8){//Final state partons radiations?
1839           fhMCOrgMass[10]->Fill(pt,mass);
1840           fhMCOrgAsym[10]->Fill(pt,asym);
1841           fhMCOrgDeltaEta[10]->Fill(pt,deta);
1842           fhMCOrgDeltaPhi[10]->Fill(pt,dphi);
1843         }
1844         // else {
1845         //   printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check ** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1846         //          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1847         // }
1848       }//status 21
1849       //else {
1850       //  printf("AliAnaPi0::FillMCVersusRecDataHistograms() - Check *** Common ancestor label %d, pdg %d, name %s, status %d; \n",
1851       //         ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
1852       // }
1853     }////Partons, colliding protons, strings, intermediate corrections
1854   }//ancLabel > -1
1855   else { //ancLabel <= -1
1856     //printf("Not related at all label = %d\n",ancLabel);
1857     fhMCOrgMass[12]->Fill(pt,mass);
1858     fhMCOrgAsym[12]->Fill(pt,asym);
1859     fhMCOrgDeltaEta[12]->Fill(pt,deta);
1860     fhMCOrgDeltaPhi[12]->Fill(pt,dphi);
1861   }
1862 }
1863
1864 //__________________________________________
1865 void AliAnaPi0::MakeAnalysisFillHistograms()
1866 {
1867   //Process one event and extract photons from AOD branch
1868   // filled with AliAnaPhoton and fill histos with invariant mass
1869
1870   //In case of simulated data, fill acceptance histograms
1871   if(IsDataMC())FillAcceptanceHistograms();
1872
1873   //if (GetReader()->GetEventNumber()%10000 == 0)
1874   // printf("--- Event %d ---\n",GetReader()->GetEventNumber());
1875
1876   if(!GetInputAODBranch())
1877   {
1878     printf("AliAnaPi0::MakeAnalysisFillHistograms() - No input aod photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1879     abort();
1880   }
1881
1882   //Init some variables
1883   Int_t   nPhot    = GetInputAODBranch()->GetEntriesFast() ;
1884
1885   if(GetDebug() > 1)
1886     printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
1887
1888   //If less than photon 2 entries in the list, skip this event
1889   if(nPhot < 2 ) {
1890
1891     if(GetDebug() > 2)
1892       printf("AliAnaPi0::MakeAnalysisFillHistograms() - nPhotons %d, cent bin %d continue to next event\n",nPhot, GetEventCentrality());
1893
1894     if(GetNCentrBin() > 1) fhCentralityNoPair->Fill(GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt());
1895
1896     return ;
1897   }
1898
1899   Int_t ncentr = GetNCentrBin();
1900   if(GetNCentrBin()==0) ncentr = 1;
1901
1902   //Init variables
1903   Int_t module1         = -1;
1904   Int_t module2         = -1;
1905   Double_t vert[]       = {0.0, 0.0, 0.0} ; //vertex
1906   Int_t evtIndex1       = 0 ;
1907   Int_t currentEvtIndex = -1;
1908   Int_t curCentrBin     = GetEventCentralityBin();
1909   //Int_t curVzBin        = GetEventVzBin();
1910   //Int_t curRPBin        = GetEventRPBin();
1911   Int_t eventbin        = GetEventMixBin();
1912
1913   //Get shower shape information of clusters
1914   TObjArray *clusters = 0;
1915   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1916   else if(fCalorimeter=="PHOS" ) clusters = GetPHOSClusters() ;
1917
1918   //---------------------------------
1919   //First loop on photons/clusters
1920   //---------------------------------
1921   for(Int_t i1=0; i1<nPhot-1; i1++)
1922   {
1923     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
1924     //printf("AliAnaPi0::MakeAnalysisFillHistograms() : cluster1 id %d\n",p1->GetCaloLabel(0));
1925
1926     // get the event index in the mixed buffer where the photon comes from
1927     // in case of mixing with analysis frame, not own mixing
1928     evtIndex1 = GetEventIndex(p1, vert) ;
1929     //printf("charge = %d\n", track->Charge());
1930     if ( evtIndex1 == -1 )
1931       return ;
1932     if ( evtIndex1 == -2 )
1933       continue ;
1934
1935     //printf("z vertex %f < %f\n",vert[2],GetZvertexCut());
1936     if(TMath::Abs(vert[2]) > GetZvertexCut()) continue ;   //vertex cut
1937
1938
1939     if (evtIndex1 != currentEvtIndex)
1940     {
1941       //Fill event bin info
1942       fhEventBin->Fill(eventbin) ;
1943       if(GetNCentrBin() > 1)
1944       {
1945         fhCentrality->Fill(curCentrBin);
1946         if(GetNRPBin() > 1 && GetEventPlane()) fhEventPlaneResolution->Fill(curCentrBin,TMath::Cos(2.*GetEventPlane()->GetQsubRes()));
1947       }
1948       currentEvtIndex = evtIndex1 ;
1949     }
1950
1951     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
1952
1953     //Get the momentum of this cluster
1954     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
1955
1956     //Get (Super)Module number of this cluster
1957     module1 = GetModuleNumber(p1);
1958
1959     //------------------------------------------
1960     //Get index in VCaloCluster array
1961     AliVCluster *cluster1 = 0;
1962     Bool_t bFound1        = kFALSE;
1963     Int_t  caloLabel1     = p1->GetCaloLabel(0);
1964     Bool_t iclus1         =-1;
1965     if(clusters)
1966     {
1967       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
1968         AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
1969         if(cluster)
1970         {
1971           if     (cluster->GetID()==caloLabel1)
1972           {
1973             bFound1  = kTRUE  ;
1974             cluster1 = cluster;
1975             iclus1   = iclus;
1976           }
1977         }
1978         if(bFound1) break;
1979       }
1980     }// calorimeter clusters loop
1981
1982     //---------------------------------
1983     //Second loop on photons/clusters
1984     //---------------------------------
1985     for(Int_t i2=i1+1; i2<nPhot; i2++)
1986     {
1987       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
1988
1989       //In case of mixing frame, check we are not in the same event as the first cluster
1990       Int_t evtIndex2 = GetEventIndex(p2, vert) ;
1991       if ( evtIndex2 == -1 )
1992         return ;
1993       if ( evtIndex2 == -2 )
1994         continue ;
1995       if (GetMixedEvent() && (evtIndex1 == evtIndex2))
1996         continue ;
1997
1998       //------------------------------------------
1999       //Get index in VCaloCluster array
2000       AliVCluster *cluster2 = 0;
2001       Bool_t bFound2        = kFALSE;
2002       Int_t caloLabel2      = p2->GetCaloLabel(0);
2003       if(clusters){
2004         for(Int_t iclus = iclus1+1; iclus < clusters->GetEntriesFast(); iclus++){
2005           AliVCluster *cluster= dynamic_cast<AliVCluster*> (clusters->At(iclus));
2006           if(cluster){
2007             if(cluster->GetID()==caloLabel2) {
2008               bFound2  = kTRUE  ;
2009               cluster2 = cluster;
2010             }
2011           }
2012           if(bFound2) break;
2013         }// calorimeter clusters loop
2014       }
2015
2016       Float_t tof1  = -1;
2017       Float_t l01   = -1;
2018       if(cluster1 && bFound1){
2019         tof1  = cluster1->GetTOF()*1e9;
2020         l01   = cluster1->GetM02();
2021       }
2022       //      else printf("cluster1 not available: calo label %d / %d, cluster ID %d\n",
2023       //                   p1->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster1->GetID());
2024
2025       Float_t tof2  = -1;
2026       Float_t l02   = -1;
2027       if(cluster2 && bFound2)
2028       {
2029         tof2  = cluster2->GetTOF()*1e9;
2030         l02   = cluster2->GetM02();
2031
2032       }
2033       //      else printf("cluster2 not available: calo label %d / %d, cluster ID %d\n",
2034       //                  p2->GetCaloLabel(0),(GetReader()->GetInputEvent())->GetNumberOfCaloClusters()-1,cluster2->GetID());
2035
2036       if(clusters)
2037       {
2038         Double_t t12diff = tof1-tof2;
2039         if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2040       }
2041       //------------------------------------------
2042
2043       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
2044
2045       //Get the momentum of this cluster
2046       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2047       //Get module number
2048       module2       = GetModuleNumber(p2);
2049
2050       //---------------------------------
2051       // Get pair kinematics
2052       //---------------------------------
2053       Double_t m    = (photon1 + photon2).M() ;
2054       Double_t pt   = (photon1 + photon2).Pt();
2055       Double_t deta = photon1.Eta() - photon2.Eta();
2056       Double_t dphi = photon1.Phi() - photon2.Phi();
2057       Double_t a    = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2058
2059       if(GetDebug() > 2)
2060         printf(" E: photon1 %f, photon2 %f; Pair: pT %f, mass %f, a %f\n", p1->E(), p2->E(), (photon1 + photon2).E(),m,a);
2061
2062       //--------------------------------
2063       // Opening angle selection
2064       //--------------------------------
2065       //Check if opening angle is too large or too small compared to what is expected
2066       Double_t angle   = photon1.Angle(photon2.Vect());
2067       if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)) {
2068         if(GetDebug() > 2)
2069           printf("AliAnaPi0::MakeAnalysisFillHistograms() -Real pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2070         continue;
2071       }
2072
2073       if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2074         if(GetDebug() > 2)
2075           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Real pair cut %f < angle %f < cut %f\n",fAngleCut, angle, fAngleMaxCut);
2076         continue;
2077       }
2078
2079       //-------------------------------------------------------------------------------------------------
2080       //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2081       //-------------------------------------------------------------------------------------------------
2082       if(a < fAsymCuts[0] && fFillSMCombinations)
2083       {
2084         if(module1==module2 && module1 >=0 && module1<fNModules)
2085           fhReMod[module1]->Fill(pt,m) ;
2086
2087         if(fCalorimeter=="EMCAL")
2088         {
2089           // Same sector
2090           Int_t j=0;
2091           for(Int_t i = 0; i < fNModules/2; i++)
2092           {
2093             j=2*i;
2094             if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhReSameSectorEMCALMod[i]->Fill(pt,m) ;
2095           }
2096
2097           // Same side
2098           for(Int_t i = 0; i < fNModules-2; i++){
2099             if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhReSameSideEMCALMod[i]->Fill(pt,m);
2100           }
2101         }//EMCAL
2102         else {//PHOS
2103           if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffPHOSMod[0]->Fill(pt,m) ;
2104           if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffPHOSMod[1]->Fill(pt,m) ;
2105           if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffPHOSMod[2]->Fill(pt,m) ;
2106         }//PHOS
2107       }
2108
2109       //In case we want only pairs in same (super) module, check their origin.
2110       Bool_t ok = kTRUE;
2111       if(fSameSM && module1!=module2) ok=kFALSE;
2112       if(ok)
2113       {
2114         //Check if one of the clusters comes from a conversion
2115         if(fCheckConversion)
2116         {
2117           if     (p1->IsTagged() && p2->IsTagged()) fhReConv2->Fill(pt,m);
2118           else if(p1->IsTagged() || p2->IsTagged()) fhReConv ->Fill(pt,m);
2119         }
2120
2121         // Fill shower shape cut histograms
2122         if(fFillSSCombinations)
2123         {
2124           if     ( l01 > 0.01 && l01 < 0.4  &&
2125                    l02 > 0.01 && l02 < 0.4 )               fhReSS[0]->Fill(pt,m); // Tight
2126           else if( l01 > 0.4  && l02 > 0.4 )               fhReSS[1]->Fill(pt,m); // Loose
2127           else if( l01 > 0.01 && l01 < 0.4  && l02 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2128           else if( l02 > 0.01 && l02 < 0.4  && l01 > 0.4 ) fhReSS[2]->Fill(pt,m); // Both
2129         }
2130
2131         //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2132         for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2133           if((p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)) && (p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))){
2134             for(Int_t iasym=0; iasym < fNAsymCuts; iasym++){
2135               if(a < fAsymCuts[iasym])
2136               {
2137                 Int_t index = ((curCentrBin*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2138                 //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);
2139
2140                 if(index < 0 || index >= ncentr*fNPIDBits*fNAsymCuts) continue ;
2141
2142                 fhRe1     [index]->Fill(pt,m);
2143                 if(fMakeInvPtPlots)fhReInvPt1[index]->Fill(pt,m,1./pt) ;
2144                 if(fFillBadDistHisto){
2145                   if(p1->DistToBad()>0 && p2->DistToBad()>0){
2146                     fhRe2     [index]->Fill(pt,m) ;
2147                     if(fMakeInvPtPlots)fhReInvPt2[index]->Fill(pt,m,1./pt) ;
2148                     if(p1->DistToBad()>1 && p2->DistToBad()>1){
2149                       fhRe3     [index]->Fill(pt,m) ;
2150                       if(fMakeInvPtPlots)fhReInvPt3[index]->Fill(pt,m,1./pt) ;
2151                     }// bad 3
2152                   }// bad2
2153                 }// Fill bad dist histos
2154               }//assymetry cut
2155             }// asymmetry cut loop
2156           }// bad 1
2157         }// pid bit loop
2158
2159         //Fill histograms with opening angle
2160         if(fFillAngleHisto)
2161         {
2162           fhRealOpeningAngle   ->Fill(pt,angle);
2163           fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2164         }
2165
2166         //Fill histograms with pair assymmetry
2167         if(fFillAsymmetryHisto)
2168         {
2169           fhRePtAsym->Fill(pt,a);
2170           if(m > 0.10 && m < 0.17) fhRePtAsymPi0->Fill(pt,a);
2171           if(m > 0.45 && m < 0.65) fhRePtAsymEta->Fill(pt,a);
2172         }
2173
2174         //-------------------------------------------------------
2175         //Get the number of cells needed for multi cut analysis.
2176         //-------------------------------------------------------
2177         Int_t ncell1 = 0;
2178         Int_t ncell2 = 0;
2179         if(fMultiCutAna || (IsDataMC() && fMultiCutAnaSim))
2180         {
2181           AliVEvent * event = GetReader()->GetInputEvent();
2182           if(event){
2183             for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++)
2184             {
2185               AliVCluster *cluster = event->GetCaloCluster(iclus);
2186
2187               Bool_t is = kFALSE;
2188               if     (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
2189               else if(fCalorimeter == "PHOS"  && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
2190
2191               if(is){
2192                 if      (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
2193                 else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
2194               } // PHOS or EMCAL cluster as requested in analysis
2195
2196               if(ncell2 > 0 &&  ncell1 > 0) break; // No need to continue the iteration
2197
2198             }
2199             //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
2200           }
2201         }
2202
2203         //---------
2204         // MC data
2205         //---------
2206         //Do some MC checks on the origin of the pair, is there any common ancestor and if there is one, who?
2207         if(IsDataMC()) {
2208
2209           if(GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2210              GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2211           {
2212             fhReMCFromConversion->Fill(pt,m);
2213           }
2214           else if(!GetMCAnalysisUtils()->CheckTagBit(p1->GetTag(),AliMCAnalysisUtils::kMCConversion) &&
2215                   !GetMCAnalysisUtils()->CheckTagBit(p2->GetTag(),AliMCAnalysisUtils::kMCConversion))
2216           {
2217             fhReMCFromNotConversion->Fill(pt,m);
2218           }
2219           else
2220           {
2221             fhReMCFromMixConversion->Fill(pt,m);
2222           }
2223
2224           FillMCVersusRecDataHistograms(p1->GetLabel(), p2->GetLabel(),p1->Pt(), p2->Pt(),ncell1, ncell2, m, pt, a,deta, dphi);
2225         }
2226
2227         //-----------------------
2228         //Multi cuts analysis
2229         //-----------------------
2230         if(fMultiCutAna)
2231         {
2232           //Histograms for different PID bits selection
2233           for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2234
2235             if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    &&
2236                p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
2237
2238             //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBitsBits+ipid]->GetName());
2239           } // pid bit cut loop
2240
2241           //Several pt,ncell and asymmetry cuts
2242           for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2243             for(Int_t icell=0; icell<fNCellNCuts; icell++){
2244               for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2245                 Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2246                 if(p1->E() >   fPtCuts[ipt]      && p2->E() > fPtCuts[ipt]        &&
2247                    a        <   fAsymCuts[iasym]                                    &&
2248                    ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell]){
2249                   fhRePtNCellAsymCuts[index]->Fill(pt,m) ;
2250                   //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2251                   if(fFillSMCombinations && module1==module2){
2252                     fhRePtNCellAsymCutsSM[module1][index]->Fill(pt,m) ;
2253                   }
2254                 }
2255               }// pid bit cut loop
2256             }// icell loop
2257           }// pt cut loop
2258           if(GetHistogramRanges()->GetHistoTrackMultiplicityBins()){
2259             for(Int_t iasym = 0; iasym < fNAsymCuts; iasym++){
2260               if(a < fAsymCuts[iasym])fhRePtMult[iasym]->Fill(pt,GetTrackMultiplicity(),m) ;
2261             }
2262           }
2263         }// multiple cuts analysis
2264       }// ok if same sm
2265     }// second same event particle
2266   }// first cluster
2267
2268   //-------------------------------------------------------------
2269   // Mixing
2270   //-------------------------------------------------------------
2271   if(DoOwnMix())
2272   {
2273     //Recover events in with same characteristics as the current event
2274
2275     //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2276     if(eventbin < 0) return ;
2277
2278     TList * evMixList=fEventsList[eventbin] ;
2279
2280     if(!evMixList)
2281     {
2282       printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mix event list not available, bin %d \n",eventbin);
2283       return;
2284     }
2285
2286     Int_t nMixed = evMixList->GetSize() ;
2287     for(Int_t ii=0; ii<nMixed; ii++)
2288     {
2289       TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
2290       Int_t nPhot2=ev2->GetEntriesFast() ;
2291       Double_t m = -999;
2292       if(GetDebug() > 1)
2293         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d, centrality bin %d\n", ii, nPhot2, GetEventCentralityBin());
2294
2295       fhEventMixBin->Fill(eventbin) ;
2296
2297       //---------------------------------
2298       //First loop on photons/clusters
2299       //---------------------------------
2300       for(Int_t i1=0; i1<nPhot; i1++){
2301         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
2302         if(fSameSM && GetModuleNumber(p1)!=module1) continue;
2303
2304         //Get kinematics of cluster and (super) module of this cluster
2305         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
2306         module1 = GetModuleNumber(p1);
2307
2308         //---------------------------------
2309         //First loop on photons/clusters
2310         //---------------------------------
2311         for(Int_t i2=0; i2<nPhot2; i2++){
2312           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
2313
2314           //Get kinematics of second cluster and calculate those of the pair
2315           TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
2316           m           = (photon1+photon2).M() ;
2317           Double_t pt = (photon1 + photon2).Pt();
2318           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
2319
2320           //Check if opening angle is too large or too small compared to what is expected
2321           Double_t angle   = photon1.Angle(photon2.Vect());
2322           if(fUseAngleEDepCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle+0.05)){
2323             if(GetDebug() > 2)
2324               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f not in E %f window\n",angle, (photon1+photon2).E());
2325             continue;
2326           }
2327           if(fUseAngleCut && (angle < fAngleCut || angle > fAngleMaxCut)) {
2328             if(GetDebug() > 2)
2329               printf("AliAnaPi0::MakeAnalysisFillHistograms() -Mix pair angle %f < cut %f\n",angle,fAngleCut);
2330             continue;
2331
2332           }
2333
2334           if(GetDebug() > 2)
2335             printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
2336                    p1->Pt(), p2->Pt(), pt,m,a);
2337
2338           //In case we want only pairs in same (super) module, check their origin.
2339           module2 = GetModuleNumber(p2);
2340
2341           //-------------------------------------------------------------------------------------------------
2342           //Fill module dependent histograms, put a cut on assymmetry on the first available cut in the array
2343           //-------------------------------------------------------------------------------------------------
2344           if(a < fAsymCuts[0] && fFillSMCombinations){
2345             if(module1==module2 && module1 >=0 && module1<fNModules)
2346               fhMiMod[module1]->Fill(pt,m) ;
2347
2348             if(fCalorimeter=="EMCAL"){
2349
2350               // Same sector
2351               Int_t j=0;
2352               for(Int_t i = 0; i < fNModules/2; i++){
2353                 j=2*i;
2354                 if((module1==j && module2==j+1) || (module1==j+1 && module2==j)) fhMiSameSectorEMCALMod[i]->Fill(pt,m) ;
2355               }
2356
2357               // Same side
2358               for(Int_t i = 0; i < fNModules-2; i++){
2359                 if((module1==i && module2==i+2) || (module1==i+2 && module2==i)) fhMiSameSideEMCALMod[i]->Fill(pt,m);
2360               }
2361             }//EMCAL
2362             else {//PHOS
2363               if((module1==0 && module2==1) || (module1==1 && module2==0)) fhMiDiffPHOSMod[0]->Fill(pt,m) ;
2364               if((module1==0 && module2==2) || (module1==2 && module2==0)) fhMiDiffPHOSMod[1]->Fill(pt,m) ;
2365               if((module1==1 && module2==2) || (module1==2 && module2==1)) fhMiDiffPHOSMod[2]->Fill(pt,m) ;
2366             }//PHOS
2367
2368
2369           }
2370
2371           Bool_t ok = kTRUE;
2372           if(fSameSM && module1!=module2) ok=kFALSE;
2373           if(ok){
2374
2375             //Check if one of the clusters comes from a conversion
2376             if(fCheckConversion){
2377               if     (p1->IsTagged() && p2->IsTagged()) fhMiConv2->Fill(pt,m);
2378               else if(p1->IsTagged() || p2->IsTagged()) fhMiConv ->Fill(pt,m);
2379             }
2380             //Fill histograms for different bad channel distance, centrality, assymmetry cut and pid bit
2381             for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
2382               if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton)))
2383               {
2384                 for(Int_t iasym=0; iasym < fNAsymCuts; iasym++)
2385                 {
2386                   if(a < fAsymCuts[iasym])
2387                   {
2388                     Int_t index = ((GetEventCentralityBin()*fNPIDBits)+ipid)*fNAsymCuts + iasym;
2389
2390                     if(index < 0 || index >= curCentrBin*fNPIDBits*fNAsymCuts) continue ;
2391
2392                     fhMi1     [index]->Fill(pt,m) ;
2393                     if(fMakeInvPtPlots)fhMiInvPt1[index]->Fill(pt,m,1./pt) ;
2394                     if(fFillBadDistHisto)
2395                     {
2396                       if(p1->DistToBad()>0 && p2->DistToBad()>0)
2397                       {
2398                         fhMi2     [index]->Fill(pt,m) ;
2399                         if(fMakeInvPtPlots)fhMiInvPt2[index]->Fill(pt,m,1./pt) ;
2400                         if(p1->DistToBad()>1 && p2->DistToBad()>1)
2401                         {
2402                           fhMi3     [index]->Fill(pt,m) ;
2403                           if(fMakeInvPtPlots)fhMiInvPt3[index]->Fill(pt,m,1./pt) ;
2404                         }
2405                       }
2406                     }// Fill bad dist histo
2407                   }//Asymmetry cut
2408                 }// Asymmetry loop
2409               }//PID cut
2410             }//loop for histograms
2411
2412             //-----------------------
2413             //Multi cuts analysis
2414             //-----------------------
2415             if(fMultiCutAna){
2416               //Several pt,ncell and asymmetry cuts
2417
2418               for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
2419                 for(Int_t icell=0; icell<fNCellNCuts; icell++){
2420                   for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
2421                     Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
2422                     if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        &&
2423                        a        <   fAsymCuts[iasym]                                    //&&
2424                        //p1->GetBtag() >=  fCellNCuts[icell] && p2->GetBtag() >= fCellNCuts[icell] // trick, correct it.
2425                        ){
2426                       fhMiPtNCellAsymCuts[index]->Fill(pt,m) ;
2427                       //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
2428                     }
2429                   }// pid bit cut loop
2430                 }// icell loop
2431               }// pt cut loop
2432             } // Multi cut ana
2433
2434             //Fill histograms with opening angle
2435             if(fFillAngleHisto){
2436               fhMixedOpeningAngle   ->Fill(pt,angle);
2437               fhMixedCosOpeningAngle->Fill(pt,TMath::Cos(angle));
2438             }
2439
2440           }//ok
2441         }// second cluster loop
2442       }//first cluster loop
2443     }//loop on mixed events
2444
2445     //--------------------------------------------------------
2446     //Add the current event to the list of events for mixing
2447     //--------------------------------------------------------
2448     TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
2449     //Add current event to buffer and Remove redundant events
2450     if(currentEvent->GetEntriesFast()>0){
2451       evMixList->AddFirst(currentEvent) ;
2452       currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
2453       if(evMixList->GetSize() >= GetNMaxEvMix())
2454       {
2455         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
2456         evMixList->RemoveLast() ;
2457         delete tmp ;
2458       }
2459     }
2460     else{ //empty event
2461       delete currentEvent ;
2462       currentEvent=0 ;
2463     }
2464   }// DoOwnMix
2465
2466 }
2467
2468 //________________________________________________________________________
2469 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
2470 {
2471   // retieves the event index and checks the vertex
2472   //    in the mixed buffer returns -2 if vertex NOK
2473   //    for normal events   returns 0 if vertex OK and -1 if vertex NOK
2474
2475   Int_t evtIndex = -1 ;
2476   if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
2477
2478     if (GetMixedEvent()){
2479
2480       evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
2481       GetVertex(vert,evtIndex);
2482
2483       if(TMath::Abs(vert[2])> GetZvertexCut())
2484         evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2485     } else {// Single event
2486
2487       GetVertex(vert);
2488
2489       if(TMath::Abs(vert[2])> GetZvertexCut())
2490         evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
2491       else
2492         evtIndex = 0 ;
2493     }
2494   }//No MC reader
2495   else {
2496     evtIndex = 0;
2497     vert[0] = 0. ;
2498     vert[1] = 0. ;
2499     vert[2] = 0. ;
2500   }
2501
2502   return evtIndex ;
2503 }
2504