PWGGA/CaloTrackCorrelations/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
  16 //_________________________________________________________________________
  17 // Class for the analysis of particle - hadron correlations
  18 // Particle (for example direct gamma) must be found in a previous analysis
  19 //-- Author: Gustavo Conesa (LNF-INFN)
  20
  21 //  Modified by Yaxian Mao:
  22 // 1. add the UE subtraction for corrlation study
  23 // 2. change the correlation variable
  24 // 3. Only use leading particle(cluster/track) as trigger for correlation (2010/07/02)
  25 // 4. Make decay photon-hadron correlations where decay contribute pi0 mass (2010/09/09)
  26 // 5. fill the pout to extract kt at the end, also to study charge asymmetry(2010/10/06)
  27 // 6. Add the possibility for event selection analysis based on vertex and multiplicity bins (10/10/2010)
  28 // 7. change the way of delta phi cut for UE study due to memory issue (reduce histograms)
  29 // 8. Add the possibility to request the absolute leading particle at the near side or not, set trigger bins, general clean-up (08/2011)
  30 //////////////////////////////////////////////////////////////////////////////
  31
  32
  33 // --- ROOT system ---
  34 //#include "TClonesArray.h"
  35 #include "TClass.h"
  36 #include "TMath.h"
  37 #include "TH3D.h"
  38 #include "TDatabasePDG.h"
  39
  40 //---- ANALYSIS system ----
  41 #include "AliNeutralMesonSelection.h"
  42 #include "AliAnaParticleHadronCorrelation.h"
  43 #include "AliCaloTrackReader.h"
  44 #include "AliCaloPID.h"
  45 #include "AliAODPWG4ParticleCorrelation.h"
  46 #include "AliFiducialCut.h"
  47 #include "AliVTrack.h"
  48 #include "AliVCluster.h"
  49 #include "AliMCAnalysisUtils.h"
  50 #include "TParticle.h"
  51 #include "AliStack.h"
  52 #include "AliAODMCParticle.h"
  53 #include "AliMixedEvent.h"
  54
  55 ClassImp(AliAnaParticleHadronCorrelation)
  56
  57
  58 //___________________________________________________________________
  59   AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation():
  60     AliAnaCaloTrackCorrBaseClass(),
  61     fMinTriggerPt(0.),
  62     fMaxAssocPt(1000.),             fMinAssocPt(0.),
  63     fDeltaPhiMaxCut(0.),            fDeltaPhiMinCut(0.),
  64     fSelectIsolated(0),             fMakeSeveralUE(0),
  65     fUeDeltaPhiMaxCut(0.),          fUeDeltaPhiMinCut(0.),
  66     fPi0AODBranchName(""),          fNeutralCorr(0),
  67     fPi0Trigger(0),                 fMakeAbsoluteLeading(0),
  68     fLeadingTriggerIndex(-1),
  69     fNAssocPtBins(0),               fAssocPtBinLimit(),
  70     //Histograms
  71     fhPtLeading(0),                 fhPhiLeading(0),
  72     fhEtaLeading(0),                fhDeltaPhiDeltaEtaCharged(0),
  73     fhPhiCharged(0),                fhEtaCharged(0),
  74     fhDeltaPhiCharged(0),           fhDeltaEtaCharged(0),
  75     fhDeltaPhiChargedPt(0),         fhDeltaPhiUeChargedPt(0),
  76     fhPtImbalanceCharged(0),        fhPtImbalanceUeCharged(0),
  77     fhPtImbalancePosCharged(0),     fhPtImbalanceNegCharged(0),
  78     fhPtHbpCharged(0),              fhPtHbpUeCharged(0),
  79     fhDeltaPhiUeLeftCharged(0),     fhDeltaPhiUeRightCharged(0),
  80     fhPtImbalanceUeLeftCharged(0),  fhPtImbalanceUeRightCharged(0),
  81     fhPtHbpUeLeftCharged(0),        fhPtHbpUeRightCharged(0),
  82     fhPtTrigPout(0),                fhPtTrigCharged(0),
  83     fhTrigDeltaPhiCharged(0x0),     fhTrigDeltaEtaCharged(0x0),
  84     fhTrigCorr(0x0),                fhTrigUeCorr(0x0),
  85     fhAssocPt(0),                   fhAssocPtBkg(0),
  86     fhDeltaPhiAssocPtBin(0),        fhDeltaPhiAssocPtBinHMPID(0),fhDeltaPhiAssocPtBinHMPIDAcc(0),
  87     fhDeltaPhiBradAssocPtBin(0),    fhDeltaPhiBrad(0),
  88     fhXEAssocPtBin(0),              fhXE(0),
  89     fhDeltaPhiDeltaEtaNeutral(0),
  90     fhPhiNeutral(0),                fhEtaNeutral(0),
  91     fhDeltaPhiNeutral(0),           fhDeltaEtaNeutral(0),
  92     fhDeltaPhiNeutralPt(0),         fhDeltaPhiUeNeutralPt(0),
  93     fhPtImbalanceNeutral(0),        fhPtImbalanceUeNeutral(0),
  94     fhPtHbpNeutral(0),              fhPtHbpUeNeutral(0),
  95     fhDeltaPhiUeLeftNeutral(0),     fhDeltaPhiUeRightNeutral(0),
  96     fhPtImbalanceUeLeftNeutral(0),  fhPtImbalanceUeRightNeutral(0),
  97     fhPtHbpUeLeftNeutral(0),        fhPtHbpUeRightNeutral(0),
  98     fhPtPi0DecayRatio(0),
  99     fhDeltaPhiDecayCharged(0),      fhPtImbalanceDecayCharged(0),
 100     fhDeltaPhiDecayNeutral(0),      fhPtImbalanceDecayNeutral(0),
 101     fh2phiLeadingParticle(0x0),
 102     fhMCLeadingCount(0),
 103     fhMCEtaCharged(0),              fhMCPhiCharged(0),
 104     fhMCDeltaEtaCharged(0),         fhMCDeltaPhiCharged(0x0),
 105     fhMCDeltaPhiDeltaEtaCharged(0), fhMCDeltaPhiChargedPt(0),
 106     fhMCPtImbalanceCharged(0),
 107     fhMCPtHbpCharged(0),
 108     fhMCPtTrigPout(0),
 109     fhMCPtAssocDeltaPhi(0)
 110 {
 111   //Default Ctor
 112
 113   //Initialize parameters
 114   InitParameters();
 115 }
 116
 117 //____________________________________________________________
 118 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
 119 {
 120   //Save parameters used for analysis
 121   TString parList ; //this will be list of parameters used for this analysis.
 122   const Int_t buffersize = 560;
 123   char onePar[buffersize] ;
 124
 125   snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;
 126   parList+=onePar ;
 127   snprintf(onePar,buffersize," Pt Trigger > %3.2f ",    fMinTriggerPt) ;
 128   parList+=onePar ;
 129   snprintf(onePar,buffersize," %3.2f < Pt associated < %3.2f ", fMinAssocPt,   fMaxAssocPt) ;
 130   parList+=onePar ;
 131   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-Hadron < %3.2f ",    fDeltaPhiMinCut,   fDeltaPhiMaxCut) ;
 132   parList+=onePar ;
 133   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-UeHadron <  %3.2f ", fUeDeltaPhiMinCut, fUeDeltaPhiMaxCut) ;
 134   parList+=onePar ;
 135   snprintf(onePar,buffersize,"Isolated Trigger?  %d\n",    fSelectIsolated) ;
 136   parList+=onePar ;
 137   snprintf(onePar,buffersize,"Several UE?  %d\n",          fMakeSeveralUE) ;
 138   parList+=onePar ;
 139   snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ", fPi0AODBranchName.Data());
 140   parList+=onePar ;
 141   snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  %d", fPi0Trigger) ;
 142   parList+=onePar ;
 143   snprintf(onePar,buffersize,"Select absolute leading for cluster triggers ?  %d\n", fMakeAbsoluteLeading) ;
 144   parList+=onePar ;
 145   snprintf(onePar,buffersize,"Associated particle pt bins  %d: ", fNAssocPtBins) ;
 146   parList+=onePar ;
 147   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
 148     snprintf(onePar,buffersize,"bin %d = [%2.1f,%2.1f];", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
 149   }
 150   parList+=onePar ;
 151
 152   //Get parameters set in base class.
 153   parList += GetBaseParametersList() ;
 154
 155   //Get parameters set in PID class.
 156   //parList += GetCaloPID()->GetPIDParametersList() ;
 157
 158   //Get parameters set in FiducialCut class (not available yet)
 159   //parlist += GetFidCut()->GetFidCutParametersList()
 160
 161   return new TObjString(parList) ;
 162
 163 }
 164
 165 //________________________________________________________________
 166 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
 167 {
 168
 169   // Create histograms to be saved in output file and
 170   // store them in fOutputContainer
 171   TList * outputContainer = new TList() ;
 172   outputContainer->SetName("CorrelationHistos") ;
 173
 174   Int_t   nptbins = GetHistogramRanges()->GetHistoPtBins(); Int_t  nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins();
 175   Float_t ptmax   = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax  = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax();
 176   Float_t ptmin   = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin  = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin();
 177
 178   fhPtLeading  = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax);
 179   fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
 180
 181   fhPhiLeading  = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 182   fhPhiLeading->SetYTitle("#phi (rad)");
 183
 184   fhEtaLeading  = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 185   fhEtaLeading->SetYTitle("#eta ");
 186
 187   outputContainer->Add(fhPtLeading);
 188   outputContainer->Add(fhPhiLeading);
 189   outputContainer->Add(fhEtaLeading);
 190
 191   //Correlation with charged hadrons
 192   if(GetReader()->IsCTSSwitchedOn()) {
 193     fhDeltaPhiDeltaEtaCharged  = new TH2F
 194     ("DeltaPhiDeltaEtaCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
 195      140,-2.,5.,200,-2,2);
 196     fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
 197     fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
 198
 199     fhPhiCharged  = new TH2F
 200     ("PhiCharged","#phi_{h^{#pm}}  vs p_{T #pm}",
 201      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 202     fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
 203     fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
 204
 205     fhEtaCharged  = new TH2F
 206     ("EtaCharged","#eta_{h^{#pm}}  vs p_{T #pm}",
 207      nptbins,ptmin,ptmax,netabins,etamin,etamax);
 208     fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
 209     fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
 210
 211     fhDeltaPhiCharged  = new TH2F
 212     ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
 213      nptbins,ptmin,ptmax,140,-2.,5.);
 214     fhDeltaPhiCharged->SetYTitle("#Delta #phi");
 215     fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
 216
 217     fhDeltaPhiChargedPt  = new TH2F
 218     ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
 219      nptbins,ptmin,ptmax,140,-2.,5.);
 220     fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
 221     fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 222
 223     fhDeltaPhiUeChargedPt  = new TH2F
 224     ("DeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
 225      nptbins,ptmin,ptmax,140,-2.,5.);
 226     fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
 227     fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 228
 229     fhDeltaEtaCharged  = new TH2F
 230     ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
 231      nptbins,ptmin,ptmax,200,-2,2);
 232     fhDeltaEtaCharged->SetYTitle("#Delta #eta");
 233     fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
 234
 235     fhPtImbalanceCharged  =
 236     new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
 237              nptbins,ptmin,ptmax,200,0.,2.);
 238     fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
 239     fhPtImbalanceCharged->SetXTitle("p_{T trigger}");
 240
 241     fhPtImbalanceUeCharged  =
 242     new TH2F("CorrelationUeCharged","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger}",
 243              nptbins,ptmin,ptmax,200,0.,2.);
 244     fhPtImbalanceUeCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
 245     fhPtImbalanceUeCharged->SetXTitle("p_{T trigger}");
 246
 247     fhPtImbalancePosCharged  =
 248     new TH2F("CorrelationPositiveCharged","z_{trigger h^{+}} = p_{T h^{+}} / p_{T trigger}",
 249              nptbins,ptmin,ptmax,200,0.,2.);
 250     fhPtImbalancePosCharged->SetYTitle("z_{trigger h^{+}}");
 251     fhPtImbalancePosCharged->SetXTitle("p_{T trigger}");
 252
 253     fhPtImbalanceNegCharged  =
 254     new TH2F("CorrelationNegativeCharged","z_{trigger h^{-}} = p_{T h^{-}} / p_{T trigger}",
 255              nptbins,ptmin,ptmax,200,0.,2.);
 256     fhPtImbalanceNegCharged->SetYTitle("z_{trigger h^{-}}");
 257     fhPtImbalanceNegCharged->SetXTitle("p_{T trigger}");
 258
 259     fhPtHbpCharged  =
 260     new TH2F("HbpCharged","#xi = ln(1/x_{E}) with charged hadrons",
 261              nptbins,ptmin,ptmax,200,0.,10.);
 262     fhPtHbpCharged->SetYTitle("ln(1/x_{E})");
 263     fhPtHbpCharged->SetXTitle("p_{T trigger}");
 264
 265     fhPtHbpUeCharged  =
 266     new TH2F("HbpUeCharged","#xi = ln(1/x_{E}) with charged hadrons",
 267              nptbins,ptmin,ptmax,200,0.,10.);
 268     fhPtHbpUeCharged->SetYTitle("ln(1/x_{E})");
 269     fhPtHbpUeCharged->SetXTitle("p_{T trigger}");
 270
 271     fhPtTrigPout  =
 272     new TH2F("PtTrigPout","Pout with triggers",
 273              nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
 274     fhPtTrigPout->SetYTitle("p_{out} (GeV/c)");
 275     fhPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
 276
 277     fhPtTrigCharged  =
 278     new TH2F("PtTrigCharged","trgger and charged tracks pt distribution",
 279              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 280     fhPtTrigCharged->SetYTitle("p_{T h^{#pm}} (GeV/c)");
 281     fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)");
 282
 283     outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
 284     outputContainer->Add(fhPhiCharged) ;
 285     outputContainer->Add(fhEtaCharged) ;
 286     outputContainer->Add(fhDeltaPhiCharged) ;
 287     outputContainer->Add(fhDeltaEtaCharged) ;
 288     outputContainer->Add(fhDeltaPhiChargedPt) ;
 289     outputContainer->Add(fhDeltaPhiUeChargedPt) ;
 290     outputContainer->Add(fhPtImbalanceCharged) ;
 291     outputContainer->Add(fhPtImbalancePosCharged) ;
 292     outputContainer->Add(fhPtImbalanceNegCharged) ;
 293     outputContainer->Add(fhPtImbalanceUeCharged) ;
 294     outputContainer->Add(fhPtHbpCharged) ;
 295     outputContainer->Add(fhPtHbpUeCharged) ;
 296     outputContainer->Add(fhPtTrigPout) ;
 297     outputContainer->Add(fhPtTrigCharged) ;
 298
 299     if(DoEventSelect()){
 300       Int_t nMultiBins = GetMultiBin();
 301       fhTrigDeltaPhiCharged = new TH2F*[nMultiBins] ;
 302       fhTrigDeltaEtaCharged = new TH2F*[nMultiBins] ;
 303       fhTrigCorr            = new TH2F*[nMultiBins];
 304       fhTrigUeCorr          = new TH2F*[nMultiBins];
 305       for(Int_t im=0; im<nMultiBins; im++){
 306         fhTrigDeltaPhiCharged[im]  = new TH2F
 307         (Form("fhTrigDeltaPhiCharged_%d",im),Form("fhTrigDeltaPhiCharged_%d",im), nptbins,ptmin,ptmax, 140,-2.,5.);
 308         fhTrigDeltaPhiCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
 309         fhTrigDeltaPhiCharged[im]->SetYTitle("#Delta #phi");
 310         fhTrigDeltaEtaCharged[im]  = new TH2F
 311         (Form("fhTrigDeltaEtaCharged_%d",im),Form("fhTrigDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, 200,-2,2);
 312         fhTrigDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
 313         fhTrigDeltaEtaCharged[im]->SetYTitle("#Delta #eta");
 314         fhTrigCorr[im]  = new TH2F
 315         (Form("fhTrigPtCorr_%d",im),Form("fhTrigPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 316         fhTrigCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
 317         fhTrigCorr[im]->SetXTitle("p_{T trigger}");
 318         fhTrigUeCorr[im]  = new TH2F
 319         (Form("fhTrigPtUeCorr_%d",im),Form("fhTrigPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 320         fhTrigUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
 321         fhTrigUeCorr[im]->SetXTitle("p_{T trigger}");
 322
 323         outputContainer->Add(fhTrigDeltaPhiCharged[im]) ;
 324         outputContainer->Add(fhTrigDeltaEtaCharged[im]) ;
 325         outputContainer->Add(fhTrigCorr[im]);
 326         outputContainer->Add(fhTrigUeCorr[im]);
 327
 328       }
 329     }
 330
 331     fhAssocPt           = new TH2F("fhAssocPt", " Trigger p_{T} vs associated hadron p_{T}",
 332                                    nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
 333     fhAssocPt->SetXTitle("p_{T trigger}");
 334     fhAssocPt->SetYTitle("p_{T associated}");
 335     outputContainer->Add(fhAssocPt) ;
 336
 337     fhAssocPtBkg        = new TH2F("fhAssocPtBkg", " Trigger p_{T} vs associated hadron p_{T} from background",
 338                                    nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
 339     fhAssocPtBkg->SetXTitle("p_{T trigger}");
 340     fhAssocPtBkg->SetYTitle("p_{T associated}");
 341     outputContainer->Add(fhAssocPtBkg) ;
 342
 343     fhDeltaPhiBrad = new TH2F("fhDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs p_{T trigger} ",
 344                               nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
 345     fhDeltaPhiBrad->SetXTitle("p_{T trigger}");
 346     fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
 347
 348     outputContainer->Add(fhDeltaPhiBrad) ;
 349
 350     fhXE       = new TH2F("fhXE", "x_{E} vs p_{T trigger}",
 351                           nptbins, ptmin, ptmax,50, 0.0, 2.0);
 352     fhXE->SetXTitle("p_{T trigger}");
 353     fhXE->SetYTitle("x_{E}");
 354
 355     outputContainer->Add(fhXE);
 356
 357
 358     fhDeltaPhiAssocPtBin     = new TH2F*[fNAssocPtBins] ;
 359     fhDeltaPhiAssocPtBinHMPID= new TH2F*[fNAssocPtBins] ;
 360     fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins] ;
 361     fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins] ;
 362     fhXEAssocPtBin           = new TH2F*[fNAssocPtBins] ;
 363     for(Int_t i = 0 ; i < fNAssocPtBins ; i++){
 364       fhDeltaPhiAssocPtBin[i] = new TH2F(Form("fhDeltaPhiPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 365                                          Form("#Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 366                                          nptbins, ptmin, ptmax,140,-2.,5.);
 367       fhDeltaPhiAssocPtBin[i]->SetXTitle("p_{T trigger}");
 368       fhDeltaPhiAssocPtBin[i]->SetYTitle("#Delta #phi");
 369
 370       fhDeltaPhiAssocPtBinHMPID[i] = new TH2F(Form("fhDeltaPhiPtAssocPt%2.1f_%2.1fHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 371                                          Form("#Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f], with track having HMPID signal", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 372                                          nptbins, ptmin, ptmax,140,-2.,5.);
 373       fhDeltaPhiAssocPtBinHMPID[i]->SetXTitle("p_{T trigger}");
 374       fhDeltaPhiAssocPtBinHMPID[i]->SetYTitle("#Delta #phi");
 375
 376       fhDeltaPhiAssocPtBinHMPIDAcc[i] = new TH2F(Form("fhDeltaPhiPtAssocPt%2.1f_%2.1fHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 377                                               Form("#Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f], with track within 5<phi<20 deg", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 378                                               nptbins, ptmin, ptmax,140,-2.,5.);
 379       fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetXTitle("p_{T trigger}");
 380       fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetYTitle("#Delta #phi");
 381
 382       fhDeltaPhiBradAssocPtBin[i] = new TH2F(Form("fhDeltaPhiBradPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 383                                              Form("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 384                                              nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
 385       fhDeltaPhiBradAssocPtBin[i]->SetXTitle("p_{T trigger}");
 386       fhDeltaPhiBradAssocPtBin[i]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
 387
 388
 389       fhXEAssocPtBin[i]       = new TH2F(Form("fhXEAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 390                                          Form("x_{E} vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 391                                          nptbins, ptmin, ptmax,50, 0.0, 2.0);
 392       fhXEAssocPtBin[i]->SetXTitle("p_{T trigger}");
 393       fhXEAssocPtBin[i]->SetYTitle("x_{E}");
 394
 395       outputContainer->Add(fhDeltaPhiAssocPtBin[i]) ;
 396       outputContainer->Add(fhDeltaPhiAssocPtBinHMPID[i]) ;
 397       outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[i]) ;
 398       outputContainer->Add(fhDeltaPhiBradAssocPtBin[i]) ;
 399       outputContainer->Add(fhXEAssocPtBin[i]);
 400     }
 401
 402     if(fPi0Trigger){
 403       fhPtPi0DecayRatio  = new TH2F
 404       ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay",
 405        nptbins,ptmin,ptmax, 100,0.,2.);
 406       fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");
 407       fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");
 408
 409       fhDeltaPhiDecayCharged  = new TH2F
 410       ("DeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",
 411        nptbins,ptmin,ptmax,140,-2.,5.);
 412       fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");
 413       fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");
 414
 415       fhPtImbalanceDecayCharged  =
 416       new TH2F("CorrelationDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",
 417                nptbins,ptmin,ptmax,200,0.,2.);
 418       fhPtImbalanceDecayCharged->SetYTitle("z_{decay h^{#pm}}");
 419       fhPtImbalanceDecayCharged->SetXTitle("p_{T decay}");
 420
 421       outputContainer->Add(fhPtPi0DecayRatio) ;
 422       outputContainer->Add(fhDeltaPhiDecayCharged) ;
 423       outputContainer->Add(fhPtImbalanceDecayCharged) ;
 424     }
 425
 426
 427     if(fMakeSeveralUE){
 428       fhDeltaPhiUeLeftCharged  = new TH2F
 429       ("DeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
 430        nptbins,ptmin,ptmax,140,-2.,5.);
 431       fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
 432       fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 433       outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
 434
 435       fhDeltaPhiUeRightCharged  = new TH2F
 436       ("DeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
 437        nptbins,ptmin,ptmax,140,-2.,5.);
 438       fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
 439       fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 440       outputContainer->Add(fhDeltaPhiUeRightCharged) ;
 441
 442       fhPtImbalanceUeLeftCharged  =
 443       new TH2F("CorrelationUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",
 444                nptbins,ptmin,ptmax,200,0.,2.);
 445       fhPtImbalanceUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
 446       fhPtImbalanceUeLeftCharged->SetXTitle("p_{T trigger}");
 447       outputContainer->Add(fhPtImbalanceUeLeftCharged) ;
 448
 449       fhPtImbalanceUeRightCharged  =
 450       new TH2F("CorrelationUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",
 451                nptbins,ptmin,ptmax,200,0.,2.);
 452       fhPtImbalanceUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
 453       fhPtImbalanceUeRightCharged->SetXTitle("p_{T trigger}");
 454       outputContainer->Add(fhPtImbalanceUeRightCharged) ;
 455
 456       fhPtHbpUeLeftCharged  =
 457       new TH2F("HbpUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",
 458                nptbins,ptmin,ptmax,200,0.,10.);
 459       fhPtHbpUeLeftCharged->SetYTitle("ln(1/x_{E})");
 460       fhPtHbpUeLeftCharged->SetXTitle("p_{T trigger}");
 461       outputContainer->Add(fhPtHbpUeLeftCharged) ;
 462
 463       fhPtHbpUeRightCharged  =
 464       new TH2F("HbpUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",
 465                nptbins,ptmin,ptmax,200,0.,10.);
 466       fhPtHbpUeRightCharged->SetYTitle("ln(1/x_{E})");
 467       fhPtHbpUeRightCharged->SetXTitle("p_{T trigger}");
 468       outputContainer->Add(fhPtHbpUeRightCharged) ;
 469
 470     }
 471   }  //Correlation with charged hadrons
 472
 473   //Correlation with neutral hadrons
 474   if(fNeutralCorr){
 475
 476     fhDeltaPhiDeltaEtaNeutral  = new TH2F
 477     ("DeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
 478      140,-2.,5.,200,-2,2);
 479     fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");
 480     fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
 481
 482     fhPhiNeutral  = new TH2F
 483     ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T #pi^{0}}",
 484      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 485     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
 486     fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
 487
 488     fhEtaNeutral  = new TH2F
 489     ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T #pi^{0}}",
 490      nptbins,ptmin,ptmax,netabins,etamin,etamax);
 491     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
 492     fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
 493
 494     fhDeltaPhiNeutral  = new TH2F
 495     ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
 496      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 497     fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
 498     fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
 499
 500     fhDeltaPhiNeutralPt  = new TH2F
 501     ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
 502      nptbins,ptmin,ptmax,140,-2.,5.);
 503     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
 504     fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
 505
 506     fhDeltaPhiUeNeutralPt  = new TH2F
 507     ("DeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
 508      nptbins,ptmin,ptmax,140,-2.,5.);
 509     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
 510     fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
 511
 512     fhDeltaEtaNeutral  = new TH2F
 513     ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
 514      nptbins,ptmin,ptmax,200,-2,2);
 515     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
 516     fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
 517
 518     fhPtImbalanceNeutral  =
 519     new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
 520              nptbins,ptmin,ptmax,200,0.,2.);
 521     fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");
 522     fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");
 523
 524     fhPtImbalanceUeNeutral  =
 525     new TH2F("CorrelationUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",
 526              nptbins,ptmin,ptmax,200,0.,2.);
 527     fhPtImbalanceUeNeutral->SetYTitle("z_{trigger #pi^{0}}");
 528     fhPtImbalanceUeNeutral->SetXTitle("p_{T trigger}");
 529
 530     fhPtHbpNeutral  =
 531     new TH2F("HbpNeutral","#xi = ln(1/x_{E}) with neutral particles",
 532              nptbins,ptmin,ptmax,200,0.,10.);
 533     fhPtHbpNeutral->SetYTitle("ln(1/x_{E})");
 534     fhPtHbpNeutral->SetXTitle("p_{T trigger}");
 535
 536     fhPtHbpUeNeutral  =
 537     new TH2F("HbpUeNeutral","#xi = ln(1/x_{E}) with neutral particles",
 538              nptbins,ptmin,ptmax,200,0.,10.);
 539     fhPtHbpUeNeutral->SetYTitle("ln(1/x_{E})");
 540     fhPtHbpUeNeutral->SetXTitle("p_{T trigger}");
 541
 542
 543     outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
 544     outputContainer->Add(fhPhiNeutral) ;
 545     outputContainer->Add(fhEtaNeutral) ;
 546     outputContainer->Add(fhDeltaPhiNeutral) ;
 547     outputContainer->Add(fhDeltaPhiNeutralPt) ;
 548     outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
 549     outputContainer->Add(fhDeltaEtaNeutral) ;
 550     outputContainer->Add(fhPtImbalanceNeutral) ;
 551     outputContainer->Add(fhPtImbalanceUeNeutral) ;
 552     outputContainer->Add(fhPtHbpNeutral) ;
 553     outputContainer->Add(fhPtHbpUeNeutral) ;
 554
 555     if(fPi0Trigger){
 556       fhDeltaPhiDecayNeutral  = new TH2F
 557       ("DeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",
 558        nptbins,ptmin,ptmax,140,-2.,5.);
 559       fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");
 560       fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");
 561
 562       fhPtImbalanceDecayNeutral  =
 563       new TH2F("CorrelationDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",
 564                nptbins,ptmin,ptmax,200,0.,2.);
 565       fhPtImbalanceDecayNeutral->SetYTitle("z_{decay h^{0}}");
 566       fhPtImbalanceDecayNeutral->SetXTitle("p_{T decay}");
 567
 568       outputContainer->Add(fhDeltaPhiDecayNeutral) ;
 569       outputContainer->Add(fhPtImbalanceDecayNeutral) ;
 570     }
 571
 572
 573     if(fMakeSeveralUE){
 574       fhDeltaPhiUeLeftNeutral  = new TH2F
 575       ("DeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",
 576        nptbins,ptmin,ptmax,140,-2.,5.);
 577       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
 578       fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
 579       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
 580
 581       fhDeltaPhiUeRightNeutral  = new TH2F
 582       ("DeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",
 583        nptbins,ptmin,ptmax,140,-2.,5.);
 584       fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
 585       fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
 586       outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
 587
 588       fhPtImbalanceUeLeftNeutral  =
 589       new TH2F("CorrelationUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
 590                nptbins,ptmin,ptmax,140,0.,2.);
 591       fhPtImbalanceUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
 592       fhPtImbalanceUeLeftNeutral->SetXTitle("p_{T trigger}");
 593       outputContainer->Add(fhPtImbalanceUeLeftNeutral) ;
 594
 595       fhPtImbalanceUeRightNeutral  =
 596       new TH2F("CorrelationUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
 597                nptbins,ptmin,ptmax,200,0.,2.);
 598       fhPtImbalanceUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
 599       fhPtImbalanceUeRightNeutral->SetXTitle("p_{T trigger}");
 600       outputContainer->Add(fhPtImbalanceUeRightNeutral) ;
 601
 602       fhPtHbpUeLeftNeutral  =
 603       new TH2F("HbpUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",
 604                nptbins,ptmin,ptmax,200,0.,10.);
 605       fhPtHbpUeLeftNeutral->SetYTitle("ln(1/x_{E})");
 606       fhPtHbpUeLeftNeutral->SetXTitle("p_{T trigger}");
 607       outputContainer->Add(fhPtHbpUeLeftNeutral) ;
 608
 609       fhPtHbpUeRightNeutral  =
 610       new TH2F("HbpUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",
 611                nptbins,ptmin,ptmax,200,0.,10.);
 612       fhPtHbpUeRightNeutral->SetYTitle("ln(1/x_{E})");
 613       fhPtHbpUeRightNeutral->SetXTitle("p_{T trigger}");
 614       outputContainer->Add(fhPtHbpUeRightNeutral) ;
 615
 616     }
 617
 618   }//Correlation with neutral hadrons
 619
 620   //if data is MC, fill more histograms
 621   if(IsDataMC()){
 622
 623     fh2phiLeadingParticle=new TH2F("fh2phiLeadingParticle","#phi resolustion for trigger particles",nptbins,ptmin,ptmax,100,-1,1);
 624     fh2phiLeadingParticle->GetXaxis()->SetTitle("p_{T gen Leading} (GeV/c)");
 625     fh2phiLeadingParticle->GetYaxis()->SetTitle("(#phi_{rec}-#phi_{gen})/#phi_{gen}");
 626
 627     fhMCLeadingCount=new TH1F("MCLeadingTriggerCount","MCLeadingTriggerCount",nptbins,ptmin,ptmax);
 628     fhMCLeadingCount->SetXTitle("p_{T trig}");
 629
 630     fhMCEtaCharged  = new TH2F
 631     ("MCEtaCharged","MC #eta_{h^{#pm}}  vs p_{T #pm}",
 632      nptbins,ptmin,ptmax,netabins,etamin,etamax);
 633     fhMCEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
 634     fhMCEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
 635
 636     fhMCPhiCharged  = new TH2F
 637     ("MCPhiCharged","#MC phi_{h^{#pm}}  vs p_{T #pm}",
 638      200,ptmin,ptmax,nphibins,phimin,phimax);
 639     fhMCPhiCharged->SetYTitle("MC #phi_{h^{#pm}} (rad)");
 640     fhMCPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
 641
 642     fhMCDeltaPhiDeltaEtaCharged  = new TH2F
 643     ("MCDeltaPhiDeltaEtaCharged","#MC phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
 644      140,-2.,5.,200,-2,2);
 645     fhMCDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
 646     fhMCDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
 647
 648     fhMCDeltaEtaCharged  = new TH2F
 649     ("MCDeltaEtaCharged","MC #eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger} and p_{T assoc}",
 650      nptbins,ptmin,ptmax,200,-2,2);
 651     fhMCDeltaEtaCharged->SetYTitle("#Delta #eta");
 652     fhMCDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
 653
 654     fhMCDeltaPhiCharged  = new TH2F
 655     ("MCDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
 656      nptbins,ptmin,ptmax,140,-2.,5.);
 657     fhMCDeltaPhiCharged->SetYTitle("#Delta #phi");
 658     fhMCDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
 659
 660     fhMCDeltaPhiChargedPt  = new TH2F
 661     ("MCDeltaPhiChargedPt","MC #phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
 662      nptbins,ptmin,ptmax,140,-2.,5.);
 663     fhMCDeltaPhiChargedPt->SetYTitle("#Delta #phi");
 664     fhMCDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 665
 666     fhMCPtImbalanceCharged  =
 667     new TH2F("MCCorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",
 668              nptbins,ptmin,ptmax,200,0.,2.);
 669     fhMCPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");
 670     fhMCPtImbalanceCharged->SetXTitle("p_{T trigger}");
 671
 672     fhMCPtHbpCharged  =
 673     new TH2F("MCHbpCharged","MC #xi = ln(1/x_{E}) with charged hadrons",
 674              nptbins,ptmin,ptmax,200,0.,10.);
 675     fhMCPtHbpCharged->SetYTitle("ln(1/x_{E})");
 676     fhMCPtHbpCharged->SetXTitle("p_{T trigger}");
 677
 678     fhMCPtTrigPout  =
 679     new TH2F("MCPtTrigPout","AOD MC Pout with triggers",
 680              nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
 681     fhMCPtTrigPout->SetYTitle("p_{out} (GeV/c)");
 682     fhMCPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
 683
 684     fhMCPtAssocDeltaPhi  =
 685     new TH2F("fhMCPtAssocDeltaPhi","AOD MC delta phi with associated charged hadrons",
 686              nptbins,ptmin,ptmax,140,-2.,5.);
 687     fhMCPtAssocDeltaPhi->SetYTitle("#Delta #phi");
 688     fhMCPtAssocDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 689
 690     outputContainer->Add(fh2phiLeadingParticle);
 691     outputContainer->Add(fhMCLeadingCount);
 692     outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged);
 693     outputContainer->Add(fhMCPhiCharged) ;
 694     outputContainer->Add(fhMCEtaCharged) ;
 695     outputContainer->Add(fhMCDeltaEtaCharged) ;
 696     outputContainer->Add(fhMCDeltaPhiCharged) ;
 697
 698     outputContainer->Add(fhMCDeltaPhiChargedPt) ;
 699     outputContainer->Add(fhMCPtImbalanceCharged) ;
 700     outputContainer->Add(fhMCPtHbpCharged) ;
 701     outputContainer->Add(fhMCPtTrigPout) ;
 702     outputContainer->Add(fhMCPtAssocDeltaPhi) ;
 703   } //for MC histogram
 704
 705
 706   return outputContainer;
 707
 708 }
 709
 710 //____________________________________________________
 711 void AliAnaParticleHadronCorrelation::InitParameters()
 712 {
 713
 714   //Initialize the parameters of the analysis.
 715   SetInputAODName("PWG4Particle");
 716   SetAODObjArrayName("Hadrons");
 717   AddToHistogramsName("AnaHadronCorr_");
 718
 719   SetPtCutRange(0.,300);
 720   fDeltaPhiMinCut       = 1.5 ;
 721   fDeltaPhiMaxCut       = 4.5 ;
 722   fSelectIsolated       = kFALSE;
 723   fMakeSeveralUE        = kFALSE;
 724   fUeDeltaPhiMinCut     = 1. ;
 725   fUeDeltaPhiMaxCut     = 1.5 ;
 726   fNeutralCorr          = kFALSE ;
 727   fPi0Trigger           = kFALSE ;
 728   fMakeAbsoluteLeading  = kTRUE;
 729
 730   fNAssocPtBins         = 7  ;
 731   fAssocPtBinLimit[0]   = 1.5  ;
 732   fAssocPtBinLimit[1]   = 3.  ;
 733   fAssocPtBinLimit[2]   = 5.  ;
 734   fAssocPtBinLimit[3]   = 7.  ;
 735   fAssocPtBinLimit[4]   = 9. ;
 736   fAssocPtBinLimit[5]   = 12. ;
 737   fAssocPtBinLimit[6]   = 15. ;
 738   fAssocPtBinLimit[7]   = 20. ;
 739   fAssocPtBinLimit[8]   = 100.;
 740   fAssocPtBinLimit[9]   = 200.;
 741
 742 }
 743
 744 //__________________________________________________________
 745 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
 746 {
 747   //Particle-Hadron Correlation Analysis, fill AODs
 748
 749   if(!GetInputAODBranch()){
 750     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
 751     abort();
 752   }
 753
 754   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
 755     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());
 756     abort();
 757   }
 758
 759   if(GetDebug() > 1){
 760     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");
 761     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
 762     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n",   GetCTSTracks()    ->GetEntriesFast());
 763     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
 764     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n",  GetPHOSClusters() ->GetEntriesFast());
 765   }
 766
 767   //Get the vertex and check it is not too large in z
 768   Double_t v[3] = {0,0,0}; //vertex ;
 769   GetReader()->GetVertex(v);
 770   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
 771
 772   //Loop on stored AOD particles, find leading trigger
 773   Double_t ptTrig      = fMinTriggerPt ;
 774   fLeadingTriggerIndex = -1 ;
 775   Int_t    naod        = GetInputAODBranch()->GetEntriesFast() ;
 776   for(Int_t iaod = 0; iaod < naod ; iaod++){
 777     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 778
 779     //vertex cut in case of mixing
 780     if (GetMixedEvent()) {
 781       Int_t evt=-1;
 782       Int_t id =-1;
 783       if     (particle->GetCaloLabel(0)  >= 0 ){
 784         id=particle->GetCaloLabel(0);
 785         if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
 786       }
 787       else if(particle->GetTrackLabel(0) >= 0 ){
 788         id=particle->GetTrackLabel(0);
 789         if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
 790       }
 791       else continue;
 792
 793       if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
 794         return ;
 795     }
 796
 797     // find the leading particles with highest momentum
 798     if (particle->Pt() > ptTrig) {
 799       ptTrig               = particle->Pt() ;
 800       fLeadingTriggerIndex = iaod ;
 801     }
 802   }// finish search of leading trigger particle
 803
 804
 805   //Do correlation with leading particle
 806   if(fLeadingTriggerIndex >= 0){
 807
 808     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(fLeadingTriggerIndex));
 809
 810     //check if the particle is isolated or if we want to take the isolation into account
 811     if(OnlyIsolated() && !particle->IsIsolated()) return;
 812
 813     //Make correlation with charged hadrons
 814     Bool_t okcharged = kTRUE;
 815     Bool_t okneutral = kTRUE;
 816     if(GetReader()->IsCTSSwitchedOn() )
 817       okcharged = MakeChargedCorrelation(particle, GetCTSTracks(),kFALSE);
 818
 819     TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
 820     if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
 821       okneutral = MakeNeutralCorrelation(particle, pi0list,kFALSE);
 822
 823   }//Correlate leading
 824
 825   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
 826
 827 }
 828
 829 //_________________________________________________________________
 830 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
 831 {
 832   //Particle-Hadron Correlation Analysis, fill histograms
 833
 834   if(!GetInputAODBranch()){
 835     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
 836     abort();
 837   }
 838
 839   if(GetDebug() > 1){
 840     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
 841     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
 842   }
 843
 844   //Get the vertex and check it is not too large in z
 845   Double_t v[3] = {0,0,0}; //vertex ;
 846   GetReader()->GetVertex(v);
 847   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
 848
 849   //Loop on stored AOD particles, find leading
 850   Double_t ptTrig    = fMinTriggerPt;
 851   if(fLeadingTriggerIndex < 0){//Search leading if not done before
 852     Int_t    naod      = GetInputAODBranch()->GetEntriesFast() ;
 853     for(Int_t iaod = 0; iaod < naod ; iaod++){   //loop on input trigger AOD file
 854       AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
 855       //vertex cut in case of mixing
 856       if (GetMixedEvent()) {
 857         Int_t evt=-1;
 858         Int_t id =-1;
 859         if     (particle->GetCaloLabel(0)  >= 0 ){
 860           id=particle->GetCaloLabel(0);
 861           if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;
 862         }
 863         else if(particle->GetTrackLabel(0) >= 0 ){
 864           id=particle->GetTrackLabel(0);
 865           if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;
 866         }
 867         else continue;
 868
 869         if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut())
 870           return ;
 871       }
 872
 873       //check if the particle is isolated or if we want to take the isolation into account
 874       if(OnlyIsolated() && !particle->IsIsolated()) continue;
 875
 876       //check if inside the vertex cut
 877       //find the leading particles with highest momentum
 878       if (particle->Pt() > ptTrig) {
 879         ptTrig               = particle->Pt() ;
 880         fLeadingTriggerIndex = iaod ;
 881       }
 882     }//finish search of leading trigger particle
 883   }//Search leading if not done before
 884
 885   if(fLeadingTriggerIndex >= 0 ){ //using trigger particle to do correlations
 886     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(fLeadingTriggerIndex));
 887
 888     //check if the particle is isolated or if we want to take the isolation into account
 889     if(OnlyIsolated() && !particle->IsIsolated()) return;
 890
 891     //Make correlation with charged hadrons
 892     Bool_t okcharged = kTRUE;
 893     Bool_t okneutral = kTRUE;
 894     if(GetReader()->IsCTSSwitchedOn() ){
 895       okcharged = MakeChargedCorrelation(particle, GetCTSTracks(),kTRUE);
 896       if(IsDataMC()){
 897         MakeMCChargedCorrelation(particle);
 898       }
 899     }
 900
 901     TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
 902     if(fNeutralCorr && pi0list){
 903       if(pi0list->GetEntriesFast() > 0)
 904         okneutral = MakeNeutralCorrelation(particle, pi0list,kTRUE);
 905     }
 906
 907     // Fill leading particle histogram if correlation went well and
 908     // no problem was found, like not absolute leading, or bad vertex in mixing.
 909     if(okcharged && okneutral){
 910       fhPtLeading->Fill(particle->Pt());
 911       Float_t phi = particle->Phi();
 912       if(phi<0)phi+=TMath::TwoPi();
 913       fhPhiLeading->Fill(particle->Pt(), phi);
 914       fhEtaLeading->Fill(particle->Pt(), particle->Eta());
 915     }//ok charged && neutral
 916   }//Aod branch loop
 917
 918   //Reinit for next event
 919   fLeadingTriggerIndex = -1;
 920
 921   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
 922 }
 923
 924 //___________________________________________________________________________________________________________
 925 Bool_t  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle,
 926                                                                 const TObjArray* pl, const Bool_t bFillHisto)
 927 {
 928   // Charged Hadron Correlation Analysis
 929   if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
 930
 931   Int_t evtIndex11   = -1 ; //cluster trigger or pi0 trigger
 932   Int_t evtIndex12   = -1 ; // pi0 trigger
 933   Int_t evtIndex13   = -1 ; // charged trigger
 934   Int_t indexPhoton1 = -1 ;
 935   Int_t indexPhoton2 = -1 ;
 936
 937   Double_t v[3]      = {0,0,0}; //vertex ;
 938   GetReader()->GetVertex(v);
 939
 940   if (GetMixedEvent()) {
 941     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
 942     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
 943     evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
 944   }
 945
 946   Double_t phiTrig = aodParticle->Phi();
 947   Double_t etaTrig = aodParticle->Eta();
 948   Double_t ptTrig  = aodParticle->Pt();
 949
 950   Double_t pt             = -100. ;
 951   Double_t px             = -100. ;
 952   Double_t py             = -100. ;
 953   Double_t rat            = -100. ;
 954   Double_t xE             = -100. ;
 955   Double_t cosi           = -100. ;
 956   Double_t phi            = -100. ;
 957   Double_t eta            = -100. ;
 958   Double_t pout           = -100. ;
 959
 960   Double_t ptDecay1       = 0. ;
 961   Double_t pxDecay1       = 0. ;
 962   Double_t pyDecay1       = 0. ;
 963   Double_t phiDecay1      = 0. ;
 964   Double_t ptDecay2       = 0. ;
 965   Double_t pxDecay2       = 0. ;
 966   Double_t pyDecay2       = 0. ;
 967   Double_t phiDecay2      = 0. ;
 968
 969   Double_t ratDecay1      = -100. ;
 970   Double_t ratDecay2      = -100. ;
 971   Double_t deltaPhi       = -100. ;
 972   Double_t deltaPhiOrg    = -100. ;
 973   Double_t deltaPhiDecay1 = -100. ;
 974   Double_t deltaPhiDecay2 = -100. ;
 975
 976   TVector3 p3;
 977   TLorentzVector photonMom ;
 978   TObjArray * clusters  = 0x0 ;
 979   TObjArray * reftracks = 0x0;
 980   Int_t nrefs           = 0;
 981   Int_t nTracks         = GetCTSTracks()->GetEntriesFast() ;
 982
 983   if(fPi0Trigger){
 984     indexPhoton1 = aodParticle->GetCaloLabel (0);
 985     indexPhoton2 = aodParticle->GetCaloLabel (1);
 986     if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
 987
 988     if(indexPhoton1!=-1 && indexPhoton2!=-1){
 989       if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
 990       else                                    clusters = GetPHOSClusters()  ;
 991       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
 992         AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));
 993         photon->GetMomentum(photonMom,GetVertex(0)) ;
 994         if(photon->GetID()==indexPhoton1) {
 995           ptDecay1  = photonMom.Pt();
 996           pxDecay1  = photonMom.Px();
 997           pyDecay1  = photonMom.Py();
 998           phiDecay1 = photonMom.Phi();
 999           if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig);
1000         }
1001         if(photon->GetID()==indexPhoton2) {
1002           ptDecay2  = photonMom.Pt();
1003           pxDecay2  = photonMom.Px();
1004           pyDecay2  = photonMom.Py();
1005           phiDecay2 = photonMom.Phi();
1006           if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay2/ptTrig);
1007         }
1008         if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
1009       } //cluster loop
1010     } //index of decay photons found
1011   } //make decay-hadron correlation
1012
1013   //Track loop, select tracks with good pt, phi and fill AODs or histograms
1014   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){
1015     AliVTrack * track = (AliVTrack *) (pl->At(ipr)) ;
1016
1017     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
1018     p3.SetXYZ(mom[0],mom[1],mom[2]);
1019     pt   = p3.Pt();
1020     px   = p3.Px();
1021     py   = p3.Py();
1022     eta  = p3.Eta();
1023     phi  = p3.Phi() ;
1024     if(phi < 0) phi+=TMath::TwoPi();
1025
1026     //Select only hadrons in pt range
1027     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
1028
1029     //remove trigger itself for correlation when use charged triggers
1030     if( track->GetID() == aodParticle->GetTrackLabel(0) || track->GetID() == aodParticle->GetTrackLabel(1) ||
1031         track->GetID() == aodParticle->GetTrackLabel(2) || track->GetID() == aodParticle->GetTrackLabel(3)   )
1032       continue ;
1033
1034     if(IsFiducialCutOn()){
1035       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
1036       if(! in ) continue ;
1037     }
1038
1039     //jump out this event if near side associated particle pt larger than trigger
1040     if (fMakeAbsoluteLeading){
1041       if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2())  return kFALSE;
1042     }
1043
1044     //Only for mixed event
1045     Int_t evtIndex2 = 0 ;
1046     if (GetMixedEvent()) {
1047       evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
1048       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
1049         continue ;
1050       //vertex cut
1051       if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
1052         return kFALSE;
1053     }
1054
1055     if(fPi0Trigger){
1056       if(indexPhoton1!=-1 && indexPhoton2!=-1){
1057         if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
1058         if(ptDecay2) ratDecay2 = pt/ptDecay2 ;
1059         deltaPhiDecay1 = phiDecay1-phi;
1060         deltaPhiDecay2 = phiDecay2-phi;
1061         if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
1062         if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
1063         if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
1064         if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
1065       }
1066     } //do decay-hadron correlation
1067
1068     //Selection within angular range
1069     deltaPhi    = phiTrig-phi;
1070     deltaPhiOrg = deltaPhi;
1071     if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
1072     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
1073
1074     rat  = pt/ptTrig ;
1075     pout = pt*TMath::Sin(deltaPhi) ;
1076     xE   =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
1077     //if(xE <0.)xE =-xE;
1078     cosi = -100;
1079     if(xE > 0 ) cosi = TMath::Log(1./xE);
1080
1081     if(GetDebug() > 2)
1082       printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Charged hadron: pt %f, phi %f, phi trigger %f. Cuts:  delta phi  %2.2f < %2.2f < %2.2f, pT Trig min %2.2f \n",
1083              pt,phi, phiTrig,fDeltaPhiMinCut, deltaPhi, fDeltaPhiMaxCut, fMinTriggerPt);
1084
1085     // Fill Histograms
1086     if(bFillHisto){
1087
1088       Int_t    assocBin   = -1;
1089       for(Int_t i = 0 ; i < fNAssocPtBins ; i++){
1090         if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i;
1091       }
1092
1093       fhAssocPt->Fill(ptTrig,pt);
1094
1095       //if(xE >= 1 ) printf("pTTrig = %f, pTHadron = %f, xE = %f\n",ptTrig,pt, xE);
1096       fhXE->Fill(ptTrig, xE) ;
1097
1098       if(TMath::Cos(deltaPhi) < 0 && assocBin >= 0 )//away side
1099         fhXEAssocPtBin[assocBin]->Fill(ptTrig, xE) ;
1100
1101       //Hardcoded values, BAD, FIXME
1102       Double_t  dphiBrad = atan2(sin(deltaPhiOrg), cos(deltaPhiOrg))/TMath::Pi();//-1 to 1
1103       if(TMath::Abs(dphiBrad)>0.325 && TMath::Abs(dphiBrad)<0.475){
1104         fhAssocPtBkg->Fill(ptTrig, pt);
1105       }
1106
1107       if(dphiBrad<-1./3) dphiBrad += 2;
1108       fhDeltaPhiBrad->Fill(ptTrig, dphiBrad);
1109
1110       if(assocBin>=0){
1111         fhDeltaPhiBradAssocPtBin[assocBin]->Fill(ptTrig, dphiBrad);
1112         fhDeltaPhiAssocPtBin    [assocBin]->Fill(ptTrig, deltaPhi);
1113         if(track->GetHMPIDsignal()>0){
1114           printf("Track pt %f with HMPID signal %f \n",pt,track->GetHMPIDsignal());
1115           fhDeltaPhiAssocPtBinHMPID[assocBin]->Fill(ptTrig, deltaPhi);
1116         }
1117
1118         if(phi > 5*TMath::DegToRad() && phi < 20*TMath::DegToRad()){
1119           //printf("Track pt %f in HMPID acceptance phi %f \n ",pt,phi*TMath::RadToDeg() );
1120           fhDeltaPhiAssocPtBinHMPIDAcc[assocBin]->Fill(ptTrig, deltaPhi);
1121
1122         }
1123
1124       }
1125
1126       fhEtaCharged     ->Fill(pt,eta);
1127       fhPhiCharged     ->Fill(pt,phi);
1128       fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);
1129       fhDeltaPhiCharged->Fill(ptTrig, deltaPhi);
1130       if(pt >2 ) fhDeltaPhiDeltaEtaCharged->Fill(deltaPhi,aodParticle->Eta()-eta);
1131
1132       if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
1133       //fill different multiplicity histogram
1134       if(DoEventSelect()){
1135         for(Int_t im=0; im<GetMultiBin(); im++){
1136           if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1)){
1137             fhTrigDeltaPhiCharged[im]->Fill(ptTrig,deltaPhi);
1138             fhTrigDeltaEtaCharged[im]->Fill(ptTrig,aodParticle->Eta()-eta);
1139           }
1140         }
1141       }
1142
1143       //delta phi cut for correlation
1144       if      ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   ) {
1145         fhDeltaPhiChargedPt ->Fill(pt,deltaPhi);
1146         fhPtImbalanceCharged->Fill(ptTrig,xE);
1147         fhPtHbpCharged ->Fill(ptTrig, cosi);
1148         fhPtTrigPout   ->Fill(ptTrig, pout) ;
1149         fhPtTrigCharged->Fill(ptTrig, pt) ;
1150         if(track->Charge() > 0) fhPtImbalancePosCharged->Fill(ptTrig,xE) ;
1151         else                    fhPtImbalanceNegCharged->Fill(ptTrig,xE) ;
1152         //fill different multiplicity histogram
1153         if(DoEventSelect()){
1154           for(Int_t im=0; im<GetMultiBin(); im++){
1155             if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
1156               fhTrigCorr[im]->Fill(ptTrig,xE);
1157           }
1158         } //multiplicity events selection
1159       } //delta phi cut for correlation
1160       else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) ) { //UE study
1161         fhDeltaPhiUeChargedPt->Fill(pt,deltaPhi);
1162         Double_t randomphi = gRandom->Uniform(TMath::Pi()/2,3*TMath::Pi()/2);
1163         Double_t uexE = -(pt/ptTrig)*TMath::Cos(randomphi);
1164         if(uexE < 0.) uexE = -uexE;
1165         if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - xe = %f, uexE = %f \n", xE, uexE);
1166         fhPtImbalanceUeCharged->Fill(ptTrig,uexE);
1167         if(uexE>0)fhPtHbpUeCharged->Fill(ptTrig,TMath::Log(1/uexE));
1168         if(DoEventSelect()){
1169           for(Int_t im=0; im<GetMultiBin(); im++){
1170             if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
1171               fhTrigUeCorr[im]->Fill(ptTrig,xE);
1172           }
1173         } //multiplicity events selection
1174
1175       } //UE study
1176
1177       if(fPi0Trigger){
1178         if(indexPhoton1!=-1 && indexPhoton2!=-1){
1179           fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaPhiDecay1);
1180           fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaPhiDecay2);
1181           if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
1182           if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
1183             fhPtImbalanceDecayCharged->Fill(ptDecay1,ratDecay1);
1184           if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
1185             fhPtImbalanceDecayCharged->Fill(ptDecay2,ratDecay2);
1186           if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
1187         } //index of decay photons found
1188       } //make decay-hadron correlation
1189
1190       //several UE calculation
1191       if(fMakeSeveralUE){
1192         if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut)){
1193           fhDeltaPhiUeLeftCharged->Fill(pt,deltaPhi);
1194           fhPtImbalanceUeLeftCharged->Fill(ptTrig,rat);
1195           fhPtHbpUeLeftCharged->Fill(ptTrig,cosi);
1196         }
1197         if((deltaPhi>fUeDeltaPhiMinCut) && (deltaPhi <fUeDeltaPhiMaxCut)){
1198           fhDeltaPhiUeRightCharged->Fill(pt,deltaPhi);
1199           fhPtImbalanceUeRightCharged->Fill(ptTrig,rat);
1200           fhPtHbpUeRightCharged->Fill(ptTrig,cosi);
1201
1202         }
1203       } //several UE calculation
1204
1205       //Fill leading particle histogram
1206       fhPtLeading->Fill(ptTrig);
1207       if(phiTrig<0)phiTrig+=TMath::TwoPi();
1208       fhPhiLeading->Fill(ptTrig, phiTrig);
1209       fhEtaLeading->Fill(ptTrig, etaTrig);
1210
1211     } //Fill histogram
1212     else{
1213       nrefs++;
1214       if(nrefs==1){
1215         reftracks = new TObjArray(0);
1216         TString trackname = Form("%s+Tracks", GetAODObjArrayName().Data());
1217         reftracks->SetName(trackname.Data());
1218         reftracks->SetOwner(kFALSE);
1219       }
1220       reftracks->Add(track);
1221     }//aod particle loop
1222   }// track loop
1223
1224   //Fill AOD with reference tracks, if not filling histograms
1225   if(!bFillHisto && reftracks) {
1226     aodParticle->AddObjArray(reftracks);
1227   }
1228
1229   return kTRUE;
1230
1231 }
1232
1233 //________________________________________________________________________________________________________________
1234 Bool_t  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle,
1235                                                                 const TObjArray* pi0list, const Bool_t bFillHisto)
1236 {
1237   // Neutral Pion Correlation Analysis
1238   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",pi0list->GetEntriesFast());
1239
1240   Int_t evtIndex11 = 0 ;
1241   Int_t evtIndex12 = 0 ;
1242   if (GetMixedEvent()) {
1243     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
1244     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
1245   }
1246
1247   Double_t pt   = -100. ;
1248   Double_t px   = -100. ;
1249   Double_t py   = -100. ;
1250   Double_t rat  = -100. ;
1251   Double_t phi  = -100. ;
1252   Double_t eta  = -100. ;
1253   Double_t xE   = -100. ;
1254   Double_t cosi = -100. ;
1255
1256   Double_t ptTrig  = aodParticle->Pt();
1257   Double_t phiTrig = aodParticle->Phi();
1258   Double_t etaTrig = aodParticle->Eta();
1259   //Double_t pxTrig  = aodParticle->Px();
1260   //Double_t pyTrig  = aodParticle->Py();
1261
1262   Int_t indexPhoton1 =-1  ;
1263   Int_t indexPhoton2 =-1  ;
1264   Double_t ptDecay1  = 0. ;
1265   Double_t pxDecay1  = 0. ;
1266   Double_t pyDecay1  = 0. ;
1267   Double_t phiDecay1 = 0. ;
1268   Double_t ptDecay2  = 0. ;
1269   Double_t pxDecay2  = 0. ;
1270   Double_t pyDecay2  = 0. ;
1271   Double_t phiDecay2 = 0. ;
1272
1273   Double_t ratDecay1      = -100. ;
1274   Double_t ratDecay2      = -100. ;
1275   Double_t deltaPhi       = -100. ;
1276   Double_t deltaPhiDecay1 = -100. ;
1277   Double_t deltaPhiDecay2 = -100. ;
1278
1279   TObjArray * clusters = 0x0 ;
1280   TLorentzVector photonMom ;
1281
1282   if(fPi0Trigger){
1283     indexPhoton1 = aodParticle->GetCaloLabel (0);
1284     indexPhoton2 = aodParticle->GetCaloLabel (1);
1285     if(GetDebug() > 1)
1286       printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
1287
1288     if(indexPhoton1!=-1 && indexPhoton2!=-1){
1289       if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
1290       else                                    clusters = GetPHOSClusters() ;
1291       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){
1292         AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));
1293         photon->GetMomentum(photonMom,GetVertex(0)) ;
1294         if(photon->GetID()==indexPhoton1) {
1295           ptDecay1  = photonMom.Pt();
1296           pxDecay1  = photonMom.Px();
1297           pyDecay1  = photonMom.Py();
1298           phiDecay1 = photonMom.Phi();
1299         }
1300         if(photon->GetID()==indexPhoton2) {
1301           ptDecay2  = photonMom.Pt();
1302           pxDecay2  = photonMom.Px();
1303           pyDecay2  = photonMom.Py();
1304           phiDecay2 = photonMom.Phi();
1305         }
1306         if(GetDebug() > 1)
1307           printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() - Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);
1308       } //photonAOD loop
1309     } //index of decay photons found
1310     if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig, ptDecay2/ptTrig);
1311   } //make decay-hadron correlation
1312
1313   TObjArray * refpi0    =0x0;
1314   Int_t nrefs = 0;
1315
1316   //Loop on stored AOD pi0
1317   Int_t naod = pi0list->GetEntriesFast();
1318   if(GetDebug() > 0)
1319     printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() -  aod branch entries %d\n", naod);
1320   for(Int_t iaod = 0; iaod < naod ; iaod++){
1321     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (pi0list->At(iaod));
1322
1323     Int_t evtIndex2 = 0 ;
1324     Int_t evtIndex3 = 0 ;
1325     if (GetMixedEvent()) {
1326       evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
1327       evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
1328
1329       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 ||
1330           evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
1331         continue ;
1332     }
1333
1334     //Int_t pdg = pi0->GetPdg();
1335     //if(pdg != AliCaloPID::kPi0) continue;
1336
1337     pt  = pi0->Pt();
1338     px  = pi0->Px();
1339     py  = pi0->Py();
1340
1341     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
1342
1343     //jumped out this event if near side associated particle pt larger than trigger
1344     if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2())  break ;
1345
1346     //Selection within angular range
1347     phi = pi0->Phi();
1348     //Float_t deltaPhi = TMath::Abs(phiTrig-phi);
1349     //if( (deltaPhi < fDeltaPhiMinCut) || ( deltaPhi > fDeltaPhiMaxCut) ) continue ;
1350
1351     if(bFillHisto){
1352
1353       deltaPhi = phiTrig-phi;
1354       if(deltaPhi<-TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
1355       if(deltaPhi>3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
1356
1357       rat = pt/ptTrig ;
1358       phi = pi0->Phi() ;
1359       eta = pi0->Eta() ;
1360
1361       rat  = pt/ptTrig ;
1362       xE   =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
1363       //if(xE <0.)xE =-xE;
1364       cosi = -100;
1365       if(xE > 0 ) cosi = TMath::Log(1./xE);
1366
1367       if(fPi0Trigger){
1368         if(indexPhoton1!=-1 && indexPhoton2!=-1){
1369           if(ptDecay1) ratDecay1 = pt/ptDecay1 ;
1370           if(ptDecay2) ratDecay2 = pt/ptDecay2 ;
1371           deltaPhiDecay1 = phiDecay1-phi;
1372           deltaPhiDecay2 = phiDecay2-phi;
1373           if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
1374           if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
1375           if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
1376           if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
1377           fhDeltaPhiDecayNeutral->Fill(ptDecay1, deltaPhiDecay1);
1378           fhDeltaPhiDecayNeutral->Fill(ptDecay2, deltaPhiDecay2);
1379           if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
1380           if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
1381             fhPtImbalanceDecayNeutral->Fill(ptDecay1,ratDecay1);
1382           if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
1383             fhPtImbalanceDecayNeutral->Fill(ptDecay2,ratDecay2);
1384           if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);
1385         }
1386       } //do decay-hadron correlation
1387
1388       fhEtaNeutral->Fill(pt,eta);
1389       fhPhiNeutral->Fill(pt,phi);
1390       fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);
1391       fhDeltaPhiNeutral->Fill(ptTrig,deltaPhi);
1392       fhDeltaPhiDeltaEtaNeutral->Fill(deltaPhi,etaTrig-eta);
1393
1394       //delta phi cut for correlation
1395       if( (deltaPhi > fDeltaPhiMinCut) && ( deltaPhi < fDeltaPhiMaxCut) ) {
1396         fhDeltaPhiNeutralPt->Fill(pt,deltaPhi);
1397         fhPtImbalanceNeutral->Fill(ptTrig,rat);
1398         fhPtHbpNeutral->Fill(ptTrig,cosi);
1399       }
1400       else {
1401         fhDeltaPhiUeNeutralPt->Fill(pt,deltaPhi);
1402         fhPtImbalanceUeNeutral->Fill(ptTrig,rat);
1403         fhPtHbpUeNeutral->Fill(ptTrig,cosi);
1404       }
1405       //several UE calculation
1406       if(fMakeSeveralUE){
1407         if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut)){
1408           fhDeltaPhiUeLeftNeutral->Fill(pt,deltaPhi);
1409           fhPtImbalanceUeLeftNeutral->Fill(ptTrig,rat);
1410           fhPtHbpUeLeftNeutral->Fill(ptTrig,cosi);
1411         }
1412         if((deltaPhi>fUeDeltaPhiMinCut) && (deltaPhi <fUeDeltaPhiMaxCut)){
1413           fhDeltaPhiUeRightNeutral->Fill(pt,deltaPhi);
1414           fhPtImbalanceUeRightNeutral->Fill(ptTrig,rat);
1415           fhPtHbpUeRightNeutral->Fill(ptTrig,cosi);
1416         }
1417       } //several UE calculation
1418           }
1419     else{
1420       nrefs++;
1421       if(nrefs==1){
1422         refpi0 = new TObjArray(0);
1423         refpi0->SetName(GetAODObjArrayName()+"Pi0s");
1424         refpi0->SetOwner(kFALSE);
1425       }
1426       refpi0->Add(pi0);
1427     }//put references in trigger AOD
1428
1429     if(GetDebug() > 2 )
1430       printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
1431
1432   }//loop
1433
1434   return kTRUE;
1435 }
1436
1437 //_________________________________________________________________________________________________________
1438 void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
1439 {
1440   // Charged Hadron Correlation Analysis with MC information
1441   if(GetDebug()>1)
1442     printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation() - Make trigger particle - charged hadron correlation in AOD MC level\n");
1443
1444   AliStack         * stack        = 0x0 ;
1445   TParticle        * primary      = 0x0 ;
1446   TClonesArray     * mcparticles0 = 0x0 ;
1447   TClonesArray     * mcparticles  = 0x0 ;
1448   AliAODMCParticle * aodprimary   = 0x0 ;
1449
1450   Double_t eprim   = 0 ;
1451   Double_t ptprim  = 0 ;
1452   Double_t phiprim = 0 ;
1453   Double_t etaprim = 0 ;
1454   Double_t pxprim  = 0 ;
1455   Double_t pyprim  = 0 ;
1456   Double_t pzprim  = 0 ;
1457   Int_t    nTracks = 0 ;
1458   Int_t iParticle  = 0 ;
1459   Double_t charge  = 0.;
1460
1461   Double_t mcrat   =-100 ;
1462   Double_t mcxE    =-100 ;
1463   Double_t mccosi  =-100 ;
1464   Double_t mcdeltaPhi = -100. ;
1465
1466   //Track loop, select tracks with good pt, phi and fill AODs or histograms
1467   //Int_t currentIndex = -1 ;
1468   Double_t mcTrackPt  = 0 ;
1469   Double_t mcTrackPhi = 0 ;
1470   Double_t mcTrackEta = 0 ;
1471   Double_t mcTrackPx  = 0 ;
1472   Double_t mcTrackPy  = 0 ;
1473   Double_t mcTrackPz  = 0 ;
1474
1475   if(GetReader()->ReadStack()){
1476     nTracks = GetMCStack()->GetNtrack() ;
1477   }
1478   else nTracks = GetReader()->GetAODMCParticles()->GetEntriesFast() ;
1479   //Int_t trackIndex[nTracks];
1480
1481   Int_t label= aodParticle->GetLabel();
1482   if(label<0){
1483     printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** bad label ***:  label %d \n", label);
1484     return;
1485   }
1486
1487   if(GetReader()->ReadStack()){
1488     stack =  GetMCStack() ;
1489     if(!stack) {
1490       printf(" AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation- Stack not available, is the MC handler called? STOP\n");
1491       abort();
1492     }
1493
1494     nTracks=stack->GetNprimary();
1495     if(label >=  stack->GetNtrack()) {
1496       if(GetDebug() > 2)  printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
1497       return ;
1498     }
1499     primary = stack->Particle(label);
1500     if(!primary){
1501       printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no primary ***:  label %d \n", label);
1502       return;
1503     }
1504
1505     eprim    = primary->Energy();
1506     ptprim   = primary->Pt();
1507     phiprim  = primary->Phi();
1508     etaprim  = primary->Eta();
1509     pxprim   = primary->Px();
1510     pyprim   = primary->Py();
1511     pzprim   = primary->Pz();
1512
1513     if(primary){
1514
1515       for (iParticle = 0 ; iParticle <  nTracks ; iParticle++) {
1516         TParticle * particle = stack->Particle(iParticle);
1517         TLorentzVector momentum;
1518         //keep only final state particles
1519         if(particle->GetStatusCode()!=1) continue ;
1520         Int_t pdg = particle->GetPdgCode();
1521         charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1522         particle->Momentum(momentum);
1523
1524         //---------- Charged particles ----------------------
1525         if(charge != 0){
1526           //Particles in CTS acceptance
1527           Bool_t inCTS =  GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
1528           if(TMath::Abs(pdg) == 11 && stack->Particle(particle->GetFirstMother())->GetPdgCode()==22) continue ;
1529           if(inCTS)
1530           {
1531             mcTrackPt  = particle->Pt();
1532             mcTrackPhi = particle->Phi();
1533             mcTrackEta = particle->Eta();
1534             mcTrackPx  = particle->Px();
1535             mcTrackPy  = particle->Py();
1536             mcTrackPz  = particle->Pz();
1537             if(mcTrackPhi < 0) mcTrackPhi+=TMath::TwoPi();
1538
1539             //Select only hadrons in pt range
1540             if(mcTrackPt < fMinAssocPt || mcTrackPt > fMaxAssocPt) continue ;
1541
1542             //remove trigger itself for correlation when use charged triggers
1543             if(label==iParticle &&
1544                TMath::Abs(mcTrackPt -ptprim ) < 1e-6 &&
1545                TMath::Abs(mcTrackPhi-phiprim) < 1e-6 &&
1546                TMath::Abs(mcTrackEta-etaprim) < 1e-6)
1547               continue ;
1548
1549             //jumped out this event if near side associated partile pt larger than trigger
1550             if( mcTrackPt > ptprim && TMath::Abs(mcTrackPhi-phiprim) < TMath::PiOver2())
1551               return ;
1552
1553             mcdeltaPhi= phiprim-mcTrackPhi;
1554             if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
1555             if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
1556
1557             mcrat     = mcTrackPt/ptprim ;
1558             mcdeltaPhi= phiprim-mcTrackPhi;
1559             mcxE      =-mcTrackPt/ptprim*TMath::Cos(mcdeltaPhi);// -(mcTrackPx*pxprim+mcTrackPy*pyprim)/(ptprim*ptprim);
1560             mccosi    =-100 ;
1561             if(mcxE > 0 ) mccosi = TMath::Log(1./mcxE);
1562
1563             // printf("rat = %f, xE = %f, cosi =%f \n", rat, xE, cosi);
1564             // printf("phi = %f \n", phi);
1565
1566             //Selection within angular range
1567             if( mcdeltaPhi< -TMath::PiOver2())  mcdeltaPhi+=TMath::TwoPi();
1568             if( mcdeltaPhi>3*TMath::PiOver2())  mcdeltaPhi-=TMath::TwoPi();
1569             Double_t mcpout = mcTrackPt*TMath::Sin(mcdeltaPhi) ;
1570             if(GetDebug()>0 )
1571               printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation() - Charged hadron: track Pt %f, track Phi %f, phi trigger %f. Cuts:  delta phi  %2.2f < %2.2f < %2.2f, pT min %2.2f \n",
1572                      mcTrackPt,mcTrackPhi, phiprim,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut, GetMinPt());
1573             // Fill Histograms
1574             fhMCEtaCharged->Fill(mcTrackPt,mcTrackEta);
1575             fhMCPhiCharged->Fill(mcTrackPt,mcTrackPhi);
1576             fhMCDeltaEtaCharged->Fill(ptprim,etaprim-mcTrackEta);
1577             fhMCDeltaPhiCharged->Fill(ptprim,mcdeltaPhi);
1578             fhMCPtAssocDeltaPhi->Fill(mcTrackPt, mcdeltaPhi);
1579             fhMCDeltaPhiDeltaEtaCharged->Fill(mcdeltaPhi,etaprim-mcTrackEta);
1580
1581             //delta phi cut for correlation
1582             if( (mcdeltaPhi > fDeltaPhiMinCut) && ( mcdeltaPhi < fDeltaPhiMaxCut) ) {
1583               fhMCDeltaPhiChargedPt->Fill(mcTrackPt,mcdeltaPhi);
1584               fhMCPtImbalanceCharged->Fill(ptprim,mcxE);
1585               fhMCPtHbpCharged->Fill(ptprim,mccosi);
1586               fhMCPtTrigPout->Fill(ptprim, mcpout) ;
1587             }//delta phi cut for correlation
1588           } //tracks after cuts
1589         }//Charged
1590       } //track loop
1591     } //when the leading particles could trace back to MC
1592   } //ESD MC
1593   else if(GetReader()->ReadAODMCParticles()){
1594     //Get the list of MC particles
1595     mcparticles0 = GetReader()->GetAODMCParticles(0);
1596     if(!mcparticles0) return;
1597     if(label >=mcparticles0->GetEntriesFast()){
1598       if(GetDebug() > 2)
1599         printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** large label ***:  label %d, n tracks %d \n", label,mcparticles0->GetEntriesFast());
1600       return;
1601     }
1602     //Get the particle
1603     aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
1604     if(!aodprimary)  {
1605       printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no AOD primary ***:  label %d \n", label);
1606       return;
1607     }
1608
1609     ptprim  = aodprimary->Pt();
1610     phiprim = aodprimary->Phi();
1611     etaprim = aodprimary->Eta();
1612     pxprim  = aodprimary->Px();
1613     pyprim  = aodprimary->Py();
1614     pzprim  = aodprimary->Pz();
1615     if(aodprimary){
1616       mcparticles= GetReader()->GetAODMCParticles();
1617       for (Int_t i=0; i<nTracks;i++) {
1618         AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(i);
1619         if (!part->IsPhysicalPrimary()) continue;
1620         Int_t pdg = part->GetPdgCode();
1621         charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
1622         TLorentzVector momentum(part->Px(),part->Py(),part->Pz(),part->E());
1623         if(charge != 0){
1624           if(part->Pt()> GetReader()->GetCTSPtMin()){
1625             //Particles in CTS acceptance
1626             Bool_t inCTS =  GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
1627             Int_t indexmother=part->GetMother();
1628             if(indexmother>-1)
1629             {
1630               Int_t mPdg = ((AliAODMCParticle*) mcparticles->At(indexmother)) ->GetPdgCode();
1631               if(TMath::Abs(pdg) == 11 && mPdg == 22) continue;
1632             }
1633             if(inCTS)
1634             {
1635               mcTrackPt  =part->Pt();
1636               mcTrackPhi =part->Phi();
1637               mcTrackEta =part->Eta();
1638               mcTrackPx  =part->Px();
1639               mcTrackPy  =part->Py();
1640               mcTrackPz  =part->Pz();
1641               if(mcTrackPhi < 0) mcTrackPhi+=TMath::TwoPi();
1642
1643               //Select only hadrons in pt range
1644               if(mcTrackPt < GetMinPt() || mcTrackPt > GetMaxPt()) continue ;
1645
1646               //remove trigger itself for correlation when use charged triggers
1647               if(label==i &&
1648                  TMath::Abs(mcTrackPt -ptprim ) < 1e-6 &&
1649                  TMath::Abs(mcTrackPhi-phiprim) < 1e-6 &&
1650                  TMath::Abs(mcTrackEta-etaprim) < 1e-6) continue ;
1651
1652               //jumped out this event if near side associated partile pt larger than trigger
1653               if( mcTrackPt> ptprim && TMath::Abs(mcTrackPhi-phiprim)<TMath::PiOver2())
1654                 return ;
1655
1656               mcdeltaPhi= phiprim-mcTrackPhi;
1657               if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
1658               if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
1659
1660               mcrat     = mcTrackPt/ptprim ;
1661               mcxE      =-mcTrackPt/ptprim*TMath::Cos(mcdeltaPhi);// -(mcTrackPx*pxprim+mcTrackPy*pyprim)/(ptprim*ptprim);
1662               mccosi    =-100 ;
1663               if(mcxE > 0 ) mccosi = TMath::Log(1./mcxE);
1664
1665               //Selection within angular range
1666               if( mcdeltaPhi< -TMath::PiOver2())  mcdeltaPhi+=TMath::TwoPi();
1667               if( mcdeltaPhi>3*TMath::PiOver2())  mcdeltaPhi-=TMath::TwoPi();
1668               Double_t mcpout = mcTrackPt*TMath::Sin(mcdeltaPhi) ;
1669               if(GetDebug()>0)
1670                 printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation() - Charged hadron: track Pt %f, track Phi %f, phi trigger %f. Cuts:  delta phi  %2.2f < %2.2f < %2.2f, pT min %2.2f \n",
1671                        mcTrackPt,mcTrackPhi, phiprim,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut, GetMinPt());
1672
1673               // Fill Histograms
1674
1675               fhMCEtaCharged     ->Fill(mcTrackPt,mcTrackEta);
1676               fhMCPhiCharged     ->Fill(mcTrackPt,mcTrackPhi);
1677               fhMCDeltaEtaCharged->Fill(ptprim,etaprim-mcTrackEta);
1678               fhMCDeltaPhiCharged->Fill(ptprim,mcdeltaPhi);
1679               fhMCPtAssocDeltaPhi->Fill(mcTrackPt, mcdeltaPhi);
1680               fhMCDeltaPhiDeltaEtaCharged->Fill(mcdeltaPhi,etaprim-mcTrackEta);
1681
1682               //delta phi cut for correlation
1683               if( (mcdeltaPhi > fDeltaPhiMinCut) && ( mcdeltaPhi < fDeltaPhiMaxCut) ) {
1684                 fhMCDeltaPhiChargedPt ->Fill(mcTrackPt,mcdeltaPhi);
1685                 fhMCPtImbalanceCharged->Fill(ptprim,mcxE);
1686                 fhMCPtHbpCharged      ->Fill(ptprim,mccosi);
1687                 fhMCPtTrigPout        ->Fill(ptprim, mcpout) ;
1688               }//delta phi cut for correlation
1689
1690             } //tracks after cuts
1691
1692           } //with minimum pt cut
1693         } //only charged particles
1694       }  //MC particle loop
1695     } //when the leading particles could trace back to MC
1696   }// AOD MC
1697 }
1698
1699 //_____________________________________________________________________
1700 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
1701 {
1702
1703   //Print some relevant parameters set for the analysis
1704   if(! opt)
1705     return;
1706
1707   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1708   AliAnaCaloTrackCorrBaseClass::Print(" ");
1709   printf("Pt trigger           >  %3.2f\n", fMinTriggerPt) ;
1710   printf("Pt associated hadron <  %3.2f\n", fMaxAssocPt) ;
1711   printf("Pt associated hadron >  %3.2f\n", fMinAssocPt) ;
1712   printf("Phi trigger particle-Hadron   <  %3.2f\n", fDeltaPhiMaxCut) ;
1713   printf("Phi trigger particle-Hadron   >  %3.2f\n", fDeltaPhiMinCut) ;
1714   printf("Phi trigger particle-UeHadron <  %3.2f\n", fUeDeltaPhiMaxCut) ;
1715   printf("Phi trigger particle-UeHadron >  %3.2f\n", fUeDeltaPhiMinCut) ;
1716   printf("Isolated Trigger?  %d\n"     , fSelectIsolated) ;
1717   printf("Several UE?  %d\n"           , fMakeSeveralUE) ;
1718   printf("Name of AOD Pi0 Branch %s \n", fPi0AODBranchName.Data());
1719   printf("Do Decay-hadron correlation ?  %d\n", fPi0Trigger) ;
1720   printf("Select absolute leading for cluster triggers ?  %d\n", fMakeAbsoluteLeading) ;
1721   printf("Trigger pt bins  %d\n", fNAssocPtBins) ;
1722   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
1723     printf("\t bin %d = [%2.1f,%2.1f]\n", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
1724   }
1725
1726 }
1727
1728 //____________________________________________________________
1729 void AliAnaParticleHadronCorrelation::SetNAssocPtBins(Int_t n)
1730 {
1731   // Set number of bins
1732
1733   fNAssocPtBins  = n ;
1734
1735
1736   if(n < 10 && n > 0)
1737   {
1738     fNAssocPtBins  = n ;
1739   }
1740   else
1741   {
1742     printf("n = larger than 9 or too small, set to 9 \n");
1743     fNAssocPtBins = 9;
1744   }
1745 }
1746
1747 //______________________________________________________________________________
1748 void AliAnaParticleHadronCorrelation::SetAssocPtBinLimit(Int_t ibin, Float_t pt)
1749 {
1750   // Set the list of limits for the trigger pt bins
1751
1752   if(ibin <= fNAssocPtBins || ibin >= 0)
1753   {
1754     fAssocPtBinLimit[ibin] = pt  ;
1755   }
1756   else {
1757     printf("AliAnaParticleHadronCorrelation::SetAssocPtBinLimit() - bin  number too large %d > %d or small, nothing done\n", ibin, fNAssocPtBins) ;
1758
1759   }
1760 }
1761