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