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