PWGGA/CaloTrackCorrelations/AliAnaParticleHadronCorrelation.cxx

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