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