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
  19 // in a previous analysis
  20 //
  21 //-- Author: Gustavo Conesa (LNF-INFN) (LPSC-IN2P3-CNRS)
  22 //           Yaxian Mao (LPSC-IN2P3-CNRS) and (CNWU) first usable implementation.
  23 //           Xiangrong Zhu (CNWU), implementtion of own mixing.
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 //#include "TClonesArray.h"
  29 #include <TClass.h>
  30 #include <TMath.h>
  31 #include <TH2F.h>
  32 #include <TDatabasePDG.h>
  33
  34 //---- ANALYSIS system ----
  35 #include "AliNeutralMesonSelection.h"
  36 #include "AliAnaParticleHadronCorrelation.h"
  37 #include "AliCaloTrackReader.h"
  38 #include "AliAODPWG4ParticleCorrelation.h"
  39 #include "AliFiducialCut.h"
  40 #include "AliVTrack.h"
  41 #include "AliVCluster.h"
  42 #include "AliMCAnalysisUtils.h"
  43 #include "TParticle.h"
  44 #include "AliStack.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliMixedEvent.h"
  47 #include "AliAnalysisManager.h"
  48 #include "AliInputEventHandler.h"
  49 #include "AliEventplane.h"
  50
  51 ClassImp(AliAnaParticleHadronCorrelation)
  52
  53
  54 //___________________________________________________________________
  55   AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation():
  56     AliAnaCaloTrackCorrBaseClass(),
  57     fFillAODWithReferences(0),      fCheckLeadingWithNeutralClusters(0),
  58     fMaxAssocPt(1000.),             fMinAssocPt(0.),
  59     fDeltaPhiMaxCut(0.),            fDeltaPhiMinCut(0.),
  60     fSelectIsolated(0),             fMakeSeveralUE(0),
  61     fUeDeltaPhiMaxCut(0.),          fUeDeltaPhiMinCut(0.),
  62     fPi0AODBranchName(""),          fNeutralCorr(0),
  63     fPi0Trigger(0),                 fDecayTrigger(0),
  64     fMakeAbsoluteLeading(0),        fMakeNearSideLeading(0),
  65     fLeadingTriggerIndex(-1),       fHMPIDCorrelation(0),  fFillBradHisto(0),
  66     fNAssocPtBins(0),               fAssocPtBinLimit(),
  67     fCorrelVzBin(0),
  68     fListMixTrackEvents(),          fListMixCaloEvents(),
  69     fUseMixStoredInReader(0),       fFillNeutralEventMixPool(0),
  70     fM02MaxCut(0),                  fM02MinCut(0),
  71     fFillPileUpHistograms(0),
  72     fSelectLeadingHadronAngle(0),
  73     fMinLeadHadPhi(0),              fMaxLeadHadPhi(0),
  74     fMinLeadHadPt(0),               fMaxLeadHadPt(0),
  75
  76     //Histograms
  77     fhPtTriggerInput(0),            fhPtTriggerSSCut(0),
  78     fhPtTriggerIsoCut(0),           fhPtTriggerFidCut(0),
  79     fhPtTrigger(0),                 fhPtTriggerVtxBC0(0),
  80     fhPtTriggerVzBin(0),            fhPtTriggerBin(0),
  81     fhPhiTrigger(0),                fhEtaTrigger(0),
  82     fhPtTriggerMC(),
  83     fhPtTriggerCentrality(0),       fhPtTriggerEventPlane(0),
  84     fhTriggerEventPlaneCentrality(0),
  85     fhPtTriggerMixed(0),            fhPtTriggerMixedVzBin(0), fhPtTriggerMixedBin(0),
  86     fhPhiTriggerMixed(0),           fhEtaTriggerMixed(0),
  87     fhPtLeadingOppositeHadron(0),   fhPtDiffPhiLeadingOppositeHadron(0), fhPtDiffEtaLeadingOppositeHadron(0),
  88     fhDeltaPhiDeltaEtaCharged(0),
  89     fhPhiCharged(0),                fhEtaCharged(0),
  90     fhDeltaPhiCharged(0),           fhDeltaEtaCharged(0),
  91     fhDeltaPhiChargedPt(0),         fhDeltaPhiUeChargedPt(0),
  92     fhUePart(0),
  93     fhXECharged(0),                 fhXECharged_Cone2(0),      fhXEUeCharged(0),
  94     fhXEPosCharged(0),              fhXENegCharged(0),
  95     fhPtHbpXECharged(0),            fhPtHbpXECharged_Cone2(0), fhPtHbpXEUeCharged(0),
  96     fhZTCharged(0),                 fhZTUeCharged(0),
  97     fhZTPosCharged(0),              fhZTNegCharged(0),
  98     fhPtHbpZTCharged(0),            fhPtHbpZTUeCharged(0),
  99     fhXEChargedMC(),                fhDeltaPhiChargedMC(),
 100     fhDeltaPhiDeltaEtaChargedPtA3GeV(0),
 101     fhDeltaPhiChargedPtA3GeV(0),    fhDeltaEtaChargedPtA3GeV(0),
 102     //Pile-Up
 103     fhDeltaPhiChargedPileUp(),      fhDeltaEtaChargedPileUp(),
 104     fhDeltaPhiChargedPtA3GeVPileUp(), fhDeltaEtaChargedPtA3GeVPileUp(),
 105     fhXEChargedPileUp(),            fhXEUeChargedPileUp(),
 106     fhZTChargedPileUp(),            fhZTUeChargedPileUp(),
 107     fhPtTrigChargedPileUp(),
 108     fhDeltaPhiChargedOtherBC(),     fhDeltaPhiChargedPtA3GeVOtherBC(),
 109     fhXEChargedOtherBC(),           fhXEUeChargedOtherBC(),
 110     fhZTChargedOtherBC(),           fhZTUeChargedOtherBC(),
 111     fhPtTrigChargedOtherBC(),
 112     fhDeltaPhiChargedBC0(),         fhDeltaPhiChargedPtA3GeVBC0(),
 113     fhXEChargedBC0(),               fhXEUeChargedBC0(),
 114     fhZTChargedBC0(),               fhZTUeChargedBC0(),
 115     fhPtTrigChargedBC0(),
 116     fhDeltaPhiChargedVtxBC0(),      fhDeltaPhiChargedPtA3GeVVtxBC0(),
 117     fhXEChargedVtxBC0(),            fhXEUeChargedVtxBC0(),
 118     fhZTChargedVtxBC0(),            fhZTUeChargedVtxBC0(),
 119     fhPtTrigChargedVtxBC0(),
 120     fhDeltaPhiUeLeftCharged(0),
 121     fhDeltaPhiUeLeftUpCharged(0),   fhDeltaPhiUeRightUpCharged(0),
 122     fhDeltaPhiUeLeftDownCharged(0), fhDeltaPhiUeRightDownCharged(0),
 123     fhXEUeLeftCharged(0),
 124     fhXEUeLeftUpCharged(0),         fhXEUeRightUpCharged(0),
 125     fhXEUeLeftDownCharged(0),       fhXEUeRightDownCharged(0),
 126     fhPtHbpXEUeLeftCharged(0),      fhZTUeLeftCharged(0),
 127     fhPtHbpZTUeLeftCharged(0),
 128     fhPtTrigPout(0),                fhPtTrigCharged(0),
 129     fhDeltaPhiChargedMult(0x0),     fhDeltaEtaChargedMult(0x0),
 130     fhXEMult(0x0),                  fhXEUeMult(0x0),
 131     fhZTMult(0x0),                  fhZTUeMult(0x0),
 132     fhAssocPtBkg(0),                fhDeltaPhiDeltaEtaAssocPtBin(0),
 133     fhDeltaPhiAssocPtBin(0),
 134     fhDeltaPhiAssocPtBinDEta08(0),  fhDeltaPhiAssocPtBinDEta0(0),
 135     fhDeltaPhiAssocPtBinHMPID(0),   fhDeltaPhiAssocPtBinHMPIDAcc(0),
 136     fhDeltaPhiBradAssocPtBin(0),    fhDeltaPhiBrad(0),
 137     fhXEAssocPtBin(0),              fhZTAssocPtBin(0),
 138     fhDeltaPhiDeltaEtaNeutral(0),
 139     fhPhiNeutral(0),                fhEtaNeutral(0),
 140     fhDeltaPhiNeutral(0),           fhDeltaEtaNeutral(0),
 141     fhDeltaPhiNeutralPt(0),         fhDeltaPhiUeNeutralPt(0),
 142     fhXENeutral(0),                 fhXEUeNeutral(0),
 143     fhPtHbpXENeutral(0),            fhPtHbpXEUeNeutral(0),
 144     fhZTNeutral(0),                 fhZTUeNeutral(0),
 145     fhPtHbpZTNeutral(0),            fhPtHbpZTUeNeutral(0),
 146     fhDeltaPhiUeLeftNeutral(0),     fhXEUeLeftNeutral(0),
 147     fhPtHbpXEUeLeftNeutral(0),      fhZTUeLeftNeutral(0),
 148     fhPtHbpZTUeLeftNeutral(0),      fhPtPi0DecayRatio(0),
 149     fhDeltaPhiDecayCharged(0),      fhXEDecayCharged(0),           fhZTDecayCharged(0),
 150     fhDeltaPhiDecayNeutral(0),      fhXEDecayNeutral(0),           fhZTDecayNeutral(0),
 151     fhDeltaPhiDecayChargedAssocPtBin(0),
 152     fh2phiTriggerParticle(0x0),
 153     fhMCPtTrigger(0),               fhMCPhiTrigger(0),             fhMCEtaTrigger(0),
 154     fhMCPtTriggerNotLeading(0),     fhMCPhiTriggerNotLeading(0),   fhMCEtaTriggerNotLeading(0),
 155     fhMCEtaCharged(0),              fhMCPhiCharged(0),
 156     fhMCDeltaEtaCharged(0),         fhMCDeltaPhiCharged(0x0),
 157     fhMCDeltaPhiDeltaEtaCharged(0), fhMCDeltaPhiChargedPt(0),
 158     fhMCPtXECharged(0),             fhMCPtXEUeCharged(0),
 159     fhMCPtXEUeLeftCharged(0),
 160     fhMCPtHbpXECharged(0),          fhMCPtHbpXEUeCharged(0),
 161     fhMCPtHbpXEUeLeftCharged(0),
 162     fhMCUePart(0),
 163     fhMCPtZTCharged(0),             fhMCPtZTUeCharged(0),
 164     fhMCPtZTUeLeftCharged(0),
 165     fhMCPtHbpZTCharged(0),          fhMCPtHbpZTUeCharged(0),
 166     fhMCPtHbpZTUeLeftCharged(0),
 167     fhMCPtTrigPout(0),              fhMCPtAssocDeltaPhi(0),
 168     //Mixing
 169     fhNEventsTrigger(0),            fhNtracksMB(0),                 fhNclustersMB(0),
 170     fhMixDeltaPhiCharged(0),        fhMixDeltaPhiDeltaEtaCharged(0),
 171     fhMixXECharged(0),              fhMixXEUeCharged(0),            fhMixHbpXECharged(0),
 172     fhMixDeltaPhiChargedAssocPtBin(),
 173     fhMixDeltaPhiChargedAssocPtBinDEta08(),
 174     fhMixDeltaPhiChargedAssocPtBinDEta0(),
 175     fhMixDeltaPhiDeltaEtaChargedAssocPtBin(),
 176     fhEventBin(0),                  fhEventMixBin(0)
 177 {
 178   //Default Ctor
 179
 180   //Initialize parameters
 181   InitParameters();
 182
 183   for(Int_t i = 0; i < 7; i++)
 184   {
 185     fhPtTriggerMC[i] = 0;
 186     fhXEChargedMC[i] = 0;
 187     fhDeltaPhiChargedMC[i] = 0;
 188   }
 189
 190   for(Int_t i = 0; i < 7; i++)
 191   {
 192     fhPtTriggerPileUp             [i] = 0 ;
 193     fhDeltaPhiChargedPileUp       [i] = 0 ; fhDeltaEtaChargedPileUp       [i] = 0 ;
 194     fhXEChargedPileUp             [i] = 0 ; fhXEUeChargedPileUp           [i] = 0 ;
 195     fhZTChargedPileUp             [i] = 0 ; fhZTUeChargedPileUp           [i] = 0 ;
 196     fhPtTrigChargedPileUp         [i] = 0 ;
 197     fhDeltaPhiChargedPtA3GeVPileUp[i] = 0 ; fhDeltaEtaChargedPtA3GeVPileUp[i] = 0 ;
 198   }
 199
 200 }
 201
 202 //_________________________________________________________________
 203 AliAnaParticleHadronCorrelation::~AliAnaParticleHadronCorrelation()
 204 {
 205   // Remove event containers
 206
 207   if(DoOwnMix())
 208   {
 209     if(fListMixTrackEvents)
 210     {
 211       for(Int_t iz=0; iz < GetNZvertBin(); iz++)
 212       {
 213         for(Int_t ic=0; ic < GetNCentrBin(); ic++)
 214         {
 215           for(Int_t irp=0; irp<GetNRPBin(); irp++)
 216           {
 217             Int_t bin = GetEventMixBin(ic, iz, irp);
 218             fListMixTrackEvents[bin]->Delete() ;
 219             delete fListMixTrackEvents[bin] ;
 220           }
 221         }
 222       }
 223     }
 224
 225     delete[] fListMixTrackEvents;
 226
 227     if(fListMixCaloEvents)
 228     {
 229       for(Int_t iz=0; iz < GetNZvertBin(); iz++)
 230       {
 231         for(Int_t ic=0; ic < GetNCentrBin(); ic++)
 232         {
 233           for(Int_t irp=0; irp<GetNRPBin(); irp++)
 234           {
 235             Int_t bin = GetEventMixBin(ic, iz, irp);
 236             fListMixCaloEvents[bin]->Delete() ;
 237             delete fListMixCaloEvents[bin] ;
 238           }
 239         }
 240       }
 241     }
 242
 243     delete[] fListMixCaloEvents;
 244
 245   }
 246 }
 247
 248 //__________________________________________________________________________________________________________________________________________
 249 void AliAnaParticleHadronCorrelation::FillChargedAngularCorrelationHistograms(Float_t ptAssoc,  Float_t ptTrig,      Int_t   bin,
 250                                                                               Float_t phiAssoc, Float_t phiTrig,     Float_t deltaPhi,
 251                                                                               Float_t etaAssoc, Float_t etaTrig,
 252                                                                               Bool_t  decay,    Float_t hmpidSignal, Int_t  outTOF,
 253                                                                               Int_t   cen,      Int_t   mcTag)
 254 {
 255   // Fill angular correlation related histograms
 256
 257   Float_t deltaEta    = etaTrig-etaAssoc;
 258   Float_t deltaPhiOrg = phiTrig-phiAssoc;
 259
 260   fhEtaCharged       ->Fill(ptAssoc,etaAssoc);
 261   fhPhiCharged       ->Fill(ptAssoc,phiAssoc);
 262   fhDeltaEtaCharged  ->Fill(ptTrig ,deltaEta);
 263   fhDeltaPhiCharged  ->Fill(ptTrig ,deltaPhi);
 264   fhDeltaPhiChargedPt->Fill(ptAssoc, deltaPhi);
 265   fhDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
 266
 267   if(ptAssoc > 3 )
 268   {
 269     fhDeltaEtaChargedPtA3GeV        ->Fill(ptTrig  ,deltaEta);
 270     fhDeltaPhiChargedPtA3GeV        ->Fill(ptTrig  ,deltaPhi);
 271     fhDeltaPhiDeltaEtaChargedPtA3GeV->Fill(deltaPhi, deltaEta);
 272   }
 273
 274   // Pile up studies
 275
 276   if(fFillPileUpHistograms)
 277   {
 278     if     (outTOF==1)
 279     {
 280       fhDeltaPhiChargedOtherBC->Fill(ptTrig ,deltaPhi) ;
 281       if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVOtherBC->Fill(ptTrig ,deltaPhi) ;
 282     }
 283     else if(outTOF==0)
 284     {
 285         fhDeltaPhiChargedBC0->Fill(ptTrig ,deltaPhi) ;
 286         if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVBC0->Fill(ptTrig ,deltaPhi) ;
 287     }
 288
 289     Int_t vtxBC = GetReader()->GetVertexBC();
 290     if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
 291     {
 292       fhDeltaPhiChargedVtxBC0->Fill(ptTrig ,deltaPhi) ;
 293       if(ptAssoc > 3 ) fhDeltaPhiChargedPtA3GeVVtxBC0->Fill(ptTrig ,deltaPhi) ;
 294     }
 295
 296     if(GetReader()->IsPileUpFromSPD())               { fhDeltaEtaChargedPileUp[0]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[0]->Fill(ptTrig ,deltaPhi) ; }
 297     if(GetReader()->IsPileUpFromEMCal())             { fhDeltaEtaChargedPileUp[1]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[1]->Fill(ptTrig ,deltaPhi) ; }
 298     if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhDeltaEtaChargedPileUp[2]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[2]->Fill(ptTrig ,deltaPhi) ; }
 299     if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhDeltaEtaChargedPileUp[3]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[3]->Fill(ptTrig ,deltaPhi) ; }
 300     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhDeltaEtaChargedPileUp[4]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[4]->Fill(ptTrig ,deltaPhi) ; }
 301     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhDeltaEtaChargedPileUp[5]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[5]->Fill(ptTrig ,deltaPhi) ; }
 302     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhDeltaEtaChargedPileUp[6]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPileUp[6]->Fill(ptTrig ,deltaPhi) ; }
 303
 304     if(ptAssoc > 3 )
 305     {
 306       if(GetReader()->IsPileUpFromSPD())               { fhDeltaEtaChargedPtA3GeVPileUp[0]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[0]->Fill(ptTrig ,deltaPhi) ; }
 307       if(GetReader()->IsPileUpFromEMCal())             { fhDeltaEtaChargedPtA3GeVPileUp[1]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[1]->Fill(ptTrig ,deltaPhi) ; }
 308       if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhDeltaEtaChargedPtA3GeVPileUp[2]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[2]->Fill(ptTrig ,deltaPhi) ; }
 309       if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhDeltaEtaChargedPtA3GeVPileUp[3]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[3]->Fill(ptTrig ,deltaPhi) ; }
 310       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhDeltaEtaChargedPtA3GeVPileUp[4]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[4]->Fill(ptTrig ,deltaPhi) ; }
 311       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhDeltaEtaChargedPtA3GeVPileUp[5]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[5]->Fill(ptTrig ,deltaPhi) ; }
 312       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhDeltaEtaChargedPtA3GeVPileUp[6]->Fill(ptTrig ,deltaEta) ; fhDeltaPhiChargedPtA3GeVPileUp[6]->Fill(ptTrig ,deltaPhi) ; }
 313     }
 314   }
 315
 316   if(IsDataMC())
 317   {
 318     Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
 319     fhDeltaPhiChargedMC[mcIndex]->Fill(ptTrig , deltaPhi);
 320   }
 321
 322   if(fDecayTrigger && decay) fhDeltaPhiDecayCharged   ->Fill(ptTrig  , deltaPhi);
 323
 324   Double_t  dphiBrad = -100;
 325   if(fFillBradHisto)
 326   {
 327     dphiBrad = atan2(sin(deltaPhiOrg), cos(deltaPhiOrg))/TMath::Pi();//-1 to 1
 328     if(TMath::Abs(dphiBrad)>0.325 && TMath::Abs(dphiBrad)<0.475)  //Hardcoded values, BAD, FIXME
 329     {
 330       fhAssocPtBkg->Fill(ptTrig, ptAssoc);
 331     }
 332
 333     if(dphiBrad<-1./3) dphiBrad += 2;
 334     fhDeltaPhiBrad->Fill(ptTrig, dphiBrad);
 335   }
 336
 337   // Fill histograms in bins of associated particle pT
 338   if(bin>=0)
 339   {
 340       fhDeltaPhiDeltaEtaAssocPtBin    [bin]->Fill(deltaPhi,deltaEta);
 341
 342       fhDeltaPhiAssocPtBin            [bin]->Fill(ptTrig, deltaPhi);
 343
 344     if(TMath::Abs(deltaEta)> 0.8)
 345       fhDeltaPhiAssocPtBinDEta08      [bin]->Fill(ptTrig, deltaPhi);
 346
 347     if(TMath::Abs(deltaEta)< 0.01)
 348       fhDeltaPhiAssocPtBinDEta0       [bin]->Fill(ptTrig, deltaPhi);
 349
 350     if (fFillBradHisto)
 351       fhDeltaPhiBradAssocPtBin        [bin]->Fill(ptTrig, dphiBrad);
 352
 353     if(fDecayTrigger && decay)
 354       fhDeltaPhiDecayChargedAssocPtBin[bin]->Fill(ptTrig, deltaPhi);
 355
 356     if(fHMPIDCorrelation)
 357     {
 358       if( hmpidSignal > 0 )
 359       {
 360         //printf("Track pt %f with HMPID signal %f \n",pt,hmpidSignal);
 361         fhDeltaPhiAssocPtBinHMPID[bin]->Fill(ptTrig, deltaPhi);
 362       }
 363
 364       if(phiAssoc > 5*TMath::DegToRad() && phiAssoc < 20*TMath::DegToRad())
 365       {
 366         //printf("Track pt %f in HMPID acceptance phi %f \n ",pt,phi*TMath::RadToDeg() );
 367         fhDeltaPhiAssocPtBinHMPIDAcc[bin]->Fill(ptTrig, deltaPhi);
 368       }
 369     }
 370   }
 371
 372   //fill different multiplicity/centrality histogram
 373   if(fFillHighMultHistograms)
 374   {
 375     fhDeltaPhiChargedMult[cen]->Fill(ptTrig,deltaPhi);
 376     fhDeltaEtaChargedMult[cen]->Fill(ptTrig,deltaEta);
 377   }
 378 }
 379
 380 //___________________________________________________________________________________________________________________________________
 381 Bool_t AliAnaParticleHadronCorrelation::FillChargedMCCorrelationHistograms(Float_t mcAssocPt, Float_t mcAssocPhi, Float_t mcAssocEta,
 382                                                                            Float_t mcTrigPt,  Float_t mcTrigPhi,  Float_t mcTrigEta)
 383 {
 384   // Fill MC histograms independently of AOD or ESD
 385
 386   Bool_t lead = kTRUE;
 387
 388   // In case we requested the trigger to be a leading particle,
 389   // check if this is true at the MC level.
 390   // Not sure if it is correct to skip or not skip this.
 391   // Absolute leading?
 392   if( fMakeAbsoluteLeading && mcAssocPt > mcTrigPt )     lead = kFALSE; // skip event
 393
 394   // Skip this event if near side associated particle pt larger than trigger
 395   if( fMakeNearSideLeading && mcAssocPt > mcTrigPt &&
 396      TMath::Abs(mcAssocPhi-mcTrigPhi)<TMath::PiOver2() ) lead = kFALSE; // skip event
 397
 398   //
 399   // Select only hadrons in pt range
 400   if ( mcAssocPt < fMinAssocPt || mcAssocPt > fMaxAssocPt ) return lead ; // exclude but continue
 401   if ( mcAssocPt < GetReader()->GetCTSPtMin())              return lead ;
 402
 403   if( mcAssocPhi < 0 ) mcAssocPhi+=TMath::TwoPi();
 404
 405   //
 406   // Remove trigger itself for correlation when use charged triggers
 407   if(TMath::Abs(mcAssocPt -mcTrigPt ) < 1e-6 &&
 408      TMath::Abs(mcAssocPhi-mcTrigPhi) < 1e-6 &&
 409      TMath::Abs(mcAssocEta-mcTrigEta) < 1e-6)            return lead ; // exclude but continue
 410
 411   Float_t mcdeltaPhi= mcTrigPhi-mcAssocPhi;
 412   if(mcdeltaPhi <= -TMath::PiOver2()) mcdeltaPhi+=TMath::TwoPi();
 413   if(mcdeltaPhi > 3*TMath::PiOver2()) mcdeltaPhi-=TMath::TwoPi();
 414
 415   Float_t mcxE    =-mcAssocPt/mcTrigPt*TMath::Cos(mcdeltaPhi);// -(mcAssocPx*pxprim+mcAssocPy*pyprim)/(mcTrigPt*mcTrigPt);
 416   Float_t mchbpXE =-100 ;
 417   if(mcxE > 0 ) mchbpXE = TMath::Log(1./mcxE);
 418
 419   Float_t mczT    = mcAssocPt/mcTrigPt ;
 420   Float_t mchbpZT =-100 ;
 421   if(mczT > 0 ) mchbpZT = TMath::Log(1./mczT);
 422
 423   Double_t mcpout = mcAssocPt*TMath::Sin(mcdeltaPhi) ;
 424
 425   if(GetDebug() > 0 )
 426   {
 427     AliInfo(Form("Charged hadron: track Pt %f, track Phi %f, phi trigger %f. Cuts:  delta phi  %2.2f < %2.2f < %2.2f \n",
 428                  mcAssocPt,mcAssocPhi, mcTrigPhi,fDeltaPhiMinCut, mcdeltaPhi, fDeltaPhiMaxCut));
 429   }
 430
 431   // Fill Histograms
 432   fhMCEtaCharged     ->Fill(mcAssocPt, mcAssocEta);
 433   fhMCPhiCharged     ->Fill(mcAssocPt, mcAssocPhi);
 434   fhMCDeltaEtaCharged->Fill(mcTrigPt , mcTrigEta-mcAssocEta);
 435   fhMCDeltaPhiCharged->Fill(mcTrigPt , mcdeltaPhi);
 436   fhMCPtAssocDeltaPhi->Fill(mcAssocPt, mcdeltaPhi);
 437
 438   fhMCDeltaPhiDeltaEtaCharged->Fill(mcdeltaPhi,mcTrigEta-mcAssocEta);
 439
 440   //delta phi cut for correlation
 441   if( (mcdeltaPhi > fDeltaPhiMinCut) && (mcdeltaPhi < fDeltaPhiMaxCut) )
 442   {
 443     fhMCDeltaPhiChargedPt->Fill(mcAssocPt,mcdeltaPhi);
 444     fhMCPtXECharged      ->Fill(mcTrigPt, mcxE);
 445     fhMCPtHbpXECharged   ->Fill(mcTrigPt, mchbpXE);
 446     fhMCPtZTCharged      ->Fill(mcTrigPt, mczT);
 447     fhMCPtHbpZTCharged   ->Fill(mcTrigPt, mchbpZT);
 448     fhMCPtTrigPout       ->Fill(mcTrigPt, mcpout) ;
 449   }
 450
 451   // Underlying event
 452
 453   // Right
 454   if ( (mcdeltaPhi > fUeDeltaPhiMinCut) && (mcdeltaPhi < fUeDeltaPhiMaxCut) )
 455   {
 456     //Double_t randomphi = gRandom->Uniform(TMath::Pi()/2,3*TMath::Pi()/2);
 457     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 458     Double_t mcUexE = -(mcAssocPt/mcTrigPt)*TMath::Cos(randomphi);
 459     Double_t mcUezT =   mcAssocPt/mcTrigPt;
 460
 461     if(mcUexE < 0.) mcUexE = -mcUexE;
 462
 463     fhMCPtXEUeCharged->Fill(mcTrigPt,mcUexE);
 464     if(mcUexE > 0) fhMCPtHbpXEUeCharged->Fill(mcTrigPt,TMath::Log(1/mcUexE));
 465
 466     fhMCPtZTUeCharged->Fill(mcTrigPt,mcUezT);
 467     if(mcUezT > 0) fhMCPtHbpZTUeCharged->Fill(mcTrigPt,TMath::Log(1/mcUezT));
 468
 469     fhMCUePart->Fill(mcTrigPt);
 470   }
 471
 472   if(fMakeSeveralUE)
 473   {
 474     // Left
 475     if((mcdeltaPhi<-fUeDeltaPhiMinCut) || (mcdeltaPhi >2*fUeDeltaPhiMaxCut))
 476     {
 477       Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 478       Double_t mcUexE = -(mcAssocPt/mcTrigPt)*TMath::Cos(randomphi);
 479       Double_t mcUezT =   mcAssocPt/mcTrigPt;
 480
 481       if(mcUexE < 0.) mcUexE = -mcUexE;
 482
 483       fhMCPtXEUeLeftCharged->Fill(mcTrigPt,mcUexE);
 484       if(mcUexE > 0) fhMCPtHbpXEUeLeftCharged->Fill(mcTrigPt,TMath::Log(1/mcUexE));
 485
 486       fhMCPtZTUeLeftCharged->Fill(mcTrigPt,mcUezT);
 487       if(mcUexE > 0) fhMCPtHbpZTUeLeftCharged->Fill(mcTrigPt,TMath::Log(1/mcUezT));
 488     }
 489   }
 490
 491   return lead;
 492 }
 493
 494 //___________________________________________________________________________________________________________________
 495 void AliAnaParticleHadronCorrelation::FillChargedMomentumImbalanceHistograms(Float_t ptTrig,   Float_t ptAssoc,
 496                                                                              Float_t deltaPhi,
 497                                                                              Int_t   cen,      Int_t   charge,
 498                                                                              Int_t   bin,      Bool_t  decay,
 499                                                                              Int_t   outTOF,   Int_t   mcTag)
 500
 501 {
 502   // Fill mostly momentum imbalance related histograms
 503
 504   Float_t zT = ptAssoc/ptTrig ;
 505   Float_t hbpZT = -100;
 506   if(zT > 0 ) hbpZT = TMath::Log(1./zT);
 507   else        hbpZT =-100;
 508
 509   Float_t xE =-ptAssoc/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 510   Float_t hbpXE = -100;
 511   //if(xE <0.)xE =-xE;
 512   if(xE > 0 ) hbpXE = TMath::Log(1./xE);
 513   else        hbpXE =-100;
 514
 515   Float_t pout = ptAssoc*TMath::Sin(deltaPhi) ;
 516
 517   fhXECharged         ->Fill(ptTrig , xE);
 518   fhPtHbpXECharged    ->Fill(ptTrig , hbpXE);
 519   fhZTCharged         ->Fill(ptTrig , zT);
 520   fhPtHbpZTCharged    ->Fill(ptTrig , hbpZT);
 521   fhPtTrigPout        ->Fill(ptTrig , pout) ;
 522   fhPtTrigCharged     ->Fill(ptTrig , ptAssoc) ;
 523   if((deltaPhi > 5*TMath::Pi()/6.)   && (deltaPhi < 7*TMath::Pi()/6.))
 524   {
 525     fhXECharged_Cone2         ->Fill(ptTrig , xE);
 526     fhPtHbpXECharged_Cone2    ->Fill(ptTrig , hbpXE);
 527   }
 528
 529   // Pile up studies
 530   if(fFillPileUpHistograms)
 531   {
 532     if     (outTOF==1)
 533     {
 534       fhXEChargedOtherBC    ->Fill(ptTrig,xE);
 535       fhZTChargedOtherBC    ->Fill(ptTrig,zT);
 536       fhPtTrigChargedOtherBC->Fill(ptTrig,ptAssoc);
 537     }
 538     else if(outTOF==0)
 539     {
 540       fhXEChargedBC0    ->Fill(ptTrig,xE);
 541       fhZTChargedBC0    ->Fill(ptTrig,zT);
 542       fhPtTrigChargedBC0->Fill(ptTrig,ptAssoc);
 543     }
 544
 545     Int_t vtxBC = GetReader()->GetVertexBC();
 546     if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
 547     {
 548       fhXEChargedVtxBC0    ->Fill(ptTrig,xE);
 549       fhZTChargedVtxBC0    ->Fill(ptTrig,zT);
 550       fhPtTrigChargedVtxBC0->Fill(ptTrig,ptAssoc);
 551     }
 552
 553     if(GetReader()->IsPileUpFromSPD())                { fhXEChargedPileUp[0]->Fill(ptTrig,xE); fhZTChargedPileUp[0]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[0]->Fill(ptTrig,ptAssoc); }
 554     if(GetReader()->IsPileUpFromEMCal())              { fhXEChargedPileUp[1]->Fill(ptTrig,xE); fhZTChargedPileUp[1]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[1]->Fill(ptTrig,ptAssoc); }
 555     if(GetReader()->IsPileUpFromSPDOrEMCal())         { fhXEChargedPileUp[2]->Fill(ptTrig,xE); fhZTChargedPileUp[2]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[2]->Fill(ptTrig,ptAssoc); }
 556     if(GetReader()->IsPileUpFromSPDAndEMCal())        { fhXEChargedPileUp[3]->Fill(ptTrig,xE); fhZTChargedPileUp[3]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[3]->Fill(ptTrig,ptAssoc); }
 557     if(GetReader()->IsPileUpFromSPDAndNotEMCal())     { fhXEChargedPileUp[4]->Fill(ptTrig,xE); fhZTChargedPileUp[4]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[4]->Fill(ptTrig,ptAssoc); }
 558     if(GetReader()->IsPileUpFromEMCalAndNotSPD())     { fhXEChargedPileUp[5]->Fill(ptTrig,xE); fhZTChargedPileUp[5]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[5]->Fill(ptTrig,ptAssoc); }
 559     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal())  { fhXEChargedPileUp[6]->Fill(ptTrig,xE); fhZTChargedPileUp[6]->Fill(ptTrig,zT); fhPtTrigChargedPileUp[6]->Fill(ptTrig,ptAssoc); }
 560   }
 561
 562   if(IsDataMC())
 563   {
 564     Int_t mcIndex = GetMCTagHistogramIndex(mcTag);
 565     fhXEChargedMC[mcIndex]->Fill(ptTrig , xE);
 566   }
 567
 568   if(fDecayTrigger && decay)
 569   {
 570     fhXEDecayCharged->Fill(ptTrig,xE);
 571     fhZTDecayCharged->Fill(ptTrig,zT);
 572   } // photon decay pi0/eta trigger
 573
 574   if(bin >= 0 )//away side
 575   {
 576     fhXEAssocPtBin[bin]->Fill(ptTrig, xE) ;
 577     fhZTAssocPtBin[bin]->Fill(ptTrig, zT) ;
 578   }
 579
 580   if(charge > 0)
 581   {
 582     fhXEPosCharged->Fill(ptTrig,xE) ;
 583     fhZTPosCharged->Fill(ptTrig,zT) ;
 584   }
 585   else
 586   {
 587     fhXENegCharged->Fill(ptTrig,xE) ;
 588     fhZTNegCharged->Fill(ptTrig,zT) ;
 589   }
 590
 591   //fill different multiplicity/centrality histogram
 592   if(fFillHighMultHistograms)
 593   {
 594     fhXEMult[cen]->Fill(ptTrig,xE);
 595     fhZTMult[cen]->Fill(ptTrig,zT);
 596   } //multiplicity/centrality events selection
 597 }
 598
 599 //_______________________________________________________________________________________________________________________
 600 void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventHistograms(Float_t ptTrig,   Float_t ptAssoc,
 601                                                                            Float_t deltaPhi, Int_t cen, Int_t outTOF)
 602 {
 603   // Fill underlying event histograms
 604
 605   fhUePart->Fill(ptTrig);
 606
 607   fhDeltaPhiUeChargedPt->Fill(ptAssoc,deltaPhi);
 608
 609   Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 610   Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 611   Double_t uezT =   ptAssoc/ptTrig;
 612
 613   if(uexE < 0.) uexE = -uexE;
 614
 615   fhXEUeCharged->Fill(ptTrig,uexE);
 616   if(uexE > 0) fhPtHbpXEUeCharged->Fill(ptTrig,TMath::Log(1/uexE));
 617
 618   fhZTUeCharged->Fill(ptTrig,uezT);
 619   if(uezT > 0) fhPtHbpZTUeCharged->Fill(ptTrig,TMath::Log(1/uezT));
 620
 621   // Pile up studies
 622
 623   if(fFillPileUpHistograms)
 624   {
 625     if     (outTOF==1)
 626     {
 627       fhXEUeChargedOtherBC->Fill(ptTrig,uexE);
 628       fhZTUeChargedOtherBC->Fill(ptTrig,uezT);
 629     }
 630     else if(outTOF==0)
 631     {
 632       fhXEUeChargedBC0->Fill(ptTrig,uexE);
 633       fhZTUeChargedBC0->Fill(ptTrig,uezT);
 634     }
 635
 636     Int_t vtxBC = GetReader()->GetVertexBC();
 637     if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)
 638     {
 639       fhXEUeChargedVtxBC0->Fill(ptTrig,uexE);
 640       fhZTUeChargedVtxBC0->Fill(ptTrig,uezT);
 641     }
 642
 643     if(GetReader()->IsPileUpFromSPD())               { fhXEUeChargedPileUp[0]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[0]->Fill(ptTrig,uezT);}
 644     if(GetReader()->IsPileUpFromEMCal())             { fhXEUeChargedPileUp[1]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[1]->Fill(ptTrig,uezT);}
 645     if(GetReader()->IsPileUpFromSPDOrEMCal())        { fhXEUeChargedPileUp[2]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[2]->Fill(ptTrig,uezT);}
 646     if(GetReader()->IsPileUpFromSPDAndEMCal())       { fhXEUeChargedPileUp[3]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[3]->Fill(ptTrig,uezT);}
 647     if(GetReader()->IsPileUpFromSPDAndNotEMCal())    { fhXEUeChargedPileUp[4]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[4]->Fill(ptTrig,uezT);}
 648     if(GetReader()->IsPileUpFromEMCalAndNotSPD())    { fhXEUeChargedPileUp[5]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[5]->Fill(ptTrig,uezT);}
 649     if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) { fhXEUeChargedPileUp[6]->Fill(ptTrig,uexE); fhZTUeChargedPileUp[6]->Fill(ptTrig,uezT);}
 650   }
 651
 652   //fill different multiplicity/centrality histogram
 653   if(fFillHighMultHistograms)
 654   {
 655     fhXEUeMult[cen]->Fill(ptTrig,uexE);
 656     fhZTUeMult[cen]->Fill(ptTrig,uezT);
 657   }
 658 }
 659
 660 //_____________________________________________________________________________________________________
 661 void AliAnaParticleHadronCorrelation::FillChargedUnderlyingEventSidesHistograms(Float_t ptTrig,
 662                                                                                 Float_t ptAssoc,
 663                                                                                 Float_t deltaPhi)
 664 {
 665   // Fill underlying event histograms to the left and right of trigger
 666   // Right cone is the default UE.
 667
 668   if((deltaPhi<-fUeDeltaPhiMinCut) || (deltaPhi >2*fUeDeltaPhiMaxCut))
 669   {
 670     fhDeltaPhiUeLeftCharged->Fill(ptAssoc,deltaPhi);
 671     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 672     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 673     Double_t uezT =   ptAssoc/ptTrig;
 674
 675     if(uexE < 0.) uexE = -uexE;
 676
 677     fhXEUeLeftCharged->Fill(ptTrig,uexE);
 678     if(uexE > 0) fhPtHbpXEUeLeftCharged->Fill(ptTrig,TMath::Log(1/uexE));
 679
 680     fhZTUeLeftCharged->Fill(ptTrig,uezT);
 681     if(uexE > 0) fhPtHbpZTUeLeftCharged->Fill(ptTrig,TMath::Log(1/uezT));
 682     fhDeltaPhiUeLeftCharged->Fill(ptAssoc, deltaPhi);
 683   }
 684
 685   if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-TMath::Pi()/2))
 686   {
 687     fhDeltaPhiUeLeftDownCharged->Fill(ptAssoc,deltaPhi);
 688     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 689     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 690
 691     if(uexE < 0.) uexE = -uexE;
 692
 693     fhXEUeLeftDownCharged->Fill(ptTrig,uexE);
 694   }
 695
 696   if((deltaPhi>2*fUeDeltaPhiMaxCut) && (deltaPhi <3*TMath::Pi()/2))
 697   {
 698     fhDeltaPhiUeLeftUpCharged->Fill(ptAssoc,deltaPhi);
 699     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 700     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 701
 702     if(uexE < 0.) uexE = -uexE;
 703
 704     fhXEUeLeftUpCharged->Fill(ptTrig,uexE);
 705   }
 706
 707   if((deltaPhi > TMath::Pi()/2) && (deltaPhi < fUeDeltaPhiMaxCut))
 708   {
 709     fhDeltaPhiUeRightUpCharged->Fill(ptAssoc,deltaPhi);
 710     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 711     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 712
 713     if(uexE < 0.) uexE = -uexE;
 714
 715     fhXEUeRightUpCharged->Fill(ptTrig,uexE);
 716   }
 717
 718   if((deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < TMath::Pi()/2))
 719   {
 720     fhDeltaPhiUeRightDownCharged->Fill(ptAssoc,deltaPhi);
 721     Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 722     Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
 723
 724     if(uexE < 0.) uexE = -uexE;
 725
 726     fhXEUeRightDownCharged->Fill(ptTrig,uexE);
 727   }
 728 }
 729
 730 //______________________________________________________________________________________________________________________________
 731 void AliAnaParticleHadronCorrelation::FillDecayPhotonCorrelationHistograms(Float_t ptAssoc,     Float_t phiAssoc,
 732                                                                            TLorentzVector mom1, TLorentzVector mom2,
 733                                                                            Bool_t bChargedOrNeutral)
 734 {
 735   // Do correlation with decay photons of triggered pi0 or eta
 736
 737   // Calculate the correlation parameters
 738   Float_t ptDecay1 = mom1.Pt();
 739   Float_t ptDecay2 = mom2.Pt();
 740
 741   Float_t zTDecay1 = -100, zTDecay2 = -100;
 742   if(ptDecay1 > 0) zTDecay1 = ptAssoc/ptDecay1 ;
 743   if(ptDecay2 > 0) zTDecay2 = ptAssoc/ptDecay2 ;
 744
 745   Float_t deltaPhiDecay1 = mom1.Phi()-phiAssoc;
 746   if(deltaPhiDecay1< -TMath::PiOver2()) deltaPhiDecay1+=TMath::TwoPi();
 747   if(deltaPhiDecay1>3*TMath::PiOver2()) deltaPhiDecay1-=TMath::TwoPi();
 748
 749   Float_t deltaPhiDecay2 = mom2.Phi()-phiAssoc;
 750   if(deltaPhiDecay2< -TMath::PiOver2()) deltaPhiDecay2+=TMath::TwoPi();
 751   if(deltaPhiDecay2>3*TMath::PiOver2()) deltaPhiDecay2-=TMath::TwoPi();
 752
 753   Float_t xEDecay1   =-zTDecay1*TMath::Cos(deltaPhiDecay1); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 754   Float_t xEDecay2   =-zTDecay2*TMath::Cos(deltaPhiDecay2); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 755
 756   if(bChargedOrNeutral) // correlate with charges
 757   {
 758     fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaPhiDecay1);
 759     fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaPhiDecay2);
 760
 761     if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::FillDecayPhotonHistograms( Charged corr) - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
 762
 763     if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
 764     {
 765       fhZTDecayCharged->Fill(ptDecay1,zTDecay1);
 766       fhXEDecayCharged->Fill(ptDecay1,xEDecay1);
 767     }
 768     if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
 769     {
 770       fhZTDecayCharged->Fill(ptDecay2,zTDecay2);
 771       fhXEDecayCharged->Fill(ptDecay2,xEDecay2);
 772     }
 773   }
 774   else // correlate with neutrals
 775   {
 776     fhDeltaPhiDecayNeutral->Fill(ptDecay1, deltaPhiDecay1);
 777     fhDeltaPhiDecayNeutral->Fill(ptDecay2, deltaPhiDecay2);
 778
 779     if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::FillDecayPhotonHistograms(Neutral corr) - deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaPhiDecay1, deltaPhiDecay2);
 780
 781     if( (deltaPhiDecay1 > fDeltaPhiMinCut) && ( deltaPhiDecay1 < fDeltaPhiMaxCut) )
 782     {
 783       fhZTDecayNeutral->Fill(ptDecay1,zTDecay1);
 784       fhXEDecayNeutral->Fill(ptDecay1,xEDecay1);
 785     }
 786     if( (deltaPhiDecay2 > fDeltaPhiMinCut) && ( deltaPhiDecay2 < fDeltaPhiMaxCut) )
 787     {
 788       fhZTDecayNeutral->Fill(ptDecay2,zTDecay2);
 789       fhXEDecayNeutral->Fill(ptDecay2,xEDecay2);
 790     }
 791   }
 792 }
 793
 794 //_____________________________________________________________________________________________________________________________
 795 void AliAnaParticleHadronCorrelation::FillNeutralUnderlyingEventSidesHistograms(Float_t ptTrig,   Float_t ptAssoc,
 796                                                                                 Float_t zT,       Float_t hbpZT,
 797                                                                                 Float_t deltaPhi)
 798 {
 799   // Fill underlying event histograms to the left of trigger
 800   // Right is the default case
 801
 802   Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
 803
 804   Float_t xE  =-ptAssoc/ptTrig*TMath::Cos(randomphi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
 805   Float_t hbpXE = -100;
 806   if(xE > 0 ) hbpXE = TMath::Log(1./xE);
 807
 808   if((deltaPhi<-fUeDeltaPhiMinCut) && (deltaPhi >-fUeDeltaPhiMaxCut))
 809   {
 810     fhDeltaPhiUeLeftNeutral->Fill(ptAssoc, deltaPhi);
 811     fhXEUeLeftNeutral      ->Fill(ptTrig , xE);
 812     fhPtHbpXEUeLeftNeutral ->Fill(ptTrig , hbpXE);
 813     fhZTUeLeftNeutral      ->Fill(ptTrig , zT);
 814     fhPtHbpZTUeLeftNeutral ->Fill(ptTrig , hbpZT);
 815   }
 816 }
 817
 818 //______________________________________________________
 819 void AliAnaParticleHadronCorrelation::FillEventMixPool()
 820 {
 821   // Fill the pool with tracks or clusters if requested
 822
 823   if ( !DoOwnMix() ) return;
 824
 825   FillChargedEventMixPool();
 826
 827   // Do the cluster pool filling only if requested
 828   // or in case of isolation cut using clusters in the cone.
 829   Bool_t isoCase = OnlyIsolated() && (GetIsolationCut()->GetParticleTypeInCone() != AliIsolationCut::kOnlyCharged);
 830
 831   if( !fFillNeutralEventMixPool && !isoCase) return;
 832
 833   FillNeutralEventMixPool();
 834 }
 835
 836 //_____________________________________________________________
 837 void AliAnaParticleHadronCorrelation::FillChargedEventMixPool()
 838 {
 839   // Mixed event pool filling for tracks
 840
 841   if(fUseMixStoredInReader && GetReader()->GetLastTracksMixedEvent() == GetEventNumber())
 842   {
 843     //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
 844     return ; // pool filled previously for another trigger
 845   }
 846
 847   Int_t nTracks = GetCTSTracks()->GetEntriesFast();
 848
 849   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
 850   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
 851
 852   if(!inputHandler) return ;
 853
 854   // Do mixing only with MB event (or the chosen mask), if not skip
 855   if( !(inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask()) ) return ;
 856
 857   fhNtracksMB->Fill(nTracks);
 858
 859   Int_t eventBin = GetEventMixBin();
 860
 861   //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
 862   if(eventBin < 0) return;
 863
 864   TObjArray * mixEventTracks = new TObjArray;
 865
 866   if(fUseMixStoredInReader)
 867   {
 868     fListMixTrackEvents[eventBin] = GetReader()->GetListWithMixedEventsForTracks(eventBin);
 869   }
 870
 871   if(!fListMixTrackEvents[eventBin]) fListMixTrackEvents[eventBin] = new TList();
 872
 873   //printf("%s ***** Pool Event bin : %d - nTracks %d\n",GetInputAODName().Data(),eventBin, GetCTSTracks()->GetEntriesFast());
 874
 875   TList * pool = fListMixTrackEvents[eventBin];
 876
 877   TVector3 p3;
 878   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
 879   {
 880     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
 881
 882     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 883     p3.SetXYZ(mom[0],mom[1],mom[2]);
 884     Float_t pt   = p3.Pt();
 885
 886     //Select only hadrons in pt range
 887     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
 888
 889     AliAODPWG4Particle * mixedTrack = new AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
 890     mixedTrack->SetDetector("CTS");
 891     mixedTrack->SetChargedBit(track->Charge()>0);
 892     mixEventTracks->Add(mixedTrack);
 893   }
 894
 895   //Set the event number where the last event was added, to avoid double pool filling
 896   GetReader()->SetLastTracksMixedEvent(GetEventNumber());
 897
 898   //printf("Add event to pool with %d tracks \n ",mixEventTracks->GetEntries());
 899   pool->AddFirst(mixEventTracks);
 900   mixEventTracks = 0;
 901
 902   //printf("Pool size %d, max %d\n",pool->GetSize(), GetNMaxEvMix());
 903
 904   if(pool->GetSize() > GetNMaxEvMix())
 905   {//Remove last event
 906     TClonesArray * tmp = static_cast<TClonesArray*>(pool->Last()) ;
 907     pool->RemoveLast() ;
 908     delete tmp ;
 909   }
 910 }
 911
 912 //_____________________________________________________________
 913 void AliAnaParticleHadronCorrelation::FillNeutralEventMixPool()
 914 {
 915   // Mixed event pool filling for neutral clusters
 916   // Right now only for EMCAL and in isolation case
 917
 918   //printf("FillNeutralEventMixPool for %s\n",GetInputAODName().Data());
 919
 920   if(fUseMixStoredInReader && GetReader()->GetLastCaloMixedEvent() == GetEventNumber())
 921   {
 922     //printf("%s : Pool already filled for this event !!!\n",GetInputAODName().Data());
 923     return ; // pool filled previously for another trigger
 924   }
 925
 926   TObjArray * pl = GetEMCALClusters();
 927   //if (GetAODObjArrayName.Contains("PHOS") )pl    = GetPHOSClusters();
 928   //else                                     pl    = GetEMCALClusters();
 929
 930   Int_t nClusters   = pl->GetEntriesFast();
 931
 932   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
 933   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
 934
 935   if(!inputHandler) return ;
 936
 937   // Do mixing only with MB event (or the chosen mask), if not skip
 938   if( !(inputHandler->IsEventSelected( ) & GetReader()->GetMixEventTriggerMask()) ) return ;
 939
 940   fhNclustersMB->Fill(nClusters);
 941
 942   Int_t eventBin = GetEventMixBin();
 943
 944   //Check that the bin exists, if not (bad determination of RP, centrality or vz bin) do nothing
 945   if(eventBin < 0) return;
 946
 947   TObjArray * mixEventCalo = new TObjArray;
 948
 949   if(fUseMixStoredInReader)
 950   {
 951     fListMixCaloEvents[eventBin] = GetReader()->GetListWithMixedEventsForCalo(eventBin);
 952   }
 953
 954   if(!fListMixCaloEvents[eventBin]) fListMixCaloEvents[eventBin] = new TList();
 955
 956   TList * poolCalo = fListMixCaloEvents[eventBin];
 957
 958   TLorentzVector mom;
 959
 960   for(Int_t ipr = 0;ipr <  nClusters ; ipr ++ )
 961   {
 962     AliVCluster * calo = (AliVCluster *) (pl->At(ipr)) ;
 963
 964     // remove matched clusters
 965     if( IsTrackMatched( calo, GetReader()->GetInputEvent() ) ) continue ;
 966
 967     //Cluster momentum calculation
 968     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 969     {
 970       calo->GetMomentum(mom,GetVertex(0)) ;
 971     }//Assume that come from vertex in straight line
 972     else
 973     {
 974       Double_t vertex[]={0,0,0};
 975       calo->GetMomentum(mom,vertex) ;
 976     }
 977
 978     Float_t pt = mom.Pt();
 979     //Select only clusters in pt range
 980     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
 981
 982     AliAODPWG4Particle * mixedCalo = new AliAODPWG4Particle(mom);
 983     mixedCalo->SetDetector("EMCAL");
 984     mixEventCalo->Add(mixedCalo);
 985   }
 986
 987   //Set the event number where the last event was added, to avoid double pool filling
 988   GetReader()->SetLastCaloMixedEvent(GetEventNumber());
 989
 990   //printf("Add event to pool with %d clusters \n ",mixEventCalo->GetEntries());
 991   poolCalo->AddFirst(mixEventCalo);
 992   mixEventCalo = 0;
 993
 994   //printf("Pool size %d, max %d\n",poolCalo->GetSize(), GetNMaxEvMix());
 995
 996   if(poolCalo->GetSize() > GetNMaxEvMix())
 997   {//Remove last event
 998     TClonesArray * tmp = static_cast<TClonesArray*>(poolCalo->Last()) ;
 999     poolCalo->RemoveLast() ;
1000     delete tmp ;
1001   }
1002 }
1003
1004 //_________________________________________________________________________________________________________________
1005 Bool_t AliAnaParticleHadronCorrelation::FindLeadingOppositeHadronInWindow(AliAODPWG4ParticleCorrelation * particle)
1006 {
1007   // Select events where the leading charged particle in the opposite hemisphere
1008   // to the trigger particle is in a window centered at 180 from the trigger
1009
1010   Float_t phiTrig    = particle->Phi();
1011   Float_t etaTrig    = particle->Eta();
1012   Float_t ptTrig     = particle->Pt();
1013   Float_t ptLeadHad  = -100 ;
1014   Float_t phiLeadHad = -100 ;
1015   Float_t etaLeadHad = -100 ;
1016   TVector3 p3;
1017
1018   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
1019   {
1020     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
1021
1022     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
1023     p3.SetXYZ(mom[0],mom[1],mom[2]);
1024
1025     Float_t pt   = p3.Pt();
1026     Float_t phi  = p3.Phi() ;
1027     if(phi < 0 ) phi+= TMath::TwoPi();
1028
1029     if(pt > ptLeadHad && TMath::Abs(phi-phiTrig) > TMath::PiOver2()) // in opposite hemisphere
1030     {
1031       ptLeadHad  = pt ;
1032       phiLeadHad = phi;
1033       etaLeadHad = p3.Eta();
1034     }
1035   }// track loop
1036
1037   fhPtLeadingOppositeHadron       ->Fill(ptTrig, ptLeadHad);
1038   fhPtDiffPhiLeadingOppositeHadron->Fill(ptTrig,phiLeadHad-phiTrig);
1039   fhPtDiffEtaLeadingOppositeHadron->Fill(ptTrig,etaLeadHad-etaTrig);
1040
1041   if(GetDebug() > 1 )
1042   {
1043     printf("AliAnaParticleHadronCorrelation::FindLeadingOppositeHadronInWindow() pT %2.2f, phi %2.2f, eta %2.2f\n",
1044            ptLeadHad,phiLeadHad*TMath::RadToDeg(),etaLeadHad);
1045
1046     printf("\t  pT trig %2.2f, Dphi (trigger-hadron) %2.2f, Deta (trigger-hadron) %2.2f\n",
1047            ptTrig, (phiLeadHad-phiTrig)*TMath::RadToDeg(), etaLeadHad-etaTrig);
1048     printf("\t cuts pT: min %2.2f, max %2.2f; DPhi: min %2.2f, max %2.2f\n",fMinLeadHadPt,fMaxLeadHadPt,fMinLeadHadPhi*TMath::RadToDeg(),fMaxLeadHadPhi*TMath::RadToDeg());
1049   }
1050
1051   if( ptLeadHad < fMinLeadHadPt ||
1052       ptLeadHad > fMaxLeadHadPt ) return kFALSE;
1053
1054   //printf("Accept leading hadron pT \n");
1055
1056   if( TMath::Abs(phiLeadHad-phiTrig) < fMinLeadHadPhi ||
1057       TMath::Abs(phiLeadHad-phiTrig) > fMaxLeadHadPhi ) return kFALSE;
1058
1059   //printf("Accept leading hadron phi \n");
1060
1061   return kTRUE ;
1062 }
1063
1064 //____________________________________________________________
1065 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()
1066 {
1067   //Save parameters used for analysis
1068   TString parList ; //this will be list of parameters used for this analysis.
1069   const Int_t buffersize = 560;
1070   char onePar[buffersize] ;
1071
1072   snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---:") ;
1073   parList+=onePar ;
1074   snprintf(onePar,buffersize," %3.2f < Pt associated < %3.2f; ", fMinAssocPt,   fMaxAssocPt) ;
1075   parList+=onePar ;
1076   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-Hadron < %3.2f; ",    fDeltaPhiMinCut,   fDeltaPhiMaxCut) ;
1077   parList+=onePar ;
1078   snprintf(onePar,buffersize," %3.2f < Phi trigger particle-UeHadron <  %3.2f; ", fUeDeltaPhiMinCut, fUeDeltaPhiMaxCut) ;
1079   parList+=onePar ;
1080   snprintf(onePar,buffersize,"Isolated Trigger?  %d;",    fSelectIsolated) ;
1081   parList+=onePar ;
1082   snprintf(onePar,buffersize,"Several UE?  %d;",          fMakeSeveralUE) ;
1083   parList+=onePar ;
1084   snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s;", fPi0AODBranchName.Data());
1085   parList+=onePar ;
1086   snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  pi0 %d, decay %d;", fPi0Trigger, fDecayTrigger) ;
1087   parList+=onePar ;
1088   snprintf(onePar,buffersize,"Select absolute leading for cluster triggers ? %d or Near Side Leading %d;",
1089            fMakeAbsoluteLeading, fMakeNearSideLeading) ;
1090   parList+=onePar ;
1091   snprintf(onePar,buffersize,"Associated particle pt bins  %d: ", fNAssocPtBins) ;
1092   parList+=onePar ;
1093   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
1094     snprintf(onePar,buffersize,"bin %d = [%2.1f,%2.1f];", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
1095   }
1096   parList+=onePar ;
1097
1098   //Get parameters set in base class.
1099   parList += GetBaseParametersList() ;
1100
1101   //Get parameters set in FiducialCut class (not available yet)
1102   //parlist += GetFidCut()->GetFidCutParametersList()
1103
1104   return new TObjString(parList) ;
1105
1106 }
1107
1108 //________________________________________________________________
1109 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()
1110 {
1111   // Create histograms to be saved in output file and
1112   // store them in fOutputContainer
1113
1114   TList * outputContainer = new TList() ;
1115   outputContainer->SetName("CorrelationHistos") ;
1116
1117   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();
1118   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();
1119   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();
1120
1121   Int_t nMixBins = GetNCentrBin()*GetNZvertBin()*GetNRPBin();
1122
1123   TString nameMC[]     = {"Photon","Pi0","Pi0Decay","EtaDecay","OtherDecay","Electron","Hadron"};
1124   TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1125
1126   // For vz dependent histograms, if option ON
1127   Int_t   nz  = 1  ;
1128   if(fCorrelVzBin) nz = GetNZvertBin();
1129   TString sz  = "" ;
1130   TString tz  = "" ;
1131
1132   fhPtTriggerInput  = new TH1F("hPtTriggerInput","Input trigger #it{p}_{T}", nptbins,ptmin,ptmax);
1133   fhPtTriggerInput->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1134   outputContainer->Add(fhPtTriggerInput);
1135
1136   if( fM02MaxCut > 0 && fM02MinCut > 0 )
1137   {
1138     fhPtTriggerSSCut  = new TH1F("hPtTriggerSSCut","Trigger #it{p}_{T} after #lambda^{2}_{0} cut", nptbins,ptmin,ptmax);
1139     fhPtTriggerSSCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1140     outputContainer->Add(fhPtTriggerSSCut);
1141   }
1142
1143   if( OnlyIsolated() )
1144   {
1145     fhPtTriggerIsoCut  = new TH1F("hPtTriggerIsoCut","Trigger #it{p}_{T} after isolation (and #lambda^{2}_{0}) cut", nptbins,ptmin,ptmax);
1146     fhPtTriggerIsoCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1147     outputContainer->Add(fhPtTriggerIsoCut);
1148   }
1149
1150   fhPtTriggerFidCut  = new TH1F("hPtTriggerFidCut","Trigger #it{p}_{T} after fiducial (isolation and #lambda^{2}_{0}) cut", nptbins,ptmin,ptmax);
1151   fhPtTriggerFidCut->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1152   outputContainer->Add(fhPtTriggerFidCut);
1153
1154   fhPtTrigger  = new TH1F("hPtTrigger","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax);
1155   fhPtTrigger->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1156   outputContainer->Add(fhPtTrigger);
1157
1158   if(IsDataMC())
1159   {
1160     for(Int_t i=0; i < 7; i++)
1161     {
1162       fhPtTriggerMC[i]  = new TH1F(Form("hPtTrigger_MC%s",nameMC[i].Data()),
1163                                    Form("#it{p}_{T} distribution of trigger particles, trigger origin is %s",nameMC[i].Data()),
1164                                    nptbins,ptmin,ptmax);
1165       fhPtTriggerMC[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1166       outputContainer->Add(fhPtTriggerMC[i]);
1167     }
1168   }
1169
1170   if(fCorrelVzBin)
1171   {
1172     fhPtTriggerVzBin  = new TH2F("hPtTriggerVzBin","#it{p}_{T} distribution of trigger particles vs vz bin", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
1173     fhPtTriggerVzBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1174     fhPtTriggerVzBin->SetYTitle("#it{v}_{#it{z}} bin");
1175     outputContainer->Add(fhPtTriggerVzBin);
1176   }
1177
1178   fhPtTriggerBin  = new TH2F ("hPtTriggerBin","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
1179   fhPtTriggerBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1180   fhPtTriggerBin->SetYTitle("Bin");
1181   outputContainer->Add(fhPtTriggerBin);
1182
1183   fhPhiTrigger  = new TH2F ("hPhiTrigger","#phi distribution of trigger Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1184   fhPhiTrigger->SetYTitle("#phi (rad)");
1185   outputContainer->Add(fhPhiTrigger);
1186
1187   fhEtaTrigger  = new TH2F ("hEtaTrigger","#eta distribution of trigger",nptbins,ptmin,ptmax, netabins,etamin,etamax);
1188   fhEtaTrigger->SetYTitle("#eta ");
1189   outputContainer->Add(fhEtaTrigger);
1190
1191   if(fFillHighMultHistograms)
1192   {
1193     fhPtTriggerCentrality   = new TH2F("hPtTriggerCentrality","Trigger particle #it{p}_{T} vs centrality",nptbins,ptmin,ptmax,100,0.,100) ;
1194     fhPtTriggerCentrality->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1195     fhPtTriggerCentrality->SetYTitle("Centrality (%)");
1196     outputContainer->Add(fhPtTriggerCentrality) ;
1197
1198     fhPtTriggerEventPlane  = new TH2F("hPtTriggerEventPlane","Trigger particle #it{p}_{T} vs event plane angle",nptbins,ptmin,ptmax, 100,0.,TMath::Pi()) ;
1199     fhPtTriggerEventPlane->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1200     fhPtTriggerEventPlane->SetXTitle("EP angle (rad)");
1201     outputContainer->Add(fhPtTriggerEventPlane) ;
1202
1203     fhTriggerEventPlaneCentrality  = new TH2F("hTriggerEventPlaneCentrality","Trigger particle centrality vs event plane angle",100,0.,100,100,0.,TMath::Pi()) ;
1204     fhTriggerEventPlaneCentrality->SetXTitle("Centrality (%)");
1205     fhTriggerEventPlaneCentrality->SetYTitle("EP angle (rad)");
1206     outputContainer->Add(fhTriggerEventPlaneCentrality) ;
1207   }
1208
1209   // Leading hadron in oposite side
1210   if(fSelectLeadingHadronAngle)
1211   {
1212     fhPtLeadingOppositeHadron  = new TH2F("hPtTriggerPtLeadingOppositeHadron","Leading hadron opposite to trigger vs trigger #it{p}_{T}",
1213                                           nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1214     fhPtLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1215     fhPtLeadingOppositeHadron->SetYTitle("#it{p}_{T}^{lead hadron} (GeV/#it{c})");
1216     outputContainer->Add(fhPtLeadingOppositeHadron);
1217
1218     fhPtDiffPhiLeadingOppositeHadron  = new TH2F("hPtTriggerDiffPhiTriggerLeadingOppositeHadron","#phi_{trigger}-#phi_{leading opposite hadron} vs #it{p}_{T}^{trig}",
1219                                                  nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1220     fhPtDiffPhiLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1221     fhPtDiffPhiLeadingOppositeHadron->SetYTitle("#phi_{trigger}-#phi_{leading opposite hadron} (rad)");
1222     outputContainer->Add(fhPtDiffPhiLeadingOppositeHadron);
1223
1224     fhPtDiffEtaLeadingOppositeHadron  = new TH2F("hPtTriggerDiffEtaTriggerPhiLeadingOppositeHadron","#eta_{trigger}-#eta_{leading opposite hadron} vs #it{p}_{T}^{trig}",
1225                                                  nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1226     fhPtDiffEtaLeadingOppositeHadron->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1227     fhPtDiffEtaLeadingOppositeHadron->SetYTitle("#eta_{trigger}-#eta_{leading opposite hadron}");
1228     outputContainer->Add(fhPtDiffEtaLeadingOppositeHadron);
1229   }
1230
1231   //Correlation with charged hadrons
1232
1233   fhDeltaPhiDeltaEtaCharged  = new TH2F
1234   ("hDeltaPhiDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}}",
1235    ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1236   fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
1237   fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
1238
1239   fhDeltaPhiDeltaEtaChargedPtA3GeV  = new TH2F
1240   ("hDeltaPhiDeltaEtaChargedPtA3GeV","#eta_{trigger} - #eta_{h^{#pm}} vs #phi_{trigger} - #phi_{h^{#pm}, #it{p}_{TA}>3 GeV/#it{c}}",
1241    ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1242   fhDeltaPhiDeltaEtaChargedPtA3GeV->SetXTitle("#Delta #phi (rad)");
1243   fhDeltaPhiDeltaEtaChargedPtA3GeV->SetYTitle("#Delta #eta");
1244
1245   fhPhiCharged  = new TH2F
1246   ("hPhiCharged","#phi_{h^{#pm}}  vs #it{p}_{T #pm}",
1247    nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1248   fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");
1249   fhPhiCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
1250
1251   fhEtaCharged  = new TH2F
1252   ("hEtaCharged","#eta_{h^{#pm}}  vs #it{p}_{T #pm}",
1253    nptbins,ptmin,ptmax,netabins,etamin,etamax);
1254   fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
1255   fhEtaCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
1256
1257   fhDeltaPhiCharged  = new TH2F
1258   ("hDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
1259    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1260   fhDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
1261   fhDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1262
1263   fhDeltaPhiChargedPtA3GeV  = new TH2F
1264   ("hDeltaPhiChargedPtA3GeV","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}",
1265    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1266   fhDeltaPhiChargedPtA3GeV->SetYTitle("#Delta #phi (rad)");
1267   fhDeltaPhiChargedPtA3GeV->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1268
1269
1270   fhDeltaPhiChargedPt  = new TH2F
1271   ("hDeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs #it{p}_{T h^{#pm}}",
1272    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1273   fhDeltaPhiChargedPt->SetYTitle("#Delta #phi (rad)");
1274   fhDeltaPhiChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1275
1276   fhDeltaEtaCharged  = new TH2F
1277   ("hDeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}",
1278    nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1279   fhDeltaEtaCharged->SetYTitle("#Delta #eta");
1280   fhDeltaEtaCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1281
1282   fhDeltaEtaChargedPtA3GeV  = new TH2F
1283   ("hDeltaEtaChargedPtA3GeV","#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}",
1284    nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1285   fhDeltaEtaChargedPtA3GeV->SetYTitle("#Delta #eta");
1286   fhDeltaEtaChargedPtA3GeV->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1287
1288   fhXECharged  =
1289   new TH2F("hXECharged","#it{x}_{#it{E}} for charged tracks",
1290            nptbins,ptmin,ptmax,200,0.,2.);
1291   fhXECharged->SetYTitle("#it{x}_{#it{E}}");
1292   fhXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1293
1294   fhXECharged_Cone2  =
1295   new TH2F("hXECharged_Cone2","#it{x}_{#it{E}} for charged tracks in cone 2 (5#pi/6-7#pi/6)",
1296            nptbins,ptmin,ptmax,200,0.,2.);
1297   fhXECharged_Cone2->SetYTitle("#it{x}_{#it{E}}");
1298   fhXECharged_Cone2->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1299
1300   fhXEPosCharged  =
1301   new TH2F("hXEPositiveCharged","#it{x}_{#it{E}} for positive charged tracks",
1302            nptbins,ptmin,ptmax,200,0.,2.);
1303   fhXEPosCharged->SetYTitle("#it{x}_{#it{E}}");
1304   fhXEPosCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1305
1306   fhXENegCharged  =
1307   new TH2F("hXENegativeCharged","#it{x}_{#it{E}} for negative charged tracks",
1308            nptbins,ptmin,ptmax,200,0.,2.);
1309   fhXENegCharged->SetYTitle("#it{x}_{#it{E}}");
1310   fhXENegCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1311
1312   fhPtHbpXECharged  =
1313   new TH2F("hHbpXECharged","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
1314            nptbins,ptmin,ptmax,200,0.,10.);
1315   fhPtHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1316   fhPtHbpXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1317
1318   fhPtHbpXECharged_Cone2  =
1319   new TH2F("hHbpXECharged_Cone2","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons in cone 2 (5#pi/6-7#pi/6)",
1320            nptbins,ptmin,ptmax,200,0.,10.);
1321   fhPtHbpXECharged_Cone2->SetYTitle("ln(1/#it{x}_{#it{E}})");
1322   fhPtHbpXECharged_Cone2->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1323
1324   fhZTCharged  =
1325   new TH2F("hZTCharged","#it{z}_{T} for charged tracks",
1326            nptbins,ptmin,ptmax,200,0.,2.);
1327   fhZTCharged->SetYTitle("#it{z}_{T}");
1328   fhZTCharged->SetXTitle("#it{p}_{T trigger}");
1329
1330   fhZTPosCharged  =
1331   new TH2F("hZTPositiveCharged","#it{z}_{T} for positive charged tracks",
1332            nptbins,ptmin,ptmax,200,0.,2.);
1333   fhZTPosCharged->SetYTitle("#it{z}_{T}");
1334   fhZTPosCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1335
1336   fhZTNegCharged  =
1337   new TH2F("hZTNegativeCharged","#it{z}_{T} for negative charged tracks",
1338            nptbins,ptmin,ptmax,200,0.,2.);
1339   fhZTNegCharged->SetYTitle("#it{z}_{T}");
1340   fhZTNegCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1341
1342   fhPtHbpZTCharged  =
1343   new TH2F("hHbpZTCharged","#xi = ln(1/#it{z}_{T}) with charged hadrons",
1344            nptbins,ptmin,ptmax,200,0.,10.);
1345   fhPtHbpZTCharged->SetYTitle("ln(1/#it{z}_{T})");
1346   fhPtHbpZTCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1347
1348   fhPtTrigPout  =
1349   new TH2F("hPtTrigPout","Pout with triggers",
1350            nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
1351   fhPtTrigPout->SetYTitle("#it{p}_{out} (GeV/#it{c})");
1352   fhPtTrigPout->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1353
1354   fhPtTrigCharged  =
1355   new TH2F("hPtTrigCharged","trigger and charged tracks pt distribution",
1356            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1357   fhPtTrigCharged->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1358   fhPtTrigCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1359
1360   outputContainer->Add(fhDeltaPhiDeltaEtaCharged);
1361   outputContainer->Add(fhDeltaPhiDeltaEtaChargedPtA3GeV);
1362   outputContainer->Add(fhPhiCharged) ;
1363   outputContainer->Add(fhEtaCharged) ;
1364   outputContainer->Add(fhDeltaPhiCharged) ;
1365   outputContainer->Add(fhDeltaPhiChargedPtA3GeV) ;
1366   outputContainer->Add(fhDeltaEtaCharged) ;
1367   outputContainer->Add(fhDeltaEtaChargedPtA3GeV) ;
1368   outputContainer->Add(fhDeltaPhiChargedPt) ;
1369
1370   outputContainer->Add(fhXECharged) ;
1371   outputContainer->Add(fhXECharged_Cone2) ;
1372
1373   if(IsDataMC())
1374   {
1375     for(Int_t i=0; i < 7; i++)
1376     {
1377
1378       fhDeltaPhiChargedMC[i]  = new TH2F(Form("hDeltaPhiCharged_MC%s",nameMC[i].Data()),
1379                                          Form("#Delta #phi for charged tracks, trigger origin is %s",nameMC[i].Data()),
1380                                          nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
1381       fhDeltaPhiChargedMC[i]->SetYTitle("#it{x}_{#it{E}}");
1382       fhDeltaPhiChargedMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1383       outputContainer->Add(fhDeltaPhiChargedMC[i]) ;
1384
1385       fhXEChargedMC[i]  = new TH2F(Form("hXECharged_MC%s",nameMC[i].Data()),
1386                                    Form("#it{x}_{#it{E}} for charged tracks, trigger origin is %s",nameMC[i].Data()),
1387                                    nptbins,ptmin,ptmax,200,0.,2.);
1388       fhXEChargedMC[i]->SetYTitle("#it{x}_{#it{E}}");
1389       fhXEChargedMC[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1390       outputContainer->Add(fhXEChargedMC[i]) ;
1391     }
1392   }
1393
1394   outputContainer->Add(fhXEPosCharged) ;
1395   outputContainer->Add(fhXENegCharged) ;
1396   outputContainer->Add(fhPtHbpXECharged) ;
1397   outputContainer->Add(fhPtHbpXECharged_Cone2) ;
1398
1399   outputContainer->Add(fhZTCharged) ;
1400   outputContainer->Add(fhZTPosCharged) ;
1401   outputContainer->Add(fhZTNegCharged) ;
1402   outputContainer->Add(fhPtHbpZTCharged) ;
1403
1404   outputContainer->Add(fhPtTrigPout) ;
1405   outputContainer->Add(fhPtTrigCharged) ;
1406
1407   TString right = "";
1408   if(fMakeSeveralUE) right = "Right";
1409
1410   fhUePart  =  new TH1F("hUePart","UE particles distribution vs pt trig",
1411                         nptbins,ptmin,ptmax);
1412   fhUePart->SetYTitle("dNch");
1413   fhUePart->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1414
1415   fhDeltaPhiUeChargedPt  = new TH2F
1416   (Form("hDeltaPhiUe%sChargedPt",right.Data()),"#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}}",
1417    nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1418   fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi (rad)");
1419   fhDeltaPhiUeChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1420
1421   fhXEUeCharged  =
1422   new TH2F(Form("hXEUeCharged%s",right.Data()),"#it{x}_{#it{E}} for Underlying Event",
1423            nptbins,ptmin,ptmax,200,0.,2.);
1424   fhXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
1425   fhXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1426
1427   fhPtHbpXEUeCharged  =
1428   new TH2F(Form("hHbpXEUeCharged%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for Underlying Event",
1429            nptbins,ptmin,ptmax,200,0.,10.);
1430   fhPtHbpXEUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1431   fhPtHbpXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1432
1433   fhZTUeCharged  =
1434   new TH2F(Form("hZTUeCharged%s",right.Data()),"#it{z}_{T} for Underlying Event",
1435            nptbins,ptmin,ptmax,200,0.,2.);
1436   fhZTUeCharged->SetYTitle("#it{z}_{T}");
1437   fhZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1438
1439   fhPtHbpZTUeCharged  =
1440   new TH2F(Form("hHbpZTUeCharged%s",right.Data()),"#xi = ln(1/#it{z}_{T}) for Underlying Event",
1441            nptbins,ptmin,ptmax,200,0.,10.);
1442   fhPtHbpZTUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1443   fhPtHbpZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1444
1445   outputContainer->Add(fhUePart);
1446   outputContainer->Add(fhDeltaPhiUeChargedPt) ;
1447   outputContainer->Add(fhXEUeCharged) ;
1448   outputContainer->Add(fhPtHbpXEUeCharged) ;
1449   outputContainer->Add(fhZTUeCharged) ;
1450   outputContainer->Add(fhPtHbpZTUeCharged) ;
1451
1452   if(fMakeSeveralUE)
1453   {
1454     fhDeltaPhiUeLeftCharged  = new TH2F
1455     ("hDeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left side range of trigger particles",
1456      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1457     fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi (rad)");
1458     fhDeltaPhiUeLeftCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1459     outputContainer->Add(fhDeltaPhiUeLeftCharged) ;
1460
1461     fhDeltaPhiUeLeftUpCharged  = new TH2F
1462     ("hDeltaPhiUeLeftUpChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left Up side range of trigger particles",
1463      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1464     fhDeltaPhiUeLeftUpCharged->SetYTitle("#Delta #phi (rad)");
1465     fhDeltaPhiUeLeftUpCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1466     outputContainer->Add(fhDeltaPhiUeLeftUpCharged) ;
1467
1468     fhDeltaPhiUeRightUpCharged  = new TH2F
1469     ("hDeltaPhiUeRightUpChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE right Up side range of trigger particles",
1470      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1471     fhDeltaPhiUeRightUpCharged->SetYTitle("#Delta #phi (rad)");
1472     fhDeltaPhiUeRightUpCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1473     outputContainer->Add(fhDeltaPhiUeRightUpCharged) ;
1474
1475     fhDeltaPhiUeLeftDownCharged  = new TH2F
1476     ("hDeltaPhiUeLeftDownChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE left Down side range of trigger particles",
1477      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1478     fhDeltaPhiUeLeftDownCharged->SetYTitle("#Delta #phi (rad)");
1479     fhDeltaPhiUeLeftDownCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1480     outputContainer->Add(fhDeltaPhiUeLeftDownCharged) ;
1481
1482     fhDeltaPhiUeRightDownCharged  = new TH2F
1483     ("hDeltaPhiUeRightDownChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs #it{p}_{T Ueh^{#pm}} with UE right Down side range of trigger particles",
1484      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1485     fhDeltaPhiUeRightDownCharged->SetYTitle("#Delta #phi (rad)");
1486     fhDeltaPhiUeRightDownCharged->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1487     outputContainer->Add(fhDeltaPhiUeRightDownCharged) ;
1488
1489     fhXEUeLeftCharged  =
1490     new TH2F("hXEUeChargedLeft","#it{x}_{#it{E}} with UE left side of trigger",
1491              nptbins,ptmin,ptmax,200,0.,2.);
1492     fhXEUeLeftCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
1493     fhXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1494     outputContainer->Add(fhXEUeLeftCharged) ;
1495
1496     fhXEUeLeftUpCharged  =
1497     new TH2F("hXEUeChargedLeftUp","#it{x}_{#it{E}} with UE left Up side of trigger",
1498              nptbins,ptmin,ptmax,200,0.,2.);
1499     fhXEUeLeftUpCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
1500     fhXEUeLeftUpCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1501     outputContainer->Add(fhXEUeLeftUpCharged) ;
1502
1503     fhXEUeRightUpCharged  =
1504     new TH2F("hXEUeChargedRightUp","#it{x}_{#it{E} h^{#pm}} with UE right Up side of trigger",
1505              nptbins,ptmin,ptmax,200,0.,2.);
1506     fhXEUeRightUpCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
1507     fhXEUeRightUpCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1508     outputContainer->Add(fhXEUeRightUpCharged) ;
1509
1510     fhXEUeLeftDownCharged  =
1511     new TH2F("hXEUeChargedLeftDown","#it{x}_{#it{E}} with UE left Down side of trigger",
1512              nptbins,ptmin,ptmax,200,0.,2.);
1513     fhXEUeLeftDownCharged->SetYTitle("#it{x}_{#it{E} Ueh^{#pm}}");
1514     fhXEUeLeftDownCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1515     outputContainer->Add(fhXEUeLeftDownCharged) ;
1516
1517     fhXEUeRightDownCharged  =
1518     new TH2F("hXEUeChargedRightDown","#it{x}_{#it{E} h^{#pm}} with UE right Down side of trigger",
1519              nptbins,ptmin,ptmax,200,0.,2.);
1520     fhXEUeRightDownCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
1521     fhXEUeRightDownCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1522     outputContainer->Add(fhXEUeRightDownCharged) ;
1523
1524     fhPtHbpXEUeLeftCharged  =
1525     new TH2F("hHbpXEUeChargedLeft","#xi = ln(1/#it{x}_{#it{E}}) with charged UE left side of trigger",
1526              nptbins,ptmin,ptmax,200,0.,10.);
1527     fhPtHbpXEUeLeftCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
1528     fhPtHbpXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1529     outputContainer->Add(fhPtHbpXEUeLeftCharged) ;
1530
1531     fhZTUeLeftCharged  =
1532     new TH2F("hZTUeChargedLeft","#it{z}_{trigger h^{#pm}} = #it{p}_{T Ueh^{#pm}} / #it{p}_{T trigger} with UE left side of trigger",
1533              nptbins,ptmin,ptmax,200,0.,2.);
1534     fhZTUeLeftCharged->SetYTitle("#it{z}_{trigger Ueh^{#pm}}");
1535     fhZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1536     outputContainer->Add(fhZTUeLeftCharged) ;
1537
1538     fhPtHbpZTUeLeftCharged  =
1539     new TH2F("hHbpZTUeChargedLeft","#xi = ln(1/#it{z}_{T}) with charged UE left side of trigger",
1540              nptbins,ptmin,ptmax,200,0.,10.);
1541     fhPtHbpZTUeLeftCharged->SetYTitle("ln(1/#it{z}_{T})");
1542     fhPtHbpZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1543     outputContainer->Add(fhPtHbpZTUeLeftCharged) ;
1544   }
1545
1546   if(fFillPileUpHistograms)
1547   {
1548     fhDeltaPhiChargedOtherBC  = new TH2F
1549     ("hDeltaPhiChargedOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC!=0",
1550      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1551     fhDeltaPhiChargedOtherBC->SetYTitle("#Delta #phi (rad)");
1552     fhDeltaPhiChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1553
1554     fhDeltaPhiChargedPtA3GeVOtherBC  = new TH2F
1555     ("hDeltaPhiChargedPtA3GeVOtherBC","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC!=0",
1556      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1557     fhDeltaPhiChargedPtA3GeVOtherBC->SetYTitle("#Delta #phi (rad)");
1558     fhDeltaPhiChargedPtA3GeVOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1559
1560     fhPtTrigChargedOtherBC  =
1561     new TH2F("hPtTrigChargedOtherBC","trigger and charged tracks pt distribution, track BC!=0",
1562              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1563     fhPtTrigChargedOtherBC->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1564     fhPtTrigChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1565
1566     fhXEChargedOtherBC  =
1567     new TH2F("hXEChargedOtherBC","#it{x}_{#it{E}} for charged tracks, track BC!=0",
1568              nptbins,ptmin,ptmax,200,0.,2.);
1569     fhXEChargedOtherBC->SetYTitle("#it{x}_{#it{E}}");
1570     fhXEChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1571
1572     fhXEUeChargedOtherBC  =
1573     new TH2F("hXEUeChargedOtherBC","#it{x}_{#it{E}} for Underlying Event, track BC!=0",
1574              nptbins,ptmin,ptmax,200,0.,2.);
1575     fhXEUeChargedOtherBC->SetYTitle("#it{x}_{#it{E}}");
1576     fhXEUeChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1577
1578     fhZTChargedOtherBC  =
1579     new TH2F("hZTChargedOtherBC","#it{z}_{T} for charged tracks, track BC!=0",
1580              nptbins,ptmin,ptmax,200,0.,2.);
1581     fhZTChargedOtherBC->SetYTitle("#it{z}_{T}");
1582     fhZTChargedOtherBC->SetXTitle("#it{p}_{T trigger}");
1583
1584     fhZTUeChargedOtherBC  =
1585     new TH2F("hZTUeChargedOtherBC","#it{z}_{T} for Underlying Event, track BC!=0",
1586              nptbins,ptmin,ptmax,200,0.,2.);
1587     fhZTUeChargedOtherBC->SetYTitle("#it{z}_{T}");
1588     fhZTUeChargedOtherBC->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1589
1590     outputContainer->Add(fhDeltaPhiChargedOtherBC) ;
1591     outputContainer->Add(fhDeltaPhiChargedPtA3GeVOtherBC) ;
1592     outputContainer->Add(fhXEChargedOtherBC) ;
1593     outputContainer->Add(fhXEUeChargedOtherBC) ;
1594     outputContainer->Add(fhZTChargedOtherBC) ;
1595     outputContainer->Add(fhZTUeChargedOtherBC) ;
1596     outputContainer->Add(fhPtTrigChargedOtherBC) ;
1597
1598     fhDeltaPhiChargedBC0  = new TH2F
1599     ("hDeltaPhiChargedBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC==0",
1600      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1601     fhDeltaPhiChargedBC0->SetYTitle("#Delta #phi (rad)");
1602     fhDeltaPhiChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1603
1604     fhDeltaPhiChargedPtA3GeVBC0  = new TH2F
1605     ("hDeltaPhiChargedPtA3GeVBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC==0",
1606      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1607     fhDeltaPhiChargedPtA3GeVBC0->SetYTitle("#Delta #phi (rad)");
1608     fhDeltaPhiChargedPtA3GeVBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1609
1610     fhPtTrigChargedBC0  =
1611     new TH2F("hPtTrigChargedBC0","trigger and charged tracks pt distribution, track BC==0",
1612              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1613     fhPtTrigChargedBC0->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1614     fhPtTrigChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1615
1616     fhXEChargedBC0  =
1617     new TH2F("hXEChargedBC0","#it{x}_{#it{E}} for charged tracks, track BC==0",
1618              nptbins,ptmin,ptmax,200,0.,2.);
1619     fhXEChargedBC0->SetYTitle("#it{x}_{#it{E}}");
1620     fhXEChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1621
1622     fhXEUeChargedBC0  =
1623     new TH2F("hXEUeChargedBC0","#it{x}_{#it{E}} for Underlying Event, track BC==0",
1624              nptbins,ptmin,ptmax,200,0.,2.);
1625     fhXEUeChargedBC0->SetYTitle("#it{x}_{#it{E}}");
1626     fhXEUeChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1627
1628     fhZTChargedBC0  =
1629     new TH2F("hZTChargedBC0","#it{z}_{T} for charged tracks, track BC==0",
1630              nptbins,ptmin,ptmax,200,0.,2.);
1631     fhZTChargedBC0->SetYTitle("#it{z}_{T}");
1632     fhZTChargedBC0->SetXTitle("#it{p}_{T trigger}");
1633
1634     fhZTUeChargedBC0  =
1635     new TH2F("hZTUeChargedBC0","#it{z}_{T} for Underlying Event, track BC==0",
1636              nptbins,ptmin,ptmax,200,0.,2.);
1637     fhZTUeChargedBC0->SetYTitle("#it{z}_{T}");
1638     fhZTUeChargedBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1639
1640     outputContainer->Add(fhDeltaPhiChargedBC0) ;
1641     outputContainer->Add(fhDeltaPhiChargedPtA3GeVBC0) ;
1642     outputContainer->Add(fhXEChargedBC0) ;
1643     outputContainer->Add(fhXEUeChargedBC0) ;
1644     outputContainer->Add(fhZTChargedBC0) ;
1645     outputContainer->Add(fhZTUeChargedBC0) ;
1646     outputContainer->Add(fhPtTrigChargedBC0) ;
1647
1648     fhPtTriggerVtxBC0  = new TH1F("hPtTriggerVtxBC0","#it{p}_{T} distribution of trigger particles", nptbins,ptmin,ptmax);
1649     fhPtTriggerVtxBC0->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1650
1651     fhDeltaPhiChargedVtxBC0  = new TH2F
1652     ("hDeltaPhiChargedVtxBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, track BC==0",
1653      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1654     fhDeltaPhiChargedVtxBC0->SetYTitle("#Delta #phi (rad)");
1655     fhDeltaPhiChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1656
1657     fhDeltaPhiChargedPtA3GeVVtxBC0  = new TH2F
1658     ("hDeltaPhiChargedPtA3GeVVtxBC0","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, track BC==0",
1659      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1660     fhDeltaPhiChargedPtA3GeVVtxBC0->SetYTitle("#Delta #phi (rad)");
1661     fhDeltaPhiChargedPtA3GeVVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1662
1663     fhPtTrigChargedVtxBC0  =
1664     new TH2F("hPtTrigChargedVtxBC0","trigger and charged tracks pt distribution, track BC==0",
1665              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1666     fhPtTrigChargedVtxBC0->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1667     fhPtTrigChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1668
1669     fhXEChargedVtxBC0  =
1670     new TH2F("hXEChargedVtxBC0","#it{x}_{#it{E}} for charged tracks, track BC==0",
1671              nptbins,ptmin,ptmax,200,0.,2.);
1672     fhXEChargedVtxBC0->SetYTitle("#it{x}_{#it{E}}");
1673     fhXEChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1674
1675     fhXEUeChargedVtxBC0  =
1676     new TH2F("hXEUeChargedVtxBC0","#it{x}_{#it{E}} for Underlying Event, track BC==0",
1677              nptbins,ptmin,ptmax,200,0.,2.);
1678     fhXEUeChargedVtxBC0->SetYTitle("#it{x}_{#it{E}}");
1679     fhXEUeChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1680
1681     fhZTChargedVtxBC0  =
1682     new TH2F("hZTChargedVtxBC0","#it{z}_{T} for charged tracks, track BC==0",
1683              nptbins,ptmin,ptmax,200,0.,2.);
1684     fhZTChargedVtxBC0->SetYTitle("#it{z}_{T}");
1685     fhZTChargedVtxBC0->SetXTitle("#it{p}_{T trigger}");
1686
1687     fhZTUeChargedVtxBC0  =
1688     new TH2F("hZTUeChargedVtxBC0","#it{z}_{T} for Underlying Event, track BC==0",
1689              nptbins,ptmin,ptmax,200,0.,2.);
1690     fhZTUeChargedVtxBC0->SetYTitle("#it{z}_{T}");
1691     fhZTUeChargedVtxBC0->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1692
1693     outputContainer->Add(fhPtTriggerVtxBC0);
1694     outputContainer->Add(fhDeltaPhiChargedVtxBC0) ;
1695     outputContainer->Add(fhDeltaPhiChargedPtA3GeVVtxBC0) ;
1696     outputContainer->Add(fhXEChargedVtxBC0) ;
1697     outputContainer->Add(fhXEUeChargedVtxBC0) ;
1698     outputContainer->Add(fhZTChargedVtxBC0) ;
1699     outputContainer->Add(fhZTUeChargedVtxBC0) ;
1700     outputContainer->Add(fhPtTrigChargedVtxBC0) ;
1701
1702     for(Int_t i = 0 ; i < 7 ; i++)
1703     {
1704       fhPtTriggerPileUp[i]  = new TH1F(Form("hPtTriggerPileUp%s",pileUpName[i].Data()),
1705                                        Form("#it{p}_{T} distribution of trigger particles, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1706       fhPtTriggerPileUp[i]->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
1707       outputContainer->Add(fhPtTriggerPileUp[i]);
1708
1709       fhDeltaPhiChargedPileUp[i]  = new TH2F(Form("hDeltaPhiChargedPileUp%s",pileUpName[i].Data()),
1710                                              Form("#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, %s Pile-Up event",pileUpName[i].Data()),
1711                                              nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1712       fhDeltaPhiChargedPileUp[i]->SetYTitle("#Delta #phi (rad)");
1713       fhDeltaPhiChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1714       outputContainer->Add(fhDeltaPhiChargedPileUp[i]) ;
1715
1716       fhDeltaPhiChargedPtA3GeVPileUp[i]  = new TH2F(Form("hDeltaPhiChargedPtA3GeVPileUp%s",pileUpName[i].Data()),
1717                                                     Form("#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, %s Pile-Up event",pileUpName[i].Data()),
1718                                                     nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1719       fhDeltaPhiChargedPtA3GeVPileUp[i]->SetYTitle("#Delta #phi (rad)");
1720       fhDeltaPhiChargedPtA3GeVPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1721       outputContainer->Add(fhDeltaPhiChargedPtA3GeVPileUp[i]) ;
1722
1723       fhDeltaEtaChargedPileUp[i]  = new TH2F(Form("hDeltaEtaChargedPileUp%s",pileUpName[i].Data()),
1724                                              Form("#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, %s Pile-Up event",pileUpName[i].Data()),
1725                                              nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1726       fhDeltaEtaChargedPileUp[i]->SetYTitle("#Delta #eta");
1727       fhDeltaEtaChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1728       outputContainer->Add(fhDeltaEtaChargedPileUp[i]) ;
1729
1730       fhDeltaEtaChargedPtA3GeVPileUp[i]  = new TH2F(Form("hDeltaEtaChargedPtA3GeVPileUp%s",pileUpName[i].Data()),
1731                                                     Form("#eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger}, #it{p}_{TA}>3 GeV/#it{c}, %s Pile-Up event",pileUpName[i].Data()),
1732                                                     nptbins,ptmin,ptmax,ndeltaetabins,deltaetamin,deltaetamax);
1733       fhDeltaEtaChargedPtA3GeVPileUp[i]->SetYTitle("#Delta #eta");
1734       fhDeltaEtaChargedPtA3GeVPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1735       outputContainer->Add(fhDeltaEtaChargedPtA3GeVPileUp[i]) ;
1736
1737       fhXEChargedPileUp[i]  = new TH2F(Form("hXEChargedPileUp%s",pileUpName[i].Data()),
1738                                        Form("#it{x}_{#it{E}} for charged tracks, %s Pile-Up event",pileUpName[i].Data()),
1739                                        nptbins,ptmin,ptmax,200,0.,2.);
1740       fhXEChargedPileUp[i]->SetYTitle("#it{x}_{#it{E}}");
1741       fhXEChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1742       outputContainer->Add(fhXEChargedPileUp[i]) ;
1743
1744       fhXEUeChargedPileUp[i]  = new TH2F(Form("hXEUeChargedPileUp%s",pileUpName[i].Data()),
1745                                          Form("#it{x}_{#it{E}} for Underlying Event, %s Pile-Up event",pileUpName[i].Data()),
1746                                          nptbins,ptmin,ptmax,200,0.,2.);
1747       fhXEUeChargedPileUp[i]->SetYTitle("#it{x}_{#it{E}}");
1748       fhXEUeChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1749       outputContainer->Add(fhXEUeChargedPileUp[i]) ;
1750
1751       fhZTChargedPileUp[i]  = new TH2F(Form("hZTChargedPileUp%s",pileUpName[i].Data()),
1752                                        Form("#it{z}_{T} for charged tracks, %s Pile-Up event",pileUpName[i].Data()),
1753                                        nptbins,ptmin,ptmax,200,0.,2.);
1754       fhZTChargedPileUp[i]->SetYTitle("#it{z}_{T}");
1755       fhZTChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1756       outputContainer->Add(fhZTChargedPileUp[i]) ;
1757
1758       fhZTUeChargedPileUp[i]  = new TH2F(Form("hZTUeChargedPileUp%s",pileUpName[i].Data()),
1759                                          Form("#it{z}_{T} for Underlying Event, %s Pile-Up event",pileUpName[i].Data()),
1760                                          nptbins,ptmin,ptmax,200,0.,2.);
1761       fhZTUeChargedPileUp[i]->SetYTitle("#it{z}_{T}");
1762       fhZTUeChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1763       outputContainer->Add(fhZTUeChargedPileUp[i]) ;
1764
1765       fhPtTrigChargedPileUp[i]  = new TH2F(Form("hPtTrigChargedPileUp%s",pileUpName[i].Data()),
1766                                            Form("trigger and charged tracks pt distribution, %s Pile-Up event",pileUpName[i].Data()),
1767                                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1768       fhPtTrigChargedPileUp[i]->SetYTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
1769       fhPtTrigChargedPileUp[i]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1770       outputContainer->Add(fhPtTrigChargedPileUp[i]) ;
1771
1772     }
1773   }
1774
1775   if(fFillHighMultHistograms)
1776   {
1777     Int_t nMultiBins = GetNCentrBin();
1778     fhDeltaPhiChargedMult = new TH2F*[nMultiBins] ;
1779     fhDeltaEtaChargedMult = new TH2F*[nMultiBins] ;
1780     fhXEMult              = new TH2F*[nMultiBins] ;
1781     fhXEUeMult            = new TH2F*[nMultiBins] ;
1782     fhZTMult              = new TH2F*[nMultiBins] ;
1783     fhZTUeMult            = new TH2F*[nMultiBins] ;
1784
1785     for(Int_t im=0; im<nMultiBins; im++)
1786     {
1787       fhDeltaPhiChargedMult[im]  = new TH2F
1788       (Form("hDeltaPhiCharged_Mult%d",im),Form("#Delta #phi charged Mult bin %d",im), nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1789       fhDeltaPhiChargedMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1790       fhDeltaPhiChargedMult[im]->SetYTitle("#Delta #phi (rad)");
1791
1792       fhDeltaEtaChargedMult[im]  = new TH2F
1793       (Form("hDeltaEtaCharged_Mult%d",im),Form("#Delta #eta charged Mult bin %d",im), nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
1794       fhDeltaEtaChargedMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1795       fhDeltaEtaChargedMult[im]->SetYTitle("#Delta #eta");
1796
1797       fhXEMult[im]  = new TH2F
1798       (Form("hXECharged_Mult%d",im),Form("#it{x}_{E} charged Mult bin %d",im), nptbins,ptmin,ptmax,200,0.,2.);
1799       fhXEMult[im]->SetYTitle("#it{x}_{E}");
1800       fhXEMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1801
1802       fhXEUeMult[im]  = new TH2F
1803       (Form("hXEUeCharged_Mult%d",im),Form("#it{x}_{E} UE charged Mult bin %d",im), nptbins,ptmin,ptmax,200,0.,2.);
1804       fhXEUeMult[im]->SetYTitle("#it{x}_{E}");
1805       fhXEUeMult[im]->SetXTitle("#it{p}_{T trigger}(GeV/#it{c})");
1806
1807       fhZTMult[im]  = new TH2F
1808       (Form("hZTCharged_Mult%d",im),Form("#it{z}_{T} charged  Mult bin %d",im), nptbins,ptmin,ptmax,200,0.,2.);
1809       fhZTMult[im]->SetYTitle("#it{z}_{T}");
1810       fhZTMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1811
1812       fhZTUeMult[im]  = new TH2F
1813       (Form("hZTUeCharged_Mult%d",im),Form("#it{z}_{T} UE charged  Mult bin %d",im), nptbins,ptmin,ptmax,200,0.,2.);
1814       fhZTUeMult[im]->SetYTitle("#it{z}_{T}");
1815       fhZTUeMult[im]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1816
1817       outputContainer->Add(fhDeltaPhiChargedMult[im]) ;
1818       outputContainer->Add(fhDeltaEtaChargedMult[im]) ;
1819       outputContainer->Add(fhXEMult  [im]);
1820       outputContainer->Add(fhXEUeMult[im]);
1821       outputContainer->Add(fhZTMult  [im]);
1822       outputContainer->Add(fhZTUeMult[im]);
1823     }
1824   }
1825
1826   if(fFillBradHisto)
1827   {
1828     fhAssocPtBkg        = new TH2F("hAssocPtBkg", " Trigger #it{p}_{T} vs associated hadron #it{p}_{T} from background",
1829                                    nptbins, ptmin, ptmax,nptbins,ptmin,ptmax);
1830     fhAssocPtBkg->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1831     fhAssocPtBkg->SetYTitle("#it{p}_{T associated} (GeV/#it{c})");
1832     outputContainer->Add(fhAssocPtBkg) ;
1833
1834     fhDeltaPhiBrad = new TH2F("hDeltaPhiBrad","atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs #it{p}_{T trigger} ",
1835                               nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
1836     fhDeltaPhiBrad->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1837     fhDeltaPhiBrad->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
1838     outputContainer->Add(fhDeltaPhiBrad) ;
1839   }
1840
1841   fhDeltaPhiDeltaEtaAssocPtBin = new TH2F*[fNAssocPtBins*nz];
1842   fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins*nz];
1843   fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
1844   fhDeltaPhiAssocPtBinDEta0  = new TH2F*[fNAssocPtBins*nz];
1845   fhXEAssocPtBin             = new TH2F*[fNAssocPtBins*nz];
1846   fhZTAssocPtBin             = new TH2F*[fNAssocPtBins*nz];
1847
1848   if(fFillBradHisto)
1849     fhDeltaPhiBradAssocPtBin = new TH2F*[fNAssocPtBins*nz];
1850
1851   if(fPi0Trigger || fDecayTrigger)
1852   {
1853     fhDeltaPhiAssocPtBin       = new TH2F*[fNAssocPtBins*nz];
1854     fhDeltaPhiAssocPtBinDEta08 = new TH2F*[fNAssocPtBins*nz];
1855     fhDeltaPhiDecayChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
1856   }
1857
1858   if(fHMPIDCorrelation)
1859   {
1860     fhDeltaPhiAssocPtBinHMPID   = new TH2F*[fNAssocPtBins*nz];
1861     fhDeltaPhiAssocPtBinHMPIDAcc= new TH2F*[fNAssocPtBins*nz];
1862   }
1863
1864   for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
1865   {
1866     for(Int_t z = 0 ; z < nz ; z++)
1867     {
1868       Int_t bin = i*nz+z;
1869
1870       if(fCorrelVzBin)
1871       {
1872         sz = Form("_vz%d",z);
1873         tz = Form(", #it{v}_{#it{z}} bin %d",z);
1874       }
1875
1876       //printf("iAssoc %d, Vz %d, bin %d - sz %s, tz %s \n",i,z,bin,sz.Data(),tz.Data());
1877
1878       fhDeltaPhiDeltaEtaAssocPtBin[bin]  = new TH2F(Form("hDeltaPhiDeltaEtaPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
1879                                                     Form("#Delta #phi vs #Delta #eta for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1880                                                     ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins,deltaetamin,deltaetamax);
1881       fhDeltaPhiDeltaEtaAssocPtBin[bin]->SetXTitle("#Delta #phi (rad)");
1882       fhDeltaPhiDeltaEtaAssocPtBin[bin]->SetYTitle("#Delta #eta");
1883
1884       fhDeltaPhiAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
1885                                            Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1886                                            nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1887       fhDeltaPhiAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1888       fhDeltaPhiAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
1889
1890       fhDeltaPhiAssocPtBinDEta08[bin] = new TH2F(Form("hDeltaPhiDeltaEta0.8PtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
1891                                                  Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta > 0.8", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1892                                                  nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1893       fhDeltaPhiAssocPtBinDEta08[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1894       fhDeltaPhiAssocPtBinDEta08[bin]->SetYTitle("#Delta #phi (rad)");
1895
1896       fhDeltaPhiAssocPtBinDEta0[bin] = new TH2F(Form("hDeltaPhiDeltaEta0PtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1897                                                 Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta = 0.", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1898                                                 nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1899       fhDeltaPhiAssocPtBinDEta0[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1900       fhDeltaPhiAssocPtBinDEta0[bin]->SetYTitle("#Delta #phi (rad)");
1901
1902       fhXEAssocPtBin[bin]       = new TH2F(Form("hXEAssocPtBin%1.f_%1.f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
1903                                            Form("#it{x}_{#it{E}} vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1904                                            nptbins, ptmin, ptmax,200, 0.0, 2.0);
1905       fhXEAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1906       fhXEAssocPtBin[bin]->SetYTitle("#it{x}_{#it{E}}");
1907
1908       fhZTAssocPtBin[bin]       = new TH2F(Form("hZTAssocPtBin%1.f_%1.f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
1909                                            Form("#it{z}_{T} vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1910                                            nptbins, ptmin, ptmax,200, 0.0, 2.0);
1911       fhZTAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1912       fhZTAssocPtBin[bin]->SetYTitle("#it{z}_{T}");
1913
1914       outputContainer->Add(fhDeltaPhiDeltaEtaAssocPtBin[bin]) ;
1915       outputContainer->Add(fhDeltaPhiAssocPtBin[bin]) ;
1916       outputContainer->Add(fhDeltaPhiAssocPtBinDEta08[bin]) ;
1917       outputContainer->Add(fhDeltaPhiAssocPtBinDEta0[bin]) ;
1918       outputContainer->Add(fhXEAssocPtBin[bin]);
1919       outputContainer->Add(fhZTAssocPtBin[bin]);
1920
1921       if(fPi0Trigger || fDecayTrigger)
1922       {
1923         fhDeltaPhiDecayChargedAssocPtBin[bin] = new TH2F(Form("hDeltaPhiPtDecayChargedAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
1924                                                          Form("#Delta #phi vs #it{p}_{T trigger} tagged as decay for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1925                                                          nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1926         fhDeltaPhiDecayChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1927         fhDeltaPhiDecayChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
1928
1929         outputContainer->Add(fhDeltaPhiDecayChargedAssocPtBin[bin]) ;
1930       }
1931
1932       if(fFillBradHisto)
1933       {
1934         fhDeltaPhiBradAssocPtBin[bin] = new TH2F(Form("hDeltaPhiBradPtAssocPt%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
1935                                                  Form("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1936                                                  nptbins, ptmin, ptmax,288, -1.0/3.0, 5.0/3.0);
1937         fhDeltaPhiBradAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1938         fhDeltaPhiBradAssocPtBin[bin]->SetYTitle("atan2(sin(#Delta #phi), cos(#Delta #phi))/#pi");
1939         outputContainer->Add(fhDeltaPhiBradAssocPtBin[bin]) ;
1940       }
1941
1942       if(fHMPIDCorrelation)
1943       {
1944         fhDeltaPhiAssocPtBinHMPID[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%sHMPID", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
1945                                                   Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, with track having HMPID signal", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1946                                                   nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1947         fhDeltaPhiAssocPtBinHMPID[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})" );
1948         fhDeltaPhiAssocPtBinHMPID[bin]->SetYTitle("#Delta #phi (rad)");
1949
1950         fhDeltaPhiAssocPtBinHMPIDAcc[bin] = new TH2F(Form("hDeltaPhiPtAssocPt%2.1f_%2.1f%sHMPIDAcc", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
1951                                                      Form("#Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, with track within 5<phi<20 deg", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
1952                                                      nptbins, ptmin, ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1953         fhDeltaPhiAssocPtBinHMPIDAcc[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
1954         fhDeltaPhiAssocPtBinHMPIDAcc[bin]->SetYTitle("#Delta #phi (rad)");
1955
1956         outputContainer->Add(fhDeltaPhiAssocPtBinHMPID   [bin]) ;
1957         outputContainer->Add(fhDeltaPhiAssocPtBinHMPIDAcc[bin]) ;
1958       }
1959     }
1960   }
1961
1962   if(fPi0Trigger || fDecayTrigger)
1963   {
1964     if(fPi0Trigger)
1965     {
1966       fhPtPi0DecayRatio  = new TH2F
1967       ("hPtPi0DecayRatio","#it{p}_{T} of #pi^{0} and the ratio of pt for two decay",
1968        nptbins,ptmin,ptmax, 100,0.,2.);
1969       fhPtPi0DecayRatio->SetXTitle("#it{p}_{T}^{#pi^{0}} (GeV/#it{c})");
1970       fhPtPi0DecayRatio->SetYTitle("#it{p}_{T}^{Decay}/#it{p}_{T}^{#pi^{0}}");
1971       outputContainer->Add(fhPtPi0DecayRatio) ;
1972     }
1973
1974     fhDeltaPhiDecayCharged  = new TH2F
1975     ("hDeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs #it{p}_{T Decay}",
1976      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
1977     fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi (rad)");
1978     fhDeltaPhiDecayCharged->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
1979
1980     fhXEDecayCharged  =
1981     new TH2F("hXEDecayCharged","#it{x}_{#it{E}}  Decay",
1982              nptbins,ptmin,ptmax,200,0.,2.);
1983     fhXEDecayCharged->SetYTitle("#it{x}_{#it{E}}");
1984     fhXEDecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
1985
1986     fhZTDecayCharged  =
1987     new TH2F("hZTDecayCharged","#it{z}_{trigger h^{#pm}} = #it{p}_{T h^{#pm}} / #it{p}_{T Decay}",
1988              nptbins,ptmin,ptmax,200,0.,2.);
1989     fhZTDecayCharged->SetYTitle("#it{z}_{decay h^{#pm}}");
1990     fhZTDecayCharged->SetXTitle("#it{p}_{T decay} (GeV/#it{c})");
1991
1992     outputContainer->Add(fhDeltaPhiDecayCharged) ;
1993     outputContainer->Add(fhXEDecayCharged) ;
1994     outputContainer->Add(fhZTDecayCharged) ;
1995   }
1996
1997   //Correlation with neutral hadrons
1998   if(fNeutralCorr)
1999   {
2000     fhDeltaPhiDeltaEtaNeutral  = new TH2F
2001     ("hDeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",
2002      ndeltaphibins ,deltaphimin,deltaphimax, ndeltaetabins ,deltaetamin,deltaetamax);
2003     fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi (rad)");
2004     fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");
2005
2006     fhPhiNeutral  = new TH2F
2007     ("hPhiNeutral","#phi_{#pi^{0}}  vs #it{p}_{T #pi^{0}}",
2008      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2009     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
2010     fhPhiNeutral->SetXTitle("#it{p}_{T #pi^{0}} (GeV/#it{c})");
2011
2012     fhEtaNeutral  = new TH2F
2013     ("hEtaNeutral","#eta_{#pi^{0}}  vs #it{p}_{T #pi^{0}}",
2014      nptbins,ptmin,ptmax,netabins,etamin,etamax);
2015     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
2016     fhEtaNeutral->SetXTitle("#it{p}_{T #pi^{0}} (GeV/#it{c})");
2017
2018     fhDeltaPhiNeutral  = new TH2F
2019     ("hDeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T trigger}",
2020      nptbins,ptmin,ptmax,nphibins,phimin,phimax);
2021     fhDeltaPhiNeutral->SetYTitle("#Delta #phi (rad)");
2022     fhDeltaPhiNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2023
2024     fhDeltaPhiNeutralPt  = new TH2F
2025     ("hDeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}}",
2026      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2027     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi (rad)");
2028     fhDeltaPhiNeutralPt->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2029
2030     fhDeltaEtaNeutral  = new TH2F
2031     ("hDeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs #it{p}_{T trigger}",
2032      nptbins,ptmin,ptmax, ndeltaetabins ,deltaetamin,deltaetamax);
2033     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");
2034     fhDeltaEtaNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2035
2036     fhXENeutral  =
2037     new TH2F("hXENeutral","#it{x}_{#it{E}} for #pi^{0} associated",
2038              nptbins,ptmin,ptmax,200,0.,2.);
2039     fhXENeutral->SetYTitle("#it{x}_{#it{E}}");
2040     fhXENeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2041
2042     fhPtHbpXENeutral  =
2043     new TH2F("hHbpXENeutral","#xi = ln(1/#it{x}_{#it{E}})for #pi^{0} associated",
2044              nptbins,ptmin,ptmax,200,0.,10.);
2045     fhPtHbpXENeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2046     fhPtHbpXENeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2047
2048     fhZTNeutral  =
2049     new TH2F("hZTNeutral","#it{z}_{trigger #pi} = #it{p}_{T #pi^{0}} / #it{p}_{T trigger} for #pi^{0} associated",
2050              nptbins,ptmin,ptmax,200,0.,2.);
2051     fhZTNeutral->SetYTitle("#it{z}_{trigger #pi^{0}}");
2052     fhZTNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2053
2054     fhPtHbpZTNeutral  =
2055     new TH2F("hHbpZTNeutral","#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2056              nptbins,ptmin,ptmax,200,0.,10.);
2057     fhPtHbpZTNeutral->SetYTitle("ln(1/#it{z}_{T})");
2058     fhPtHbpZTNeutral->SetXTitle("#it{p}_{T trigger}");
2059
2060     outputContainer->Add(fhDeltaPhiDeltaEtaNeutral);
2061     outputContainer->Add(fhPhiNeutral) ;
2062     outputContainer->Add(fhEtaNeutral) ;
2063     outputContainer->Add(fhDeltaPhiNeutral) ;
2064     outputContainer->Add(fhDeltaPhiNeutralPt) ;
2065     outputContainer->Add(fhDeltaEtaNeutral) ;
2066     outputContainer->Add(fhXENeutral) ;
2067     outputContainer->Add(fhPtHbpXENeutral) ;
2068     outputContainer->Add(fhZTNeutral) ;
2069     outputContainer->Add(fhPtHbpZTNeutral) ;
2070
2071     fhDeltaPhiUeNeutralPt  = new TH2F
2072     (Form("hDeltaPhiUe%sNeutralPt",right.Data()),"#phi_{trigger} - #phi_{#pi^{0}} vs #it{p}_{T #pi^{0}}}",
2073      nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2074     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi (rad)");
2075     fhDeltaPhiUeNeutralPt->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2076
2077     fhXEUeNeutral  =
2078     new TH2F(Form("hXEUeNeutral%s",right.Data()),"#it{x}_{#it{E}} for #pi^{0} associated",
2079              nptbins,ptmin,ptmax,200,0.,2.);
2080     fhXEUeNeutral->SetYTitle("#it{x}_{#it{E}}");
2081     fhXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2082
2083     fhPtHbpXEUeNeutral  =
2084     new TH2F(Form("hHbpXEUeNeutral%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2085              nptbins,ptmin,ptmax,200,0.,10.);
2086     fhPtHbpXEUeNeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2087     fhPtHbpXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2088
2089     fhZTUeNeutral  =
2090     new TH2F(Form("hZTUeNeutral%s",right.Data()),"#it{z}_{trigger #pi} = #it{p}_{T #pi^{0}} / #it{p}_{T trigger} for #pi^{0} associated",
2091              nptbins,ptmin,ptmax,200,0.,2.);
2092     fhZTUeNeutral->SetYTitle("#it{z}_{trigger #pi^{0}}");
2093     fhZTUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2094
2095     fhPtHbpZTUeNeutral  =
2096     new TH2F(Form("hHbpZTUeNeutral%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) for #pi^{0} associated",
2097              nptbins,ptmin,ptmax,200,0.,10.);
2098     fhPtHbpXEUeNeutral->SetYTitle("ln(1/#it{z}_{T})");
2099     fhPtHbpXEUeNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2100
2101     outputContainer->Add(fhDeltaPhiUeNeutralPt) ;
2102     outputContainer->Add(fhXEUeNeutral) ;
2103     outputContainer->Add(fhPtHbpXEUeNeutral) ;
2104     outputContainer->Add(fhZTUeNeutral) ;
2105     outputContainer->Add(fhPtHbpZTUeNeutral) ;
2106
2107     if(fMakeSeveralUE)
2108     {
2109       fhDeltaPhiUeLeftNeutral  = new TH2F
2110       ("hDeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs #it{p}_{T h^{0}} with neutral UE left side range of trigger particles",
2111        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2112       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi (rad)");
2113       fhDeltaPhiUeLeftNeutral->SetXTitle("#it{p}_{T h^{0}} (GeV/#it{c})");
2114       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;
2115
2116       fhXEUeLeftNeutral  =
2117       new TH2F("hXEUeNeutralLeft","#it{x}_{#it{E}} = #it{p}_{T Ueh^{0}} / #it{p}_{T trigger} with neutral UE left side of trigger",
2118                nptbins,ptmin,ptmax,140,0.,2.);
2119       fhXEUeLeftNeutral->SetYTitle("#it{z}_{trigger Ueh^{0}}");
2120       fhXEUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2121       outputContainer->Add(fhXEUeLeftNeutral) ;
2122
2123       fhPtHbpXEUeLeftNeutral  =
2124       new TH2F("hHbpXEUeNeutralLeft","#xi = ln(1/#it{x}_{#it{E}}) with neutral UE left side of trigger",
2125                nptbins,ptmin,ptmax,200,0.,10.);
2126       fhPtHbpXEUeLeftNeutral->SetYTitle("ln(1/#it{x}_{#it{E}})");
2127       fhPtHbpXEUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2128       outputContainer->Add(fhPtHbpXEUeLeftNeutral) ;
2129
2130       fhZTUeLeftNeutral  =
2131       new TH2F("hZTUeNeutralLeft","#it{z}_{trigger h^{0}} = #it{p}_{T Ueh^{0}} / #it{p}_{T trigger} with neutral UE left side of trigger",
2132                nptbins,ptmin,ptmax,140,0.,2.);
2133       fhZTUeLeftNeutral->SetYTitle("#it{z}_{trigger Ueh^{0}}");
2134       fhZTUeLeftNeutral->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2135       outputContainer->Add(fhZTUeLeftNeutral) ;
2136
2137       fhPtHbpZTUeLeftNeutral  =
2138       new TH2F("hHbpZTUeNeutralLeft","#xi = ln(1/#it{z}_{T}) with neutral UE left side of trigger",
2139                nptbins,ptmin,ptmax,200,0.,10.);
2140       fhPtHbpZTUeLeftNeutral->SetYTitle("ln(1/#it{z}_{T})");
2141       fhPtHbpZTUeLeftNeutral->SetXTitle("#it{p}_{T trigger}");
2142       outputContainer->Add(fhPtHbpZTUeLeftNeutral) ;
2143     }
2144
2145     if(fPi0Trigger || fDecayTrigger)
2146     {
2147       fhDeltaPhiDecayNeutral  = new TH2F
2148       ("hDeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs #it{p}_{T Decay}",
2149        nptbins,ptmin,ptmax, ndeltaphibins ,deltaphimin,deltaphimax);
2150       fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi (rad)");
2151       fhDeltaPhiDecayNeutral->SetXTitle("#it{p}_{T Decay} (GeV/#it{c})");
2152
2153       fhXEDecayNeutral  =
2154       new TH2F("hXEDecayNeutral","#it{x}_{#it{E}} for decay trigger",
2155                nptbins,ptmin,ptmax,200,0.,2.);
2156       fhXEDecayNeutral->SetYTitle("#it{x}_{#it{E}}");
2157       fhXEDecayNeutral->SetXTitle("#it{p}_{T decay}");
2158
2159       fhZTDecayNeutral  =
2160       new TH2F("hZTDecayNeutral","#it{z}_{trigger h^{0}} = #it{p}_{T h^{0}} / #it{p}_{T Decay}",
2161                nptbins,ptmin,ptmax,200,0.,2.);
2162       fhZTDecayNeutral->SetYTitle("#it{z}_{h^{0}}");
2163       fhZTDecayNeutral->SetXTitle("#it{p}_{T decay}");
2164
2165       outputContainer->Add(fhDeltaPhiDecayNeutral) ;
2166       outputContainer->Add(fhXEDecayNeutral) ;
2167       outputContainer->Add(fhZTDecayNeutral) ;
2168
2169     }
2170   }//Correlation with neutral hadrons
2171
2172   //if data is MC, fill more histograms
2173   if(IsDataMC())
2174   {
2175     fh2phiTriggerParticle = new TH2F("h2PhiTriggerParticle","#phi resolution for trigger particles",nptbins,ptmin,ptmax,100,-1,1);
2176     fh2phiTriggerParticle->GetXaxis()->SetTitle("#it{p}_{T gen Trigger} (GeV/#it{c})");
2177     fh2phiTriggerParticle->GetYaxis()->SetTitle("(#phi_{rec}-#phi_{gen})/#phi_{gen}");
2178
2179     fhMCPtTrigger  = new TH1F ("hMCPtTrigger","MC : trigger #it{p}_{T}", nptbins,ptmin,ptmax);
2180     fhMCPtTrigger->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2181
2182     fhMCPhiTrigger  = new TH2F ("hMCPhiTrigger","MC : trigger #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2183     fhMCPhiTrigger->SetYTitle("#phi (rad)");
2184     fhMCPhiTrigger->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2185
2186     fhMCEtaTrigger  = new TH2F ("hMCEtaTrigger","MC : trigger #eta",nptbins,ptmin,ptmax, netabins,etamin,etamax);
2187     fhMCEtaTrigger->SetYTitle("#eta");
2188     fhMCEtaTrigger->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2189
2190     if(fMakeAbsoluteLeading || fMakeNearSideLeading)
2191     {
2192       fhMCPtTriggerNotLeading  = new TH1F ("hMCPtTriggerNotLeading","MC : trigger #it{p}_{T}, when not leading of primaries",
2193                                            nptbins,ptmin,ptmax);
2194       fhMCPtTriggerNotLeading->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2195
2196       fhMCPhiTriggerNotLeading  = new TH2F ("hMCPhiTriggerNotLeading","MC : trigger #phi, when not leading of primaries",
2197                                             nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2198       fhMCPhiTriggerNotLeading->SetYTitle("#phi (rad)");
2199       fhMCPhiTriggerNotLeading->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2200
2201
2202       fhMCEtaTriggerNotLeading  = new TH2F ("hMCEtaTriggerNotLeading","MC : triogger #eta, when not leading of primaries",
2203                                             nptbins,ptmin,ptmax, netabins,etamin,etamax);
2204       fhMCEtaTriggerNotLeading->SetYTitle("#eta ");
2205       fhMCEtaTriggerNotLeading->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2206     }
2207
2208     fhMCEtaCharged  = new TH2F
2209     ("hMCEtaCharged","MC #eta_{h^{#pm}}  vs #it{p}_{T #pm}",
2210      nptbins,ptmin,ptmax,netabins,etamin,etamax);
2211     fhMCEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");
2212     fhMCEtaCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
2213
2214     fhMCPhiCharged  = new TH2F
2215     ("hMCPhiCharged","#MC phi_{h^{#pm}}  vs #it{p}_{T #pm}",
2216      200,ptmin,ptmax,nphibins,phimin,phimax);
2217     fhMCPhiCharged->SetYTitle("MC #phi_{h^{#pm}} (rad)");
2218     fhMCPhiCharged->SetXTitle("#it{p}_{T #pm} (GeV/#it{c})");
2219
2220     fhMCDeltaPhiDeltaEtaCharged  = new TH2F
2221     ("hMCDeltaPhiDeltaEtaCharged","#MC phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
2222      140,-2.,5.,200,-2,2);
2223     fhMCDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
2224     fhMCDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
2225
2226     fhMCDeltaEtaCharged  = new TH2F
2227     ("hMCDeltaEtaCharged","MC #eta_{trigger} - #eta_{h^{#pm}} vs #it{p}_{T trigger} and #it{p}_{T assoc}",
2228      nptbins,ptmin,ptmax,200,-2,2);
2229     fhMCDeltaEtaCharged->SetYTitle("#Delta #eta");
2230     fhMCDeltaEtaCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2231
2232     fhMCDeltaPhiCharged  = new TH2F
2233     ("hMCDeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
2234      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2235     fhMCDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
2236     fhMCDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2237
2238     fhMCDeltaPhiChargedPt  = new TH2F
2239     ("hMCDeltaPhiChargedPt","MC #phi_{trigger} - #phi_{#h^{#pm}} vs #it{p}_{T h^{#pm}}",
2240      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2241     fhMCDeltaPhiChargedPt->SetYTitle("#Delta #phi (rad)");
2242     fhMCDeltaPhiChargedPt->SetXTitle("#it{p}_{T h^{#pm}} (GeV/#it{c})");
2243
2244     fhMCPtXECharged  =
2245     new TH2F("hMCPtXECharged","#it{x}_{#it{E}} with charged hadrons",
2246              nptbins,ptmin,ptmax,200,0.,2.);
2247     fhMCPtXECharged->SetYTitle("#it{x}_{#it{E}}");
2248     fhMCPtXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2249
2250     fhMCPtHbpXECharged  =
2251     new TH2F("hMCHbpXECharged","MC #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
2252              nptbins,ptmin,ptmax,200,0.,10.);
2253     fhMCPtHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
2254     fhMCPtHbpXECharged->SetXTitle("#it{p}_{T trigger}");
2255
2256     fhMCPtZTCharged  =
2257     new TH2F("hMCPtZTCharged","#it{z}_{T} with charged hadrons",
2258              nptbins,ptmin,ptmax,200,0.,2.);
2259     fhMCPtZTCharged->SetYTitle("#it{z}_{T}");
2260     fhMCPtZTCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2261
2262     fhMCPtHbpZTCharged  =
2263     new TH2F("hMCHbpZTCharged","MC #xi = ln(1/#it{z}_{T}) with charged hadrons",
2264              nptbins,ptmin,ptmax,200,0.,10.);
2265     fhMCPtHbpZTCharged->SetYTitle("ln(1/#it{z}_{T})");
2266     fhMCPtHbpZTCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2267
2268     fhMCPtTrigPout  =
2269     new TH2F("hMCPtTrigPout","AOD MC Pout with triggers",
2270              nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax);
2271     fhMCPtTrigPout->SetYTitle("#it{p}_{out} (GeV/#it{c})");
2272     fhMCPtTrigPout->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2273
2274     fhMCPtAssocDeltaPhi  =
2275     new TH2F("hMCPtAssocDeltaPhi","AOD MC delta phi with associated charged hadrons",
2276              nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2277     fhMCPtAssocDeltaPhi->SetYTitle("#Delta #phi (rad)");
2278     fhMCPtAssocDeltaPhi->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2279
2280     outputContainer->Add(fh2phiTriggerParticle);
2281     outputContainer->Add(fhMCPtTrigger);
2282     outputContainer->Add(fhMCPhiTrigger);
2283     outputContainer->Add(fhMCEtaTrigger);
2284
2285     if(fMakeAbsoluteLeading || fMakeNearSideLeading)
2286     {
2287       outputContainer->Add(fhMCPtTriggerNotLeading);
2288       outputContainer->Add(fhMCPhiTriggerNotLeading);
2289       outputContainer->Add(fhMCEtaTriggerNotLeading);
2290     }
2291
2292     outputContainer->Add(fhMCDeltaPhiDeltaEtaCharged);
2293     outputContainer->Add(fhMCPhiCharged) ;
2294     outputContainer->Add(fhMCEtaCharged) ;
2295     outputContainer->Add(fhMCDeltaEtaCharged) ;
2296     outputContainer->Add(fhMCDeltaPhiCharged) ;
2297
2298     outputContainer->Add(fhMCDeltaPhiChargedPt) ;
2299     outputContainer->Add(fhMCPtXECharged) ;
2300     outputContainer->Add(fhMCPtZTCharged) ;
2301     outputContainer->Add(fhMCPtHbpXECharged) ;
2302     outputContainer->Add(fhMCPtHbpZTCharged) ;
2303     outputContainer->Add(fhMCPtTrigPout) ;
2304     outputContainer->Add(fhMCPtAssocDeltaPhi) ;
2305
2306     // Underlying event
2307
2308     fhMCUePart  =
2309     new TH1F("hMCUePart","MC UE particles distribution vs pt trig",
2310              nptbins,ptmin,ptmax);
2311     fhMCUePart->SetYTitle("#it{dN}^{ch}");
2312     fhMCUePart->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2313
2314     fhMCPtXEUeCharged  =
2315     new TH2F(Form("hMCPtXEUeCharged%s",right.Data()),"#it{x}_{#it{E}} with charged hadrons, Underlying Event",
2316              nptbins,ptmin,ptmax,200,0.,2.);
2317     fhMCPtXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
2318     fhMCPtXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2319
2320     fhMCPtHbpXEUeCharged =
2321     new TH2F(Form("hMCPtHbpXEUeCharged%s",right.Data()),"#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, Underlying Event",
2322              nptbins,ptmin,ptmax,200,0.,10.);
2323     fhMCPtHbpXEUeCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
2324     fhMCPtHbpXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2325
2326     fhMCPtZTUeCharged =
2327     new TH2F(Form("hMCPtZTUeCharged%s",right.Data()),"#it{z}_{T} with charged hadrons, Underlying Event",
2328              nptbins,ptmin,ptmax,200,0.,2.);
2329     fhMCPtZTUeCharged->SetYTitle("#it{z}_{T}");
2330     fhMCPtZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2331
2332     fhMCPtHbpZTUeCharged =
2333     new TH2F(Form("hMCPtHbpZTUeCharged%s",right.Data()),"#xi = ln(1/#it{z}_{T}) with charged hadrons, Underlying Event",
2334              nptbins,ptmin,ptmax,200,0.,10.);
2335     fhMCPtHbpZTUeCharged->SetYTitle("ln(1/#it{z}_{T})");
2336     fhMCPtHbpZTUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2337
2338     outputContainer->Add(fhMCUePart);
2339     outputContainer->Add(fhMCPtXEUeCharged) ;
2340     outputContainer->Add(fhMCPtZTUeCharged) ;
2341     outputContainer->Add(fhMCPtHbpZTUeCharged);
2342     outputContainer->Add(fhMCPtHbpXEUeCharged);
2343
2344     if(fMakeSeveralUE)
2345     {
2346       fhMCPtXEUeLeftCharged  =
2347       new TH2F("hMCPtXEUeChargedLeft","#it{x}_{#it{E}} with charged hadrons, with UE left side range of trigger particles",
2348                nptbins,ptmin,ptmax,200,0.,2.);
2349       fhMCPtXEUeLeftCharged->SetYTitle("#it{x}_{#it{E}}");
2350       fhMCPtXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2351
2352
2353       fhMCPtHbpXEUeLeftCharged =
2354       new TH2F("hMCPtHbpXEUeChargedLeft","#xi = ln(1/#it{x}_{#it{E}}) with charged hadrons, with UE left side range of trigger particles",
2355                nptbins,ptmin,ptmax,200,0.,10.);
2356       fhMCPtHbpXEUeLeftCharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
2357       fhMCPtHbpXEUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2358
2359       fhMCPtZTUeLeftCharged  =
2360       new TH2F("hMCPtZTUeChargedLeft","#it{z}_{T} with charged hadrons, with UE left side range of trigger particles",
2361                nptbins,ptmin,ptmax,200,0.,2.);
2362       fhMCPtZTUeLeftCharged->SetYTitle("#it{z}_{T}");
2363       fhMCPtZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2364
2365       fhMCPtHbpZTUeLeftCharged =
2366       new TH2F("hMCPtHbpZTUeChargedLeft","#xi = ln(1/#it{z}_{T}) with charged hadrons, with UE left side range of trigger particles",
2367                nptbins,ptmin,ptmax,200,0.,10.);
2368       fhMCPtHbpZTUeLeftCharged->SetYTitle("ln(1/#it{z}_{T})");
2369       fhMCPtHbpZTUeLeftCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2370
2371       outputContainer->Add(fhMCPtXEUeLeftCharged) ;
2372       outputContainer->Add(fhMCPtZTUeLeftCharged) ;
2373       outputContainer->Add(fhMCPtHbpXEUeLeftCharged);
2374       outputContainer->Add(fhMCPtHbpZTUeLeftCharged) ;
2375
2376     }
2377   } //for MC histogram
2378
2379   if(DoOwnMix())
2380   {
2381     //create event containers
2382
2383     if(!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForTracksExists()))
2384     {
2385       Int_t nvz = GetNZvertBin();
2386       Int_t nrp = GetNRPBin();
2387       Int_t nce = GetNCentrBin();
2388
2389       fListMixTrackEvents= new TList*[nvz*nrp*nce] ;
2390
2391       for( Int_t ice = 0 ; ice < nce ; ice++ )
2392       {
2393         for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
2394         {
2395           for( Int_t irp = 0 ; irp < nrp ; irp++ )
2396           {
2397             Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
2398
2399             //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
2400             //       ic,iz, irp, bin);
2401
2402             fListMixTrackEvents[bin] = new TList() ;
2403             fListMixTrackEvents[bin]->SetOwner(kFALSE);
2404           }
2405         }
2406       }
2407     }
2408
2409     fhPtTriggerMixed  = new TH1F ("hPtTriggerMixed","#it{p}_{T} distribution of trigger particles, used for mixing", nptbins,ptmin,ptmax);
2410     fhPtTriggerMixed->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2411
2412     if(fCorrelVzBin)
2413     {
2414       fhPtTriggerMixedVzBin  = new TH2F ("hPtTriggerMixedVzBin","#it{p}_{T} distribution of trigger particles, used for mixing", nptbins,ptmin,ptmax,GetNZvertBin(),0,GetNZvertBin());
2415       fhPtTriggerMixedVzBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2416       fhPtTriggerMixedVzBin->SetYTitle("#it{v}_{#it{z}} bin");
2417       outputContainer->Add(fhPtTriggerMixedVzBin);
2418     }
2419
2420     fhPtTriggerMixedBin  = new TH2F ("hPtTriggerMixedBin","#it{p}_{T} distribution of trigger particles vs mixing bin", nptbins,ptmin,ptmax,nMixBins,0,nMixBins);
2421     fhPtTriggerMixedBin->SetXTitle("#it{p}_{T}^{trig} (GeV/#it{c})");
2422     fhPtTriggerMixedBin->SetYTitle("Bin");
2423
2424     fhPhiTriggerMixed  = new TH2F ("hPhiTriggerMixed","#phi distribution of trigger Particles, used for mixing",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
2425     fhPhiTriggerMixed->SetYTitle("#phi (rad)");
2426
2427     fhEtaTriggerMixed  = new TH2F ("hEtaTriggerMixed","#eta distribution of trigger, used for mixing",nptbins,ptmin,ptmax, netabins,etamin,etamax);
2428     fhEtaTriggerMixed->SetYTitle("#eta ");
2429
2430     outputContainer->Add(fhPtTriggerMixed);
2431     outputContainer->Add(fhPtTriggerMixedBin);
2432     outputContainer->Add(fhPhiTriggerMixed);
2433     outputContainer->Add(fhEtaTriggerMixed);
2434
2435     // Fill the cluster pool only in isolation analysis or if requested
2436     if( ( OnlyIsolated()        ||  fFillNeutralEventMixPool) &&
2437        (!fUseMixStoredInReader || (fUseMixStoredInReader && !GetReader()->ListWithMixedEventsForCaloExists())))
2438     {
2439       Int_t nvz = GetNZvertBin();
2440       Int_t nrp = GetNRPBin();
2441       Int_t nce = GetNCentrBin();
2442
2443       fListMixCaloEvents= new TList*[nvz*nrp*nce] ;
2444
2445       for( Int_t ice = 0 ; ice < nce ; ice++ )
2446       {
2447         for( Int_t ivz = 0 ; ivz < nvz ; ivz++ )
2448         {
2449           for( Int_t irp = 0 ; irp < nrp ; irp++ )
2450           {
2451             Int_t bin = GetEventMixBin(ice,ivz,irp); //ic*nvz*nrp+iz*nrp+irp;
2452
2453             //printf("GetCreateOutputObjects - Bins : cent %d, vz %d, RP %d, event %d/%d\n",
2454             //       ic,iz, irp, bin);
2455
2456             fListMixCaloEvents[bin] = new TList() ;
2457             fListMixCaloEvents[bin]->SetOwner(kFALSE);
2458           }
2459         }
2460       }
2461     }
2462
2463     //Init the list in the reader if not done previously
2464     if(fUseMixStoredInReader)
2465     {
2466       if( !GetReader()->ListWithMixedEventsForTracksExists() )
2467         GetReader()->SetListWithMixedEventsForTracks(fListMixTrackEvents);
2468
2469       if( !GetReader()->ListWithMixedEventsForCaloExists()   )
2470         GetReader()->SetListWithMixedEventsForCalo  (fListMixCaloEvents );
2471     }
2472
2473     fhEventBin=new TH1I("hEventBin","Number of real events per bin(cen,vz,rp)",
2474                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2475                         GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2476     fhEventBin->SetXTitle("bin");
2477     outputContainer->Add(fhEventBin) ;
2478
2479     fhEventMixBin=new TH1I("hEventMixBin","Number of events  per bin(cen,vz,rp)",
2480                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1,0,
2481                            GetNCentrBin()*GetNZvertBin()*GetNRPBin()+1) ;
2482     fhEventMixBin->SetXTitle("bin");
2483     outputContainer->Add(fhEventMixBin) ;
2484
2485     fhNtracksMB=new TH1F("hNtracksMBEvent","Number of tracks w/ event trigger kMB",2000,0,2000);
2486     outputContainer->Add(fhNtracksMB);
2487
2488     if(fFillNeutralEventMixPool || OnlyIsolated())
2489     {
2490       fhNclustersMB=new TH1F("hNclustersMBEvent","Number of clusters w/ event trigger kMB",2000,0,2000);
2491       outputContainer->Add(fhNclustersMB);
2492     }
2493
2494     fhMixDeltaPhiCharged  = new TH2F
2495     ("hMixDeltaPhiCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #it{p}_{T trigger}",
2496      nptbins,ptmin,ptmax,ndeltaphibins ,deltaphimin,deltaphimax);
2497     fhMixDeltaPhiCharged->SetYTitle("#Delta #phi (rad)");
2498     fhMixDeltaPhiCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2499     outputContainer->Add(fhMixDeltaPhiCharged);
2500
2501     fhMixDeltaPhiDeltaEtaCharged  = new TH2F
2502     ("hMixDeltaPhiDeltaEtaCharged","Mixed event : #phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",
2503      ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
2504     fhMixDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi (rad)");
2505     fhMixDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");
2506     outputContainer->Add(fhMixDeltaPhiDeltaEtaCharged);
2507
2508     fhMixXECharged  =
2509     new TH2F("hMixXECharged","Mixed event : #it{x}_{#it{E}} for charged tracks",
2510              nptbins,ptmin,ptmax,200,0.,2.);
2511     fhMixXECharged->SetYTitle("#it{x}_{#it{E}}");
2512     fhMixXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2513     outputContainer->Add(fhMixXECharged);
2514
2515     fhMixXEUeCharged  =
2516     new TH2F("hMixXEUeCharged","Mixed event : #it{x}_{#it{E}} for charged tracks in Ue region",
2517              nptbins,ptmin,ptmax,200,0.,2.);
2518     fhMixXEUeCharged->SetYTitle("#it{x}_{#it{E}}");
2519     fhMixXEUeCharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2520     outputContainer->Add(fhMixXEUeCharged);
2521
2522     fhMixHbpXECharged  =
2523     new TH2F("hMixHbpXECharged","mixed event : #xi = ln(1/#it{x}_{#it{E}}) with charged hadrons",
2524              nptbins,ptmin,ptmax,200,0.,10.);
2525     fhMixHbpXECharged->SetYTitle("ln(1/#it{x}_{#it{E}})");
2526     fhMixHbpXECharged->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2527     outputContainer->Add(fhMixHbpXECharged);
2528
2529     fhMixDeltaPhiChargedAssocPtBin         = new TH2F*[fNAssocPtBins*nz];
2530     fhMixDeltaPhiChargedAssocPtBinDEta08   = new TH2F*[fNAssocPtBins*nz];
2531     fhMixDeltaPhiChargedAssocPtBinDEta0    = new TH2F*[fNAssocPtBins*nz];
2532     fhMixDeltaPhiDeltaEtaChargedAssocPtBin = new TH2F*[fNAssocPtBins*nz];
2533
2534     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
2535     {
2536       for(Int_t z = 0 ; z < nz ; z++)
2537       {
2538         Int_t bin = i*nz+z;
2539
2540         if(fCorrelVzBin)
2541         {
2542           sz = Form("_vz%d",z);
2543           tz = Form(", #it{v}_{#it{z}} bin %d",z);
2544         }
2545
2546         //printf("MIX : iAssoc %d, Vz %d, bin %d - sz %s, tz %s \n",i,z,bin,sz.Data(),tz.Data());
2547
2548         fhMixDeltaPhiChargedAssocPtBin[bin] = new TH2F(Form("hMixDeltaPhiChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2549                                                        Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2550                                                        nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
2551         fhMixDeltaPhiChargedAssocPtBin[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2552         fhMixDeltaPhiChargedAssocPtBin[bin]->SetYTitle("#Delta #phi (rad)");
2553
2554         fhMixDeltaPhiChargedAssocPtBinDEta08[bin] = new TH2F(Form("hMixDeltaPhiDeltaEta0.8ChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2555                                                              Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta > 0.8", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2556                                                              nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
2557         fhMixDeltaPhiChargedAssocPtBinDEta08[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2558         fhMixDeltaPhiChargedAssocPtBinDEta08[bin]->SetYTitle("#Delta #phi (rad)");
2559
2560         fhMixDeltaPhiChargedAssocPtBinDEta0[bin] = new TH2F(Form("hMixDeltaPhiDeltaEta0ChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2561                                                             Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s, for #Delta #eta = 0", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2562                                                             nptbins, ptmin, ptmax,  ndeltaphibins ,deltaphimin,deltaphimax);
2563         fhMixDeltaPhiChargedAssocPtBinDEta0[bin]->SetXTitle("#it{p}_{T trigger} (GeV/#it{c})");
2564         fhMixDeltaPhiChargedAssocPtBinDEta0[bin]->SetYTitle("#Delta #phi (rad)");
2565
2566         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin] = new TH2F(Form("hMixDeltaPhiDeltaEtaChargedAssocPtBin%2.1f_%2.1f%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],sz.Data()),
2567                                                                Form("Mixed event #Delta #phi vs #it{p}_{T trigger} for associated #it{p}_{T} bin [%2.1f,%2.1f]%s", fAssocPtBinLimit[i], fAssocPtBinLimit[i+1],tz.Data()),
2568                                                                ndeltaphibins ,deltaphimin,deltaphimax,ndeltaetabins ,deltaetamin,deltaetamax);
2569         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->SetXTitle("#Delta #phi (rad)");
2570         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->SetYTitle("#Delta #eta");
2571
2572         outputContainer->Add(fhMixDeltaPhiChargedAssocPtBin[bin]);
2573         outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta08[bin]);
2574         outputContainer->Add(fhMixDeltaPhiChargedAssocPtBinDEta0[bin]);
2575         outputContainer->Add(fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]);
2576       }
2577     }
2578   }
2579
2580   return outputContainer;
2581
2582 }
2583
2584 //_________________________________________________________________________________________________
2585 Bool_t AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum(AliAODPWG4Particle* trigger,
2586                                                                TLorentzVector & mom1,
2587                                                                TLorentzVector & mom2)
2588 {
2589   // Get the momentum of the pi0/eta assigned decay photons
2590   // In case of pi0/eta trigger, we may want to check their decay correlation,
2591   // get their decay children
2592
2593   Int_t indexPhoton1 = trigger->GetCaloLabel(0);
2594   Int_t indexPhoton2 = trigger->GetCaloLabel(1);
2595
2596   if(indexPhoton1!=-1 || indexPhoton2!=-1) return kFALSE;
2597
2598   if(GetDebug() > 1)
2599     printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);
2600
2601   TObjArray * clusters  = 0x0 ;
2602   if(trigger->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;
2603   else                                clusters = GetPHOSClusters()  ;
2604
2605   for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++)
2606   {
2607     AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));
2608
2609     if(photon->GetID()==indexPhoton1) photon->GetMomentum(mom1,GetVertex(0)) ;
2610     if(photon->GetID()==indexPhoton2) photon->GetMomentum(mom1,GetVertex(0)) ;
2611
2612     if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::GetDecayPhotonMomentum() - Photon1 = %f, Photon2 = %f \n", mom1.Pt(), mom2.Pt());
2613
2614   } //cluster loop
2615
2616   return kTRUE;
2617
2618 }
2619
2620 //_____________________________________________________________
2621 Int_t AliAnaParticleHadronCorrelation::GetMCTagHistogramIndex(Int_t mcTag)
2622 {
2623   // Index of MC histograms depending on MC origin
2624
2625   if     ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPrompt) ||
2626            GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation)) return 0;
2627   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0))           return 1;
2628   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0Decay))      return 2;
2629   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEtaDecay))      return 3;
2630   else if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCOtherDecay))    return 4;
2631   else if(!GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron))      return 5;
2632   else                                                                                    return 6;
2633
2634 }
2635
2636 //_________________________________________
2637 void AliAnaParticleHadronCorrelation::Init()
2638 {
2639   //Init
2640   //Do some checks
2641
2642   if(!GetReader()->IsCTSSwitchedOn())
2643     AliFatal("STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!\n");
2644 }
2645
2646 //____________________________________________________
2647 void AliAnaParticleHadronCorrelation::InitParameters()
2648 {
2649
2650   //Initialize the parameters of the analysis.
2651   SetInputAODName("Particle");
2652   SetAODObjArrayName("Hadrons");
2653   AddToHistogramsName("AnaHadronCorr_");
2654
2655   SetPtCutRange(0.,300);
2656   fDeltaPhiMinCut       = TMath::DegToRad()*120.;
2657   fDeltaPhiMaxCut       = TMath::DegToRad()*240. ;
2658   fSelectIsolated       = kFALSE;
2659   fMakeSeveralUE        = kFALSE;
2660   fUeDeltaPhiMinCut     = TMath::DegToRad()*60.;
2661   fUeDeltaPhiMaxCut     = TMath::DegToRad()*120 ;
2662
2663   fNeutralCorr          = kFALSE ;
2664   fPi0Trigger           = kFALSE ;
2665   fDecayTrigger         = kFALSE ;
2666   fHMPIDCorrelation     = kFALSE ;
2667
2668   fMakeAbsoluteLeading  = kTRUE;
2669   fMakeNearSideLeading  = kFALSE;
2670
2671   fNAssocPtBins         = 9   ;
2672   fAssocPtBinLimit[0]   = 0.2 ;
2673   fAssocPtBinLimit[1]   = 0.5 ;
2674   fAssocPtBinLimit[2]   = 1.0 ;
2675   fAssocPtBinLimit[3]   = 2.0 ;
2676   fAssocPtBinLimit[4]   = 3.0 ;
2677   fAssocPtBinLimit[5]   = 4.0 ;
2678   fAssocPtBinLimit[6]   = 5.0 ;
2679   fAssocPtBinLimit[7]   = 6.0 ;
2680   fAssocPtBinLimit[8]   = 7.0 ;
2681   fAssocPtBinLimit[9]   = 8.0 ;
2682   fAssocPtBinLimit[10]  = 9.0 ;
2683   fAssocPtBinLimit[11]  = 10.0 ;
2684   fAssocPtBinLimit[12]  = 12.0 ;
2685   fAssocPtBinLimit[13]  = 14.0 ;
2686   fAssocPtBinLimit[14]  = 16.0 ;
2687   fAssocPtBinLimit[15]  = 20.0 ;
2688   fAssocPtBinLimit[16]  = 30.0 ;
2689   fAssocPtBinLimit[17]  = 40.0 ;
2690   fAssocPtBinLimit[18]  = 50.0 ;
2691   fAssocPtBinLimit[19]  = 200.0 ;
2692
2693   fUseMixStoredInReader = kTRUE;
2694
2695   fM02MinCut   = -1 ;
2696   fM02MaxCut   = -1 ;
2697
2698   fSelectLeadingHadronAngle = kFALSE;
2699   fMinLeadHadPhi = 150*TMath::DegToRad();
2700   fMaxLeadHadPhi = 210*TMath::DegToRad();
2701
2702   fMinLeadHadPt  = 1;
2703   fMaxLeadHadPt  = 100;
2704
2705 }
2706
2707 //_________________________________________________________________________
2708 Bool_t AliAnaParticleHadronCorrelation::IsTriggerTheEventLeadingParticle()
2709 {
2710   // Check if the what of the selected triggers is leading particle comparing
2711   // with all the triggers, all the tracks or all the clusters (if requested for the clusters).
2712
2713   Double_t ptTrig      = GetMinPt();
2714   Double_t phiTrig     = 0 ;
2715   fLeadingTriggerIndex =-1 ;
2716   Int_t index          =-1 ;
2717   AliAODPWG4ParticleCorrelation* pLeading = 0;
2718
2719   // Loop on stored AOD particles, find leading trigger on the selected list, with at least min pT.
2720
2721   for(Int_t iaod = 0; iaod < GetInputAODBranch()->GetEntriesFast() ; iaod++)
2722   {
2723     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
2724     particle->SetLeadingParticle(kFALSE); // set it later
2725
2726     // Vertex cut in case of mixing
2727     Int_t check = CheckMixedEventVertex(particle->GetCaloLabel(0), particle->GetTrackLabel(0));
2728     if(check ==  0) continue;
2729     if(check == -1) return kFALSE; // not sure if it is correct.
2730
2731     // find the leading particles with highest momentum
2732     if (particle->Pt() > ptTrig)
2733     {
2734       ptTrig   = particle->Pt() ;
2735       phiTrig  = particle->Phi();
2736       index    = iaod     ;
2737       pLeading = particle ;
2738     }
2739   }// finish search of leading trigger particle on the AOD branch.
2740
2741   if(index < 0) return kFALSE;
2742
2743   //printf("AOD leading pT %2.2f, ID %d\n", pLeading->Pt(),pLeading->GetCaloLabel(0));
2744
2745   if(phiTrig < 0 ) phiTrig += TMath::TwoPi();
2746
2747   // Compare if it is the leading of all tracks
2748
2749   TVector3 p3;
2750   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
2751   {
2752     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
2753
2754     if(track->GetID() == pLeading->GetTrackLabel(0) || track->GetID() == pLeading->GetTrackLabel(1) ||
2755        track->GetID() == pLeading->GetTrackLabel(2) || track->GetID() == pLeading->GetTrackLabel(3)   ) continue ;
2756
2757     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
2758     p3.SetXYZ(mom[0],mom[1],mom[2]);
2759     Float_t pt   = p3.Pt();
2760     Float_t phi  = p3.Phi() ;
2761     if(phi < 0) phi+=TMath::TwoPi();
2762
2763     //jump out this event if near side associated particle pt larger than trigger
2764     if (fMakeNearSideLeading)
2765     {
2766       if(pt > ptTrig && TMath::Abs(phi-phiTrig) < TMath::PiOver2())  return kFALSE;
2767     }
2768     //jump out this event if there is any other particle with pt larger than trigger
2769     else
2770     {
2771       if(pt > ptTrig)  return kFALSE ;
2772     }
2773    }// track loop
2774
2775   // Compare if it is leading of all calorimeter clusters
2776
2777   if(fCheckLeadingWithNeutralClusters)
2778   {
2779     // Select the calorimeter cluster list
2780     TObjArray * nePl = 0x0;
2781     if      (pLeading->GetDetector() == "PHOS" )
2782       nePl = GetPHOSClusters();
2783     else
2784       nePl = GetEMCALClusters();
2785
2786     if(!nePl) return kTRUE; // Do the selection just with the tracks if no calorimeter is available.
2787
2788     TLorentzVector lv;
2789     for(Int_t ipr = 0;ipr < nePl->GetEntriesFast() ; ipr ++ )
2790     {
2791       AliVCluster * cluster = (AliVCluster *) (nePl->At(ipr)) ;
2792
2793       if(cluster->GetID() == pLeading->GetCaloLabel(0) || cluster->GetID() == pLeading->GetCaloLabel(1) ) continue ;
2794
2795       cluster->GetMomentum(lv,GetVertex(0));
2796
2797       Float_t pt   = lv.Pt();
2798       Float_t phi  = lv.Phi() ;
2799       if(phi < 0) phi+=TMath::TwoPi();
2800
2801       if(IsTrackMatched(cluster,GetReader()->GetInputEvent())) continue ; // avoid charged clusters, already covered by tracks, or cluster merging with track.
2802
2803       //jump out this event if near side associated particle pt larger than trigger
2804       // not really needed for calorimeter, unless DCal is included
2805       if (fMakeNearSideLeading)
2806       {
2807         if(pt > ptTrig && TMath::Abs(phi-phiTrig) < TMath::PiOver2()) return kFALSE ;
2808       }
2809       //jump out this event if there is any other particle with pt larger than trigger
2810       else
2811       {
2812         if(pt > ptTrig)  return kFALSE ;
2813       }
2814     }// cluster loop
2815   } // check neutral clusters
2816
2817   fLeadingTriggerIndex = index ;
2818   pLeading->SetLeadingParticle(kTRUE);
2819
2820   if( GetDebug() > 1 ) printf("\t particle AOD with index %d is leading with pT %2.2f\n", fLeadingTriggerIndex, pLeading->Pt());
2821
2822   return kTRUE;
2823
2824 }
2825
2826 //_________________________________________________________________
2827 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()
2828 {
2829   //Particle-Hadron Correlation Analysis, fill histograms
2830
2831   if(!GetInputAODBranch())
2832   {
2833     AliFatal(Form("No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data()));
2834     return ; // coverity
2835   }
2836
2837   Int_t naod = GetInputAODBranch()->GetEntriesFast();
2838   if( naod == 0 )
2839   {
2840     if(GetDebug() > 1)
2841       printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No particle AOD found! \n");
2842
2843     return ; // no trigger particles found.
2844   }
2845
2846   if(GetDebug() > 1)
2847   {
2848     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");
2849     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", naod);
2850     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In CTS aod entries %d\n",   GetCTSTracks()->GetEntriesFast());
2851   }
2852
2853   //------------------------------------------------------
2854   // Find leading trigger if analysis request only leading,
2855   // if there is no leading trigger, then skip the event
2856
2857   Int_t iaod = 0 ;
2858   if( fMakeAbsoluteLeading || fMakeNearSideLeading )
2859   {
2860     Bool_t leading = IsTriggerTheEventLeadingParticle();
2861
2862     if(GetDebug() > 1)
2863       printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - AOD Leading trigger? %d, with index %d\n",leading,fLeadingTriggerIndex);
2864
2865     if(!leading)
2866     {
2867       if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Leading was requested and not found\n");
2868       return ;
2869     }
2870     else
2871     {
2872       // Select only the leading in the trigger AOD loop
2873       naod = 1 ;
2874       iaod = fLeadingTriggerIndex;
2875     }
2876   }
2877
2878   //------------------------------------------------------
2879   // Get event multiplicity and bins
2880
2881   Float_t cen         = GetEventCentrality();
2882   Float_t ep          = GetEventPlaneAngle();
2883   if(fFillHighMultHistograms) fhTriggerEventPlaneCentrality->Fill(cen,ep);
2884
2885   Int_t   mixEventBin = GetEventMixBin();
2886   Int_t   vzbin       = GetEventVzBin();
2887
2888   //------------------------------------------------------
2889   // Loop on trigger AOD
2890
2891   for( iaod = 0; iaod < naod; iaod++ )
2892   {
2893     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
2894
2895     //
2896     // Trigger particle selection criteria:
2897     //
2898     Float_t pt = particle->Pt();
2899
2900     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
2901
2902     fhPtTriggerInput->Fill(pt);
2903
2904     //
2905     // check if it was a calorimeter cluster
2906     // and if the shower shape cut was requested apply it.
2907     // Not needed if already done at the particle identification level,
2908     // but for isolation studies, it is preferred not to remove so we do it here
2909     //
2910     Int_t clID1  = particle->GetCaloLabel(0) ;
2911     Int_t clID2  = particle->GetCaloLabel(1) ; // for photon clusters should not be set.
2912     if( GetDebug() > 1 ) printf("%s Trigger : id1 %d, id2 %d, min %f, max %f, det %s\n",
2913            GetInputAODName().Data(),clID1,clID2,fM02MinCut,fM02MaxCut,(particle->GetDetector()).Data());
2914
2915     if(clID1 > 0 && clID2 < 0 && fM02MaxCut > 0 && fM02MinCut > 0)
2916     {
2917       Int_t iclus = -1;
2918       TObjArray* clusters = 0x0;
2919       if     (particle->GetDetector() == "EMCAL") clusters = GetEMCALClusters();
2920       else if(particle->GetDetector() == "PHOS" ) clusters = GetPHOSClusters();
2921
2922       if(clusters)
2923       {
2924         AliVCluster *cluster = FindCluster(clusters,clID1,iclus);
2925         Float_t m02 = cluster->GetM02();
2926         if(m02 > fM02MaxCut || m02 < fM02MinCut) continue ;
2927       }
2928
2929       fhPtTriggerSSCut->Fill(pt);
2930     }
2931
2932     //
2933     // Check if the particle is isolated or if we want to take the isolation into account
2934     // This bool is set in AliAnaParticleIsolation
2935     //
2936     if(OnlyIsolated())
2937     {
2938       if( !particle->IsIsolated() ) continue;
2939       fhPtTriggerIsoCut->Fill(pt);
2940     }
2941
2942     //
2943     // Check if trigger is in fiducial region
2944     //
2945     if(IsFiducialCutOn())
2946     {
2947       Bool_t in = GetFiducialCut()->IsInFiducialCut(*particle->Momentum(),particle->GetDetector()) ;
2948       if(! in ) continue ;
2949     }
2950
2951     fhPtTriggerFidCut->Fill(pt);
2952
2953     //---------------------------------------
2954     // Make correlation
2955
2956     // Find the leading hadron in the opposite hemisphere to the triggeer
2957     // and accept the trigger if leading is in defined window.
2958     Bool_t okLeadHad = kTRUE;
2959     if(fSelectLeadingHadronAngle)
2960       okLeadHad = FindLeadingOppositeHadronInWindow(particle);
2961     if(!okLeadHad) continue;
2962
2963     //
2964     // Charged particles correlation
2965     //
2966     MakeChargedCorrelation(particle);
2967
2968     // MC
2969     if(IsDataMC())
2970       MakeMCChargedCorrelation(particle->GetLabel());
2971
2972     // Do own mixed event with charged,
2973     // add event and remove previous or fill the mixed histograms
2974     if(DoOwnMix())
2975       MakeChargedMixCorrelation(particle);
2976
2977     //
2978     // Neutral particles correlation
2979     //
2980     if(fNeutralCorr)
2981       MakeNeutralCorrelation(particle);
2982
2983     //----------------------------------------------------------------
2984     // Fill trigger pT related histograms if not absolute leading
2985
2986     //
2987     // pT of the trigger, vs trigger origin if MC
2988     //
2989     fhPtTrigger->Fill(pt);
2990     if(IsDataMC())
2991     {
2992       Int_t mcIndex = GetMCTagHistogramIndex(particle->GetTag());
2993       fhPtTriggerMC[mcIndex]->Fill(pt);
2994     }
2995
2996     //
2997     // Acceptance of the trigger
2998     //
2999     Float_t phi = particle->Phi();
3000     if( phi<0 ) phi+=TMath::TwoPi();
3001     fhPhiTrigger->Fill(pt, phi);
3002
3003     fhEtaTrigger->Fill(pt, particle->Eta());
3004     //printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Trigger particle : pt %f, eta %f, phi %f\n",particle->Pt(),particle->Eta(),phi);
3005
3006     //----------------------------------
3007     // Trigger particle pT vs event bins
3008
3009     fhPtTriggerBin->Fill(pt,mixEventBin);
3010     if(fCorrelVzBin)
3011       fhPtTriggerVzBin->Fill(pt,vzbin);
3012
3013     if(fFillHighMultHistograms)
3014     {
3015       fhPtTriggerCentrality->Fill(pt,cen);
3016       fhPtTriggerEventPlane->Fill(pt,ep);
3017     }
3018
3019     //----------------------------------
3020     // Trigger particle pT vs pile-up
3021
3022     if(fFillPileUpHistograms)
3023     {
3024       Int_t vtxBC = GetReader()->GetVertexBC();
3025       if(vtxBC == 0 || vtxBC==AliVTrack::kTOFBCNA)     fhPtTriggerVtxBC0->Fill(pt);
3026
3027       if(GetReader()->IsPileUpFromSPD())               fhPtTriggerPileUp[0]->Fill(pt);
3028       if(GetReader()->IsPileUpFromEMCal())             fhPtTriggerPileUp[1]->Fill(pt);
3029       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtTriggerPileUp[2]->Fill(pt);
3030       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtTriggerPileUp[3]->Fill(pt);
3031       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtTriggerPileUp[4]->Fill(pt);
3032       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtTriggerPileUp[5]->Fill(pt);
3033       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtTriggerPileUp[6]->Fill(pt);
3034     }
3035   } // AOD trigger loop
3036
3037   //Reinit for next event
3038   fLeadingTriggerIndex = -1;
3039
3040   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");
3041 }
3042
3043 //_______________________________________________________________________________________________________
3044 void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
3045 {
3046   // Charged Hadron Correlation Analysis
3047   if(GetDebug() > 1)
3048     printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");
3049
3050   Float_t phiTrig = aodParticle->Phi();
3051   Float_t etaTrig = aodParticle->Eta();
3052   Float_t ptTrig  = aodParticle->Pt();
3053   Bool_t   decay  = aodParticle->IsTagged();
3054   Int_t    mcTag  = aodParticle->GetTag();
3055   Double_t bz     = GetReader()->GetInputEvent()->GetMagneticField();
3056
3057   Float_t pt       = -100. ;
3058   Float_t phi      = -100. ;
3059   Float_t eta      = -100. ;
3060   Float_t deltaPhi = -100. ;
3061
3062   TVector3 p3;
3063   TLorentzVector photonMom ;
3064   TObjArray * reftracks = 0x0;
3065   Int_t nrefs           = 0;
3066
3067   // Mixed event settings
3068   Int_t evtIndex11   = -1 ; // cluster trigger or pi0 trigger
3069   Int_t evtIndex12   = -1 ; // pi0 trigger
3070   Int_t evtIndex13   = -1 ; // charged trigger
3071
3072   if (GetMixedEvent())
3073   {
3074     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
3075     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
3076     evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;
3077   }
3078
3079   //
3080   // In case of pi0/eta trigger, we may want to check their decay correlation,
3081   // get their decay children
3082   //
3083   TLorentzVector decayMom1;
3084   TLorentzVector decayMom2;
3085   Bool_t decayFound = kFALSE;
3086   if( fPi0Trigger )
3087   {
3088     decayFound = GetDecayPhotonMomentum(aodParticle,decayMom1, decayMom2);
3089     if(decayFound)
3090     {
3091       fhPtPi0DecayRatio->Fill(ptTrig, decayMom1.Pt()/ptTrig);
3092       fhPtPi0DecayRatio->Fill(ptTrig, decayMom2.Pt()/ptTrig);
3093     }
3094   }
3095
3096   //-----------------------------------------------------------------------
3097   // Track loop, select tracks with good pt, phi and fill AODs or histograms
3098   //-----------------------------------------------------------------------
3099
3100   for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ )
3101   {
3102     AliVTrack * track = (AliVTrack *) (GetCTSTracks()->At(ipr)) ;
3103
3104     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
3105     p3.SetXYZ(mom[0],mom[1],mom[2]);
3106     pt   = p3.Pt();
3107     eta  = p3.Eta();
3108     phi  = p3.Phi() ;
3109     if(phi < 0) phi+=TMath::TwoPi();
3110
3111     //Select only hadrons in pt range
3112     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
3113
3114     //remove trigger itself for correlation when use charged triggers
3115     if( track->GetID() == aodParticle->GetTrackLabel(0) || track->GetID() == aodParticle->GetTrackLabel(1) ||
3116         track->GetID() == aodParticle->GetTrackLabel(2) || track->GetID() == aodParticle->GetTrackLabel(3)   )
3117       continue ;
3118
3119     //Only for mixed event frame
3120     Int_t evtIndex2 = 0 ;
3121     if (GetMixedEvent())
3122     {
3123       evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;
3124       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events
3125         continue ;
3126       //vertex cut
3127       if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut())
3128         continue;
3129     }
3130
3131      if(GetDebug() > 2 )
3132       printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
3133
3134     // ------------------------------
3135     // Track type bin or bits setting
3136     //
3137
3138     //
3139     // * Set the pt associated bin for the defined bins *
3140     //
3141     Int_t assocBin   = -1;
3142     for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
3143     {
3144       if(pt > fAssocPtBinLimit[i] && pt < fAssocPtBinLimit[i+1]) assocBin= i;
3145     }
3146
3147     //
3148     // * Assign to the histogram array a bin corresponding
3149     // to a combination of pTa and vz bins *
3150     //
3151     Int_t nz = 1;
3152     Int_t vz = 0;
3153
3154     if(fCorrelVzBin)
3155     {
3156       nz = GetNZvertBin();
3157       vz = GetEventVzBin();
3158     }
3159
3160     Int_t bin = assocBin*nz+vz;
3161
3162     //printf("assoc Bin = %d, vZ bin  = %d, bin = %d \n", assocBin,GetEventVzBin(),bin);
3163
3164     //
3165     // * Get the status of the TOF bit *
3166     //
3167     ULong_t status = track->GetStatus();
3168     Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) ;
3169     //Double32_t tof = track->GetTOFsignal()*1e-3;
3170     Int_t trackBC = track->GetTOFBunchCrossing(bz);
3171
3172     Int_t outTOF = -1;
3173     if     (okTOF && trackBC!=0) outTOF = 1;
3174     else if(okTOF && trackBC==0) outTOF = 0;
3175
3176     // Track multiplicity or cent bin
3177     Int_t cenbin = GetEventCentralityBin(); // combine with vz assoc bin???
3178
3179     //----------------
3180     // Fill Histograms
3181
3182     //
3183     // Azimuthal Angle histograms
3184     //
3185     // Calculate deltaPhi shift
3186     deltaPhi = phiTrig-phi;
3187     if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3188     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3189
3190     FillChargedAngularCorrelationHistograms(pt,  ptTrig,  bin, phi, phiTrig,  deltaPhi,
3191                                             eta, etaTrig, decay, track->GetHMPIDsignal(),
3192                                             outTOF, cenbin, mcTag);
3193     //
3194     // Imbalance zT/xE/pOut histograms
3195     //
3196
3197     //
3198     // Delta phi cut for momentum imbalance correlation
3199     //
3200     if  ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   )
3201       FillChargedMomentumImbalanceHistograms(ptTrig, pt, deltaPhi, cenbin, track->Charge(),
3202                                              bin, decay, outTOF, mcTag);
3203
3204     //
3205     // Underlying event, right side, default case
3206     //
3207     if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
3208       FillChargedUnderlyingEventHistograms(ptTrig, pt, deltaPhi, cenbin, outTOF);
3209
3210     //
3211     // Several UE calculation,  in different perpendicular regions, up to 6:
3212     // left, right, upper-left, lower left, upper-right, lower-right
3213     //
3214     if(fMakeSeveralUE)
3215       FillChargedUnderlyingEventSidesHistograms(ptTrig,pt,deltaPhi);
3216
3217     //
3218     if(fPi0Trigger && decayFound)
3219       FillDecayPhotonCorrelationHistograms(pt, phi, decayMom1,decayMom2, kTRUE) ;
3220
3221     //
3222     // Add track reference to array
3223     //
3224     if(fFillAODWithReferences)
3225     {
3226       nrefs++;
3227       if(nrefs==1)
3228       {
3229         reftracks = new TObjArray(0);
3230         TString trackname = Form("%sTracks", GetAODObjArrayName().Data());
3231         reftracks->SetName(trackname.Data());
3232         reftracks->SetOwner(kFALSE);
3233       }
3234
3235       reftracks->Add(track);
3236     }// reference track to AOD
3237   }// track loop
3238
3239   //Fill AOD with reference tracks, if not filling histograms
3240   if(fFillAODWithReferences && reftracks)
3241   {
3242     aodParticle->AddObjArray(reftracks);
3243   }
3244
3245 }
3246
3247 //_________________________________________________________________________________________________________
3248 void AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation(AliAODPWG4ParticleCorrelation *aodParticle)
3249 {
3250   // Mix current trigger with tracks in another MB event
3251
3252   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Make trigger particle - charged hadron mixed event correlation \n");
3253
3254   if(GetMixedEvent()) return;  // This is not the mixed event from general mixing frame
3255
3256   // Get the event with similar caracteristics
3257   //printf("MakeChargedMixCorrelation for %s\n",GetInputAODName().Data());
3258
3259   AliAnalysisManager   * manager      = AliAnalysisManager::GetAnalysisManager();
3260
3261   AliInputEventHandler * inputHandler = dynamic_cast<AliInputEventHandler*>(manager->GetInputEventHandler());
3262
3263   if(!inputHandler) return;
3264
3265   if(!(inputHandler->IsEventSelected( ) & GetReader()->GetEventTriggerMask())) return;
3266
3267   // Get the pool, check if it exits
3268   Int_t eventBin = GetEventMixBin();
3269
3270   fhEventBin->Fill(eventBin);
3271
3272   //Check that the bin exists, if not (bad determination of RP, ntrality or vz bin) do nothing
3273   if(eventBin < 0) return;
3274
3275   TList * pool     = 0;
3276   TList * poolCalo = 0;
3277   if(fUseMixStoredInReader)
3278   {
3279     pool     = GetReader()->GetListWithMixedEventsForTracks(eventBin);
3280     if(OnlyIsolated() || fFillNeutralEventMixPool) poolCalo = GetReader()->GetListWithMixedEventsForCalo  (eventBin);
3281   }
3282   else
3283   {
3284     pool     = fListMixTrackEvents[eventBin];
3285     if(OnlyIsolated()  || fFillNeutralEventMixPool) poolCalo = fListMixCaloEvents [eventBin];
3286   }
3287
3288   if(!pool) return ;
3289
3290   if((OnlyIsolated()  || fFillNeutralEventMixPool ) && !poolCalo &&
3291      (GetIsolationCut()->GetParticleTypeInCone()!=AliIsolationCut::AliIsolationCut::kOnlyCharged))
3292     printf("AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation() - Careful, cluster pool not available\n");
3293
3294   Double_t ptTrig  = aodParticle->Pt();
3295   Double_t etaTrig = aodParticle->Eta();
3296   Double_t phiTrig = aodParticle->Phi();
3297   if(phiTrig < 0.) phiTrig+=TMath::TwoPi();
3298
3299   if(GetDebug() > 1)
3300     printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Pool bin %d size %d, trigger trigger pt=%f, phi=%f, eta=%f\n",
3301            eventBin,pool->GetSize(), ptTrig,phiTrig,etaTrig);
3302
3303   Double_t ptAssoc  = -999.;
3304   Double_t phiAssoc = -999.;
3305   Double_t etaAssoc = -999.;
3306   Double_t deltaPhi = -999.;
3307   Double_t deltaEta = -999.;
3308   Double_t xE = -999.;
3309   Double_t hbpXE = -999.;
3310
3311   //Start from first event in pool except if in this same event the pool was filled
3312   Int_t ev0 = 0;
3313   if(GetReader()->GetLastTracksMixedEvent() == GetEventNumber()) ev0 = 1;
3314
3315   for(Int_t ev=ev0; ev < pool->GetSize(); ev++)
3316   {
3317     TObjArray* bgTracks = static_cast<TObjArray*>(pool->At(ev));
3318     TObjArray* bgCalo   = 0;
3319
3320     // Check if the particle is isolated in the mixed event, it not, do not fill the histograms
3321     if((OnlyIsolated() || fFillNeutralEventMixPool) && poolCalo)
3322     {
3323       if(pool->GetSize()!=poolCalo->GetSize())
3324         printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Different size of calo and track pools\n");
3325
3326       bgCalo = static_cast<TObjArray*>(poolCalo->At(ev));
3327
3328       if(!bgCalo)
3329         printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation() - Event %d in calo pool not available?\n",ev);
3330
3331       if(OnlyIsolated())
3332       {
3333         Int_t n=0; Int_t nfrac = 0; Bool_t isolated = kFALSE; Float_t coneptsum = 0;
3334         GetIsolationCut()->MakeIsolationCut(bgTracks,bgCalo,
3335                                             GetReader(), GetCaloPID(),
3336                                             kFALSE, aodParticle, "",
3337                                             n,nfrac,coneptsum, isolated);
3338
3339         //printf("AliAnaParticleHadronCorrelation::MakeChargedMixCorrelation() - Isolated? %d - cone %f, ptthres %f",
3340         //       isolated,GetIsolationCut()->GetConeSize(),GetIsolationCut()->GetPtThreshold());
3341         //if(bgTracks)printf(" - n track %d", bgTracks->GetEntriesFast());
3342         //printf("\n");
3343
3344         if(!isolated) continue ;
3345       }
3346     }
3347
3348     fhEventMixBin->Fill(eventBin);
3349
3350     Int_t nTracks=bgTracks->GetEntriesFast();
3351     //printf("\t Read Pool event %d, nTracks %d\n",ev,nTracks);
3352
3353     //Check if it is leading if mixed event
3354     if(fMakeNearSideLeading || fMakeAbsoluteLeading)
3355     {
3356       Bool_t leading = kTRUE;
3357       for(Int_t jlead = 0;jlead <nTracks; jlead++ )
3358       {
3359         AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(jlead) ;
3360
3361         ptAssoc  = track->Pt();
3362         phiAssoc = track->Phi() ;
3363
3364         if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
3365         if (fMakeNearSideLeading)
3366         {
3367           if(ptAssoc > ptTrig && TMath::Abs(phiAssoc-phiTrig) < TMath::PiOver2())
3368           {
3369             leading = kFALSE;
3370             break;
3371           }
3372         }
3373         //jump out this event if there is any other particle with pt larger than trigger
3374         else if(fMakeAbsoluteLeading)
3375         {
3376           if(ptAssoc > ptTrig)
3377           {
3378             leading = kFALSE;
3379             break;
3380           }
3381         }
3382       }
3383
3384       if(fFillNeutralEventMixPool && bgCalo)
3385       {
3386         Int_t nClusters=bgCalo->GetEntriesFast();
3387         TLorentzVector mom ;
3388         for(Int_t jlead = 0;jlead <nClusters; jlead++ )
3389         {
3390           AliAODPWG4Particle *cluster= (AliAODPWG4Particle*) bgCalo->At(jlead) ;
3391
3392           ptAssoc  = cluster->Pt();
3393           phiAssoc = cluster->Phi() ;
3394
3395           if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
3396           if (fMakeNearSideLeading)
3397           {
3398             if(ptAssoc > ptTrig && TMath::Abs(phiAssoc-phiTrig) < TMath::PiOver2())
3399             {
3400               leading = kFALSE;
3401               break;
3402             }
3403           }
3404           //jump out this event if there is any other particle with pt larger than trigger
3405           else if(fMakeAbsoluteLeading)
3406           {
3407             if(ptAssoc > ptTrig)
3408             {
3409               leading = kFALSE;
3410               break;
3411             }
3412           }
3413         }
3414       }
3415
3416       if(!leading) continue; // not leading, check the next event in pool
3417
3418     }
3419
3420     fhPtTriggerMixed   ->Fill(ptTrig);
3421     fhPhiTriggerMixed  ->Fill(ptTrig, phiTrig);
3422     fhEtaTriggerMixed  ->Fill(ptTrig, etaTrig);
3423     fhPtTriggerMixedBin->Fill(ptTrig,eventBin);
3424     if(fCorrelVzBin)fhPtTriggerMixedVzBin->Fill(ptTrig, GetEventVzBin());
3425
3426     for(Int_t j1 = 0;j1 <nTracks; j1++ )
3427     {
3428       AliAODPWG4Particle *track = (AliAODPWG4Particle*) bgTracks->At(j1) ;
3429
3430       if(!track) continue;
3431
3432       ptAssoc  = track->Pt();
3433       etaAssoc = track->Eta();
3434       phiAssoc = track->Phi() ;
3435       if(phiAssoc < 0) phiAssoc+=TMath::TwoPi();
3436
3437       if(IsFiducialCutOn())
3438       {
3439         Bool_t in = GetFiducialCut()->IsInFiducialCut(*aodParticle->Momentum(),"CTS") ;
3440         if(!in) continue ;
3441       }
3442
3443       deltaPhi = phiTrig-phiAssoc;
3444       if(deltaPhi < -TMath::PiOver2())  deltaPhi+=TMath::TwoPi();
3445       if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3446       deltaEta = etaTrig-etaAssoc;
3447
3448       if(GetDebug()>0)
3449         printf("AliAnaParticleHadronCorrelationNew::MakeChargedMixCorrelation(): deltaPhi= %f, deltaEta=%f\n",deltaPhi, deltaEta);
3450
3451       // Set the pt associated bin for the defined bins
3452       Int_t assocBin   = -1;
3453       for(Int_t i = 0 ; i < fNAssocPtBins ; i++)
3454       {
3455         if(ptAssoc > fAssocPtBinLimit[i] && ptAssoc < fAssocPtBinLimit[i+1]) assocBin= i;
3456       }
3457
3458       // Assign to the histogram array a bin corresponding to a combination of pTa and vz bins
3459       Int_t nz = 1;
3460       Int_t vz = 0;
3461
3462       if(fCorrelVzBin)
3463       {
3464         nz = GetNZvertBin();
3465         vz = GetEventVzBin();
3466       }
3467
3468       Int_t bin = assocBin*nz+vz;
3469
3470       fhMixDeltaPhiCharged        ->Fill(ptTrig,  deltaPhi);
3471       fhMixDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
3472
3473       fhMixDeltaPhiCharged        ->Fill(ptTrig,  deltaPhi);
3474       fhMixDeltaPhiDeltaEtaCharged->Fill(deltaPhi, deltaEta);
3475
3476       xE   =-ptAssoc/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
3477       //if(xE <0.)xE =-xE;
3478       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
3479       else        hbpXE =-100;
3480
3481       if ( (deltaPhi > fDeltaPhiMinCut)   && (deltaPhi < fDeltaPhiMaxCut)   )
3482       {
3483         fhMixXECharged->Fill(ptTrig,xE);
3484         fhMixHbpXECharged->Fill(ptTrig,hbpXE);
3485       }
3486
3487       if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
3488       {
3489         //Underlying event region
3490         Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
3491         Double_t uexE = -(ptAssoc/ptTrig)*TMath::Cos(randomphi);
3492
3493         if(uexE < 0.) uexE = -uexE;
3494
3495         fhMixXEUeCharged->Fill(ptTrig,uexE);
3496       }
3497
3498       if(bin < 0) continue ; // this pt bin was not considered
3499
3500       if(TMath::Abs(deltaEta) > 0.8)
3501         fhMixDeltaPhiChargedAssocPtBinDEta08  [bin]->Fill(ptTrig,   deltaPhi);
3502       if(TMath::Abs(deltaEta) < 0.01)
3503         fhMixDeltaPhiChargedAssocPtBinDEta0   [bin]->Fill(ptTrig,   deltaPhi);
3504
3505         fhMixDeltaPhiChargedAssocPtBin        [bin]->Fill(ptTrig,   deltaPhi);
3506         fhMixDeltaPhiDeltaEtaChargedAssocPtBin[bin]->Fill(deltaPhi, deltaEta);
3507
3508     } // track loop
3509   } // mixed event loop
3510 }
3511
3512
3513 //_______________________________________________________________________________________________________
3514 void AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * aodParticle)
3515 {
3516   // Neutral Pion Correlation Analysis
3517
3518   TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); // For the future, foresee more possible pi0 lists
3519   if(!pi0list) return ;
3520
3521   Int_t npi0 = pi0list->GetEntriesFast();
3522   if(npi0 == 0) return ;
3523
3524   if(GetDebug() > 1)
3525     printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Particle - pi0 correlation, %d pi0's\n",npi0);
3526
3527   Int_t evtIndex11 = 0 ;
3528   Int_t evtIndex12 = 0 ;
3529   if (GetMixedEvent())
3530   {
3531     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;
3532     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;
3533   }
3534
3535   Float_t pt   = -100. ;
3536   Float_t zT   = -100. ;
3537   Float_t phi  = -100. ;
3538   Float_t eta  = -100. ;
3539   Float_t xE   = -100. ;
3540   Float_t hbpXE= -100. ;
3541   Float_t hbpZT= -100. ;
3542
3543   Float_t ptTrig  = aodParticle->Pt();
3544   Float_t phiTrig = aodParticle->Phi();
3545   Float_t etaTrig = aodParticle->Eta();
3546   Float_t deltaPhi= -100. ;
3547   Float_t deltaEta= -100. ;
3548
3549   TLorentzVector photonMom ;
3550
3551   // In case of pi0/eta trigger, we may want to check their decay correlation,
3552   // get their decay children
3553   TLorentzVector decayMom1;
3554   TLorentzVector decayMom2;
3555   Bool_t decayFound = kFALSE;
3556   if(fPi0Trigger) decayFound = GetDecayPhotonMomentum(aodParticle,decayMom1, decayMom2);
3557
3558   TObjArray * refpi0 = 0x0;
3559   Int_t nrefs        = 0;
3560
3561   //Loop on stored AOD pi0
3562
3563   for(Int_t iaod = 0; iaod < npi0 ; iaod++)
3564   {
3565     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (pi0list->At(iaod));
3566
3567     Int_t evtIndex2 = 0 ;
3568     Int_t evtIndex3 = 0 ;
3569     if (GetMixedEvent())
3570     {
3571       evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;
3572       evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;
3573
3574       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 ||
3575           evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events
3576         continue ;
3577     }
3578
3579     pt  = pi0->Pt();
3580
3581     if(pt < fMinAssocPt || pt > fMaxAssocPt) continue ;
3582
3583     //remove trigger itself for correlation when use charged triggers
3584     if(aodParticle->GetCaloLabel(0) >= 0 &&
3585        (pi0->GetCaloLabel(0) == aodParticle->GetCaloLabel(0) || pi0->GetCaloLabel(1) == aodParticle->GetCaloLabel(0))) continue ;
3586
3587     if( aodParticle->GetCaloLabel(1) >= 0 &&
3588        (pi0->GetCaloLabel(0) == aodParticle->GetCaloLabel(1) || pi0->GetCaloLabel(1) == aodParticle->GetCaloLabel(1))) continue ;
3589
3590     //
3591     // Angular correlations
3592     //
3593     phi      = pi0->Phi() ;
3594     eta      = pi0->Eta() ;
3595     deltaEta = etaTrig-eta;
3596     deltaPhi = phiTrig-phi;
3597     if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
3598     if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
3599
3600     fhEtaNeutral     ->Fill(pt    , eta     );
3601     fhPhiNeutral     ->Fill(pt    , phi     );
3602     fhDeltaEtaNeutral->Fill(ptTrig, deltaEta);
3603     fhDeltaPhiNeutral->Fill(ptTrig, deltaPhi);
3604
3605     if(pt > 2 ) fhDeltaPhiDeltaEtaNeutral->Fill(deltaPhi, deltaEta);
3606
3607     //
3608     // Momentum imbalance
3609     //
3610     zT  = pt/ptTrig ;
3611
3612     hbpZT = -100;
3613     hbpXE = -100;
3614
3615     if(zT > 0 ) hbpZT = TMath::Log(1./zT);
3616
3617     //delta phi cut for correlation
3618     if( (deltaPhi > fDeltaPhiMinCut) && ( deltaPhi < fDeltaPhiMaxCut) )
3619     {
3620       xE  =-pt/ptTrig*TMath::Cos(deltaPhi); // -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);
3621       //if(xE <0.)xE =-xE;
3622       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
3623
3624       fhDeltaPhiNeutralPt->Fill(pt,deltaPhi);
3625       fhXENeutral        ->Fill(ptTrig,xE);
3626       fhPtHbpXENeutral   ->Fill(ptTrig,hbpXE);
3627       fhZTNeutral        ->Fill(ptTrig,zT);
3628       fhPtHbpZTNeutral   ->Fill(ptTrig,hbpZT);
3629     }
3630     else if ( (deltaPhi > fUeDeltaPhiMinCut) && (deltaPhi < fUeDeltaPhiMaxCut) )
3631     {
3632       // Randomize angle for xE calculation
3633       Double_t randomphi = gRandom->Uniform(fDeltaPhiMinCut,fDeltaPhiMaxCut);
3634
3635       xE = -(pt/ptTrig)*TMath::Cos(randomphi);
3636       if(xE > 0 ) hbpXE = TMath::Log(1./xE);
3637
3638       fhDeltaPhiUeNeutralPt->Fill(pt,deltaPhi);
3639       fhZTUeNeutral        ->Fill(ptTrig,zT);
3640       fhPtHbpZTUeNeutral   ->Fill(ptTrig,hbpZT);
3641       fhXEUeNeutral        ->Fill(ptTrig,xE);
3642       fhPtHbpXEUeNeutral   ->Fill(ptTrig,hbpXE);
3643     }
3644
3645     // Several UE calculation, not sure it is useful
3646     // with partical calorimter acceptance
3647     if(fMakeSeveralUE) FillNeutralUnderlyingEventSidesHistograms(ptTrig,pt,zT,hbpZT,deltaPhi);
3648
3649     //
3650     // Decay photon correlations
3651     //
3652     if(fPi0Trigger && decayFound)
3653       FillDecayPhotonCorrelationHistograms(pt, phi, decayMom1,decayMom2,kFALSE) ;
3654
3655     if(fFillAODWithReferences)
3656     {
3657       nrefs++;
3658       if(nrefs==1)
3659       {
3660         refpi0 = new TObjArray(0);
3661         refpi0->SetName(GetAODObjArrayName()+"Pi0s");
3662         refpi0->SetOwner(kFALSE);
3663       }
3664       refpi0->Add(pi0);
3665     }//put references in trigger AOD
3666
3667     if(GetDebug() > 2 )
3668       printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected pi0: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);
3669
3670   }//loop
3671
3672   //Fill AOD with reference tracks, if not filling histograms
3673   if(fFillAODWithReferences && refpi0)
3674   {
3675     aodParticle->AddObjArray(refpi0);
3676   }
3677 }
3678
3679 //__________________________________________________________________________
3680 void  AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(Int_t label)
3681 {
3682   // Charged Hadron Correlation Analysis with MC information
3683
3684   if ( GetDebug() > 1 )
3685     AliInfo("Make trigger particle - charged hadron correlation in AOD MC level");
3686
3687   if( label < 0 )
3688   {
3689     if( GetDebug() > 0 ) AliInfo(Form(" *** bad label ***:  label %d", label));
3690     return;
3691   }
3692
3693   AliStack         * stack        = 0x0 ;
3694   TParticle        * primary      = 0x0 ;
3695   TClonesArray     * mcparticles  = 0x0 ;
3696   AliAODMCParticle * aodprimary   = 0x0 ;
3697
3698   Double_t eprim   = 0 ;
3699   Double_t ptprim  = 0 ;
3700   Double_t phiprim = 0 ;
3701   Double_t etaprim = 0 ;
3702   Int_t    nTracks = 0 ;
3703   Int_t iParticle  = 0 ;
3704
3705   Bool_t leadTrig = kTRUE;
3706
3707   if( GetReader()->ReadStack() )
3708   {
3709     stack =  GetMCStack() ;
3710     if(!stack)
3711     {
3712       AliFatal("Stack not available, is the MC handler called? STOP");
3713       return;
3714     }
3715
3716     //nTracks = stack->GetNtrack() ;
3717     nTracks = stack->GetNprimary();
3718     if( label >=  stack->GetNtrack() )
3719     {
3720       if(GetDebug() > 2)
3721         AliInfo(Form("*** large label ***:  label %d, n tracks %d", label, stack->GetNtrack()));
3722       return ;
3723     }
3724
3725     primary = stack->Particle(label);
3726     if ( !primary )
3727     {
3728       AliInfo(Form(" *** no primary ***:  label %d", label));
3729       return;
3730     }
3731
3732     eprim    = primary->Energy();
3733     ptprim   = primary->Pt();
3734     etaprim  = primary->Eta();
3735     phiprim  = primary->Phi();
3736     if(phiprim < 0) phiprim+=TMath::TwoPi();
3737
3738     if(ptprim < 0.01 || eprim < 0.01) return ;
3739
3740     for (iParticle = 0 ; iParticle <  nTracks ; iParticle++)
3741     {
3742       TParticle * particle = stack->Particle(iParticle);
3743       TLorentzVector momentum;
3744
3745       //keep only final state particles
3746       if( particle->GetStatusCode() != 1 ) continue ;
3747
3748       //---------- Charged particles ----------------------
3749       Int_t pdg    = particle->GetPdgCode();
3750       Int_t charge = (Int_t) TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
3751       if(charge == 0) continue;
3752
3753       particle->Momentum(momentum);
3754
3755       //Particles in CTS acceptance, make sure to use the same selection as in the reader
3756       Bool_t inCTS =  GetReader()->GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
3757       //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,momentum.Pt(),momentum.Eta(),momentum.Phi()*TMath::RadToDeg());
3758       if( !inCTS ) continue;
3759
3760       // Remove conversions
3761       if ( TMath::Abs(pdg) == 11 && stack->Particle(particle->GetFirstMother())->GetPdgCode() == 22 ) continue ;
3762
3763       if ( label == iParticle ) continue; // avoid trigger particle
3764
3765       Float_t phi = particle->Phi();
3766       if(phi < 0) phi+=TMath::TwoPi();
3767
3768       Bool_t lead = FillChargedMCCorrelationHistograms(particle->Pt(),phi,particle->Eta(),ptprim,phiprim,etaprim);
3769       if(!lead) leadTrig = kFALSE;
3770       //if ( !lead && (fMakeAbsoluteLeading || fMakeNearSideLeading) ) return;
3771
3772     } //track loop
3773
3774   } //ESD MC
3775
3776   else if( GetReader()->ReadAODMCParticles() )
3777   {
3778     //Get the list of MC particles
3779     mcparticles = GetReader()->GetAODMCParticles();
3780     if( !mcparticles ) return;
3781
3782     nTracks = mcparticles->GetEntriesFast() ;
3783
3784     if( label >= nTracks )
3785     {
3786       if(GetDebug() > 2)
3787         AliInfo(Form(" *** large label ***:  label %d, n tracks %d", label,nTracks));
3788       return;
3789     }
3790
3791     //Get the particle
3792     aodprimary = (AliAODMCParticle*) mcparticles->At(label);
3793     if( !aodprimary )
3794     {
3795       AliInfo(Form(" *** no AOD primary ***:  label %d", label));
3796       return;
3797     }
3798
3799     eprim   = aodprimary->E();
3800     ptprim  = aodprimary->Pt();
3801     etaprim = aodprimary->Eta();
3802     phiprim = aodprimary->Phi();
3803     if(phiprim < 0) phiprim+=TMath::TwoPi();
3804
3805     if(ptprim < 0.01 || eprim < 0.01) return ;
3806
3807     for (iParticle = 0; iParticle < nTracks; iParticle++)
3808     {
3809       AliAODMCParticle *part = (AliAODMCParticle*) mcparticles->At(iParticle);
3810
3811       if (!part->IsPhysicalPrimary() ) continue; // same as part->GetStatus() !=1
3812
3813       if ( part->Charge() == 0 ) continue;
3814
3815       TLorentzVector momentum(part->Px(),part->Py(),part->Pz(),part->E());
3816
3817       //Particles in CTS acceptance, make sure to use the same selection as in the reader
3818       Bool_t inCTS =  GetReader()->GetFiducialCut()->IsInFiducialCut(momentum,"CTS");
3819       //printf("Accepted? %d; pt %2.2f, eta %2.2f, phi %2.2f\n",inCTS,momentum.Pt(),momentum.Eta(),momentum.Phi()*TMath::RadToDeg());
3820       if( !inCTS ) continue;
3821
3822       // Remove conversions
3823       Int_t indexmother = part->GetMother();
3824       if ( indexmother > -1 )
3825       {
3826         Int_t pdg = part->GetPdgCode();
3827         Int_t mPdg = ((AliAODMCParticle*) mcparticles->At(indexmother)) ->GetPdgCode();
3828         if (TMath::Abs(pdg) == 11 && mPdg == 22) continue;
3829       }
3830
3831       if ( label == iParticle ) continue; // avoid trigger particle
3832
3833       Float_t phi = part->Phi();
3834       if(phi < 0) phi+=TMath::TwoPi();
3835
3836       Bool_t lead = FillChargedMCCorrelationHistograms(part->Pt(),phi,part->Eta(),ptprim,phiprim,etaprim);
3837       if(!lead) leadTrig = kFALSE;
3838       //if ( !lead && (fMakeAbsoluteLeading || fMakeNearSideLeading)) return;
3839
3840     }  //MC particle loop
3841   }// AOD MC
3842
3843   // Trigger MC particle histograms
3844   //if (!lead  && (fMakeAbsoluteLeading || fMakeNearSideLeading)) return;
3845
3846   fhMCPtTrigger ->Fill(ptprim);
3847   fhMCPhiTrigger->Fill(ptprim,phiprim);
3848   fhMCEtaTrigger->Fill(ptprim,etaprim);
3849
3850   if(!leadTrig && (fMakeAbsoluteLeading || fMakeNearSideLeading) )
3851   {
3852     if(GetDebug() > 1)
3853       printf("AliAnaParticleHadronCorrelation::MakeMCChargedCorrelation(): Not leading primary trigger: pT %2.2f, phi %2.2f, eta %2.2f\n",
3854              ptprim,phiprim*TMath::RadToDeg(),etaprim);
3855
3856     fhMCPtTriggerNotLeading ->Fill(ptprim);
3857     fhMCPhiTriggerNotLeading->Fill(ptprim,phiprim);
3858     fhMCEtaTriggerNotLeading->Fill(ptprim,etaprim);
3859   }
3860 }
3861
3862 //_____________________________________________________________________
3863 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const
3864 {
3865
3866   //Print some relevant parameters set for the analysis
3867   if(! opt)
3868     return;
3869
3870   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
3871   AliAnaCaloTrackCorrBaseClass::Print(" ");
3872   printf("Pt trigger > %2.2f; < %2.2f\n", GetMinPt() , GetMaxPt()) ;
3873   printf("Pt associa > %2.2f; < %2.2f\n", fMinAssocPt, fMaxAssocPt) ;
3874   printf("Phi trigger particle-Hadron   <  %3.2f\n", fDeltaPhiMaxCut) ;
3875   printf("Phi trigger particle-Hadron   >  %3.2f\n", fDeltaPhiMinCut) ;
3876   printf("Phi trigger particle-UeHadron <  %3.2f\n", fUeDeltaPhiMaxCut) ;
3877   printf("Phi trigger particle-UeHadron >  %3.2f\n", fUeDeltaPhiMinCut) ;
3878   printf("Isolated Trigger?  %d\n"     , fSelectIsolated) ;
3879   printf("Several UE?  %d\n"           , fMakeSeveralUE) ;
3880   printf("Name of AOD Pi0 Branch %s \n", fPi0AODBranchName.Data());
3881   printf("Do Decay-hadron correlation ?  %d\n", fPi0Trigger) ;
3882   printf("Select absolute leading for cluster triggers ?  %d or Near Side %d\n",
3883          fMakeAbsoluteLeading, fMakeNearSideLeading) ;
3884   printf("Trigger pt bins  %d\n", fNAssocPtBins) ;
3885   for (Int_t ibin = 0; ibin<fNAssocPtBins; ibin++) {
3886     printf("\t bin %d = [%2.1f,%2.1f]\n", ibin, fAssocPtBinLimit[ibin], fAssocPtBinLimit[ibin+1]) ;
3887   }
3888
3889 }
3890
3891 //____________________________________________________________
3892 void AliAnaParticleHadronCorrelation::SetNAssocPtBins(Int_t n)
3893 {
3894   // Set number of bins
3895
3896   fNAssocPtBins  = n ;
3897
3898   if(n < 20 && n > 0)
3899   {
3900     fNAssocPtBins  = n ;
3901   }
3902   else
3903   {
3904     printf("n = larger than 19 or too small, set to 19 \n");
3905     fNAssocPtBins = 19;
3906   }
3907 }
3908
3909 //______________________________________________________________________________
3910 void AliAnaParticleHadronCorrelation::SetAssocPtBinLimit(Int_t ibin, Float_t pt)
3911 {
3912   // Set the list of limits for the trigger pt bins
3913
3914   if(ibin <= fNAssocPtBins || ibin >= 0)
3915   {
3916     fAssocPtBinLimit[ibin] = pt  ;
3917   }
3918   else
3919   {
3920     printf("AliAnaParticleHadronCorrelation::SetAssocPtBinLimit() - bin  number too large %d > %d or small, nothing done\n", ibin, fNAssocPtBins) ;
3921   }
3922 }
3923