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