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 <TH2F.h>
  38 #include <TDatabasePDG.h>
  39
  40 //---- ANALYSIS system ----
  41 #include "AliNeutralMesonSelection.h"
  42 #include "AliAnaParticleHadronCorrelation.h"
  43 #include "AliCaloTrackReader.h"
  44 #include "AliAODPWG4ParticleCorrelation.h"
  45 #include "AliFiducialCut.h"
  46 #include "AliVTrack.h"
  47 #include "AliVCluster.h"
  48 #include "AliMCAnalysisUtils.h"
  49 #include "TParticle.h"
  50 #include "AliStack.h"
  51 #include "AliAODMCParticle.h"
  52 #include "AliMixedEvent.h"
  53 #include "AliAnalysisManager.h"
  54 #include "AliInputEventHandler.h"
  55 #include "AliEventplane.h"
  56
  57 ClassImp(AliAnaParticleHadronCorrelation)
  58
  59
  60 //___________________________________________________________________
  61   AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation():
  62     AliAnaCaloTrackCorrBaseClass(),
  63     fMinTriggerPt(0.),
  64     fMaxAssocPt(1000.),             fMinAssocPt(0.),
  65     fDeltaPhiMaxCut(0.),            fDeltaPhiMinCut(0.),
  66     fSelectIsolated(0),             fMakeSeveralUE(0),
  67     fUeDeltaPhiMaxCut(0.),          fUeDeltaPhiMinCut(0.),
  68     fPi0AODBranchName(""),          fNeutralCorr(0),
  69     fPi0Trigger(0),                 fDecayTrigger(0),
  70     fMakeAbsoluteLeading(0),        fMakeNearSideLeading(0),
  71     fLeadingTriggerIndex(-1),       fHMPIDCorrelation(0),  fFillBradHisto(0),
  72     fNAssocPtBins(0),               fAssocPtBinLimit(),
  73     fListMixTrackEvents(),          fListMixCaloEvents(),  fUseMixStoredInReader(0),
  74     fM02MaxCut(0),                  fM02MinCut(0),
  75     //Histograms
  76     fhPtLeading(0),                 fhPtLeadingBin(0),
  77     fhPhiLeading(0),                fhEtaLeading(0),
  78     fhPtLeadingCentrality(0),       fhPtLeadingEventPlane(0),
  79     fhLeadingEventPlaneCentrality(0),
  80     fhPtLeadingMixed(0),            fhPtLeadingMixedBin(0),
  81     fhPhiLeadingMixed(0),           fhEtaLeadingMixed(0),
  82     fhDeltaPhiDeltaEtaCharged(0),
  83     fhPhiCharged(0),                fhEtaCharged(0),
  84     fhDeltaPhiCharged(0),           fhDeltaEtaCharged(0),
  85     fhDeltaPhiChargedPt(0),         fhDeltaPhiUeChargedPt(0),
  86     fhUePart(0),
  87     fhXECharged(0),                 fhXEUeCharged(0),
  88     fhXEPosCharged(0),              fhXENegCharged(0),
  89     fhPtHbpXECharged(0),            fhPtHbpXEUeCharged(0),
  90     fhZTCharged(0),                 fhZTUeCharged(0),
  91     fhZTPosCharged(0),              fhZTNegCharged(0),
  92     fhPtHbpZTCharged(0),            fhPtHbpZTUeCharged(0),
  93     fhDeltaPhiUeLeftCharged(0),     fhDeltaPhiUeRightCharged(0),
  94     fhXEUeLeftCharged(0),           fhXEUeRightCharged(0),
  95     fhPtHbpXEUeLeftCharged(0),      fhPtHbpXEUeRightCharged(0),
  96     fhZTUeLeftCharged(0),           fhZTUeRightCharged(0),
  97     fhPtHbpZTUeLeftCharged(0),      fhPtHbpZTUeRightCharged(0),
  98     fhPtTrigPout(0),                fhPtTrigCharged(0),
  99     fhTrigDeltaPhiCharged(0x0),     fhTrigDeltaEtaCharged(0x0),
 100     fhTrigXECorr(0x0),              fhTrigXEUeCorr(0x0),
 101     fhTrigZTCorr(0x0),              fhTrigZTUeCorr(0x0),
 102     fhAssocPtBkg(0),
 103     fhDeltaPhiAssocPtBin(0),        fhDeltaPhiAssocPtBinDEta08(0),
 104     fhDeltaPhiAssocPtBinHMPID(0),   fhDeltaPhiAssocPtBinHMPIDAcc(0),
 105     fhDeltaPhiBradAssocPtBin(0),    fhDeltaPhiBrad(0),
 106     fhXEAssocPtBin(0),              fhZTAssocPtBin(0),
 107     fhDeltaPhiDeltaEtaNeutral(0),
 108     fhPhiNeutral(0),                fhEtaNeutral(0),
 109     fhDeltaPhiNeutral(0),           fhDeltaEtaNeutral(0),
 110     fhDeltaPhiNeutralPt(0),         fhDeltaPhiUeNeutralPt(0),
 111     fhXENeutral(0),                 fhXEUeNeutral(0),
 112     fhPtHbpXENeutral(0),            fhPtHbpXEUeNeutral(0),
 113     fhZTNeutral(0),                 fhZTUeNeutral(0),
 114     fhPtHbpZTNeutral(0),            fhPtHbpZTUeNeutral(0),
 115     fhDeltaPhiUeLeftNeutral(0),     fhDeltaPhiUeRightNeutral(0),
 116     fhXEUeLeftNeutral(0),           fhXEUeRightNeutral(0),
 117     fhPtHbpXEUeLeftNeutral(0),      fhPtHbpXEUeRightNeutral(0),
 118     fhZTUeLeftNeutral(0),           fhZTUeRightNeutral(0),
 119     fhPtHbpZTUeLeftNeutral(0),      fhPtHbpZTUeRightNeutral(0),
 120     fhPtPi0DecayRatio(0),
 121     fhDeltaPhiDecayCharged(0),      fhXEDecayCharged(0), fhZTDecayCharged(0),
 122     fhDeltaPhiDecayNeutral(0),      fhXEDecayNeutral(0), fhZTDecayNeutral(0),
 123     fhDeltaPhiDecayChargedAssocPtBin(0),
 124     fhXEDecayChargedAssocPtBin(0),  fhZTDecayChargedAssocPtBin(0),
 125     fh2phiLeadingParticle(0x0),     fhMCPtLeading(0),
 126     fhMCPhiLeading(0),              fhMCEtaLeading(0),
 127     fhMCEtaCharged(0),              fhMCPhiCharged(0),
 128     fhMCDeltaEtaCharged(0),         fhMCDeltaPhiCharged(0x0),
 129     fhMCDeltaPhiDeltaEtaCharged(0), fhMCDeltaPhiChargedPt(0),
 130     fhMCPtXECharged(0),             fhMCPtXEUeCharged(0),
 131     fhMCPtHbpXECharged(0),          fhMCPtHbpXEUeCharged(0),
 132     fhMCUePart(0),
 133     fhMCPtZTCharged(0),             fhMCPtHbpZTCharged(0),
 134     fhMCPtTrigPout(0),              fhMCPtAssocDeltaPhi(0),
 135     //Mixing
 136     fhNEventsTrigger(0),
 137     fhNtracksAll(0),                fhNtracksTrigger(0),
 138     fhNtracksMB(0),
 139     fhMixDeltaPhiCharged(0),        fhMixDeltaPhiDeltaEtaCharged(0),
 140     fhMixXECharged(0), fhMixHbpXECharged(0),
 141     fhMixDeltaPhiChargedAssocPtBin(),
 142     fhMixDeltaPhiChargedAssocPtBinDEta08(),
 143     fhMixDeltaPhiDeltaEtaChargedAssocPtBin(),
 144     fhEventBin(0),                  fhEventMixBin(0)
 145 {
 146   //Default Ctor
 147
 148   //Initialize parameters
 149   InitParameters();
 150 }
 151
 152 //_________________________________________________________________
 153 AliAnaParticleHadronCorrelation::~AliAnaParticleHadronCorrelation()
 154 {
 155   // Remove event containers
 156
 157   if(DoOwnMix())
 158   {
 159     if(fListMixTrackEvents)
 160     {
 161       for(Int_t iz=0; iz < GetNZvertBin(); iz++)
 162       {
 163         for(Int_t ic=0; ic < GetNCentrBin(); ic++)
 164         {
 165           for(Int_t irp=0; irp<GetNRPBin(); irp++)
 166           {
 167             Int_t bin = GetEventMixBin(ic, iz, irp);
 168             fListMixTrackEvents[bin]->Delete() ;
 169             delete fListMixTrackEvents[bin] ;
 170           }
 171         }
 172       }
 173     }
 174
 175     delete[] fListMixTrackEvents;
 176
 177     if(fListMixCaloEvents)
 178     {
 179       for(Int_t iz=0; iz < GetNZvertBin(); iz++)
 180       {
 181         for(Int_t ic=0; ic < GetNCentrBin(); ic++)
 182         {
 183           for(Int_t irp=0; irp<GetNRPBin(); irp++)
 184           {
 185             Int_t bin = GetEventMixBin(ic, iz, irp);
 186             fListMixCaloEvents[bin]->Delete() ;
 187             delete fListMixCaloEvents[bin] ;
 188           }
 189         }
 190       }
 191     }
 192
 193     delete[] fListMixCaloEvents;
 194
 195   }
 196 }
 197
 198 //______________________________________________________________________________________________________________________________________________________
 199 void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(const Float_t ptAssoc,  const Float_t ptTrig,   const Int_t   assocBin,
 200                                                                               const Float_t phiAssoc, const Float_t phiTrig,  Float_t &     deltaPhi,
 201                                                                               const Float_t etaAssoc, const Float_t etaTrig,
 202                                                                               const Bool_t  decay,    const Float_t hmpidSignal, const Int_t nTracks)
 203 {
 204   // Fill angular correlation related histograms
 205
 206   Float_t deltaEta    = etaTrig-etaAssoc;
 207           deltaPhi    = phiTrig-phiAssoc;
 208   Float_t deltaPhiOrg = deltaPhi;
 209
 210   if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
 211   if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
 212
 213   fhEtaCharged     ->Fill(ptAssoc,etaAssoc);
 214   fhPhiCharged     ->Fill(ptAssoc,phiAssoc);
 215   fhDeltaEtaCharged->Fill(ptTrig ,deltaEta);
 216   fhDeltaPhiCharged->Fill(ptTrig ,deltaPhi);
 217
 218   if(ptAssoc > 2 )           fhDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
 219
 220   if(fDecayTrigger && decay) fhDeltaPhiDecayCharged   ->Fill(ptTrig  , deltaPhi);
 221
 222   Double_t  dphiBrad = -100;
 223   if(fFillBradHisto)
 224   {
 225     dphiBrad = atan2(sin(deltaPhiOrg), cos(deltaPhiOrg))/TMath::Pi();//-1 to 1
 226     if(TMath::Abs(dphiBrad)>0.325 && TMath::Abs(dphiBrad)<0.475)  //Hardcoded values, BAD, FIXME
 227     {
 228       fhAssocPtBkg->Fill(ptTrig, ptAssoc);
 229     }
 230
 231     if(dphiBrad<-1./3) dphiBrad += 2;
 232     fhDeltaPhiBrad->Fill(ptTrig, dphiBrad);
 233   }
 234
 235   // Fill histograms in bins of associated particle pT
 236   if(assocBin>=0)
 237   {
 238       fhDeltaPhiAssocPtBin            [assocBin]->Fill(ptTrig, deltaPhi);
 239
 240     if(TMath::Abs(deltaEta)> 0.8)
 241       fhDeltaPhiAssocPtBinDEta08      [assocBin]->Fill(ptTrig, deltaPhi);
 242
 243     if (fFillBradHisto)
 244       fhDeltaPhiBradAssocPtBin        [assocBin]->Fill(ptTrig, dphiBrad);
 245
 246     if(fDecayTrigger && decay)
 247       fhDeltaPhiDecayChargedAssocPtBin[assocBin]->Fill(ptTrig, deltaPhi);
 248
 249     if(fHMPIDCorrelation)
 250     {
 251       if( hmpidSignal > 0 )
 252       {
 253         //printf("Track pt %f with HMPID signal %f \n",pt,hmpidSignal);
 254         fhDeltaPhiAssocPtBinHMPID[assocBin]->Fill(ptTrig, deltaPhi);
 255       }
 256
 257       if(phiAssoc > 5*TMath::DegToRad() && phiAssoc < 20*TMath::DegToRad())
 258       {
 259         //printf("Track pt %f in HMPID acceptance phi %f \n ",pt,phi*TMath::RadToDeg() );
 260         fhDeltaPhiAssocPtBinHMPIDAcc[assocBin]->Fill(ptTrig, deltaPhi);
 261       }
 262     }
 263   }
 264
 265   //fill different multiplicity histogram
 266   if(DoEventSelect())
 267   {
 268     for(Int_t im = 0; im<GetMultiBin(); im++)
 269     {
 270       if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
 271       {
 272         fhTrigDeltaPhiCharged[im]->Fill(ptTrig,deltaPhi);
 273         fhTrigDeltaEtaCharged[im]->Fill(ptTrig,deltaEta);
 274       }
 275     }
 276   }
 277 }
 278
 279 //____________________________________________________________________________________________________________________________________________________
 280 Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(const Float_t mcAssocPt,      Float_t mcAssocPhi, const Float_t mcAssocEta,
 281                                                                            const Float_t mcTrigPt, const Float_t mcTrigPhi,  const Float_t mcTrigEta)
 282 {
 283   // Fill MC histograms independently of AOD or ESD
 284
 285   //Select only hadrons in pt range
 286   if(mcAssocPt < fMinAssocPt || mcAssocPt > fMaxAssocPt) return kTRUE ; // exclude but continue
 287
 288   if(mcAssocPhi < 0) mcAssocPhi+=TMath::TwoPi();
 289
 290   //remove trigger itself for correlation when use charged triggers
 291   if(TMath::Abs(mcAssocPt -mcTrigPt ) < 1e-6 &&
 292      TMath::Abs(mcAssocPhi-mcTrigPhi) < 1e-6 &&
 293      TMath::Abs(mcAssocEta-mcTrigEta) < 1e-6)            return kTRUE ; // exclude but continue
 294
 295   // Absolute leading?
 296   if( fMakeAbsoluteLeading && mcAssocPt > mcTrigPt )     return kFALSE; // jump event
 297
 298   //jump out this event if near side associated partile pt larger than trigger
 299   if( fMakeNearSideLeading && mcAssocPt > mcTrigPt &&
 300      TMath::Abs(mcAssocPhi-mcTrigPhi)<TMath::PiOver2() ) return kFALSE; // jump event
 301
 302   Float_t mcdeltaPhi= mcTrigPhi-mcAssocPhi;
 303   if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
 304   if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
 305
 306   Float_t mcxE    =-mcAssocPt/mcTrigPt*TMath::Cos(mcdeltaPhi);// -(mcAssocPx*pxprim+mcAssocPy*pyprim)/(mcTrigPt*mcTrigPt);
 307   Float_t mchbpXE =-100 ;
 308   if(mcxE > 0 ) mchbpXE = TMath::Log(1./mcxE);
 309
 310   Float_t mczT    = mcAssocPt/mcTrigPt ;
 311   Float_t mchbpZT =-100 ;
 312   if(mczT > 0 ) mchbpZT = TMath::Log(1./mczT);
 313
 314   //Selection within angular range
 315   if( mcdeltaPhi< -TMath::PiOver2())  mcdeltaPhi+=TMath::TwoPi();
 316   if( mcdeltaPhi>3*TMath::PiOver2())  mcdeltaPhi-=TMath::TwoPi();
 317
 318   Double_t mcpout = mcAssocPt*TMath::Sin(mcdeltaPhi) ;
 319
 320   if(GetDebug() > 0 )
 321     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",
 322            mcAssocPt,mcAssocPhi, mcTrigPhi,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut, GetMinPt());
 323
 324   // Fill Histograms
 325   fhMCEtaCharged     ->Fill(mcAssocPt, mcAssocEta);
 326   fhMCPhiCharged     ->Fill(mcAssocPt, mcAssocPhi);
 327   fhMCDeltaEtaCharged->Fill(mcTrigPt,    mcTrigEta-mcAssocEta);
 328   fhMCDeltaPhiCharged->Fill(mcTrigPt,    mcdeltaPhi);
 329   fhMCPtAssocDeltaPhi->Fill(mcAssocPt, mcdeltaPhi);
 330
 331   fhMCDeltaPhiDeltaEtaCharged->Fill(mcdeltaPhi,mcTrigEta-mcAssocEta);
 332
 333   //delta phi cut for correlation
 334   if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) )
 335   {
 336     fhMCDeltaPhiChargedPt->Fill(mcAssocPt,mcdeltaPhi);
 337     fhMCPtXECharged      ->Fill(mcTrigPt,mcxE);
 338     fhMCPtHbpXECharged   ->Fill(mcTrigPt,mchbpXE);
 339     fhMCPtZTCharged      ->Fill(mcTrigPt,mczT);
 340     fhMCPtHbpZTCharged   ->Fill(mcTrigPt,mchbpZT);
 341     fhMCPtTrigPout       ->Fill(mcTrigPt, mcpout) ;
 342   }
 343
 344   //underlying event
 345   if ( (mcdeltaPhi > fUeDeltaPhiMinCut) && (mcdeltaPhi < fUeDeltaPhiMaxCut) )
 346     {
 347        Double_t randomphi = gRandom->Uniform(TMath::Pi()/2,3*TMath::Pi()/2);
 348        Double_t mcUexE = -(mcAssocPt/mcTrigPt)*TMath::Cos(randomphi);
 349
 350        if(mcUexE < 0.) mcUexE = -mcUexE;
 351
 352        fhMCPtXEUeCharged->Fill(mcTrigPt,mcUexE);
 353        if(mcUexE > 0) fhMCPtHbpXEUeCharged->Fill(mcTrigPt,TMath::Log(1/mcUexE));
 354
 355        fhMCUePart->Fill(mcTrigPt);
 356     }
 357
 358   return kTRUE;
 359 }
 360
 361 //___________________________________________________________________________________________________________________
 362 void AliAnaParticleHadronCorrelation::FillChargedMomentumImbalanceHistograms(const Float_t ptTrig,   const Float_t ptAssoc,
 363                                                                              const Float_t xE,       const Float_t hbpXE,
 364                                                                              const Float_t zT,       const Float_t hbpZT,
 365                                                                              const Float_t pout,     const Float_t deltaPhi,
 366                                                                              const Int_t   nTracks,  const Int_t   charge,
 367                                                                              const Int_t   assocBin, const Bool_t  decay    )
 368
 369 {
 370   // Fill mostly momentum imbalance related histograms
 371
 372   fhDeltaPhiChargedPt ->Fill(ptAssoc, deltaPhi);
 373   fhXECharged         ->Fill(ptTrig , xE);
 374   fhPtHbpXECharged    ->Fill(ptTrig , hbpXE);
 375   fhZTCharged         ->Fill(ptTrig , zT);
 376   fhPtHbpZTCharged    ->Fill(ptTrig , hbpZT);
 377   fhPtTrigPout        ->Fill(ptTrig , pout) ;
 378   fhPtTrigCharged     ->Fill(ptTrig , ptAssoc) ;
 379
 380   if(fDecayTrigger && decay)
 381   {
 382     fhXEDecayCharged->Fill(ptTrig,xE);
 383     fhZTDecayCharged->Fill(ptTrig,zT);
 384   } // photon decay pi0/eta trigger
 385
 386   if(assocBin >= 0 )//away side
 387   {
 388     fhXEAssocPtBin[assocBin]->Fill(ptTrig, xE) ;
 389     fhZTAssocPtBin[assocBin]->Fill(ptTrig, zT) ;
 390
 391     if(fDecayTrigger && decay)
 392     {
 393       fhXEDecayChargedAssocPtBin[assocBin]->Fill(ptTrig, xE);
 394       fhZTDecayChargedAssocPtBin[assocBin]->Fill(ptTrig, zT);
 395     }
 396   }
 397
 398   if(charge > 0)
 399   {
 400     fhXEPosCharged->Fill(ptTrig,xE) ;
 401     fhZTPosCharged->Fill(ptTrig,zT) ;
 402   }
 403   else
 404   {
 405     fhXENegCharged->Fill(ptTrig,xE) ;
 406     fhZTNegCharged->Fill(ptTrig,zT) ;
 407   }
 408
 409   //fill different multiplicity histogram
 410   if(DoEventSelect())
 411   {
 412     for(Int_t im=0; im<GetMultiBin(); im++)
 413     {
 414       if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
 415       {
 416         fhTrigXECorr[im]->Fill(ptTrig,xE);
 417         fhTrigZTCorr[im]->Fill(ptTrig,zT);
 418       }
 419     }
 420   } //multiplicity events selection
 421 }
 422
 423 //_______________________________________________________________________________________________________________________
 424 void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventHistograms(const Float_t ptTrig,   const Float_t ptAssoc,
 425                                                                            const Float_t deltaPhi, const Int_t nTracks)
 426 {
 427   // Fill underlying event histograms
 428
 429   fhDeltaPhiUeChargedPt->Fill(ptAssoc,deltaPhi);
 430
 431   Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 432   Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 433   Double_t uezT =   ptAssoc/ptTrig;
 434
 435   if(uexE < 0.) uexE = -uexE;
 436
 437   fhXEUeCharged->Fill(ptTrig,uexE);
 438   if(uexE > 0) fhPtHbpXEUeCharged->Fill(ptTrig,TMath::Log(1/uexE));
 439
 440   fhZTUeCharged->Fill(ptTrig,uezT);
 441   if(uexE > 0) fhPtHbpZTUeCharged->Fill(ptTrig,TMath::Log(1/uezT));
 442
 443   if(DoEventSelect())
 444   {
 445     for(Int_t im=0; im<GetMultiBin(); im++)
 446     {
 447       if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))
 448       {
 449         fhTrigXEUeCorr[im]->Fill(ptTrig,uexE); // xE? CHECK
 450         fhTrigZTUeCorr[im]->Fill(ptTrig,uezT); // zT? CHECK
 451       }
 452     }
 453   } //multiplicity events selection
 454 }
 455
 456 //___________________________________________________________________________________________________________________________
 457 void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventSidesHistograms(const Float_t ptTrig,   const Float_t ptAssoc,
 458                                                                                 const Float_t xE,       const Float_t hbpXE,
 459                                                                                 const Float_t zT,       const Float_t hbpZT,
 460                                                                                 const Float_t deltaPhi)
 461 {
 462  // Fill underlying event histograms to the left and right of trigger
 463   if((deltaPhi<-fUeDeltaPhiMinCut) || (deltaPhi >2*fUeDeltaPhiMaxCut))
 464   {
 465     fhDeltaPhiUeLeftCharged->Fill(ptAssoc,deltaPhi);
 466     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 467     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 468     Double_t uezT =   ptAssoc/ptTrig;
 469
 470     if(uexE < 0.) uexE = -uexE;
 471
 472     fhXEUeLeftCharged->Fill(ptTrig,uexE);
 473     if(uexE > 0) fhPtHbpXEUeLeftCharged->Fill(ptTrig,TMath::Log(1/uexE));
 474
 475     fhZTUeLeftCharged->Fill(ptTrig,uezT);
 476     if(uexE > 0) fhPtHbpZTUeLeftCharged->Fill(ptTrig,TMath::Log(1/uezT));
 477     fhDeltaPhiUeLeftCharged->Fill(ptAssoc, deltaPhi);
 478   }
 479
 480   if((deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut))
 481   {
 482     fhDeltaPhiUeRightCharged->Fill(ptAssoc,deltaPhi);
 483     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 484     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 485     Double_t uezT =   ptAssoc/ptTrig;
 486
 487     if(uexE < 0.) uexE = -uexE;
 488
 489     fhXEUeRightCharged->Fill(ptTrig,uexE);
 490     if(uexE > 0) fhPtHbpXEUeRightCharged->Fill(ptTrig,TMath::Log(1/uexE));
 491
 492     fhZTUeRightCharged->Fill(ptTrig,uezT);
 493     if(uexE > 0) fhPtHbpZTUeRightCharged->Fill(ptTrig,TMath::Log(1/uezT));
 494     fhDeltaPhiUeRightCharged->Fill(ptAssoc, deltaPhi);
 495   }
 496 }
 497
 498 //______________________________________________________________________________________________________________________________
 499 void AliAnaParticleHadronCorrelation::FillDecayPhotonCorrelationHistograms(const Float_t ptAssoc,     const Float_t phiAssoc,
 500                                                                            const TLorentzVector mom1, const TLorentzVector mom2,
 501                                                                            const Bool_t bChargedOrNeutral)
 502 {
 503   // Do correlation with decay photons of triggered pi0 or eta
 504
 505   // Calculate the correlation parameters
 506   Float_t ptDecay1 = mom1.Pt();
 507   Float_t ptDecay2 = mom2.Pt();
 508
 509   Float_t zTDecay1 = -100, zTDecay2 = -100;
 510   if(ptDecay1) zTDecay1 = ptAssoc/ptDecay1 ;
 511   if(ptDecay2) zTDecay2 = ptAssoc/ptDecay2 ;
 512
 513   Float_t deltaPhiDecay1 = mom1.Phi()-phiAssoc;
 514   if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
 515   if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
 516
 517   Float_t deltaPhiDecay2 = mom2.Phi()-phiAssoc;
 518   if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
 519   if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
 520
 521   Float_t xEDecay1   =-zTDecay1*TMath::Cos(deltaPhiDecay1); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 522   Float_t xEDecay2   =-zTDecay2*TMath::Cos(deltaPhiDecay2); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 523
 524   if(bChargedOrNeutral) // correlate with charges
 525   {
 526     fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaPhiDecay1);
 527     fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaPhiDecay2);
 528
 529     if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::FillDecayPhotonHistograms( Charged corr) - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
 530
 531     if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
 532     {
 533       fhZTDecayCharged->Fill(ptDecay1,zTDecay1);
 534       fhXEDecayCharged->Fill(ptDecay1,xEDecay1);
 535     }
 536     if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
 537     {
 538       fhZTDecayCharged->Fill(ptDecay2,zTDecay2);
 539       fhXEDecayCharged->Fill(ptDecay2,xEDecay2);
 540     }
 541   }
 542   else // correlate with neutrals
 543   {
 544     fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaPhiDecay1);
 545     fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaPhiDecay2);
 546
 547     if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::FillDecayPhotonHistograms(Neutral corr) - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
 548
 549     if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
 550     {
 551       fhZTDecayCharged->Fill(ptDecay1,zTDecay1);
 552       fhXEDecayCharged->Fill(ptDecay1,xEDecay1);
 553     }
 554     if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
 555     {
 556       fhZTDecayCharged->Fill(ptDecay2,zTDecay2);
 557       fhXEDecayCharged->Fill(ptDecay2,xEDecay2);
 558     }
 559   }
 560 }
 561
 562 //______________________________________________________________________________________________________________________________________________________
 563 void AliAnaParticleHadronCorrelation::FillNeutralAngularCorrelationHistograms(const Float_t ptAssoc,  const Float_t ptTrig,
 564                                                                               const Float_t phiAssoc, const Float_t phiTrig,  Float_t &     deltaPhi,
 565                                                                               const Float_t etaAssoc, const Float_t etaTrig)
 566 {
 567   // Fill angular correlation related histograms
 568
 569   Float_t deltaEta    = etaTrig-etaAssoc;
 570   deltaPhi    = phiTrig-phiAssoc;
 571
 572   if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
 573   if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
 574
 575   fhEtaNeutral     ->Fill(ptAssoc,etaAssoc);
 576   fhPhiNeutral     ->Fill(ptAssoc,phiAssoc);
 577   fhDeltaEtaNeutral->Fill(ptTrig ,deltaEta);
 578   fhDeltaPhiNeutral->Fill(ptTrig ,deltaPhi);
 579
 580   if(ptAssoc > 2 ) fhDeltaPhiDeltaEtaNeutral->Fill(deltaPhi, deltaEta);
 581
 582 }
 583
 584 //_____________________________________________________________________________________________________________________________
 585 void AliAnaParticleHadronCorrelation::FillNeutralUnderlyingEventSidesHistograms(const Float_t ptTrig,   const Float_t ptAssoc,
 586                                                                                 const Float_t xE,       const Float_t hbpXE,
 587                                                                                 const Float_t zT,       const Float_t hbpZT,
 588                                                                                 const Float_t deltaPhi)
 589 {
 590   // Fill underlying event histograms to the left and right of trigger
 591
 592   if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut))
 593   {
 594     fhDeltaPhiUeLeftNeutral->Fill(ptAssoc, deltaPhi);
 595     fhXEUeLeftNeutral      ->Fill(ptTrig , xE);
 596     fhPtHbpXEUeLeftNeutral ->Fill(ptTrig , hbpXE);
 597     fhZTUeLeftNeutral      ->Fill(ptTrig , zT);
 598     fhPtHbpZTUeLeftNeutral ->Fill(ptTrig , hbpZT);
 599   }
 600
 601   if((deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut))
 602   {
 603     fhDeltaPhiUeRightNeutral->Fill(ptAssoc, deltaPhi);
 604     fhXEUeRightNeutral      ->Fill(ptTrig , xE);
 605     fhPtHbpXEUeRightNeutral ->Fill(ptTrig , hbpXE);
 606     fhZTUeRightNeutral      ->Fill(ptTrig , zT);
 607     fhPtHbpZTUeRightNeutral ->Fill(ptTrig , hbpZT);
 608   }
 609 }
 610
 611 //_____________________________________________________________
 612 void AliAnaParticleHadronCorrelation::FillChargedEventMixPool()
 613 {
 614   // Mixed event pool filling for tracks
 615
 616   //printf("FillChargedEventMixPool for %s\n",GetInputAODName().Data());
 617
 618   if(fUseMixStoredInReader && GetReader()->GetLastTracksMixedEvent() == GetEventNumber())
 619   {
 620     //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
 621     return ; // pool filled previously for another trigger
 622   }
 623
 624   Int_t nTracks = GetCTSTracks()->GetEntriesFast();
 625
 626   fhNtracksAll->Fill(nTracks);
 627
 628   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
 629   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
 630
 631   if(!inputHandler) return ;
 632
 633   if( inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask()    )
 634   {
 635     fhNtracksTrigger->Fill(nTracks);
 636   }
 637
 638   // Do mixing only with MB event (or the chosen mask), if not skip
 639   if( !(inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask()) ) return ;
 640
 641   fhNtracksMB->Fill(nTracks);
 642
 643   Int_t eventBin = GetEventMixBin();
 644
 645   //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
 646   if(eventBin < 0) return;
 647
 648   TObjArray * mixEventTracks = new TObjArray;
 649
 650   if(fUseMixStoredInReader)
 651   {
 652     fListMixTrackEvents[eventBin] = GetReader()->GetListWithMixedEventsForTracks(eventBin);
 653   }
 654
 655   if(!fListMixTrackEvents[eventBin]) fListMixTrackEvents[eventBin] = new TList();
 656
 657   //printf("%s ***** Pool Event bin : %d - nTracks %d\n",GetInputAODName().Data(),eventBin, GetCTSTracks()->GetEntriesFast());
 658
 659   TList * pool = fListMixTrackEvents[eventBin];
 660
 661   TVector3 p3;
 662   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
 663   {
 664     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
 665
 666     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 667     p3.SetXYZ(mom[0],mom[1],mom[2]);
 668     Float_t pt   = p3.Pt();
 669
 670     //Select only hadrons in pt range
 671     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
 672
 673     AliAODPWG4Particle * mixedTrack = new AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
 674     mixedTrack->SetDetector("CTS");
 675     mixedTrack->SetChargedBit(track->Charge()>0);
 676     mixEventTracks->Add(mixedTrack);
 677   }
 678
 679   //Set the event number where the last event was added, to avoid double pool filling
 680   GetReader()->SetLastTracksMixedEvent(GetEventNumber());
 681
 682   //printf("Add event to pool with %d tracks \n ",mixEventTracks->GetEntries());
 683   pool->AddFirst(mixEventTracks);
 684   mixEventTracks = 0;
 685
 686   //printf("Pool size %d, max %d\n",pool->GetSize(), GetNMaxEvMix());
 687
 688   if(pool->GetSize() > GetNMaxEvMix())
 689   {//Remove last event
 690     TClonesArray * tmp = static_cast<TClonesArray*>(pool->Last()) ;
 691     pool->RemoveLast() ;
 692     delete tmp ;
 693   }
 694 }
 695
 696 //_____________________________________________________________
 697 void AliAnaParticleHadronCorrelation::FillNeutralEventMixPool()
 698 {
 699   // Mixed event pool filling for neutral clusters
 700   // Right now only for EMCAL and in isolation case
 701
 702   //printf("FillNeutralEventMixPool for %s\n",GetInputAODName().Data());
 703
 704   if(fUseMixStoredInReader && GetReader()->GetLastCaloMixedEvent() == GetEventNumber())
 705   {
 706     //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
 707     return ; // pool filled previously for another trigger
 708   }
 709
 710   //  Int_t nClusters = GetEMCALClusters()->GetEntriesFast();
 711   //
 712   //  fhNclustersAll->Fill(nClusters);
 713
 714   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
 715   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
 716
 717   if(!inputHandler) return ;
 718
 719   //  if( inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask()    )
 720   //  {
 721   //    fhNclustersTrigger->Fill(nClusters);
 722   //  }
 723
 724   // Do mixing only with MB event (or the chosen mask), if not skip
 725   if( !(inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask()) ) return ;
 726
 727   //  fhNClustersMB->Fill(nCluster);
 728
 729   Int_t eventBin = GetEventMixBin();
 730
 731   //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
 732   if(eventBin < 0) return;
 733
 734   TObjArray * mixEventCalo = new TObjArray;
 735
 736   if(fUseMixStoredInReader)
 737   {
 738     fListMixCaloEvents[eventBin] = GetReader()->GetListWithMixedEventsForCalo(eventBin);
 739   }
 740
 741   if(!fListMixCaloEvents[eventBin]) fListMixCaloEvents[eventBin] = new TList();
 742
 743   TList * poolCalo = fListMixCaloEvents[eventBin];
 744
 745   TObjArray * pl = GetEMCALClusters();
 746   //if (GetAODObjArrayName.Contains("PHOS") )pl    = GetPHOSClusters();
 747   //else                                     pl    = GetEMCALClusters();
 748
 749   TLorentzVector mom;
 750   //printf("NClusters before selection %d\n",pl->GetEntriesFast());
 751   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ )
 752   {
 753     AliVCluster * calo = (AliVCluster *) (pl->At(ipr)) ;
 754
 755     // remove matched clusters
 756     if( IsTrackMatched( calo, GetReader()->GetInputEvent() ) ) continue ;
 757
 758     //Cluster momentum calculation
 759     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 760     {
 761       calo->GetMomentum(mom,GetVertex(0)) ;
 762     }//Assume that come from vertex in straight line
 763     else
 764     {
 765       Double_t vertex[]={0,0,0};
 766       calo->GetMomentum(mom,vertex) ;
 767     }
 768
 769     Float_t pt = mom.Pt();
 770     //Select only clusters in pt range
 771     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
 772
 773     AliAODPWG4Particle * mixedCalo = new AliAODPWG4Particle(mom);
 774     mixedCalo->SetDetector("EMCAL");
 775     mixEventCalo->Add(mixedCalo);
 776   }
 777
 778   //Set the event number where the last event was added, to avoid double pool filling
 779   GetReader()->SetLastCaloMixedEvent(GetEventNumber());
 780
 781   //printf("Add event to pool with %d clusters \n ",mixEventCalo->GetEntries());
 782   poolCalo->AddFirst(mixEventCalo);
 783   mixEventCalo = 0;
 784
 785   //printf("Pool size %d, max %d\n",poolCalo->GetSize(), GetNMaxEvMix());
 786
 787   if(poolCalo->GetSize() > GetNMaxEvMix())
 788   {//Remove last event
 789     TClonesArray * tmp = static_cast<TClonesArray*>(poolCalo->Last()) ;
 790     poolCalo->RemoveLast() ;
 791     delete tmp ;
 792   }
 793 }
 794
 795 //____________________________________________________________
 796 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
 797 {
 798   //Save parameters used for analysis
 799   TString parList ; //this will be list of parameters used for this analysis.
 800   const Int_t buffersize = 560;
 801   char onePar[buffersize] ;
 802
 803   snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;
 804   parList+=onePar ;
 805   snprintf(onePar,buffersize," Pt Trigger > %3.2f ",    fMinTriggerPt) ;
 806   parList+=onePar ;
 807   snprintf(onePar,buffersize," %3.2f < Pt associated < %3.2f ", fMinAssocPt,   fMaxAssocPt) ;
 808   parList+=onePar ;
 809   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-Hadron < %3.2f ",    fDeltaPhiMinCut,   fDeltaPhiMaxCut) ;
 810   parList+=onePar ;
 811   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-UeHadron <  %3.2f ", fUeDeltaPhiMinCut, fUeDeltaPhiMaxCut) ;
 812   parList+=onePar ;
 813   snprintf(onePar,buffersize,"Isolated Trigger?  %d\n",    fSelectIsolated) ;
 814   parList+=onePar ;
 815   snprintf(onePar,buffersize,"Several UE?  %d\n",          fMakeSeveralUE) ;
 816   parList+=onePar ;
 817   snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ", fPi0AODBranchName.Data());
 818   parList+=onePar ;
 819   snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  pi0 %d, decay %d", fPi0Trigger, fDecayTrigger) ;
 820   parList+=onePar ;
 821   snprintf(onePar,buffersize,"Select absolute leading for cluster triggers ? %d or Near Side Leading %d \n",
 822            fMakeAbsoluteLeading, fMakeNearSideLeading) ;
 823   parList+=onePar ;
 824   snprintf(onePar,buffersize,"Associated particle pt bins  %d: ", fNAssocPtBins) ;
 825   parList+=onePar ;
 826   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
 827     snprintf(onePar,buffersize,"bin %d = [%2.1f,%2.1f];", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
 828   }
 829   parList+=onePar ;
 830
 831   //Get parameters set in base class.
 832   parList += GetBaseParametersList() ;
 833
 834   //Get parameters set in FiducialCut class (not available yet)
 835   //parlist += GetFidCut()->GetFidCutParametersList()
 836
 837   return new TObjString(parList) ;
 838
 839 }
 840
 841 //________________________________________________________________
 842 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
 843 {
 844
 845   // Create histograms to be saved in output file and
 846   // store them in fOutputContainer
 847
 848   TList * outputContainer = new TList() ;
 849   outputContainer->SetName("CorrelationHistos") ;
 850
 851   Int_t   nptbins = GetHistogramRanges()->GetHistoPtBins(); Int_t  nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins(); Int_t  ndeltaphibins = GetHistogramRanges()->GetHistoDeltaPhiBins(); Int_t   ndeltaetabins = GetHistogramRanges()->GetHistoDeltaEtaBins();
 852   Float_t ptmax   = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax  = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax(); Float_t deltaphimax  = GetHistogramRanges()->GetHistoDeltaPhiMax();  Float_t deltaetamax   = GetHistogramRanges()->GetHistoDeltaEtaMax();
 853   Float_t ptmin   = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin  = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin(); Float_t deltaphimin  = GetHistogramRanges()->GetHistoDeltaPhiMin();  Float_t deltaetamin   = GetHistogramRanges()->GetHistoDeltaEtaMin();
 854
 855   Int_t nMixBins = GetNCentrBin()*GetNZvertBin()*GetNRPBin();
 856
 857   fhPtLeading  = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax);
 858   fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
 859
 860   fhPtLeadingBin  = new TH2F ("hPtLeadingBin","p_T distribution of leading particles", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
 861   fhPtLeadingBin->SetXTitle("p_{T}^{trig} (GeV/c)");
 862   fhPtLeadingBin->SetYTitle("Bin");
 863
 864   fhPhiLeading  = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 865   fhPhiLeading->SetYTitle("#phi (rad)");
 866
 867   fhEtaLeading  = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 868   fhEtaLeading->SetYTitle("#eta ");
 869
 870   outputContainer->Add(fhPtLeading);
 871   outputContainer->Add(fhPtLeadingBin);
 872   outputContainer->Add(fhPhiLeading);
 873   outputContainer->Add(fhEtaLeading);
 874
 875   fhPtLeadingCentrality   = new TH2F("hPtLeadingCentrality","Leading particle p_{T} vs centrality",nptbins,ptmin,ptmax,100,0.,100) ;
 876   fhPtLeadingCentrality->SetXTitle("p_{T}^{trig} (GeV/c)");
 877   fhPtLeadingCentrality->SetYTitle("Centrality (%)");
 878   outputContainer->Add(fhPtLeadingCentrality) ;
 879
 880   fhPtLeadingEventPlane  = new TH2F("hPtLeadingEventPlane","Leading particle p_{T} vs event plane angle",nptbins,ptmin,ptmax, 100,0.,TMath::Pi()) ;
 881   fhPtLeadingEventPlane->SetXTitle("p_{T}^{trig} (GeV/c)");
 882   fhPtLeadingEventPlane->SetXTitle("EP angle (rad)");
 883   outputContainer->Add(fhPtLeadingEventPlane) ;
 884
 885   fhLeadingEventPlaneCentrality  = new TH2F("hLeadingEventPlane","Leading particle centrality vs event plane angle",100,0.,100,100,0.,TMath::Pi()) ;
 886   fhLeadingEventPlaneCentrality->SetXTitle("Centrality (%)");
 887   fhLeadingEventPlaneCentrality->SetYTitle("EP angle (rad)");
 888   outputContainer->Add(fhLeadingEventPlaneCentrality) ;
 889
 890   //Correlation with charged hadrons
 891   if(GetReader()->IsCTSSwitchedOn())
 892   {
 893     fhDeltaPhiDeltaEtaCharged  = new TH2F
 894     ("hDeltaPhiDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}}",
 895     ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
 896     fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
 897     fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
 898
 899     fhPhiCharged  = new TH2F
 900     ("hPhiCharged","#phi_{h^{#pm}}  vs p_{T #pm}",
 901      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 902     fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
 903     fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
 904
 905     fhEtaCharged  = new TH2F
 906     ("hEtaCharged","#eta_{h^{#pm}}  vs p_{T #pm}",
 907      nptbins,ptmin,ptmax,netabins,etamin,etamax);
 908     fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
 909     fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
 910
 911     fhDeltaPhiCharged  = new TH2F
 912     ("hDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
 913      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
 914     fhDeltaPhiCharged->SetYTitle("#Delta #phi");
 915     fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
 916
 917     fhDeltaPhiChargedPt  = new TH2F
 918     ("hDeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
 919      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
 920     fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
 921     fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 922
 923     fhDeltaPhiUeChargedPt  = new TH2F
 924     ("hDeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
 925      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
 926     fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
 927     fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 928
 929     fhUePart  =  new TH1F("hUePart","UE particles distribution vs pt trig",
 930              nptbins,ptmin,ptmax);
 931     fhUePart->SetYTitle("dNch");
 932     fhUePart->SetXTitle("p_{T trigger}");
 933
 934
 935     fhDeltaEtaCharged  = new TH2F
 936     ("hDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
 937      nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
 938     fhDeltaEtaCharged->SetYTitle("#Delta #eta");
 939     fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
 940
 941     fhXECharged  =
 942     new TH2F("hXECharged","x_{E} for charged tracks",
 943              nptbins,ptmin,ptmax,200,0.,2.);
 944     fhXECharged->SetYTitle("x_{E}");
 945     fhXECharged->SetXTitle("p_{T trigger}");
 946
 947     fhXEUeCharged  =
 948     new TH2F("hXEUeCharged","x_{E} for Underlying Event",
 949              nptbins,ptmin,ptmax,200,0.,2.);
 950     fhXEUeCharged->SetYTitle("x_{E}");
 951     fhXEUeCharged->SetXTitle("p_{T trigger}");
 952
 953     fhXEPosCharged  =
 954     new TH2F("hXEPositiveCharged","x_{E} for positive charged tracks",
 955              nptbins,ptmin,ptmax,200,0.,2.);
 956     fhXEPosCharged->SetYTitle("x_{E}");
 957     fhXEPosCharged->SetXTitle("p_{T trigger}");
 958
 959     fhXENegCharged  =
 960     new TH2F("hXENegativeCharged","x_{E} for negative charged tracks",
 961              nptbins,ptmin,ptmax,200,0.,2.);
 962     fhXENegCharged->SetYTitle("x_{E}");
 963     fhXENegCharged->SetXTitle("p_{T trigger}");
 964
 965     fhPtHbpXECharged  =
 966     new TH2F("hHbpXECharged","#xi = ln(1/x_{E}) with charged hadrons",
 967              nptbins,ptmin,ptmax,200,0.,10.);
 968     fhPtHbpXECharged->SetYTitle("ln(1/x_{E})");
 969     fhPtHbpXECharged->SetXTitle("p_{T trigger}");
 970
 971     fhPtHbpXEUeCharged  =
 972     new TH2F("hHbpXEUeCharged","#xi = ln(1/x_{E}) with charged hadrons,Underlying Event",
 973              nptbins,ptmin,ptmax,200,0.,10.);
 974     fhPtHbpXEUeCharged->SetYTitle("ln(1/x_{E})");
 975     fhPtHbpXEUeCharged->SetXTitle("p_{T trigger}");
 976
 977     fhZTCharged  =
 978     new TH2F("hZTCharged","z_{T} for charged tracks",
 979              nptbins,ptmin,ptmax,200,0.,2.);
 980     fhZTCharged->SetYTitle("z_{T}");
 981     fhZTCharged->SetXTitle("p_{T trigger}");
 982
 983     fhZTUeCharged  =
 984     new TH2F("hZTUeCharged","z_{T} for Underlying Event",
 985              nptbins,ptmin,ptmax,200,0.,2.);
 986     fhZTUeCharged->SetYTitle("z_{T}");
 987     fhZTUeCharged->SetXTitle("p_{T trigger}");
 988
 989     fhZTPosCharged  =
 990     new TH2F("hZTPositiveCharged","z_{T} for positive charged tracks",
 991              nptbins,ptmin,ptmax,200,0.,2.);
 992     fhZTPosCharged->SetYTitle("z_{T}");
 993     fhZTPosCharged->SetXTitle("p_{T trigger}");
 994
 995     fhZTNegCharged  =
 996     new TH2F("hZTNegativeCharged","z_{T} for negative charged tracks",
 997              nptbins,ptmin,ptmax,200,0.,2.);
 998     fhZTNegCharged->SetYTitle("z_{T}");
 999     fhZTNegCharged->SetXTitle("p_{T trigger}");
1000
1001     fhPtHbpZTCharged  =
1002     new TH2F("hHbpZTCharged","#xi = ln(1/z_{T}) with charged hadrons",
1003              nptbins,ptmin,ptmax,200,0.,10.);
1004     fhPtHbpZTCharged->SetYTitle("ln(1/z_{T})");
1005     fhPtHbpZTCharged->SetXTitle("p_{T trigger}");
1006
1007     fhPtHbpZTUeCharged  =
1008     new TH2F("hHbpZTUeCharged","#xi = ln(1/z_{T}) with charged hadrons,Underlying Event",
1009              nptbins,ptmin,ptmax,200,0.,10.);
1010     fhPtHbpZTUeCharged->SetYTitle("ln(1/x_{E})");
1011     fhPtHbpZTUeCharged->SetXTitle("p_{T trigger}");
1012
1013     fhPtTrigPout  =
1014     new TH2F("hPtTrigPout","Pout with triggers",
1015              nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
1016     fhPtTrigPout->SetYTitle("p_{out} (GeV/c)");
1017     fhPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
1018
1019     fhPtTrigCharged  =
1020     new TH2F("hPtTrigCharged","trgger and charged tracks pt distribution",
1021              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1022     fhPtTrigCharged->SetYTitle("p_{T h^{#pm}} (GeV/c)");
1023     fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)");
1024
1025     outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
1026     outputContainer->Add(fhPhiCharged) ;
1027     outputContainer->Add(fhEtaCharged) ;
1028     outputContainer->Add(fhDeltaPhiCharged) ;
1029     outputContainer->Add(fhDeltaEtaCharged) ;
1030     outputContainer->Add(fhDeltaPhiChargedPt) ;
1031     outputContainer->Add(fhDeltaPhiUeChargedPt) ;
1032     outputContainer->Add(fhUePart);
1033
1034     outputContainer->Add(fhXECharged) ;
1035     outputContainer->Add(fhXEPosCharged) ;
1036     outputContainer->Add(fhXENegCharged) ;
1037     outputContainer->Add(fhXEUeCharged) ;
1038     outputContainer->Add(fhPtHbpXECharged) ;
1039     outputContainer->Add(fhPtHbpXEUeCharged) ;
1040
1041     outputContainer->Add(fhZTCharged) ;
1042     outputContainer->Add(fhZTPosCharged) ;
1043     outputContainer->Add(fhZTNegCharged) ;
1044     outputContainer->Add(fhZTUeCharged) ;
1045     outputContainer->Add(fhPtHbpZTCharged) ;
1046     outputContainer->Add(fhPtHbpZTUeCharged) ;
1047
1048     outputContainer->Add(fhPtTrigPout) ;
1049     outputContainer->Add(fhPtTrigCharged) ;
1050
1051     if(DoEventSelect())
1052     {
1053       Int_t nMultiBins = GetMultiBin();
1054       fhTrigDeltaPhiCharged = new TH2F*[nMultiBins] ;
1055       fhTrigDeltaEtaCharged = new TH2F*[nMultiBins] ;
1056       fhTrigXECorr          = new TH2F*[nMultiBins] ;
1057       fhTrigXEUeCorr        = new TH2F*[nMultiBins] ;
1058       fhTrigZTCorr          = new TH2F*[nMultiBins] ;
1059       fhTrigZTUeCorr        = new TH2F*[nMultiBins] ;
1060
1061       for(Int_t im=0; im<nMultiBins; im++)
1062       {
1063         fhTrigDeltaPhiCharged[im]  = new TH2F
1064         (Form("hTrigDeltaPhiCharged_%d",im),Form("hTrigDeltaPhiCharged_%d",im), nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1065         fhTrigDeltaPhiCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
1066         fhTrigDeltaPhiCharged[im]->SetYTitle("#Delta #phi");
1067
1068         fhTrigDeltaEtaCharged[im]  = new TH2F
1069         (Form("hTrigDeltaEtaCharged_%d",im),Form("hTrigDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
1070         fhTrigDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
1071         fhTrigDeltaEtaCharged[im]->SetYTitle("#Delta #eta");
1072
1073         fhTrigXECorr[im]  = new TH2F
1074         (Form("hTrigXEPtCorr_%d",im),Form("hTrigXEPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
1075         fhTrigXECorr[im]->SetYTitle("x_{E trigger h^{#pm}}");
1076         fhTrigXECorr[im]->SetXTitle("p_{T trigger}");
1077
1078         fhTrigXEUeCorr[im]  = new TH2F
1079         (Form("hTrigXEPtUeCorr_%d",im),Form("hTrigXEPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
1080         fhTrigXEUeCorr[im]->SetYTitle("x_{E trigger h^{#pm}}");
1081         fhTrigXEUeCorr[im]->SetXTitle("p_{T trigger}");
1082
1083         fhTrigZTCorr[im]  = new TH2F
1084         (Form("hTrigZTPtCorr_%d",im),Form("hTrigZTPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
1085         fhTrigZTCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
1086         fhTrigZTCorr[im]->SetXTitle("p_{T trigger}");
1087
1088         fhTrigZTUeCorr[im]  = new TH2F
1089         (Form("hTrigZTPtUeCorr_%d",im),Form("hTrigZTPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
1090         fhTrigZTUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
1091         fhTrigZTUeCorr[im]->SetXTitle("p_{T trigger}");
1092
1093         outputContainer->Add(fhTrigDeltaPhiCharged[im]) ;
1094         outputContainer->Add(fhTrigDeltaEtaCharged[im]) ;
1095         outputContainer->Add(fhTrigXECorr[im]);
1096         outputContainer->Add(fhTrigXEUeCorr[im]);
1097         outputContainer->Add(fhTrigZTCorr[im]);
1098         outputContainer->Add(fhTrigZTUeCorr[im]);
1099       }
1100     }
1101
1102     if(fFillBradHisto)
1103     {
1104       fhAssocPtBkg        = new TH2F("hAssocPtBkg", " Trigger p_{T} vs associated hadron p_{T} from background",
1105                                      nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
1106       fhAssocPtBkg->SetXTitle("p_{T trigger}");
1107       fhAssocPtBkg->SetYTitle("p_{T associated}");
1108       outputContainer->Add(fhAssocPtBkg) ;
1109
1110       fhDeltaPhiBrad = new TH2F("hDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs p_{T trigger} ",
1111                                 nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
1112       fhDeltaPhiBrad->SetXTitle("p_{T trigger}");
1113       fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
1114       outputContainer->Add(fhDeltaPhiBrad) ;
1115     }
1116
1117     fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins] ;
1118     fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins] ;
1119     fhXEAssocPtBin             = new TH2F*[fNAssocPtBins] ;
1120     fhZTAssocPtBin             = new TH2F*[fNAssocPtBins] ;
1121
1122     if(fFillBradHisto)
1123       fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins] ;
1124
1125     if(fPi0Trigger || fDecayTrigger)
1126     {
1127       fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins] ;
1128       fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins] ;
1129       fhXEAssocPtBin             = new TH2F*[fNAssocPtBins] ;
1130       fhZTAssocPtBin             = new TH2F*[fNAssocPtBins] ;
1131       fhXEDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
1132       fhZTDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
1133       fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
1134     }
1135
1136     if(fHMPIDCorrelation)
1137     {
1138       fhDeltaPhiAssocPtBinHMPID   = new TH2F*[fNAssocPtBins] ;
1139       fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins] ;
1140     }
1141
1142     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
1143     {
1144       fhDeltaPhiAssocPtBin[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1145                                          Form("#Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1146                                          nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1147       fhDeltaPhiAssocPtBin[i]->SetXTitle("p_{T trigger}");
1148       fhDeltaPhiAssocPtBin[i]->SetYTitle("#Delta #phi");
1149
1150       fhDeltaPhiAssocPtBinDEta08[i] = new TH2F(Form("hDeltaPhiDeltaEta0.8PtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1151                                                Form("#Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f], for #Delta #eta > 0.8", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1152                                                nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1153       fhDeltaPhiAssocPtBinDEta08[i]->SetXTitle("p_{T trigger}");
1154       fhDeltaPhiAssocPtBinDEta08[i]->SetYTitle("#Delta #phi");
1155
1156       fhXEAssocPtBin[i]       = new TH2F(Form("hXEAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1157                                          Form("x_{E} vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1158                                          nptbins, ptmin, ptmax,200, 0.0, 2.0);
1159       fhXEAssocPtBin[i]->SetXTitle("p_{T trigger}");
1160       fhXEAssocPtBin[i]->SetYTitle("x_{E}");
1161
1162       fhZTAssocPtBin[i]       = new TH2F(Form("hZTAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1163                                          Form("z_{T} vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1164                                          nptbins, ptmin, ptmax,200, 0.0, 2.0);
1165       fhZTAssocPtBin[i]->SetXTitle("p_{T trigger}");
1166       fhZTAssocPtBin[i]->SetYTitle("z_{T}");
1167
1168       outputContainer->Add(fhDeltaPhiAssocPtBin[i]) ;
1169       outputContainer->Add(fhDeltaPhiAssocPtBinDEta08[i]) ;
1170       outputContainer->Add(fhXEAssocPtBin[i]);
1171       outputContainer->Add(fhZTAssocPtBin[i]);
1172
1173       if(fPi0Trigger || fDecayTrigger)
1174       {
1175         fhDeltaPhiDecayChargedAssocPtBin[i] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1176                                            Form("#Delta #phi vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1177                                            nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1178         fhDeltaPhiDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
1179         fhDeltaPhiDecayChargedAssocPtBin[i]->SetYTitle("#Delta #phi");
1180
1181         fhXEDecayChargedAssocPtBin[i]       = new TH2F(Form("hXEDecayChargedAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1182                                            Form("x_{E} vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1183                                            nptbins, ptmin, ptmax,200, 0.0, 2.0);
1184         fhXEDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
1185         fhXEDecayChargedAssocPtBin[i]->SetYTitle("x_{E}");
1186
1187         fhZTDecayChargedAssocPtBin[i]       = new TH2F(Form("hZTDecayChargedAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1188                                            Form("z_{T} vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1189                                            nptbins, ptmin, ptmax,200, 0.0, 2.0);
1190         fhZTDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
1191         fhZTDecayChargedAssocPtBin[i]->SetYTitle("z_{T}");
1192
1193         outputContainer->Add(fhDeltaPhiDecayChargedAssocPtBin[i]) ;
1194         outputContainer->Add(fhXEDecayChargedAssocPtBin[i]);
1195         outputContainer->Add(fhZTDecayChargedAssocPtBin[i]);
1196
1197       }
1198
1199       if(fFillBradHisto)
1200       {
1201         fhDeltaPhiBradAssocPtBin[i] = new TH2F(Form("hDeltaPhiBradPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1202                                                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]),
1203                                                nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
1204         fhDeltaPhiBradAssocPtBin[i]->SetXTitle("p_{T trigger}");
1205         fhDeltaPhiBradAssocPtBin[i]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
1206         outputContainer->Add(fhDeltaPhiBradAssocPtBin[i]) ;
1207       }
1208
1209       if(fHMPIDCorrelation)
1210       {
1211         fhDeltaPhiAssocPtBinHMPID[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1fHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1212                                                 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]),
1213                                                 nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1214         fhDeltaPhiAssocPtBinHMPID[i]->SetXTitle("p_{T trigger}");
1215         fhDeltaPhiAssocPtBinHMPID[i]->SetYTitle("#Delta #phi");
1216
1217         fhDeltaPhiAssocPtBinHMPIDAcc[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1fHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1218                                                    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]),
1219                                                    nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1220         fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetXTitle("p_{T trigger}");
1221         fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetYTitle("#Delta #phi");
1222
1223         outputContainer->Add(fhDeltaPhiAssocPtBinHMPID[i]) ;
1224         outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[i]) ;
1225
1226       }
1227     }
1228
1229     if(fPi0Trigger || fDecayTrigger)
1230     {
1231       if(fPi0Trigger)
1232       {
1233         fhPtPi0DecayRatio  = new TH2F
1234         ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay",
1235          nptbins,ptmin,ptmax, 100,0.,2.);
1236         fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");
1237         fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");
1238         outputContainer->Add(fhPtPi0DecayRatio) ;
1239       }
1240
1241       fhDeltaPhiDecayCharged  = new TH2F
1242       ("hDeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",
1243        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1244       fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");
1245       fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");
1246
1247       fhXEDecayCharged  =
1248       new TH2F("hXEDecayCharged","x_{E}  Decay",
1249                nptbins,ptmin,ptmax,200,0.,2.);
1250       fhXEDecayCharged->SetYTitle("x_{E}");
1251       fhXEDecayCharged->SetXTitle("p_{T decay}");
1252
1253       fhZTDecayCharged  =
1254       new TH2F("hZTDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",
1255                nptbins,ptmin,ptmax,200,0.,2.);
1256       fhZTDecayCharged->SetYTitle("z_{decay h^{#pm}}");
1257       fhZTDecayCharged->SetXTitle("p_{T decay}");
1258
1259       outputContainer->Add(fhDeltaPhiDecayCharged) ;
1260       outputContainer->Add(fhXEDecayCharged) ;
1261       outputContainer->Add(fhZTDecayCharged) ;
1262     }
1263
1264     if(fMakeSeveralUE)
1265     {
1266       fhDeltaPhiUeLeftCharged  = new TH2F
1267       ("hDeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
1268        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1269       fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
1270       fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
1271       outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
1272
1273       fhDeltaPhiUeRightCharged  = new TH2F
1274       ("hDeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
1275        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1276       fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
1277       fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
1278       outputContainer->Add(fhDeltaPhiUeRightCharged) ;
1279
1280       fhXEUeLeftCharged  =
1281       new TH2F("hXEUeChargedLeft","x_{E} with UE left side of trigger",
1282                nptbins,ptmin,ptmax,200,0.,2.);
1283       fhXEUeLeftCharged->SetYTitle("x_{E Ueh^{#pm}}");
1284       fhXEUeLeftCharged->SetXTitle("p_{T trigger}");
1285       outputContainer->Add(fhXEUeLeftCharged) ;
1286
1287       fhXEUeRightCharged  =
1288       new TH2F("hXEUeChargedRight","x_{E h^{#pm}} with UE right side of trigger",
1289                nptbins,ptmin,ptmax,200,0.,2.);
1290       fhXEUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
1291       fhXEUeRightCharged->SetXTitle("p_{T trigger}");
1292       outputContainer->Add(fhXEUeRightCharged) ;
1293
1294       fhPtHbpXEUeLeftCharged  =
1295       new TH2F("hHbpXEUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",
1296                nptbins,ptmin,ptmax,200,0.,10.);
1297       fhPtHbpXEUeLeftCharged->SetYTitle("ln(1/x_{E})");
1298       fhPtHbpXEUeLeftCharged->SetXTitle("p_{T trigger}");
1299       outputContainer->Add(fhPtHbpXEUeLeftCharged) ;
1300
1301       fhPtHbpXEUeRightCharged  =
1302       new TH2F("hHbpXEUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",
1303                nptbins,ptmin,ptmax,200,0.,10.);
1304       fhPtHbpXEUeRightCharged->SetYTitle("ln(1/x_{E})");
1305       fhPtHbpXEUeRightCharged->SetXTitle("p_{T trigger}");
1306       outputContainer->Add(fhPtHbpXEUeRightCharged) ;
1307
1308       fhZTUeLeftCharged  =
1309       new TH2F("hZTUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",
1310                nptbins,ptmin,ptmax,200,0.,2.);
1311       fhZTUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
1312       fhZTUeLeftCharged->SetXTitle("p_{T trigger}");
1313       outputContainer->Add(fhZTUeLeftCharged) ;
1314
1315       fhZTUeRightCharged  =
1316       new TH2F("hZTUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",
1317                nptbins,ptmin,ptmax,200,0.,2.);
1318       fhZTUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
1319       fhZTUeRightCharged->SetXTitle("p_{T trigger}");
1320       outputContainer->Add(fhZTUeRightCharged) ;
1321
1322       fhPtHbpZTUeLeftCharged  =
1323       new TH2F("hHbpZTUeChargedLeft","#xi = ln(1/z_{T}) with charged UE left side of trigger",
1324                nptbins,ptmin,ptmax,200,0.,10.);
1325       fhPtHbpZTUeLeftCharged->SetYTitle("ln(1/z_{T})");
1326       fhPtHbpZTUeLeftCharged->SetXTitle("p_{T trigger}");
1327       outputContainer->Add(fhPtHbpZTUeLeftCharged) ;
1328
1329       fhPtHbpZTUeRightCharged  =
1330       new TH2F("hHbpZTUeChargedRight","#xi = ln(1/z_{T}) with charged UE right side of trigger",
1331                nptbins,ptmin,ptmax,200,0.,10.);
1332       fhPtHbpZTUeRightCharged->SetYTitle("ln(1/z_{T})");
1333       fhPtHbpZTUeRightCharged->SetXTitle("p_{T trigger}");
1334       outputContainer->Add(fhPtHbpZTUeRightCharged) ;
1335
1336     }
1337   }  //Correlation with charged hadrons
1338
1339   //Correlation with neutral hadrons
1340   if(fNeutralCorr)
1341   {
1342     fhDeltaPhiDeltaEtaNeutral  = new TH2F
1343     ("hDeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
1344      ndeltaphibins ,deltaphimin,deltaphimax, ndeltaetabins ,deltaetamin,deltaetamax);
1345     fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");
1346     fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
1347
1348     fhPhiNeutral  = new TH2F
1349     ("hPhiNeutral","#phi_{#pi^{0}}  vs p_{T #pi^{0}}",
1350      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1351     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
1352     fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
1353
1354     fhEtaNeutral  = new TH2F
1355     ("hEtaNeutral","#eta_{#pi^{0}}  vs p_{T #pi^{0}}",
1356      nptbins,ptmin,ptmax,netabins,etamin,etamax);
1357     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
1358     fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
1359
1360     fhDeltaPhiNeutral  = new TH2F
1361     ("hDeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
1362      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1363     fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
1364     fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
1365
1366     fhDeltaPhiNeutralPt  = new TH2F
1367     ("hDeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
1368      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1369     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
1370     fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
1371
1372     fhDeltaPhiUeNeutralPt  = new TH2F
1373     ("hDeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
1374      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1375     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
1376     fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
1377
1378     fhDeltaEtaNeutral  = new TH2F
1379     ("hDeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
1380      nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
1381     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
1382     fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
1383
1384     fhXENeutral  =
1385     new TH2F("hXENeutral","x_{E} for #pi^{0} associated",
1386              nptbins,ptmin,ptmax,200,0.,2.);
1387     fhXENeutral->SetYTitle("x_{E}");
1388     fhXENeutral->SetXTitle("p_{T trigger}");
1389
1390     fhXEUeNeutral  =
1391     new TH2F("hXEUeNeutral","x_{E} for #pi^{0} associated",
1392              nptbins,ptmin,ptmax,200,0.,2.);
1393     fhXEUeNeutral->SetYTitle("x_{E}");
1394     fhXEUeNeutral->SetXTitle("p_{T trigger}");
1395
1396     fhPtHbpXENeutral  =
1397     new TH2F("hHbpXENeutral","#xi = ln(1/x_{E})for #pi^{0} associated",
1398              nptbins,ptmin,ptmax,200,0.,10.);
1399     fhPtHbpXENeutral->SetYTitle("ln(1/x_{E})");
1400     fhPtHbpXENeutral->SetXTitle("p_{T trigger}");
1401
1402     fhPtHbpXEUeNeutral  =
1403     new TH2F("hHbpXEUeNeutral","#xi = ln(1/x_{E}) for #pi^{0} associated",
1404              nptbins,ptmin,ptmax,200,0.,10.);
1405     fhPtHbpXEUeNeutral->SetYTitle("ln(1/x_{E})");
1406     fhPtHbpXEUeNeutral->SetXTitle("p_{T trigger}");
1407
1408     fhZTNeutral  =
1409     new TH2F("hZTNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger} for #pi^{0} associated",
1410              nptbins,ptmin,ptmax,200,0.,2.);
1411     fhZTNeutral->SetYTitle("z_{trigger #pi^{0}}");
1412     fhZTNeutral->SetXTitle("p_{T trigger}");
1413
1414     fhZTUeNeutral  =
1415     new TH2F("hZTUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger} for #pi^{0} associated",
1416              nptbins,ptmin,ptmax,200,0.,2.);
1417     fhZTUeNeutral->SetYTitle("z_{trigger #pi^{0}}");
1418     fhZTUeNeutral->SetXTitle("p_{T trigger}");
1419
1420     fhPtHbpZTNeutral  =
1421     new TH2F("hHbpZTNeutral","#xi = ln(1/x_{E}) for #pi^{0} associated",
1422              nptbins,ptmin,ptmax,200,0.,10.);
1423     fhPtHbpZTNeutral->SetYTitle("ln(1/z_{T})");
1424     fhPtHbpZTNeutral->SetXTitle("p_{T trigger}");
1425
1426     fhPtHbpZTUeNeutral  =
1427     new TH2F("hHbpZTUeNeutral","#xi = ln(1/x_{E}) for #pi^{0} associated",
1428              nptbins,ptmin,ptmax,200,0.,10.);
1429     fhPtHbpXEUeNeutral->SetYTitle("ln(1/z_{T})");
1430     fhPtHbpXEUeNeutral->SetXTitle("p_{T trigger}");
1431
1432     outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
1433     outputContainer->Add(fhPhiNeutral) ;
1434     outputContainer->Add(fhEtaNeutral) ;
1435     outputContainer->Add(fhDeltaPhiNeutral) ;
1436     outputContainer->Add(fhDeltaPhiNeutralPt) ;
1437     outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
1438     outputContainer->Add(fhDeltaEtaNeutral) ;
1439     outputContainer->Add(fhXENeutral) ;
1440     outputContainer->Add(fhXEUeNeutral) ;
1441     outputContainer->Add(fhPtHbpXENeutral) ;
1442     outputContainer->Add(fhPtHbpXEUeNeutral) ;
1443     outputContainer->Add(fhZTNeutral) ;
1444     outputContainer->Add(fhZTUeNeutral) ;
1445     outputContainer->Add(fhPtHbpZTNeutral) ;
1446     outputContainer->Add(fhPtHbpZTUeNeutral) ;
1447
1448     if(fPi0Trigger || fDecayTrigger)
1449     {
1450       fhDeltaPhiDecayNeutral  = new TH2F
1451       ("hDeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",
1452        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1453       fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");
1454       fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");
1455
1456       fhXEDecayNeutral  =
1457       new TH2F("hXEDecayNeutral","x_{E} for decay trigger",
1458                nptbins,ptmin,ptmax,200,0.,2.);
1459       fhXEDecayNeutral->SetYTitle("x_{E}");
1460       fhXEDecayNeutral->SetXTitle("p_{T decay}");
1461
1462       fhZTDecayNeutral  =
1463       new TH2F("hZTDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",
1464                nptbins,ptmin,ptmax,200,0.,2.);
1465       fhZTDecayNeutral->SetYTitle("z_{h^{0}}");
1466       fhZTDecayNeutral->SetXTitle("p_{T decay}");
1467
1468       outputContainer->Add(fhDeltaPhiDecayNeutral) ;
1469       outputContainer->Add(fhXEDecayNeutral) ;
1470       outputContainer->Add(fhZTDecayNeutral) ;
1471
1472     }
1473
1474     if(fMakeSeveralUE)
1475     {
1476       fhDeltaPhiUeLeftNeutral  = new TH2F
1477       ("hDeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",
1478        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1479       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
1480       fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
1481       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
1482
1483       fhDeltaPhiUeRightNeutral  = new TH2F
1484       ("hDeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",
1485        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1486       fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
1487       fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
1488       outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
1489
1490       fhXEUeLeftNeutral  =
1491       new TH2F("hXEUeNeutralLeft","x_{E} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
1492                nptbins,ptmin,ptmax,140,0.,2.);
1493       fhXEUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1494       fhXEUeLeftNeutral->SetXTitle("p_{T trigger}");
1495       outputContainer->Add(fhXEUeLeftNeutral) ;
1496
1497       fhXEUeRightNeutral  =
1498       new TH2F("hXEUeNeutralRight","x_{E} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
1499                nptbins,ptmin,ptmax,200,0.,2.);
1500       fhXEUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1501       fhXEUeRightNeutral->SetXTitle("p_{T trigger}");
1502       outputContainer->Add(fhXEUeRightNeutral) ;
1503
1504       fhPtHbpXEUeLeftNeutral  =
1505       new TH2F("hHbpXEUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",
1506                nptbins,ptmin,ptmax,200,0.,10.);
1507       fhPtHbpXEUeLeftNeutral->SetYTitle("ln(1/x_{E})");
1508       fhPtHbpXEUeLeftNeutral->SetXTitle("p_{T trigger}");
1509       outputContainer->Add(fhPtHbpXEUeLeftNeutral) ;
1510
1511       fhPtHbpXEUeRightNeutral  =
1512       new TH2F("hHbpXEUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",
1513                nptbins,ptmin,ptmax,200,0.,10.);
1514       fhPtHbpXEUeRightNeutral->SetYTitle("ln(1/x_{E})");
1515       fhPtHbpXEUeRightNeutral->SetXTitle("p_{T trigger}");
1516       outputContainer->Add(fhPtHbpXEUeRightNeutral) ;
1517
1518       fhZTUeLeftNeutral  =
1519       new TH2F("hZTUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
1520                nptbins,ptmin,ptmax,140,0.,2.);
1521       fhZTUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1522       fhZTUeLeftNeutral->SetXTitle("p_{T trigger}");
1523       outputContainer->Add(fhZTUeLeftNeutral) ;
1524
1525       fhZTUeRightNeutral  =
1526       new TH2F("hZTUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
1527                nptbins,ptmin,ptmax,200,0.,2.);
1528       fhZTUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1529       fhZTUeRightNeutral->SetXTitle("p_{T trigger}");
1530       outputContainer->Add(fhZTUeRightNeutral) ;
1531
1532       fhPtHbpZTUeLeftNeutral  =
1533       new TH2F("hHbpZTUeNeutralLeft","#xi = ln(1/z_{T}) with neutral UE left side of trigger",
1534                nptbins,ptmin,ptmax,200,0.,10.);
1535       fhPtHbpZTUeLeftNeutral->SetYTitle("ln(1/z_{T})");
1536       fhPtHbpZTUeLeftNeutral->SetXTitle("p_{T trigger}");
1537       outputContainer->Add(fhPtHbpZTUeLeftNeutral) ;
1538
1539       fhPtHbpZTUeRightNeutral  =
1540       new TH2F("hHbpZTUeNeutralRight","#xi = ln(1/z_{T}) with neutral UE right side of trigger",
1541                nptbins,ptmin,ptmax,200,0.,10.);
1542       fhPtHbpZTUeRightNeutral->SetYTitle("ln(1/z_{T})");
1543       fhPtHbpZTUeRightNeutral->SetXTitle("p_{T trigger}");
1544       outputContainer->Add(fhPtHbpZTUeRightNeutral) ;
1545
1546     }
1547
1548   }//Correlation with neutral hadrons
1549
1550   //if data is MC, fill more histograms
1551   if(IsDataMC())
1552   {
1553     fh2phiLeadingParticle=new TH2F("h2phiLeadingParticle","#phi resolustion for trigger particles",nptbins,ptmin,ptmax,100,-1,1);
1554     fh2phiLeadingParticle->GetXaxis()->SetTitle("p_{T gen Leading} (GeV/c)");
1555     fh2phiLeadingParticle->GetYaxis()->SetTitle("(#phi_{rec}-#phi_{gen})/#phi_{gen}");
1556
1557     fhMCPtLeading  = new TH1F ("hMCPtLeading","MC : p_T distribution of leading particles", nptbins,ptmin,ptmax);
1558     fhMCPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
1559
1560     fhMCPhiLeading  = new TH2F ("hMCPhiLeading","MC : #phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1561     fhMCPhiLeading->SetYTitle("#phi (rad)");
1562
1563     fhMCEtaLeading  = new TH2F ("hMCEtaLeading","MC : #eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax);
1564     fhMCEtaLeading->SetYTitle("#eta ");
1565
1566
1567     fhMCEtaCharged  = new TH2F
1568     ("hMCEtaCharged","MC #eta_{h^{#pm}}  vs p_{T #pm}",
1569      nptbins,ptmin,ptmax,netabins,etamin,etamax);
1570     fhMCEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
1571     fhMCEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
1572
1573     fhMCPhiCharged  = new TH2F
1574     ("hMCPhiCharged","#MC phi_{h^{#pm}}  vs p_{T #pm}",
1575      200,ptmin,ptmax,nphibins,phimin,phimax);
1576     fhMCPhiCharged->SetYTitle("MC #phi_{h^{#pm}} (rad)");
1577     fhMCPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
1578
1579     fhMCDeltaPhiDeltaEtaCharged  = new TH2F
1580     ("hMCDeltaPhiDeltaEtaCharged","#MC phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
1581      140,-2.,5.,200,-2,2);
1582     fhMCDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
1583     fhMCDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1584
1585     fhMCDeltaEtaCharged  = new TH2F
1586     ("hMCDeltaEtaCharged","MC #eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger} and p_{T assoc}",
1587      nptbins,ptmin,ptmax,200,-2,2);
1588     fhMCDeltaEtaCharged->SetYTitle("#Delta #eta");
1589     fhMCDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
1590
1591     fhMCDeltaPhiCharged  = new TH2F
1592     ("hMCDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
1593      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1594     fhMCDeltaPhiCharged->SetYTitle("#Delta #phi");
1595     fhMCDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
1596
1597     fhMCDeltaPhiChargedPt  = new TH2F
1598     ("hMCDeltaPhiChargedPt","MC #phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
1599      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1600     fhMCDeltaPhiChargedPt->SetYTitle("#Delta #phi");
1601     fhMCDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
1602
1603     fhMCPtXECharged  =
1604     new TH2F("hMCPtXECharged","x_{E}",
1605              nptbins,ptmin,ptmax,200,0.,2.);
1606     fhMCPtXECharged->SetYTitle("x_{E}");
1607     fhMCPtXECharged->SetXTitle("p_{T trigger}");
1608
1609     fhMCPtXEUeCharged  =
1610     new TH2F("hMCPtXEUeCharged","x_{E}",
1611              nptbins,ptmin,ptmax,200,0.,2.);
1612     fhMCPtXEUeCharged->SetYTitle("x_{E}");
1613     fhMCPtXEUeCharged->SetXTitle("p_{T trigger}");
1614
1615     fhMCPtHbpXECharged  =
1616     new TH2F("hMCHbpXECharged","MC #xi = ln(1/x_{E}) with charged hadrons",
1617              nptbins,ptmin,ptmax,200,0.,10.);
1618     fhMCPtHbpXECharged->SetYTitle("ln(1/x_{E})");
1619     fhMCPtHbpXECharged->SetXTitle("p_{T trigger}");
1620
1621     fhMCPtHbpXEUeCharged =
1622     new TH2F("hMCPtHbpXEUeCharged","#xi = ln(1/x_{E}) with charged hadrons,Underlying Event",
1623              nptbins,ptmin,ptmax,200,0.,10.);
1624     fhMCPtHbpXEUeCharged->SetYTitle("ln(1/x_{E})");
1625     fhMCPtHbpXEUeCharged->SetXTitle("p_{T trigger}");
1626
1627     fhMCUePart  =
1628     new TH1F("hMCUePart","MC UE particles distribution vs pt trig",
1629              nptbins,ptmin,ptmax);
1630     fhMCUePart->SetYTitle("dNch");
1631     fhMCUePart->SetXTitle("p_{T trigger}");
1632
1633     fhMCPtZTCharged  =
1634     new TH2F("hMCPtZTCharged","z_{T}",
1635              nptbins,ptmin,ptmax,200,0.,2.);
1636     fhMCPtZTCharged->SetYTitle("z_{T}");
1637     fhMCPtZTCharged->SetXTitle("p_{T trigger}");
1638
1639     fhMCPtHbpZTCharged  =
1640     new TH2F("hMCHbpZTCharged","MC #xi = ln(1/z_{T}) with charged hadrons",
1641              nptbins,ptmin,ptmax,200,0.,10.);
1642     fhMCPtHbpZTCharged->SetYTitle("ln(1/z_{T})");
1643     fhMCPtHbpZTCharged->SetXTitle("p_{T trigger}");
1644
1645     fhMCPtTrigPout  =
1646     new TH2F("hMCPtTrigPout","AOD MC Pout with triggers",
1647              nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
1648     fhMCPtTrigPout->SetYTitle("p_{out} (GeV/c)");
1649     fhMCPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
1650
1651     fhMCPtAssocDeltaPhi  =
1652     new TH2F("hMCPtAssocDeltaPhi","AOD MC delta phi with associated charged hadrons",
1653              nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1654     fhMCPtAssocDeltaPhi->SetYTitle("#Delta #phi");
1655     fhMCPtAssocDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
1656
1657     outputContainer->Add(fh2phiLeadingParticle);
1658     outputContainer->Add(fhMCPtLeading);
1659     outputContainer->Add(fhMCPhiLeading);
1660     outputContainer->Add(fhMCEtaLeading);
1661     outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged);
1662     outputContainer->Add(fhMCPhiCharged) ;
1663     outputContainer->Add(fhMCEtaCharged) ;
1664     outputContainer->Add(fhMCDeltaEtaCharged) ;
1665     outputContainer->Add(fhMCDeltaPhiCharged) ;
1666
1667     outputContainer->Add(fhMCDeltaPhiChargedPt) ;
1668     outputContainer->Add(fhMCPtXECharged) ;
1669     outputContainer->Add(fhMCPtXEUeCharged) ;
1670     outputContainer->Add(fhMCPtZTCharged) ;
1671     outputContainer->Add(fhMCPtHbpXECharged) ;
1672     outputContainer->Add(fhMCPtHbpXEUeCharged);
1673     outputContainer->Add(fhMCUePart);
1674     outputContainer->Add(fhMCPtHbpZTCharged) ;
1675     outputContainer->Add(fhMCPtTrigPout) ;
1676     outputContainer->Add(fhMCPtAssocDeltaPhi) ;
1677   } //for MC histogram
1678
1679   if(DoOwnMix())
1680   {
1681     //create event containers
1682
1683     if(!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForTracksExists()))
1684     {
1685       Int_t nvz = GetNZvertBin();
1686       Int_t nrp = GetNRPBin();
1687       Int_t nce = GetNCentrBin();
1688
1689       fListMixTrackEvents= new TList*[nvz*nrp*nce] ;
1690
1691       for( Int_t ice = 0 ; ice < nce ; ice++ )
1692       {
1693         for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
1694         {
1695           for( Int_t irp = 0 ; irp < nrp ; irp++ )
1696           {
1697             Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
1698
1699             //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
1700             //       ic,iz, irp, bin);
1701
1702             fListMixTrackEvents[bin] = new TList() ;
1703             fListMixTrackEvents[bin]->SetOwner(kFALSE);
1704           }
1705         }
1706       }
1707     }
1708
1709     fhPtLeadingMixed  = new TH1F ("hPtLeadingMixed","p_T distribution of leading particles, used for mixing", nptbins,ptmin,ptmax);
1710     fhPtLeadingMixed->SetXTitle("p_{T}^{trig} (GeV/c)");
1711
1712     fhPtLeadingMixedBin  = new TH2F ("hPtLeadingMixedBin","p_T distribution of leading particles vs mixing bin", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
1713     fhPtLeadingMixedBin->SetXTitle("p_{T}^{trig} (GeV/c)");
1714     fhPtLeadingMixedBin->SetYTitle("Bin");
1715
1716     fhPhiLeadingMixed  = new TH2F ("hPhiLeadingMixed","#phi distribution of leading Particles, used for mixing",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1717     fhPhiLeadingMixed->SetYTitle("#phi (rad)");
1718
1719     fhEtaLeadingMixed  = new TH2F ("hEtaLeadingMixed","#eta distribution of leading, used for mixing",nptbins,ptmin,ptmax, netabins,etamin,etamax);
1720     fhEtaLeadingMixed->SetYTitle("#eta ");
1721
1722     outputContainer->Add(fhPtLeadingMixed);
1723     outputContainer->Add(fhPtLeadingMixedBin);
1724     outputContainer->Add(fhPhiLeadingMixed);
1725     outputContainer->Add(fhEtaLeadingMixed);
1726
1727     // Fill the cluster pool only in isolation analysis
1728     if( OnlyIsolated() && (!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForCaloExists())))
1729     {
1730       Int_t nvz = GetNZvertBin();
1731       Int_t nrp = GetNRPBin();
1732       Int_t nce = GetNCentrBin();
1733
1734       fListMixCaloEvents= new TList*[nvz*nrp*nce] ;
1735
1736       for( Int_t ice = 0 ; ice < nce ; ice++ )
1737       {
1738         for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
1739         {
1740           for( Int_t irp = 0 ; irp < nrp ; irp++ )
1741           {
1742             Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
1743
1744             //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
1745             //       ic,iz, irp, bin);
1746
1747             fListMixCaloEvents[bin] = new TList() ;
1748             fListMixCaloEvents[bin]->SetOwner(kFALSE);
1749           }
1750         }
1751       }
1752     }
1753
1754     //Init the list in the reader if not done previously
1755     if(fUseMixStoredInReader)
1756     {
1757       if( !GetReader()->ListWithMixedEventsForTracksExists() )
1758         GetReader()->SetListWithMixedEventsForTracks(fListMixTrackEvents);
1759
1760       if( !GetReader()->ListWithMixedEventsForCaloExists()   )
1761         GetReader()->SetListWithMixedEventsForCalo  (fListMixCaloEvents );
1762     }
1763
1764     fhEventBin=new TH1I("hEventBin","Number of real events per bin(cen,vz,rp)",
1765                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
1766                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
1767     fhEventBin->SetXTitle("bin");
1768     outputContainer->Add(fhEventBin) ;
1769
1770     fhEventMixBin=new TH1I("hEventMixBin","Number of events  per bin(cen,vz,rp)",
1771                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
1772                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
1773     fhEventMixBin->SetXTitle("bin");
1774     outputContainer->Add(fhEventMixBin) ;
1775
1776     fhNtracksAll=new TH1F("hNtracksAll","Number of tracks w/o event trigger",2000,0,2000);
1777     outputContainer->Add(fhNtracksAll);
1778
1779     fhNtracksTrigger=new TH1F("hNtracksTriggerEvent","Number of tracks w/ event trigger",2000,0,2000);
1780     outputContainer->Add(fhNtracksTrigger);
1781
1782     fhNtracksMB=new TH1F("hNtracksMBEvent","Number of tracks w/ event trigger kMB",2000,0,2000);
1783     outputContainer->Add(fhNtracksMB);
1784
1785     fhMixDeltaPhiCharged  = new TH2F
1786     ("hMixDeltaPhiCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
1787      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1788     fhMixDeltaPhiCharged->SetYTitle("#Delta #phi");
1789     fhMixDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
1790     outputContainer->Add(fhMixDeltaPhiCharged);
1791
1792     fhMixDeltaPhiDeltaEtaCharged  = new TH2F
1793     ("hMixDeltaPhiDeltaEtaCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
1794      ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
1795     fhMixDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
1796     fhMixDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1797     outputContainer->Add(fhMixDeltaPhiDeltaEtaCharged);
1798
1799     fhMixXECharged  =
1800     new TH2F("hMixXECharged","Mixed event : x_{E} for charged tracks",
1801              nptbins,ptmin,ptmax,200,0.,2.);
1802     fhMixXECharged->SetYTitle("x_{E}");
1803     fhMixXECharged->SetXTitle("p_{T trigger}");
1804     outputContainer->Add(fhMixXECharged);
1805
1806     fhMixHbpXECharged  =
1807     new TH2F("hMixHbpXECharged","mixed event : #xi = ln(1/x_{E}) with charged hadrons",
1808              nptbins,ptmin,ptmax,200,0.,10.);
1809     fhMixHbpXECharged->SetYTitle("ln(1/x_{E})");
1810     fhMixHbpXECharged->SetXTitle("p_{T trigger}");
1811     outputContainer->Add(fhMixHbpXECharged);
1812
1813     fhMixDeltaPhiChargedAssocPtBin         = new TH2F*[fNAssocPtBins] ;
1814     fhMixDeltaPhiChargedAssocPtBinDEta08   = new TH2F*[fNAssocPtBins] ;
1815     fhMixDeltaPhiDeltaEtaChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
1816
1817     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
1818     {
1819       fhMixDeltaPhiChargedAssocPtBin[i] = new TH2F(Form("hMixDeltaPhiChargedAssocPtBin%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1820                                                    Form("Mixed event #Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1821                                                    nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
1822       fhMixDeltaPhiChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
1823       fhMixDeltaPhiChargedAssocPtBin[i]->SetYTitle("#Delta #phi");
1824
1825       fhMixDeltaPhiChargedAssocPtBinDEta08[i] = new TH2F(Form("hMixDeltaPhiDeltaEta0.8ChargedAssocPtBin%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1826                                                          Form("Mixed event #Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f] for #Delta #eta > 0.8", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1827                                                          nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
1828       fhMixDeltaPhiChargedAssocPtBinDEta08[i]->SetXTitle("p_{T trigger}");
1829       fhMixDeltaPhiChargedAssocPtBinDEta08[i]->SetYTitle("#Delta #phi");
1830
1831       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i] = new TH2F(Form("hMixDeltaPhiDeltaEtaChargedAssocPtBin%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1832                                                            Form("Mixed event #Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1833                                                            ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
1834       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]->SetXTitle("#Delta #phi");
1835       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]->SetYTitle("#Delta #eta");
1836
1837       outputContainer->Add(fhMixDeltaPhiChargedAssocPtBin[i]);
1838       outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta08[i]);
1839       outputContainer->Add(fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]);
1840
1841     }
1842   }
1843
1844   return outputContainer;
1845
1846 }
1847
1848 //_________________________________________________________________________________________________
1849 Bool_t AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum(const AliAODPWG4Particle* trigger,
1850                                                              TLorentzVector & mom1,
1851                                                              TLorentzVector & mom2)
1852 {
1853   // Get the momentum of the pi0/eta assigned decay photons
1854   // In case of pi0/eta trigger, we may want to check their decay correlation,
1855   // get their decay children
1856
1857   Int_t indexPhoton1 = trigger->GetCaloLabel(0);
1858   Int_t indexPhoton2 = trigger->GetCaloLabel(1);
1859   Float_t ptTrig     = trigger->Pt();
1860
1861   if(indexPhoton1!=-1 || indexPhoton2!=-1) return kFALSE;
1862
1863   if(GetDebug() > 1)
1864     printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
1865
1866   TObjArray * clusters  = 0x0 ;
1867   if(trigger->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
1868   else                                clusters = GetPHOSClusters()  ;
1869
1870   for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++)
1871   {
1872     AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));
1873     if(photon->GetID()==indexPhoton1)
1874     {
1875       photon->GetMomentum(mom1,GetVertex(0)) ;
1876       if(ptTrig) fhPtPi0DecayRatio->Fill(ptTrig, mom1.Pt()/ptTrig);
1877     }
1878     if(photon->GetID()==indexPhoton2)
1879     {
1880       photon->GetMomentum(mom1,GetVertex(0)) ;
1881       if(ptTrig > 0) fhPtPi0DecayRatio->Fill(ptTrig, mom2.Pt()/ptTrig);
1882     }
1883
1884     if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - Photon1 = %f, Photon2 = %f \n", mom1.Pt(), mom2.Pt());
1885
1886   } //cluster loop
1887
1888   return kTRUE;
1889
1890 }
1891
1892 //____________________________________________________
1893 void AliAnaParticleHadronCorrelation::InitParameters()
1894 {
1895
1896   //Initialize the parameters of the analysis.
1897   SetInputAODName("Particle");
1898   SetAODObjArrayName("Hadrons");
1899   AddToHistogramsName("AnaHadronCorr_");
1900
1901   SetPtCutRange(0.,300);
1902   fDeltaPhiMinCut       = 1.5 ;
1903   fDeltaPhiMaxCut       = 4.5 ;
1904   fSelectIsolated       = kFALSE;
1905   fMakeSeveralUE        = kFALSE;
1906   fUeDeltaPhiMinCut     = 1. ;
1907   fUeDeltaPhiMaxCut     = 1.5 ;
1908
1909   fNeutralCorr          = kFALSE ;
1910   fPi0Trigger           = kFALSE ;
1911   fDecayTrigger         = kFALSE ;
1912   fHMPIDCorrelation     = kFALSE ;
1913
1914   fMakeAbsoluteLeading  = kTRUE;
1915   fMakeNearSideLeading  = kFALSE;
1916
1917   fNAssocPtBins         = 9   ;
1918   fAssocPtBinLimit[0]   = 0.2 ;
1919   fAssocPtBinLimit[1]   = 2.0 ;
1920   fAssocPtBinLimit[2]   = 4.0 ;
1921   fAssocPtBinLimit[3]   = 6.0 ;
1922   fAssocPtBinLimit[4]   = 8.0 ;
1923   fAssocPtBinLimit[5]   = 10. ;
1924   fAssocPtBinLimit[6]   = 12. ;
1925   fAssocPtBinLimit[7]   = 15. ;
1926   fAssocPtBinLimit[8]   = 25. ;
1927   fAssocPtBinLimit[9]   = 50. ;
1928
1929   fUseMixStoredInReader = kTRUE;
1930
1931   fM02MinCut   = -1 ;
1932   fM02MaxCut   = -1 ;
1933
1934 }
1935
1936 //__________________________________________________________
1937 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
1938 {
1939   //Particle-Hadron Correlation Analysis, fill AODs
1940
1941   if(!GetInputAODBranch())
1942   {
1943     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1944     abort();
1945   }
1946
1947   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
1948   {
1949     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());
1950     abort();
1951   }
1952
1953   if(GetDebug() > 1)
1954   {
1955     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");
1956     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
1957     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n",   GetCTSTracks()    ->GetEntriesFast());
1958     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
1959     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n",  GetPHOSClusters() ->GetEntriesFast());
1960   }
1961
1962   //Get the vertex and check it is not too large in z
1963   Double_t v[3] = {0,0,0}; //vertex ;
1964   GetReader()->GetVertex(v);
1965   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
1966
1967   // Fill the pool with tracks if requested
1968   if(DoOwnMix())
1969   {
1970     FillChargedEventMixPool();
1971     if(OnlyIsolated())
1972       FillNeutralEventMixPool();
1973   }
1974
1975   //Loop on stored AOD particles, find leading trigger
1976   Double_t ptTrig      = fMinTriggerPt ;
1977   fLeadingTriggerIndex = -1 ;
1978   Int_t    naod        = GetInputAODBranch()->GetEntriesFast() ;
1979   for(Int_t iaod = 0; iaod < naod ; iaod++)
1980   {
1981     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1982
1983     // Vertex cut in case of mixing
1984     Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
1985     if(check ==  0) continue;
1986     if(check == -1) return;
1987
1988     // find the leading particles with highest momentum
1989     if (particle->Pt() > ptTrig)
1990     {
1991       ptTrig               = particle->Pt() ;
1992       fLeadingTriggerIndex = iaod ;
1993     }
1994   }// finish search of leading trigger particle
1995
1996
1997   //Do correlation with leading particle
1998   if(fLeadingTriggerIndex >= 0)
1999   {
2000
2001     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(fLeadingTriggerIndex));
2002
2003     //check if the particle is isolated or if we want to take the isolation into account
2004     if(OnlyIsolated() && !particle->IsIsolated()) return;
2005
2006     //Make correlation with charged hadrons
2007     Bool_t okcharged = kTRUE;
2008     Bool_t okneutral = kTRUE;
2009     if(GetReader()->IsCTSSwitchedOn() )
2010       okcharged = MakeChargedCorrelation(particle, GetCTSTracks(),kFALSE);
2011
2012     TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
2013     if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
2014       okneutral = MakeNeutralCorrelation(particle, pi0list,kFALSE);
2015
2016   }//Correlate leading
2017
2018   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
2019
2020 }
2021
2022 //_________________________________________________________________
2023 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
2024 {
2025   //Particle-Hadron Correlation Analysis, fill histograms
2026
2027   if(!GetInputAODBranch())
2028   {
2029     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
2030     abort();
2031   }
2032
2033   if(GetDebug() > 1)
2034   {
2035     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
2036     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
2037   }
2038
2039   //Get the vertex and check it is not too large in z
2040   Double_t v[3] = {0,0,0}; //vertex ;
2041   GetReader()->GetVertex(v);
2042   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
2043
2044   //Loop on stored AOD particles, find leading
2045   Double_t ptTrig    = fMinTriggerPt;
2046   if(fLeadingTriggerIndex < 0)
2047   {
2048     //Search leading if not done before
2049     Int_t    naod      = GetInputAODBranch()->GetEntriesFast() ;
2050     for(Int_t iaod = 0; iaod < naod ; iaod++)
2051     {    //loop on input trigger AOD file
2052       AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
2053
2054       // Vertex cut in case of mixing
2055       Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
2056       if(check ==  0) continue;
2057       if(check == -1) return;
2058
2059       //check if the particle is isolated or if we want to take the isolation into account
2060       if(OnlyIsolated() && !particle->IsIsolated()) continue;
2061
2062       //find the leading particles with highest momentum
2063       if (particle->Pt() > ptTrig)
2064       {
2065         ptTrig               = particle->Pt() ;
2066         fLeadingTriggerIndex = iaod ;
2067       }
2068
2069     }// Finish search of leading trigger particle
2070   }// Search leading if not done before
2071
2072   if(fLeadingTriggerIndex >= 0 )
2073   { //using trigger particle to do correlations
2074
2075     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(fLeadingTriggerIndex));
2076
2077     // check if it was a calorimeter cluster and if the SS cut was requested, if so, apply it
2078     Int_t clID1  = particle->GetCaloLabel(0) ;
2079     Int_t clID2  = particle->GetCaloLabel(1) ; // for photon clusters should not be set.
2080     //printf("Leading for for %s: id1 %d, id2 %d, min %f, max %f, det %s\n",
2081     //       GetInputAODName().Data(),clID1,clID2,fM02MinCut,fM02MaxCut,(particle->GetDetector()).Data());
2082
2083     if(clID1 > 0 && clID2 < 0 && fM02MaxCut > 0 && fM02MinCut > 0)
2084     {
2085       Int_t iclus = -1;
2086       TObjArray* clusters = 0x0;
2087       if     (particle->GetDetector() == "EMCAL") clusters = GetEMCALClusters();
2088       else if(particle->GetDetector() == "PHOS" ) clusters = GetPHOSClusters();
2089
2090       if(clusters)
2091       {
2092         AliVCluster *cluster = FindCluster(clusters,clID1,iclus);
2093         Float_t m02 = cluster->GetM02();
2094         //printf("\t Check m02 = %2.2f\n",m02);
2095         if(m02 > fM02MaxCut || m02 < fM02MinCut)
2096         {
2097           //printf("\t \t Not accepted\n");
2098           return;
2099         }
2100       }
2101     }
2102
2103     // Check if trigger is in fiducial region
2104     if(IsFiducialCutOn())
2105     {
2106       Bool_t in = GetFiducialCut()->IsInFiducialCut(*particle->Momentum(),particle->GetDetector()) ;
2107       if(! in ) return ;
2108     }
2109
2110     // Check if the particle is isolated or if we want to take the isolation into account
2111     if(OnlyIsolated() && !particle->IsIsolated()) return;
2112
2113     // Make correlation with charged hadrons
2114     Bool_t okcharged = kTRUE;
2115     Bool_t okneutral = kTRUE;
2116     if(GetReader()->IsCTSSwitchedOn() )
2117     {
2118       okcharged = MakeChargedCorrelation(particle, GetCTSTracks(),kTRUE);
2119       if(IsDataMC())
2120       {
2121         MakeMCChargedCorrelation(particle);
2122       }
2123     }
2124
2125     TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
2126     if(fNeutralCorr && pi0list)
2127     {
2128       if(pi0list->GetEntriesFast() > 0)
2129         okneutral = MakeNeutralCorrelation(particle, pi0list,kTRUE);
2130     }
2131
2132     // Fill leading particle histogram if correlation went well and
2133     // no problem was found, like not absolute leading, or bad vertex in mixing.
2134     if(okcharged && okneutral)
2135     {
2136       Float_t pt = particle->Pt();
2137       fhPtLeading->Fill(pt);
2138       fhPtLeadingBin->Fill(pt,GetEventMixBin());
2139
2140       Float_t phi = particle->Phi();
2141       if(phi<0)phi+=TMath::TwoPi();
2142       fhPhiLeading->Fill(pt, phi);
2143
2144       fhEtaLeading->Fill(pt, particle->Eta());
2145       //printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Leading particle : pt %f, eta %f, phi %f\n",particle->Pt(),particle->Eta(),phi);
2146
2147       Float_t cen = GetEventCentrality();
2148       Float_t ep  = GetEventPlaneAngle();
2149
2150       fhPtLeadingCentrality        ->Fill(pt,cen);
2151       fhPtLeadingEventPlane        ->Fill(pt,ep);
2152       fhLeadingEventPlaneCentrality->Fill(cen,ep);
2153
2154     }//ok charged && neutral
2155   }//Aod branch loop
2156
2157   //Reinit for next event
2158   fLeadingTriggerIndex = -1;
2159
2160   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
2161 }
2162
2163 //___________________________________________________________________________________________________________
2164 Bool_t  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle,
2165                                                                 const TObjArray* pl, const Bool_t bFillHisto)
2166 {
2167   // Charged Hadron Correlation Analysis
2168   if(GetDebug() > 1)
2169     printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
2170
2171   Float_t phiTrig = aodParticle->Phi();
2172   Float_t etaTrig = aodParticle->Eta();
2173   Float_t ptTrig  = aodParticle->Pt();
2174   Bool_t   decay  = aodParticle->IsTagged();
2175
2176   Float_t pt       = -100. ;
2177   Float_t zT       = -100. ;
2178   Float_t xE       = -100. ;
2179   Float_t hbpXE    = -100. ;
2180   Float_t hbpZT    = -100. ;
2181   Float_t phi      = -100. ;
2182   Float_t eta      = -100. ;
2183   Float_t pout     = -100. ;
2184   Float_t deltaPhi = -100. ;
2185
2186   TVector3 p3;
2187   TLorentzVector photonMom ;
2188   TObjArray * reftracks = 0x0;
2189   Int_t nrefs           = 0;
2190   Int_t nTracks         = GetCTSTracks()->GetEntriesFast() ;
2191
2192   // Mixed event settings
2193   Int_t evtIndex11   = -1 ; // cluster trigger or pi0 trigger
2194   Int_t evtIndex12   = -1 ; // pi0 trigger
2195   Int_t evtIndex13   = -1 ; // charged trigger
2196
2197   Double_t v[3]      = {0,0,0}; //vertex ;
2198   GetReader()->GetVertex(v);
2199
2200   if (GetMixedEvent())
2201   {
2202     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
2203     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
2204     evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
2205   }
2206
2207   // In case of pi0/eta trigger, we may want to check their decay correlation,
2208   // get their decay children
2209   TLorentzVector decayMom1;
2210   TLorentzVector decayMom2;
2211   Bool_t decayFound = kFALSE;
2212   if(fPi0Trigger && bFillHisto) decayFound = GetDecayPhotonMomentum(aodParticle,decayMom1, decayMom2);
2213
2214   //-----------------------------------------------------------------------
2215   //Track loop, select tracks with good pt, phi and fill AODs or histograms
2216   //-----------------------------------------------------------------------
2217
2218   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ )
2219   {
2220     AliVTrack * track = (AliVTrack *) (pl->At(ipr)) ;
2221
2222     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
2223     p3.SetXYZ(mom[0],mom[1],mom[2]);
2224     pt   = p3.Pt();
2225     eta  = p3.Eta();
2226     phi  = p3.Phi() ;
2227     if(phi < 0) phi+=TMath::TwoPi();
2228
2229     //Select only hadrons in pt range
2230     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
2231
2232     //remove trigger itself for correlation when use charged triggers
2233     if( track->GetID() == aodParticle->GetTrackLabel(0) || track->GetID() == aodParticle->GetTrackLabel(1) ||
2234         track->GetID() == aodParticle->GetTrackLabel(2) || track->GetID() == aodParticle->GetTrackLabel(3)   )
2235       continue ;
2236
2237     //jump out this event if near side associated particle pt larger than trigger
2238     if (fMakeNearSideLeading)
2239     {
2240       if(pt > ptTrig && TMath::Abs(phi-phiTrig) < TMath::PiOver2())  return kFALSE;
2241     }
2242     //jump out this event if there is any other particle with pt larger than trigger
2243     else if(fMakeAbsoluteLeading)
2244     {
2245       if(pt > ptTrig)  return kFALSE;
2246     }
2247
2248     //Only for mixed event
2249     Int_t evtIndex2 = 0 ;
2250     if (GetMixedEvent())
2251     {
2252       evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
2253       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
2254         continue ;
2255       //vertex cut
2256       if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
2257         return kFALSE;
2258     }
2259
2260     // Fill Histograms
2261     if(bFillHisto)
2262     {
2263
2264       if(GetDebug() > 2 )
2265         printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
2266
2267       // Set the pt associated bin for the defined bins
2268       Int_t assocBin   = -1;
2269       for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2270       {
2271         if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i;
2272       }
2273
2274       // Azimuthal Angle
2275       // calculate deltaPhi for later, shift when needed
2276       FillChargedAngularCorrelationHistograms(pt,  ptTrig,  assocBin, phi, phiTrig,  deltaPhi,
2277                                               eta, etaTrig, decay, track->GetHMPIDsignal(),nTracks);
2278
2279       // Imbalance zT/xE/pOut
2280       zT = pt/ptTrig ;
2281       if(zT > 0 ) hbpZT = TMath::Log(1./zT);
2282       else        hbpZT =-100;
2283
2284       xE   =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
2285       //if(xE <0.)xE =-xE;
2286       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
2287       else        hbpXE =-100;
2288
2289       pout = pt*TMath::Sin(deltaPhi) ;
2290
2291       //delta phi cut for momentum imbalance correlation
2292       if      ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   )
2293       {
2294
2295         FillChargedMomentumImbalanceHistograms(ptTrig, pt, xE, hbpXE, zT, hbpZT, pout, deltaPhi,
2296                                                nTracks, track->Charge(), assocBin, decay);
2297
2298       }
2299       if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
2300       { //UE study
2301
2302         FillChargedUnderlyingEventHistograms(ptTrig, pt, deltaPhi, nTracks);
2303
2304         fhUePart->Fill(ptTrig);
2305
2306       }
2307
2308       if(fPi0Trigger && decayFound)
2309         FillDecayPhotonCorrelationHistograms(pt, phi, decayMom1,decayMom2, kTRUE) ;
2310
2311       //several UE calculation
2312       if(fMakeSeveralUE) FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,xE,hbpXE,zT,hbpZT,deltaPhi);
2313
2314     } //Fill histogram
2315     else
2316     {
2317       nrefs++;
2318       if(nrefs==1)
2319       {
2320         reftracks = new TObjArray(0);
2321         TString trackname = Form("%s+Tracks", GetAODObjArrayName().Data());
2322         reftracks->SetName(trackname.Data());
2323         reftracks->SetOwner(kFALSE);
2324       }
2325
2326       reftracks->Add(track);
2327
2328     }//aod particle loop
2329   }// track loop
2330
2331   //Fill AOD with reference tracks, if not filling histograms
2332   if(!bFillHisto && reftracks)
2333   {
2334     aodParticle->AddObjArray(reftracks);
2335   }
2336
2337   //Own mixed event, add event and remove previous or fill the mixed histograms
2338   if(DoOwnMix() && bFillHisto)
2339   {
2340       MakeChargedMixCorrelation(aodParticle);
2341   }
2342
2343   return kTRUE;
2344
2345 }
2346
2347
2348 //_________________________________________________________________________________________________________
2349 void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
2350 {
2351   // Mix current trigger with tracks in another MB event
2352
2353   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Make trigger particle - charged hadron mixed event correlation \n");
2354
2355   if(GetMixedEvent()) return;  // This is not the mixed event from general mixing frame
2356
2357   // Get the event with similar caracteristics
2358   //printf("MakeChargedMixCorrelation for %s\n",GetInputAODName().Data());
2359
2360   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
2361
2362   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
2363
2364   if(!inputHandler) return;
2365
2366   if(!(inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask())) return;
2367
2368   // Get the pool, check if it exits
2369   Int_t eventBin = GetEventMixBin();
2370
2371   fhEventBin->Fill(eventBin);
2372
2373   //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2374   if(eventBin < 0) return;
2375
2376   TList * pool     = 0;
2377   TList * poolCalo = 0;
2378   if(fUseMixStoredInReader)
2379   {
2380     pool     = GetReader()->GetListWithMixedEventsForTracks(eventBin);
2381     poolCalo = GetReader()->GetListWithMixedEventsForCalo  (eventBin);
2382   }
2383   else
2384   {
2385     pool     = fListMixTrackEvents[eventBin];
2386     poolCalo = fListMixCaloEvents [eventBin];
2387   }
2388
2389   if(!pool) return ;
2390
2391   if(OnlyIsolated() && !poolCalo &&
2392      (GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::AliIsolationCut::kOnlyCharged))
2393     printf("AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation() - Careful, cluster pool not available\n");
2394
2395   Double_t ptTrig  = aodParticle->Pt();
2396   Double_t etaTrig = aodParticle->Eta();
2397   Double_t phiTrig = aodParticle->Phi();
2398   if(phiTrig < 0.) phiTrig+=TMath::TwoPi();
2399
2400   if(GetDebug() > 1)
2401     printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Pool bin %d size %d, leading trigger pt=%f, phi=%f, eta=%f\n",
2402            eventBin,pool->GetSize(), ptTrig,phiTrig,etaTrig);
2403
2404   Double_t ptAssoc  = -999.;
2405   Double_t phiAssoc = -999.;
2406   Double_t etaAssoc = -999.;
2407   Double_t deltaPhi = -999.;
2408   Double_t deltaEta = -999.;
2409   Double_t xE = -999.;
2410   Double_t hbpXE = -999.;
2411
2412   //Start from first event in pool except if in this same event the pool was filled
2413   Int_t ev0 = 0;
2414   if(GetReader()->GetLastTracksMixedEvent() == GetEventNumber()) ev0 = 1;
2415
2416   for(Int_t ev=ev0; ev < pool->GetSize(); ev++)
2417   {
2418     TObjArray* bgTracks = static_cast<TObjArray*>(pool->At(ev));
2419
2420     // Check if the particle is isolated in the mixed event, it not, do not fill the histograms
2421     if(OnlyIsolated())
2422     {
2423       if(pool->GetSize()!=poolCalo->GetSize())
2424         printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Different size of calo and track pools\n");
2425
2426       TObjArray* bgCalo = static_cast<TObjArray*>(poolCalo->At(ev));
2427
2428       if(!bgCalo)
2429         printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Event %d in calo pool not available?\n",ev);
2430
2431       Int_t n=0; Int_t nfrac = 0; Bool_t isolated = kFALSE; Float_t coneptsum = 0;
2432       GetIsolationCut()->MakeIsolationCut(bgTracks,bgCalo,
2433                                           GetReader(), GetCaloPID(),
2434                                           kFALSE, aodParticle, "",
2435                                           n,nfrac,coneptsum, isolated);
2436
2437       //printf("AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation() - Isolated? %d - cone %f, ptthres %f",
2438       //       isolated,GetIsolationCut()->GetConeSize(),GetIsolationCut()->GetPtThreshold());
2439       //if(bgTracks)printf(" - n track %d", bgTracks->GetEntriesFast());
2440       //printf("\n");
2441
2442       if(!isolated) continue ;
2443     }
2444
2445     fhEventMixBin->Fill(eventBin);
2446
2447     Int_t nTracks=bgTracks->GetEntriesFast();
2448     //printf("\t Read Pool event %d, nTracks %d\n",ev,nTracks);
2449
2450     //Check if it is leading if mixed event
2451     if(fMakeNearSideLeading || fMakeAbsoluteLeading)
2452     {
2453       Bool_t leading = kTRUE;
2454       for(Int_t jlead = 0;jlead <nTracks; jlead++ )
2455       {
2456         AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(jlead) ;
2457
2458         ptAssoc  = track->Pt();
2459         phiAssoc = track->Phi() ;
2460
2461         if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
2462         if (fMakeNearSideLeading)
2463         {
2464           if(ptAssoc > ptTrig && TMath::Abs(phiAssoc-phiTrig) < TMath::PiOver2())
2465           {
2466             leading = kFALSE;
2467             break;
2468           }
2469         }
2470         //jump out this event if there is any other particle with pt larger than trigger
2471         else if(fMakeAbsoluteLeading)
2472         {
2473           if(ptAssoc > ptTrig)
2474           {
2475             leading = kFALSE;
2476             break;
2477           }
2478         }
2479       }
2480       if(!leading) continue; // not leading, check the next event in pool
2481     }
2482
2483     fhPtLeadingMixed ->Fill(ptTrig);
2484     fhPhiLeadingMixed->Fill(ptTrig, phiTrig);
2485     fhEtaLeadingMixed->Fill(ptTrig, etaTrig);
2486     fhPtLeadingMixedBin->Fill(ptTrig,eventBin);
2487
2488     for(Int_t j1 = 0;j1 <nTracks; j1++ )
2489     {
2490       AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(j1) ;
2491
2492       if(!track) continue;
2493
2494       ptAssoc  = track->Pt();
2495       etaAssoc = track->Eta();
2496       phiAssoc = track->Phi() ;
2497       if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
2498
2499       if(IsFiducialCutOn())
2500       {
2501         Bool_t in = GetFiducialCut()->IsInFiducialCut(*aodParticle->Momentum(),"CTS") ;
2502         if(!in) continue ;
2503       }
2504
2505       deltaPhi = phiTrig-phiAssoc;
2506       if(deltaPhi < -TMath::PiOver2())  deltaPhi+=TMath::TwoPi();
2507       if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
2508       deltaEta = etaTrig-etaAssoc;
2509
2510       if(GetDebug()>0)
2511         printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation(): deltaPhi= %f, deltaEta=%f\n",deltaPhi, deltaEta);
2512
2513       // Set the pt associated bin for the defined bins
2514       Int_t assocBin   = -1;
2515       for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2516       {
2517         if(ptAssoc > fAssocPtBinLimit[i] && ptAssoc < fAssocPtBinLimit[i+1]) assocBin= i;
2518       }
2519
2520       fhMixDeltaPhiCharged        ->Fill(ptTrig,  deltaPhi);
2521       fhMixDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
2522
2523       fhMixDeltaPhiCharged        ->Fill(ptTrig,  deltaPhi);
2524       fhMixDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
2525
2526       xE   =-ptAssoc/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
2527       //if(xE <0.)xE =-xE;
2528       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
2529       else        hbpXE =-100;
2530
2531       if      ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   )
2532       {
2533         fhMixXECharged->Fill(ptTrig,xE);
2534         fhMixHbpXECharged->Fill(ptTrig,hbpXE);
2535       }
2536
2537       if(assocBin < 0) continue ; // this pt bin was not considered
2538
2539       if(TMath::Abs(deltaEta) > 0.8)
2540         fhMixDeltaPhiChargedAssocPtBinDEta08  [assocBin]->Fill(ptTrig,   deltaPhi);
2541
2542         fhMixDeltaPhiChargedAssocPtBin        [assocBin]->Fill(ptTrig,   deltaPhi);
2543         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[assocBin]->Fill(deltaPhi, deltaEta);
2544
2545     } // track loop
2546   } // mixed event loop
2547 }
2548
2549
2550 //________________________________________________________________________________________________________________
2551 Bool_t  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle,
2552                                                                 const TObjArray* pi0list, const Bool_t bFillHisto)
2553 {
2554   // Neutral Pion Correlation Analysis
2555   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",
2556                             pi0list->GetEntriesFast());
2557
2558   Int_t evtIndex11 = 0 ;
2559   Int_t evtIndex12 = 0 ;
2560   if (GetMixedEvent())
2561   {
2562     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
2563     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
2564   }
2565
2566   Float_t pt   = -100. ;
2567   Float_t zT   = -100. ;
2568   Float_t phi  = -100. ;
2569   Float_t eta  = -100. ;
2570   Float_t xE   = -100. ;
2571   Float_t hbpXE= -100. ;
2572   Float_t hbpZT= -100. ;
2573
2574   Float_t ptTrig  = aodParticle->Pt();
2575   Float_t phiTrig = aodParticle->Phi();
2576   Float_t etaTrig = aodParticle->Eta();
2577   Float_t deltaPhi= -100. ;
2578
2579   TLorentzVector photonMom ;
2580
2581   // In case of pi0/eta trigger, we may want to check their decay correlation,
2582   // get their decay children
2583   TLorentzVector decayMom1;
2584   TLorentzVector decayMom2;
2585   Bool_t decayFound = kFALSE;
2586   if(fPi0Trigger && bFillHisto) decayFound = GetDecayPhotonMomentum(aodParticle,decayMom1, decayMom2);
2587
2588   TObjArray * refpi0 = 0x0;
2589   Int_t nrefs        = 0;
2590
2591   //Loop on stored AOD pi0
2592
2593   Int_t naod = pi0list->GetEntriesFast();
2594   if(GetDebug() > 0)
2595     printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() -  aod branch entries %d\n", naod);
2596
2597   for(Int_t iaod = 0; iaod < naod ; iaod++)
2598   {
2599     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (pi0list->At(iaod));
2600
2601     Int_t evtIndex2 = 0 ;
2602     Int_t evtIndex3 = 0 ;
2603     if (GetMixedEvent())
2604     {
2605       evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
2606       evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
2607
2608       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 ||
2609           evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
2610         continue ;
2611     }
2612
2613     pt  = pi0->Pt();
2614
2615     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
2616
2617     //jump out this event if near side associated particle pt larger than trigger
2618     if (fMakeNearSideLeading)
2619     {
2620       if(pt > ptTrig && TMath::Abs(phi-phiTrig) < TMath::PiOver2())  return kFALSE;
2621     }
2622     //jump out this event if there is any other particle with pt larger than trigger
2623     else if(fMakeAbsoluteLeading)
2624     {
2625       if(pt > ptTrig)  return kFALSE;
2626     }
2627
2628     if(bFillHisto)
2629     {
2630       phi = pi0->Phi() ;
2631       eta = pi0->Eta() ;
2632
2633       FillNeutralAngularCorrelationHistograms(pt, ptTrig, phi, phiTrig, deltaPhi, eta, etaTrig);
2634
2635       zT  = pt/ptTrig ;
2636       xE  =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
2637
2638       //if(xE <0.)xE =-xE;
2639
2640       hbpXE = -100;
2641       hbpZT = -100;
2642
2643       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
2644       if(zT > 0 ) hbpZT = TMath::Log(1./zT);
2645
2646       if(fPi0Trigger && decayFound)
2647         FillDecayPhotonCorrelationHistograms(pt, phi, decayMom1,decayMom2,kFALSE) ;
2648
2649       //delta phi cut for correlation
2650       if( (deltaPhi > fDeltaPhiMinCut) && ( deltaPhi < fDeltaPhiMaxCut) )
2651       {
2652         fhDeltaPhiNeutralPt->Fill(pt,deltaPhi);
2653         fhXENeutral        ->Fill(ptTrig,xE);
2654         fhPtHbpXENeutral   ->Fill(ptTrig,hbpXE);
2655       }
2656       else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
2657       {
2658         fhDeltaPhiUeNeutralPt->Fill(pt,deltaPhi);
2659         fhXEUeNeutral        ->Fill(ptTrig,xE);
2660         fhPtHbpXEUeNeutral   ->Fill(ptTrig,hbpXE);
2661       }
2662
2663       //several UE calculation
2664       if(fMakeSeveralUE) FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,xE,hbpXE,zT,hbpZT,deltaPhi);
2665
2666           }
2667     else
2668     {
2669       nrefs++;
2670       if(nrefs==1)
2671       {
2672         refpi0 = new TObjArray(0);
2673         refpi0->SetName(GetAODObjArrayName()+"Pi0s");
2674         refpi0->SetOwner(kFALSE);
2675       }
2676       refpi0->Add(pi0);
2677     }//put references in trigger AOD
2678
2679     if(GetDebug() > 2 )
2680       printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
2681
2682   }//loop
2683
2684   return kTRUE;
2685 }
2686
2687 //_________________________________________________________________________________________________________
2688 void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
2689 {
2690   // Charged Hadron Correlation Analysis with MC information
2691
2692   if(GetDebug()>1)
2693     printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation() - Make trigger particle - charged hadron correlation in AOD MC level\n");
2694
2695   AliStack         * stack        = 0x0 ;
2696   TParticle        * primary      = 0x0 ;
2697   TClonesArray     * mcparticles0 = 0x0 ;
2698   TClonesArray     * mcparticles  = 0x0 ;
2699   AliAODMCParticle * aodprimary   = 0x0 ;
2700
2701   Double_t eprim   = 0 ;
2702   Double_t ptprim  = 0 ;
2703   Double_t phiprim = 0 ;
2704   Double_t etaprim = 0 ;
2705   Int_t    nTracks = 0 ;
2706   Int_t iParticle  = 0 ;
2707   Double_t charge  = 0.;
2708
2709
2710   if(GetReader()->ReadStack())
2711   {
2712     nTracks = GetMCStack()->GetNtrack() ;
2713   }
2714   else
2715   {
2716     nTracks = GetReader()->GetAODMCParticles()->GetEntriesFast() ;
2717   }
2718   //Int_t trackIndex[nTracks];
2719
2720   Int_t label= aodParticle->GetLabel();
2721   if(label < 0)
2722   {
2723     if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** bad label ***:  label %d \n", label);
2724     return;
2725   }
2726
2727
2728   if(GetReader()->ReadStack())
2729   {
2730     stack =  GetMCStack() ;
2731     if(!stack)
2732     {
2733       printf(" AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation- Stack not available, is the MC handler called? STOP\n");
2734       abort();
2735     }
2736
2737     nTracks=stack->GetNprimary();
2738     if(label >=  stack->GetNtrack())
2739     {
2740       if(GetDebug() > 2)  printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
2741       return ;
2742     }
2743
2744     primary = stack->Particle(label);
2745     if(!primary)
2746     {
2747       printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no primary ***:  label %d \n", label);
2748       return;
2749     }
2750
2751     if(primary)
2752     {
2753       eprim    = primary->Energy();
2754       ptprim   = primary->Pt();
2755       phiprim  = primary->Phi();
2756       etaprim  = primary->Eta();
2757
2758       if(ptprim < 0.01 || eprim < 0.01) return ;
2759
2760       for (iParticle = 0 ; iParticle <  nTracks ; iParticle++)
2761       {
2762         TParticle * particle = stack->Particle(iParticle);
2763         TLorentzVector momentum;
2764
2765         //keep only final state particles
2766         if(particle->GetStatusCode()!=1) continue ;
2767
2768         Int_t pdg = particle->GetPdgCode();
2769
2770         charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
2771
2772         particle->Momentum(momentum);
2773
2774         //---------- Charged particles ----------------------
2775         if(charge != 0)
2776         {
2777           //Particles in CTS acceptance
2778           Bool_t inCTS =  GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
2779
2780           if(TMath::Abs(pdg) == 11 && stack->Particle(particle->GetFirstMother())->GetPdgCode()==22) continue ;
2781
2782           if(inCTS)
2783           {
2784             if( label!=iParticle) // avoid trigger particle
2785             {
2786               if(!FillChargedMCCorrelationHistograms(particle->Pt(),particle->Phi(),particle->Eta(),ptprim,phiprim,etaprim)) return;
2787             }
2788           }// in CTS acceptance
2789         }// charged
2790       } //track loop
2791     } //when the leading particles could trace back to MC
2792   } //ESD MC
2793   else if(GetReader()->ReadAODMCParticles())
2794   {
2795     //Get the list of MC particles
2796     mcparticles0 = GetReader()->GetAODMCParticles(0);
2797     if(!mcparticles0) return;
2798     if(label >=mcparticles0->GetEntriesFast())
2799     {
2800       if(GetDebug() > 2)
2801         printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** large label ***:  label %d, n tracks %d \n", label,mcparticles0->GetEntriesFast());
2802       return;
2803     }
2804     //Get the particle
2805     aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
2806     if(!aodprimary)
2807     {
2808       printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no AOD primary ***:  label %d \n", label);
2809       return;
2810     }
2811
2812
2813     if(aodprimary)
2814     {
2815       ptprim  = aodprimary->Pt();
2816       phiprim = aodprimary->Phi();
2817       etaprim = aodprimary->Eta();
2818       eprim   = aodprimary->E();
2819
2820       Bool_t lead = kFALSE;
2821
2822       if(ptprim < 0.01 || eprim < 0.01) return ;
2823
2824       mcparticles= GetReader()->GetAODMCParticles();
2825       for (Int_t i = 0; i < nTracks; i++)
2826       {
2827         AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(i);
2828
2829         if (!part->IsPhysicalPrimary()) continue;
2830         Int_t pdg = part->GetPdgCode();
2831         charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
2832         TLorentzVector momentum(part->Px(),part->Py(),part->Pz(),part->E());
2833         if(charge != 0)
2834         {
2835           if(part->Pt()> GetReader()->GetCTSPtMin())
2836           {
2837             //Particles in CTS acceptance
2838             Bool_t inCTS =  GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
2839             Int_t indexmother=part->GetMother();
2840             if(indexmother>-1)
2841             {
2842               Int_t mPdg = ((AliAODMCParticle*) mcparticles->At(indexmother)) ->GetPdgCode();
2843               if(TMath::Abs(pdg) == 11 && mPdg == 22) continue;
2844             }
2845
2846             if(inCTS)
2847             {
2848               if( label!=iParticle) // avoid trigger particle
2849               {
2850                 if(!FillChargedMCCorrelationHistograms(part->Pt(),part->Phi(),part->Eta(),ptprim,phiprim,etaprim)) return;
2851                 else lead = kTRUE;
2852               }
2853             } // in acceptance
2854           } // min pt cut
2855         } //only charged particles
2856       }  //MC particle loop
2857       if (lead)
2858         {
2859           fhMCPtLeading->Fill(ptprim);
2860           fhMCPhiLeading->Fill(ptprim);
2861           fhMCEtaLeading->Fill(ptprim);
2862         }
2863     } //when the leading particles could trace back to MC
2864   }// AOD MC
2865 }
2866
2867 //_____________________________________________________________________
2868 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
2869 {
2870
2871   //Print some relevant parameters set for the analysis
2872   if(! opt)
2873     return;
2874
2875   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2876   AliAnaCaloTrackCorrBaseClass::Print(" ");
2877   printf("Pt trigger           >  %3.2f\n", fMinTriggerPt) ;
2878   printf("Pt associated hadron <  %3.2f\n", fMaxAssocPt) ;
2879   printf("Pt associated hadron >  %3.2f\n", fMinAssocPt) ;
2880   printf("Phi trigger particle-Hadron   <  %3.2f\n", fDeltaPhiMaxCut) ;
2881   printf("Phi trigger particle-Hadron   >  %3.2f\n", fDeltaPhiMinCut) ;
2882   printf("Phi trigger particle-UeHadron <  %3.2f\n", fUeDeltaPhiMaxCut) ;
2883   printf("Phi trigger particle-UeHadron >  %3.2f\n", fUeDeltaPhiMinCut) ;
2884   printf("Isolated Trigger?  %d\n"     , fSelectIsolated) ;
2885   printf("Several UE?  %d\n"           , fMakeSeveralUE) ;
2886   printf("Name of AOD Pi0 Branch %s \n", fPi0AODBranchName.Data());
2887   printf("Do Decay-hadron correlation ?  %d\n", fPi0Trigger) ;
2888   printf("Select absolute leading for cluster triggers ?  %d or Near Side %d\n",
2889          fMakeAbsoluteLeading, fMakeNearSideLeading) ;
2890   printf("Trigger pt bins  %d\n", fNAssocPtBins) ;
2891   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
2892     printf("\t bin %d = [%2.1f,%2.1f]\n", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
2893   }
2894
2895 }
2896
2897 //____________________________________________________________
2898 void AliAnaParticleHadronCorrelation::SetNAssocPtBins(Int_t n)
2899 {
2900   // Set number of bins
2901
2902   fNAssocPtBins  = n ;
2903
2904
2905   if(n < 10 && n > 0)
2906   {
2907     fNAssocPtBins  = n ;
2908   }
2909   else
2910   {
2911     printf("n = larger than 9 or too small, set to 9 \n");
2912     fNAssocPtBins = 9;
2913   }
2914 }
2915
2916 //______________________________________________________________________________
2917 void AliAnaParticleHadronCorrelation::SetAssocPtBinLimit(Int_t ibin, Float_t pt)
2918 {
2919   // Set the list of limits for the trigger pt bins
2920
2921   if(ibin <= fNAssocPtBins || ibin >= 0)
2922   {
2923     fAssocPtBinLimit[ibin] = pt  ;
2924   }
2925   else
2926   {
2927     printf("AliAnaParticleHadronCorrelation::SetAssocPtBinLimit() - bin  number too large %d > %d or small, nothing done\n", ibin, fNAssocPtBins) ;
2928   }
2929 }
2930