PWGGA/CaloTrackCorrelations/AliAnaPi0.cxx

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