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