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->IsEventSelected( ) & GetReader()->GetEventTriggerMask())
 548   {
 549     fhNtracksTrigger->Fill(nTracks);
 550   }
 551
 552   if(inputHandler->IsEventSelected( ) & AliVEvent::kAnyINT)
 553   {
 554
 555     fhNtracksINT->Fill(nTracks);
 556
 557     //Get vertex z bin
 558     Double_t v[3] = {0,0,0}; //vertex ;
 559     GetReader()->GetVertex(v);
 560
 561     Int_t curZvertBin = (Int_t)(0.5*GetNZvertBin()*(v[2]+GetZvertexCut())/GetZvertexCut()) ;
 562
 563     // centrality or tracks bin
 564     Int_t curCentrBin = 0;
 565     if(fUseTrackMultBins)
 566     { // Track multiplicity bins
 567       //curCentrBin = (GetTrackMultiplicity()-1)/5;
 568       //if(curCentrBin > GetNCentrBin()-1) curCentrBin=GetNCentrBin()-1;
 569       Int_t trackMult = GetReader()->GetTrackMultiplicity();
 570       if(trackMult<=5)
 571         curCentrBin=8;
 572       else if(trackMult<=10)
 573         curCentrBin=7;
 574       else if(trackMult<=15)
 575         curCentrBin=6;
 576       else if(trackMult<=20)
 577         curCentrBin=5;
 578       else if(trackMult<=30)
 579         curCentrBin=4;
 580       else if(trackMult<=40)
 581         curCentrBin=3;
 582       else if(trackMult<=55)
 583         curCentrBin=2;
 584       else if(trackMult<=70)
 585         curCentrBin=1 ;
 586       else curCentrBin=0 ;
 587     }
 588     else // Set centrality based on centrality task
 589     {
 590       curCentrBin = GetEventCentrality() * GetNCentrBin() / GetReader()->GetCentralityOpt();
 591       if(GetDebug() > 0 )printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - curCentrBin %d, centrality %d, n bins %d, max bin from centrality %d\n",
 592                                 curCentrBin, GetEventCentrality(), GetNCentrBin(), GetReader()->GetCentralityOpt());
 593     }
 594
 595     TObjArray * mixEventTracks = new TObjArray;
 596
 597     //printf("curCen %d, curZ %d, bin %d\n",curCentrBin,curZvertBin,curCentrBin*GetNZvertBin()+curZvertBin);
 598
 599     if(!fListMixEvents[curCentrBin*GetNZvertBin()+curZvertBin])
 600       fListMixEvents[curCentrBin*GetNZvertBin()+curZvertBin]=new TList();
 601
 602     TList * pool = fListMixEvents[curCentrBin*GetNZvertBin()+curZvertBin];
 603
 604     TVector3 p3;
 605     for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
 606     {
 607       AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
 608
 609       Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 610       p3.SetXYZ(mom[0],mom[1],mom[2]);
 611       Float_t pt   = p3.Pt();
 612
 613       //Select only hadrons in pt range
 614       if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
 615
 616       AliAODPWG4Particle * mixedTrack = new AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
 617       mixedTrack->SetDetector("CTS");
 618       mixedTrack->SetChargedBit(track->Charge()>0);
 619
 620       mixEventTracks->Add(mixedTrack);
 621     }
 622
 623     pool->AddFirst(mixEventTracks);
 624     mixEventTracks = 0;
 625     if(pool->GetSize() >= GetNMaxEvMix())
 626     {//Remove last event
 627       TClonesArray * tmp = static_cast<TClonesArray*>(pool->Last()) ;
 628       pool->RemoveLast() ;
 629       delete tmp ;
 630     }
 631
 632   } // MB event
 633
 634 }
 635
 636 //____________________________________________________________
 637 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
 638 {
 639   //Save parameters used for analysis
 640   TString parList ; //this will be list of parameters used for this analysis.
 641   const Int_t buffersize = 560;
 642   char onePar[buffersize] ;
 643
 644   snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;
 645   parList+=onePar ;
 646   snprintf(onePar,buffersize," Pt Trigger > %3.2f ",    fMinTriggerPt) ;
 647   parList+=onePar ;
 648   snprintf(onePar,buffersize," %3.2f < Pt associated < %3.2f ", fMinAssocPt,   fMaxAssocPt) ;
 649   parList+=onePar ;
 650   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-Hadron < %3.2f ",    fDeltaPhiMinCut,   fDeltaPhiMaxCut) ;
 651   parList+=onePar ;
 652   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-UeHadron <  %3.2f ", fUeDeltaPhiMinCut, fUeDeltaPhiMaxCut) ;
 653   parList+=onePar ;
 654   snprintf(onePar,buffersize,"Isolated Trigger?  %d\n",    fSelectIsolated) ;
 655   parList+=onePar ;
 656   snprintf(onePar,buffersize,"Several UE?  %d\n",          fMakeSeveralUE) ;
 657   parList+=onePar ;
 658   snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ", fPi0AODBranchName.Data());
 659   parList+=onePar ;
 660   snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  pi0 %d, decay %d", fPi0Trigger, fDecayTrigger) ;
 661   parList+=onePar ;
 662   snprintf(onePar,buffersize,"Select absolute leading for cluster triggers ? %d or Near Side Leading %d \n",
 663            fMakeAbsoluteLeading, fMakeNearSideLeading) ;
 664   parList+=onePar ;
 665   snprintf(onePar,buffersize,"Associated particle pt bins  %d: ", fNAssocPtBins) ;
 666   parList+=onePar ;
 667   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
 668     snprintf(onePar,buffersize,"bin %d = [%2.1f,%2.1f];", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
 669   }
 670   parList+=onePar ;
 671
 672   //Get parameters set in base class.
 673   parList += GetBaseParametersList() ;
 674
 675   //Get parameters set in FiducialCut class (not available yet)
 676   //parlist += GetFidCut()->GetFidCutParametersList()
 677
 678   return new TObjString(parList) ;
 679
 680 }
 681
 682 //________________________________________________________________
 683 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
 684 {
 685
 686   // Create histograms to be saved in output file and
 687   // store them in fOutputContainer
 688
 689
 690   TList * outputContainer = new TList() ;
 691   outputContainer->SetName("CorrelationHistos") ;
 692
 693   Int_t   nptbins = GetHistogramRanges()->GetHistoPtBins(); Int_t  nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t   netabins = GetHistogramRanges()->GetHistoEtaBins();
 694   Float_t ptmax   = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax  = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax   = GetHistogramRanges()->GetHistoEtaMax();
 695   Float_t ptmin   = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin  = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin   = GetHistogramRanges()->GetHistoEtaMin();
 696
 697   fhPtLeading  = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax);
 698   fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
 699
 700   fhPhiLeading  = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 701   fhPhiLeading->SetYTitle("#phi (rad)");
 702
 703   fhEtaLeading  = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 704   fhEtaLeading->SetYTitle("#eta ");
 705
 706   outputContainer->Add(fhPtLeading);
 707   outputContainer->Add(fhPhiLeading);
 708   outputContainer->Add(fhEtaLeading);
 709
 710   //Correlation with charged hadrons
 711   if(GetReader()->IsCTSSwitchedOn())
 712   {
 713     fhDeltaPhiDeltaEtaCharged  = new TH2F
 714     ("hDeltaPhiDeltaEtaCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
 715      140,-2.,5.,200,-2,2);
 716     fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
 717     fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
 718
 719     fhPhiCharged  = new TH2F
 720     ("hPhiCharged","#phi_{h^{#pm}}  vs p_{T #pm}",
 721      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 722     fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
 723     fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
 724
 725     fhEtaCharged  = new TH2F
 726     ("hEtaCharged","#eta_{h^{#pm}}  vs p_{T #pm}",
 727      nptbins,ptmin,ptmax,netabins,etamin,etamax);
 728     fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
 729     fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
 730
 731     fhDeltaPhiCharged  = new TH2F
 732     ("hDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
 733      nptbins,ptmin,ptmax,140,-2.,5.);
 734     fhDeltaPhiCharged->SetYTitle("#Delta #phi");
 735     fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
 736
 737     fhDeltaPhiChargedPt  = new TH2F
 738     ("hDeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
 739      nptbins,ptmin,ptmax,140,-2.,5.);
 740     fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
 741     fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 742
 743     fhDeltaPhiUeChargedPt  = new TH2F
 744     ("hDeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
 745      nptbins,ptmin,ptmax,140,-2.,5.);
 746     fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
 747     fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 748
 749     fhDeltaEtaCharged  = new TH2F
 750     ("hDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
 751      nptbins,ptmin,ptmax,200,-2,2);
 752     fhDeltaEtaCharged->SetYTitle("#Delta #eta");
 753     fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
 754
 755     fhXECharged  =
 756     new TH2F("hXECharged","x_{E} for charged tracks",
 757              nptbins,ptmin,ptmax,200,0.,2.);
 758     fhXECharged->SetYTitle("x_{E}");
 759     fhXECharged->SetXTitle("p_{T trigger}");
 760
 761     fhXEUeCharged  =
 762     new TH2F("hXEUeCharged","x_{E} for Underlying Event",
 763              nptbins,ptmin,ptmax,200,0.,2.);
 764     fhXEUeCharged->SetYTitle("x_{E}");
 765     fhXEUeCharged->SetXTitle("p_{T trigger}");
 766
 767     fhXEPosCharged  =
 768     new TH2F("hXEPositiveCharged","x_{E} for positive charged tracks",
 769              nptbins,ptmin,ptmax,200,0.,2.);
 770     fhXEPosCharged->SetYTitle("x_{E}");
 771     fhXEPosCharged->SetXTitle("p_{T trigger}");
 772
 773     fhXENegCharged  =
 774     new TH2F("hXENegativeCharged","x_{E} for negative charged tracks",
 775              nptbins,ptmin,ptmax,200,0.,2.);
 776     fhXENegCharged->SetYTitle("x_{E}");
 777     fhXENegCharged->SetXTitle("p_{T trigger}");
 778
 779     fhPtHbpXECharged  =
 780     new TH2F("hHbpXECharged","#xi = ln(1/x_{E}) with charged hadrons",
 781              nptbins,ptmin,ptmax,200,0.,10.);
 782     fhPtHbpXECharged->SetYTitle("ln(1/x_{E})");
 783     fhPtHbpXECharged->SetXTitle("p_{T trigger}");
 784
 785     fhPtHbpXEUeCharged  =
 786     new TH2F("hHbpXEUeCharged","#xi = ln(1/x_{E}) with charged hadrons,Underlying Event",
 787              nptbins,ptmin,ptmax,200,0.,10.);
 788     fhPtHbpXEUeCharged->SetYTitle("ln(1/x_{E})");
 789     fhPtHbpXEUeCharged->SetXTitle("p_{T trigger}");
 790
 791     fhZTCharged  =
 792     new TH2F("hZTCharged","z_{T} for charged tracks",
 793              nptbins,ptmin,ptmax,200,0.,2.);
 794     fhZTCharged->SetYTitle("z_{T}");
 795     fhZTCharged->SetXTitle("p_{T trigger}");
 796
 797     fhZTUeCharged  =
 798     new TH2F("hZTUeCharged","z_{T} for Underlying Event",
 799              nptbins,ptmin,ptmax,200,0.,2.);
 800     fhZTUeCharged->SetYTitle("z_{T}");
 801     fhZTUeCharged->SetXTitle("p_{T trigger}");
 802
 803     fhZTPosCharged  =
 804     new TH2F("hZTPositiveCharged","z_{T} for positive charged tracks",
 805              nptbins,ptmin,ptmax,200,0.,2.);
 806     fhZTPosCharged->SetYTitle("z_{T}");
 807     fhZTPosCharged->SetXTitle("p_{T trigger}");
 808
 809     fhZTNegCharged  =
 810     new TH2F("hZTNegativeCharged","z_{T} for negative charged tracks",
 811              nptbins,ptmin,ptmax,200,0.,2.);
 812     fhZTNegCharged->SetYTitle("z_{T}");
 813     fhZTNegCharged->SetXTitle("p_{T trigger}");
 814
 815     fhPtHbpZTCharged  =
 816     new TH2F("hHbpZTCharged","#xi = ln(1/z_{T}) with charged hadrons",
 817              nptbins,ptmin,ptmax,200,0.,10.);
 818     fhPtHbpZTCharged->SetYTitle("ln(1/z_{T})");
 819     fhPtHbpZTCharged->SetXTitle("p_{T trigger}");
 820
 821     fhPtHbpZTUeCharged  =
 822     new TH2F("hHbpZTUeCharged","#xi = ln(1/z_{T}) with charged hadrons,Underlying Event",
 823              nptbins,ptmin,ptmax,200,0.,10.);
 824     fhPtHbpZTUeCharged->SetYTitle("ln(1/x_{E})");
 825     fhPtHbpZTUeCharged->SetXTitle("p_{T trigger}");
 826
 827     fhPtTrigPout  =
 828     new TH2F("hPtTrigPout","Pout with triggers",
 829              nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
 830     fhPtTrigPout->SetYTitle("p_{out} (GeV/c)");
 831     fhPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
 832
 833     fhPtTrigCharged  =
 834     new TH2F("hPtTrigCharged","trgger and charged tracks pt distribution",
 835              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 836     fhPtTrigCharged->SetYTitle("p_{T h^{#pm}} (GeV/c)");
 837     fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)");
 838
 839     outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
 840     outputContainer->Add(fhPhiCharged) ;
 841     outputContainer->Add(fhEtaCharged) ;
 842     outputContainer->Add(fhDeltaPhiCharged) ;
 843     outputContainer->Add(fhDeltaEtaCharged) ;
 844     outputContainer->Add(fhDeltaPhiChargedPt) ;
 845     outputContainer->Add(fhDeltaPhiUeChargedPt) ;
 846
 847     outputContainer->Add(fhXECharged) ;
 848     outputContainer->Add(fhXEPosCharged) ;
 849     outputContainer->Add(fhXENegCharged) ;
 850     outputContainer->Add(fhXEUeCharged) ;
 851     outputContainer->Add(fhPtHbpXECharged) ;
 852     outputContainer->Add(fhPtHbpXEUeCharged) ;
 853
 854     outputContainer->Add(fhZTCharged) ;
 855     outputContainer->Add(fhZTPosCharged) ;
 856     outputContainer->Add(fhZTNegCharged) ;
 857     outputContainer->Add(fhZTUeCharged) ;
 858     outputContainer->Add(fhPtHbpZTCharged) ;
 859     outputContainer->Add(fhPtHbpZTUeCharged) ;
 860
 861     outputContainer->Add(fhPtTrigPout) ;
 862     outputContainer->Add(fhPtTrigCharged) ;
 863
 864     if(DoEventSelect())
 865     {
 866       Int_t nMultiBins = GetMultiBin();
 867       fhTrigDeltaPhiCharged = new TH2F*[nMultiBins] ;
 868       fhTrigDeltaEtaCharged = new TH2F*[nMultiBins] ;
 869       fhTrigXECorr          = new TH2F*[nMultiBins] ;
 870       fhTrigXEUeCorr        = new TH2F*[nMultiBins] ;
 871       fhTrigZTCorr          = new TH2F*[nMultiBins] ;
 872       fhTrigZTUeCorr        = new TH2F*[nMultiBins] ;
 873
 874       for(Int_t im=0; im<nMultiBins; im++)
 875       {
 876         fhTrigDeltaPhiCharged[im]  = new TH2F
 877         (Form("hTrigDeltaPhiCharged_%d",im),Form("hTrigDeltaPhiCharged_%d",im), nptbins,ptmin,ptmax, 140,-2.,5.);
 878         fhTrigDeltaPhiCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
 879         fhTrigDeltaPhiCharged[im]->SetYTitle("#Delta #phi");
 880
 881         fhTrigDeltaEtaCharged[im]  = new TH2F
 882         (Form("hTrigDeltaEtaCharged_%d",im),Form("hTrigDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, 200,-2,2);
 883         fhTrigDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
 884         fhTrigDeltaEtaCharged[im]->SetYTitle("#Delta #eta");
 885
 886         fhTrigXECorr[im]  = new TH2F
 887         (Form("hTrigXEPtCorr_%d",im),Form("hTrigXEPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 888         fhTrigXECorr[im]->SetYTitle("x_{E trigger h^{#pm}}");
 889         fhTrigXECorr[im]->SetXTitle("p_{T trigger}");
 890
 891         fhTrigXEUeCorr[im]  = new TH2F
 892         (Form("hTrigXEPtUeCorr_%d",im),Form("hTrigXEPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 893         fhTrigXEUeCorr[im]->SetYTitle("x_{E trigger h^{#pm}}");
 894         fhTrigXEUeCorr[im]->SetXTitle("p_{T trigger}");
 895
 896         fhTrigZTCorr[im]  = new TH2F
 897         (Form("hTrigZTPtCorr_%d",im),Form("hTrigZTPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 898         fhTrigZTCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
 899         fhTrigZTCorr[im]->SetXTitle("p_{T trigger}");
 900
 901         fhTrigZTUeCorr[im]  = new TH2F
 902         (Form("hTrigZTPtUeCorr_%d",im),Form("hTrigZTPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 903         fhTrigZTUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
 904         fhTrigZTUeCorr[im]->SetXTitle("p_{T trigger}");
 905
 906         outputContainer->Add(fhTrigDeltaPhiCharged[im]) ;
 907         outputContainer->Add(fhTrigDeltaEtaCharged[im]) ;
 908         outputContainer->Add(fhTrigXECorr[im]);
 909         outputContainer->Add(fhTrigXEUeCorr[im]);
 910         outputContainer->Add(fhTrigZTCorr[im]);
 911         outputContainer->Add(fhTrigZTUeCorr[im]);
 912       }
 913     }
 914
 915     if(fFillBradHisto)
 916     {
 917       fhAssocPtBkg        = new TH2F("hAssocPtBkg", " Trigger p_{T} vs associated hadron p_{T} from background",
 918                                      nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
 919       fhAssocPtBkg->SetXTitle("p_{T trigger}");
 920       fhAssocPtBkg->SetYTitle("p_{T associated}");
 921       outputContainer->Add(fhAssocPtBkg) ;
 922
 923       fhDeltaPhiBrad = new TH2F("hDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs p_{T trigger} ",
 924                                 nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
 925       fhDeltaPhiBrad->SetXTitle("p_{T trigger}");
 926       fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
 927       outputContainer->Add(fhDeltaPhiBrad) ;
 928     }
 929
 930     fhDeltaPhiAssocPtBin     = new TH2F*[fNAssocPtBins] ;
 931     fhXEAssocPtBin           = new TH2F*[fNAssocPtBins] ;
 932     fhZTAssocPtBin           = new TH2F*[fNAssocPtBins] ;
 933
 934     if(fFillBradHisto)
 935       fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins] ;
 936
 937     if(fPi0Trigger || fDecayTrigger)
 938     {
 939       fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins] ;
 940       fhXEAssocPtBin             = new TH2F*[fNAssocPtBins] ;
 941       fhZTAssocPtBin             = new TH2F*[fNAssocPtBins] ;
 942       fhXEDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
 943       fhZTDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
 944       fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
 945     }
 946
 947     if(fHMPIDCorrelation)
 948     {
 949       fhDeltaPhiAssocPtBinHMPID   = new TH2F*[fNAssocPtBins] ;
 950       fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins] ;
 951     }
 952
 953     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
 954     {
 955       fhDeltaPhiAssocPtBin[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 956                                          Form("#Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 957                                          nptbins, ptmin, ptmax,140,-2.,5.);
 958       fhDeltaPhiAssocPtBin[i]->SetXTitle("p_{T trigger}");
 959       fhDeltaPhiAssocPtBin[i]->SetYTitle("#Delta #phi");
 960
 961       fhXEAssocPtBin[i]       = new TH2F(Form("hXEAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 962                                          Form("x_{E} vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 963                                          nptbins, ptmin, ptmax,200, 0.0, 2.0);
 964       fhXEAssocPtBin[i]->SetXTitle("p_{T trigger}");
 965       fhXEAssocPtBin[i]->SetYTitle("x_{E}");
 966
 967       fhZTAssocPtBin[i]       = new TH2F(Form("hZTAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 968                                          Form("z_{T} vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 969                                          nptbins, ptmin, ptmax,200, 0.0, 2.0);
 970       fhZTAssocPtBin[i]->SetXTitle("p_{T trigger}");
 971       fhZTAssocPtBin[i]->SetYTitle("z_{T}");
 972
 973       outputContainer->Add(fhDeltaPhiAssocPtBin[i]) ;
 974       outputContainer->Add(fhXEAssocPtBin[i]);
 975       outputContainer->Add(fhZTAssocPtBin[i]);
 976
 977       if(fPi0Trigger || fDecayTrigger)
 978       {
 979         fhDeltaPhiDecayChargedAssocPtBin[i] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 980                                            Form("#Delta #phi vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 981                                            nptbins, ptmin, ptmax,140,-2.,5.);
 982         fhDeltaPhiDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
 983         fhDeltaPhiDecayChargedAssocPtBin[i]->SetYTitle("#Delta #phi");
 984
 985         fhXEDecayChargedAssocPtBin[i]       = new TH2F(Form("hXEDecayChargedAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 986                                            Form("x_{E} vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 987                                            nptbins, ptmin, ptmax,200, 0.0, 2.0);
 988         fhXEDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
 989         fhXEDecayChargedAssocPtBin[i]->SetYTitle("x_{E}");
 990
 991         fhZTDecayChargedAssocPtBin[i]       = new TH2F(Form("hZTDecayChargedAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 992                                            Form("z_{T} vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
 993                                            nptbins, ptmin, ptmax,200, 0.0, 2.0);
 994         fhZTDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
 995         fhZTDecayChargedAssocPtBin[i]->SetYTitle("z_{T}");
 996
 997         outputContainer->Add(fhDeltaPhiDecayChargedAssocPtBin[i]) ;
 998         outputContainer->Add(fhXEDecayChargedAssocPtBin[i]);
 999         outputContainer->Add(fhZTDecayChargedAssocPtBin[i]);
1000
1001       }
1002
1003       if(fFillBradHisto)
1004       {
1005         fhDeltaPhiBradAssocPtBin[i] = new TH2F(Form("hDeltaPhiBradPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1006                                                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]),
1007                                                nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
1008         fhDeltaPhiBradAssocPtBin[i]->SetXTitle("p_{T trigger}");
1009         fhDeltaPhiBradAssocPtBin[i]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
1010         outputContainer->Add(fhDeltaPhiBradAssocPtBin[i]) ;
1011       }
1012
1013       if(fHMPIDCorrelation)
1014       {
1015         fhDeltaPhiAssocPtBinHMPID[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1fHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1016                                                 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]),
1017                                                 nptbins, ptmin, ptmax,140,-2.,5.);
1018         fhDeltaPhiAssocPtBinHMPID[i]->SetXTitle("p_{T trigger}");
1019         fhDeltaPhiAssocPtBinHMPID[i]->SetYTitle("#Delta #phi");
1020
1021         fhDeltaPhiAssocPtBinHMPIDAcc[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1fHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1022                                                    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]),
1023                                                    nptbins, ptmin, ptmax,140,-2.,5.);
1024         fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetXTitle("p_{T trigger}");
1025         fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetYTitle("#Delta #phi");
1026
1027         outputContainer->Add(fhDeltaPhiAssocPtBinHMPID[i]) ;
1028         outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[i]) ;
1029
1030       }
1031     }
1032
1033     if(fPi0Trigger || fDecayTrigger)
1034     {
1035       if(fPi0Trigger)
1036       {
1037         fhPtPi0DecayRatio  = new TH2F
1038         ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay",
1039          nptbins,ptmin,ptmax, 100,0.,2.);
1040         fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");
1041         fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");
1042       }
1043
1044       fhDeltaPhiDecayCharged  = new TH2F
1045       ("hDeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",
1046        nptbins,ptmin,ptmax,140,-2.,5.);
1047       fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");
1048       fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");
1049
1050       fhXEDecayCharged  =
1051       new TH2F("hXEDecayCharged","x_{E}  Decay",
1052                nptbins,ptmin,ptmax,200,0.,2.);
1053       fhXEDecayCharged->SetYTitle("x_{E}");
1054       fhXEDecayCharged->SetXTitle("p_{T decay}");
1055
1056       fhZTDecayCharged  =
1057       new TH2F("hZTDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",
1058                nptbins,ptmin,ptmax,200,0.,2.);
1059       fhZTDecayCharged->SetYTitle("z_{decay h^{#pm}}");
1060       fhZTDecayCharged->SetXTitle("p_{T decay}");
1061
1062       outputContainer->Add(fhPtPi0DecayRatio) ;
1063       outputContainer->Add(fhDeltaPhiDecayCharged) ;
1064       outputContainer->Add(fhXEDecayCharged) ;
1065       outputContainer->Add(fhZTDecayCharged) ;
1066     }
1067
1068     if(fMakeSeveralUE)
1069     {
1070       fhDeltaPhiUeLeftCharged  = new TH2F
1071       ("hDeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
1072        nptbins,ptmin,ptmax,140,-2.,5.);
1073       fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
1074       fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
1075       outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
1076
1077       fhDeltaPhiUeRightCharged  = new TH2F
1078       ("hDeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
1079        nptbins,ptmin,ptmax,140,-2.,5.);
1080       fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
1081       fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
1082       outputContainer->Add(fhDeltaPhiUeRightCharged) ;
1083
1084       fhXEUeLeftCharged  =
1085       new TH2F("hXEUeChargedLeft","x_{E} with UE left side of trigger",
1086                nptbins,ptmin,ptmax,200,0.,2.);
1087       fhXEUeLeftCharged->SetYTitle("x_{E Ueh^{#pm}}");
1088       fhXEUeLeftCharged->SetXTitle("p_{T trigger}");
1089       outputContainer->Add(fhXEUeLeftCharged) ;
1090
1091       fhXEUeRightCharged  =
1092       new TH2F("hXEUeChargedRight","x_{E h^{#pm}} with UE right side of trigger",
1093                nptbins,ptmin,ptmax,200,0.,2.);
1094       fhXEUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
1095       fhXEUeRightCharged->SetXTitle("p_{T trigger}");
1096       outputContainer->Add(fhXEUeRightCharged) ;
1097
1098       fhPtHbpXEUeLeftCharged  =
1099       new TH2F("hHbpXEUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",
1100                nptbins,ptmin,ptmax,200,0.,10.);
1101       fhPtHbpXEUeLeftCharged->SetYTitle("ln(1/x_{E})");
1102       fhPtHbpXEUeLeftCharged->SetXTitle("p_{T trigger}");
1103       outputContainer->Add(fhPtHbpXEUeLeftCharged) ;
1104
1105       fhPtHbpXEUeRightCharged  =
1106       new TH2F("hHbpXEUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",
1107                nptbins,ptmin,ptmax,200,0.,10.);
1108       fhPtHbpXEUeRightCharged->SetYTitle("ln(1/x_{E})");
1109       fhPtHbpXEUeRightCharged->SetXTitle("p_{T trigger}");
1110       outputContainer->Add(fhPtHbpXEUeRightCharged) ;
1111
1112       fhZTUeLeftCharged  =
1113       new TH2F("hZTUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",
1114                nptbins,ptmin,ptmax,200,0.,2.);
1115       fhZTUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
1116       fhZTUeLeftCharged->SetXTitle("p_{T trigger}");
1117       outputContainer->Add(fhZTUeLeftCharged) ;
1118
1119       fhZTUeRightCharged  =
1120       new TH2F("hZTUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",
1121                nptbins,ptmin,ptmax,200,0.,2.);
1122       fhZTUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
1123       fhZTUeRightCharged->SetXTitle("p_{T trigger}");
1124       outputContainer->Add(fhZTUeRightCharged) ;
1125
1126       fhPtHbpZTUeLeftCharged  =
1127       new TH2F("hHbpZTUeChargedLeft","#xi = ln(1/z_{T}) with charged UE left side of trigger",
1128                nptbins,ptmin,ptmax,200,0.,10.);
1129       fhPtHbpZTUeLeftCharged->SetYTitle("ln(1/z_{T})");
1130       fhPtHbpZTUeLeftCharged->SetXTitle("p_{T trigger}");
1131       outputContainer->Add(fhPtHbpZTUeLeftCharged) ;
1132
1133       fhPtHbpZTUeRightCharged  =
1134       new TH2F("hHbpZTUeChargedRight","#xi = ln(1/z_{T}) with charged UE right side of trigger",
1135                nptbins,ptmin,ptmax,200,0.,10.);
1136       fhPtHbpZTUeRightCharged->SetYTitle("ln(1/z_{T})");
1137       fhPtHbpZTUeRightCharged->SetXTitle("p_{T trigger}");
1138       outputContainer->Add(fhPtHbpZTUeRightCharged) ;
1139
1140     }
1141   }  //Correlation with charged hadrons
1142
1143   //Correlation with neutral hadrons
1144   if(fNeutralCorr)
1145   {
1146     fhDeltaPhiDeltaEtaNeutral  = new TH2F
1147     ("hDeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
1148      140,-2.,5.,200,-2,2);
1149     fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");
1150     fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
1151
1152     fhPhiNeutral  = new TH2F
1153     ("hPhiNeutral","#phi_{#pi^{0}}  vs p_{T #pi^{0}}",
1154      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1155     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
1156     fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
1157
1158     fhEtaNeutral  = new TH2F
1159     ("hEtaNeutral","#eta_{#pi^{0}}  vs p_{T #pi^{0}}",
1160      nptbins,ptmin,ptmax,netabins,etamin,etamax);
1161     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
1162     fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
1163
1164     fhDeltaPhiNeutral  = new TH2F
1165     ("hDeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
1166      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1167     fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
1168     fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
1169
1170     fhDeltaPhiNeutralPt  = new TH2F
1171     ("hDeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
1172      nptbins,ptmin,ptmax,140,-2.,5.);
1173     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
1174     fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
1175
1176     fhDeltaPhiUeNeutralPt  = new TH2F
1177     ("hDeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
1178      nptbins,ptmin,ptmax,140,-2.,5.);
1179     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
1180     fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
1181
1182     fhDeltaEtaNeutral  = new TH2F
1183     ("hDeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
1184      nptbins,ptmin,ptmax,200,-2,2);
1185     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
1186     fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
1187
1188     fhXENeutral  =
1189     new TH2F("hXENeutral","x_{E} for #pi^{0} associated",
1190              nptbins,ptmin,ptmax,200,0.,2.);
1191     fhXENeutral->SetYTitle("x_{E}");
1192     fhXENeutral->SetXTitle("p_{T trigger}");
1193
1194     fhXEUeNeutral  =
1195     new TH2F("hXEUeNeutral","x_{E} for #pi^{0} associated",
1196              nptbins,ptmin,ptmax,200,0.,2.);
1197     fhXEUeNeutral->SetYTitle("x_{E}");
1198     fhXEUeNeutral->SetXTitle("p_{T trigger}");
1199
1200     fhPtHbpXENeutral  =
1201     new TH2F("hHbpXENeutral","#xi = ln(1/x_{E})for #pi^{0} associated",
1202              nptbins,ptmin,ptmax,200,0.,10.);
1203     fhPtHbpXENeutral->SetYTitle("ln(1/x_{E})");
1204     fhPtHbpXENeutral->SetXTitle("p_{T trigger}");
1205
1206     fhPtHbpXEUeNeutral  =
1207     new TH2F("hHbpXEUeNeutral","#xi = ln(1/x_{E}) for #pi^{0} associated",
1208              nptbins,ptmin,ptmax,200,0.,10.);
1209     fhPtHbpXEUeNeutral->SetYTitle("ln(1/x_{E})");
1210     fhPtHbpXEUeNeutral->SetXTitle("p_{T trigger}");
1211
1212     fhZTNeutral  =
1213     new TH2F("hZTNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger} for #pi^{0} associated",
1214              nptbins,ptmin,ptmax,200,0.,2.);
1215     fhZTNeutral->SetYTitle("z_{trigger #pi^{0}}");
1216     fhZTNeutral->SetXTitle("p_{T trigger}");
1217
1218     fhZTUeNeutral  =
1219     new TH2F("hZTUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger} for #pi^{0} associated",
1220              nptbins,ptmin,ptmax,200,0.,2.);
1221     fhZTUeNeutral->SetYTitle("z_{trigger #pi^{0}}");
1222     fhZTUeNeutral->SetXTitle("p_{T trigger}");
1223
1224     fhPtHbpZTNeutral  =
1225     new TH2F("hHbpZTNeutral","#xi = ln(1/x_{E}) for #pi^{0} associated",
1226              nptbins,ptmin,ptmax,200,0.,10.);
1227     fhPtHbpZTNeutral->SetYTitle("ln(1/z_{T})");
1228     fhPtHbpZTNeutral->SetXTitle("p_{T trigger}");
1229
1230     fhPtHbpZTUeNeutral  =
1231     new TH2F("hHbpZTUeNeutral","#xi = ln(1/x_{E}) for #pi^{0} associated",
1232              nptbins,ptmin,ptmax,200,0.,10.);
1233     fhPtHbpXEUeNeutral->SetYTitle("ln(1/z_{T})");
1234     fhPtHbpXEUeNeutral->SetXTitle("p_{T trigger}");
1235
1236     outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
1237     outputContainer->Add(fhPhiNeutral) ;
1238     outputContainer->Add(fhEtaNeutral) ;
1239     outputContainer->Add(fhDeltaPhiNeutral) ;
1240     outputContainer->Add(fhDeltaPhiNeutralPt) ;
1241     outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
1242     outputContainer->Add(fhDeltaEtaNeutral) ;
1243     outputContainer->Add(fhXENeutral) ;
1244     outputContainer->Add(fhXEUeNeutral) ;
1245     outputContainer->Add(fhPtHbpXENeutral) ;
1246     outputContainer->Add(fhPtHbpXEUeNeutral) ;
1247     outputContainer->Add(fhZTNeutral) ;
1248     outputContainer->Add(fhZTUeNeutral) ;
1249     outputContainer->Add(fhPtHbpZTNeutral) ;
1250     outputContainer->Add(fhPtHbpZTUeNeutral) ;
1251
1252     if(fPi0Trigger || fDecayTrigger)
1253     {
1254       fhDeltaPhiDecayNeutral  = new TH2F
1255       ("hDeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",
1256        nptbins,ptmin,ptmax,140,-2.,5.);
1257       fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");
1258       fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");
1259
1260       fhXEDecayNeutral  =
1261       new TH2F("hXEDecayNeutral","x_{E} for decay trigger",
1262                nptbins,ptmin,ptmax,200,0.,2.);
1263       fhXEDecayNeutral->SetYTitle("x_{E}");
1264       fhXEDecayNeutral->SetXTitle("p_{T decay}");
1265
1266       fhZTDecayNeutral  =
1267       new TH2F("hZTDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",
1268                nptbins,ptmin,ptmax,200,0.,2.);
1269       fhZTDecayNeutral->SetYTitle("z_{h^{0}}");
1270       fhZTDecayNeutral->SetXTitle("p_{T decay}");
1271
1272       outputContainer->Add(fhDeltaPhiDecayNeutral) ;
1273       outputContainer->Add(fhXEDecayNeutral) ;
1274       outputContainer->Add(fhZTDecayNeutral) ;
1275
1276     }
1277
1278     if(fMakeSeveralUE)
1279     {
1280       fhDeltaPhiUeLeftNeutral  = new TH2F
1281       ("hDeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",
1282        nptbins,ptmin,ptmax,140,-2.,5.);
1283       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
1284       fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
1285       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
1286
1287       fhDeltaPhiUeRightNeutral  = new TH2F
1288       ("hDeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",
1289        nptbins,ptmin,ptmax,140,-2.,5.);
1290       fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
1291       fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
1292       outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
1293
1294       fhXEUeLeftNeutral  =
1295       new TH2F("hXEUeNeutralLeft","x_{E} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
1296                nptbins,ptmin,ptmax,140,0.,2.);
1297       fhXEUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1298       fhXEUeLeftNeutral->SetXTitle("p_{T trigger}");
1299       outputContainer->Add(fhXEUeLeftNeutral) ;
1300
1301       fhXEUeRightNeutral  =
1302       new TH2F("hXEUeNeutralRight","x_{E} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
1303                nptbins,ptmin,ptmax,200,0.,2.);
1304       fhXEUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1305       fhXEUeRightNeutral->SetXTitle("p_{T trigger}");
1306       outputContainer->Add(fhXEUeRightNeutral) ;
1307
1308       fhPtHbpXEUeLeftNeutral  =
1309       new TH2F("hHbpXEUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",
1310                nptbins,ptmin,ptmax,200,0.,10.);
1311       fhPtHbpXEUeLeftNeutral->SetYTitle("ln(1/x_{E})");
1312       fhPtHbpXEUeLeftNeutral->SetXTitle("p_{T trigger}");
1313       outputContainer->Add(fhPtHbpXEUeLeftNeutral) ;
1314
1315       fhPtHbpXEUeRightNeutral  =
1316       new TH2F("hHbpXEUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",
1317                nptbins,ptmin,ptmax,200,0.,10.);
1318       fhPtHbpXEUeRightNeutral->SetYTitle("ln(1/x_{E})");
1319       fhPtHbpXEUeRightNeutral->SetXTitle("p_{T trigger}");
1320       outputContainer->Add(fhPtHbpXEUeRightNeutral) ;
1321
1322       fhZTUeLeftNeutral  =
1323       new TH2F("hZTUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
1324                nptbins,ptmin,ptmax,140,0.,2.);
1325       fhZTUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1326       fhZTUeLeftNeutral->SetXTitle("p_{T trigger}");
1327       outputContainer->Add(fhZTUeLeftNeutral) ;
1328
1329       fhZTUeRightNeutral  =
1330       new TH2F("hZTUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
1331                nptbins,ptmin,ptmax,200,0.,2.);
1332       fhZTUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1333       fhZTUeRightNeutral->SetXTitle("p_{T trigger}");
1334       outputContainer->Add(fhZTUeRightNeutral) ;
1335
1336       fhPtHbpZTUeLeftNeutral  =
1337       new TH2F("hHbpZTUeNeutralLeft","#xi = ln(1/z_{T}) with neutral UE left side of trigger",
1338                nptbins,ptmin,ptmax,200,0.,10.);
1339       fhPtHbpZTUeLeftNeutral->SetYTitle("ln(1/z_{T})");
1340       fhPtHbpZTUeLeftNeutral->SetXTitle("p_{T trigger}");
1341       outputContainer->Add(fhPtHbpZTUeLeftNeutral) ;
1342
1343       fhPtHbpZTUeRightNeutral  =
1344       new TH2F("hHbpZTUeNeutralRight","#xi = ln(1/z_{T}) with neutral UE right side of trigger",
1345                nptbins,ptmin,ptmax,200,0.,10.);
1346       fhPtHbpZTUeRightNeutral->SetYTitle("ln(1/z_{T})");
1347       fhPtHbpZTUeRightNeutral->SetXTitle("p_{T trigger}");
1348       outputContainer->Add(fhPtHbpZTUeRightNeutral) ;
1349
1350     }
1351
1352   }//Correlation with neutral hadrons
1353
1354   //if data is MC, fill more histograms
1355   if(IsDataMC())
1356   {
1357     fh2phiLeadingParticle=new TH2F("h2phiLeadingParticle","#phi resolustion for trigger particles",nptbins,ptmin,ptmax,100,-1,1);
1358     fh2phiLeadingParticle->GetXaxis()->SetTitle("p_{T gen Leading} (GeV/c)");
1359     fh2phiLeadingParticle->GetYaxis()->SetTitle("(#phi_{rec}-#phi_{gen})/#phi_{gen}");
1360
1361     fhMCEtaCharged  = new TH2F
1362     ("hMCEtaCharged","MC #eta_{h^{#pm}}  vs p_{T #pm}",
1363      nptbins,ptmin,ptmax,netabins,etamin,etamax);
1364     fhMCEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
1365     fhMCEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
1366
1367     fhMCPhiCharged  = new TH2F
1368     ("hMCPhiCharged","#MC phi_{h^{#pm}}  vs p_{T #pm}",
1369      200,ptmin,ptmax,nphibins,phimin,phimax);
1370     fhMCPhiCharged->SetYTitle("MC #phi_{h^{#pm}} (rad)");
1371     fhMCPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
1372
1373     fhMCDeltaPhiDeltaEtaCharged  = new TH2F
1374     ("hMCDeltaPhiDeltaEtaCharged","#MC phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
1375      140,-2.,5.,200,-2,2);
1376     fhMCDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
1377     fhMCDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1378
1379     fhMCDeltaEtaCharged  = new TH2F
1380     ("hMCDeltaEtaCharged","MC #eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger} and p_{T assoc}",
1381      nptbins,ptmin,ptmax,200,-2,2);
1382     fhMCDeltaEtaCharged->SetYTitle("#Delta #eta");
1383     fhMCDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
1384
1385     fhMCDeltaPhiCharged  = new TH2F
1386     ("hMCDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
1387      nptbins,ptmin,ptmax,140,-2.,5.);
1388     fhMCDeltaPhiCharged->SetYTitle("#Delta #phi");
1389     fhMCDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
1390
1391     fhMCDeltaPhiChargedPt  = new TH2F
1392     ("hMCDeltaPhiChargedPt","MC #phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
1393      nptbins,ptmin,ptmax,140,-2.,5.);
1394     fhMCDeltaPhiChargedPt->SetYTitle("#Delta #phi");
1395     fhMCDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
1396
1397     fhMCPtXECharged  =
1398     new TH2F("hMCPtXECharged","x_{E}",
1399              nptbins,ptmin,ptmax,200,0.,2.);
1400     fhMCPtXECharged->SetYTitle("x_{E}");
1401     fhMCPtXECharged->SetXTitle("p_{T trigger}");
1402
1403     fhMCPtHbpXECharged  =
1404     new TH2F("hMCHbpXECharged","MC #xi = ln(1/x_{E}) with charged hadrons",
1405              nptbins,ptmin,ptmax,200,0.,10.);
1406     fhMCPtHbpXECharged->SetYTitle("ln(1/x_{E})");
1407     fhMCPtHbpXECharged->SetXTitle("p_{T trigger}");
1408
1409     fhMCPtZTCharged  =
1410     new TH2F("hMCPtZTCharged","z_{T}",
1411              nptbins,ptmin,ptmax,200,0.,2.);
1412     fhMCPtZTCharged->SetYTitle("z_{T}");
1413     fhMCPtZTCharged->SetXTitle("p_{T trigger}");
1414
1415     fhMCPtHbpZTCharged  =
1416     new TH2F("hMCHbpZTCharged","MC #xi = ln(1/z_{T}) with charged hadrons",
1417              nptbins,ptmin,ptmax,200,0.,10.);
1418     fhMCPtHbpZTCharged->SetYTitle("ln(1/z_{T})");
1419     fhMCPtHbpZTCharged->SetXTitle("p_{T trigger}");
1420
1421     fhMCPtTrigPout  =
1422     new TH2F("hMCPtTrigPout","AOD MC Pout with triggers",
1423              nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
1424     fhMCPtTrigPout->SetYTitle("p_{out} (GeV/c)");
1425     fhMCPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
1426
1427     fhMCPtAssocDeltaPhi  =
1428     new TH2F("hMCPtAssocDeltaPhi","AOD MC delta phi with associated charged hadrons",
1429              nptbins,ptmin,ptmax,140,-2.,5.);
1430     fhMCPtAssocDeltaPhi->SetYTitle("#Delta #phi");
1431     fhMCPtAssocDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
1432
1433     outputContainer->Add(fh2phiLeadingParticle);
1434     outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged);
1435     outputContainer->Add(fhMCPhiCharged) ;
1436     outputContainer->Add(fhMCEtaCharged) ;
1437     outputContainer->Add(fhMCDeltaEtaCharged) ;
1438     outputContainer->Add(fhMCDeltaPhiCharged) ;
1439
1440     outputContainer->Add(fhMCDeltaPhiChargedPt) ;
1441     outputContainer->Add(fhMCPtXECharged) ;
1442     outputContainer->Add(fhMCPtZTCharged) ;
1443     outputContainer->Add(fhMCPtHbpXECharged) ;
1444     outputContainer->Add(fhMCPtHbpZTCharged) ;
1445     outputContainer->Add(fhMCPtTrigPout) ;
1446     outputContainer->Add(fhMCPtAssocDeltaPhi) ;
1447   } //for MC histogram
1448
1449   if(fDoOwnMix)
1450   {
1451     //create event containers
1452     fListMixEvents= new TList*[GetNCentrBin()*GetNZvertBin()] ;
1453
1454     for(Int_t ic=0; ic<GetNCentrBin(); ic++){
1455       for(Int_t iz=0; iz<GetNZvertBin(); iz++){
1456         fListMixEvents[ic*GetNZvertBin()+iz] = new TList() ;
1457         fListMixEvents[ic*GetNZvertBin()+iz]->SetOwner(kFALSE);
1458       }
1459     }
1460
1461     fhNtracksAll=new TH1F("hNtracksAll","Number of tracks w/o event trigger",2000,0,2000);
1462     outputContainer->Add(fhNtracksAll);
1463
1464     fhNtracksTrigger=new TH1F("hNtracksTriggerEvent","Number of tracks w/ event trigger",2000,0,2000);
1465     outputContainer->Add(fhNtracksTrigger);
1466
1467     fhNtracksINT=new TH1F("hNtracksMBEvent","Number of tracks w/ event trigger kAnyINT",2000,0,2000);
1468     outputContainer->Add(fhNtracksINT);
1469
1470     fhMixDeltaPhiCharged  = new TH2F
1471     ("hMixDeltaPhiCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
1472      nptbins,ptmin,ptmax,140,-2.,5.);
1473     fhMixDeltaPhiCharged->SetYTitle("#Delta #phi");
1474     fhMixDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
1475     outputContainer->Add(fhMixDeltaPhiCharged);
1476
1477     fhMixDeltaPhiDeltaEtaCharged  = new TH2F
1478     ("hMixDeltaPhiDeltaEtaCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
1479      140,-2.,5.,200,-2,2);
1480     fhMixDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
1481     fhMixDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1482     outputContainer->Add(fhMixDeltaPhiDeltaEtaCharged);
1483
1484     fhMixDeltaPhiChargedAssocPtBin         = new TH2F*[fNAssocPtBins] ;
1485     fhMixDeltaPhiDeltaEtaChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
1486
1487     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
1488     {
1489       fhMixDeltaPhiChargedAssocPtBin[i] = new TH2F(Form("hMixDeltaPhiChargedAssocPtBin%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1490                                          Form("Mixed event #Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1491                                          nptbins, ptmin, ptmax,140,-2.,5.);
1492       fhMixDeltaPhiChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
1493       fhMixDeltaPhiChargedAssocPtBin[i]->SetYTitle("#Delta #phi");
1494
1495       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i] = new TH2F(Form("hMixDeltaPhiDeltaEtaChargedAssocPtBin%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1496                                                    Form("Mixed event #Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1497                                                    140,-2.,5.,200,-2,2);
1498       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]->SetXTitle("#Delta #phi");
1499       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]->SetYTitle("#Delta #eta");
1500
1501       outputContainer->Add(fhMixDeltaPhiChargedAssocPtBin[i]);
1502       outputContainer->Add(fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]);
1503
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         = 7  ;
1583   fAssocPtBinLimit[0]   = 1.5  ;
1584   fAssocPtBinLimit[1]   = 3.  ;
1585   fAssocPtBinLimit[2]   = 5.  ;
1586   fAssocPtBinLimit[3]   = 7.  ;
1587   fAssocPtBinLimit[4]   = 9. ;
1588   fAssocPtBinLimit[5]   = 12. ;
1589   fAssocPtBinLimit[6]   = 15. ;
1590   fAssocPtBinLimit[7]   = 20. ;
1591   fAssocPtBinLimit[8]   = 100.;
1592   fAssocPtBinLimit[9]   = 200.;
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   if(GetDebug()>1)
2221     printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation() - Make trigger particle - charged hadron correlation in AOD MC level\n");
2222
2223   AliStack         * stack        = 0x0 ;
2224   TParticle        * primary      = 0x0 ;
2225   TClonesArray     * mcparticles0 = 0x0 ;
2226   TClonesArray     * mcparticles  = 0x0 ;
2227   AliAODMCParticle * aodprimary   = 0x0 ;
2228
2229   Double_t eprim   = 0 ;
2230   Double_t ptprim  = 0 ;
2231   Double_t phiprim = 0 ;
2232   Double_t etaprim = 0 ;
2233   Int_t    nTracks = 0 ;
2234   Int_t iParticle  = 0 ;
2235   Double_t charge  = 0.;
2236
2237   if(GetReader()->ReadStack())
2238   {
2239     nTracks = GetMCStack()->GetNtrack() ;
2240   }
2241   else
2242   {
2243     nTracks = GetReader()->GetAODMCParticles()->GetEntriesFast() ;
2244   }
2245   //Int_t trackIndex[nTracks];
2246
2247   Int_t label= aodParticle->GetLabel();
2248   if(label<0){
2249     printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** bad label ***:  label %d \n", label);
2250     return;
2251   }
2252
2253   if(GetReader()->ReadStack())
2254   {
2255     stack =  GetMCStack() ;
2256     if(!stack) {
2257       printf(" AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation- Stack not available, is the MC handler called? STOP\n");
2258       abort();
2259     }
2260
2261     nTracks=stack->GetNprimary();
2262     if(label >=  stack->GetNtrack())
2263     {
2264       if(GetDebug() > 2)  printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
2265       return ;
2266     }
2267
2268     primary = stack->Particle(label);
2269     if(!primary)
2270     {
2271       printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no primary ***:  label %d \n", label);
2272       return;
2273     }
2274
2275     if(primary)
2276     {
2277       eprim    = primary->Energy();
2278       ptprim   = primary->Pt();
2279       phiprim  = primary->Phi();
2280       etaprim  = primary->Eta();
2281
2282       if(ptprim < 0.01 || eprim < 0.01) return ;
2283
2284       for (iParticle = 0 ; iParticle <  nTracks ; iParticle++)
2285       {
2286         TParticle * particle = stack->Particle(iParticle);
2287         TLorentzVector momentum;
2288
2289         //keep only final state particles
2290         if(particle->GetStatusCode()!=1) continue ;
2291
2292         Int_t pdg = particle->GetPdgCode();
2293
2294         charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
2295
2296         particle->Momentum(momentum);
2297
2298         //---------- Charged particles ----------------------
2299         if(charge != 0)
2300         {
2301           //Particles in CTS acceptance
2302           Bool_t inCTS =  GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
2303
2304           if(TMath::Abs(pdg) == 11 && stack->Particle(particle->GetFirstMother())->GetPdgCode()==22) continue ;
2305
2306           if(inCTS)
2307           {
2308             if( label!=iParticle) // avoid trigger particle
2309             {
2310               if(!FillChargedMCCorrelationHistograms(particle->Pt(),particle->Phi(),particle->Eta(),ptprim,phiprim,etaprim)) return;
2311             }
2312           }// in CTS acceptance
2313         }// charged
2314       } //track loop
2315     } //when the leading particles could trace back to MC
2316   } //ESD MC
2317   else if(GetReader()->ReadAODMCParticles())
2318   {
2319     //Get the list of MC particles
2320     mcparticles0 = GetReader()->GetAODMCParticles(0);
2321     if(!mcparticles0) return;
2322     if(label >=mcparticles0->GetEntriesFast())
2323     {
2324       if(GetDebug() > 2)
2325         printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** large label ***:  label %d, n tracks %d \n", label,mcparticles0->GetEntriesFast());
2326       return;
2327     }
2328     //Get the particle
2329     aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
2330     if(!aodprimary)
2331     {
2332       printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no AOD primary ***:  label %d \n", label);
2333       return;
2334     }
2335
2336     if(aodprimary)
2337     {
2338       ptprim  = aodprimary->Pt();
2339       phiprim = aodprimary->Phi();
2340       etaprim = aodprimary->Eta();
2341
2342       if(ptprim < 0.01 || eprim < 0.01) return ;
2343
2344       mcparticles= GetReader()->GetAODMCParticles();
2345       for (Int_t i = 0; i < nTracks; i++)
2346       {
2347         AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(i);
2348         if (!part->IsPhysicalPrimary()) continue;
2349         Int_t pdg = part->GetPdgCode();
2350         charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
2351         TLorentzVector momentum(part->Px(),part->Py(),part->Pz(),part->E());
2352         if(charge != 0)
2353         {
2354           if(part->Pt()> GetReader()->GetCTSPtMin())
2355           {
2356             //Particles in CTS acceptance
2357             Bool_t inCTS =  GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
2358             Int_t indexmother=part->GetMother();
2359             if(indexmother>-1)
2360             {
2361               Int_t mPdg = ((AliAODMCParticle*) mcparticles->At(indexmother)) ->GetPdgCode();
2362               if(TMath::Abs(pdg) == 11 && mPdg == 22) continue;
2363             }
2364
2365             if(inCTS)
2366             {
2367               if( label!=iParticle) // avoid trigger particle
2368               {
2369                 if(!FillChargedMCCorrelationHistograms(part->Pt(),part->Phi(),part->Eta(),ptprim,phiprim,etaprim)) return;
2370               }
2371             } // in acceptance
2372           } // min pt cut
2373         } //only charged particles
2374       }  //MC particle loop
2375     } //when the leading particles could trace back to MC
2376   }// AOD MC
2377 }
2378
2379 //_____________________________________________________________________
2380 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
2381 {
2382
2383   //Print some relevant parameters set for the analysis
2384   if(! opt)
2385     return;
2386
2387   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2388   AliAnaCaloTrackCorrBaseClass::Print(" ");
2389   printf("Pt trigger           >  %3.2f\n", fMinTriggerPt) ;
2390   printf("Pt associated hadron <  %3.2f\n", fMaxAssocPt) ;
2391   printf("Pt associated hadron >  %3.2f\n", fMinAssocPt) ;
2392   printf("Phi trigger particle-Hadron   <  %3.2f\n", fDeltaPhiMaxCut) ;
2393   printf("Phi trigger particle-Hadron   >  %3.2f\n", fDeltaPhiMinCut) ;
2394   printf("Phi trigger particle-UeHadron <  %3.2f\n", fUeDeltaPhiMaxCut) ;
2395   printf("Phi trigger particle-UeHadron >  %3.2f\n", fUeDeltaPhiMinCut) ;
2396   printf("Isolated Trigger?  %d\n"     , fSelectIsolated) ;
2397   printf("Several UE?  %d\n"           , fMakeSeveralUE) ;
2398   printf("Name of AOD Pi0 Branch %s \n", fPi0AODBranchName.Data());
2399   printf("Do Decay-hadron correlation ?  %d\n", fPi0Trigger) ;
2400   printf("Select absolute leading for cluster triggers ?  %d or Near Side %d\n",
2401          fMakeAbsoluteLeading, fMakeNearSideLeading) ;
2402   printf("Trigger pt bins  %d\n", fNAssocPtBins) ;
2403   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
2404     printf("\t bin %d = [%2.1f,%2.1f]\n", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
2405   }
2406
2407 }
2408
2409 //____________________________________________________________
2410 void AliAnaParticleHadronCorrelation::SetNAssocPtBins(Int_t n)
2411 {
2412   // Set number of bins
2413
2414   fNAssocPtBins  = n ;
2415
2416
2417   if(n < 10 && n > 0)
2418   {
2419     fNAssocPtBins  = n ;
2420   }
2421   else
2422   {
2423     printf("n = larger than 9 or too small, set to 9 \n");
2424     fNAssocPtBins = 9;
2425   }
2426 }
2427
2428 //______________________________________________________________________________
2429 void AliAnaParticleHadronCorrelation::SetAssocPtBinLimit(Int_t ibin, Float_t pt)
2430 {
2431   // Set the list of limits for the trigger pt bins
2432
2433   if(ibin <= fNAssocPtBins || ibin >= 0)
2434   {
2435     fAssocPtBinLimit[ibin] = pt  ;
2436   }
2437   else {
2438     printf("AliAnaParticleHadronCorrelation::SetAssocPtBinLimit() - bin  number too large %d > %d or small, nothing done\n", ibin, fNAssocPtBins) ;
2439
2440   }
2441 }
2442