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