PWGGA/CaloTrackCorrelations/AliAnaParticleHadronCorrelation.cxx

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