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