PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.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 /* $Id: $ */
  16
  17 //_________________________________________________________________________
  18 // Class for the analysis of particle - hadron correlations
  19 // Particle (for example direct gamma) must be found in a previous analysis
  20 //-- Author: Gustavo Conesa (LNF-INFN)
  21
  22 //  Modified by Yaxian Mao:
  23 // 1. add the UE subtraction for corrlation study
  24 // 2. change the correlation variable
  25 //////////////////////////////////////////////////////////////////////////////
  26
  27
  28 // --- ROOT system ---
  29 #include "TH2F.h"
  30 #include "TClonesArray.h"
  31 #include "TClass.h"
  32 #include "TMath.h"
  33
  34 //---- ANALYSIS system ----
  35 #include "AliNeutralMesonSelection.h"
  36 #include "AliAnaParticleHadronCorrelation.h"
  37 #include "AliCaloTrackReader.h"
  38 #include "AliCaloPID.h"
  39 #include "AliAODPWG4ParticleCorrelation.h"
  40 #include "AliFidutialCut.h"
  41 #include "AliAODTrack.h"
  42 #include "AliAODCaloCluster.h"
  43 #include "AliMCAnalysisUtils.h"
  44 #include "TParticle.h"
  45 #include "AliStack.h"
  46 #include "AliAODMCParticle.h"
  47
  48 ClassImp(AliAnaParticleHadronCorrelation)
  49
  50
  51 //____________________________________________________________________________
  52   AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation():
  53     AliAnaPartCorrBaseClass(),
  54     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fSelectIsolated(0),
  55     fMakeSeveralUE(0),  fUeDeltaPhiMaxCut(0.), fUeDeltaPhiMinCut(0.),
  56     fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0),
  57     fhDeltaPhiCharged(0), fhDeltaPhiNeutral(0),
  58     fhDeltaEtaCharged(0), fhDeltaEtaNeutral(0),
  59     fhDeltaPhiChargedPt(0), fhDeltaPhiNeutralPt(0),
  60     fhDeltaPhiUeChargedPt(0), fhDeltaPhiUeNeutralPt(0),
  61     fhPtImbalanceNeutral(0), fhPtImbalanceCharged(0),
  62     fhPtImbalanceUeCharged(0),fhPtImbalanceUeNeutral(0),
  63     fhPtHbpCharged(0), fhPtHbpUeCharged(0),
  64     fhPtHbpNeutral(0), fhPtHbpUeNeutral(0),
  65     fhDeltaPhiUeLeftCharged(0),fhDeltaPhiUeRightCharged(0),
  66     fhDeltaPhiUeLeftNeutral(0),fhDeltaPhiUeRightNeutral(0),
  67     fhPtImbalanceUeLeftCharged(0),fhPtImbalanceUeRightCharged(0),
  68     fhPtImbalanceUeLeftNeutral(0),fhPtImbalanceUeRightNeutral(0),
  69     fhPtHbpUeLeftCharged(0),fhPtHbpUeRightCharged(0),
  70     fhPtHbpUeLeftNeutral(0),fhPtHbpUeRightNeutral(0)
  71 {
  72   //Default Ctor
  73
  74   //Initialize parameters
  75   InitParameters();
  76 }
  77
  78 //____________________________________________________________________________
  79 AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation(const AliAnaParticleHadronCorrelation & g):
  80   AliAnaPartCorrBaseClass(g),
  81   fDeltaPhiMaxCut(g.fDeltaPhiMaxCut), fDeltaPhiMinCut(g.fDeltaPhiMinCut),
  82   fSelectIsolated(g.fSelectIsolated),
  83   fMakeSeveralUE(g.fMakeSeveralUE),  fUeDeltaPhiMaxCut(g.fUeDeltaPhiMaxCut),
  84   fUeDeltaPhiMinCut(g.fUeDeltaPhiMinCut),
  85   fhPhiCharged(g.fhPhiCharged), fhPhiNeutral(g.fhPhiNeutral),
  86   fhEtaCharged(g.fhEtaCharged), fhEtaNeutral(g.fhEtaNeutral),
  87   fhDeltaPhiCharged(g.fhDeltaPhiCharged),
  88   fhDeltaPhiNeutral(g.fhDeltaPhiNeutral),
  89   fhDeltaEtaCharged(g.fhDeltaEtaCharged),
  90   fhDeltaEtaNeutral(g.fhDeltaEtaNeutral),
  91   fhDeltaPhiChargedPt(g.fhDeltaPhiChargedPt),
  92   fhDeltaPhiNeutralPt(g.fhDeltaPhiNeutralPt),
  93   fhDeltaPhiUeChargedPt(g.fhDeltaPhiUeChargedPt),
  94   fhDeltaPhiUeNeutralPt(g.fhDeltaPhiUeNeutralPt),
  95   fhPtImbalanceNeutral(g.fhPtImbalanceNeutral),
  96   fhPtImbalanceCharged(g.fhPtImbalanceCharged),
  97   fhPtImbalanceUeCharged(g.fhPtImbalanceUeCharged),
  98   fhPtImbalanceUeNeutral(g.fhPtImbalanceUeNeutral),
  99   fhPtHbpCharged(g.fhPtHbpCharged),
 100   fhPtHbpUeCharged(g.fhPtHbpUeCharged),
 101   fhPtHbpNeutral(g.fhPtHbpNeutral),
 102   fhPtHbpUeNeutral(g.fhPtHbpUeNeutral),
 103   fhDeltaPhiUeLeftCharged(g.fhDeltaPhiUeLeftCharged),
 104   fhDeltaPhiUeRightCharged(g.fhDeltaPhiUeRightCharged),
 105   fhDeltaPhiUeLeftNeutral(g.fhDeltaPhiUeLeftNeutral),
 106   fhDeltaPhiUeRightNeutral(g.fhDeltaPhiUeRightNeutral),
 107   fhPtImbalanceUeLeftCharged(g.fhPtImbalanceUeLeftCharged),
 108   fhPtImbalanceUeRightCharged(g.fhPtImbalanceUeRightCharged),
 109   fhPtImbalanceUeLeftNeutral(g.fhPtImbalanceUeLeftNeutral),
 110   fhPtImbalanceUeRightNeutral(g.fhPtImbalanceUeRightNeutral),
 111   fhPtHbpUeLeftCharged(g.fhPtHbpUeLeftCharged),
 112   fhPtHbpUeRightCharged(g.fhPtHbpUeRightCharged),
 113   fhPtHbpUeLeftNeutral(g.fhPtHbpUeLeftNeutral),
 114   fhPtHbpUeRightNeutral(g.fhPtHbpUeRightNeutral)
 115 {
 116   // cpy ctor
 117
 118 }
 119
 120 //_________________________________________________________________________
 121 AliAnaParticleHadronCorrelation & AliAnaParticleHadronCorrelation::operator = (const AliAnaParticleHadronCorrelation & source)
 122 {
 123   // assignment operator
 124
 125   if(this == &source)return *this;
 126   ((AliAnaPartCorrBaseClass *)this)->operator=(source);
 127
 128   fDeltaPhiMaxCut = source.fDeltaPhiMaxCut ;
 129   fDeltaPhiMinCut = source.fDeltaPhiMinCut ;
 130   fSelectIsolated = source.fSelectIsolated ;
 131   fMakeSeveralUE = source.fMakeSeveralUE ;
 132   fUeDeltaPhiMaxCut = source.fUeDeltaPhiMaxCut ;
 133   fUeDeltaPhiMinCut = source.fUeDeltaPhiMinCut ;
 134   fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ;
 135   fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ;
 136   fhDeltaPhiCharged = source.fhDeltaPhiCharged ;
 137   fhDeltaPhiNeutral = source.fhDeltaPhiNeutral ;
 138   fhDeltaEtaCharged = source.fhDeltaEtaCharged ;
 139   fhDeltaEtaNeutral = source.fhDeltaEtaNeutral ;
 140   fhDeltaPhiChargedPt = source.fhDeltaPhiChargedPt ;
 141   fhDeltaPhiNeutralPt = source.fhDeltaPhiNeutralPt ;
 142   fhDeltaPhiUeChargedPt = source.fhDeltaPhiUeChargedPt ;
 143   fhDeltaPhiUeNeutralPt = source.fhDeltaPhiUeNeutralPt ;
 144   fhPtImbalanceNeutral = source.fhPtImbalanceNeutral ;
 145   fhPtImbalanceCharged = source.fhPtImbalanceCharged ;
 146   fhPtImbalanceUeCharged = source.fhPtImbalanceUeCharged ;
 147   fhPtImbalanceUeNeutral = source.fhPtImbalanceUeNeutral ;
 148   fhPtHbpCharged = source.fhPtHbpCharged ;
 149   fhPtHbpUeCharged = source.fhPtHbpUeCharged;
 150   fhPtHbpNeutral = source.fhPtHbpNeutral ;
 151   fhPtHbpUeNeutral = source.fhPtHbpUeNeutral;
 152   fhDeltaPhiUeLeftCharged = source.fhDeltaPhiUeLeftCharged ;
 153   fhDeltaPhiUeRightCharged = source.fhDeltaPhiUeRightCharged ;
 154   fhDeltaPhiUeLeftNeutral = source.fhDeltaPhiUeLeftNeutral ;
 155   fhDeltaPhiUeRightNeutral = source.fhDeltaPhiUeRightNeutral ;
 156   fhPtImbalanceUeLeftCharged = source.fhPtImbalanceUeLeftCharged ;
 157   fhPtImbalanceUeRightCharged = source.fhPtImbalanceUeRightCharged ;
 158   fhPtImbalanceUeLeftNeutral = source.fhPtImbalanceUeLeftNeutral ;
 159   fhPtImbalanceUeRightNeutral = source.fhPtImbalanceUeRightNeutral ;
 160   fhPtHbpUeLeftCharged = source.fhPtHbpUeLeftCharged ;
 161   fhPtHbpUeRightCharged = source.fhPtHbpUeRightCharged ;
 162   fhPtHbpUeLeftNeutral = source.fhPtHbpUeLeftNeutral ;
 163   fhPtHbpUeRightNeutral = source.fhPtHbpUeRightNeutral ;
 164
 165   return *this;
 166
 167 }
 168
 169 //________________________________________________________________________
 170 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
 171 {
 172
 173   // Create histograms to be saved in output file and
 174   // store them in fOutputContainer
 175   TList * outputContainer = new TList() ;
 176   outputContainer->SetName("CorrelationHistos") ;
 177
 178   Int_t nptbins  = GetHistoNPtBins();
 179   Int_t nphibins = GetHistoNPhiBins();
 180   Int_t netabins = GetHistoNEtaBins();
 181   Float_t ptmax  = GetHistoPtMax();
 182   Float_t phimax = GetHistoPhiMax();
 183   Float_t etamax = GetHistoEtaMax();
 184   Float_t ptmin  = GetHistoPtMin();
 185   Float_t phimin = GetHistoPhiMin();
 186   Float_t etamin = GetHistoEtaMin();
 187
 188   //Correlation with charged hadrons
 189   if(GetReader()->IsCTSSwitchedOn()) {
 190     fhPhiCharged  = new TH2F
 191       ("PhiCharged","#phi_{h^{#pm}}  vs p_{T #pm}",
 192        nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 193     fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
 194     fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
 195
 196     fhEtaCharged  = new TH2F
 197       ("EtaCharged","#eta_{h^{#pm}}  vs p_{T #pm}",
 198        nptbins,ptmin,ptmax,netabins,etamin,etamax);
 199     fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
 200     fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
 201
 202     fhDeltaPhiCharged  = new TH2F
 203       ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
 204        nptbins,ptmin,ptmax,700,-2.,5.);
 205     fhDeltaPhiCharged->SetYTitle("#Delta #phi");
 206     fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
 207
 208     fhDeltaPhiChargedPt  = new TH2F
 209       ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
 210        nptbins,ptmin,ptmax,700,-2.,5.);
 211     fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
 212     fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 213
 214     fhDeltaPhiUeChargedPt  = new TH2F
 215       ("DeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
 216        nptbins,ptmin,ptmax,700,-2.,5.);
 217     fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
 218     fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 219
 220     fhDeltaEtaCharged  = new TH2F
 221       ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
 222        nptbins,ptmin,ptmax,200,-2,2);
 223     fhDeltaEtaCharged->SetYTitle("#Delta #eta");
 224     fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
 225
 226     fhPtImbalanceCharged  =
 227       new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
 228                nptbins,ptmin,ptmax,1000,0.,2.);
 229     fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
 230     fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
 231
 232     fhPtImbalanceUeCharged  =
 233       new TH2F("CorrelationUeCharged","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger}",
 234                nptbins,ptmin,ptmax,1000,0.,2.);
 235     fhPtImbalanceUeCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
 236     fhPtImbalanceUeCharged->SetXTitle("p_{T trigger}");
 237
 238    fhPtHbpCharged  =
 239       new TH2F("HbpCharged","#xi = ln(1/x_{E}) with charged hadrons",
 240                nptbins,ptmin,ptmax,1000,0.,10.);
 241     fhPtHbpCharged->SetYTitle("ln(1/x_{E})");
 242     fhPtHbpCharged->SetXTitle("p_{T trigger}");
 243
 244     fhPtHbpUeCharged  =
 245       new TH2F("HbpUeCharged","#xi = ln(1/x_{E}) with charged hadrons",
 246                nptbins,ptmin,ptmax,1000,0.,10.);
 247     fhPtHbpUeCharged->SetYTitle("ln(1/x_{E})");
 248     fhPtHbpUeCharged->SetXTitle("p_{T trigger}");
 249
 250     outputContainer->Add(fhPhiCharged) ;
 251     outputContainer->Add(fhEtaCharged) ;
 252     outputContainer->Add(fhDeltaPhiCharged) ;
 253     outputContainer->Add(fhDeltaEtaCharged) ;
 254     outputContainer->Add(fhDeltaPhiChargedPt) ;
 255     outputContainer->Add(fhDeltaPhiUeChargedPt) ;
 256     outputContainer->Add(fhPtImbalanceCharged) ;
 257     outputContainer->Add(fhPtImbalanceUeCharged) ;
 258     outputContainer->Add(fhPtHbpCharged) ;
 259     outputContainer->Add(fhPtHbpUeCharged) ;
 260     if(fMakeSeveralUE){
 261       fhDeltaPhiUeLeftCharged  = new TH2F
 262         ("DeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
 263          nptbins,ptmin,ptmax,700,-2.,5.);
 264       fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
 265       fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 266       outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
 267
 268       fhDeltaPhiUeRightCharged  = new TH2F
 269         ("DeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
 270          nptbins,ptmin,ptmax,700,-2.,5.);
 271       fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
 272       fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 273       outputContainer->Add(fhDeltaPhiUeRightCharged) ;
 274
 275       fhPtImbalanceUeLeftCharged  =
 276         new TH2F("CorrelationUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",
 277                  nptbins,ptmin,ptmax,1000,0.,2.);
 278       fhPtImbalanceUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
 279       fhPtImbalanceUeLeftCharged->SetXTitle("p_{T trigger}");
 280       outputContainer->Add(fhPtImbalanceUeLeftCharged) ;
 281
 282       fhPtImbalanceUeRightCharged  =
 283         new TH2F("CorrelationUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",
 284                  nptbins,ptmin,ptmax,1000,0.,2.);
 285       fhPtImbalanceUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
 286       fhPtImbalanceUeRightCharged->SetXTitle("p_{T trigger}");
 287       outputContainer->Add(fhPtImbalanceUeRightCharged) ;
 288
 289                 fhPtHbpUeLeftCharged  =
 290                 new TH2F("HbpUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",
 291                                  nptbins,ptmin,ptmax,1000,0.,10.);
 292                 fhPtHbpUeLeftCharged->SetYTitle("ln(1/x_{E})");
 293                 fhPtHbpUeLeftCharged->SetXTitle("p_{T trigger}");
 294                 outputContainer->Add(fhPtHbpUeLeftCharged) ;
 295
 296                 fhPtHbpUeRightCharged  =
 297                 new TH2F("HbpUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",
 298                                  nptbins,ptmin,ptmax,1000,0.,10.);
 299                 fhPtHbpUeRightCharged->SetYTitle("ln(1/x_{E})");
 300                 fhPtHbpUeRightCharged->SetXTitle("p_{T trigger}");
 301                 outputContainer->Add(fhPtHbpUeRightCharged) ;
 302
 303     }
 304   }  //Correlation with charged hadrons
 305
 306   //Correlation with neutral hadrons
 307   if(GetReader()->IsEMCALSwitchedOn() || GetReader()->IsPHOSSwitchedOn()){
 308
 309     fhPhiNeutral  = new TH2F
 310       ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T #pi^{0}}",
 311        nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 312     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
 313     fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
 314
 315     fhEtaNeutral  = new TH2F
 316       ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T #pi^{0}}",
 317        nptbins,ptmin,ptmax,netabins,etamin,etamax);
 318     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
 319     fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
 320
 321     fhDeltaPhiNeutral  = new TH2F
 322       ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
 323        nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 324     fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
 325     fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
 326
 327     fhDeltaPhiNeutralPt  = new TH2F
 328       ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
 329        nptbins,ptmin,ptmax,700,-2.,5.);
 330     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
 331     fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0} (GeV/c)");
 332
 333     fhDeltaPhiUeNeutralPt  = new TH2F
 334       ("DeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
 335        nptbins,ptmin,ptmax,700,-2.,5.);
 336     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
 337     fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0} (GeV/c)");
 338
 339     fhDeltaEtaNeutral  = new TH2F
 340       ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
 341        nptbins,ptmin,ptmax,200,-2,2);
 342     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
 343     fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
 344
 345     fhPtImbalanceNeutral  =
 346       new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
 347                nptbins,ptmin,ptmax,1000,0.,2.);
 348     fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
 349     fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
 350
 351     fhPtImbalanceUeNeutral  =
 352       new TH2F("CorrelationUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
 353                nptbins,ptmin,ptmax,1000,0.,2.);
 354     fhPtImbalanceUeNeutral->SetYTitle("z_{trigger #pi^{0}}");
 355     fhPtImbalanceUeNeutral->SetXTitle("p_{T trigger}");
 356
 357         fhPtHbpNeutral  =
 358       new TH2F("HbpNeutral","#xi = ln(1/x_{E}) with neutral particles",
 359                            nptbins,ptmin,ptmax,1000,0.,10.);
 360         fhPtHbpNeutral->SetYTitle("ln(1/x_{E})");
 361         fhPtHbpNeutral->SetXTitle("p_{T trigger}");
 362
 363         fhPtHbpUeNeutral  =
 364       new TH2F("HbpUeNeutral","#xi = ln(1/x_{E}) with neutral particles",
 365                            nptbins,ptmin,ptmax,1000,0.,10.);
 366         fhPtHbpUeNeutral->SetYTitle("ln(1/x_{E})");
 367         fhPtHbpUeNeutral->SetXTitle("p_{T trigger}");
 368
 369
 370
 371     outputContainer->Add(fhPhiNeutral) ;
 372     outputContainer->Add(fhEtaNeutral) ;
 373     outputContainer->Add(fhDeltaPhiNeutral) ;
 374     outputContainer->Add(fhDeltaPhiNeutralPt) ;
 375     outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
 376     outputContainer->Add(fhDeltaEtaNeutral) ;
 377     outputContainer->Add(fhPtImbalanceNeutral) ;
 378     outputContainer->Add(fhPtImbalanceUeNeutral) ;
 379         outputContainer->Add(fhPtHbpNeutral) ;
 380         outputContainer->Add(fhPtHbpUeNeutral) ;
 381
 382     if(fMakeSeveralUE){
 383       fhDeltaPhiUeLeftNeutral  = new TH2F
 384         ("DeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",
 385          nptbins,ptmin,ptmax,700,-2.,5.);
 386       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
 387       fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
 388       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
 389
 390       fhDeltaPhiUeRightNeutral  = new TH2F
 391         ("DeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",
 392          nptbins,ptmin,ptmax,700,-2.,5.);
 393       fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
 394       fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
 395       outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
 396
 397       fhPtImbalanceUeLeftNeutral  =
 398         new TH2F("CorrelationUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
 399                  nptbins,ptmin,ptmax,1000,0.,2.);
 400       fhPtImbalanceUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
 401       fhPtImbalanceUeLeftNeutral->SetXTitle("p_{T trigger}");
 402       outputContainer->Add(fhPtImbalanceUeLeftNeutral) ;
 403
 404       fhPtImbalanceUeRightNeutral  =
 405         new TH2F("CorrelationUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
 406                  nptbins,ptmin,ptmax,1000,0.,2.);
 407       fhPtImbalanceUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
 408       fhPtImbalanceUeRightNeutral->SetXTitle("p_{T trigger}");
 409       outputContainer->Add(fhPtImbalanceUeRightNeutral) ;
 410
 411         fhPtHbpUeLeftNeutral  =
 412         new TH2F("HbpUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",
 413                         nptbins,ptmin,ptmax,1000,0.,10.);
 414         fhPtHbpUeLeftNeutral->SetYTitle("ln(1/x_{E})");
 415         fhPtHbpUeLeftNeutral->SetXTitle("p_{T trigger}");
 416         outputContainer->Add(fhPtHbpUeLeftNeutral) ;
 417
 418         fhPtHbpUeRightNeutral  =
 419         new TH2F("HbpUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",
 420                         nptbins,ptmin,ptmax,1000,0.,10.);
 421         fhPtHbpUeRightNeutral->SetYTitle("ln(1/x_{E})");
 422         fhPtHbpUeRightNeutral->SetXTitle("p_{T trigger}");
 423         outputContainer->Add(fhPtHbpUeRightNeutral) ;
 424
 425     }
 426
 427     //Keep neutral meson selection histograms if requiered
 428     //Setting done in AliNeutralMesonSelection
 429     if(GetNeutralMesonSelection()){
 430       TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
 431       if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
 432         for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
 433     }
 434
 435   }//Correlation with neutral hadrons
 436
 437   return outputContainer;
 438
 439 }
 440
 441 //____________________________________________________________________________
 442 void AliAnaParticleHadronCorrelation::InitParameters()
 443 {
 444
 445   //Initialize the parameters of the analysis.
 446   SetInputAODName("PWG4Particle");
 447   SetAODObjArrayName("Hadrons");
 448   AddToHistogramsName("AnaHadronCorr_");
 449
 450   //Correlation with neutrals
 451   //SetOutputAODClassName("AliAODPWG4Particle");
 452   //SetOutputAODName("Pi0Correlated");
 453
 454   SetPtCutRange(0.,300);
 455   fDeltaPhiMinCut = 1.5 ;
 456   fDeltaPhiMaxCut = 4.5 ;
 457   fSelectIsolated = kFALSE;
 458   fMakeSeveralUE = kFALSE;
 459   fUeDeltaPhiMinCut = 1. ;
 460   fUeDeltaPhiMaxCut = 1.5 ;
 461 }
 462
 463 //__________________________________________________________________
 464 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
 465 {
 466
 467   //Print some relevant parameters set for the analysis
 468   if(! opt)
 469     return;
 470
 471   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 472   AliAnaPartCorrBaseClass::Print(" ");
 473
 474   printf("Phi trigger particle-Hadron      <     %3.2f\n", fDeltaPhiMaxCut) ;
 475   printf("Phi trigger particle-Hadron      >     %3.2f\n", fDeltaPhiMinCut) ;
 476   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
 477   printf("Phi trigger particle-UeHadron    <    %3.2f\n", fUeDeltaPhiMaxCut) ;
 478   printf("Phi trigger particle-UeHadron    >    %3.2f\n", fUeDeltaPhiMinCut) ;
 479   printf("Several UE?  %d\n", fMakeSeveralUE) ;
 480 }
 481
 482 //____________________________________________________________________________
 483 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
 484 {
 485   //Particle-Hadron Correlation Analysis, fill AODs
 486
 487   if(!GetInputAODBranch()){
 488     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
 489     abort();
 490   }
 491
 492   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
 493         printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
 494         abort();
 495   }
 496
 497   if(GetDebug() > 1){
 498     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");
 499     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
 500     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
 501     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
 502     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
 503   }
 504
 505   //Loop on stored AOD particles, trigger
 506   Int_t naod = GetInputAODBranch()->GetEntriesFast();
 507   for(Int_t iaod = 0; iaod < naod ; iaod++){
 508     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 509
 510         //Make correlation with charged hadrons
 511     if(GetReader()->IsCTSSwitchedOn() )
 512       MakeChargedCorrelation(particle, GetAODCTS(),kFALSE);
 513
 514     //Make correlation with neutral pions
 515     //Trigger particle in PHOS, correlation with EMCAL
 516     if(particle->GetDetector()=="PHOS" && GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntriesFast() > 0)
 517       MakeNeutralCorrelationFillAOD(particle, GetAODEMCAL(),"EMCAL");
 518     //Trigger particle in EMCAL, correlation with PHOS
 519     else if(particle->GetDetector()=="EMCAL" && GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntriesFast() > 0)
 520       MakeNeutralCorrelationFillAOD(particle, GetAODPHOS(),"PHOS");
 521     //Trigger particle in CTS, correlation with PHOS, EMCAL and CTS
 522     else if(particle->GetDetector()=="CTS" ){
 523       if(GetReader()->IsPHOSSwitchedOn() && GetAODPHOS()->GetEntriesFast() > 0)
 524         MakeNeutralCorrelationFillAOD(particle, GetAODPHOS(),"PHOS");
 525       if(GetReader()->IsEMCALSwitchedOn() && GetAODEMCAL()->GetEntriesFast() > 0)
 526         MakeNeutralCorrelationFillAOD(particle, GetAODEMCAL(),"EMCAL");
 527     }
 528
 529
 530   }//Aod branch loop
 531
 532   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
 533
 534 }
 535
 536 //____________________________________________________________________________
 537 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
 538 {
 539   //Particle-Hadron Correlation Analysis, fill histograms
 540
 541   if(!GetInputAODBranch()){
 542     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
 543     abort();
 544   }
 545
 546   if(GetDebug() > 1){
 547     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
 548     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
 549   }
 550
 551   //Loop on stored AOD particles
 552   Int_t naod = GetInputAODBranch()->GetEntriesFast();
 553   for(Int_t iaod = 0; iaod < naod ; iaod++){
 554     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 555
 556     //check if the particle is isolated or if we want to take the isolation into account
 557     if(OnlyIsolated() && !particle->IsIsolated()) continue;
 558
 559     //Make correlation with charged hadrons
 560     TObjArray * reftracks   = particle->GetObjArray(GetAODObjArrayName()+"Tracks");
 561     if(reftracks){
 562       if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Particle %d, In Track Refs  entries %d\n", iaod, reftracks->GetEntriesFast());
 563       if(reftracks->GetEntriesFast() > 0) MakeChargedCorrelation(particle, reftracks,kTRUE);
 564     }
 565
 566     //Make correlation with neutral pions
 567     if(GetOutputAODBranch() && GetOutputAODBranch()->GetEntriesFast() > 0){
 568       if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Particle %d, In Cluster Refs entries %d\n",iaod, GetOutputAODBranch()->GetEntriesFast());
 569       MakeNeutralCorrelationFillHistograms(particle);
 570     }
 571
 572   }//Aod branch loop
 573
 574   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
 575
 576 }
 577
 578 //____________________________________________________________________________
 579 void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle, TObjArray* const pl, const Bool_t bFillHisto)
 580 {
 581    // Charged Hadron Correlation Analysis
 582    if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
 583
 584    Double_t ptTrig  = aodParticle->Pt();
 585    Double_t pxTrig  = aodParticle->Px();
 586    Double_t pyTrig  = aodParticle->Py();
 587
 588    Double_t phiTrig = aodParticle->Phi();
 589
 590    Double_t pt   = -100.;
 591    Double_t px   = -100.;
 592    Double_t py   = -100.;
 593    Double_t rat  = -100.;
 594    Double_t xE  = -100.;
 595    Double_t cosi  = -100.;
 596    Double_t phi  = -100. ;
 597    Double_t eta  = -100. ;
 598    Double_t p[3];
 599
 600     TObjArray * reftracks    =0x0;
 601    if(!bFillHisto)
 602      reftracks    = new TObjArray;
 603
 604    //Track loop, select tracks with good pt, phi and fill AODs or histograms
 605    for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
 606      AliAODTrack * track = (AliAODTrack *) (pl->At(ipr)) ;
 607      track->GetPxPyPz(p) ;
 608      TLorentzVector mom(p[0],p[1],p[2],0);
 609      pt   = mom.Pt();
 610      px   = mom.Px();
 611      py   = mom.Py();
 612      eta  = mom.Eta();
 613      phi  = mom.Phi() ;
 614      if(phi < 0) phi+=TMath::TwoPi();
 615      rat   = pt/ptTrig ;
 616      xE   = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 617      if(xE <0.)xE =-xE;
 618      cosi = TMath::Log(1/xE);
 619     // printf("rat = %f, xE = %f, cosi =%f \n", rat, xE, cosi);
 620      // printf("phi = %f \n", phi);
 621
 622      if(IsFidutialCutOn()){
 623        Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
 624        if(! in ) continue ;
 625      }
 626
 627      //Select only hadrons in pt range
 628      if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
 629
 630      //Selection within angular range
 631      Float_t deltaphi = phiTrig-phi;
 632      if(deltaphi<-TMath::PiOver2()) deltaphi+=TMath::TwoPi();
 633      if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();
 634
 635      if(GetDebug() > 2)
 636        printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Charged hadron: pt %f, phi %f, phi trigger %f. Cuts:  delta phi  %2.2f < %2.2f < %2.2f, pT min %2.2f \n",
 637              pt,phi, phiTrig,fDeltaPhiMinCut, deltaphi, fDeltaPhiMaxCut, GetMinPt());
 638
 639      if(bFillHisto){
 640        // Fill Histograms
 641        fhEtaCharged->Fill(pt,eta);
 642        fhPhiCharged->Fill(pt,phi);
 643        fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
 644        fhDeltaPhiCharged->Fill(ptTrig, deltaphi);
 645
 646        if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
 647
 648        //delta phi cut for correlation
 649        if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {
 650          fhDeltaPhiChargedPt->Fill(pt,deltaphi);
 651          fhPtImbalanceCharged->Fill(ptTrig,rat);
 652          fhPtHbpCharged->Fill(ptTrig,cosi);
 653        }
 654       else {
 655                 fhDeltaPhiUeChargedPt->Fill(pt,deltaphi);
 656                 fhPtImbalanceUeCharged->Fill(ptTrig,rat);
 657                 fhPtHbpUeCharged->Fill(ptTrig,cosi);
 658       }
 659        //several UE calculation
 660       if(fMakeSeveralUE){
 661          if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){
 662            fhDeltaPhiUeLeftCharged->Fill(pt,deltaphi);
 663            fhPtImbalanceUeLeftCharged->Fill(ptTrig,rat);
 664            fhPtHbpUeLeftCharged->Fill(ptTrig,cosi);
 665         }
 666          if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){
 667            fhDeltaPhiUeRightCharged->Fill(pt,deltaphi);
 668            fhPtImbalanceUeRightCharged->Fill(ptTrig,rat);
 669            fhPtHbpUeRightCharged->Fill(ptTrig,cosi);
 670
 671          }
 672       } //several UE calculation
 673
 674      }
 675      else{
 676        reftracks->Add(track);
 677      }//aod particle loop
 678    }// track loop
 679
 680    //Fill AOD with reference tracks, if not filling histograms
 681    if(!bFillHisto && reftracks->GetEntriesFast() > 0) {
 682      reftracks->SetName(GetAODObjArrayName()+"Tracks");
 683      aodParticle->AddObjArray(reftracks);
 684    }
 685
 686 }
 687
 688 //____________________________________________________________________________
 689 void  AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD(AliAODPWG4ParticleCorrelation* const aodParticle,TObjArray* const pl, TString detector)
 690 {
 691   // Neutral Pion Correlation Analysis, find pi0, put them in new output aod, if correlation cuts passed
 692   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Make trigger particle - neutral hadron correlation \n");
 693
 694   if(!NewOutputAOD()){
 695     printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Output aod not created, set AOD class name and branch name in the configuration file, STOP! \n");
 696     abort();
 697   }
 698
 699   Double_t phiTrig = aodParticle->Phi();
 700   Int_t tag = 0;
 701   TLorentzVector gammai;
 702   TLorentzVector gammaj;
 703
 704   //Get vertex for photon momentum calculation
 705   Double_t vertex[]  = {0,0,0} ; //vertex
 706   Double_t vertex2[] = {0,0,0} ; //vertex of second input aod
 707   if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
 708   {
 709          GetReader()->GetVertex(vertex);
 710          if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
 711   }
 712
 713   //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
 714   //Int_t iEvent= GetReader()->GetEventNumber() ;
 715   Int_t nclus = pl->GetEntriesFast();
 716   for(Int_t iclus = 0;iclus < nclus ; iclus ++ ){
 717     AliAODCaloCluster * calo = (AliAODCaloCluster *) (pl->At(iclus)) ;
 718
 719         //Input from second AOD?
 720         Int_t inputi = 0;
 721         if     (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus) inputi = 1 ;
 722         else if(aodParticle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus) inputi = 1;
 723
 724         //Cluster selection, not charged, with photon or pi0 id and in fidutial cut
 725     Int_t pdg=0;
 726         if     (inputi == 0 && !SelectCluster(calo, vertex,  gammai, pdg))  continue ;
 727         else if(inputi == 1 && !SelectCluster(calo, vertex2, gammai, pdg))  continue ;
 728
 729         if(GetDebug() > 2)
 730       printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral cluster in %s: pt %f, phi %f, phi trigger %f. Cuts:  delta phi min %2.2f,  max %2.2f, pT min %2.2f \n",
 731              detector.Data(), gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
 732
 733     //2 gamma overlapped, found with PID
 734     if(pdg == AliCaloPID::kPi0){
 735
 736       //Select only hadrons in pt range
 737       if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt()) continue ;
 738
 739       //Selection within angular range
 740       Float_t phi = gammai.Phi();
 741       if(phi < 0) phi+=TMath::TwoPi();
 742       //Float_t deltaphi = TMath::Abs(phiTrig-phi);
 743       //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
 744
 745       AliAODPWG4Particle pi0 = AliAODPWG4Particle(gammai);
 746       //pi0.SetLabel(calo->GetLabel(0));
 747       pi0.SetPdg(AliCaloPID::kPi0);
 748       pi0.SetDetector(detector);
 749
 750       if(IsDataMC()){
 751             pi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(0),GetReader(),inputi));
 752             if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of candidate %d\n",pi0.GetTag());
 753           }//Work with stack also
 754       //Set the indeces of the original caloclusters
 755       pi0.SetCaloLabel(calo->GetID(),-1);
 756       AddAODParticle(pi0);
 757
 758       if(GetDebug() > 2)
 759         printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Correlated with selected pi0 (pid): pt %f, phi %f\n",pi0.Pt(),pi0.Phi());
 760
 761     }// pdg = 111
 762
 763     //Make invariant mass analysis
 764     else if(pdg == AliCaloPID::kPhoton){
 765       //Search the photon companion in case it comes from  a Pi0 decay
 766       //Apply several cuts to select the good pair;
 767       for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){
 768         AliAODCaloCluster * calo2 = (AliAODCaloCluster *) (pl->At(jclus)) ;
 769
 770         //Input from second AOD?
 771         Int_t inputj = 0;
 772         if     (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= jclus) inputj = 1;
 773         else if(aodParticle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= jclus) inputj = 1;
 774
 775         //Cluster selection, not charged with photon or pi0 id and in fidutial cut
 776         Int_t pdgj=0;
 777         if     (inputj == 0 && !SelectCluster(calo2, vertex,  gammaj, pdgj))  continue ;
 778         else if(inputj == 1 && !SelectCluster(calo2, vertex2, gammaj, pdgj))  continue ;
 779
 780         if(!SelectCluster(calo2,vertex, gammaj, pdgj)) continue ;
 781
 782         if(pdgj == AliCaloPID::kPhoton ){
 783
 784           if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt()) continue ;
 785
 786           //Selection within angular range
 787           Float_t phi = (gammai+gammaj).Phi();
 788           if(phi < 0) phi+=TMath::TwoPi();
 789           //Float_t deltaphi = TMath::Abs(phiTrig-phi);
 790           //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
 791
 792           //Select good pair (aperture and invariant mass)
 793           if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
 794
 795             if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral Hadron Correlation: AOD Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f\n",
 796                                        (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());
 797
 798             TLorentzVector pi0mom = gammai+gammaj;
 799             AliAODPWG4Particle pi0 = AliAODPWG4Particle(pi0mom);
 800             //pi0.SetLabel(calo->GetLabel(0));
 801             pi0.SetPdg(AliCaloPID::kPi0);
 802             pi0.SetDetector(detector);
 803             if(IsDataMC()){
 804               //Check origin of the candidates
 805
 806                   Int_t label1 = calo->GetLabel(0);
 807                   Int_t label2 = calo2->GetLabel(0);
 808               Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), inputi);
 809               Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), inputj);
 810
 811               if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
 812               if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
 813
 814                           //Check if pi0 mother is the same
 815                           if(GetReader()->ReadStack()){
 816                                         TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
 817                                         label1 = mother1->GetFirstMother();
 818                                         //mother1 = GetMCStack()->Particle(label1);//pi0
 819
 820                                         TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
 821                                         label2 = mother2->GetFirstMother();
 822                                         //mother2 = GetMCStack()->Particle(label2);//pi0
 823                           }
 824                           else if(GetReader()->ReadAODMCParticles()){
 825                                         AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputi))->At(label1);//photon in kine tree
 826                                         label1 = mother1->GetMother();
 827                                         //mother1 = GetMCStack()->Particle(label1);//pi0
 828                                         AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputj))->At(label2);//photon in kine tree
 829                                         label2 = mother2->GetMother();
 830                                         //mother2 = GetMCStack()->Particle(label2);//pi0
 831                           }
 832
 833                           //printf("mother1 %d, mother2 %d\n",label1,label2);
 834                           if(label1 == label2)
 835                                   GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
 836               }
 837             }//Work with mc information also
 838             pi0.SetTag(tag);
 839             //Set the indeces of the original caloclusters
 840             pi0.SetCaloLabel(calo->GetID(), calo2->GetID());
 841             AddAODParticle(pi0);
 842
 843
 844           }//Pair selected
 845         }//if pair of gammas
 846       }//2nd loop
 847     }// if pdg = 22
 848   }//1st loop
 849
 850   if(GetDebug() > 1)
 851     printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - End, %d pi0's found \n",GetOutputAODBranch()->GetEntriesFast());
 852 }
 853
 854 //____________________________________________________________________________
 855 void  AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms(AliAODPWG4ParticleCorrelation * const aodParticle)
 856 {
 857   // Neutral Pion Correlation Analysis
 858   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistogramS() - Make trigger particle - pi0 correlation, %d pi0's \n",GetOutputAODBranch()->GetEntriesFast());
 859
 860   Double_t pt  = -100.;
 861   Double_t px   = -100.;
 862   Double_t py   = -100.;
 863   Double_t rat = -100.;
 864   Double_t phi = -100.;
 865   Double_t eta = -100.;
 866   Double_t xE  = -100.;
 867   Double_t cosi  = -100.;
 868
 869   Double_t ptTrig  = aodParticle->Pt();
 870   Double_t phiTrig = aodParticle->Phi();
 871   Double_t etaTrig = aodParticle->Eta();
 872   Double_t pxTrig  = aodParticle->Px();
 873   Double_t pyTrig  = aodParticle->Py();
 874
 875
 876   if(!GetOutputAODBranch()){
 877     printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
 878     abort();
 879   }
 880
 881   //Loop on stored AOD pi0
 882   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 883   if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() -  aod branch entries %d\n", naod);
 884   for(Int_t iaod = 0; iaod < naod ; iaod++){
 885     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 886     Int_t pdg = pi0->GetPdg();
 887
 888     if(pdg != AliCaloPID::kPi0) continue;
 889     pt  = pi0->Pt();
 890     px   = pi0->Px();
 891         py   = pi0->Py();
 892     if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
 893
 894     //Selection within angular range
 895     phi = pi0->Phi();
 896     //Float_t deltaphi = TMath::Abs(phiTrig-phi);
 897     //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
 898     Float_t deltaphi = phiTrig-phi;
 899     if(deltaphi<-TMath::PiOver2()) deltaphi+=TMath::TwoPi();
 900     if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();
 901
 902     rat = pt/ptTrig ;
 903     phi = pi0->Phi() ;
 904     eta = pi0->Eta() ;
 905         xE   = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 906         if(xE <0.)xE =-xE;
 907         cosi = TMath::Log(1/xE);
 908
 909     fhEtaNeutral->Fill(pt,eta);
 910     fhPhiNeutral->Fill(pt,phi);
 911     fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);
 912     fhDeltaPhiNeutral->Fill(ptTrig,deltaphi);
 913
 914        //delta phi cut for correlation
 915        if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {
 916          fhDeltaPhiNeutralPt->Fill(pt,deltaphi);
 917          fhPtImbalanceNeutral->Fill(ptTrig,rat);
 918          fhPtHbpNeutral->Fill(ptTrig,cosi);
 919        }
 920       else {
 921         fhDeltaPhiUeNeutralPt->Fill(pt,deltaphi);
 922         fhPtImbalanceUeNeutral->Fill(ptTrig,rat);
 923         fhPtHbpUeNeutral->Fill(ptTrig,cosi);
 924       }
 925        //several UE calculation
 926       if(fMakeSeveralUE){
 927          if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){
 928            fhDeltaPhiUeLeftNeutral->Fill(pt,deltaphi);
 929            fhPtImbalanceUeLeftNeutral->Fill(ptTrig,rat);
 930            fhPtHbpUeLeftNeutral->Fill(ptTrig,cosi);
 931         }
 932          if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){
 933            fhDeltaPhiUeRightNeutral->Fill(pt,deltaphi);
 934            fhPtImbalanceUeRightNeutral->Fill(ptTrig,rat);
 935            fhPtHbpUeRightNeutral->Fill(ptTrig,cosi);
 936          }
 937       } //several UE calculation
 938
 939
 940     if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
 941
 942   }//loop
 943 }
 944
 945
 946 //____________________________________________________________________________
 947 Bool_t  AliAnaParticleHadronCorrelation::SelectCluster(AliAODCaloCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) const {
 948   //Select cluster depending on its pid and acceptance selections
 949
 950   //Skip matched clusters with tracks
 951   if(calo->GetNTracksMatched() > 0) return kFALSE;
 952
 953   TString detector = "";
 954   if(calo->IsPHOSCluster()) detector= "PHOS";
 955   else if(calo->IsEMCALCluster()) detector= "EMCAL";
 956
 957   //Check PID
 958   calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
 959   pdg = AliCaloPID::kPhoton;
 960   if(IsCaloPIDOn()){
 961     //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
 962     //or redo PID, recommended option for EMCal.
 963
 964     if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
 965       pdg = GetCaloPID()->GetPdg(detector,calo->PID(),mom.E());//PID with weights
 966     else
 967       pdg = GetCaloPID()->GetPdg(detector,mom,calo);//PID recalculated
 968
 969     if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);
 970
 971     //If it does not pass pid, skip
 972     if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) {
 973       return kFALSE ;
 974     }
 975   }//PID on
 976
 977   //Check acceptance selection
 978   if(IsFidutialCutOn()){
 979     Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,detector) ;
 980     if(! in ) return kFALSE ;
 981   }
 982
 983   if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
 984
 985   return kTRUE;
 986
 987 }