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