PWG4/PartCorrDep/AliAnaParticleHadronCorrelation.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 /* $Id: $ */\r
  16 \r
  17 //_________________________________________________________________________\r
  18 // Class for the analysis of particle - hadron correlations\r
  19 // Particle (for example direct gamma) must be found in a previous analysis \r
  20 //-- Author: Gustavo Conesa (LNF-INFN) \r
  21 \r
  22 //  Modified by Yaxian Mao:\r
  23 // 1. add the UE subtraction for corrlation study\r
  24 // 2. change the correlation variable\r
  25 // 3. Only use leading particle(cluster/track) as trigger for correlation (2010/07/02)\r
  26 // 4. Make decay photon-hadron correlations where decay contribute pi0 mass (2010/09/09)\r
  27 // 5. fill the pout to extract kt at the end, also to study charge asymmetry(2010/10/06) \r
  28 //////////////////////////////////////////////////////////////////////////////\r
  29 \r
  30 \r
  31 // --- ROOT system ---\r
  32 //#include "TClonesArray.h"\r
  33 #include "TClass.h"\r
  34 #include "TMath.h"\r
  35 #include "TH3D.h"\r
  36 \r
  37 //---- ANALYSIS system ----\r
  38 #include "AliNeutralMesonSelection.h" \r
  39 #include "AliAnaParticleHadronCorrelation.h" \r
  40 #include "AliCaloTrackReader.h"\r
  41 #include "AliCaloPID.h"\r
  42 #include "AliAODPWG4ParticleCorrelation.h"\r
  43 #include "AliFiducialCut.h"\r
  44 #include "AliVTrack.h"\r
  45 #include "AliVCluster.h"\r
  46 #include "AliMCAnalysisUtils.h"\r
  47 #include "TParticle.h"\r
  48 #include "AliStack.h"\r
  49 #include "AliAODMCParticle.h"\r
  50 #include "AliMixedEvent.h"\r
  51 \r
  52 ClassImp(AliAnaParticleHadronCorrelation)\r
  53 \r
  54 \r
  55 //____________________________________________________________________________\r
  56   AliAnaParticleHadronCorrelation::AliAnaParticleHadronCorrelation(): \r
  57     AliAnaPartCorrBaseClass(),\r
  58     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.), fSelectIsolated(0),\r
  59     fMakeSeveralUE(0),  \r
  60     fUeDeltaPhiMaxCut(0.), fUeDeltaPhiMinCut(0.), \r
  61     fPi0AODBranchName(""),fNeutralCorr(0), fPi0Trigger(0),\r
  62    // fMultiBin(0),fNZvertBin(0),fNrpBin(0),fZvtxCut(0.),\r
  63    // fUseSelectEvent(kFALSE), \r
  64    // fhNclustersNtracks(0), //fhVertex(0),\r
  65     fhPtLeading(0),fhPhiLeading(0),fhEtaLeading(0), \r
  66     fhDeltaPhiDeltaEtaCharged(0),\r
  67     fhPhiCharged(0), fhEtaCharged(0), \r
  68     fhDeltaPhiCharged(0), \r
  69     fhDeltaEtaCharged(0), \r
  70     fhDeltaPhiChargedPt(0), \r
  71     fhDeltaPhiUeChargedPt(0), \r
  72     fhPtImbalanceCharged(0), \r
  73     fhPtImbalanceUeCharged(0),\r
  74     fhPtImbalancePosCharged(0),fhPtImbalanceNegCharged(0),\r
  75     fhPtHbpCharged(0), fhPtHbpUeCharged(0),\r
  76     fhDeltaPhiUeLeftCharged(0),fhDeltaPhiUeRightCharged(0),\r
  77     fhPtImbalanceUeLeftCharged(0),fhPtImbalanceUeRightCharged(0),\r
  78     fhPtHbpUeLeftCharged(0),fhPtHbpUeRightCharged(0), \r
  79     fhPtTrigPout(0), fhPtAssocDeltaPhi(0),\r
  80     fhPtTrigCharged(0),\r
  81     fhTrigDeltaPhiCharged(0x0), fhTrigDeltaEtaCharged(0x0),fhTrigCorr(0x0),fhTrigUeCorr(0x0),\r
  82     fhDeltaPhiDeltaEtaNeutral(0), \r
  83     fhPhiNeutral(0), fhEtaNeutral(0), \r
  84     fhDeltaPhiNeutral(0), fhDeltaEtaNeutral(0),\r
  85     fhDeltaPhiNeutralPt(0),fhDeltaPhiUeNeutralPt(0), \r
  86     fhPtImbalanceNeutral(0),fhPtImbalanceUeNeutral(0),\r
  87     fhPtHbpNeutral(0), fhPtHbpUeNeutral(0),\r
  88     fhDeltaPhiUeLeftNeutral(0),fhDeltaPhiUeRightNeutral(0),\r
  89     fhPtImbalanceUeLeftNeutral(0),fhPtImbalanceUeRightNeutral(0),\r
  90     fhPtHbpUeLeftNeutral(0),fhPtHbpUeRightNeutral(0),\r
  91     fhPtPi0DecayRatio(0),\r
  92     fhDeltaPhiDecayCharged(0),\r
  93     fhPtImbalanceDecayCharged(0), \r
  94     fhDeltaPhiDecayNeutral(0),\r
  95     fhPtImbalanceDecayNeutral(0)\r
  96 {\r
  97   //Default Ctor\r
  98   \r
  99   //Initialize parameters\r
 100   InitParameters();\r
 101 }\r
 102 \r
 103 //________________________________________________________________________\r
 104 TList *  AliAnaParticleHadronCorrelation::GetCreateOutputObjects()\r
 105 {  \r
 106   \r
 107   // Create histograms to be saved in output file and \r
 108   // store them in fOutputContainer\r
 109   TList * outputContainer = new TList() ; \r
 110   outputContainer->SetName("CorrelationHistos") ; \r
 111   \r
 112   Int_t nptbins  = GetHistoPtBins();\r
 113   Int_t nphibins = GetHistoPhiBins();\r
 114   Int_t netabins = GetHistoEtaBins();\r
 115   Float_t ptmax  = GetHistoPtMax();\r
 116   Float_t phimax = GetHistoPhiMax();\r
 117   Float_t etamax = GetHistoEtaMax();\r
 118   Float_t ptmin  = GetHistoPtMin();\r
 119   Float_t phimin = GetHistoPhiMin();\r
 120   Float_t etamin = GetHistoEtaMin();    \r
 121 \r
 122   \r
 123 //  fhNclustersNtracks  = new TH2F ("Multiplicity","Neutral cluster and charged track multiplicity",1000, 0., 1000.,1000, 0., 1000.); \r
 124 //  fhNclustersNtracks->SetYTitle("# of tracks");\r
 125 //  fhNclustersNtracks->SetXTitle("# of clusters");\r
 126 \r
 127   fhPtLeading  = new TH1F ("hPtLeading","p_T distribution of leading particles", nptbins,ptmin,ptmax); \r
 128   fhPtLeading->SetXTitle("p_{T}^{trig} (GeV/c)");\r
 129   \r
 130   fhPhiLeading  = new TH2F ("hPhiLeading","#phi distribution of leading Particles",nptbins,ptmin,ptmax, nphibins,phimin,phimax); \r
 131   fhPhiLeading->SetYTitle("#phi (rad)");\r
 132   \r
 133   fhEtaLeading  = new TH2F ("hEtaLeading","#eta distribution of leading",nptbins,ptmin,ptmax, netabins,etamin,etamax); \r
 134   fhEtaLeading->SetYTitle("#eta ");  \r
 135  // outputContainer->Add(fhNclustersNtracks);\r
 136   outputContainer->Add(fhPtLeading);\r
 137   outputContainer->Add(fhPhiLeading);\r
 138   outputContainer->Add(fhEtaLeading);\r
 139   \r
 140   //Correlation with charged hadrons\r
 141   if(GetReader()->IsCTSSwitchedOn()) {\r
 142     fhDeltaPhiDeltaEtaCharged  = new TH2F\r
 143     ("DeltaPhiDeltaEtaCharged","#phi_{trigger} - #phi_{h^{#pm}} vs #eta_{trigger} - #eta_{h^{#pm}}",\r
 144      140,-2.,5.,200,-2,2); \r
 145     fhDeltaPhiDeltaEtaCharged->SetXTitle("#Delta #phi");\r
 146     fhDeltaPhiDeltaEtaCharged->SetYTitle("#Delta #eta");    \r
 147     \r
 148     fhPhiCharged  = new TH2F\r
 149     ("PhiCharged","#phi_{h^{#pm}}  vs p_{T #pm}",\r
 150      nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
 151     fhPhiCharged->SetYTitle("#phi_{h^{#pm}} (rad)");\r
 152     fhPhiCharged->SetXTitle("p_{T #pm} (GeV/c)");\r
 153     \r
 154     fhEtaCharged  = new TH2F\r
 155     ("EtaCharged","#eta_{h^{#pm}}  vs p_{T #pm}",\r
 156      nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
 157     fhEtaCharged->SetYTitle("#eta_{h^{#pm}} (rad)");\r
 158     fhEtaCharged->SetXTitle("p_{T #pm} (GeV/c)");\r
 159     \r
 160     fhDeltaPhiCharged  = new TH2F\r
 161     ("DeltaPhiCharged","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",\r
 162      nptbins,ptmin,ptmax,140,-2.,5.); \r
 163     fhDeltaPhiCharged->SetYTitle("#Delta #phi");\r
 164     fhDeltaPhiCharged->SetXTitle("p_{T trigger} (GeV/c)");\r
 165     \r
 166     fhDeltaPhiChargedPt  = new TH2F\r
 167     ("DeltaPhiChargedPt","#phi_{trigger} - #phi_{#h^{#pm}} vs p_{T h^{#pm}}",\r
 168      nptbins,ptmin,ptmax,140,-2.,5.);\r
 169     fhDeltaPhiChargedPt->SetYTitle("#Delta #phi");\r
 170     fhDeltaPhiChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");\r
 171     \r
 172     fhDeltaPhiUeChargedPt  = new TH2F\r
 173     ("DeltaPhiUeChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}}",\r
 174      nptbins,ptmin,ptmax,140,-2.,5.);\r
 175     fhDeltaPhiUeChargedPt->SetYTitle("#Delta #phi");\r
 176     fhDeltaPhiUeChargedPt->SetXTitle("p_{T h^{#pm}} (GeV/c)");\r
 177     \r
 178     fhDeltaEtaCharged  = new TH2F\r
 179     ("DeltaEtaCharged","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",\r
 180      nptbins,ptmin,ptmax,200,-2,2); \r
 181     fhDeltaEtaCharged->SetYTitle("#Delta #eta");\r
 182     fhDeltaEtaCharged->SetXTitle("p_{T trigger} (GeV/c)");\r
 183     \r
 184     fhPtImbalanceCharged  = \r
 185     new TH2F("CorrelationCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T trigger}",\r
 186              nptbins,ptmin,ptmax,200,0.,2.); \r
 187     fhPtImbalanceCharged->SetYTitle("z_{trigger h^{#pm}}");\r
 188     fhPtImbalanceCharged->SetXTitle("p_{T trigger}");\r
 189     \r
 190     fhPtImbalanceUeCharged  = \r
 191     new TH2F("CorrelationUeCharged","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger}",\r
 192              nptbins,ptmin,ptmax,200,0.,2.); \r
 193     fhPtImbalanceUeCharged->SetYTitle("z_{trigger Ueh^{#pm}}");\r
 194     fhPtImbalanceUeCharged->SetXTitle("p_{T trigger}");\r
 195     \r
 196     fhPtImbalancePosCharged  = \r
 197     new TH2F("CorrelationPositiveCharged","z_{trigger h^{+}} = p_{T h^{+}} / p_{T trigger}",\r
 198              nptbins,ptmin,ptmax,200,0.,2.); \r
 199     fhPtImbalancePosCharged->SetYTitle("z_{trigger h^{+}}");\r
 200     fhPtImbalancePosCharged->SetXTitle("p_{T trigger}");\r
 201     \r
 202     fhPtImbalanceNegCharged  = \r
 203     new TH2F("CorrelationNegativeCharged","z_{trigger h^{-}} = p_{T h^{-}} / p_{T trigger}",\r
 204              nptbins,ptmin,ptmax,200,0.,2.); \r
 205     fhPtImbalanceNegCharged->SetYTitle("z_{trigger h^{-}}");\r
 206     fhPtImbalanceNegCharged->SetXTitle("p_{T trigger}");\r
 207     \r
 208     fhPtHbpCharged  = \r
 209     new TH2F("HbpCharged","#xi = ln(1/x_{E}) with charged hadrons",\r
 210              nptbins,ptmin,ptmax,200,0.,10.); \r
 211     fhPtHbpCharged->SetYTitle("ln(1/x_{E})");\r
 212     fhPtHbpCharged->SetXTitle("p_{T trigger}");\r
 213     \r
 214     fhPtHbpUeCharged  = \r
 215     new TH2F("HbpUeCharged","#xi = ln(1/x_{E}) with charged hadrons",\r
 216              nptbins,ptmin,ptmax,200,0.,10.); \r
 217     fhPtHbpUeCharged->SetYTitle("ln(1/x_{E})");\r
 218     fhPtHbpUeCharged->SetXTitle("p_{T trigger}");\r
 219     \r
 220     fhPtTrigPout  = \r
 221     new TH2F("PtTrigPout","Pout with triggers",\r
 222              nptbins,ptmin,ptmax,2*nptbins,-ptmax,ptmax); \r
 223     fhPtTrigPout->SetYTitle("p_{out} (GeV/c)");\r
 224     fhPtTrigPout->SetXTitle("p_{T trigger} (GeV/c)"); \r
 225     \r
 226     fhPtAssocDeltaPhi  = \r
 227     new TH2F("fhPtAssocDeltaPhi"," charged hadrons vs. delta phi",\r
 228              nptbins,ptmin,ptmax,140,-2.,5.); \r
 229     fhPtAssocDeltaPhi->SetXTitle("p_{T h^{#pm}} (GeV/c)");  \r
 230     fhPtAssocDeltaPhi->SetYTitle("#Delta #phi (GeV/c)"); \r
 231     \r
 232     fhPtTrigCharged  = \r
 233     new TH2F("PtTrigCharged","trgger and charged tracks pt distribution",\r
 234              nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); \r
 235     fhPtTrigCharged->SetYTitle("p_{T h^{#pm}} (GeV/c)");\r
 236     fhPtTrigCharged->SetXTitle("p_{T trigger} (GeV/c)");    \r
 237           \r
 238     outputContainer->Add(fhDeltaPhiDeltaEtaCharged);\r
 239     outputContainer->Add(fhPhiCharged) ;\r
 240     outputContainer->Add(fhEtaCharged) ;\r
 241     outputContainer->Add(fhDeltaPhiCharged) ; \r
 242     outputContainer->Add(fhDeltaEtaCharged) ;\r
 243     outputContainer->Add(fhDeltaPhiChargedPt) ;\r
 244     outputContainer->Add(fhDeltaPhiUeChargedPt) ;\r
 245     outputContainer->Add(fhPtImbalanceCharged) ;\r
 246     outputContainer->Add(fhPtImbalancePosCharged) ;\r
 247     outputContainer->Add(fhPtImbalanceNegCharged) ;\r
 248     outputContainer->Add(fhPtImbalanceUeCharged) ;\r
 249     outputContainer->Add(fhPtHbpCharged) ;\r
 250     outputContainer->Add(fhPtHbpUeCharged) ;\r
 251     outputContainer->Add(fhPtTrigPout) ;\r
 252     outputContainer->Add(fhPtAssocDeltaPhi) ;\r
 253     outputContainer->Add(fhPtTrigCharged) ;\r
 254     \r
 255     if(DoEventSelect()){ \r
 256       Int_t nMultiBins = GetMultiBin();\r
 257       fhTrigDeltaPhiCharged = new TH2F*[nMultiBins] ;\r
 258       fhTrigDeltaEtaCharged = new TH2F*[nMultiBins] ;\r
 259       fhTrigCorr = new TH2F*[nMultiBins];\r
 260       fhTrigUeCorr = new TH2F*[nMultiBins];\r
 261       for(Int_t im=0; im<nMultiBins; im++){\r
 262         fhTrigDeltaPhiCharged[im]  = new TH2F \r
 263         (Form("fhTrigDeltaPhiCharged_%d",im),Form("fhTrigDeltaPhiCharged_%d",im), nptbins,ptmin,ptmax, 140,-2.,5.); \r
 264         fhTrigDeltaPhiCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");\r
 265         fhTrigDeltaPhiCharged[im]->SetYTitle("#Delta #phi");\r
 266         fhTrigDeltaEtaCharged[im]  = new TH2F \r
 267         (Form("fhTrigDeltaEtaCharged_%d",im),Form("fhTrigDeltaEtaCharged_%d",im), nptbins,ptmin,ptmax, 200,-2,2); \r
 268         fhTrigDeltaEtaCharged[im]->SetXTitle("p_{T trigger} (GeV/c)");\r
 269         fhTrigDeltaEtaCharged[im]->SetYTitle("#Delta #eta");\r
 270         fhTrigCorr[im]  = new TH2F\r
 271         (Form("fhTrigPtCorr_%d",im),Form("fhTrigPtCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.); \r
 272         fhTrigCorr[im]->SetYTitle("z_{trigger h^{#pm}}");\r
 273         fhTrigCorr[im]->SetXTitle("p_{T trigger}");\r
 274         fhTrigUeCorr[im]  = new TH2F\r
 275         (Form("fhTrigPtUeCorr_%d",im),Form("fhTrigPtUeCorr_%d",im), nptbins,ptmin,ptmax,200,0.,2.); \r
 276         fhTrigUeCorr[im]->SetYTitle("z_{trigger h^{#pm}}");\r
 277         fhTrigUeCorr[im]->SetXTitle("p_{T trigger}");       \r
 278 \r
 279         outputContainer->Add(fhTrigDeltaPhiCharged[im]) ;\r
 280         outputContainer->Add(fhTrigDeltaEtaCharged[im]) ;\r
 281         outputContainer->Add(fhTrigCorr[im]);\r
 282         outputContainer->Add(fhTrigUeCorr[im]);\r
 283         \r
 284         }\r
 285     }\r
 286    \r
 287     if(fPi0Trigger){\r
 288       fhPtPi0DecayRatio  = new TH2F\r
 289       ("hPtPi0DecayRatio","p_T of #pi^{0} and the ratio of pt for two decay", \r
 290        nptbins,ptmin,ptmax, 100,0.,2.); \r
 291       fhPtPi0DecayRatio->SetXTitle("p_{T}^{#pi^{0}} (GeV/c)");\r
 292       fhPtPi0DecayRatio->SetYTitle("p_{T}^{Decay}/p_{T}^{#pi^{0}}");\r
 293       \r
 294       fhDeltaPhiDecayCharged  = new TH2F\r
 295       ("DeltaPhiDecayCharged","#phi_{Decay} - #phi_{h^{#pm}} vs p_{T Decay}",\r
 296        nptbins,ptmin,ptmax,140,-2.,5.); \r
 297       fhDeltaPhiDecayCharged->SetYTitle("#Delta #phi");\r
 298       fhDeltaPhiDecayCharged->SetXTitle("p_{T Decay} (GeV/c)");\r
 299       \r
 300       fhPtImbalanceDecayCharged  = \r
 301       new TH2F("CorrelationDecayCharged","z_{trigger h^{#pm}} = p_{T h^{#pm}} / p_{T Decay}",\r
 302                nptbins,ptmin,ptmax,200,0.,2.); \r
 303       fhPtImbalanceDecayCharged->SetYTitle("z_{decay h^{#pm}}");\r
 304       fhPtImbalanceDecayCharged->SetXTitle("p_{T decay}");\r
 305       \r
 306       outputContainer->Add(fhPtPi0DecayRatio) ; \r
 307       outputContainer->Add(fhDeltaPhiDecayCharged) ; \r
 308       outputContainer->Add(fhPtImbalanceDecayCharged) ;\r
 309     }    \r
 310     \r
 311     \r
 312     if(fMakeSeveralUE){ \r
 313       fhDeltaPhiUeLeftCharged  = new TH2F\r
 314       ("DeltaPhiUeLeftChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE left side range of trigger particles",\r
 315        nptbins,ptmin,ptmax,140,-2.,5.);\r
 316       fhDeltaPhiUeLeftCharged->SetYTitle("#Delta #phi");\r
 317       fhDeltaPhiUeLeftCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");\r
 318       outputContainer->Add(fhDeltaPhiUeLeftCharged) ;\r
 319       \r
 320       fhDeltaPhiUeRightCharged  = new TH2F\r
 321       ("DeltaPhiUeRightChargedPt","#phi_{trigger} - #phi_{#Ueh^{#pm}} vs p_{T Ueh^{#pm}} with UE right side range of trigger particles",\r
 322        nptbins,ptmin,ptmax,140,-2.,5.);\r
 323       fhDeltaPhiUeRightCharged->SetYTitle("#Delta #phi");\r
 324       fhDeltaPhiUeRightCharged->SetXTitle("p_{T h^{#pm}} (GeV/c)");\r
 325       outputContainer->Add(fhDeltaPhiUeRightCharged) ;\r
 326       \r
 327       fhPtImbalanceUeLeftCharged  = \r
 328       new TH2F("CorrelationUeChargedLeft","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE left side of trigger",\r
 329                nptbins,ptmin,ptmax,200,0.,2.); \r
 330       fhPtImbalanceUeLeftCharged->SetYTitle("z_{trigger Ueh^{#pm}}");\r
 331       fhPtImbalanceUeLeftCharged->SetXTitle("p_{T trigger}");\r
 332       outputContainer->Add(fhPtImbalanceUeLeftCharged) ;\r
 333       \r
 334       fhPtImbalanceUeRightCharged  = \r
 335       new TH2F("CorrelationUeChargedRight","z_{trigger h^{#pm}} = p_{T Ueh^{#pm}} / p_{T trigger} with UE right side of trigger",\r
 336                nptbins,ptmin,ptmax,200,0.,2.); \r
 337       fhPtImbalanceUeRightCharged->SetYTitle("z_{trigger Ueh^{#pm}}");\r
 338       fhPtImbalanceUeRightCharged->SetXTitle("p_{T trigger}");\r
 339       outputContainer->Add(fhPtImbalanceUeRightCharged) ;\r
 340       \r
 341       fhPtHbpUeLeftCharged  = \r
 342       new TH2F("HbpUeChargedLeft","#xi = ln(1/x_{E}) with charged UE left side of trigger",\r
 343                nptbins,ptmin,ptmax,200,0.,10.); \r
 344       fhPtHbpUeLeftCharged->SetYTitle("ln(1/x_{E})");\r
 345       fhPtHbpUeLeftCharged->SetXTitle("p_{T trigger}");\r
 346       outputContainer->Add(fhPtHbpUeLeftCharged) ;\r
 347       \r
 348       fhPtHbpUeRightCharged  = \r
 349       new TH2F("HbpUeChargedRight","#xi = ln(1/x_{E}) with charged UE right side of trigger",\r
 350                nptbins,ptmin,ptmax,200,0.,10.); \r
 351       fhPtHbpUeRightCharged->SetYTitle("ln(1/x_{E})");\r
 352       fhPtHbpUeRightCharged->SetXTitle("p_{T trigger}");\r
 353       outputContainer->Add(fhPtHbpUeRightCharged) ;\r
 354       \r
 355     }  \r
 356   }  //Correlation with charged hadrons\r
 357   \r
 358   //Correlation with neutral hadrons\r
 359   if(fNeutralCorr){\r
 360     \r
 361     fhDeltaPhiDeltaEtaNeutral  = new TH2F\r
 362     ("DeltaPhiDeltaEtaNeutral","#phi_{trigger} - #phi_{h^{0}} vs #eta_{trigger} - #eta_{h^{0}}",\r
 363      140,-2.,5.,200,-2,2); \r
 364     fhDeltaPhiDeltaEtaNeutral->SetXTitle("#Delta #phi");\r
 365     fhDeltaPhiDeltaEtaNeutral->SetYTitle("#Delta #eta");   \r
 366           \r
 367     fhPhiNeutral  = new TH2F\r
 368     ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T #pi^{0}}",\r
 369      nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
 370     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");\r
 371     fhPhiNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
 372     \r
 373     fhEtaNeutral  = new TH2F\r
 374     ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T #pi^{0}}",\r
 375      nptbins,ptmin,ptmax,netabins,etamin,etamax); \r
 376     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");\r
 377     fhEtaNeutral->SetXTitle("p_{T #pi^{0}} (GeV/c)");\r
 378     \r
 379     fhDeltaPhiNeutral  = new TH2F\r
 380     ("DeltaPhiNeutral","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",\r
 381      nptbins,ptmin,ptmax,nphibins,phimin,phimax); \r
 382     fhDeltaPhiNeutral->SetYTitle("#Delta #phi");\r
 383     fhDeltaPhiNeutral->SetXTitle("p_{T trigger} (GeV/c)");\r
 384     \r
 385     fhDeltaPhiNeutralPt  = new TH2F\r
 386     ("DeltaPhiNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",\r
 387      nptbins,ptmin,ptmax,140,-2.,5.); \r
 388     fhDeltaPhiNeutralPt->SetYTitle("#Delta #phi");\r
 389     fhDeltaPhiNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");\r
 390     \r
 391     fhDeltaPhiUeNeutralPt  = new TH2F\r
 392     ("DeltaPhiUeNeutralPt","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T #pi^{0}}}",\r
 393      nptbins,ptmin,ptmax,140,-2.,5.); \r
 394     fhDeltaPhiUeNeutralPt->SetYTitle("#Delta #phi");\r
 395     fhDeltaPhiUeNeutralPt->SetXTitle("p_{T h^{0}} (GeV/c)");\r
 396     \r
 397     fhDeltaEtaNeutral  = new TH2F\r
 398     ("DeltaEtaNeutral","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",\r
 399      nptbins,ptmin,ptmax,200,-2,2); \r
 400     fhDeltaEtaNeutral->SetYTitle("#Delta #eta");\r
 401     fhDeltaEtaNeutral->SetXTitle("p_{T trigger} (GeV/c)");\r
 402     \r
 403     fhPtImbalanceNeutral  = \r
 404     new TH2F("CorrelationNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",\r
 405              nptbins,ptmin,ptmax,200,0.,2.); \r
 406     fhPtImbalanceNeutral->SetYTitle("z_{trigger #pi^{0}}");\r
 407     fhPtImbalanceNeutral->SetXTitle("p_{T trigger}");\r
 408     \r
 409     fhPtImbalanceUeNeutral  = \r
 410     new TH2F("CorrelationUeNeutral","z_{trigger #pi} = p_{T #pi^{0}} / p_{T trigger}",\r
 411              nptbins,ptmin,ptmax,200,0.,2.); \r
 412     fhPtImbalanceUeNeutral->SetYTitle("z_{trigger #pi^{0}}");\r
 413     fhPtImbalanceUeNeutral->SetXTitle("p_{T trigger}");\r
 414     \r
 415     fhPtHbpNeutral  = \r
 416     new TH2F("HbpNeutral","#xi = ln(1/x_{E}) with neutral particles",\r
 417              nptbins,ptmin,ptmax,200,0.,10.); \r
 418     fhPtHbpNeutral->SetYTitle("ln(1/x_{E})");\r
 419     fhPtHbpNeutral->SetXTitle("p_{T trigger}");\r
 420     \r
 421     fhPtHbpUeNeutral  = \r
 422     new TH2F("HbpUeNeutral","#xi = ln(1/x_{E}) with neutral particles",\r
 423              nptbins,ptmin,ptmax,200,0.,10.); \r
 424     fhPtHbpUeNeutral->SetYTitle("ln(1/x_{E})");\r
 425     fhPtHbpUeNeutral->SetXTitle("p_{T trigger}");\r
 426     \r
 427     \r
 428     outputContainer->Add(fhDeltaPhiDeltaEtaNeutral); \r
 429     outputContainer->Add(fhPhiNeutral) ;  \r
 430     outputContainer->Add(fhEtaNeutral) ;   \r
 431     outputContainer->Add(fhDeltaPhiNeutral) ; \r
 432     outputContainer->Add(fhDeltaPhiNeutralPt) ; \r
 433     outputContainer->Add(fhDeltaPhiUeNeutralPt) ; \r
 434     outputContainer->Add(fhDeltaEtaNeutral) ; \r
 435     outputContainer->Add(fhPtImbalanceNeutral) ;\r
 436     outputContainer->Add(fhPtImbalanceUeNeutral) ;  \r
 437     outputContainer->Add(fhPtHbpNeutral) ;\r
 438     outputContainer->Add(fhPtHbpUeNeutral) ;    \r
 439     \r
 440     if(fPi0Trigger){\r
 441       fhDeltaPhiDecayNeutral  = new TH2F\r
 442       ("DeltaPhiDecayNeutral","#phi_{Decay} - #phi_{h^{0}} vs p_{T Decay}",\r
 443        nptbins,ptmin,ptmax,140,-2.,5.); \r
 444       fhDeltaPhiDecayNeutral->SetYTitle("#Delta #phi");\r
 445       fhDeltaPhiDecayNeutral->SetXTitle("p_{T Decay} (GeV/c)");\r
 446       \r
 447       fhPtImbalanceDecayNeutral  = \r
 448       new TH2F("CorrelationDecayNeutral","z_{trigger h^{0}} = p_{T h^{0}} / p_{T Decay}",\r
 449                nptbins,ptmin,ptmax,200,0.,2.); \r
 450       fhPtImbalanceDecayNeutral->SetYTitle("z_{decay h^{0}}");\r
 451       fhPtImbalanceDecayNeutral->SetXTitle("p_{T decay}");\r
 452       \r
 453       outputContainer->Add(fhDeltaPhiDecayNeutral) ; \r
 454       outputContainer->Add(fhPtImbalanceDecayNeutral) ;\r
 455     }\r
 456     \r
 457     \r
 458     if(fMakeSeveralUE){ \r
 459       fhDeltaPhiUeLeftNeutral  = new TH2F\r
 460       ("DeltaPhiUeLeftNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T h^{0}} with neutral UE left side range of trigger particles",\r
 461        nptbins,ptmin,ptmax,140,-2.,5.);\r
 462       fhDeltaPhiUeLeftNeutral->SetYTitle("#Delta #phi");\r
 463       fhDeltaPhiUeLeftNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");\r
 464       outputContainer->Add(fhDeltaPhiUeLeftNeutral) ;\r
 465       \r
 466       fhDeltaPhiUeRightNeutral  = new TH2F\r
 467       ("DeltaPhiUeRightNeutralPt","#phi_{trigger} - #phi_{#Ueh^{0}} vs p_{T Ueh^{0}} with neutral UE right side range of trigger particles",\r
 468        nptbins,ptmin,ptmax,140,-2.,5.);\r
 469       fhDeltaPhiUeRightNeutral->SetYTitle("#Delta #phi");\r
 470       fhDeltaPhiUeRightNeutral->SetXTitle("p_{T h^{0}} (GeV/c)");\r
 471       outputContainer->Add(fhDeltaPhiUeRightNeutral) ;\r
 472       \r
 473       fhPtImbalanceUeLeftNeutral  = \r
 474       new TH2F("CorrelationUeNeutralLeft","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE left side of trigger",\r
 475                nptbins,ptmin,ptmax,140,0.,2.); \r
 476       fhPtImbalanceUeLeftNeutral->SetYTitle("z_{trigger Ueh^{0}}");\r
 477       fhPtImbalanceUeLeftNeutral->SetXTitle("p_{T trigger}");\r
 478       outputContainer->Add(fhPtImbalanceUeLeftNeutral) ;\r
 479       \r
 480       fhPtImbalanceUeRightNeutral  = \r
 481       new TH2F("CorrelationUeNeutralRight","z_{trigger h^{0}} = p_{T Ueh^{0}} / p_{T trigger} with neutral UE right side of trigger",\r
 482                nptbins,ptmin,ptmax,200,0.,2.); \r
 483       fhPtImbalanceUeRightNeutral->SetYTitle("z_{trigger Ueh^{0}}");\r
 484       fhPtImbalanceUeRightNeutral->SetXTitle("p_{T trigger}");\r
 485       outputContainer->Add(fhPtImbalanceUeRightNeutral) ;\r
 486       \r
 487       fhPtHbpUeLeftNeutral  = \r
 488       new TH2F("HbpUeNeutralLeft","#xi = ln(1/x_{E}) with neutral UE left side of trigger",\r
 489                nptbins,ptmin,ptmax,200,0.,10.); \r
 490       fhPtHbpUeLeftNeutral->SetYTitle("ln(1/x_{E})");\r
 491       fhPtHbpUeLeftNeutral->SetXTitle("p_{T trigger}");\r
 492       outputContainer->Add(fhPtHbpUeLeftNeutral) ;\r
 493       \r
 494       fhPtHbpUeRightNeutral  = \r
 495       new TH2F("HbpUeNeutralRight","#xi = ln(1/x_{E}) with neutral UE right side of trigger",\r
 496                nptbins,ptmin,ptmax,200,0.,10.); \r
 497       fhPtHbpUeRightNeutral->SetYTitle("ln(1/x_{E})");\r
 498       fhPtHbpUeRightNeutral->SetXTitle("p_{T trigger}");\r
 499       outputContainer->Add(fhPtHbpUeRightNeutral) ;\r
 500       \r
 501     }  \r
 502     \r
 503     //Keep neutral meson selection histograms if requiered\r
 504     //Setting done in AliNeutralMesonSelection\r
 505     if(GetNeutralMesonSelection()){\r
 506       TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;\r
 507       if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())\r
 508         for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;\r
 509       delete nmsHistos;\r
 510     }\r
 511     \r
 512   }//Correlation with neutral hadrons\r
 513   \r
 514   return outputContainer;\r
 515   \r
 516 }\r
 517 \r
 518 //____________________________________________________________________________\r
 519 void AliAnaParticleHadronCorrelation::InitParameters()\r
 520 {\r
 521   \r
 522   //Initialize the parameters of the analysis.\r
 523   SetInputAODName("PWG4Particle");\r
 524   SetAODObjArrayName("Hadrons");  \r
 525   AddToHistogramsName("AnaHadronCorr_");\r
 526 \r
 527   SetPtCutRange(0.,300);\r
 528   fDeltaPhiMinCut = 1.5 ;\r
 529   fDeltaPhiMaxCut = 4.5 ;\r
 530   fSelectIsolated = kFALSE;\r
 531   fMakeSeveralUE = kFALSE;\r
 532   fUeDeltaPhiMinCut = 1. ;\r
 533   fUeDeltaPhiMaxCut = 1.5 ;\r
 534   fNeutralCorr = kFALSE ;\r
 535   fPi0Trigger = kFALSE ;\r
 536 \r
 537 }\r
 538 \r
 539 //__________________________________________________________________\r
 540 void AliAnaParticleHadronCorrelation::Print(const Option_t * opt) const\r
 541 {\r
 542 \r
 543   //Print some relevant parameters set for the analysis\r
 544   if(! opt)\r
 545     return;\r
 546   \r
 547   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
 548   AliAnaPartCorrBaseClass::Print(" ");\r
 549 \r
 550   printf("Phi trigger particle-Hadron      <     %3.2f\n", fDeltaPhiMaxCut) ; \r
 551   printf("Phi trigger particle-Hadron      >     %3.2f\n", fDeltaPhiMinCut) ;\r
 552   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;\r
 553   printf("Phi trigger particle-UeHadron    <    %3.2f\n", fUeDeltaPhiMaxCut) ; \r
 554   printf("Phi trigger particle-UeHadron    >    %3.2f\n", fUeDeltaPhiMinCut) ;\r
 555   printf("Several UE?  %d\n", fMakeSeveralUE) ;\r
 556   printf("Name of AOD Pi0 Branch %s \n",fPi0AODBranchName.Data());\r
 557   printf("Do Decay-hadron correlation ?  %d\n", fPi0Trigger) ;\r
 558 \r
 559   \r
 560 } \r
 561 \r
 562 //__________________________________________________________________\r
 563 TObjString* AliAnaParticleHadronCorrelation::GetAnalysisCuts()\r
 564 {\r
 565   //Save parameters used for analysis\r
 566   TString parList ; //this will be list of parameters used for this analysis.\r
 567   const Int_t buffersize = 255;\r
 568   char onePar[buffersize] ;\r
 569  \r
 570   snprintf(onePar,buffersize,"--- AliAnaPaticleHadronCorrelation ---\n") ;\r
 571   parList+=onePar ;     \r
 572   snprintf(onePar,buffersize,"Phi trigger particle-Hadron      <     %3.2f ", fDeltaPhiMaxCut) ; \r
 573   parList+=onePar ;\r
 574   snprintf(onePar,buffersize,"Phi trigger particle-Hadron      >     %3.2f ", fDeltaPhiMinCut) ;\r
 575   parList+=onePar ;\r
 576   snprintf(onePar,buffersize,"Isolated Trigger?  %d\n", fSelectIsolated) ;\r
 577   parList+=onePar ;\r
 578   snprintf(onePar,buffersize,"Phi trigger particle-UeHadron    <    %3.2f ", fUeDeltaPhiMaxCut) ; \r
 579   parList+=onePar ;\r
 580   snprintf(onePar,buffersize,"Phi trigger particle-UeHadron    >    %3.2f ", fUeDeltaPhiMinCut) ;\r
 581   parList+=onePar ;\r
 582   snprintf(onePar,buffersize,"Several UE?  %d\n", fMakeSeveralUE) ;\r
 583   parList+=onePar ;\r
 584   snprintf(onePar,buffersize,"Name of AOD Pi0 Branch %s ",fPi0AODBranchName.Data());\r
 585   parList+=onePar ;\r
 586   snprintf(onePar,buffersize,"Do Decay-hadron correlation ?  %d", fPi0Trigger) ;\r
 587   parList+=onePar ;\r
 588 \r
 589   //Get parameters set in base class.\r
 590   parList += GetBaseParametersList() ;\r
 591   \r
 592   //Get parameters set in PID class.\r
 593   //parList += GetCaloPID()->GetPIDParametersList() ;\r
 594   \r
 595   //Get parameters set in FiducialCut class (not available yet)\r
 596   //parlist += GetFidCut()->GetFidCutParametersList() \r
 597   \r
 598   return new TObjString(parList) ;  \r
 599 \r
 600 } \r
 601 \r
 602 \r
 603 //____________________________________________________________________________\r
 604 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD()  \r
 605 {  \r
 606   //Particle-Hadron Correlation Analysis, fill AODs\r
 607   \r
 608   if(!GetInputAODBranch()){\r
 609     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());\r
 610     abort();\r
 611   }\r
 612         \r
 613   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){\r
 614     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());\r
 615     abort();\r
 616   }\r
 617         \r
 618   if(GetDebug() > 1){\r
 619     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - Begin hadron correlation analysis, fill AODs \n");\r
 620     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());\r
 621     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetCTSTracks()->GetEntriesFast());\r
 622     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());\r
 623     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetPHOSClusters()->GetEntriesFast());\r
 624   }\r
 625   \r
 626   //Loop on stored AOD particles, trigger\r
 627   Double_t ptTrig    = 0.;\r
 628   Int_t    trigIndex = -1;\r
 629   Int_t naod = GetInputAODBranch()->GetEntriesFast();\r
 630   //fhNclustersNtracks->Fill(naod, GetCTSTracks()->GetEntriesFast());\r
 631   for(Int_t iaod = 0; iaod < naod ; iaod++){\r
 632     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));\r
 633     //find the leading particles with highest momentum\r
 634     if (particle->Pt()>ptTrig) {\r
 635       ptTrig = particle->Pt() ;\r
 636       trigIndex = iaod ;\r
 637     }\r
 638   }//Aod branch loop\r
 639         \r
 640   //Do correlation with leading particle\r
 641   if(trigIndex!=-1){\r
 642           \r
 643     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));\r
 644     //Make correlation with charged hadrons\r
 645     if(GetReader()->IsCTSSwitchedOn() )\r
 646       MakeChargedCorrelation(particle, GetCTSTracks(),kFALSE);\r
 647     \r
 648     TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists\r
 649     if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)\r
 650       MakeNeutralCorrelation(particle, pi0list,kFALSE);\r
 651     \r
 652   }//Correlate leading\r
 653   \r
 654   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillAOD() - End fill AODs \n");\r
 655   \r
 656 }\r
 657 \r
 658 //____________________________________________________________________________\r
 659 void  AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms()  \r
 660 {  \r
 661   //Particle-Hadron Correlation Analysis, fill histograms\r
 662   \r
 663   if(!GetInputAODBranch()){\r
 664     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());\r
 665     abort();\r
 666   }\r
 667   \r
 668   if(GetDebug() > 1){\r
 669     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - Begin hadron correlation analysis, fill histograms \n");\r
 670     printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());\r
 671   }\r
 672   \r
 673   //Get the vertex and check it is not too large in z\r
 674   Double_t v[3] = {0,0,0}; //vertex ;\r
 675   GetReader()->GetVertex(v);\r
 676   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;  \r
 677   \r
 678   //Loop on stored AOD particles, find leading\r
 679   Int_t naod = GetInputAODBranch()->GetEntriesFast();\r
 680  // if(naod!=0)fhVertex->Fill(v[0],v[1],v[2]);\r
 681   Double_t ptTrig = 0.;\r
 682   Int_t trigIndex = -1 ;\r
 683   for(Int_t iaod = 0; iaod < naod ; iaod++){     //loop on input trigger AOD file \r
 684     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));\r
 685     //vertex cut in case of mixing\r
 686     if (GetMixedEvent()) {\r
 687       Int_t evt=-1;\r
 688       Int_t id =-1;\r
 689       if     (particle->GetCaloLabel(0)  >= 0 ){\r
 690         id=particle->GetCaloLabel(0); \r
 691         if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;\r
 692       }\r
 693       else if(particle->GetTrackLabel(0) >= 0 ){\r
 694         id=particle->GetTrackLabel(0);\r
 695         if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;\r
 696       }\r
 697       else continue;\r
 698       \r
 699       if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut()) \r
 700         return ;\r
 701     }\r
 702 \r
 703     //check if the particle is isolated or if we want to take the isolation into account\r
 704     if(OnlyIsolated() && !particle->IsIsolated()) continue;\r
 705     //check if inside the vertex cut\r
 706     //find the leading particles with highest momentum\r
 707     if (particle->Pt()>ptTrig) {\r
 708       ptTrig = particle->Pt() ;\r
 709                   trigIndex = iaod ;\r
 710     }\r
 711   }//finish searching for leading trigger particle\r
 712   if(trigIndex!=-1){ //using trigger partilce to do correlations\r
 713     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(trigIndex));\r
 714           \r
 715     if (GetMixedEvent()) {\r
 716       Int_t evt=-1;\r
 717       Int_t id = 0;\r
 718       if     (particle->GetCaloLabel(0)  >= 0 ){\r
 719         id=particle->GetCaloLabel(0); \r
 720         if(id >= 0 )evt= GetMixedEvent()-> EventIndexForCaloCluster(id) ;\r
 721       }\r
 722       else if(particle->GetTrackLabel(0) >= 0 ){\r
 723         id=particle->GetTrackLabel(0);\r
 724         if(id >= 0 )evt= GetMixedEvent()->EventIndex(id) ;\r
 725       }\r
 726       else return;\r
 727       \r
 728       if (TMath::Abs(GetVertex(evt)[2]) > GetZvertexCut()) \r
 729         return ;\r
 730     }\r
 731         \r
 732     //Fill leading particle histogram   \r
 733     fhPtLeading->Fill(particle->Pt());\r
 734     Float_t phi = particle->Phi();\r
 735     if(phi<0)phi+=TMath::TwoPi();\r
 736     fhPhiLeading->Fill(particle->Pt(), phi);\r
 737     fhEtaLeading->Fill(particle->Pt(), particle->Eta());\r
 738     \r
 739     //Make correlation with charged hadrons\r
 740     if(GetReader()->IsCTSSwitchedOn() )\r
 741       MakeChargedCorrelation(particle, GetCTSTracks(),kTRUE);\r
 742     \r
 743     TObjArray * pi0list = (TObjArray*) GetAODBranch(fPi0AODBranchName); //For the future, foresee more possible pi0 lists\r
 744     if(fNeutralCorr && pi0list && pi0list->GetEntriesFast() > 0)\r
 745       MakeNeutralCorrelation(particle, pi0list,kTRUE);\r
 746     \r
 747   }//Aod branch loop\r
 748   \r
 749   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeAnalysisFillHistograms() - End fill histograms \n");\r
 750   \r
 751 }\r
 752 \r
 753 //____________________________________________________________________________\r
 754 void  AliAnaParticleHadronCorrelation::MakeChargedCorrelation(AliAODPWG4ParticleCorrelation *aodParticle, TObjArray* const pl, const Bool_t bFillHisto)\r
 755 {  \r
 756   // Charged Hadron Correlation Analysis\r
 757   if(GetDebug() > 1)printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Make trigger particle - charged hadron correlation \n");\r
 758   \r
 759   Int_t evtIndex11 = -1 ; //cluster trigger or pi0 trigger \r
 760   Int_t evtIndex12 = -1 ; // pi0 trigger\r
 761   Int_t evtIndex13 = -1 ; // charged trigger\r
 762   Int_t indexPhoton1 = -1 ;\r
 763   Int_t indexPhoton2 = -1 ;  \r
 764   Double_t v[3] = {0,0,0}; //vertex ;\r
 765   GetReader()->GetVertex(v);\r
 766   if(!GetMixedEvent() && TMath::Abs(v[2]) > GetZvertexCut()) return ;  \r
 767 \r
 768   Int_t nTracks = GetCTSTracks()->GetEntriesFast() ;\r
 769   \r
 770   if (GetMixedEvent()) {\r
 771     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;\r
 772     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;    \r
 773     evtIndex13 = GetMixedEvent()->EventIndex(aodParticle->GetTrackLabel(0)) ;\r
 774   }\r
 775   \r
 776   Double_t ptTrig  = aodParticle->Pt();\r
 777   Double_t pxTrig  = aodParticle->Px();\r
 778   Double_t pyTrig  = aodParticle->Py();\r
 779   \r
 780   Double_t phiTrig = aodParticle->Phi();\r
 781   Double_t etaTrig = aodParticle->Eta();\r
 782   \r
 783   Double_t pt   = -100.;\r
 784   Double_t px   = -100.;\r
 785   Double_t py   = -100.;\r
 786   Double_t rat  = -100.; \r
 787   Double_t xE   = -100.; \r
 788   Double_t cosi = -100.; \r
 789   Double_t phi  = -100. ;\r
 790   Double_t eta  = -100. ;\r
 791   TVector3 p3;  \r
 792         \r
 793   TObjArray * reftracks = 0x0;\r
 794   Int_t nrefs = 0;\r
 795         \r
 796   Double_t ptDecay1  = 0. ;\r
 797   Double_t pxDecay1  = 0. ;\r
 798   Double_t pyDecay1  = 0. ;\r
 799   Double_t phiDecay1 = 0. ;\r
 800   Double_t ptDecay2  = 0. ;\r
 801   Double_t pxDecay2  = 0. ;\r
 802   Double_t pyDecay2  = 0. ;\r
 803   Double_t phiDecay2 = 0. ;\r
 804   \r
 805   Double_t ratDecay1  = -100.;  \r
 806   Double_t ratDecay2  = -100.; \r
 807   Float_t deltaphi    = -100. ;\r
 808   Float_t deltaphiDecay1 = -100. ;\r
 809   Float_t deltaphiDecay2 = -100. ;\r
 810   TObjArray * clusters   = 0x0 ;  \r
 811   TLorentzVector photonMom ;    \r
 812   \r
 813   if(fPi0Trigger){\r
 814     indexPhoton1 = aodParticle->GetCaloLabel (0);\r
 815     indexPhoton2 = aodParticle->GetCaloLabel (1);\r
 816     if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);\r
 817     \r
 818     if(indexPhoton1!=-1 && indexPhoton2!=-1){\r
 819       if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;\r
 820       else  clusters = GetPHOSClusters() ;\r
 821       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){\r
 822         AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));   \r
 823         photon->GetMomentum(photonMom,GetVertex(0)) ;\r
 824         if(photon->GetID()==indexPhoton1) {\r
 825           ptDecay1  = photonMom.Pt();\r
 826           pxDecay1  = photonMom.Px();\r
 827           pyDecay1  = photonMom.Py();\r
 828           phiDecay1 = photonMom.Phi();\r
 829           if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig);\r
 830         }\r
 831         if(photon->GetID()==indexPhoton2) {\r
 832           ptDecay2  = photonMom.Pt();\r
 833           pxDecay2  = photonMom.Px();\r
 834           pyDecay2  = photonMom.Py();\r
 835           phiDecay2 = photonMom.Phi();\r
 836           if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay2/ptTrig);\r
 837         } \r
 838         if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);\r
 839       } //cluster loop        \r
 840     } //index of decay photons found\r
 841     \r
 842   } //make decay-hadron correlation\r
 843   \r
 844   //Track loop, select tracks with good pt, phi and fill AODs or histograms\r
 845   //Int_t currentIndex = -1 ; \r
 846   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){\r
 847     AliVTrack * track = (AliVTrack *) (pl->At(ipr)) ;\r
 848 \r
 849     //check if inside the vertex cut\r
 850     //printf("charge = %d\n", track->Charge());\r
 851     Int_t evtIndex2 = 0 ; \r
 852     if (GetMixedEvent()) {\r
 853       evtIndex2 = GetMixedEvent()->EventIndex(track->GetID()) ;\r
 854       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex13 == evtIndex2 ) // photon and track from different events\r
 855         continue ; \r
 856       //vertex cut\r
 857       if (TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) \r
 858          return ;\r
 859       //      if(currentIndex == evtIndex2) // tracks from different event \r
 860       //        continue ;\r
 861       //      currentIndex = evtIndex2 ;\r
 862     }\r
 863     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};\r
 864     p3.SetXYZ(mom[0],mom[1],mom[2]);\r
 865     pt   = p3.Pt();\r
 866     px   = p3.Px();\r
 867     py   = p3.Py();\r
 868     eta  = p3.Eta();\r
 869     phi  = p3.Phi() ;\r
 870     if(phi < 0) phi+=TMath::TwoPi();\r
 871 \r
 872     //Select only hadrons in pt range\r
 873     if(pt < GetMinPt() || pt > GetMaxPt()) continue ;\r
 874     //remove trigger itself for correlation when use charged triggers    \r
 875     if(track->GetID()==aodParticle->GetTrackLabel(0) && pt==ptTrig && phi==phiTrig && eta==etaTrig) \r
 876       continue ;    \r
 877     if(IsFiducialCutOn()){\r
 878       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;\r
 879       if(! in ) continue ;\r
 880     }    \r
 881     //jumped out this event if near side associated partile pt larger than trigger\r
 882     if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2())  break ;\r
 883     rat   = pt/ptTrig ;\r
 884     xE    = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);\r
 885     if(xE <0.)xE =-xE;\r
 886     cosi = TMath::Log(1/xE);\r
 887     // printf("rat = %f, xE = %f, cosi =%f \n", rat, xE, cosi);\r
 888     // printf("phi = %f \n", phi);\r
 889     \r
 890      if(fPi0Trigger){\r
 891       if(indexPhoton1!=-1 && indexPhoton2!=-1){\r
 892         if(ptDecay1) ratDecay1 = pt/ptDecay1 ;\r
 893         if(ptDecay2) ratDecay2 = pt/ptDecay2 ; \r
 894         deltaphiDecay1 = phiDecay1-phi;\r
 895         deltaphiDecay2 = phiDecay2-phi;\r
 896         if(deltaphiDecay1< -TMath::PiOver2()) deltaphiDecay1+=TMath::TwoPi();\r
 897         if(deltaphiDecay1>3*TMath::PiOver2()) deltaphiDecay1-=TMath::TwoPi();\r
 898         if(deltaphiDecay2< -TMath::PiOver2()) deltaphiDecay2+=TMath::TwoPi();\r
 899         if(deltaphiDecay2>3*TMath::PiOver2()) deltaphiDecay2-=TMath::TwoPi();    \r
 900       }\r
 901     } //do decay-hadron correlation    \r
 902             \r
 903     //Selection within angular range\r
 904     deltaphi = phiTrig-phi;\r
 905     if(deltaphi< -TMath::PiOver2()) deltaphi+=TMath::TwoPi();\r
 906     if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();\r
 907 \r
 908     Double_t pout = pt*TMath::Sin(deltaphi) ;\r
 909     \r
 910     if(GetDebug() > 2)\r
 911       printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Charged hadron: pt %f, phi %f, phi trigger %f. Cuts:  delta phi  %2.2f < %2.2f < %2.2f, pT min %2.2f \n",\r
 912              pt,phi, phiTrig,fDeltaPhiMinCut, deltaphi, fDeltaPhiMaxCut, GetMinPt());\r
 913     \r
 914     if(bFillHisto){\r
 915       // Fill Histograms\r
 916       fhEtaCharged->Fill(pt,eta);\r
 917       fhPhiCharged->Fill(pt,phi);\r
 918       fhDeltaEtaCharged->Fill(ptTrig,aodParticle->Eta()-eta);\r
 919       fhDeltaPhiCharged->Fill(ptTrig, deltaphi);\r
 920       fhDeltaPhiDeltaEtaCharged->Fill(deltaphi,aodParticle->Eta()-eta);\r
 921       fhPtAssocDeltaPhi->Fill(pt, deltaphi);\r
 922       \r
 923       if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeChargedCorrelation() - Selected charge for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);\r
 924       //fill different multiplicity histogram\r
 925       if(DoEventSelect()){\r
 926         for(Int_t im=0; im<GetMultiBin(); im++){\r
 927           if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1)){\r
 928             fhTrigDeltaPhiCharged[im]->Fill(ptTrig,deltaphi);\r
 929             fhTrigDeltaEtaCharged[im]->Fill(ptTrig,aodParticle->Eta()-eta);\r
 930          }\r
 931         }\r
 932       }\r
 933       //delta phi cut for correlation\r
 934       if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {\r
 935         fhDeltaPhiChargedPt->Fill(pt,deltaphi);\r
 936         fhPtImbalanceCharged->Fill(ptTrig,xE); \r
 937         fhPtHbpCharged->Fill(ptTrig,cosi);\r
 938         fhPtTrigPout->Fill(ptTrig, pout) ;\r
 939         fhPtTrigCharged->Fill(ptTrig, pt) ;\r
 940         if(track->Charge()>0) fhPtImbalancePosCharged->Fill(ptTrig,xE) ;\r
 941         else fhPtImbalanceNegCharged->Fill(ptTrig,xE) ;\r
 942         //fill different multiplicity histogram\r
 943         if(DoEventSelect()){\r
 944           for(Int_t im=0; im<GetMultiBin(); im++){\r
 945             if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))\r
 946               fhTrigCorr[im]->Fill(ptTrig,xE);\r
 947           }\r
 948         } //multiplicity events selection\r
 949       } //delta phi cut for correlation\r
 950       else if ((deltaphi > fUeDeltaPhiMinCut) && ( deltaphi < fUeDeltaPhiMaxCut)) { //UE study\r
 951         fhDeltaPhiUeChargedPt->Fill(pt,deltaphi);\r
 952         //fhUePoutPtTrigPtAssoc->Fill(ptTrig, pt, pout) ;\r
 953         fhPtImbalanceUeCharged->Fill(ptTrig,xE);\r
 954         fhPtHbpUeCharged->Fill(ptTrig,cosi);\r
 955         if(DoEventSelect()){\r
 956           for(Int_t im=0; im<GetMultiBin(); im++){\r
 957             if(nTracks < ( GetMaxMulti() - GetMinMulti() )/GetMultiBin()*(im+1))\r
 958               fhTrigUeCorr[im]->Fill(ptTrig,xE);\r
 959           }\r
 960         } //multiplicity events selection\r
 961         \r
 962       } //UE study\r
 963       \r
 964       if(fPi0Trigger){\r
 965         if(indexPhoton1!=-1 && indexPhoton2!=-1){\r
 966           fhDeltaPhiDecayCharged->Fill(ptDecay1, deltaphiDecay1);\r
 967           fhDeltaPhiDecayCharged->Fill(ptDecay2, deltaphiDecay2);\r
 968           if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);\r
 969           if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )\r
 970             fhPtImbalanceDecayCharged->Fill(ptDecay1,ratDecay1); \r
 971           if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )\r
 972             fhPtImbalanceDecayCharged->Fill(ptDecay2,ratDecay2);\r
 973           if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);\r
 974         } //index of decay photons found\r
 975       } //make decay-hadron correlation          \r
 976       \r
 977       //several UE calculation \r
 978       if(fMakeSeveralUE){\r
 979         if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){  \r
 980           fhDeltaPhiUeLeftCharged->Fill(pt,deltaphi);\r
 981           fhPtImbalanceUeLeftCharged->Fill(ptTrig,rat);\r
 982           fhPtHbpUeLeftCharged->Fill(ptTrig,cosi);\r
 983         }\r
 984         if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){  \r
 985           fhDeltaPhiUeRightCharged->Fill(pt,deltaphi);\r
 986           fhPtImbalanceUeRightCharged->Fill(ptTrig,rat);\r
 987           fhPtHbpUeRightCharged->Fill(ptTrig,cosi);\r
 988           \r
 989         }\r
 990       } //several UE calculation\r
 991       \r
 992     } //Fill histogram \r
 993     else{\r
 994       nrefs++;\r
 995       if(nrefs==1){\r
 996         reftracks = new TObjArray(0);\r
 997         reftracks->SetName(GetAODObjArrayName()+"Tracks");\r
 998         reftracks->SetOwner(kFALSE);\r
 999       }\r
1000       reftracks->Add(track);\r
1001     }//aod particle loop\r
1002   }// track loop\r
1003   \r
1004   //Fill AOD with reference tracks, if not filling histograms\r
1005   if(!bFillHisto && reftracks) {\r
1006     aodParticle->AddObjArray(reftracks);\r
1007   }\r
1008   \r
1009   //delete reftracks;\r
1010   \r
1011 }  \r
1012 \r
1013 \r
1014 \r
1015 //____________________________________________________________________________\r
1016 //void  AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD(AliAODPWG4ParticleCorrelation* const aodParticle,TObjArray* const pl, TString detector)  \r
1017 //{  \r
1018 //    // Neutral Pion Correlation Analysis, find pi0, put them in new output aod, if correlation cuts passed\r
1019 //  if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Make trigger particle - neutral hadron correlation \n");\r
1020 //  \r
1021 //  if(!NewOutputAOD()){\r
1022 //    printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Output aod not created, set AOD class name and branch name in the configuration file, STOP! \n");\r
1023 //    abort();\r
1024 //  }\r
1025 //  \r
1026 //  Double_t phiTrig = aodParticle->Phi();\r
1027 //  Int_t       tag = 0;\r
1028 //  TLorentzVector gammai;\r
1029 //  TLorentzVector gammaj;\r
1030 //  \r
1031 //  //Get vertex for photon momentum calculation\r
1032 //  \r
1033 //  if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) \r
1034 //  {\r
1035 //    for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) {\r
1036 //      if (!GetMixedEvent()) \r
1037 //        GetReader()->GetVertex(GetVertex(iev));\r
1038 //      else \r
1039 //        GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev)); \r
1040 //    }\r
1041 //  }\r
1042 //  Double_t vertex2[] = {0.0,0.0,0.0} ; //vertex of second input aod\r
1043 //  if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) \r
1044 //  {\r
1045 //     if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);\r
1046 //  }\r
1047 //      \r
1048 //    //Cluster loop, select pairs with good pt, phi and fill AODs or histograms\r
1049 //    //Int_t iEvent= GetReader()->GetEventNumber() ;\r
1050 //  Int_t nclus = pl->GetEntriesFast();\r
1051 //  for(Int_t iclus = 0;iclus < nclus ; iclus ++ ){\r
1052 //    AliVCluster * calo = (AliVCluster *) (pl->At(iclus)) ;\r
1053 //    \r
1054 //    Int_t evtIndex1 = 0 ; \r
1055 //    if (GetMixedEvent()) {\r
1056 //      evtIndex1=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; \r
1057 //    }\r
1058 //    \r
1059 //\r
1060 //      //Input from second AOD?\r
1061 //    Int_t inputi = 0;\r
1062 //    if     (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= iclus) \r
1063 //      inputi = 1 ;\r
1064 //    else if(aodParticle->GetDetector() == "PHOS"  && GetReader()->GetPHOSClustersNormalInputEntries()  <= iclus) \r
1065 //      inputi = 1;\r
1066 //    \r
1067 //      //Cluster selection, not charged, with photon or pi0 id and in fiducial cut\r
1068 //    //FIXME\r
1069 //    Int_t pdg=0;\r
1070 //    //if     (inputi == 0 && !SelectCluster(calo, GetVertex(evtIndex1),  gammai, pdg))  \r
1071 //      continue ;\r
1072 //    //MEFIX\r
1073 //    else if(inputi == 1 && !SelectCluster(calo, vertex2, gammai, pdg))  \r
1074 //      continue ;\r
1075 //    \r
1076 //    if(GetDebug() > 2)\r
1077 //      printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral cluster in %s: pt %f, phi %f, phi trigger %f. Cuts:  delta phi min %2.2f,  max %2.2f, pT min %2.2f \n",\r
1078 //             detector.Data(), gammai.Pt(),gammai.Phi(),phiTrig,fDeltaPhiMinCut, fDeltaPhiMaxCut, GetMinPt());\r
1079 //    \r
1080 //      //2 gamma overlapped, found with PID\r
1081 //    if(pdg == AliCaloPID::kPi0){ \r
1082 //      \r
1083 //        //Select only hadrons in pt range\r
1084 //      if(gammai.Pt() < GetMinPt() || gammai.Pt() > GetMaxPt()) \r
1085 //        continue ;\r
1086 //      \r
1087 //        //Selection within angular range\r
1088 //      Float_t phi = gammai.Phi();\r
1089 //      if(phi < 0) phi+=TMath::TwoPi();\r
1090 //        //Float_t deltaphi = TMath::Abs(phiTrig-phi);\r
1091 //        //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;\r
1092 //      \r
1093 //      AliAODPWG4Particle pi0 = AliAODPWG4Particle(gammai);\r
1094 //        //pi0.SetLabel(calo->GetLabel());\r
1095 //      pi0.SetPdg(AliCaloPID::kPi0);\r
1096 //      pi0.SetDetector(detector);\r
1097 //      \r
1098 //      if(IsDataMC()){\r
1099 //        pi0.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(),inputi));\r
1100 //        if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of candidate %d\n",pi0.GetTag());\r
1101 //      }//Work with stack also \r
1102 //       //Set the indeces of the original caloclusters  \r
1103 //      pi0.SetCaloLabel(calo->GetID(),-1);\r
1104 //      AddAODParticle(pi0);\r
1105 //      \r
1106 //      if(GetDebug() > 2) \r
1107 //        printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Correlated with selected pi0 (pid): pt %f, phi %f\n",pi0.Pt(),pi0.Phi());\r
1108 //      \r
1109 //    }// pdg = 111\r
1110 //    \r
1111 //      //Make invariant mass analysis\r
1112 //    else if(pdg == AliCaloPID::kPhoton){      \r
1113 //        //Search the photon companion in case it comes from  a Pi0 decay\r
1114 //        //Apply several cuts to select the good pair;\r
1115 //      for(Int_t jclus = iclus+1; jclus < pl->GetEntries() ; jclus ++ ){\r
1116 //        AliVCluster * calo2 = (AliVCluster *) (pl->At(jclus)) ;\r
1117 //        Int_t evtIndex2 = 0 ; \r
1118 //        if (GetMixedEvent()) {\r
1119 //          evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; \r
1120 //        }\r
1121 //        if (GetMixedEvent() && (evtIndex1 == evtIndex2))\r
1122 //          continue ;\r
1123 //        \r
1124 //          //Input from second AOD?\r
1125 //        Int_t inputj = 0;\r
1126 //        if     (aodParticle->GetDetector() == "EMCAL" && GetReader()->GetEMCALClustersNormalInputEntries() <= jclus) \r
1127 //          inputj = 1;\r
1128 //        else if(aodParticle->GetDetector() == "PHOS"  && GetReader()->GetPHOSClustersNormalInputEntries()  <= jclus) \r
1129 //          inputj = 1;\r
1130 //        \r
1131 //          //Cluster selection, not charged with photon or pi0 id and in fiducial cut\r
1132 //        Int_t pdgj=0;\r
1133 //        //FIXME\r
1134 //        //if     (inputj == 0 && !SelectCluster(calo2, GetVertex(evtIndex2),  gammaj, pdgj))  \r
1135 //          continue ;\r
1136 //        //MEFIX\r
1137 //\r
1138 //        else if(inputj == 1 && !SelectCluster(calo2, vertex2, gammaj, pdgj))  \r
1139 //          continue ;\r
1140 //        //FIXME\r
1141 //        //if(!SelectCluster(calo2,GetVertex(evtIndex2), gammaj, pdgj)) \r
1142 //        //MEFIX\r
1143 //          continue ;\r
1144 //        \r
1145 //        if(pdgj == AliCaloPID::kPhoton ){\r
1146 //          \r
1147 //          if((gammai+gammaj).Pt() < GetMinPt() || (gammai+gammaj).Pt() > GetMaxPt()) \r
1148 //            continue ;\r
1149 //          \r
1150 //            //Selection within angular range\r
1151 //          Float_t phi = (gammai+gammaj).Phi();\r
1152 //          if(phi < 0) phi+=TMath::TwoPi();\r
1153 //            //Float_t deltaphi = TMath::Abs(phiTrig-phi);\r
1154 //            //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;\r
1155 //          \r
1156 //            //Select good pair (aperture and invariant mass)\r
1157 //          if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){\r
1158 //            \r
1159 //            if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Neutral Hadron Correlation: AOD Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f\n",\r
1160 //                                       (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());\r
1161 //            \r
1162 //            TLorentzVector pi0mom = gammai+gammaj;\r
1163 //            AliAODPWG4Particle pi0 = AliAODPWG4Particle(pi0mom);\r
1164 //              //pi0.SetLabel(calo->GetLabel());\r
1165 //            pi0.SetPdg(AliCaloPID::kPi0);\r
1166 //            pi0.SetDetector(detector);        \r
1167 //            if(IsDataMC()){\r
1168 //                //Check origin of the candidates\r
1169 //              \r
1170 //              Int_t label1 = calo->GetLabel();\r
1171 //              Int_t label2 = calo2->GetLabel();\r
1172 //              Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), inputi);\r
1173 //              Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), inputj);\r
1174 //              \r
1175 //              if(GetDebug() > 0) \r
1176 //                printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);\r
1177 //              if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) && GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)){\r
1178 //                \r
1179 //                  //Check if pi0 mother is the same\r
1180 //                if(GetReader()->ReadStack()){ \r
1181 //                  TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree\r
1182 //                  label1 = mother1->GetFirstMother();\r
1183 //                    //mother1 = GetMCStack()->Particle(label1);//pi0\r
1184 //                  \r
1185 //                  TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree\r
1186 //                  label2 = mother2->GetFirstMother();\r
1187 //                    //mother2 = GetMCStack()->Particle(label2);//pi0\r
1188 //                }\r
1189 //                else if(GetReader()->ReadAODMCParticles()){\r
1190 //                  AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputi))->At(label1);//photon in kine tree\r
1191 //                  label1 = mother1->GetMother();\r
1192 //                    //mother1 = GetMCStack()->Particle(label1);//pi0\r
1193 //                  AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(inputj))->At(label2);//photon in kine tree\r
1194 //                  label2 = mother2->GetMother();\r
1195 //                    //mother2 = GetMCStack()->Particle(label2);//pi0\r
1196 //                }\r
1197 //                \r
1198 //                  //printf("mother1 %d, mother2 %d\n",label1,label2);\r
1199 //                if(label1 == label2)\r
1200 //                  GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);\r
1201 //              }\r
1202 //            }//Work with mc information also   \r
1203 //            pi0.SetTag(tag);\r
1204 //              //Set the indeces of the original caloclusters  \r
1205 //            pi0.SetCaloLabel(calo->GetID(), calo2->GetID());\r
1206 //            AddAODParticle(pi0);\r
1207 //            \r
1208 //            \r
1209 //          }//Pair selected\r
1210 //        }//if pair of gammas\r
1211 //      }//2nd loop\r
1212 //    }// if pdg = 22\r
1213 //  }//1st loop\r
1214 //  \r
1215 //  if(GetDebug() > 1) \r
1216 //    printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillAOD() - End, %d pi0's found \n",GetOutputAODBranch()->GetEntriesFast());\r
1217 //}\r
1218 \r
1219 //____________________________________________________________________________\r
1220 void  AliAnaParticleHadronCorrelation::MakeNeutralCorrelation(AliAODPWG4ParticleCorrelation * const aodParticle, TObjArray* pi0list, const Bool_t bFillHisto)  \r
1221 {  \r
1222   // Neutral Pion Correlation Analysis\r
1223   if(GetDebug() > 1) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Make trigger particle - pi0 correlation, %d pi0's \n",pi0list->GetEntriesFast());\r
1224   \r
1225   Int_t evtIndex11 = 0 ; \r
1226   Int_t evtIndex12 = 0 ; \r
1227   if (GetMixedEvent()) {\r
1228     evtIndex11 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(0)) ;\r
1229     evtIndex12 = GetMixedEvent()->EventIndexForCaloCluster(aodParticle->GetCaloLabel(1)) ;    \r
1230   }  \r
1231   \r
1232   Double_t pt   = -100.;\r
1233   Double_t px   = -100.;\r
1234   Double_t py   = -100.;\r
1235   Double_t rat = -100.; \r
1236   Double_t phi = -100.;\r
1237   Double_t eta = -100.;\r
1238   Double_t xE  = -100.; \r
1239   Double_t cosi  = -100.; \r
1240   \r
1241   Double_t ptTrig  = aodParticle->Pt();\r
1242   Double_t phiTrig = aodParticle->Phi();\r
1243   Double_t etaTrig = aodParticle->Eta();\r
1244   Double_t pxTrig  = aodParticle->Px();\r
1245   Double_t pyTrig  = aodParticle->Py();\r
1246   \r
1247   \r
1248   Int_t indexPhoton1 = -1 ;\r
1249   Int_t indexPhoton2 = -1 ;    \r
1250   Double_t ptDecay1 = 0. ;\r
1251   Double_t pxDecay1  = 0. ;\r
1252   Double_t pyDecay1  = 0. ;\r
1253   Double_t phiDecay1  = 0. ;\r
1254   Double_t ptDecay2 = 0. ;\r
1255   Double_t pxDecay2  = 0. ;\r
1256   Double_t pyDecay2  = 0. ;\r
1257   Double_t phiDecay2  = 0. ;\r
1258   \r
1259   Double_t ratDecay1  = -100.;  \r
1260   Double_t ratDecay2  = -100.; \r
1261   Float_t deltaphi = -100. ;\r
1262   Float_t deltaphiDecay1 = -100. ;\r
1263   Float_t deltaphiDecay2 = -100. ;\r
1264   TObjArray * clusters = 0x0 ;  \r
1265   TLorentzVector photonMom ;    \r
1266   if(fPi0Trigger){\r
1267     indexPhoton1 = aodParticle->GetCaloLabel (0);\r
1268     indexPhoton2 = aodParticle->GetCaloLabel (1);\r
1269     if(GetDebug() > 1)printf("indexPhoton1 = %d, indexPhoton2 = %d \n", indexPhoton1, indexPhoton2);\r
1270     \r
1271     if(indexPhoton1!=-1 && indexPhoton2!=-1){\r
1272       if(aodParticle->GetDetector()=="EMCAL") clusters = GetEMCALClusters() ;\r
1273       else                                    clusters = GetPHOSClusters() ;\r
1274       for(Int_t iclus = 0; iclus < clusters->GetEntriesFast(); iclus++){\r
1275         AliVCluster * photon =  (AliVCluster*) (clusters->At(iclus));   \r
1276         photon->GetMomentum(photonMom,GetVertex(0)) ;\r
1277         if(photon->GetID()==indexPhoton1) {\r
1278           ptDecay1  = photonMom.Pt();\r
1279           pxDecay1  = photonMom.Px();\r
1280           pyDecay1  = photonMom.Py();\r
1281           phiDecay1 = photonMom.Phi();\r
1282         }\r
1283         if(photon->GetID()==indexPhoton2) {\r
1284           ptDecay2  = photonMom.Pt();\r
1285           pxDecay2  = photonMom.Px();\r
1286           pyDecay2  = photonMom.Py();\r
1287           phiDecay2 = photonMom.Phi();\r
1288         } \r
1289         if(GetDebug() > 1)printf("Photon1 = %f, Photon2 = %f \n", ptDecay1, ptDecay2);\r
1290       } //photonAOD loop        \r
1291     } //index of decay photons found\r
1292     if(ptTrig && bFillHisto) fhPtPi0DecayRatio->Fill(ptTrig, ptDecay1/ptTrig, ptDecay2/ptTrig);\r
1293   } //make decay-hadron correlation\r
1294   \r
1295   TObjArray * refpi0    =0x0;\r
1296   Int_t nrefs = 0;\r
1297   \r
1298   //Loop on stored AOD pi0\r
1299   Int_t naod = pi0list->GetEntriesFast();\r
1300   if(GetDebug() > 0) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelationFillHistograms() -  aod branch entries %d\n", naod);\r
1301   for(Int_t iaod = 0; iaod < naod ; iaod++){\r
1302     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (pi0list->At(iaod));\r
1303     \r
1304     Int_t evtIndex2 = 0 ; \r
1305     Int_t evtIndex3 = 0 ; \r
1306     if (GetMixedEvent()) {\r
1307       evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(0)) ;\r
1308       evtIndex3 = GetMixedEvent()->EventIndexForCaloCluster(pi0->GetCaloLabel(1)) ;\r
1309       \r
1310       if (evtIndex11 == evtIndex2 || evtIndex12 == evtIndex2 || evtIndex11 == evtIndex3 || evtIndex12 == evtIndex3) // trigger and pi0 are not from different events\r
1311         continue ; \r
1312     }      \r
1313     \r
1314     //Int_t pdg = pi0->GetPdg();\r
1315     //if(pdg != AliCaloPID::kPi0) continue;  \r
1316     \r
1317     pt  = pi0->Pt();\r
1318     px  = pi0->Px();\r
1319     py  = pi0->Py();    \r
1320     if(pt < GetMinPt() || pt > GetMaxPt()) continue ;\r
1321     //jumped out this event if near side associated partile pt larger than trigger\r
1322     if(pt > ptTrig && TMath::Abs(phi-phiTrig)<TMath::PiOver2())  break ;\r
1323 \r
1324     //Selection within angular range\r
1325     phi = pi0->Phi();\r
1326     //Float_t deltaphi = TMath::Abs(phiTrig-phi);\r
1327     //if( (deltaphi < fDeltaPhiMinCut) || ( deltaphi > fDeltaPhiMaxCut) ) continue ;\r
1328     \r
1329     if(bFillHisto){\r
1330       \r
1331       deltaphi = phiTrig-phi;\r
1332       if(deltaphi<-TMath::PiOver2()) deltaphi+=TMath::TwoPi();\r
1333       if(deltaphi>3*TMath::PiOver2()) deltaphi-=TMath::TwoPi();\r
1334       \r
1335       rat = pt/ptTrig ;\r
1336       phi = pi0->Phi() ;\r
1337       eta = pi0->Eta() ;\r
1338       xE   = -(px*pxTrig+py*pyTrig)/(ptTrig*ptTrig);\r
1339       if(xE <0.)xE =-xE;\r
1340       cosi = TMath::Log(1/xE);\r
1341       \r
1342       if(fPi0Trigger){\r
1343         if(indexPhoton1!=-1 && indexPhoton2!=-1){\r
1344           if(ptDecay1) ratDecay1 = pt/ptDecay1 ;\r
1345           if(ptDecay2) ratDecay2 = pt/ptDecay2 ; \r
1346           deltaphiDecay1 = phiDecay1-phi;\r
1347           deltaphiDecay2 = phiDecay2-phi;\r
1348           if(deltaphiDecay1< -TMath::PiOver2()) deltaphiDecay1+=TMath::TwoPi();\r
1349           if(deltaphiDecay1>3*TMath::PiOver2()) deltaphiDecay1-=TMath::TwoPi();\r
1350           if(deltaphiDecay2< -TMath::PiOver2()) deltaphiDecay2+=TMath::TwoPi();\r
1351           if(deltaphiDecay2>3*TMath::PiOver2()) deltaphiDecay2-=TMath::TwoPi();   \r
1352           fhDeltaPhiDecayNeutral->Fill(ptDecay1, deltaphiDecay1);\r
1353           fhDeltaPhiDecayNeutral->Fill(ptDecay2, deltaphiDecay2);\r
1354           if(GetDebug() > 1)printf("deltaPhoton1 = %f, deltaPhoton2 = %f \n", deltaphiDecay1, deltaphiDecay2);\r
1355           if( (deltaphiDecay1 > fDeltaPhiMinCut) && ( deltaphiDecay1 < fDeltaPhiMaxCut) )\r
1356             fhPtImbalanceDecayNeutral->Fill(ptDecay1,ratDecay1); \r
1357           if( (deltaphiDecay2 > fDeltaPhiMinCut) && ( deltaphiDecay2 < fDeltaPhiMaxCut) )\r
1358             fhPtImbalanceDecayNeutral->Fill(ptDecay2,ratDecay2);\r
1359           if(GetDebug() > 1)printf("ratPhoton1 = %f, ratPhoton2 = %f \n", pt/ptDecay1, pt/ptDecay2);\r
1360         }\r
1361       } //do decay-hadron correlation\r
1362       \r
1363       fhEtaNeutral->Fill(pt,eta);\r
1364       fhPhiNeutral->Fill(pt,phi);\r
1365       fhDeltaEtaNeutral->Fill(ptTrig,etaTrig-eta);\r
1366       fhDeltaPhiNeutral->Fill(ptTrig,deltaphi);\r
1367       fhDeltaPhiDeltaEtaNeutral->Fill(deltaphi,etaTrig-eta);\r
1368       \r
1369       //delta phi cut for correlation\r
1370       if( (deltaphi > fDeltaPhiMinCut) && ( deltaphi < fDeltaPhiMaxCut) ) {\r
1371         fhDeltaPhiNeutralPt->Fill(pt,deltaphi);\r
1372         fhPtImbalanceNeutral->Fill(ptTrig,rat); \r
1373         fhPtHbpNeutral->Fill(ptTrig,cosi); \r
1374       }\r
1375       else {\r
1376         fhDeltaPhiUeNeutralPt->Fill(pt,deltaphi);\r
1377         fhPtImbalanceUeNeutral->Fill(ptTrig,rat);\r
1378         fhPtHbpUeNeutral->Fill(ptTrig,cosi); \r
1379       }\r
1380       //several UE calculation \r
1381       if(fMakeSeveralUE){\r
1382         if((deltaphi<-fUeDeltaPhiMinCut) && (deltaphi >-fUeDeltaPhiMaxCut)){  \r
1383           fhDeltaPhiUeLeftNeutral->Fill(pt,deltaphi);\r
1384           fhPtImbalanceUeLeftNeutral->Fill(ptTrig,rat);\r
1385           fhPtHbpUeLeftNeutral->Fill(ptTrig,cosi);\r
1386         }\r
1387         if((deltaphi>fUeDeltaPhiMinCut) && (deltaphi <fUeDeltaPhiMaxCut)){  \r
1388           fhDeltaPhiUeRightNeutral->Fill(pt,deltaphi);\r
1389           fhPtImbalanceUeRightNeutral->Fill(ptTrig,rat);\r
1390           fhPtHbpUeRightNeutral->Fill(ptTrig,cosi);\r
1391         }\r
1392       } //several UE calculation\r
1393           }\r
1394     else{\r
1395       nrefs++;\r
1396       if(nrefs==1){\r
1397         refpi0 = new TObjArray(0);\r
1398         refpi0->SetName(GetAODObjArrayName()+"Pi0s");\r
1399         refpi0->SetOwner(kFALSE);\r
1400       }\r
1401       refpi0->Add(pi0);\r
1402     }//put references in trigger AOD \r
1403       \r
1404      //if(GetDebug() > 2 ) printf("AliAnaParticleHadronCorrelation::MakeNeutralCorrelation() - Selected neutral for momentum imbalance: pt %2.2f, phi %2.2f, eta %2.2f \n",pt,phi,eta);\r
1405       \r
1406     }//loop\r
1407 }\r
1408   \r
1409 \r
1410 //____________________________________________________________________________\r
1411 //Bool_t  AliAnaParticleHadronCorrelation::SelectCluster(AliVCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg) {\r
1412 //  //Select cluster depending on its pid and acceptance selections\r
1413 //  \r
1414 //  //Skip matched clusters with tracks\r
1415 //  if(IsTrackMatched(calo)) return kFALSE;\r
1416 //  \r
1417 //  TString detector = "";\r
1418 //  if     (calo->IsPHOS())  detector= "PHOS";\r
1419 //  else if(calo->IsEMCAL()) detector= "EMCAL";\r
1420 //              \r
1421 //  //Check PID\r
1422 //  calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line\r
1423 //  pdg = AliCaloPID::kPhoton;   \r
1424 //  if(IsCaloPIDOn()){\r
1425 //    //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)\r
1426 //    //or redo PID, recommended option for EMCal.\r
1427 //    \r
1428 //    if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )\r
1429 //      pdg = GetCaloPID()->GetPdg(detector,calo->GetPID(),mom.E());//PID with weights\r
1430 //    else\r
1431 //      pdg = GetCaloPID()->GetPdg(detector,mom,calo);//PID recalculated\r
1432 //    \r
1433 //    if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - PDG of identified particle %d\n",pdg);\r
1434 //    \r
1435 //    //If it does not pass pid, skip\r
1436 //    if(pdg != AliCaloPID::kPhoton && pdg != AliCaloPID::kPi0) {\r
1437 //      return kFALSE ;\r
1438 //    }\r
1439 //  }//PID on\r
1440 //  \r
1441 //  //Check acceptance selection\r
1442 //  if(IsFiducialCutOn()){\r
1443 //    Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,detector) ;\r
1444 //    if(! in ) return kFALSE ;\r
1445 //  }\r
1446 //  \r
1447 //  if(GetDebug() > 5) printf("AliAnaParticleHadronCorrelation::SelectCluster() - Correlation photon selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);\r
1448 //  \r
1449 //  return kTRUE;\r
1450 //  \r
1451 //}\r