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