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 // 4. Make decay photon-hadron correlations where decay contribute pi0 mass (2010/09/09)
  27 // 5. fill the pout to extract kt at the end, also to study charge asymmetry(2010/10/06)
  28 //////////////////////////////////////////////////////////////////////////////
  29
  30
  31 // --- ROOT system ---
  32 //#include "TClonesArray.h"
  33 #include "TClass.h"
  34 #include "TMath.h"
  35 #include "TH3D.h"
  36
  37 //---- ANALYSIS system ----
  38 #include "AliNeutralMesonSelection.h"
  39 #include "AliAnaParticleHadronCorrelation.h"
  40 #include "AliCaloTrackReader.h"
  41 #include "AliCaloPID.h"
  42 #include "AliAODPWG4ParticleCorrelation.h"
  43 #include "AliFiducialCut.h"
  44 #include "AliAODTrack.h"
  45 #include "AliVCluster.h"
  46 #include "AliMCAnalysisUtils.h"
  47 #include "TParticle.h"
  48 #include "AliStack.h"
  49 #include "AliAODMCParticle.h"
  50 #include "AliMixedEvent.h"
  51
  52 ClassImp(AliAnaParticleHadronCorrelation)
  53
  54
  55 //____________________________________________________________________________
  56   AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation():
  57     AliAnaPartCorrBaseClass(),
  58     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fSelectIsolated(0),
  59     fMakeSeveralUE(0),  fUeDeltaPhiMaxCut(0.), fUeDeltaPhiMinCut(0.),
  60     fPi0AODBranchName(""),fNeutralCorr(0), fPi0Trigger(0),
  61    // fMultiBin(0),fNZvertBin(0),fNrpBin(0),fZvtxCut(0.),
  62    // fUseSelectEvent(kFALSE),
  63    // fhNclustersNtracks(0), //fhVertex(0),
  64     fhPtLeading(0),fhPhiLeading(0),fhEtaLeading(0),
  65     fhDeltaPhiDeltaEtaCharged(0),
  66     fhPhiCharged(0), fhEtaCharged(0),
  67     fhDeltaPhiCharged(0),
  68     fhDeltaEtaCharged(0),
  69     fhDeltaPhiChargedPt(0),
  70     fhDeltaPhiUeChargedPt(0),
  71     fhPtImbalanceCharged(0),
  72     fhPtImbalanceUeCharged(0),
  73     fhPtImbalancePosCharged(0),fhPtImbalanceNegCharged(0),
  74     fhPtHbpCharged(0), fhPtHbpUeCharged(0),
  75     fhDeltaPhiUeLeftCharged(0),fhDeltaPhiUeRightCharged(0),
  76     fhPtImbalanceUeLeftCharged(0),fhPtImbalanceUeRightCharged(0),
  77     fhPtHbpUeLeftCharged(0),fhPtHbpUeRightCharged(0),
  78     fhPoutPtTrigPtAssoc(0), fhUePoutPtTrigPtAssoc(0),
  79     fhPtTrigCharged(0),
  80     fhTrigDeltaPhiDeltaEtaCharged(0x0), fhTrigCorr(0x0),fhTrigUeCorr(0x0),
  81     fhDeltaPhiDeltaEtaNeutral(0),
  82     fhPhiNeutral(0), fhEtaNeutral(0),
  83     fhDeltaPhiNeutral(0), fhDeltaEtaNeutral(0),
  84     fhDeltaPhiNeutralPt(0),fhDeltaPhiUeNeutralPt(0),
  85     fhPtImbalanceNeutral(0),fhPtImbalanceUeNeutral(0),
  86     fhPtHbpNeutral(0), fhPtHbpUeNeutral(0),
  87     fhDeltaPhiUeLeftNeutral(0),fhDeltaPhiUeRightNeutral(0),
  88     fhPtImbalanceUeLeftNeutral(0),fhPtImbalanceUeRightNeutral(0),
  89     fhPtHbpUeLeftNeutral(0),fhPtHbpUeRightNeutral(0),
  90     fhPtPi0DecayRatio(0),
  91     fhDeltaPhiDecayCharged(0),
  92     fhPtImbalanceDecayCharged(0),
  93     fhDeltaPhiDecayNeutral(0),
  94     fhPtImbalanceDecayNeutral(0)
  95 {
  96   //Default Ctor
  97
  98   //Initialize parameters
  99   InitParameters();
 100 }
 101
 102 //________________________________________________________________________
 103 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
 104 {
 105
 106   // Create histograms to be saved in output file and
 107   // store them in fOutputContainer
 108   TList * outputContainer = new TList() ;
 109   outputContainer->SetName("CorrelationHistos") ;
 110
 111   Int_t nptbins  = GetHistoPtBins();
 112   Int_t nphibins = GetHistoPhiBins();
 113   Int_t netabins = GetHistoEtaBins();
 114   Float_t ptmax  = GetHistoPtMax();
 115   Float_t phimax = GetHistoPhiMax();
 116   Float_t etamax = GetHistoEtaMax();
 117   Float_t ptmin  = GetHistoPtMin();
 118   Float_t phimin = GetHistoPhiMin();
 119   Float_t etamin = GetHistoEtaMin();
 120
 121
 122 //  fhNclustersNtracks  = new TH2F ("Multiplicity","Neutral cluster and charged track multiplicity",1000, 0., 1000.,1000, 0., 1000.);
 123 //  fhNclustersNtracks->SetYTitle("# of tracks");
 124 //  fhNclustersNtracks->SetXTitle("# of clusters");
 125
 126   fhPtLeading  = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax);
 127   fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
 128
 129   fhPhiLeading  = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 130   fhPhiLeading->SetYTitle("#phi (rad)");
 131
 132   fhEtaLeading  = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 133   fhEtaLeading->SetYTitle("#eta ");
 134  // outputContainer->Add(fhNclustersNtracks);
 135   outputContainer->Add(fhPtLeading);
 136   outputContainer->Add(fhPhiLeading);
 137   outputContainer->Add(fhEtaLeading);
 138
 139   //Correlation with charged hadrons
 140   if(GetReader()->IsCTSSwitchedOn()) {
 141     fhDeltaPhiDeltaEtaCharged  = new TH2F
 142     ("DeltaPhiDeltaEtaCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
 143      140,-2.,5.,200,-2,2);
 144     fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
 145     fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
 146
 147     fhPhiCharged  = new TH2F
 148     ("PhiCharged","#phi_{h^{#pm}}  vs p_{T #pm}",
 149      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 150     fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
 151     fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
 152
 153     fhEtaCharged  = new TH2F
 154     ("EtaCharged","#eta_{h^{#pm}}  vs p_{T #pm}",
 155      nptbins,ptmin,ptmax,netabins,etamin,etamax);
 156     fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
 157     fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
 158
 159     fhDeltaPhiCharged  = new TH2F
 160     ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
 161      nptbins,ptmin,ptmax,140,-2.,5.);
 162     fhDeltaPhiCharged->SetYTitle("#Delta #phi");
 163     fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
 164
 165     fhDeltaPhiChargedPt  = new TH2F
 166     ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
 167      nptbins,ptmin,ptmax,140,-2.,5.);
 168     fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
 169     fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 170
 171     fhDeltaPhiUeChargedPt  = new TH2F
 172     ("DeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
 173      nptbins,ptmin,ptmax,140,-2.,5.);
 174     fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
 175     fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 176
 177     fhDeltaEtaCharged  = new TH2F
 178     ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
 179      nptbins,ptmin,ptmax,200,-2,2);
 180     fhDeltaEtaCharged->SetYTitle("#Delta #eta");
 181     fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
 182
 183     fhPtImbalanceCharged  =
 184     new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
 185              nptbins,ptmin,ptmax,200,0.,2.);
 186     fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
 187     fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
 188
 189     fhPtImbalanceUeCharged  =
 190     new TH2F("CorrelationUeCharged","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger}",
 191              nptbins,ptmin,ptmax,200,0.,2.);
 192     fhPtImbalanceUeCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
 193     fhPtImbalanceUeCharged->SetXTitle("p_{T trigger}");
 194
 195     fhPtImbalancePosCharged  =
 196     new TH2F("CorrelationPositiveCharged","z_{trigger h^{+}} = p_{T h^{+}} / p_{T trigger}",
 197              nptbins,ptmin,ptmax,200,0.,2.);
 198     fhPtImbalancePosCharged->SetYTitle("z_{trigger h^{+}}");
 199     fhPtImbalancePosCharged->SetXTitle("p_{T trigger}");
 200
 201     fhPtImbalanceNegCharged  =
 202     new TH2F("CorrelationNegativeCharged","z_{trigger h^{-}} = p_{T h^{-}} / p_{T trigger}",
 203              nptbins,ptmin,ptmax,200,0.,2.);
 204     fhPtImbalanceNegCharged->SetYTitle("z_{trigger h^{-}}");
 205     fhPtImbalanceNegCharged->SetXTitle("p_{T trigger}");
 206
 207     fhPtHbpCharged  =
 208     new TH2F("HbpCharged","#xi = ln(1/x_{E}) with charged hadrons",
 209              nptbins,ptmin,ptmax,200,0.,10.);
 210     fhPtHbpCharged->SetYTitle("ln(1/x_{E})");
 211     fhPtHbpCharged->SetXTitle("p_{T trigger}");
 212
 213     fhPtHbpUeCharged  =
 214     new TH2F("HbpUeCharged","#xi = ln(1/x_{E}) with charged hadrons",
 215              nptbins,ptmin,ptmax,200,0.,10.);
 216     fhPtHbpUeCharged->SetYTitle("ln(1/x_{E})");
 217     fhPtHbpUeCharged->SetXTitle("p_{T trigger}");
 218
 219     fhPoutPtTrigPtAssoc  =
 220     new TH3D("PoutPtTrigPtAssoc","Pout with charged hadrons",
 221              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
 222     fhPoutPtTrigPtAssoc->SetZTitle("p_{out} (GeV/c)");
 223     fhPoutPtTrigPtAssoc->SetYTitle("p_{T associated} (GeV/c)");
 224     fhPoutPtTrigPtAssoc->SetXTitle("p_{T trigger} (GeV/c)");
 225
 226     fhUePoutPtTrigPtAssoc  =
 227     new TH3D("UePoutPtTrigPtAssoc"," UE Pout with charged hadrons",
 228              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
 229     fhUePoutPtTrigPtAssoc->SetZTitle("p_{out} (GeV/c)");
 230     fhUePoutPtTrigPtAssoc->SetYTitle("p_{T associated} (GeV/c)");
 231     fhUePoutPtTrigPtAssoc->SetXTitle("p_{T trigger} (GeV/c)");
 232
 233     fhPtTrigCharged  =
 234     new TH2F("PtTrigCharged","trgger and charged tracks pt distribution",
 235              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 236     fhPtTrigCharged->SetYTitle("p_{T asso} (GeV/c)");
 237     fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)");
 238
 239     outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
 240     outputContainer->Add(fhPhiCharged) ;
 241     outputContainer->Add(fhEtaCharged) ;
 242     outputContainer->Add(fhDeltaPhiCharged) ;
 243     outputContainer->Add(fhDeltaEtaCharged) ;
 244     outputContainer->Add(fhDeltaPhiChargedPt) ;
 245     outputContainer->Add(fhDeltaPhiUeChargedPt) ;
 246     outputContainer->Add(fhPtImbalanceCharged) ;
 247     outputContainer->Add(fhPtImbalancePosCharged) ;
 248     outputContainer->Add(fhPtImbalanceNegCharged) ;
 249     outputContainer->Add(fhPtImbalanceUeCharged) ;
 250     outputContainer->Add(fhPtHbpCharged) ;
 251     outputContainer->Add(fhPtHbpUeCharged) ;
 252     outputContainer->Add(fhPoutPtTrigPtAssoc) ;
 253     outputContainer->Add(fhUePoutPtTrigPtAssoc) ;
 254     outputContainer->Add(fhPtTrigCharged) ;
 255
 256     if(DoEventSelect()){
 257       Int_t nMultiBins = GetMultiBin();
 258       fhTrigDeltaPhiDeltaEtaCharged = new TH3D*[nMultiBins] ;
 259       fhTrigCorr = new TH2F*[nMultiBins];
 260       fhTrigUeCorr = new TH2F*[nMultiBins];
 261       for(Int_t im=0; im<nMultiBins; im++){
 262         fhTrigDeltaPhiDeltaEtaCharged[im]  = new TH3D
 263         (Form("fhTrigDeltaPhiDeltaEtaCharged_%d",im),Form("fhTrigDeltaPhiDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, 140,-2.,5., 200,-2,2);
 264         fhTrigDeltaPhiDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
 265         fhTrigDeltaPhiDeltaEtaCharged[im]->SetYTitle("#Delta #phi");
 266         fhTrigDeltaPhiDeltaEtaCharged[im]->SetZTitle("#Delta #eta");
 267         fhTrigCorr[im]  = new TH2F
 268         (Form("fhTrigPtCorr_%d",im),Form("fhTrigPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 269         fhTrigCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
 270         fhTrigCorr[im]->SetXTitle("p_{T trigger}");
 271         fhTrigUeCorr[im]  = new TH2F
 272         (Form("fhTrigPtUeCorr_%d",im),Form("fhTrigPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 273         fhTrigUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
 274         fhTrigUeCorr[im]->SetXTitle("p_{T trigger}");
 275
 276         outputContainer->Add(fhTrigDeltaPhiDeltaEtaCharged[im]) ;
 277         outputContainer->Add(fhTrigCorr[im]);
 278         outputContainer->Add(fhTrigUeCorr[im]);
 279
 280         }
 281     }
 282
 283     if(fPi0Trigger){
 284       fhPtPi0DecayRatio  = new TH2F
 285       ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay",
 286        nptbins,ptmin,ptmax, 100,0.,2.);
 287       fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");
 288       fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");
 289
 290       fhDeltaPhiDecayCharged  = new TH2F
 291       ("DeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",
 292        nptbins,ptmin,ptmax,140,-2.,5.);
 293       fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");
 294       fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");
 295
 296       fhPtImbalanceDecayCharged  =
 297       new TH2F("CorrelationDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",
 298                nptbins,ptmin,ptmax,200,0.,2.);
 299       fhPtImbalanceDecayCharged->SetYTitle("z_{decay h^{#pm}}");
 300       fhPtImbalanceDecayCharged->SetXTitle("p_{T decay}");
 301
 302       outputContainer->Add(fhPtPi0DecayRatio) ;
 303       outputContainer->Add(fhDeltaPhiDecayCharged) ;
 304       outputContainer->Add(fhPtImbalanceDecayCharged) ;
 305     }
 306
 307
 308     if(fMakeSeveralUE){
 309       fhDeltaPhiUeLeftCharged  = new TH2F
 310       ("DeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
 311        nptbins,ptmin,ptmax,140,-2.,5.);
 312       fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
 313       fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 314       outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
 315
 316       fhDeltaPhiUeRightCharged  = new TH2F
 317       ("DeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
 318        nptbins,ptmin,ptmax,140,-2.,5.);
 319       fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
 320       fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 321       outputContainer->Add(fhDeltaPhiUeRightCharged) ;
 322
 323       fhPtImbalanceUeLeftCharged  =
 324       new TH2F("CorrelationUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",
 325                nptbins,ptmin,ptmax,200,0.,2.);
 326       fhPtImbalanceUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
 327       fhPtImbalanceUeLeftCharged->SetXTitle("p_{T trigger}");
 328       outputContainer->Add(fhPtImbalanceUeLeftCharged) ;
 329
 330       fhPtImbalanceUeRightCharged  =
 331       new TH2F("CorrelationUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",
 332                nptbins,ptmin,ptmax,200,0.,2.);
 333       fhPtImbalanceUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
 334       fhPtImbalanceUeRightCharged->SetXTitle("p_{T trigger}");
 335       outputContainer->Add(fhPtImbalanceUeRightCharged) ;
 336
 337       fhPtHbpUeLeftCharged  =
 338       new TH2F("HbpUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",
 339                nptbins,ptmin,ptmax,200,0.,10.);
 340       fhPtHbpUeLeftCharged->SetYTitle("ln(1/x_{E})");
 341       fhPtHbpUeLeftCharged->SetXTitle("p_{T trigger}");
 342       outputContainer->Add(fhPtHbpUeLeftCharged) ;
 343
 344       fhPtHbpUeRightCharged  =
 345       new TH2F("HbpUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",
 346                nptbins,ptmin,ptmax,200,0.,10.);
 347       fhPtHbpUeRightCharged->SetYTitle("ln(1/x_{E})");
 348       fhPtHbpUeRightCharged->SetXTitle("p_{T trigger}");
 349       outputContainer->Add(fhPtHbpUeRightCharged) ;
 350
 351     }
 352   }  //Correlation with charged hadrons
 353
 354   //Correlation with neutral hadrons
 355   if(fNeutralCorr){
 356
 357     fhDeltaPhiDeltaEtaNeutral  = new TH2F
 358     ("DeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
 359      140,-2.,5.,200,-2,2);
 360     fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");
 361     fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
 362
 363     fhPhiNeutral  = new TH2F
 364     ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T #pi^{0}}",
 365      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 366     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
 367     fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
 368
 369     fhEtaNeutral  = new TH2F
 370     ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T #pi^{0}}",
 371      nptbins,ptmin,ptmax,netabins,etamin,etamax);
 372     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
 373     fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
 374
 375     fhDeltaPhiNeutral  = new TH2F
 376     ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
 377      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 378     fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
 379     fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
 380
 381     fhDeltaPhiNeutralPt  = new TH2F
 382     ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
 383      nptbins,ptmin,ptmax,140,-2.,5.);
 384     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
 385     fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
 386
 387     fhDeltaPhiUeNeutralPt  = new TH2F
 388     ("DeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
 389      nptbins,ptmin,ptmax,140,-2.,5.);
 390     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
 391     fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
 392
 393     fhDeltaEtaNeutral  = new TH2F
 394     ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
 395      nptbins,ptmin,ptmax,200,-2,2);
 396     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
 397     fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
 398
 399     fhPtImbalanceNeutral  =
 400     new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
 401              nptbins,ptmin,ptmax,200,0.,2.);
 402     fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
 403     fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
 404
 405     fhPtImbalanceUeNeutral  =
 406     new TH2F("CorrelationUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
 407              nptbins,ptmin,ptmax,200,0.,2.);
 408     fhPtImbalanceUeNeutral->SetYTitle("z_{trigger #pi^{0}}");
 409     fhPtImbalanceUeNeutral->SetXTitle("p_{T trigger}");
 410
 411     fhPtHbpNeutral  =
 412     new TH2F("HbpNeutral","#xi = ln(1/x_{E}) with neutral particles",
 413              nptbins,ptmin,ptmax,200,0.,10.);
 414     fhPtHbpNeutral->SetYTitle("ln(1/x_{E})");
 415     fhPtHbpNeutral->SetXTitle("p_{T trigger}");
 416
 417     fhPtHbpUeNeutral  =
 418     new TH2F("HbpUeNeutral","#xi = ln(1/x_{E}) with neutral particles",
 419              nptbins,ptmin,ptmax,200,0.,10.);
 420     fhPtHbpUeNeutral->SetYTitle("ln(1/x_{E})");
 421     fhPtHbpUeNeutral->SetXTitle("p_{T trigger}");
 422
 423
 424     outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
 425     outputContainer->Add(fhPhiNeutral) ;
 426     outputContainer->Add(fhEtaNeutral) ;
 427     outputContainer->Add(fhDeltaPhiNeutral) ;
 428     outputContainer->Add(fhDeltaPhiNeutralPt) ;
 429     outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
 430     outputContainer->Add(fhDeltaEtaNeutral) ;
 431     outputContainer->Add(fhPtImbalanceNeutral) ;
 432     outputContainer->Add(fhPtImbalanceUeNeutral) ;
 433     outputContainer->Add(fhPtHbpNeutral) ;
 434     outputContainer->Add(fhPtHbpUeNeutral) ;
 435
 436     if(fPi0Trigger){
 437       fhDeltaPhiDecayNeutral  = new TH2F
 438       ("DeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",
 439        nptbins,ptmin,ptmax,140,-2.,5.);
 440       fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");
 441       fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");
 442
 443       fhPtImbalanceDecayNeutral  =
 444       new TH2F("CorrelationDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",
 445                nptbins,ptmin,ptmax,200,0.,2.);
 446       fhPtImbalanceDecayNeutral->SetYTitle("z_{decay h^{0}}");
 447       fhPtImbalanceDecayNeutral->SetXTitle("p_{T decay}");
 448
 449       outputContainer->Add(fhDeltaPhiDecayNeutral) ;
 450       outputContainer->Add(fhPtImbalanceDecayNeutral) ;
 451     }
 452
 453
 454     if(fMakeSeveralUE){
 455       fhDeltaPhiUeLeftNeutral  = new TH2F
 456       ("DeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",
 457        nptbins,ptmin,ptmax,140,-2.,5.);
 458       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
 459       fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
 460       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
 461
 462       fhDeltaPhiUeRightNeutral  = new TH2F
 463       ("DeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",
 464        nptbins,ptmin,ptmax,140,-2.,5.);
 465       fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
 466       fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
 467       outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
 468
 469       fhPtImbalanceUeLeftNeutral  =
 470       new TH2F("CorrelationUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
 471                nptbins,ptmin,ptmax,140,0.,2.);
 472       fhPtImbalanceUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
 473       fhPtImbalanceUeLeftNeutral->SetXTitle("p_{T trigger}");
 474       outputContainer->Add(fhPtImbalanceUeLeftNeutral) ;
 475
 476       fhPtImbalanceUeRightNeutral  =
 477       new TH2F("CorrelationUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
 478                nptbins,ptmin,ptmax,200,0.,2.);
 479       fhPtImbalanceUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
 480       fhPtImbalanceUeRightNeutral->SetXTitle("p_{T trigger}");
 481       outputContainer->Add(fhPtImbalanceUeRightNeutral) ;
 482
 483       fhPtHbpUeLeftNeutral  =
 484       new TH2F("HbpUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",
 485                nptbins,ptmin,ptmax,200,0.,10.);
 486       fhPtHbpUeLeftNeutral->SetYTitle("ln(1/x_{E})");
 487       fhPtHbpUeLeftNeutral->SetXTitle("p_{T trigger}");
 488       outputContainer->Add(fhPtHbpUeLeftNeutral) ;
 489
 490       fhPtHbpUeRightNeutral  =
 491       new TH2F("HbpUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",
 492                nptbins,ptmin,ptmax,200,0.,10.);
 493       fhPtHbpUeRightNeutral->SetYTitle("ln(1/x_{E})");
 494       fhPtHbpUeRightNeutral->SetXTitle("p_{T trigger}");
 495       outputContainer->Add(fhPtHbpUeRightNeutral) ;
 496
 497     }
 498
 499     //Keep neutral meson selection histograms if requiered
 500     //Setting done in AliNeutralMesonSelection
 501     if(GetNeutralMesonSelection()){
 502       TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
 503       if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
 504         for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
 505       delete nmsHistos;
 506     }
 507
 508   }//Correlation with neutral hadrons
 509
 510   return outputContainer;
 511
 512 }
 513
 514 //____________________________________________________________________________
 515 void AliAnaParticleHadronCorrelation::InitParameters()
 516 {
 517
 518   //Initialize the parameters of the analysis.
 519   SetInputAODName("PWG4Particle");
 520   SetAODObjArrayName("Hadrons");
 521   AddToHistogramsName("AnaHadronCorr_");
 522
 523   SetPtCutRange(0.,300);
 524   fDeltaPhiMinCut = 1.5 ;
 525   fDeltaPhiMaxCut = 4.5 ;
 526   fSelectIsolated = kFALSE;
 527   fMakeSeveralUE = kFALSE;
 528   fUeDeltaPhiMinCut = 1. ;
 529   fUeDeltaPhiMaxCut = 1.5 ;
 530   fNeutralCorr = kFALSE ;
 531   fPi0Trigger = kFALSE ;
 532
 533 }
 534
 535 //__________________________________________________________________
 536 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
 537 {
 538
 539   //Print some relevant parameters set for the analysis
 540   if(! opt)
 541     return;
 542
 543   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 544   AliAnaPartCorrBaseClass::Print(" ");
 545
 546   printf("Phi trigger particle-Hadron      <     %3.2f\n", fDeltaPhiMaxCut) ;
 547   printf("Phi trigger particle-Hadron      >     %3.2f\n", fDeltaPhiMinCut) ;
 548   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
 549   printf("Phi trigger particle-UeHadron    <    %3.2f\n", fUeDeltaPhiMaxCut) ;
 550   printf("Phi trigger particle-UeHadron    >    %3.2f\n", fUeDeltaPhiMinCut) ;
 551   printf("Several UE?  %d\n", fMakeSeveralUE) ;
 552   printf("Name of AOD Pi0 Branch %s \n",fPi0AODBranchName.Data());
 553   printf("Do Decay-hadron correlation ?  %d\n", fPi0Trigger) ;
 554
 555
 556 }
 557
 558 //__________________________________________________________________
 559 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
 560 {
 561   //Save parameters used for analysis
 562   TString parList ; //this will be list of parameters used for this analysis.
 563   const Int_t buffersize = 255;
 564   char onePar[buffersize] ;
 565
 566   snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;
 567   parList+=onePar ;
 568   snprintf(onePar,buffersize,"Phi trigger particle-Hadron      <     %3.2f ", fDeltaPhiMaxCut) ;
 569   parList+=onePar ;
 570   snprintf(onePar,buffersize,"Phi trigger particle-Hadron      >     %3.2f ", fDeltaPhiMinCut) ;
 571   parList+=onePar ;
 572   snprintf(onePar,buffersize,"Isolated Trigger?  %d\n", fSelectIsolated) ;
 573   parList+=onePar ;
 574   snprintf(onePar,buffersize,"Phi trigger particle-UeHadron    <    %3.2f ", fUeDeltaPhiMaxCut) ;
 575   parList+=onePar ;
 576   snprintf(onePar,buffersize,"Phi trigger particle-UeHadron    >    %3.2f ", fUeDeltaPhiMinCut) ;
 577   parList+=onePar ;
 578   snprintf(onePar,buffersize,"Several UE?  %d\n", fMakeSeveralUE) ;
 579   parList+=onePar ;
 580   snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ",fPi0AODBranchName.Data());
 581   parList+=onePar ;
 582   snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  %d", fPi0Trigger) ;
 583   parList+=onePar ;
 584
 585   //Get parameters set in base class.
 586   parList += GetBaseParametersList() ;
 587
 588   //Get parameters set in PID class.
 589   //parList += GetCaloPID()->GetPIDParametersList() ;
 590
 591   //Get parameters set in FiducialCut class (not available yet)
 592   //parlist += GetFidCut()->GetFidCutParametersList()
 593
 594   return new TObjString(parList) ;
 595
 596 }
 597
 598
 599 //____________________________________________________________________________
 600 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
 601 {
 602   //Particle-Hadron Correlation Analysis, fill AODs
 603
 604   if(!GetInputAODBranch()){
 605     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
 606     abort();
 607   }
 608
 609   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
 610     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());
 611     abort();
 612   }
 613
 614   if(GetDebug() > 1){
 615     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");
 616     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
 617     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
 618     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
 619     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
 620   }
 621
 622   //Loop on stored AOD particles, trigger
 623   Double_t ptTrig    = 0.;
 624   Int_t    trigIndex = -1;
 625   Int_t naod = GetInputAODBranch()->GetEntriesFast();
 626   //fhNclustersNtracks->Fill(naod, GetAODCTS()->GetEntriesFast());
 627   for(Int_t iaod = 0; iaod < naod ; iaod++){
 628     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 629     //find the leading particles with highest momentum
 630     if (particle->Pt()>ptTrig) {
 631       ptTrig = particle->Pt() ;
 632       trigIndex = iaod ;
 633     }
 634   }//Aod branch loop
 635
 636   //Do correlation with leading particle
 637   if(trigIndex!=-1){
 638
 639     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));
 640     //Make correlation with charged hadrons
 641     if(GetReader()->IsCTSSwitchedOn() )
 642       MakeChargedCorrelation(particle, GetAODCTS(),kFALSE);
 643
 644     TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
 645     if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
 646       MakeNeutralCorrelation(particle, pi0list,kFALSE);
 647
 648   }//Correlate leading
 649
 650   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
 651
 652 }
 653
 654 //____________________________________________________________________________
 655 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
 656 {
 657   //Particle-Hadron Correlation Analysis, fill histograms
 658
 659   if(!GetInputAODBranch()){
 660     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
 661     abort();
 662   }
 663
 664   if(GetDebug() > 1){
 665     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
 666     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
 667   }
 668
 669   //Get the vertex and check it is not too large in z
 670   Double_t v[3] = {0,0,0}; //vertex ;
 671   GetReader()->GetVertex(v);
 672   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
 673
 674   //Loop on stored AOD particles, find leading
 675   Int_t naod = GetInputAODBranch()->GetEntriesFast();
 676  // if(naod!=0)fhVertex->Fill(v[0],v[1],v[2]);
 677   Double_t ptTrig = 0.;
 678   Int_t trigIndex = -1 ;
 679   for(Int_t iaod = 0; iaod < naod ; iaod++){     //loop on input trigger AOD file
 680     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 681     //vertex cut in case of mixing
 682     if (GetMixedEvent()) {
 683       Int_t evt=-1;
 684       Int_t id =-1;
 685       if     (particle->GetCaloLabel(0)  >= 0 ){
 686         id=particle->GetCaloLabel(0);
 687         if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
 688       }
 689       else if(particle->GetTrackLabel(0) >= 0 ){
 690         id=particle->GetTrackLabel(0);
 691         if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
 692       }
 693       else continue;
 694
 695       if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
 696         return ;
 697     }
 698
 699     //check if the particle is isolated or if we want to take the isolation into account
 700     if(OnlyIsolated() && !particle->IsIsolated()) continue;
 701     //check if inside the vertex cut
 702     //find the leading particles with highest momentum
 703     if (particle->Pt()>ptTrig) {
 704       ptTrig = particle->Pt() ;
 705                   trigIndex = iaod ;
 706     }
 707   }//finish searching for leading trigger particle
 708   if(trigIndex!=-1){ //using trigger partilce to do correlations
 709     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));
 710
 711     if (GetMixedEvent()) {
 712       Int_t evt=-1;
 713       Int_t id = 0;
 714       if     (particle->GetCaloLabel(0)  >= 0 ){
 715         id=particle->GetCaloLabel(0);
 716         if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
 717       }
 718       else if(particle->GetTrackLabel(0) >= 0 ){
 719         id=particle->GetTrackLabel(0);
 720         if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
 721       }
 722       else return;
 723
 724       if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
 725         return ;
 726     }
 727
 728     //Fill leading particle histogram
 729     fhPtLeading->Fill(particle->Pt());
 730     Float_t phi = particle->Phi();
 731     if(phi<0)phi+=TMath::TwoPi();
 732     fhPhiLeading->Fill(particle->Pt(), phi);
 733     fhEtaLeading->Fill(particle->Pt(), particle->Eta());
 734
 735     //Make correlation with charged hadrons
 736     if(GetReader()->IsCTSSwitchedOn() )
 737       MakeChargedCorrelation(particle, GetAODCTS(),kTRUE);
 738
 739     TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
 740     if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
 741       MakeNeutralCorrelation(particle, pi0list,kTRUE);
 742
 743   }//Aod branch loop
 744
 745   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
 746
 747 }
 748
 749 //____________________________________________________________________________
 750 void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle, TObjArray* const pl, const Bool_t bFillHisto)
 751 {
 752   // Charged Hadron Correlation Analysis
 753   if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
 754
 755   Int_t evtIndex11 = 0 ;
 756   Int_t evtIndex12 = 0 ;
 757   Int_t indexPhoton1 = -1 ;
 758   Int_t indexPhoton2 = -1 ;
 759   Double_t v[3] = {0,0,0}; //vertex ;
 760   GetReader()->GetVertex(v);
 761   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
 762
 763   Int_t nTracks = GetAODCTS()->GetEntriesFast() ;
 764
 765   if (GetMixedEvent()) {
 766     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
 767     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
 768   }
 769
 770   Double_t ptTrig  = aodParticle->Pt();
 771   Double_t pxTrig  = aodParticle->Px();
 772   Double_t pyTrig  = aodParticle->Py();
 773
 774   Double_t phiTrig = aodParticle->Phi();
 775   Double_t etaTrig = aodParticle->Eta();
 776
 777   Double_t pt   = -100.;
 778   Double_t px   = -100.;
 779   Double_t py   = -100.;
 780   Double_t rat  = -100.;
 781   Double_t xE   = -100.;
 782   Double_t cosi = -100.;
 783   Double_t phi  = -100. ;
 784   Double_t eta  = -100. ;
 785   TVector3 p3;
 786
 787   TObjArray * reftracks = 0x0;
 788   Int_t nrefs = 0;
 789
 790   Double_t ptDecay1  = 0. ;
 791   Double_t pxDecay1  = 0. ;
 792   Double_t pyDecay1  = 0. ;
 793   Double_t phiDecay1 = 0. ;
 794   Double_t ptDecay2  = 0. ;
 795   Double_t pxDecay2  = 0. ;
 796   Double_t pyDecay2  = 0. ;
 797   Double_t phiDecay2 = 0. ;
 798
 799   Double_t ratDecay1  = -100.;
 800   Double_t ratDecay2  = -100.;
 801   Float_t deltaphi    = -100. ;
 802   Float_t deltaphiDecay1 = -100. ;
 803   Float_t deltaphiDecay2 = -100. ;
 804   TObjArray * clusters   = 0x0 ;
 805   TLorentzVector photonMom ;
 806
 807   if(fPi0Trigger){
 808     indexPhoton1 = aodParticle->GetCaloLabel (0);
 809     indexPhoton2 = aodParticle->GetCaloLabel (1);
 810     if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
 811
 812     if(indexPhoton1!=-1 && indexPhoton2!=-1){
 813       if(aodParticle->GetDetector()=="EMCAL") clusters = GetAODEMCAL() ;
 814       else  clusters = GetAODPHOS() ;
 815       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
 816         AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));
 817         photon->GetMomentum(photonMom,GetVertex(0)) ;
 818         if(photon->GetID()==indexPhoton1) {
 819           ptDecay1  = photonMom.Pt();
 820           pxDecay1  = photonMom.Px();
 821           pyDecay1  = photonMom.Py();
 822           phiDecay1 = photonMom.Phi();
 823           if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig);
 824         }
 825         if(photon->GetID()==indexPhoton2) {
 826           ptDecay2  = photonMom.Pt();
 827           pxDecay2  = photonMom.Px();
 828           pyDecay2  = photonMom.Py();
 829           phiDecay2 = photonMom.Phi();
 830           if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay2/ptTrig);
 831         }
 832         if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
 833       } //cluster loop
 834     } //index of decay photons found
 835
 836   } //make decay-hadron correlation
 837
 838   //Track loop, select tracks with good pt, phi and fill AODs or histograms
 839   //Int_t currentIndex = -1 ;
 840   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){
 841     AliAODTrack * track = (AliAODTrack *) (pl->At(ipr)) ;
 842
 843     //check if inside the vertex cut
 844     //printf("charge = %d\n", track->Charge());
 845     Int_t evtIndex2 = 0 ;
 846     if (GetMixedEvent()) {
 847       evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
 848       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2) // photon and track from different events
 849         continue ;
 850       //vertex cut
 851       if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
 852          return ;
 853       //      if(currentIndex == evtIndex2) // tracks from different event
 854       //        continue ;
 855       //      currentIndex = evtIndex2 ;
 856     }
 857     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 858     p3.SetXYZ(mom[0],mom[1],mom[2]);
 859     pt   = p3.Pt();
 860     px   = p3.Px();
 861     py   = p3.Py();
 862     eta  = p3.Eta();
 863     phi  = p3.Phi() ;
 864     if(phi < 0) phi+=TMath::TwoPi();
 865
 866     //Select only hadrons in pt range
 867     if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
 868     //remove trigger itself for correlation when use charged triggers
 869     if(track->GetID()==aodParticle->GetTrackLabel(0) && pt==ptTrig && phi==phiTrig && eta==etaTrig)
 870       continue ;
 871     if(IsFiducialCutOn()){
 872       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
 873       if(! in ) continue ;
 874     }
 875     //jumped out this event if near side associated partile pt larger than trigger
 876     if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2())  break ;
 877     rat   = pt/ptTrig ;
 878     xE    = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 879     if(xE <0.)xE =-xE;
 880     cosi = TMath::Log(1/xE);
 881     // printf("rat = %f, xE = %f, cosi =%f \n", rat, xE, cosi);
 882     // printf("phi = %f \n", phi);
 883
 884      if(fPi0Trigger){
 885       if(indexPhoton1!=-1 && indexPhoton2!=-1){
 886         if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
 887         if(ptDecay2) ratDecay2 = pt/ptDecay2 ;
 888         deltaphiDecay1 = phiDecay1-phi;
 889         deltaphiDecay2 = phiDecay2-phi;
 890         if(deltaphiDecay1< -TMath::PiOver2()) deltaphiDecay1+=TMath::TwoPi();
 891         if(deltaphiDecay1>3*TMath::PiOver2()) deltaphiDecay1-=TMath::TwoPi();
 892         if(deltaphiDecay2< -TMath::PiOver2()) deltaphiDecay2+=TMath::TwoPi();
 893         if(deltaphiDecay2>3*TMath::PiOver2()) deltaphiDecay2-=TMath::TwoPi();
 894       }
 895     } //do decay-hadron correlation
 896
 897     //Selection within angular range
 898     deltaphi = phiTrig-phi;
 899     if(deltaphi< -TMath::PiOver2()) deltaphi+=TMath::TwoPi();
 900     if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();
 901
 902     Double_t pout = pt*TMath::Sin(deltaphi) ;
 903
 904     if(GetDebug() > 2)
 905       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",
 906              pt,phi, phiTrig,fDeltaPhiMinCut, deltaphi, fDeltaPhiMaxCut, GetMinPt());
 907
 908     if(bFillHisto){
 909       // Fill Histograms
 910       fhEtaCharged->Fill(pt,eta);
 911       fhPhiCharged->Fill(pt,phi);
 912       fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
 913       fhDeltaPhiCharged->Fill(ptTrig, deltaphi);
 914       fhDeltaPhiDeltaEtaCharged->Fill(deltaphi,aodParticle->Eta()-eta);
 915
 916       if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
 917       //fill different multiplicity histogram
 918       if(DoEventSelect()){
 919         for(Int_t im=0; im<GetMultiBin(); im++){
 920           if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))fhTrigDeltaPhiDeltaEtaCharged[im]->Fill(ptTrig,deltaphi,aodParticle->Eta()-eta);
 921         }
 922       }
 923       //delta phi cut for correlation
 924       if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {
 925         fhDeltaPhiChargedPt->Fill(pt,deltaphi);
 926         fhPtImbalanceCharged->Fill(ptTrig,xE);
 927         fhPtHbpCharged->Fill(ptTrig,cosi);
 928         fhPoutPtTrigPtAssoc->Fill(ptTrig, pt, pout) ;
 929         fhPtTrigCharged->Fill(ptTrig, pt) ;
 930         if(track->Charge()>0) fhPtImbalancePosCharged->Fill(ptTrig,xE) ;
 931         else fhPtImbalanceNegCharged->Fill(ptTrig,xE) ;
 932         //fill different multiplicity histogram
 933         if(DoEventSelect()){
 934           for(Int_t im=0; im<GetMultiBin(); im++){
 935             if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
 936               fhTrigCorr[im]->Fill(ptTrig,xE);
 937           }
 938         } //multiplicity events selection
 939       } //delta phi cut for correlation
 940       else { //UE study
 941         fhDeltaPhiUeChargedPt->Fill(pt,deltaphi);
 942         fhUePoutPtTrigPtAssoc->Fill(ptTrig, pt, pout) ;
 943         fhPtImbalanceUeCharged->Fill(ptTrig,xE);
 944         fhPtHbpUeCharged->Fill(ptTrig,cosi);
 945         if(DoEventSelect()){
 946           for(Int_t im=0; im<GetMultiBin(); im++){
 947             if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
 948               fhTrigUeCorr[im]->Fill(ptTrig,xE);
 949           }
 950         } //multiplicity events selection
 951
 952       } //UE study
 953
 954       if(fPi0Trigger){
 955         if(indexPhoton1!=-1 && indexPhoton2!=-1){
 956           fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaphiDecay1);
 957           fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaphiDecay2);
 958           if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);
 959           if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )
 960             fhPtImbalanceDecayCharged->Fill(ptDecay1,ratDecay1);
 961           if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )
 962             fhPtImbalanceDecayCharged->Fill(ptDecay2,ratDecay2);
 963           if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
 964         } //index of decay photons found
 965       } //make decay-hadron correlation
 966
 967       //several UE calculation
 968       if(fMakeSeveralUE){
 969         if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){
 970           fhDeltaPhiUeLeftCharged->Fill(pt,deltaphi);
 971           fhPtImbalanceUeLeftCharged->Fill(ptTrig,rat);
 972           fhPtHbpUeLeftCharged->Fill(ptTrig,cosi);
 973         }
 974         if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){
 975           fhDeltaPhiUeRightCharged->Fill(pt,deltaphi);
 976           fhPtImbalanceUeRightCharged->Fill(ptTrig,rat);
 977           fhPtHbpUeRightCharged->Fill(ptTrig,cosi);
 978
 979         }
 980       } //several UE calculation
 981
 982     } //Fill histogram
 983     else{
 984       nrefs++;
 985       if(nrefs==1){
 986         reftracks = new TObjArray(0);
 987         reftracks->SetName(GetAODObjArrayName()+"Tracks");
 988         reftracks->SetOwner(kFALSE);
 989       }
 990       reftracks->Add(track);
 991     }//aod particle loop
 992   }// track loop
 993
 994   //Fill AOD with reference tracks, if not filling histograms
 995   if(!bFillHisto && reftracks) {
 996     aodParticle->AddObjArray(reftracks);
 997   }
 998
 999   //delete reftracks;
