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   if(fUseMixStoredInReader && GetReader()->GetLastTracksMixedEvent() == GetEventNumber())
 578   {
 579     //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
 580     return ; // pool filled previously for another trigger
 581   }
 582
 583   Int_t nTracks = GetCTSTracks()->GetEntriesFast();
 584
 585   fhNtracksAll->Fill(nTracks);
 586
 587   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
 588   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
 589
 590   if(!inputHandler) return ;
 591
 592   if( inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask()    )
 593   {
 594     fhNtracksTrigger->Fill(nTracks);
 595   }
 596
 597   if( inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask() )
 598   {
 599
 600     fhNtracksMB->Fill(nTracks);
 601
 602     Int_t eventBin = GetEventMixBin();
 603
 604     //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
 605     if(eventBin < 0) return;
 606
 607     TObjArray * mixEventTracks = new TObjArray;
 608
 609     if(fUseMixStoredInReader)
 610     {
 611       fListMixEvents[eventBin] = GetReader()->GetListWithMixedEventsForTracks(eventBin);
 612     }
 613
 614     if(!fListMixEvents[eventBin]) fListMixEvents[eventBin] = new TList();
 615
 616     //printf("%s ***** Pool Event bin : %d - nTracks %d\n",GetInputAODName().Data(),eventBin, GetCTSTracks()->GetEntriesFast());
 617
 618     TList * pool = fListMixEvents[eventBin];
 619
 620     TVector3 p3;
 621     for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
 622     {
 623       AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
 624
 625       Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 626       p3.SetXYZ(mom[0],mom[1],mom[2]);
 627       Float_t pt   = p3.Pt();
 628
 629       //Select only hadrons in pt range
 630       if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
 631
 632       AliAODPWG4Particle * mixedTrack = new AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
 633       mixedTrack->SetDetector("CTS");
 634       mixedTrack->SetChargedBit(track->Charge()>0);
 635
 636       mixEventTracks->Add(mixedTrack);
 637     }
 638
 639     //Set the event number where the last event was added, to avoid double pool filling
 640     GetReader()->SetLastTracksMixedEvent(GetEventNumber());
 641
 642     pool->AddFirst(mixEventTracks);
 643     mixEventTracks = 0;
 644     //printf("Pool size %d, max %d\n",pool->GetSize(), GetNMaxEvMix());
 645
 646     if(pool->GetSize() > GetNMaxEvMix())
 647     {//Remove last event
 648       TClonesArray * tmp = static_cast<TClonesArray*>(pool->Last()) ;
 649       pool->RemoveLast() ;
 650       delete tmp ;
 651     }
 652
 653   } // MB event
 654
 655 }
 656
 657 //____________________________________________________________
 658 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
 659 {
 660   //Save parameters used for analysis
 661   TString parList ; //this will be list of parameters used for this analysis.
 662   const Int_t buffersize = 560;
 663   char onePar[buffersize] ;
 664
 665   snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;
 666   parList+=onePar ;
 667   snprintf(onePar,buffersize," Pt Trigger > %3.2f ",    fMinTriggerPt) ;
 668   parList+=onePar ;
 669   snprintf(onePar,buffersize," %3.2f < Pt associated < %3.2f ", fMinAssocPt,   fMaxAssocPt) ;
 670   parList+=onePar ;
 671   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-Hadron < %3.2f ",    fDeltaPhiMinCut,   fDeltaPhiMaxCut) ;
 672   parList+=onePar ;
 673   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-UeHadron <  %3.2f ", fUeDeltaPhiMinCut, fUeDeltaPhiMaxCut) ;
 674   parList+=onePar ;
 675   snprintf(onePar,buffersize,"Isolated Trigger?  %d\n",    fSelectIsolated) ;
 676   parList+=onePar ;
 677   snprintf(onePar,buffersize,"Several UE?  %d\n",          fMakeSeveralUE) ;
 678   parList+=onePar ;
 679   snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ", fPi0AODBranchName.Data());
 680   parList+=onePar ;
 681   snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  pi0 %d, decay %d", fPi0Trigger, fDecayTrigger) ;
 682   parList+=onePar ;
 683   snprintf(onePar,buffersize,"Select absolute leading for cluster triggers ? %d or Near Side Leading %d \n",
 684            fMakeAbsoluteLeading, fMakeNearSideLeading) ;
 685   parList+=onePar ;
 686   snprintf(onePar,buffersize,"Associated particle pt bins  %d: ", fNAssocPtBins) ;
 687   parList+=onePar ;
 688   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
 689     snprintf(onePar,buffersize,"bin %d = [%2.1f,%2.1f];", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
 690   }
 691   parList+=onePar ;
 692
 693   //Get parameters set in base class.
 694   parList += GetBaseParametersList() ;
 695
 696   //Get parameters set in FiducialCut class (not available yet)
 697   //parlist += GetFidCut()->GetFidCutParametersList()
 698
 699   return new TObjString(parList) ;
 700
 701 }
 702
 703 //________________________________________________________________
 704 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
 705 {
 706
 707   // Create histograms to be saved in output file and
 708   // store them in fOutputContainer
 709
 710   TList * outputContainer = new TList() ;
 711   outputContainer->SetName("CorrelationHistos") ;
 712
 713   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();
 714   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();
 715   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();
 716
 717   fhPtLeading  = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax);
 718   fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
 719
 720   fhPhiLeading  = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 721   fhPhiLeading->SetYTitle("#phi (rad)");
 722
 723   fhEtaLeading  = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 724   fhEtaLeading->SetYTitle("#eta ");
 725
 726   outputContainer->Add(fhPtLeading);
 727   outputContainer->Add(fhPhiLeading);
 728   outputContainer->Add(fhEtaLeading);
 729
 730   fhPtLeadingCentrality   = new TH2F("hPtLeadingCentrality","Leading particle p_{T} vs centrality",nptbins,ptmin,ptmax,100,0.,100) ;
 731   fhPtLeadingCentrality->SetXTitle("p_{T}^{trig} (GeV/c)");
 732   fhPtLeadingCentrality->SetYTitle("Centrality (%)");
 733   outputContainer->Add(fhPtLeadingCentrality) ;
 734
 735   fhPtLeadingEventPlane  = new TH2F("hPtLeadingEventPlane","Leading particle p_{T} vs event plane angle",nptbins,ptmin,ptmax, 100,0.,TMath::Pi()) ;
 736   fhPtLeadingEventPlane->SetXTitle("p_{T}^{trig} (GeV/c)");
 737   fhPtLeadingEventPlane->SetXTitle("EP angle (rad)");
 738   outputContainer->Add(fhPtLeadingEventPlane) ;
 739
 740   fhLeadingEventPlaneCentrality  = new TH2F("hLeadingEventPlane","Leading particle centrality vs event plane angle",100,0.,100,100,0.,TMath::Pi()) ;
 741   fhLeadingEventPlaneCentrality->SetXTitle("Centrality (%)");
 742   fhLeadingEventPlaneCentrality->SetYTitle("EP angle (rad)");
 743   outputContainer->Add(fhLeadingEventPlaneCentrality) ;
 744
 745   //Correlation with charged hadrons
 746   if(GetReader()->IsCTSSwitchedOn())
 747   {
 748     fhDeltaPhiDeltaEtaCharged  = new TH2F
 749     ("hDeltaPhiDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}}",
 750     ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
 751     fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
 752     fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
 753
 754     fhPhiCharged  = new TH2F
 755     ("hPhiCharged","#phi_{h^{#pm}}  vs p_{T #pm}",
 756      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 757     fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
 758     fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
 759
 760     fhEtaCharged  = new TH2F
 761     ("hEtaCharged","#eta_{h^{#pm}}  vs p_{T #pm}",
 762      nptbins,ptmin,ptmax,netabins,etamin,etamax);
 763     fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
 764     fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
 765
 766     fhDeltaPhiCharged  = new TH2F
 767     ("hDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
 768      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
 769     fhDeltaPhiCharged->SetYTitle("#Delta #phi");
 770     fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
 771
 772     fhDeltaPhiChargedPt  = new TH2F
 773     ("hDeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
 774      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
 775     fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");
 776     fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 777
 778     fhDeltaPhiUeChargedPt  = new TH2F
 779     ("hDeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",
 780      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
 781     fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");
 782     fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
 783
 784     fhUePart  =  new TH1F("hUePart","UE particles distribution vs pt trig",
 785              nptbins,ptmin,ptmax);
 786     fhUePart->SetYTitle("dNch");
 787     fhUePart->SetXTitle("p_{T trigger}");
 788
 789
 790     fhDeltaEtaCharged  = new TH2F
 791     ("hDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",
 792      nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
 793     fhDeltaEtaCharged->SetYTitle("#Delta #eta");
 794     fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
 795
 796     fhXECharged  =
 797     new TH2F("hXECharged","x_{E} for charged tracks",
 798              nptbins,ptmin,ptmax,200,0.,2.);
 799     fhXECharged->SetYTitle("x_{E}");
 800     fhXECharged->SetXTitle("p_{T trigger}");
 801
 802     fhXEUeCharged  =
 803     new TH2F("hXEUeCharged","x_{E} for Underlying Event",
 804              nptbins,ptmin,ptmax,200,0.,2.);
 805     fhXEUeCharged->SetYTitle("x_{E}");
 806     fhXEUeCharged->SetXTitle("p_{T trigger}");
 807
 808     fhXEPosCharged  =
 809     new TH2F("hXEPositiveCharged","x_{E} for positive charged tracks",
 810              nptbins,ptmin,ptmax,200,0.,2.);
 811     fhXEPosCharged->SetYTitle("x_{E}");
 812     fhXEPosCharged->SetXTitle("p_{T trigger}");
 813
 814     fhXENegCharged  =
 815     new TH2F("hXENegativeCharged","x_{E} for negative charged tracks",
 816              nptbins,ptmin,ptmax,200,0.,2.);
 817     fhXENegCharged->SetYTitle("x_{E}");
 818     fhXENegCharged->SetXTitle("p_{T trigger}");
 819
 820     fhPtHbpXECharged  =
 821     new TH2F("hHbpXECharged","#xi = ln(1/x_{E}) with charged hadrons",
 822              nptbins,ptmin,ptmax,200,0.,10.);
 823     fhPtHbpXECharged->SetYTitle("ln(1/x_{E})");
 824     fhPtHbpXECharged->SetXTitle("p_{T trigger}");
 825
 826     fhPtHbpXEUeCharged  =
 827     new TH2F("hHbpXEUeCharged","#xi = ln(1/x_{E}) with charged hadrons,Underlying Event",
 828              nptbins,ptmin,ptmax,200,0.,10.);
 829     fhPtHbpXEUeCharged->SetYTitle("ln(1/x_{E})");
 830     fhPtHbpXEUeCharged->SetXTitle("p_{T trigger}");
 831
 832     fhZTCharged  =
 833     new TH2F("hZTCharged","z_{T} for charged tracks",
 834              nptbins,ptmin,ptmax,200,0.,2.);
 835     fhZTCharged->SetYTitle("z_{T}");
 836     fhZTCharged->SetXTitle("p_{T trigger}");
 837
 838     fhZTUeCharged  =
 839     new TH2F("hZTUeCharged","z_{T} for Underlying Event",
 840              nptbins,ptmin,ptmax,200,0.,2.);
 841     fhZTUeCharged->SetYTitle("z_{T}");
 842     fhZTUeCharged->SetXTitle("p_{T trigger}");
 843
 844     fhZTPosCharged  =
 845     new TH2F("hZTPositiveCharged","z_{T} for positive charged tracks",
 846              nptbins,ptmin,ptmax,200,0.,2.);
 847     fhZTPosCharged->SetYTitle("z_{T}");
 848     fhZTPosCharged->SetXTitle("p_{T trigger}");
 849
 850     fhZTNegCharged  =
 851     new TH2F("hZTNegativeCharged","z_{T} for negative charged tracks",
 852              nptbins,ptmin,ptmax,200,0.,2.);
 853     fhZTNegCharged->SetYTitle("z_{T}");
 854     fhZTNegCharged->SetXTitle("p_{T trigger}");
 855
 856     fhPtHbpZTCharged  =
 857     new TH2F("hHbpZTCharged","#xi = ln(1/z_{T}) with charged hadrons",
 858              nptbins,ptmin,ptmax,200,0.,10.);
 859     fhPtHbpZTCharged->SetYTitle("ln(1/z_{T})");
 860     fhPtHbpZTCharged->SetXTitle("p_{T trigger}");
 861
 862     fhPtHbpZTUeCharged  =
 863     new TH2F("hHbpZTUeCharged","#xi = ln(1/z_{T}) with charged hadrons,Underlying Event",
 864              nptbins,ptmin,ptmax,200,0.,10.);
 865     fhPtHbpZTUeCharged->SetYTitle("ln(1/x_{E})");
 866     fhPtHbpZTUeCharged->SetXTitle("p_{T trigger}");
 867
 868     fhPtTrigPout  =
 869     new TH2F("hPtTrigPout","Pout with triggers",
 870              nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
 871     fhPtTrigPout->SetYTitle("p_{out} (GeV/c)");
 872     fhPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
 873
 874     fhPtTrigCharged  =
 875     new TH2F("hPtTrigCharged","trgger and charged tracks pt distribution",
 876              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 877     fhPtTrigCharged->SetYTitle("p_{T h^{#pm}} (GeV/c)");
 878     fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)");
 879
 880     outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
 881     outputContainer->Add(fhPhiCharged) ;
 882     outputContainer->Add(fhEtaCharged) ;
 883     outputContainer->Add(fhDeltaPhiCharged) ;
 884     outputContainer->Add(fhDeltaEtaCharged) ;
 885     outputContainer->Add(fhDeltaPhiChargedPt) ;
 886     outputContainer->Add(fhDeltaPhiUeChargedPt) ;
 887     outputContainer->Add(fhUePart);
 888
 889     outputContainer->Add(fhXECharged) ;
 890     outputContainer->Add(fhXEPosCharged) ;
 891     outputContainer->Add(fhXENegCharged) ;
 892     outputContainer->Add(fhXEUeCharged) ;
 893     outputContainer->Add(fhPtHbpXECharged) ;
 894     outputContainer->Add(fhPtHbpXEUeCharged) ;
 895
 896     outputContainer->Add(fhZTCharged) ;
 897     outputContainer->Add(fhZTPosCharged) ;
 898     outputContainer->Add(fhZTNegCharged) ;
 899     outputContainer->Add(fhZTUeCharged) ;
 900     outputContainer->Add(fhPtHbpZTCharged) ;
 901     outputContainer->Add(fhPtHbpZTUeCharged) ;
 902
 903     outputContainer->Add(fhPtTrigPout) ;
 904     outputContainer->Add(fhPtTrigCharged) ;
 905
 906     if(DoEventSelect())
 907     {
 908       Int_t nMultiBins = GetMultiBin();
 909       fhTrigDeltaPhiCharged = new TH2F*[nMultiBins] ;
 910       fhTrigDeltaEtaCharged = new TH2F*[nMultiBins] ;
 911       fhTrigXECorr          = new TH2F*[nMultiBins] ;
 912       fhTrigXEUeCorr        = new TH2F*[nMultiBins] ;
 913       fhTrigZTCorr          = new TH2F*[nMultiBins] ;
 914       fhTrigZTUeCorr        = new TH2F*[nMultiBins] ;
 915
 916       for(Int_t im=0; im<nMultiBins; im++)
 917       {
 918         fhTrigDeltaPhiCharged[im]  = new TH2F
 919         (Form("hTrigDeltaPhiCharged_%d",im),Form("hTrigDeltaPhiCharged_%d",im), nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
 920         fhTrigDeltaPhiCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
 921         fhTrigDeltaPhiCharged[im]->SetYTitle("#Delta #phi");
 922
 923         fhTrigDeltaEtaCharged[im]  = new TH2F
 924         (Form("hTrigDeltaEtaCharged_%d",im),Form("hTrigDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
 925         fhTrigDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");
 926         fhTrigDeltaEtaCharged[im]->SetYTitle("#Delta #eta");
 927
 928         fhTrigXECorr[im]  = new TH2F
 929         (Form("hTrigXEPtCorr_%d",im),Form("hTrigXEPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 930         fhTrigXECorr[im]->SetYTitle("x_{E trigger h^{#pm}}");
 931         fhTrigXECorr[im]->SetXTitle("p_{T trigger}");
 932
 933         fhTrigXEUeCorr[im]  = new TH2F
 934         (Form("hTrigXEPtUeCorr_%d",im),Form("hTrigXEPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 935         fhTrigXEUeCorr[im]->SetYTitle("x_{E trigger h^{#pm}}");
 936         fhTrigXEUeCorr[im]->SetXTitle("p_{T trigger}");
 937
 938         fhTrigZTCorr[im]  = new TH2F
 939         (Form("hTrigZTPtCorr_%d",im),Form("hTrigZTPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 940         fhTrigZTCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
 941         fhTrigZTCorr[im]->SetXTitle("p_{T trigger}");
 942
 943         fhTrigZTUeCorr[im]  = new TH2F
 944         (Form("hTrigZTPtUeCorr_%d",im),Form("hTrigZTPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.);
 945         fhTrigZTUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");
 946         fhTrigZTUeCorr[im]->SetXTitle("p_{T trigger}");
 947
 948         outputContainer->Add(fhTrigDeltaPhiCharged[im]) ;
 949         outputContainer->Add(fhTrigDeltaEtaCharged[im]) ;
 950         outputContainer->Add(fhTrigXECorr[im]);
 951         outputContainer->Add(fhTrigXEUeCorr[im]);
 952         outputContainer->Add(fhTrigZTCorr[im]);
 953         outputContainer->Add(fhTrigZTUeCorr[im]);
 954       }
 955     }
 956
 957     if(fFillBradHisto)
 958     {
 959       fhAssocPtBkg        = new TH2F("hAssocPtBkg", " Trigger p_{T} vs associated hadron p_{T} from background",
 960                                      nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
 961       fhAssocPtBkg->SetXTitle("p_{T trigger}");
 962       fhAssocPtBkg->SetYTitle("p_{T associated}");
 963       outputContainer->Add(fhAssocPtBkg) ;
 964
 965       fhDeltaPhiBrad = new TH2F("hDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs p_{T trigger} ",
 966                                 nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
 967       fhDeltaPhiBrad->SetXTitle("p_{T trigger}");
 968       fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
 969       outputContainer->Add(fhDeltaPhiBrad) ;
 970     }
 971
 972     fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins] ;
 973     fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins] ;
 974     fhXEAssocPtBin             = new TH2F*[fNAssocPtBins] ;
 975     fhZTAssocPtBin             = new TH2F*[fNAssocPtBins] ;
 976
 977     if(fFillBradHisto)
 978       fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins] ;
 979
 980     if(fPi0Trigger || fDecayTrigger)
 981     {
 982       fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins] ;
 983       fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins] ;
 984       fhXEAssocPtBin             = new TH2F*[fNAssocPtBins] ;
 985       fhZTAssocPtBin             = new TH2F*[fNAssocPtBins] ;
 986       fhXEDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
 987       fhZTDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
 988       fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
 989     }
 990
 991     if(fHMPIDCorrelation)
 992     {
 993       fhDeltaPhiAssocPtBinHMPID   = new TH2F*[fNAssocPtBins] ;
 994       fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins] ;
 995     }
 996
 997     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
 998     {
 999       fhDeltaPhiAssocPtBin[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1000                                          Form("#Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1001                                          nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1002       fhDeltaPhiAssocPtBin[i]->SetXTitle("p_{T trigger}");
1003       fhDeltaPhiAssocPtBin[i]->SetYTitle("#Delta #phi");
1004
1005       fhDeltaPhiAssocPtBinDEta08[i] = new TH2F(Form("hDeltaPhiDeltaEta0.8PtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1006                                                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]),
1007                                                nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1008       fhDeltaPhiAssocPtBinDEta08[i]->SetXTitle("p_{T trigger}");
1009       fhDeltaPhiAssocPtBinDEta08[i]->SetYTitle("#Delta #phi");
1010
1011       fhXEAssocPtBin[i]       = new TH2F(Form("hXEAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1012                                          Form("x_{E} vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1013                                          nptbins, ptmin, ptmax,200, 0.0, 2.0);
1014       fhXEAssocPtBin[i]->SetXTitle("p_{T trigger}");
1015       fhXEAssocPtBin[i]->SetYTitle("x_{E}");
1016
1017       fhZTAssocPtBin[i]       = new TH2F(Form("hZTAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1018                                          Form("z_{T} vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1019                                          nptbins, ptmin, ptmax,200, 0.0, 2.0);
1020       fhZTAssocPtBin[i]->SetXTitle("p_{T trigger}");
1021       fhZTAssocPtBin[i]->SetYTitle("z_{T}");
1022
1023       outputContainer->Add(fhDeltaPhiAssocPtBin[i]) ;
1024       outputContainer->Add(fhDeltaPhiAssocPtBinDEta08[i]) ;
1025       outputContainer->Add(fhXEAssocPtBin[i]);
1026       outputContainer->Add(fhZTAssocPtBin[i]);
1027
1028       if(fPi0Trigger || fDecayTrigger)
1029       {
1030         fhDeltaPhiDecayChargedAssocPtBin[i] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1031                                            Form("#Delta #phi vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1032                                            nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1033         fhDeltaPhiDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
1034         fhDeltaPhiDecayChargedAssocPtBin[i]->SetYTitle("#Delta #phi");
1035
1036         fhXEDecayChargedAssocPtBin[i]       = new TH2F(Form("hXEDecayChargedAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1037                                            Form("x_{E} vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1038                                            nptbins, ptmin, ptmax,200, 0.0, 2.0);
1039         fhXEDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
1040         fhXEDecayChargedAssocPtBin[i]->SetYTitle("x_{E}");
1041
1042         fhZTDecayChargedAssocPtBin[i]       = new TH2F(Form("hZTDecayChargedAssocPtBin%1.f_%1.f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1043                                            Form("z_{T} vs p_{T trigger} tagged as decay for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1044                                            nptbins, ptmin, ptmax,200, 0.0, 2.0);
1045         fhZTDecayChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
1046         fhZTDecayChargedAssocPtBin[i]->SetYTitle("z_{T}");
1047
1048         outputContainer->Add(fhDeltaPhiDecayChargedAssocPtBin[i]) ;
1049         outputContainer->Add(fhXEDecayChargedAssocPtBin[i]);
1050         outputContainer->Add(fhZTDecayChargedAssocPtBin[i]);
1051
1052       }
1053
1054       if(fFillBradHisto)
1055       {
1056         fhDeltaPhiBradAssocPtBin[i] = new TH2F(Form("hDeltaPhiBradPtAssocPt%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1057                                                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]),
1058                                                nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
1059         fhDeltaPhiBradAssocPtBin[i]->SetXTitle("p_{T trigger}");
1060         fhDeltaPhiBradAssocPtBin[i]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
1061         outputContainer->Add(fhDeltaPhiBradAssocPtBin[i]) ;
1062       }
1063
1064       if(fHMPIDCorrelation)
1065       {
1066         fhDeltaPhiAssocPtBinHMPID[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1fHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1067                                                 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]),
1068                                                 nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1069         fhDeltaPhiAssocPtBinHMPID[i]->SetXTitle("p_{T trigger}");
1070         fhDeltaPhiAssocPtBinHMPID[i]->SetYTitle("#Delta #phi");
1071
1072         fhDeltaPhiAssocPtBinHMPIDAcc[i] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1fHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1073                                                    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]),
1074                                                    nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1075         fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetXTitle("p_{T trigger}");
1076         fhDeltaPhiAssocPtBinHMPIDAcc[i]->SetYTitle("#Delta #phi");
1077
1078         outputContainer->Add(fhDeltaPhiAssocPtBinHMPID[i]) ;
1079         outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[i]) ;
1080
1081       }
1082     }
1083
1084     if(fPi0Trigger || fDecayTrigger)
1085     {
1086       if(fPi0Trigger)
1087       {
1088         fhPtPi0DecayRatio  = new TH2F
1089         ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay",
1090          nptbins,ptmin,ptmax, 100,0.,2.);
1091         fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");
1092         fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");
1093         outputContainer->Add(fhPtPi0DecayRatio) ;
1094       }
1095
1096       fhDeltaPhiDecayCharged  = new TH2F
1097       ("hDeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",
1098        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1099       fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");
1100       fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");
1101
1102       fhXEDecayCharged  =
1103       new TH2F("hXEDecayCharged","x_{E}  Decay",
1104                nptbins,ptmin,ptmax,200,0.,2.);
1105       fhXEDecayCharged->SetYTitle("x_{E}");
1106       fhXEDecayCharged->SetXTitle("p_{T decay}");
1107
1108       fhZTDecayCharged  =
1109       new TH2F("hZTDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",
1110                nptbins,ptmin,ptmax,200,0.,2.);
1111       fhZTDecayCharged->SetYTitle("z_{decay h^{#pm}}");
1112       fhZTDecayCharged->SetXTitle("p_{T decay}");
1113
1114       outputContainer->Add(fhDeltaPhiDecayCharged) ;
1115       outputContainer->Add(fhXEDecayCharged) ;
1116       outputContainer->Add(fhZTDecayCharged) ;
1117     }
1118
1119     if(fMakeSeveralUE)
1120     {
1121       fhDeltaPhiUeLeftCharged  = new TH2F
1122       ("hDeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",
1123        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1124       fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");
1125       fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
1126       outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
1127
1128       fhDeltaPhiUeRightCharged  = new TH2F
1129       ("hDeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",
1130        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1131       fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");
1132       fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");
1133       outputContainer->Add(fhDeltaPhiUeRightCharged) ;
1134
1135       fhXEUeLeftCharged  =
1136       new TH2F("hXEUeChargedLeft","x_{E} with UE left side of trigger",
1137                nptbins,ptmin,ptmax,200,0.,2.);
1138       fhXEUeLeftCharged->SetYTitle("x_{E Ueh^{#pm}}");
1139       fhXEUeLeftCharged->SetXTitle("p_{T trigger}");
1140       outputContainer->Add(fhXEUeLeftCharged) ;
1141
1142       fhXEUeRightCharged  =
1143       new TH2F("hXEUeChargedRight","x_{E h^{#pm}} with UE right side of trigger",
1144                nptbins,ptmin,ptmax,200,0.,2.);
1145       fhXEUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
1146       fhXEUeRightCharged->SetXTitle("p_{T trigger}");
1147       outputContainer->Add(fhXEUeRightCharged) ;
1148
1149       fhPtHbpXEUeLeftCharged  =
1150       new TH2F("hHbpXEUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",
1151                nptbins,ptmin,ptmax,200,0.,10.);
1152       fhPtHbpXEUeLeftCharged->SetYTitle("ln(1/x_{E})");
1153       fhPtHbpXEUeLeftCharged->SetXTitle("p_{T trigger}");
1154       outputContainer->Add(fhPtHbpXEUeLeftCharged) ;
1155
1156       fhPtHbpXEUeRightCharged  =
1157       new TH2F("hHbpXEUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",
1158                nptbins,ptmin,ptmax,200,0.,10.);
1159       fhPtHbpXEUeRightCharged->SetYTitle("ln(1/x_{E})");
1160       fhPtHbpXEUeRightCharged->SetXTitle("p_{T trigger}");
1161       outputContainer->Add(fhPtHbpXEUeRightCharged) ;
1162
1163       fhZTUeLeftCharged  =
1164       new TH2F("hZTUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",
1165                nptbins,ptmin,ptmax,200,0.,2.);
1166       fhZTUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
1167       fhZTUeLeftCharged->SetXTitle("p_{T trigger}");
1168       outputContainer->Add(fhZTUeLeftCharged) ;
1169
1170       fhZTUeRightCharged  =
1171       new TH2F("hZTUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",
1172                nptbins,ptmin,ptmax,200,0.,2.);
1173       fhZTUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");
1174       fhZTUeRightCharged->SetXTitle("p_{T trigger}");
1175       outputContainer->Add(fhZTUeRightCharged) ;
1176
1177       fhPtHbpZTUeLeftCharged  =
1178       new TH2F("hHbpZTUeChargedLeft","#xi = ln(1/z_{T}) with charged UE left side of trigger",
1179                nptbins,ptmin,ptmax,200,0.,10.);
1180       fhPtHbpZTUeLeftCharged->SetYTitle("ln(1/z_{T})");
1181       fhPtHbpZTUeLeftCharged->SetXTitle("p_{T trigger}");
1182       outputContainer->Add(fhPtHbpZTUeLeftCharged) ;
1183
1184       fhPtHbpZTUeRightCharged  =
1185       new TH2F("hHbpZTUeChargedRight","#xi = ln(1/z_{T}) with charged UE right side of trigger",
1186                nptbins,ptmin,ptmax,200,0.,10.);
1187       fhPtHbpZTUeRightCharged->SetYTitle("ln(1/z_{T})");
1188       fhPtHbpZTUeRightCharged->SetXTitle("p_{T trigger}");
1189       outputContainer->Add(fhPtHbpZTUeRightCharged) ;
1190
1191     }
1192   }  //Correlation with charged hadrons
1193
1194   //Correlation with neutral hadrons
1195   if(fNeutralCorr)
1196   {
1197     fhDeltaPhiDeltaEtaNeutral  = new TH2F
1198     ("hDeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
1199      ndeltaphibins ,deltaphimin,deltaphimax, ndeltaetabins ,deltaetamin,deltaetamax);
1200     fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");
1201     fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
1202
1203     fhPhiNeutral  = new TH2F
1204     ("hPhiNeutral","#phi_{#pi^{0}}  vs p_{T #pi^{0}}",
1205      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1206     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
1207     fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
1208
1209     fhEtaNeutral  = new TH2F
1210     ("hEtaNeutral","#eta_{#pi^{0}}  vs p_{T #pi^{0}}",
1211      nptbins,ptmin,ptmax,netabins,etamin,etamax);
1212     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
1213     fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");
1214
1215     fhDeltaPhiNeutral  = new TH2F
1216     ("hDeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",
1217      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1218     fhDeltaPhiNeutral->SetYTitle("#Delta #phi");
1219     fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");
1220
1221     fhDeltaPhiNeutralPt  = new TH2F
1222     ("hDeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
1223      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1224     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");
1225     fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
1226
1227     fhDeltaPhiUeNeutralPt  = new TH2F
1228     ("hDeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",
1229      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1230     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");
1231     fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");
1232
1233     fhDeltaEtaNeutral  = new TH2F
1234     ("hDeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",
1235      nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
1236     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
1237     fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");
1238
1239     fhXENeutral  =
1240     new TH2F("hXENeutral","x_{E} for #pi^{0} associated",
1241              nptbins,ptmin,ptmax,200,0.,2.);
1242     fhXENeutral->SetYTitle("x_{E}");
1243     fhXENeutral->SetXTitle("p_{T trigger}");
1244
1245     fhXEUeNeutral  =
1246     new TH2F("hXEUeNeutral","x_{E} for #pi^{0} associated",
1247              nptbins,ptmin,ptmax,200,0.,2.);
1248     fhXEUeNeutral->SetYTitle("x_{E}");
1249     fhXEUeNeutral->SetXTitle("p_{T trigger}");
1250
1251     fhPtHbpXENeutral  =
1252     new TH2F("hHbpXENeutral","#xi = ln(1/x_{E})for #pi^{0} associated",
1253              nptbins,ptmin,ptmax,200,0.,10.);
1254     fhPtHbpXENeutral->SetYTitle("ln(1/x_{E})");
1255     fhPtHbpXENeutral->SetXTitle("p_{T trigger}");
1256
1257     fhPtHbpXEUeNeutral  =
1258     new TH2F("hHbpXEUeNeutral","#xi = ln(1/x_{E}) for #pi^{0} associated",
1259              nptbins,ptmin,ptmax,200,0.,10.);
1260     fhPtHbpXEUeNeutral->SetYTitle("ln(1/x_{E})");
1261     fhPtHbpXEUeNeutral->SetXTitle("p_{T trigger}");
1262
1263     fhZTNeutral  =
1264     new TH2F("hZTNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger} for #pi^{0} associated",
1265              nptbins,ptmin,ptmax,200,0.,2.);
1266     fhZTNeutral->SetYTitle("z_{trigger #pi^{0}}");
1267     fhZTNeutral->SetXTitle("p_{T trigger}");
1268
1269     fhZTUeNeutral  =
1270     new TH2F("hZTUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger} for #pi^{0} associated",
1271              nptbins,ptmin,ptmax,200,0.,2.);
1272     fhZTUeNeutral->SetYTitle("z_{trigger #pi^{0}}");
1273     fhZTUeNeutral->SetXTitle("p_{T trigger}");
1274
1275     fhPtHbpZTNeutral  =
1276     new TH2F("hHbpZTNeutral","#xi = ln(1/x_{E}) for #pi^{0} associated",
1277              nptbins,ptmin,ptmax,200,0.,10.);
1278     fhPtHbpZTNeutral->SetYTitle("ln(1/z_{T})");
1279     fhPtHbpZTNeutral->SetXTitle("p_{T trigger}");
1280
1281     fhPtHbpZTUeNeutral  =
1282     new TH2F("hHbpZTUeNeutral","#xi = ln(1/x_{E}) for #pi^{0} associated",
1283              nptbins,ptmin,ptmax,200,0.,10.);
1284     fhPtHbpXEUeNeutral->SetYTitle("ln(1/z_{T})");
1285     fhPtHbpXEUeNeutral->SetXTitle("p_{T trigger}");
1286
1287     outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
1288     outputContainer->Add(fhPhiNeutral) ;
1289     outputContainer->Add(fhEtaNeutral) ;
1290     outputContainer->Add(fhDeltaPhiNeutral) ;
1291     outputContainer->Add(fhDeltaPhiNeutralPt) ;
1292     outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
1293     outputContainer->Add(fhDeltaEtaNeutral) ;
1294     outputContainer->Add(fhXENeutral) ;
1295     outputContainer->Add(fhXEUeNeutral) ;
1296     outputContainer->Add(fhPtHbpXENeutral) ;
1297     outputContainer->Add(fhPtHbpXEUeNeutral) ;
1298     outputContainer->Add(fhZTNeutral) ;
1299     outputContainer->Add(fhZTUeNeutral) ;
1300     outputContainer->Add(fhPtHbpZTNeutral) ;
1301     outputContainer->Add(fhPtHbpZTUeNeutral) ;
1302
1303     if(fPi0Trigger || fDecayTrigger)
1304     {
1305       fhDeltaPhiDecayNeutral  = new TH2F
1306       ("hDeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",
1307        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1308       fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");
1309       fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");
1310
1311       fhXEDecayNeutral  =
1312       new TH2F("hXEDecayNeutral","x_{E} for decay trigger",
1313                nptbins,ptmin,ptmax,200,0.,2.);
1314       fhXEDecayNeutral->SetYTitle("x_{E}");
1315       fhXEDecayNeutral->SetXTitle("p_{T decay}");
1316
1317       fhZTDecayNeutral  =
1318       new TH2F("hZTDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",
1319                nptbins,ptmin,ptmax,200,0.,2.);
1320       fhZTDecayNeutral->SetYTitle("z_{h^{0}}");
1321       fhZTDecayNeutral->SetXTitle("p_{T decay}");
1322
1323       outputContainer->Add(fhDeltaPhiDecayNeutral) ;
1324       outputContainer->Add(fhXEDecayNeutral) ;
1325       outputContainer->Add(fhZTDecayNeutral) ;
1326
1327     }
1328
1329     if(fMakeSeveralUE)
1330     {
1331       fhDeltaPhiUeLeftNeutral  = new TH2F
1332       ("hDeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",
1333        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1334       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");
1335       fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
1336       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
1337
1338       fhDeltaPhiUeRightNeutral  = new TH2F
1339       ("hDeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",
1340        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1341       fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");
1342       fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");
1343       outputContainer->Add(fhDeltaPhiUeRightNeutral) ;
1344
1345       fhXEUeLeftNeutral  =
1346       new TH2F("hXEUeNeutralLeft","x_{E} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
1347                nptbins,ptmin,ptmax,140,0.,2.);
1348       fhXEUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1349       fhXEUeLeftNeutral->SetXTitle("p_{T trigger}");
1350       outputContainer->Add(fhXEUeLeftNeutral) ;
1351
1352       fhXEUeRightNeutral  =
1353       new TH2F("hXEUeNeutralRight","x_{E} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
1354                nptbins,ptmin,ptmax,200,0.,2.);
1355       fhXEUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1356       fhXEUeRightNeutral->SetXTitle("p_{T trigger}");
1357       outputContainer->Add(fhXEUeRightNeutral) ;
1358
1359       fhPtHbpXEUeLeftNeutral  =
1360       new TH2F("hHbpXEUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",
1361                nptbins,ptmin,ptmax,200,0.,10.);
1362       fhPtHbpXEUeLeftNeutral->SetYTitle("ln(1/x_{E})");
1363       fhPtHbpXEUeLeftNeutral->SetXTitle("p_{T trigger}");
1364       outputContainer->Add(fhPtHbpXEUeLeftNeutral) ;
1365
1366       fhPtHbpXEUeRightNeutral  =
1367       new TH2F("hHbpXEUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",
1368                nptbins,ptmin,ptmax,200,0.,10.);
1369       fhPtHbpXEUeRightNeutral->SetYTitle("ln(1/x_{E})");
1370       fhPtHbpXEUeRightNeutral->SetXTitle("p_{T trigger}");
1371       outputContainer->Add(fhPtHbpXEUeRightNeutral) ;
1372
1373       fhZTUeLeftNeutral  =
1374       new TH2F("hZTUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",
1375                nptbins,ptmin,ptmax,140,0.,2.);
1376       fhZTUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1377       fhZTUeLeftNeutral->SetXTitle("p_{T trigger}");
1378       outputContainer->Add(fhZTUeLeftNeutral) ;
1379
1380       fhZTUeRightNeutral  =
1381       new TH2F("hZTUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",
1382                nptbins,ptmin,ptmax,200,0.,2.);
1383       fhZTUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");
1384       fhZTUeRightNeutral->SetXTitle("p_{T trigger}");
1385       outputContainer->Add(fhZTUeRightNeutral) ;
1386
1387       fhPtHbpZTUeLeftNeutral  =
1388       new TH2F("hHbpZTUeNeutralLeft","#xi = ln(1/z_{T}) with neutral UE left side of trigger",
1389                nptbins,ptmin,ptmax,200,0.,10.);
1390       fhPtHbpZTUeLeftNeutral->SetYTitle("ln(1/z_{T})");
1391       fhPtHbpZTUeLeftNeutral->SetXTitle("p_{T trigger}");
1392       outputContainer->Add(fhPtHbpZTUeLeftNeutral) ;
1393
1394       fhPtHbpZTUeRightNeutral  =
1395       new TH2F("hHbpZTUeNeutralRight","#xi = ln(1/z_{T}) with neutral UE right side of trigger",
1396                nptbins,ptmin,ptmax,200,0.,10.);
1397       fhPtHbpZTUeRightNeutral->SetYTitle("ln(1/z_{T})");
1398       fhPtHbpZTUeRightNeutral->SetXTitle("p_{T trigger}");
1399       outputContainer->Add(fhPtHbpZTUeRightNeutral) ;
1400
1401     }
1402
1403   }//Correlation with neutral hadrons
1404
1405   //if data is MC, fill more histograms
1406   if(IsDataMC())
1407   {
1408     fh2phiLeadingParticle=new TH2F("h2phiLeadingParticle","#phi resolustion for trigger particles",nptbins,ptmin,ptmax,100,-1,1);
1409     fh2phiLeadingParticle->GetXaxis()->SetTitle("p_{T gen Leading} (GeV/c)");
1410     fh2phiLeadingParticle->GetYaxis()->SetTitle("(#phi_{rec}-#phi_{gen})/#phi_{gen}");
1411
1412     fhMCPtLeading  = new TH1F ("hMCPtLeading","MC : p_T distribution of leading particles", nptbins,ptmin,ptmax);
1413     fhMCPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");
1414
1415     fhMCEtaCharged  = new TH2F
1416     ("hMCEtaCharged","MC #eta_{h^{#pm}}  vs p_{T #pm}",
1417      nptbins,ptmin,ptmax,netabins,etamin,etamax);
1418     fhMCEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
1419     fhMCEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");
1420
1421     fhMCPhiCharged  = new TH2F
1422     ("hMCPhiCharged","#MC phi_{h^{#pm}}  vs p_{T #pm}",
1423      200,ptmin,ptmax,nphibins,phimin,phimax);
1424     fhMCPhiCharged->SetYTitle("MC #phi_{h^{#pm}} (rad)");
1425     fhMCPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");
1426
1427     fhMCDeltaPhiDeltaEtaCharged  = new TH2F
1428     ("hMCDeltaPhiDeltaEtaCharged","#MC phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
1429      140,-2.,5.,200,-2,2);
1430     fhMCDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
1431     fhMCDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1432
1433     fhMCDeltaEtaCharged  = new TH2F
1434     ("hMCDeltaEtaCharged","MC #eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger} and p_{T assoc}",
1435      nptbins,ptmin,ptmax,200,-2,2);
1436     fhMCDeltaEtaCharged->SetYTitle("#Delta #eta");
1437     fhMCDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");
1438
1439     fhMCDeltaPhiCharged  = new TH2F
1440     ("hMCDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
1441      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1442     fhMCDeltaPhiCharged->SetYTitle("#Delta #phi");
1443     fhMCDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
1444
1445     fhMCDeltaPhiChargedPt  = new TH2F
1446     ("hMCDeltaPhiChargedPt","MC #phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",
1447      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1448     fhMCDeltaPhiChargedPt->SetYTitle("#Delta #phi");
1449     fhMCDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");
1450
1451     fhMCPtXECharged  =
1452     new TH2F("hMCPtXECharged","x_{E}",
1453              nptbins,ptmin,ptmax,200,0.,2.);
1454     fhMCPtXECharged->SetYTitle("x_{E}");
1455     fhMCPtXECharged->SetXTitle("p_{T trigger}");
1456
1457     fhMCPtXEUeCharged  =
1458     new TH2F("hMCPtXEUeCharged","x_{E}",
1459              nptbins,ptmin,ptmax,200,0.,2.);
1460     fhMCPtXEUeCharged->SetYTitle("x_{E}");
1461     fhMCPtXEUeCharged->SetXTitle("p_{T trigger}");
1462
1463     fhMCPtHbpXECharged  =
1464     new TH2F("hMCHbpXECharged","MC #xi = ln(1/x_{E}) with charged hadrons",
1465              nptbins,ptmin,ptmax,200,0.,10.);
1466     fhMCPtHbpXECharged->SetYTitle("ln(1/x_{E})");
1467     fhMCPtHbpXECharged->SetXTitle("p_{T trigger}");
1468
1469     fhMCPtHbpXEUeCharged =
1470     new TH2F("hMCPtHbpXEUeCharged","#xi = ln(1/x_{E}) with charged hadrons,Underlying Event",
1471              nptbins,ptmin,ptmax,200,0.,10.);
1472     fhMCPtHbpXEUeCharged->SetYTitle("ln(1/x_{E})");
1473     fhMCPtHbpXEUeCharged->SetXTitle("p_{T trigger}");
1474
1475     fhMCUePart  =
1476     new TH1F("hMCUePart","MC UE particles distribution vs pt trig",
1477              nptbins,ptmin,ptmax);
1478     fhMCUePart->SetYTitle("dNch");
1479     fhMCUePart->SetXTitle("p_{T trigger}");
1480
1481     fhMCPtZTCharged  =
1482     new TH2F("hMCPtZTCharged","z_{T}",
1483              nptbins,ptmin,ptmax,200,0.,2.);
1484     fhMCPtZTCharged->SetYTitle("z_{T}");
1485     fhMCPtZTCharged->SetXTitle("p_{T trigger}");
1486
1487     fhMCPtHbpZTCharged  =
1488     new TH2F("hMCHbpZTCharged","MC #xi = ln(1/z_{T}) with charged hadrons",
1489              nptbins,ptmin,ptmax,200,0.,10.);
1490     fhMCPtHbpZTCharged->SetYTitle("ln(1/z_{T})");
1491     fhMCPtHbpZTCharged->SetXTitle("p_{T trigger}");
1492
1493     fhMCPtTrigPout  =
1494     new TH2F("hMCPtTrigPout","AOD MC Pout with triggers",
1495              nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
1496     fhMCPtTrigPout->SetYTitle("p_{out} (GeV/c)");
1497     fhMCPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)");
1498
1499     fhMCPtAssocDeltaPhi  =
1500     new TH2F("hMCPtAssocDeltaPhi","AOD MC delta phi with associated charged hadrons",
1501              nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1502     fhMCPtAssocDeltaPhi->SetYTitle("#Delta #phi");
1503     fhMCPtAssocDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
1504
1505     outputContainer->Add(fh2phiLeadingParticle);
1506     outputContainer->Add(fhMCPtLeading);
1507     outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged);
1508     outputContainer->Add(fhMCPhiCharged) ;
1509     outputContainer->Add(fhMCEtaCharged) ;
1510     outputContainer->Add(fhMCDeltaEtaCharged) ;
1511     outputContainer->Add(fhMCDeltaPhiCharged) ;
1512
1513     outputContainer->Add(fhMCDeltaPhiChargedPt) ;
1514     outputContainer->Add(fhMCPtXECharged) ;
1515     outputContainer->Add(fhMCPtXEUeCharged) ;
1516     outputContainer->Add(fhMCPtZTCharged) ;
1517     outputContainer->Add(fhMCPtHbpXECharged) ;
1518     outputContainer->Add(fhMCPtHbpXEUeCharged);
1519     outputContainer->Add(fhMCUePart);
1520     outputContainer->Add(fhMCPtHbpZTCharged) ;
1521     outputContainer->Add(fhMCPtTrigPout) ;
1522     outputContainer->Add(fhMCPtAssocDeltaPhi) ;
1523   } //for MC histogram
1524
1525   if(DoOwnMix())
1526   {
1527     //create event containers
1528
1529     if(!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForTracksExists()))
1530     {
1531
1532       Int_t nvz = GetNZvertBin();
1533       Int_t nrp = GetNRPBin();
1534       Int_t nce = GetNCentrBin();
1535
1536       fListMixEvents= new TList*[nvz*nrp*nce] ;
1537
1538       for( Int_t ice = 0 ; ice < nce ; ice++ )
1539       {
1540         for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
1541         {
1542           for( Int_t irp = 0 ; irp < nrp ; irp++ )
1543           {
1544             Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
1545
1546             //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
1547             //       ic,iz, irp, bin);
1548
1549             fListMixEvents[bin] = new TList() ;
1550             fListMixEvents[bin]->SetOwner(kFALSE);
1551           }
1552         }
1553       }
1554     }
1555
1556     //Init the list in the reader if not done previously
1557     if(!GetReader()->ListWithMixedEventsForTracksExists() && fUseMixStoredInReader)
1558     {
1559       //printf("%s : Set the list of events \n",GetInputAODName().Data());
1560       GetReader()->SetListWithMixedEventsForTracks(fListMixEvents);
1561     }
1562
1563     fhEventBin=new TH1I("hEventBin","Number of real events per bin(cen,vz,rp)",
1564                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
1565                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
1566     fhEventBin->SetXTitle("bin");
1567     outputContainer->Add(fhEventBin) ;
1568
1569     fhEventMixBin=new TH1I("hEventMixBin","Number of events  per bin(cen,vz,rp)",
1570                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
1571                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
1572     fhEventMixBin->SetXTitle("bin");
1573     outputContainer->Add(fhEventMixBin) ;
1574
1575     fhNtracksAll=new TH1F("hNtracksAll","Number of tracks w/o event trigger",2000,0,2000);
1576     outputContainer->Add(fhNtracksAll);
1577
1578     fhNtracksTrigger=new TH1F("hNtracksTriggerEvent","Number of tracks w/ event trigger",2000,0,2000);
1579     outputContainer->Add(fhNtracksTrigger);
1580
1581     fhNtracksMB=new TH1F("hNtracksMBEvent","Number of tracks w/ event trigger kMB",2000,0,2000);
1582     outputContainer->Add(fhNtracksMB);
1583
1584     fhMixDeltaPhiCharged  = new TH2F
1585     ("hMixDeltaPhiCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",
1586      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1587     fhMixDeltaPhiCharged->SetYTitle("#Delta #phi");
1588     fhMixDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");
1589     outputContainer->Add(fhMixDeltaPhiCharged);
1590
1591     fhMixDeltaPhiDeltaEtaCharged  = new TH2F
1592     ("hMixDeltaPhiDeltaEtaCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
1593      ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
1594     fhMixDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");
1595     fhMixDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1596     outputContainer->Add(fhMixDeltaPhiDeltaEtaCharged);
1597
1598     fhMixDeltaPhiChargedAssocPtBin         = new TH2F*[fNAssocPtBins] ;
1599     fhMixDeltaPhiChargedAssocPtBinDEta08   = new TH2F*[fNAssocPtBins] ;
1600     fhMixDeltaPhiDeltaEtaChargedAssocPtBin = new TH2F*[fNAssocPtBins] ;
1601
1602     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
1603     {
1604       fhMixDeltaPhiChargedAssocPtBin[i] = new TH2F(Form("hMixDeltaPhiChargedAssocPtBin%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1605                                                    Form("Mixed event #Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1606                                                    nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
1607       fhMixDeltaPhiChargedAssocPtBin[i]->SetXTitle("p_{T trigger}");
1608       fhMixDeltaPhiChargedAssocPtBin[i]->SetYTitle("#Delta #phi");
1609
1610       fhMixDeltaPhiChargedAssocPtBinDEta08[i] = new TH2F(Form("hMixDeltaPhiDeltaEta0.8ChargedAssocPtBinDEta08%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1611                                                          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]),
1612                                                          nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
1613       fhMixDeltaPhiChargedAssocPtBinDEta08[i]->SetXTitle("p_{T trigger}");
1614       fhMixDeltaPhiChargedAssocPtBinDEta08[i]->SetYTitle("#Delta #phi");
1615
1616       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i] = new TH2F(Form("hMixDeltaPhiDeltaEtaChargedAssocPtBin%2.1f_%2.1f", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1617                                                            Form("Mixed event #Delta #phi vs p_{T trigger} for associated p_{T} bin [%2.1f,%2.1f]", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1]),
1618                                                            ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
1619       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]->SetXTitle("#Delta #phi");
1620       fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]->SetYTitle("#Delta #eta");
1621
1622       outputContainer->Add(fhMixDeltaPhiChargedAssocPtBin[i]);
1623       outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta08[i]);
1624       outputContainer->Add(fhMixDeltaPhiDeltaEtaChargedAssocPtBin[i]);
1625
1626     }
1627   }
1628
1629   return outputContainer;
1630
1631 }
1632
1633 //_________________________________________________________________________________________________
1634 Bool_t AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum(const AliAODPWG4Particle* trigger,
1635                                                              TLorentzVector & mom1,
1636                                                              TLorentzVector & mom2)
1637 {
1638   // Get the momentum of the pi0/eta assigned decay photons
1639   // In case of pi0/eta trigger, we may want to check their decay correlation,
1640   // get their decay children
1641
1642   Int_t indexPhoton1 = trigger->GetCaloLabel(0);
1643   Int_t indexPhoton2 = trigger->GetCaloLabel(1);
1644   Float_t ptTrig     = trigger->Pt();
1645
1646   if(indexPhoton1!=-1 || indexPhoton2!=-1) return kFALSE;
1647
1648   if(GetDebug() > 1)
1649     printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
1650
1651   TObjArray * clusters  = 0x0 ;
1652   if(trigger->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
1653   else                                clusters = GetPHOSClusters()  ;
1654
1655   for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++)
1656   {
1657     AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));
1658     if(photon->GetID()==indexPhoton1)
1659     {
1660       photon->GetMomentum(mom1,GetVertex(0)) ;
1661       if(ptTrig) fhPtPi0DecayRatio->Fill(ptTrig, mom1.Pt()/ptTrig);
1662     }
1663     if(photon->GetID()==indexPhoton2)
1664     {
1665       photon->GetMomentum(mom1,GetVertex(0)) ;
1666       if(ptTrig > 0) fhPtPi0DecayRatio->Fill(ptTrig, mom2.Pt()/ptTrig);
1667     }
1668
1669     if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - Photon1 = %f, Photon2 = %f \n", mom1.Pt(), mom2.Pt());
1670
1671   } //cluster loop
1672
1673   return kTRUE;
1674
1675 }
1676
1677 //____________________________________________________
1678 void AliAnaParticleHadronCorrelation::InitParameters()
1679 {
1680
1681   //Initialize the parameters of the analysis.
1682   SetInputAODName("Particle");
1683   SetAODObjArrayName("Hadrons");
1684   AddToHistogramsName("AnaHadronCorr_");
1685
1686   SetPtCutRange(0.,300);
1687   fDeltaPhiMinCut       = 1.5 ;
1688   fDeltaPhiMaxCut       = 4.5 ;
1689   fSelectIsolated       = kFALSE;
1690   fMakeSeveralUE        = kFALSE;
1691   fUeDeltaPhiMinCut     = 1. ;
1692   fUeDeltaPhiMaxCut     = 1.5 ;
1693
1694   fNeutralCorr          = kFALSE ;
1695   fPi0Trigger           = kFALSE ;
1696   fDecayTrigger         = kFALSE ;
1697   fHMPIDCorrelation     = kFALSE ;
1698
1699   fMakeAbsoluteLeading  = kTRUE;
1700   fMakeNearSideLeading  = kFALSE;
1701
1702   fNAssocPtBins         = 9   ;
1703   fAssocPtBinLimit[0]   = 0.2 ;
1704   fAssocPtBinLimit[1]   = 2.0 ;
1705   fAssocPtBinLimit[2]   = 4.0 ;
1706   fAssocPtBinLimit[3]   = 6.0 ;
1707   fAssocPtBinLimit[4]   = 8.0 ;
1708   fAssocPtBinLimit[5]   = 10. ;
1709   fAssocPtBinLimit[6]   = 12. ;
1710   fAssocPtBinLimit[7]   = 15. ;
1711   fAssocPtBinLimit[8]   = 25. ;
1712   fAssocPtBinLimit[9]   = 50. ;
1713
1714   fUseMixStoredInReader = kTRUE;
1715
1716   fM02MinCut   = -1 ;
1717   fM02MaxCut   = -1 ;
1718
1719 }
1720
1721 //__________________________________________________________
1722 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()
1723 {
1724   //Particle-Hadron Correlation Analysis, fill AODs
1725
1726   if(!GetInputAODBranch())
1727   {
1728     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1729     abort();
1730   }
1731
1732   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation"))
1733   {
1734     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());
1735     abort();
1736   }
1737
1738   if(GetDebug() > 1)
1739   {
1740     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");
1741     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
1742     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n",   GetCTSTracks()    ->GetEntriesFast());
1743     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
1744     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n",  GetPHOSClusters() ->GetEntriesFast());
1745   }
1746
1747   //Get the vertex and check it is not too large in z
1748   Double_t v[3] = {0,0,0}; //vertex ;
1749   GetReader()->GetVertex(v);
1750   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
1751
1752   // Fill the pool with tracks if requested
1753   if(DoOwnMix()) FillChargedEventMixPool();
1754
1755   //Loop on stored AOD particles, find leading trigger
1756   Double_t ptTrig      = fMinTriggerPt ;
1757   fLeadingTriggerIndex = -1 ;
1758   Int_t    naod        = GetInputAODBranch()->GetEntriesFast() ;
1759   for(Int_t iaod = 0; iaod < naod ; iaod++)
1760   {
1761     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1762
1763     // Vertex cut in case of mixing
1764     Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
1765     if(check ==  0) continue;
1766     if(check == -1) return;
1767
1768     // find the leading particles with highest momentum
1769     if (particle->Pt() > ptTrig)
1770     {
1771       ptTrig               = particle->Pt() ;
1772       fLeadingTriggerIndex = iaod ;
1773     }
1774   }// finish search of leading trigger particle
1775
1776
1777   //Do correlation with leading particle
1778   if(fLeadingTriggerIndex >= 0)
1779   {
1780
1781     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(fLeadingTriggerIndex));
1782
1783     //check if the particle is isolated or if we want to take the isolation into account
1784     if(OnlyIsolated() && !particle->IsIsolated()) return;
1785
1786     //Make correlation with charged hadrons
1787     Bool_t okcharged = kTRUE;
1788     Bool_t okneutral = kTRUE;
1789     if(GetReader()->IsCTSSwitchedOn() )
1790       okcharged = MakeChargedCorrelation(particle, GetCTSTracks(),kFALSE);
1791
1792     TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
1793     if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)
1794       okneutral = MakeNeutralCorrelation(particle, pi0list,kFALSE);
1795
1796   }//Correlate leading
1797
1798   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");
1799
1800 }
1801
1802 //_________________________________________________________________
1803 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
1804 {
1805   //Particle-Hadron Correlation Analysis, fill histograms
1806
1807   if(!GetInputAODBranch())
1808   {
1809     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1810     abort();
1811   }
1812
1813   if(GetDebug() > 1)
1814   {
1815     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
1816     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
1817   }
1818
1819   //Get the vertex and check it is not too large in z
1820   Double_t v[3] = {0,0,0}; //vertex ;
1821   GetReader()->GetVertex(v);
1822   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;
1823
1824   //Loop on stored AOD particles, find leading
1825   Double_t ptTrig    = fMinTriggerPt;
1826   if(fLeadingTriggerIndex < 0)
1827   {
1828     //Search leading if not done before
1829     Int_t    naod      = GetInputAODBranch()->GetEntriesFast() ;
1830     for(Int_t iaod = 0; iaod < naod ; iaod++)
1831     {    //loop on input trigger AOD file
1832       AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1833
1834       // Vertex cut in case of mixing
1835       Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
1836       if(check ==  0) continue;
1837       if(check == -1) return;
1838
1839       //check if the particle is isolated or if we want to take the isolation into account
1840       if(OnlyIsolated() && !particle->IsIsolated()) continue;
1841
1842       //find the leading particles with highest momentum
1843       if (particle->Pt() > ptTrig)
1844       {
1845         ptTrig               = particle->Pt() ;
1846         fLeadingTriggerIndex = iaod ;
1847       }
1848
1849     }// Finish search of leading trigger particle
1850   }// Search leading if not done before
1851
1852   if(fLeadingTriggerIndex >= 0 )
1853   { //using trigger particle to do correlations
1854
1855     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(fLeadingTriggerIndex));
1856
1857     // check if it was a calorimeter cluster and if the SS cut was requested, if so, apply it
1858     Int_t clID1  = particle->GetCaloLabel(0) ;
1859     Int_t clID2  = particle->GetCaloLabel(1) ; // for photon clusters should not be set.
1860     //printf("Leading for for %s: id1 %d, id2 %d, min %f, max %f, det %s\n",
1861     //       GetInputAODName().Data(),clID1,clID2,fM02MinCut,fM02MaxCut,(particle->GetDetector()).Data());
1862
1863     if(clID1 > 0 && clID2 < 0 && fM02MaxCut > 0 && fM02MinCut > 0)
1864     {
1865       Int_t iclus = -1;
1866       TObjArray* clusters = 0x0;
1867       if     (particle->GetDetector() == "EMCAL") clusters = GetEMCALClusters();
1868       else if(particle->GetDetector() == "PHOS" ) clusters = GetPHOSClusters();
1869
1870       if(clusters)
1871       {
1872         AliVCluster *cluster = FindCluster(clusters,clID1,iclus);
1873         Float_t m02 = cluster->GetM02();
1874         //printf("\t Check m02 = %2.2f\n",m02);
1875         if(m02 > fM02MaxCut || m02 < fM02MinCut)
1876         {
1877           //printf("\t \t Not accepted\n");
1878           return;
1879         }
1880       }
1881     }
1882
1883     // Check if trigger is in fiducial region
1884     if(IsFiducialCutOn())
1885     {
1886       Bool_t in = GetFiducialCut()->IsInFiducialCut(*particle->Momentum(),particle->GetDetector()) ;
1887       if(! in ) return ;
1888     }
1889
1890     // Check if the particle is isolated or if we want to take the isolation into account
1891     if(OnlyIsolated() && !particle->IsIsolated()) return;
1892
1893     // Make correlation with charged hadrons
1894     Bool_t okcharged = kTRUE;
1895     Bool_t okneutral = kTRUE;
1896     if(GetReader()->IsCTSSwitchedOn() )
1897     {
1898       okcharged = MakeChargedCorrelation(particle, GetCTSTracks(),kTRUE);
1899       if(IsDataMC())
1900       {
1901         MakeMCChargedCorrelation(particle);
1902       }
1903     }
1904
1905     TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists
1906     if(fNeutralCorr && pi0list)
1907     {
1908       if(pi0list->GetEntriesFast() > 0)
1909         okneutral = MakeNeutralCorrelation(particle, pi0list,kTRUE);
1910     }
1911
1912     // Fill leading particle histogram if correlation went well and
1913     // no problem was found, like not absolute leading, or bad vertex in mixing.
1914     if(okcharged && okneutral)
1915     {
1916       Float_t pt = particle->Pt();
1917       fhPtLeading->Fill(pt);
1918
1919       Float_t phi = particle->Phi();
1920       if(phi<0)phi+=TMath::TwoPi();
1921       fhPhiLeading->Fill(pt, phi);
1922
1923       fhEtaLeading->Fill(pt, particle->Eta());
1924       //printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Leading particle : pt %f, eta %f, phi %f\n",particle->Pt(),particle->Eta(),phi);
1925
1926       Float_t cen = GetEventCentrality();
1927       Float_t ep  = GetEventPlaneAngle();
1928
1929       fhPtLeadingCentrality        ->Fill(pt,cen);
1930       fhPtLeadingEventPlane        ->Fill(pt,ep);
1931       fhLeadingEventPlaneCentrality->Fill(cen,ep);
1932
1933     }//ok charged && neutral
1934   }//Aod branch loop
1935
1936   //Reinit for next event
1937   fLeadingTriggerIndex = -1;
1938
1939   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
1940 }
1941
1942 //___________________________________________________________________________________________________________
1943 Bool_t  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle,
1944                                                                 const TObjArray* pl, const Bool_t bFillHisto)
1945 {
1946   // Charged Hadron Correlation Analysis
1947   if(GetDebug() > 1)
1948     printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
1949
1950   Float_t phiTrig = aodParticle->Phi();
1951   Float_t etaTrig = aodParticle->Eta();
1952   Float_t ptTrig  = aodParticle->Pt();
1953   Bool_t   decay  = aodParticle->IsTagged();
1954
1955   Float_t pt       = -100. ;
1956   Float_t zT       = -100. ;
1957   Float_t xE       = -100. ;
1958   Float_t hbpXE    = -100. ;
1959   Float_t hbpZT    = -100. ;
1960   Float_t phi      = -100. ;
1961   Float_t eta      = -100. ;
1962   Float_t pout     = -100. ;
1963   Float_t deltaPhi = -100. ;
1964
1965   TVector3 p3;
1966   TLorentzVector photonMom ;
1967   TObjArray * reftracks = 0x0;
1968   Int_t nrefs           = 0;
1969   Int_t nTracks         = GetCTSTracks()->GetEntriesFast() ;
1970
1971   // Mixed event settings
1972   Int_t evtIndex11   = -1 ; // cluster trigger or pi0 trigger
1973   Int_t evtIndex12   = -1 ; // pi0 trigger
1974   Int_t evtIndex13   = -1 ; // charged trigger
1975
1976   Double_t v[3]      = {0,0,0}; //vertex ;
1977   GetReader()->GetVertex(v);
1978
1979   if (GetMixedEvent())
1980   {
1981     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
1982     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
1983     evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
1984   }
1985
1986   // In case of pi0/eta trigger, we may want to check their decay correlation,
1987   // get their decay children
1988   TLorentzVector decayMom1;
1989   TLorentzVector decayMom2;
1990   Bool_t decayFound = kFALSE;
1991   if(fPi0Trigger && bFillHisto) decayFound = GetDecayPhotonMomentum(aodParticle,decayMom1, decayMom2);
1992
1993   //-----------------------------------------------------------------------
1994   //Track loop, select tracks with good pt, phi and fill AODs or histograms
1995   //-----------------------------------------------------------------------
1996
1997   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ )
1998   {
1999     AliVTrack * track = (AliVTrack *) (pl->At(ipr)) ;
2000
2001     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
2002     p3.SetXYZ(mom[0],mom[1],mom[2]);
2003     pt   = p3.Pt();
2004     eta  = p3.Eta();
2005     phi  = p3.Phi() ;
2006     if(phi < 0) phi+=TMath::TwoPi();
2007
2008     //Select only hadrons in pt range
2009     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
2010
2011     //remove trigger itself for correlation when use charged triggers
2012     if( track->GetID() == aodParticle->GetTrackLabel(0) || track->GetID() == aodParticle->GetTrackLabel(1) ||
2013         track->GetID() == aodParticle->GetTrackLabel(2) || track->GetID() == aodParticle->GetTrackLabel(3)   )
2014       continue ;
2015
2016     //jump out this event if near side associated particle pt larger than trigger
2017     if (fMakeNearSideLeading)
2018     {
2019       if(pt > ptTrig && TMath::Abs(phi-phiTrig) < TMath::PiOver2())  return kFALSE;
2020     }
2021     //jump out this event if there is any other particle with pt larger than trigger
2022     else if(fMakeAbsoluteLeading)
2023     {
2024       if(pt > ptTrig)  return kFALSE;
2025     }
2026
2027     //Only for mixed event
2028     Int_t evtIndex2 = 0 ;
2029     if (GetMixedEvent())
2030     {
2031       evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
2032       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
2033         continue ;
2034       //vertex cut
2035       if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
2036         return kFALSE;
2037     }
2038
2039     // Fill Histograms
2040     if(bFillHisto)
2041     {
2042
2043       if(GetDebug() > 2 )
2044         printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
2045
2046       // Set the pt associated bin for the defined bins
2047       Int_t assocBin   = -1;
2048       for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2049       {
2050         if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i;
2051       }
2052
2053       // Azimuthal Angle
2054       // calculate deltaPhi for later, shift when needed
2055       FillChargedAngularCorrelationHistograms(pt,  ptTrig,  assocBin, phi, phiTrig,  deltaPhi,
2056                                               eta, etaTrig, decay, track->GetHMPIDsignal(),nTracks);
2057
2058       // Imbalance zT/xE/pOut
2059       zT = pt/ptTrig ;
2060       if(zT > 0 ) hbpZT = TMath::Log(1./zT);
2061       else        hbpZT =-100;
2062
2063       xE   =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
2064       //if(xE <0.)xE =-xE;
2065       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
2066       else        hbpXE =-100;
2067
2068       pout = pt*TMath::Sin(deltaPhi) ;
2069
2070       //delta phi cut for momentum imbalance correlation
2071       if      ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   )
2072       {
2073
2074         FillChargedMomentumImbalanceHistograms(ptTrig, pt, xE, hbpXE, zT, hbpZT, pout, deltaPhi,
2075                                                nTracks, track->Charge(), assocBin, decay);
2076
2077       }
2078       if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
2079       { //UE study
2080
2081         FillChargedUnderlyingEventHistograms(ptTrig, pt, deltaPhi, nTracks);
2082
2083         fhUePart->Fill(ptTrig);
2084
2085       }
2086
2087       if(fPi0Trigger && decayFound)
2088         FillDecayPhotonCorrelationHistograms(pt, phi, decayMom1,decayMom2, kTRUE) ;
2089
2090       //several UE calculation
2091       if(fMakeSeveralUE) FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,xE,hbpXE,zT,hbpZT,deltaPhi);
2092
2093     } //Fill histogram
2094     else
2095     {
2096       nrefs++;
2097       if(nrefs==1)
2098       {
2099         reftracks = new TObjArray(0);
2100         TString trackname = Form("%s+Tracks", GetAODObjArrayName().Data());
2101         reftracks->SetName(trackname.Data());
2102         reftracks->SetOwner(kFALSE);
2103       }
2104
2105       reftracks->Add(track);
2106
2107     }//aod particle loop
2108   }// track loop
2109
2110   //Fill AOD with reference tracks, if not filling histograms
2111   if(!bFillHisto && reftracks)
2112   {
2113     aodParticle->AddObjArray(reftracks);
2114   }
2115
2116   //Own mixed event, add event and remove previous or fill the mixed histograms
2117   if(DoOwnMix() && bFillHisto)
2118   {
2119       MakeChargedMixCorrelation(aodParticle);
2120   }
2121
2122   return kTRUE;
2123
2124 }
2125
2126
2127 //_________________________________________________________________________________________________________
2128 void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
2129 {
2130   // Mix current trigger with tracks in another MB event
2131
2132   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Make trigger particle - charged hadron mixed event correlation \n");
2133
2134   if(GetMixedEvent()) return;  // This is not the mixed event from general mixing frame
2135
2136   // Get the event with similar caracteristics
2137   //printf("MakeChargedMixCorrelation for %s\n",GetInputAODName().Data());
2138
2139   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
2140
2141   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
2142
2143   if(!inputHandler) return;
2144
2145   if(!(inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask())) return;
2146
2147   // Get the pool, check if it exits
2148   Int_t eventBin = GetEventMixBin();
2149
2150   fhEventBin->Fill(eventBin);
2151
2152   //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
2153   if(eventBin < 0) return;
2154
2155   TList * pool = 0;
2156   if(fUseMixStoredInReader) pool = GetReader()->GetListWithMixedEventsForTracks(eventBin);
2157   else                      pool = fListMixEvents[eventBin];
2158
2159   if(!pool) return ;
2160
2161   Double_t ptTrig  = aodParticle->Pt();
2162   Double_t etaTrig = aodParticle->Eta();
2163   Double_t phiTrig = aodParticle->Phi();
2164   if(phiTrig < 0.) phiTrig+=TMath::TwoPi();
2165
2166   if(GetDebug() > 1)
2167     printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Pool bin %d size %d, leading trigger pt=%f, phi=%f, eta=%f\n",
2168            eventBin,pool->GetSize(), ptTrig,phiTrig,etaTrig);
2169
2170   Double_t ptAssoc  = -999.;
2171   Double_t phiAssoc = -999.;
2172   Double_t etaAssoc = -999.;
2173   Double_t deltaPhi = -999.;
2174   Double_t deltaEta = -999.;
2175
2176   // Get the clusters array, needed for isolation
2177   TObjArray * caloList    = 0x0; ;
2178   if      (aodParticle->GetDetector() == "PHOS" )
2179     caloList = GetPHOSClusters();
2180   else if (aodParticle->GetDetector() == "EMCAL")
2181     caloList = GetEMCALClusters();
2182
2183   //Start from first event in pool except if in this same event the pool was filled
2184   Int_t ev0 = 0;
2185   if(GetReader()->GetLastTracksMixedEvent() == GetEventNumber()) ev0 = 1;
2186
2187   for(Int_t ev=ev0; ev < pool->GetSize(); ev++)
2188   {
2189     TObjArray* bgTracks = static_cast<TObjArray*>(pool->At(ev));
2190
2191     // Check if the particle is isolated in the mixed event, it not, do not fill the histograms
2192     if(OnlyIsolated())
2193     {
2194       Int_t n=0; Int_t nfrac = 0; Bool_t isolated = kFALSE; Float_t coneptsum = 0;
2195       GetIsolationCut()->MakeIsolationCut(bgTracks,caloList,
2196                                           GetReader(), GetCaloPID(),
2197                                           kTRUE, aodParticle, GetAODObjArrayName(),
2198                                           n,nfrac,coneptsum, isolated);
2199       printf("Isolated? %d - cone %f, ptthres %f",
2200              isolated,GetIsolationCut()->GetConeSize(),GetIsolationCut()->GetPtThreshold());
2201       if(caloList)printf(" - n clus %d",caloList->GetEntriesFast());
2202       if(bgTracks)printf(" - n track %d", bgTracks->GetEntriesFast());
2203       printf("\n");
2204       if(!isolated) return ;
2205     }
2206
2207     fhEventMixBin->Fill(eventBin);
2208
2209     Int_t nTracks=bgTracks->GetEntriesFast();
2210     //printf("\t Read Pool event %d, nTracks %d\n",ev,nTracks);
2211
2212     for(Int_t j1 = 0;j1 <nTracks; j1++ )
2213     {
2214       AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(j1) ;
2215
2216       if(!track) continue;
2217
2218       ptAssoc  = track->Pt();
2219       etaAssoc = track->Eta();
2220       phiAssoc = track->Phi() ;
2221       if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
2222
2223       if(IsFiducialCutOn())
2224       {
2225         Bool_t in = GetFiducialCut()->IsInFiducialCut(*aodParticle->Momentum(),"CTS") ;
2226         if(!in) continue ;
2227       }
2228
2229       deltaPhi = phiTrig-phiAssoc;
2230       if(deltaPhi < -TMath::PiOver2())  deltaPhi+=TMath::TwoPi();
2231       if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
2232       deltaEta = etaTrig-etaAssoc;
2233
2234       if(GetDebug()>0)
2235         printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation(): deltaPhi= %f, deltaEta=%f\n",deltaPhi, deltaEta);
2236
2237       // Set the pt associated bin for the defined bins
2238       Int_t assocBin   = -1;
2239       for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2240       {
2241         if(ptAssoc > fAssocPtBinLimit[i] && ptAssoc < fAssocPtBinLimit[i+1]) assocBin= i;
2242       }
2243
2244       fhMixDeltaPhiCharged        ->Fill(ptTrig,  deltaPhi);
2245       fhMixDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
2246
2247       if(assocBin < 0) continue ; // this pt bin was not considered
2248
2249       if(TMath::Abs(deltaEta) > 0.8)
2250         fhMixDeltaPhiChargedAssocPtBinDEta08  [assocBin]->Fill(ptTrig,   deltaPhi);
2251
2252         fhMixDeltaPhiChargedAssocPtBin        [assocBin]->Fill(ptTrig,   deltaPhi);
2253         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[assocBin]->Fill(deltaPhi, deltaEta);
2254
2255     } // track loop
2256   } // mixed event loop
2257 }
2258
2259
2260 //________________________________________________________________________________________________________________
2261 Bool_t  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle,
2262                                                                 const TObjArray* pi0list, const Bool_t bFillHisto)
2263 {
2264   // Neutral Pion Correlation Analysis
2265   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",
2266                             pi0list->GetEntriesFast());
2267
2268   Int_t evtIndex11 = 0 ;
2269   Int_t evtIndex12 = 0 ;
2270   if (GetMixedEvent())
2271   {
2272     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
2273     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
2274   }
2275
2276   Float_t pt   = -100. ;
2277   Float_t zT   = -100. ;
2278   Float_t phi  = -100. ;
2279   Float_t eta  = -100. ;
2280   Float_t xE   = -100. ;
2281   Float_t hbpXE= -100. ;
2282   Float_t hbpZT= -100. ;
2283
2284   Float_t ptTrig  = aodParticle->Pt();
2285   Float_t phiTrig = aodParticle->Phi();
2286   Float_t etaTrig = aodParticle->Eta();
2287   Float_t deltaPhi= -100. ;
2288
2289   TLorentzVector photonMom ;
2290
2291   // In case of pi0/eta trigger, we may want to check their decay correlation,
2292   // get their decay children
2293   TLorentzVector decayMom1;
2294   TLorentzVector decayMom2;
2295   Bool_t decayFound = kFALSE;
2296   if(fPi0Trigger && bFillHisto) decayFound = GetDecayPhotonMomentum(aodParticle,decayMom1, decayMom2);
2297
2298   TObjArray * refpi0 = 0x0;
2299   Int_t nrefs        = 0;
2300
2301   //Loop on stored AOD pi0
2302
2303   Int_t naod = pi0list->GetEntriesFast();
2304   if(GetDebug() > 0)
2305     printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() -  aod branch entries %d\n", naod);
2306
2307   for(Int_t iaod = 0; iaod < naod ; iaod++)
2308   {
2309     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (pi0list->At(iaod));
2310
2311     Int_t evtIndex2 = 0 ;
2312     Int_t evtIndex3 = 0 ;
2313     if (GetMixedEvent())
2314     {
2315       evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
2316       evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
2317
2318       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 ||
2319           evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
2320         continue ;
2321     }
2322
2323     pt  = pi0->Pt();
2324
2325     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
2326
2327     //jump out this event if near side associated particle pt larger than trigger
2328     if (fMakeNearSideLeading)
2329     {
2330       if(pt > ptTrig && TMath::Abs(phi-phiTrig) < TMath::PiOver2())  return kFALSE;
2331     }
2332     //jump out this event if there is any other particle with pt larger than trigger
2333     else if(fMakeAbsoluteLeading)
2334     {
2335       if(pt > ptTrig)  return kFALSE;
2336     }
2337
2338     if(bFillHisto)
2339     {
2340       phi = pi0->Phi() ;
2341       eta = pi0->Eta() ;
2342
2343       FillNeutralAngularCorrelationHistograms(pt, ptTrig, phi, phiTrig, deltaPhi, eta, etaTrig);
2344
2345       zT  = pt/ptTrig ;
2346       xE  =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
2347
2348       //if(xE <0.)xE =-xE;
2349
2350       hbpXE = -100;
2351       hbpZT = -100;
2352
2353       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
2354       if(zT > 0 ) hbpZT = TMath::Log(1./zT);
2355
2356       if(fPi0Trigger && decayFound)
2357         FillDecayPhotonCorrelationHistograms(pt, phi, decayMom1,decayMom2,kFALSE) ;
2358
2359       //delta phi cut for correlation
2360       if( (deltaPhi > fDeltaPhiMinCut) && ( deltaPhi < fDeltaPhiMaxCut) )
2361       {
2362         fhDeltaPhiNeutralPt->Fill(pt,deltaPhi);
2363         fhXENeutral        ->Fill(ptTrig,xE);
2364         fhPtHbpXENeutral   ->Fill(ptTrig,hbpXE);
2365       }
2366       else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
2367       {
2368         fhDeltaPhiUeNeutralPt->Fill(pt,deltaPhi);
2369         fhXEUeNeutral        ->Fill(ptTrig,xE);
2370         fhPtHbpXEUeNeutral   ->Fill(ptTrig,hbpXE);
2371       }
2372
2373       //several UE calculation
2374       if(fMakeSeveralUE) FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,xE,hbpXE,zT,hbpZT,deltaPhi);
2375
2376           }
2377     else
2378     {
2379       nrefs++;
2380       if(nrefs==1)
2381       {
2382         refpi0 = new TObjArray(0);
2383         refpi0->SetName(GetAODObjArrayName()+"Pi0s");
2384         refpi0->SetOwner(kFALSE);
2385       }
2386       refpi0->Add(pi0);
2387     }//put references in trigger AOD
2388
2389     if(GetDebug() > 2 )
2390       printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
2391
2392   }//loop
2393
2394   return kTRUE;
2395 }
2396
2397 //_________________________________________________________________________________________________________
2398 void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
2399 {
2400   // Charged Hadron Correlation Analysis with MC information
2401
2402   if(GetDebug()>1)
2403     printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation() - Make trigger particle - charged hadron correlation in AOD MC level\n");
2404
2405   AliStack         * stack        = 0x0 ;
2406   TParticle        * primary      = 0x0 ;
2407   TClonesArray     * mcparticles0 = 0x0 ;
2408   TClonesArray     * mcparticles  = 0x0 ;
2409   AliAODMCParticle * aodprimary   = 0x0 ;
2410
2411   Double_t eprim   = 0 ;
2412   Double_t ptprim  = 0 ;
2413   Double_t phiprim = 0 ;
2414   Double_t etaprim = 0 ;
2415   Int_t    nTracks = 0 ;
2416   Int_t iParticle  = 0 ;
2417   Double_t charge  = 0.;
2418
2419
2420   if(GetReader()->ReadStack())
2421   {
2422     nTracks = GetMCStack()->GetNtrack() ;
2423   }
2424   else
2425   {
2426     nTracks = GetReader()->GetAODMCParticles()->GetEntriesFast() ;
2427   }
2428   //Int_t trackIndex[nTracks];
2429
2430   Int_t label= aodParticle->GetLabel();
2431   if(label < 0)
2432   {
2433     if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** bad label ***:  label %d \n", label);
2434     return;
2435   }
2436
2437
2438   if(GetReader()->ReadStack())
2439   {
2440     stack =  GetMCStack() ;
2441     if(!stack)
2442     {
2443       printf(" AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation- Stack not available, is the MC handler called? STOP\n");
2444       abort();
2445     }
2446
2447     nTracks=stack->GetNprimary();
2448     if(label >=  stack->GetNtrack())
2449     {
2450       if(GetDebug() > 2)  printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
2451       return ;
2452     }
2453
2454     primary = stack->Particle(label);
2455     if(!primary)
2456     {
2457       printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no primary ***:  label %d \n", label);
2458       return;
2459     }
2460
2461     if(primary)
2462     {
2463       eprim    = primary->Energy();
2464       ptprim   = primary->Pt();
2465       phiprim  = primary->Phi();
2466       etaprim  = primary->Eta();
2467
2468       if(ptprim < 0.01 || eprim < 0.01) return ;
2469
2470       for (iParticle = 0 ; iParticle <  nTracks ; iParticle++)
2471       {
2472         TParticle * particle = stack->Particle(iParticle);
2473         TLorentzVector momentum;
2474
2475         //keep only final state particles
2476         if(particle->GetStatusCode()!=1) continue ;
2477
2478         Int_t pdg = particle->GetPdgCode();
2479
2480         charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
2481
2482         particle->Momentum(momentum);
2483
2484         //---------- Charged particles ----------------------
2485         if(charge != 0)
2486         {
2487           //Particles in CTS acceptance
2488           Bool_t inCTS =  GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
2489
2490           if(TMath::Abs(pdg) == 11 && stack->Particle(particle->GetFirstMother())->GetPdgCode()==22) continue ;
2491
2492           if(inCTS)
2493           {
2494             if( label!=iParticle) // avoid trigger particle
2495             {
2496               if(!FillChargedMCCorrelationHistograms(particle->Pt(),particle->Phi(),particle->Eta(),ptprim,phiprim,etaprim)) return;
2497             }
2498           }// in CTS acceptance
2499         }// charged
2500       } //track loop
2501     } //when the leading particles could trace back to MC
2502   } //ESD MC
2503   else if(GetReader()->ReadAODMCParticles())
2504   {
2505     //Get the list of MC particles
2506     mcparticles0 = GetReader()->GetAODMCParticles(0);
2507     if(!mcparticles0) return;
2508     if(label >=mcparticles0->GetEntriesFast())
2509     {
2510       if(GetDebug() > 2)
2511         printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** large label ***:  label %d, n tracks %d \n", label,mcparticles0->GetEntriesFast());
2512       return;
2513     }
2514     //Get the particle
2515     aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
2516     if(!aodprimary)
2517     {
2518       printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation *** no AOD primary ***:  label %d \n", label);
2519       return;
2520     }
2521
2522
2523     if(aodprimary)
2524     {
2525       ptprim  = aodprimary->Pt();
2526       phiprim = aodprimary->Phi();
2527       etaprim = aodprimary->Eta();
2528       eprim   = aodprimary->E();
2529
2530       Bool_t lead = kFALSE;
2531
2532       if(ptprim < 0.01 || eprim < 0.01) return ;
2533
2534       mcparticles= GetReader()->GetAODMCParticles();
2535       for (Int_t i = 0; i < nTracks; i++)
2536       {
2537         AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(i);
2538
2539         if (!part->IsPhysicalPrimary()) continue;
2540         Int_t pdg = part->GetPdgCode();
2541         charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
2542         TLorentzVector momentum(part->Px(),part->Py(),part->Pz(),part->E());
2543         if(charge != 0)
2544         {
2545           if(part->Pt()> GetReader()->GetCTSPtMin())
2546           {
2547             //Particles in CTS acceptance
2548             Bool_t inCTS =  GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
2549             Int_t indexmother=part->GetMother();
2550             if(indexmother>-1)
2551             {
2552               Int_t mPdg = ((AliAODMCParticle*) mcparticles->At(indexmother)) ->GetPdgCode();
2553               if(TMath::Abs(pdg) == 11 && mPdg == 22) continue;
2554             }
2555
2556             if(inCTS)
2557             {
2558               if( label!=iParticle) // avoid trigger particle
2559               {
2560                 if(!FillChargedMCCorrelationHistograms(part->Pt(),part->Phi(),part->Eta(),ptprim,phiprim,etaprim)) return;
2561                 else lead = kTRUE;
2562               }
2563             } // in acceptance
2564           } // min pt cut
2565         } //only charged particles
2566       }  //MC particle loop
2567       if (lead) fhMCPtLeading->Fill(ptprim);
2568     } //when the leading particles could trace back to MC
2569   }// AOD MC
2570 }
2571
2572 //_____________________________________________________________________
2573 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
2574 {
2575
2576   //Print some relevant parameters set for the analysis
2577   if(! opt)
2578     return;
2579
2580   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2581   AliAnaCaloTrackCorrBaseClass::Print(" ");
2582   printf("Pt trigger           >  %3.2f\n", fMinTriggerPt) ;
2583   printf("Pt associated hadron <  %3.2f\n", fMaxAssocPt) ;
2584   printf("Pt associated hadron >  %3.2f\n", fMinAssocPt) ;
2585   printf("Phi trigger particle-Hadron   <  %3.2f\n", fDeltaPhiMaxCut) ;
2586   printf("Phi trigger particle-Hadron   >  %3.2f\n", fDeltaPhiMinCut) ;
2587   printf("Phi trigger particle-UeHadron <  %3.2f\n", fUeDeltaPhiMaxCut) ;
2588   printf("Phi trigger particle-UeHadron >  %3.2f\n", fUeDeltaPhiMinCut) ;
2589   printf("Isolated Trigger?  %d\n"     , fSelectIsolated) ;
2590   printf("Several UE?  %d\n"           , fMakeSeveralUE) ;
2591   printf("Name of AOD Pi0 Branch %s \n", fPi0AODBranchName.Data());
2592   printf("Do Decay-hadron correlation ?  %d\n", fPi0Trigger) ;
2593   printf("Select absolute leading for cluster triggers ?  %d or Near Side %d\n",
2594          fMakeAbsoluteLeading, fMakeNearSideLeading) ;
2595   printf("Trigger pt bins  %d\n", fNAssocPtBins) ;
2596   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
2597     printf("\t bin %d = [%2.1f,%2.1f]\n", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
2598   }
2599
2600 }
2601
2602 //____________________________________________________________
2603 void AliAnaParticleHadronCorrelation::SetNAssocPtBins(Int_t n)
2604 {
2605   // Set number of bins
2606
2607   fNAssocPtBins  = n ;
2608
2609
2610   if(n < 10 && n > 0)
2611   {
2612     fNAssocPtBins  = n ;
2613   }
2614   else
2615   {
2616     printf("n = larger than 9 or too small, set to 9 \n");
2617     fNAssocPtBins = 9;
2618   }
2619 }
2620
2621 //______________________________________________________________________________
2622 void AliAnaParticleHadronCorrelation::SetAssocPtBinLimit(Int_t ibin, Float_t pt)
2623 {
2624   // Set the list of limits for the trigger pt bins
2625
2626   if(ibin <= fNAssocPtBins || ibin >= 0)
2627   {
2628     fAssocPtBinLimit[ibin] = pt  ;
2629   }
2630   else
2631   {
2632     printf("AliAnaParticleHadronCorrelation::SetAssocPtBinLimit() - bin  number too large %d > %d or small, nothing done\n", ibin, fNAssocPtBins) ;
2633   }
2634 }
2635