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