1000
1001 }
1002
1003
1004
1005 //____________________________________________________________________________
1006 //void  AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD(AliAODPWG4ParticleCorrelation* const aodParticle,TObjArray* const pl, TString detector)
1007 //{
1008 //    // Neutral Pion Correlation Analysis, find pi0, put them in new output aod, if correlation cuts passed
1009 //  if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Make trigger particle - neutral hadron correlation \n");
1010 //
1011 //  if(!NewOutputAOD()){
1012 //    printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Output aod not created, set AOD class name and branch name in the configuration file, STOP! \n");
1013 //    abort();
1014 //  }
1015 //
1016 //  Double_t phiTrig = aodParticle->Phi();
1017 //  Int_t       tag = 0;
1018 //  TLorentzVector gammai;
1019 //  TLorentzVector gammaj;
1020 //
1021 //  //Get vertex for photon momentum calculation
1022 //
1023 //  if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
1024 //  {
1025 //    for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) {
1026 //      if (!GetMixedEvent())
1027 //        GetReader()->GetVertex(GetVertex(iev));
1028 //      else
1029 //        GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev));
1030 //    }
1031 //  }
1032 //  Double_t vertex2[] = {0.0,0.0,0.0} ; //vertex of second input aod
1033 //  if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
1034 //  {
1035 //     if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
1036 //  }
1037 //
1038 //    //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
1039 //    //Int_t iEvent= GetReader()->GetEventNumber() ;
1040 //  Int_t nclus = pl->GetEntriesFast();
1041 //  for(Int_t iclus = 0;iclus < nclus ; iclus ++ ){
1042 //    AliVCluster * calo = (AliVCluster *) (pl->At(iclus)) ;
1043 //
1044 //    Int_t evtIndex1 = 0 ;
1045 //    if (GetMixedEvent()) {
1046 //      evtIndex1=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1047 //    }
1048 //
1049 //
1050 //      //Input from second AOD?
1051 //    Int_t inputi = 0;
1052 //    if     (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= iclus)
1053 //      inputi = 1 ;
1054 //    else if(aodParticle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= iclus)
1055 //      inputi = 1;
1056 //
1057 //      //Cluster selection, not charged, with photon or pi0 id and in fiducial cut
1058 //    //FIXME
1059 //    Int_t pdg=0;
1060 //    //if     (inputi == 0 && !SelectCluster(calo, GetVertex(evtIndex1),  gammai, pdg))
1061 //      continue ;
1062 //    //MEFIX
1063 //    else if(inputi == 1 && !SelectCluster(calo, vertex2, gammai, pdg))
1064 //      continue ;
1065 //
1066 //    if(GetDebug() > 2)
1067 //      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",
1068 //             detector.Data(), gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());
1069 //
1070 //      //2 gamma overlapped, found with PID
1071 //    if(pdg == AliCaloPID::kPi0){
1072 //
1073 //        //Select only hadrons in pt range
1074 //      if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt())
1075 //        continue ;
1076 //
1077 //        //Selection within angular range
1078 //      Float_t phi = gammai.Phi();
1079 //      if(phi < 0) phi+=TMath::TwoPi();
1080 //        //Float_t deltaphi = TMath::Abs(phiTrig-phi);
1081 //        //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
1082 //
1083 //      AliAODPWG4Particle pi0 = AliAODPWG4Particle(gammai);
1084 //        //pi0.SetLabel(calo->GetLabel());
1085 //      pi0.SetPdg(AliCaloPID::kPi0);
1086 //      pi0.SetDetector(detector);
1087 //
1088 //      if(IsDataMC()){
1089 //        pi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(),inputi));
1090 //        if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of candidate %d\n",pi0.GetTag());
1091 //      }//Work with stack also
1092 //       //Set the indeces of the original caloclusters
1093 //      pi0.SetCaloLabel(calo->GetID(),-1);
1094 //      AddAODParticle(pi0);
1095 //
1096 //      if(GetDebug() > 2)
1097 //        printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Correlated with selected pi0 (pid): pt %f, phi %f\n",pi0.Pt(),pi0.Phi());
1098 //
1099 //    }// pdg = 111
1100 //
1101 //      //Make invariant mass analysis
1102 //    else if(pdg == AliCaloPID::kPhoton){
1103 //        //Search the photon companion in case it comes from  a Pi0 decay
1104 //        //Apply several cuts to select the good pair;
1105 //      for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){
1106 //        AliVCluster * calo2 = (AliVCluster *) (pl->At(jclus)) ;
1107 //        Int_t evtIndex2 = 0 ;
1108 //        if (GetMixedEvent()) {
1109 //          evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1110 //        }
1111 //        if (GetMixedEvent() && (evtIndex1 == evtIndex2))
1112 //          continue ;
1113 //
1114 //          //Input from second AOD?
1115 //        Int_t inputj = 0;
1116 //        if     (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetAODEMCALNormalInputEntries() <= jclus)
1117 //          inputj = 1;
1118 //        else if(aodParticle->GetDetector() == "PHOS"  && GetReader()->GetAODPHOSNormalInputEntries()  <= jclus)
1119 //          inputj = 1;
1120 //
1121 //          //Cluster selection, not charged with photon or pi0 id and in fiducial cut
1122 //        Int_t pdgj=0;
1123 //        //FIXME
1124 //        //if     (inputj == 0 && !SelectCluster(calo2, GetVertex(evtIndex2),  gammaj, pdgj))
1125 //          continue ;
1126 //        //MEFIX
1127 //
1128 //        else if(inputj == 1 && !SelectCluster(calo2, vertex2, gammaj, pdgj))
1129 //          continue ;
1130 //        //FIXME
1131 //        //if(!SelectCluster(calo2,GetVertex(evtIndex2), gammaj, pdgj))
1132 //        //MEFIX
1133 //          continue ;
1134 //
1135 //        if(pdgj == AliCaloPID::kPhoton ){
1136 //
1137 //          if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt())
1138 //            continue ;
1139 //
1140 //            //Selection within angular range
1141 //          Float_t phi = (gammai+gammaj).Phi();
1142 //          if(phi < 0) phi+=TMath::TwoPi();
1143 //            //Float_t deltaphi = TMath::Abs(phiTrig-phi);
1144 //            //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
1145 //
1146 //            //Select good pair (aperture and invariant mass)
1147 //          if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
1148 //
1149 //            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",
1150 //                                       (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());
1151 //
1152 //            TLorentzVector pi0mom = gammai+gammaj;
1153 //            AliAODPWG4Particle pi0 = AliAODPWG4Particle(pi0mom);
1154 //              //pi0.SetLabel(calo->GetLabel());
1155 //            pi0.SetPdg(AliCaloPID::kPi0);
1156 //            pi0.SetDetector(detector);
1157 //            if(IsDataMC()){
1158 //                //Check origin of the candidates
1159 //
1160 //              Int_t label1 = calo->GetLabel();
1161 //              Int_t label2 = calo2->GetLabel();
1162 //              Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), inputi);
1163 //              Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), inputj);
1164 //
1165 //              if(GetDebug() > 0)
1166 //                printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
1167 //              if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){
1168 //
1169 //                  //Check if pi0 mother is the same
1170 //                if(GetReader()->ReadStack()){
1171 //                  TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
1172 //                  label1 = mother1->GetFirstMother();
1173 //                    //mother1 = GetMCStack()->Particle(label1);//pi0
1174 //
1175 //                  TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
1176 //                  label2 = mother2->GetFirstMother();
1177 //                    //mother2 = GetMCStack()->Particle(label2);//pi0
1178 //                }
1179 //                else if(GetReader()->ReadAODMCParticles()){
1180 //                  AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputi))->At(label1);//photon in kine tree
1181 //                  label1 = mother1->GetMother();
1182 //                    //mother1 = GetMCStack()->Particle(label1);//pi0
1183 //                  AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputj))->At(label2);//photon in kine tree
1184 //                  label2 = mother2->GetMother();
1185 //                    //mother2 = GetMCStack()->Particle(label2);//pi0
1186 //                }
1187 //
1188 //                  //printf("mother1 %d, mother2 %d\n",label1,label2);
1189 //                if(label1 == label2)
1190 //                  GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
1191 //              }
1192 //            }//Work with mc information also
1193 //            pi0.SetTag(tag);
1194 //              //Set the indeces of the original caloclusters
1195 //            pi0.SetCaloLabel(calo->GetID(), calo2->GetID());
1196 //            AddAODParticle(pi0);
1197 //
1198 //
1199 //          }//Pair selected
1200 //        }//if pair of gammas
1201 //      }//2nd loop
1202 //    }// if pdg = 22
1203 //  }//1st loop
1204 //
1205 //  if(GetDebug() > 1)
1206 //    printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - End, %d pi0's found \n",GetOutputAODBranch()->GetEntriesFast());
1207 //}
1208
1209 //____________________________________________________________________________
1210 void  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle, TObjArray* pi0list, const Bool_t bFillHisto)
1211 {
1212   // Neutral Pion Correlation Analysis
1213   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",pi0list->GetEntriesFast());
1214
1215   Int_t evtIndex11 = 0 ;
1216   Int_t evtIndex12 = 0 ;
1217   if (GetMixedEvent()) {
1218     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
1219     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
1220   }
1221
1222   Double_t pt   = -100.;
1223   Double_t px   = -100.;
1224   Double_t py   = -100.;
1225   Double_t rat = -100.;
1226   Double_t phi = -100.;
1227   Double_t eta = -100.;
1228   Double_t xE  = -100.;
1229   Double_t cosi  = -100.;
1230
1231   Double_t ptTrig  = aodParticle->Pt();
1232   Double_t phiTrig = aodParticle->Phi();
1233   Double_t etaTrig = aodParticle->Eta();
1234   Double_t pxTrig  = aodParticle->Px();
1235   Double_t pyTrig  = aodParticle->Py();
1236
1237
1238   Int_t indexPhoton1 = -1 ;
1239   Int_t indexPhoton2 = -1 ;
1240   Double_t ptDecay1 = 0. ;
1241   Double_t pxDecay1  = 0. ;
1242   Double_t pyDecay1  = 0. ;
1243   Double_t phiDecay1  = 0. ;
1244   Double_t ptDecay2 = 0. ;
1245   Double_t pxDecay2  = 0. ;
1246   Double_t pyDecay2  = 0. ;
1247   Double_t phiDecay2  = 0. ;
1248
1249   Double_t ratDecay1  = -100.;
1250   Double_t ratDecay2  = -100.;
1251   Float_t deltaphi = -100. ;
1252   Float_t deltaphiDecay1 = -100. ;
1253   Float_t deltaphiDecay2 = -100. ;
1254   TObjArray * clusters = 0x0 ;
1255   TLorentzVector photonMom ;
1256   if(fPi0Trigger){
1257     indexPhoton1 = aodParticle->GetCaloLabel (0);
1258     indexPhoton2 = aodParticle->GetCaloLabel (1);
1259     if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
1260
1261     if(indexPhoton1!=-1 && indexPhoton2!=-1){
1262       if(aodParticle->GetDetector()=="EMCAL") clusters = GetAODEMCAL() ;
1263       else                                    clusters = GetAODPHOS() ;
1264       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
1265         AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));
1266         photon->GetMomentum(photonMom,GetVertex(0)) ;
1267         if(photon->GetID()==indexPhoton1) {
1268           ptDecay1  = photonMom.Pt();
1269           pxDecay1  = photonMom.Px();
1270           pyDecay1  = photonMom.Py();
1271           phiDecay1 = photonMom.Phi();
1272         }
1273         if(photon->GetID()==indexPhoton2) {
1274           ptDecay2  = photonMom.Pt();
1275           pxDecay2  = photonMom.Px();
1276           pyDecay2  = photonMom.Py();
1277           phiDecay2 = photonMom.Phi();
1278         }
1279         if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
1280       } //photonAOD loop
1281     } //index of decay photons found
1282     if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig, ptDecay2/ptTrig);
1283   } //make decay-hadron correlation
1284
1285   TObjArray * refpi0    =0x0;
1286   Int_t nrefs = 0;
1287
1288   //Loop on stored AOD pi0
1289   Int_t naod = pi0list->GetEntriesFast();
1290   if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() -  aod branch entries %d\n", naod);
1291   for(Int_t iaod = 0; iaod < naod ; iaod++){
1292     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (pi0list->At(iaod));
1293
1294     Int_t evtIndex2 = 0 ;
1295     Int_t evtIndex3 = 0 ;
1296     if (GetMixedEvent()) {
1297       evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
1298       evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
1299
1300       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
1301         continue ;
1302     }
1303
1304     //Int_t pdg = pi0->GetPdg();
1305     //if(pdg != AliCaloPID::kPi0) continue;
1306
1307     pt  = pi0->Pt();
1308     px  = pi0->Px();
1309     py  = pi0->Py();
1310     if(pt < GetMinPt() || pt > GetMaxPt()) continue ;
1311     //jumped out this event if near side associated partile pt larger than trigger
1312     if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2())  break ;
1313
1314     //Selection within angular range
1315     phi = pi0->Phi();
1316     //Float_t deltaphi = TMath::Abs(phiTrig-phi);
1317     //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;
1318
1319     if(bFillHisto){
1320
1321       deltaphi = phiTrig-phi;
1322       if(deltaphi<-TMath::PiOver2()) deltaphi+=TMath::TwoPi();
1323       if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();
1324
1325       rat = pt/ptTrig ;
1326       phi = pi0->Phi() ;
1327       eta = pi0->Eta() ;
1328       xE   = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
1329       if(xE <0.)xE =-xE;
1330       cosi = TMath::Log(1/xE);
1331
1332       if(fPi0Trigger){
1333         if(indexPhoton1!=-1 && indexPhoton2!=-1){
1334           if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
1335           if(ptDecay2) ratDecay2 = pt/ptDecay2 ;
1336           deltaphiDecay1 = phiDecay1-phi;
1337           deltaphiDecay2 = phiDecay2-phi;
1338           if(deltaphiDecay1< -TMath::PiOver2()) deltaphiDecay1+=TMath::TwoPi();
1339           if(deltaphiDecay1>3*TMath::PiOver2()) deltaphiDecay1-=TMath::TwoPi();
1340           if(deltaphiDecay2< -TMath::PiOver2()) deltaphiDecay2+=TMath::TwoPi();
1341           if(deltaphiDecay2>3*TMath::PiOver2()) deltaphiDecay2-=TMath::TwoPi();
1342           fhDeltaPhiDecayNeutral->Fill(ptDecay1, deltaphiDecay1);
1343           fhDeltaPhiDecayNeutral->Fill(ptDecay2, deltaphiDecay2);
1344           if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);
1345           if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )
1346             fhPtImbalanceDecayNeutral->Fill(ptDecay1,ratDecay1);
1347           if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )
1348             fhPtImbalanceDecayNeutral->Fill(ptDecay2,ratDecay2);
1349           if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
1350         }
1351       } //do decay-hadron correlation
1352
1353       fhEtaNeutral->Fill(pt,eta);
1354       fhPhiNeutral->Fill(pt,phi);
1355       fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);
1356       fhDeltaPhiNeutral->Fill(ptTrig,deltaphi);
1357       fhDeltaPhiDeltaEtaNeutral->Fill(deltaphi,etaTrig-eta);
1358
1359       //delta phi cut for correlation
1360       if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {
1361         fhDeltaPhiNeutralPt->Fill(pt,deltaphi);
1362         fhPtImbalanceNeutral->Fill(ptTrig,rat);
1363         fhPtHbpNeutral->Fill(ptTrig,cosi);
1364       }
1365       else {
1366         fhDeltaPhiUeNeutralPt->Fill(pt,deltaphi);
1367         fhPtImbalanceUeNeutral->Fill(ptTrig,rat);
1368         fhPtHbpUeNeutral->Fill(ptTrig,cosi);
1369       }
1370       //several UE calculation
1371       if(fMakeSeveralUE){
1372         if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){
1373           fhDeltaPhiUeLeftNeutral->Fill(pt,deltaphi);
1374           fhPtImbalanceUeLeftNeutral->Fill(ptTrig,rat);
1375           fhPtHbpUeLeftNeutral->Fill(ptTrig,cosi);
1376         }
1377         if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){
1378           fhDeltaPhiUeRightNeutral->Fill(pt,deltaphi);
1379           fhPtImbalanceUeRightNeutral->Fill(ptTrig,rat);
1380           fhPtHbpUeRightNeutral->Fill(ptTrig,cosi);
1381         }
1382       } //several UE calculation
1383           }
1384     else{
1385       nrefs++;
1386       if(nrefs==1){
1387         refpi0 = new TObjArray(0);
1388         refpi0->SetName(GetAODObjArrayName()+"Pi0s");
1389         refpi0->SetOwner(kFALSE);
1390       }
1391       refpi0->Add(pi0);
1392     }//put references in trigger AOD
1393
1394      //if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
1395
1396     }//loop
1397 }
1398
1399
1400 //____________________________________________________________________________
1401 //Bool_t  AliAnaParticleHadronCorrelation::SelectCluster(AliVCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) {
1402 //  //Select cluster depending on its pid and acceptance selections
1403 //
1404 //  //Skip matched clusters with tracks
1405 //  if(IsTrackMatched(calo)) return kFALSE;
1406 //
1407 //  TString detector = "";
1408 //  if     (calo->IsPHOS())  detector= "PHOS";
1409 //  else if(calo->IsEMCAL()) detector= "EMCAL";
1410 //
1411 //  //Check PID
1412 //  calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
1413 //  pdg = AliCaloPID::kPhoton;
1414 //  if(IsCaloPIDOn()){
1415 //    //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
1416 //    //or redo PID, recommended option for EMCal.
1417 //
1418 //    if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
1419 //      pdg = GetCaloPID()->GetPdg(detector,calo->GetPID(),mom.E());//PID with weights
1420 //    else
1421 //      pdg = GetCaloPID()->GetPdg(detector,mom,calo);//PID recalculated
1422 //
1423 //    if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);
1424 //
1425 //    //If it does not pass pid, skip
1426 //    if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) {
1427 //      return kFALSE ;
1428 //    }
1429 //  }//PID on
1430 //
1431 //  //Check acceptance selection
1432 //  if(IsFiducialCutOn()){
1433 //    Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,detector) ;
1434 //    if(! in ) return kFALSE ;
1435 //  }
1436 //
1437 //  if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
1438 //
1439 //  return kTRUE;
1440 //
1441 //}