PWG4/PartCorrDep/AliAnaParticleJetLeadingConeCorrelation.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 /* $Id: $ */
  16
  17 //_________________________________________________________________________
  18 // Class that contains the algorithm for the reconstruction of jet, cone around leading particle
  19 // The seed is a backward particle (direct photon)
  20 // 1) Take the trigger particle stored in AliAODPWG4ParticleCorrelation,
  21 // 2) Search for the highest pt leading particle opposite to the photon within a phi, pt window
  22 // 3) Take all particles around leading in a cone R with pt larger than threshold and construct the jet
  23 //
  24 //  Class created from old AliPHOSGammaJet
  25 //  (see AliRoot versions previous Release 4-09)
  26 //
  27 //*-- Author: Gustavo Conesa (LNF-INFN)
  28 //////////////////////////////////////////////////////////////////////////////
  29
  30
  31 // --- ROOT system ---
  32 #include "TH2F.h"
  33 #include "TClonesArray.h"
  34 #include "TClass.h"
  35 //#include "Riostream.h"
  36
  37 //---- Analysis system ----
  38 #include "AliVTrack.h"
  39 #include "AliVCluster.h"
  40 #include "AliCaloTrackReader.h"
  41 #include "AliNeutralMesonSelection.h"
  42 #include "AliAnaParticleJetLeadingConeCorrelation.h"
  43 #include "AliCaloPID.h"
  44 #include "AliAODPWG4ParticleCorrelation.h"
  45 #include "AliFiducialCut.h"
  46
  47 ClassImp(AliAnaParticleJetLeadingConeCorrelation)
  48
  49
  50 //____________________________________________________________________________
  51   AliAnaParticleJetLeadingConeCorrelation::AliAnaParticleJetLeadingConeCorrelation() :
  52     AliAnaPartCorrBaseClass(), fJetsOnlyInCTS(kFALSE), fPbPb(kFALSE),
  53     fSeveralConeAndPtCuts(0),  fReMakeJet(0),
  54     fDeltaPhiMaxCut(0.), fDeltaPhiMinCut(0.),
  55     fLeadingRatioMaxCut(0.),  fLeadingRatioMinCut(0.),
  56     fJetCTSRatioMaxCut(0.), fJetCTSRatioMinCut(0.),
  57     fJetRatioMaxCut(0.),  fJetRatioMinCut(0.),
  58     fJetNCone(0),fJetNPt(0), fJetCone(0),
  59     fJetPtThreshold(0),fJetPtThresPbPb(0),
  60     fPtTriggerSelectionCut(0.0), fSelect(0),
  61     fSelectIsolated(0),
  62     //Histograms
  63     fOutCont(0x0),
  64     //Leading
  65     fhChargedLeadingPt(0),fhChargedLeadingPhi(0),fhChargedLeadingEta(0),
  66     fhChargedLeadingDeltaPt(0),fhChargedLeadingDeltaPhi(0),fhChargedLeadingDeltaEta(0),
  67     fhChargedLeadingRatioPt(0),
  68     fhNeutralLeadingPt(0),fhNeutralLeadingPhi(0),fhNeutralLeadingEta(0),
  69     fhNeutralLeadingDeltaPt(0),fhNeutralLeadingDeltaPhi(0),fhNeutralLeadingDeltaEta(0),
  70     fhNeutralLeadingRatioPt(0),fhChargedLeadingXi(0), fhNeutralLeadingXi(0),
  71     fhChargedLeadingDeltaPhiRatioPt30(0), fhNeutralLeadingDeltaPhiRatioPt30(0),
  72     fhChargedLeadingDeltaPhiRatioPt50(0), fhNeutralLeadingDeltaPhiRatioPt50(0),
  73     //Jet
  74     fhJetPt(0),fhJetRatioPt(0),fhJetDeltaPhi(0), fhJetDeltaEta(0),
  75     fhJetLeadingRatioPt(0),fhJetLeadingDeltaPhi(0),fhJetLeadingDeltaEta(0),
  76     fhJetFFz(0),fhJetFFxi(0),fhJetFFpt(0),fhJetNTracksInCone(0),
  77     fhBkgPt(0),fhBkgRatioPt(0),fhBkgDeltaPhi(0), fhBkgDeltaEta(0),
  78     fhBkgLeadingRatioPt(0),fhBkgLeadingDeltaPhi(0),fhBkgLeadingDeltaEta(0),
  79     fhBkgFFz(0),fhBkgFFxi(0),fhBkgFFpt(0),fhBkgNTracksInCone(0),
  80     //Several cones and thres histograms
  81     fhJetPts(),fhJetRatioPts(),fhJetDeltaPhis(), fhJetDeltaEtas(),
  82     fhJetLeadingRatioPts(),fhJetLeadingDeltaPhis(),fhJetLeadingDeltaEtas(),
  83     fhJetFFzs(),fhJetFFxis(),fhJetFFpts(),fhJetNTracksInCones(),
  84     fhBkgPts(),fhBkgRatioPts(),fhBkgDeltaPhis(), fhBkgDeltaEtas(),
  85     fhBkgLeadingRatioPts(),fhBkgLeadingDeltaPhis(),fhBkgLeadingDeltaEtas(),
  86     fhBkgFFzs(),fhBkgFFxis(),fhBkgFFpts(),fhBkgNTracksInCones()
  87 {
  88   //Default Ctor
  89
  90   //Initialize parameters
  91
  92   for(Int_t i = 0; i<6; i++){
  93     fJetXMin1[i]     = 0.0 ;
  94     fJetXMin2[i]     = 0.0 ;
  95     fJetXMax1[i]     = 0.0 ;
  96     fJetXMax2[i]     = 0.0 ;
  97     fBkgMean[i]      = 0.0 ;
  98     fBkgRMS[i]       = 0.0 ;
  99     if( i < 2 ){
 100       fJetE1[i]        = 0.0 ;
 101       fJetE2[i]        = 0.0 ;
 102       fJetSigma1[i]    = 0.0 ;
 103       fJetSigma2[i]    = 0.0 ;
 104     }
 105   }
 106
 107   //Several cones and thres histograms
 108   for(Int_t i = 0; i<5; i++){
 109     fJetCones[i]         = 0.0 ;
 110     fJetNameCones[i]     = ""  ;
 111     fJetPtThres[i]      = 0.0 ;
 112     fJetNamePtThres[i]  = ""  ;
 113     for(Int_t j = 0; j<5; j++){
 114       fhJetPts[i][j]=0 ;
 115       fhJetRatioPts[i][j]=0 ;
 116       fhJetDeltaPhis[i][j]=0 ;
 117       fhJetDeltaEtas[i][j]=0 ;
 118       fhJetLeadingRatioPts[i][j]=0 ;
 119       fhJetLeadingDeltaPhis[i][j]=0 ;
 120       fhJetLeadingDeltaEtas[i][j]=0 ;
 121       fhJetFFzs[i][j]=0 ;
 122       fhJetFFxis[i][j]=0 ;
 123       fhJetFFpts[i][j]=0 ;
 124       fhJetNTracksInCones[i][j]=0 ;
 125       fhBkgPts[i][j]=0 ;
 126       fhBkgRatioPts[i][j]=0 ;
 127       fhBkgDeltaPhis[i][j]=0 ;
 128       fhBkgDeltaEtas[i][j]=0 ;
 129       fhBkgLeadingRatioPts[i][j]=0 ;
 130       fhBkgLeadingDeltaPhis[i][j]=0 ;
 131       fhBkgLeadingDeltaEtas[i][j]=0 ;
 132       fhBkgFFzs[i][j]=0 ;
 133       fhBkgFFxis[i][j]=0 ;
 134       fhBkgFFpts[i][j]=0 ;
 135       fhBkgNTracksInCones[i][j]=0 ;
 136     }
 137   }
 138
 139   InitParameters();
 140
 141 }
 142 /*
 143 //____________________________________________________________________________
 144 AliAnaParticleJetLeadingConeCorrelation::AliAnaParticleJetLeadingConeCorrelation(const AliAnaParticleJetLeadingConeCorrelation & jetlc) :
 145   AliAnaPartCorrBaseClass(jetlc), fJetsOnlyInCTS(jetlc.fJetsOnlyInCTS), fPbPb(jetlc.fPbPb),
 146   fSeveralConeAndPtCuts(jetlc.fSeveralConeAndPtCuts),  fReMakeJet(jetlc. fReMakeJet),
 147   fDeltaPhiMaxCut(jetlc. fDeltaPhiMaxCut), fDeltaPhiMinCut(jetlc.fDeltaPhiMinCut),
 148   fLeadingRatioMaxCut(jetlc.fLeadingRatioMaxCut),  fLeadingRatioMinCut(jetlc.fLeadingRatioMinCut),
 149   fJetCTSRatioMaxCut(jetlc.fJetCTSRatioMaxCut),
 150   fJetCTSRatioMinCut(jetlc.fJetCTSRatioMinCut), fJetRatioMaxCut(jetlc.fJetRatioMaxCut),
 151   fJetRatioMinCut(jetlc.fJetRatioMinCut),  fJetNCone(jetlc.fJetNCone),
 152   fJetNPt(jetlc.fJetNPt), fJetCone(jetlc.fJetCone),
 153   fJetPtThreshold(jetlc.fJetPtThreshold),fJetPtThresPbPb(jetlc.fJetPtThresPbPb),
 154   fPtTriggerSelectionCut(jetlc.fPtTriggerSelectionCut), fSelect(jetlc.fSelect),
 155   fSelectIsolated(jetlc.fSelectIsolated),
 156   //Histograms
 157   fOutCont(jetlc. fOutCont),
 158   //Leading
 159   fhChargedLeadingPt(jetlc.fhChargedLeadingPt), fhChargedLeadingPhi(jetlc.fhChargedLeadingPhi),
 160   fhChargedLeadingEta(jetlc.fhChargedLeadingEta), fhChargedLeadingDeltaPt(jetlc.fhChargedLeadingDeltaPt),
 161   fhChargedLeadingDeltaPhi(jetlc.fhChargedLeadingDeltaPhi),fhChargedLeadingDeltaEta(jetlc.fhChargedLeadingDeltaEta),
 162   fhChargedLeadingRatioPt(jetlc.fhChargedLeadingRatioPt),
 163   fhNeutralLeadingPt(jetlc.fhNeutralLeadingPt),fhNeutralLeadingPhi(jetlc.fhNeutralLeadingPhi),
 164   fhNeutralLeadingEta(jetlc.fhNeutralLeadingEta), fhNeutralLeadingDeltaPt(jetlc.fhNeutralLeadingDeltaPt),
 165   fhNeutralLeadingDeltaPhi(jetlc.fhNeutralLeadingDeltaPhi),fhNeutralLeadingDeltaEta(jetlc.fhNeutralLeadingDeltaEta),
 166   fhNeutralLeadingRatioPt(jetlc.fhNeutralLeadingRatioPt),
 167   fhChargedLeadingXi(jetlc.fhChargedLeadingXi), fhNeutralLeadingXi(jetlc.fhNeutralLeadingXi),
 168   fhChargedLeadingDeltaPhiRatioPt30(jetlc.fhChargedLeadingDeltaPhiRatioPt30), fhNeutralLeadingDeltaPhiRatioPt30(jetlc.fhNeutralLeadingDeltaPhiRatioPt30),
 169   fhChargedLeadingDeltaPhiRatioPt50(jetlc.fhChargedLeadingDeltaPhiRatioPt50), fhNeutralLeadingDeltaPhiRatioPt50(jetlc.fhNeutralLeadingDeltaPhiRatioPt50),
 170   //Jet
 171   fhJetPt(jetlc.fhJetPt),fhJetRatioPt(jetlc.fhJetRatioPt),fhJetDeltaPhi(jetlc.fhJetDeltaPhi),
 172   fhJetDeltaEta(jetlc.fhJetDeltaEta), fhJetLeadingRatioPt(jetlc.fhJetLeadingRatioPt),
 173   fhJetLeadingDeltaPhi(jetlc.fhJetLeadingDeltaPhi),fhJetLeadingDeltaEta(jetlc.fhJetLeadingDeltaEta),
 174   fhJetFFz(jetlc.fhJetFFz),fhJetFFxi(jetlc.fhJetFFxi),fhJetFFpt(jetlc.fhJetFFpt),
 175   fhJetNTracksInCone(jetlc.fhJetNTracksInCone),
 176   fhBkgPt(jetlc.fhBkgPt),fhBkgRatioPt(jetlc.fhBkgRatioPt),fhBkgDeltaPhi(jetlc.fhBkgDeltaPhi),
 177   fhBkgDeltaEta(jetlc.fhBkgDeltaEta), fhBkgLeadingRatioPt(jetlc.fhBkgLeadingRatioPt),
 178   fhBkgLeadingDeltaPhi(jetlc.fhBkgLeadingDeltaPhi),fhBkgLeadingDeltaEta(jetlc.fhBkgLeadingDeltaEta),
 179   fhBkgFFz(jetlc.fhBkgFFz),fhBkgFFxi(jetlc.fhBkgFFxi),fhBkgFFpt(jetlc.fhBkgFFpt),
 180   fhBkgNTracksInCone(jetlc.fhBkgNTracksInCone),
 181   //Several cones and thres histograms
 182   fhJetPts(),fhJetRatioPts(),fhJetDeltaPhis(), fhJetDeltaEtas(),
 183   fhJetLeadingRatioPts(),fhJetLeadingDeltaPhis(),fhJetLeadingDeltaEtas(),
 184   fhJetFFzs(),fhJetFFxis(),fhJetFFpts(),fhJetNTracksInCones(),
 185   fhBkgPts(),fhBkgRatioPts(),fhBkgDeltaPhis(), fhBkgDeltaEtas(),
 186   fhBkgLeadingRatioPts(),fhBkgLeadingDeltaPhis(),fhBkgLeadingDeltaEtas(),
 187   fhBkgFFzs(),fhBkgFFxis(),fhBkgFFpts(),fhBkgNTracksInCones()
 188 {
 189   // cpy ctor
 190
 191   for(Int_t i = 0; i<6; i++){
 192     fJetXMin1[i]       = jetlc.fJetXMin1[i] ;
 193     fJetXMin2[i]       = jetlc.fJetXMin2[i] ;
 194     fJetXMax1[i]       = jetlc.fJetXMax1[i] ;
 195     fJetXMax2[i]       = jetlc.fJetXMax2[i] ;
 196     fBkgMean[i]        = jetlc.fBkgMean[i] ;
 197     fBkgRMS[i]         = jetlc.fBkgRMS[i] ;
 198     if( i < 2 ){
 199       fJetE1[i]        = jetlc.fJetE1[i] ;
 200       fJetE2[i]        = jetlc.fJetE2[i] ;
 201       fJetSigma1[i]    = jetlc.fJetSigma1[i] ;
 202       fJetSigma2[i]    = jetlc.fJetSigma2[i] ;
 203     }
 204   }
 205
 206   //Several cones and thres histograms
 207   for(Int_t i = 0; i<5; i++){
 208     fJetCones[i]        = jetlc.fJetCones[i] ;
 209     fJetNameCones[i]    = jetlc.fJetNameCones[i] ;
 210     fJetPtThres[i]      = jetlc.fJetPtThres[i] ;
 211     fJetNamePtThres[i]  = jetlc.fJetNamePtThres[i] ;
 212     for(Int_t j = 0; j<5; j++){
 213       fhJetPts[i][j] = jetlc.fhJetPts[i][j] ;
 214       fhJetRatioPts[i][j] = jetlc.fhJetRatioPts[i][j] ;
 215       fhJetDeltaPhis[i][j] = jetlc.fhJetDeltaPhis[i][j] ;
 216       fhJetDeltaEtas[i][j] = jetlc.fhJetDeltaEtas[i][j] ;
 217       fhJetLeadingRatioPts[i][j] = jetlc.fhJetLeadingRatioPts[i][j] ;
 218       fhJetLeadingDeltaPhis[i][j] = jetlc.fhJetLeadingDeltaPhis[i][j] ;
 219       fhJetLeadingDeltaEtas[i][j] = jetlc.fhJetLeadingDeltaEtas[i][j] ;
 220       fhJetFFzs[i][j] = jetlc.fhJetFFzs[i][j] ;
 221       fhJetFFxis[i][j] = jetlc.fhJetFFxis[i][j] ;
 222       fhJetFFpts[i][j] = jetlc.fhJetFFpts[i][j] ;
 223       fhJetNTracksInCones[i][j] = fhJetNTracksInCones[i][j] ;
 224       fhBkgPts[i][j] = jetlc.fhBkgPts[i][j] ;
 225       fhBkgRatioPts[i][j] = jetlc.fhBkgRatioPts[i][j] ;
 226       fhBkgDeltaPhis[i][j] = jetlc.fhBkgDeltaPhis[i][j] ;
 227       fhBkgDeltaEtas[i][j] = jetlc.fhBkgDeltaEtas[i][j] ;
 228       fhBkgLeadingRatioPts[i][j] = jetlc.fhBkgLeadingRatioPts[i][j] ;
 229       fhBkgLeadingDeltaPhis[i][j] = jetlc.fhBkgLeadingDeltaPhis[i][j] ;
 230       fhBkgLeadingDeltaEtas[i][j] = jetlc.fhBkgLeadingDeltaEtas[i][j] ;
 231       fhBkgFFzs[i][j] = jetlc.fhBkgFFzs[i][j] ;
 232       fhBkgFFxis[i][j] = jetlc.fhBkgFFxis[i][j] ;
 233       fhBkgFFpts[i][j] = jetlc.fhBkgFFpts[i][j] ;
 234       fhBkgNTracksInCones[i][j] = jetlc.fhBkgNTracksInCones[i][j] ;
 235     }
 236   }
 237 }
 238
 239 //_________________________________________________________________________
 240 AliAnaParticleJetLeadingConeCorrelation & AliAnaParticleJetLeadingConeCorrelation::operator = (const AliAnaParticleJetLeadingConeCorrelation & jetlc)
 241 {
 242   // assignment operator
 243
 244   if(this == &jetlc)return *this;
 245   ((AliAnaPartCorrBaseClass *)this)->operator=(jetlc);
 246
 247   fSeveralConeAndPtCuts  = jetlc.fSeveralConeAndPtCuts ;
 248   fPbPb                  = jetlc.fPbPb ;
 249   fReMakeJet             = jetlc.fReMakeJet ;
 250   fJetsOnlyInCTS         = jetlc.fJetsOnlyInCTS;
 251
 252   fDeltaPhiMaxCut        = jetlc.fDeltaPhiMaxCut ;
 253   fDeltaPhiMinCut        = jetlc.fDeltaPhiMinCut ;
 254   fLeadingRatioMaxCut    = jetlc.fLeadingRatioMaxCut ;
 255   fLeadingRatioMinCut    = jetlc.fLeadingRatioMinCut ;
 256
 257   fJetCTSRatioMaxCut     = jetlc.fJetCTSRatioMaxCut ;
 258   fJetCTSRatioMinCut     = jetlc.fJetCTSRatioMinCut ;
 259   fJetRatioMaxCut        = jetlc.fJetRatioMaxCut ;
 260   fJetRatioMinCut        = jetlc.fJetRatioMinCut ;
 261
 262   fJetNCone              = jetlc.fJetNCone ;
 263   fJetNPt                = jetlc.fJetNPt ; fJetCone = jetlc.fJetCone ;
 264   fJetPtThreshold        = jetlc.fJetPtThreshold ;
 265   fJetPtThresPbPb        = jetlc.fJetPtThresPbPb ;
 266   fPtTriggerSelectionCut = jetlc.fPtTriggerSelectionCut ;
 267   fSelect                = jetlc.fSelect ;
 268   fSelectIsolated        = jetlc.fSelectIsolated ;
 269
 270   for(Int_t i = 0; i<6; i++){
 271     fJetXMin1[i]       = jetlc.fJetXMin1[i] ;
 272     fJetXMin2[i]       = jetlc.fJetXMin2[i] ;
 273     fJetXMax1[i]       = jetlc.fJetXMax1[i] ;
 274     fJetXMax2[i]       = jetlc.fJetXMax2[i] ;
 275     fBkgMean[i]        = jetlc.fBkgMean[i] ;
 276     fBkgRMS[i]         = jetlc.fBkgRMS[i] ;
 277     if( i < 2 ){
 278       fJetE1[i]        = jetlc.fJetE1[i] ;
 279       fJetE2[i]        = jetlc.fJetE2[i] ;
 280       fJetSigma1[i]    = jetlc.fJetSigma1[i] ;
 281       fJetSigma2[i]    = jetlc.fJetSigma2[i] ;
 282     }
 283   }
 284
 285   //Histograms
 286   fOutCont = jetlc. fOutCont ;
 287   fhChargedLeadingPt = jetlc.fhChargedLeadingPt; fhChargedLeadingPhi = jetlc.fhChargedLeadingPhi;
 288   fhChargedLeadingEta = jetlc.fhChargedLeadingEta; fhChargedLeadingDeltaPt = jetlc.fhChargedLeadingDeltaPt;
 289   fhChargedLeadingDeltaPhi = jetlc.fhChargedLeadingDeltaPhi;fhChargedLeadingDeltaEta = jetlc.fhChargedLeadingDeltaEta;
 290   fhChargedLeadingRatioPt = jetlc.fhChargedLeadingRatioPt;
 291   fhNeutralLeadingPt = jetlc.fhNeutralLeadingPt;fhNeutralLeadingPhi = jetlc.fhNeutralLeadingPhi;
 292   fhNeutralLeadingEta = jetlc.fhNeutralLeadingEta; fhNeutralLeadingDeltaPt = jetlc.fhNeutralLeadingDeltaPt;
 293   fhNeutralLeadingDeltaPhi = jetlc.fhNeutralLeadingDeltaPhi;fhNeutralLeadingDeltaEta = jetlc.fhNeutralLeadingDeltaEta;
 294   fhNeutralLeadingRatioPt = jetlc.fhNeutralLeadingRatioPt;
 295   fhChargedLeadingXi = jetlc.fhChargedLeadingXi;
 296   fhNeutralLeadingXi = jetlc.fhNeutralLeadingXi;
 297
 298   fhChargedLeadingDeltaPhiRatioPt30 = jetlc.fhChargedLeadingDeltaPhiRatioPt30;
 299   fhNeutralLeadingDeltaPhiRatioPt30 = jetlc.fhNeutralLeadingDeltaPhiRatioPt30;
 300   fhChargedLeadingDeltaPhiRatioPt50 = jetlc.fhChargedLeadingDeltaPhiRatioPt50;
 301   fhNeutralLeadingDeltaPhiRatioPt50 = jetlc.fhNeutralLeadingDeltaPhiRatioPt50;
 302
 303   fhJetPt = jetlc.fhJetPt;fhJetRatioPt = jetlc.fhJetRatioPt;fhJetDeltaPhi = jetlc.fhJetDeltaPhi;
 304   fhJetDeltaEta = jetlc.fhJetDeltaEta; fhJetLeadingRatioPt = jetlc.fhJetLeadingRatioPt;
 305   fhJetLeadingDeltaPhi = jetlc.fhJetLeadingDeltaPhi;fhJetLeadingDeltaEta = jetlc.fhJetLeadingDeltaEta;
 306   fhJetFFz = jetlc.fhJetFFz;fhJetFFxi = jetlc.fhJetFFxi;fhJetFFpt = jetlc.fhJetFFpt;
 307   fhJetNTracksInCone = jetlc.fhJetNTracksInCone;
 308   fhBkgPt = jetlc.fhBkgPt;fhBkgRatioPt = jetlc.fhBkgRatioPt;fhBkgDeltaPhi = jetlc.fhBkgDeltaPhi;
 309   fhBkgDeltaEta = jetlc.fhBkgDeltaEta; fhBkgLeadingRatioPt = jetlc.fhBkgLeadingRatioPt;
 310   fhBkgLeadingDeltaPhi = jetlc.fhBkgLeadingDeltaPhi;fhBkgLeadingDeltaEta = jetlc.fhBkgLeadingDeltaEta;
 311   fhBkgFFz = jetlc.fhBkgFFz;fhBkgFFxi = jetlc.fhBkgFFxi;fhBkgFFpt = jetlc.fhBkgFFpt;
 312   fhBkgNTracksInCone = jetlc.fhBkgNTracksInCone;
 313
 314
 315   //Several cones and thres histograms
 316   for(Int_t i = 0; i<5; i++){
 317     fJetCones[i]        = jetlc.fJetCones[i] ;
 318     fJetNameCones[i]    = jetlc.fJetNameCones[i] ;
 319     fJetPtThres[i]      = jetlc.fJetPtThres[i] ;
 320     fJetNamePtThres[i]  = jetlc.fJetNamePtThres[i] ;
 321
 322     for(Int_t j = 0; j<5; j++){
 323       fhJetPts[i][j] = jetlc.fhJetPts[i][j] ;
 324       fhJetRatioPts[i][j] = jetlc.fhJetRatioPts[i][j] ;
 325       fhJetDeltaPhis[i][j] = jetlc.fhJetDeltaPhis[i][j] ;
 326       fhJetDeltaEtas[i][j] = jetlc.fhJetDeltaEtas[i][j] ;
 327       fhJetLeadingRatioPts[i][j] = jetlc.fhJetLeadingRatioPts[i][j] ;
 328       fhJetLeadingDeltaPhis[i][j] = jetlc.fhJetLeadingDeltaPhis[i][j] ;
 329       fhJetLeadingDeltaEtas[i][j] = jetlc.fhJetLeadingDeltaEtas[i][j] ;
 330       fhJetFFzs[i][j] = jetlc.fhJetFFzs[i][j] ;
 331       fhJetFFxis[i][j] = jetlc.fhJetFFxis[i][j] ;
 332       fhJetFFpts[i][j] = jetlc.fhJetFFpts[i][j] ;
 333       fhJetNTracksInCones[i][j] = fhJetNTracksInCones[i][j] ;
 334       fhBkgPts[i][j] = jetlc.fhBkgPts[i][j] ;
 335       fhBkgRatioPts[i][j] = jetlc.fhBkgRatioPts[i][j] ;
 336       fhBkgDeltaPhis[i][j] = jetlc.fhBkgDeltaPhis[i][j] ;
 337       fhBkgDeltaEtas[i][j] = jetlc.fhBkgDeltaEtas[i][j] ;
 338       fhBkgLeadingRatioPts[i][j] = jetlc.fhBkgLeadingRatioPts[i][j] ;
 339       fhBkgLeadingDeltaPhis[i][j] = jetlc.fhBkgLeadingDeltaPhis[i][j] ;
 340       fhBkgLeadingDeltaEtas[i][j] = jetlc.fhBkgLeadingDeltaEtas[i][j] ;
 341       fhBkgFFzs[i][j] = jetlc.fhBkgFFzs[i][j] ;
 342       fhBkgFFxis[i][j] = jetlc.fhBkgFFxis[i][j] ;
 343       fhBkgFFpts[i][j] = jetlc.fhBkgFFpts[i][j] ;
 344       fhBkgNTracksInCones[i][j] = jetlc.fhBkgNTracksInCones[i][j] ;
 345     }
 346   }
 347
 348   return *this;
 349
 350 }
 351 */
 352 //____________________________________________________________________________
 353 AliAnaParticleJetLeadingConeCorrelation::~AliAnaParticleJetLeadingConeCorrelation()
 354 {
 355    // Remove all pointers except analysis output pointers.
 356 }
 357
 358 //____________________________________________________________________________
 359 Double_t AliAnaParticleJetLeadingConeCorrelation::CalculateJetRatioLimit(const Double_t ptg, const Double_t *par, const Double_t *x) const {
 360   //Calculate the ratio of the jet and trigger particle limit for the selection
 361   //WARNING: need to check what it does
 362   //printf("CalculateLimit: x1 %2.3f, x2%2.3f\n",x[0],x[1]);
 363   Double_t ePP = par[0] + par[1] * ptg ;
 364   Double_t sPP = par[2] + par[3] * ptg ;
 365   Double_t f   = x[0]   + x[1]   * ptg ;
 366   Double_t ePbPb = ePP + par[4] ;
 367   Double_t sPbPb = TMath::Sqrt(sPP*sPP+ par[5]*par[5]) ;
 368   Double_t rat = (ePbPb - sPbPb * f) / ptg ;
 369   //printf("CalculateLimit: ePP %2.3f, sPP %2.3f, f %2.3f\n", ePP, sPP, f);
 370   //printf("CalculateLimit: ePbPb %2.3f, sPbPb %2.3f, rat %2.3f\n", ePbPb, sPbPb, rat);
 371   return rat ;
 372 }
 373
 374 //____________________________________________________________________________
 375 void AliAnaParticleJetLeadingConeCorrelation::FillJetHistos(AliAODPWG4ParticleCorrelation * particle, const TLorentzVector  leading, const TLorentzVector jet, const TString type, const TString lastname)
 376 {
 377   //Fill jet and background histograms
 378   Double_t ptTrig = particle->Pt();
 379   Double_t ptJet = jet.Pt();
 380   Double_t ptLead = leading.Pt();
 381   Double_t phiTrig = particle->Phi();
 382   Double_t phiJet = jet.Phi();
 383   if(phiJet < 0) phiJet+=TMath::TwoPi();
 384   Double_t phiLead = leading.Phi();
 385   if(phiLead < 0) phiLead+=TMath::TwoPi();
 386   Double_t etaTrig = particle->Eta();
 387   Double_t etaJet = jet.Eta();
 388   Double_t etaLead = leading.Eta();
 389
 390   TH2F *h1 = 0x0;
 391   h1 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 392   if(h1)h1->Fill(ptTrig,ptJet);
 393
 394   TH2F *h2 = 0x0;
 395   h2 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sRatioPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 396   if(h2) h2->Fill(ptTrig,ptJet/ptTrig);
 397
 398   TH2F *h3 = 0x0;
 399   h3 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingRatioPt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 400   if(h3)h3->Fill(ptTrig,ptLead/ptJet);
 401
 402   //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))->
 403   //     Fill(ptTrig,phiJet);
 404   TH2F *h4 = 0x0;
 405   h4 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sDeltaPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 406   if(h4) h4->Fill(ptTrig,phiJet-phiTrig);
 407   TH2F *h5 = 0x0;
 408   h5 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingDeltaPhi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 409   if(h5) h5->Fill(ptTrig,phiJet-phiLead);
 410
 411   //   dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())))->
 412   //     Fill(ptTrig,etaJet);
 413   TH2F *h6 = 0x0;
 414   h6 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sDeltaEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 415   if(h6) h6->Fill(ptTrig,etaJet-etaTrig);
 416   TH2F *h7 = 0x0;
 417   h7 = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sLeadingDeltaEta%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 418   if(h7) h7->Fill(ptTrig,etaJet-etaLead);
 419
 420   //Construct fragmentation function
 421   TObjArray * pl = new TObjArray;
 422
 423   if(type == "Jet") pl = particle->GetObjArray(Form("%sTracks",GetAODObjArrayName().Data()));
 424   else if(type == "Bkg") particle->GetObjArray(Form("%sTracksBkg",GetAODObjArrayName().Data()));
 425
 426   if(!pl) return ;
 427
 428   //Different pt cut for jet particles in different collisions systems
 429   //Only needed when jet is recalculated from AODs
 430   Float_t ptcut = fJetPtThreshold;
 431   if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
 432
 433   TVector3 p3;
 434   Int_t nTracksInCone = 0;
 435
 436   for(Int_t ipr = 0;ipr < pl->GetEntriesFast() ; ipr ++ ){
 437     AliVTrack* track = dynamic_cast<AliVTrack *>(pl->At(ipr)) ;
 438     if(track)p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 439     else printf("AliAnaParticleJetLeadingConeCorrelation::FillJetHistos() - Track not available\n");
 440
 441     //Recheck if particle is in jet cone
 442     if(fReMakeJet || fSeveralConeAndPtCuts)
 443       if(!IsParticleInJetCone(p3.Eta(), p3.Phi(), leading.Eta(), leading.Phi()) ) continue ;
 444
 445     nTracksInCone++;
 446
 447     TH2F *ha =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFz%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 448     if(ha) ha->Fill(ptTrig,p3.Pt()/ptTrig);
 449     TH2F *hb  =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFxi%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 450     if(hb) hb->Fill(ptTrig,TMath::Log(ptTrig/p3.Pt()));
 451     TH2F *hc =dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sFFpt%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 452     if(hc) hc->Fill(ptTrig,p3.Pt());
 453
 454   }//track loop
 455
 456   if(nTracksInCone > 0) {
 457     TH2F *hd = dynamic_cast<TH2F*>(GetOutputContainer()->FindObject(Form("%s%sNTracksInCone%s",GetAddedHistogramsStringToName().Data(),type.Data(),lastname.Data())));
 458     if(hd)hd->Fill(ptTrig, nTracksInCone);
 459   }
 460
 461 }
 462
 463 //________________________________________________________________________
 464 TList *  AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects()
 465 {
 466   // Create histograms to be saved in output file and
 467   // store them in fOutCont
 468
 469   if(GetDebug()>1) printf("AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects() - Init histograms \n");
 470
 471   fOutCont = new TList() ;
 472   fOutCont->SetName("ParticleJetLeadingInConeCorrelationHistograms") ;
 473
 474   Int_t nptbins  = GetHistoPtBins();
 475   Int_t nphibins = GetHistoPhiBins();
 476   Int_t netabins = GetHistoEtaBins();
 477   Float_t ptmax  = GetHistoPtMax();
 478   Float_t phimax = GetHistoPhiMax();
 479   Float_t etamax = GetHistoEtaMax();
 480   Float_t ptmin  = GetHistoPtMin();
 481   Float_t phimin = GetHistoPhiMin();
 482   Float_t etamin = GetHistoEtaMin();
 483
 484   fhChargedLeadingPt  = new TH2F("ChargedLeadingPt","p_{T leading charge} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 485   fhChargedLeadingPt->SetYTitle("p_{T leading charge}");
 486   fhChargedLeadingPt->SetXTitle("p_{T trigger} (GeV/c)");
 487
 488   fhChargedLeadingPhi  = new TH2F("ChargedLeadingPhi","#phi_{h^{#pm}}  vs p_{T trigger}", nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 489   fhChargedLeadingPhi->SetYTitle("#phi_{h^{#pm}} (rad)");
 490   fhChargedLeadingPhi->SetXTitle("p_{T trigger} (GeV/c)");
 491
 492   fhChargedLeadingEta  = new TH2F("ChargedLeadingEta","#eta_{h^{#pm}}  vs p_{T trigger}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 493   fhChargedLeadingEta->SetYTitle("#eta_{h^{#pm}} ");
 494   fhChargedLeadingEta->SetXTitle("p_{T trigger} (GeV/c)");
 495
 496   fhChargedLeadingDeltaPt  = new TH2F("ChargedLeadingDeltaPt","p_{T trigger} - p_{T h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 497   fhChargedLeadingDeltaPt->SetYTitle("#Delta p_{T} (GeV/c)");
 498   fhChargedLeadingDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
 499
 500   fhChargedLeadingDeltaPhi  = new TH2F("ChargedLeadingDeltaPhi","#phi_{trigger} - #phi_{h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 501   fhChargedLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
 502   fhChargedLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 503
 504   fhChargedLeadingDeltaEta  = new TH2F("ChargedLeadingDeltaEta","#eta_{trigger} - #eta_{h^{#pm}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 505   fhChargedLeadingDeltaEta->SetYTitle("#Delta #eta");
 506   fhChargedLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 507
 508   fhChargedLeadingRatioPt  = new TH2F("ChargedLeadingRatioPt","p_{T leading charge} /p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 509   fhChargedLeadingRatioPt->SetYTitle("p_{T lead charge} /p_{T trigger}");
 510   fhChargedLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 511
 512   fhChargedLeadingXi  = new TH2F("ChargedLeadingXi","ln(p_{T trigger} / p_{T leading charge} ) vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,10);
 513   fhChargedLeadingXi->SetYTitle("#xi");
 514   fhChargedLeadingXi->SetXTitle("p_{T trigger} (GeV/c)");
 515
 516   fOutCont->Add(fhChargedLeadingPt) ;
 517   fOutCont->Add(fhChargedLeadingPhi) ;
 518   fOutCont->Add(fhChargedLeadingEta) ;
 519   fOutCont->Add(fhChargedLeadingDeltaPt) ;
 520   fOutCont->Add(fhChargedLeadingDeltaPhi) ;
 521   fOutCont->Add(fhChargedLeadingDeltaEta) ;
 522   fOutCont->Add(fhChargedLeadingRatioPt) ;
 523   fOutCont->Add(fhChargedLeadingXi) ;
 524
 525   fhChargedLeadingDeltaPhiRatioPt30  = new TH2F("ChargedLeadingDeltaPhiRatioPt30","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, charged leading, p_{T trigger} > 30 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1);
 526   fhChargedLeadingDeltaPhiRatioPt30->SetXTitle("#Delta #phi (rad)");
 527   fhChargedLeadingDeltaPhiRatioPt30->SetYTitle("p_{T leading} / p_{T trigger}");
 528
 529   fhChargedLeadingDeltaPhiRatioPt50  = new TH2F("ChargedLeadingDeltaPhiRatioPt50","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, charged leading, p_{T trigger} > 50 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1);
 530   fhChargedLeadingDeltaPhiRatioPt50->SetXTitle("#Delta #phi (rad)");
 531   fhChargedLeadingDeltaPhiRatioPt50->SetYTitle("p_{T leading} / p_{T trigger}");
 532
 533   fOutCont->Add(fhChargedLeadingDeltaPhiRatioPt30) ;
 534   fOutCont->Add(fhChargedLeadingDeltaPhiRatioPt50) ;
 535
 536   if(!fJetsOnlyInCTS){
 537
 538     fhNeutralLeadingPt  = new TH2F("NeutralLeadingPt","p_{T leading #pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 539     fhNeutralLeadingPt->SetYTitle("p_{T leading #pi^{0}}");
 540     fhNeutralLeadingPt->SetXTitle("p_{T trigger} (GeV/c)");
 541
 542     fhNeutralLeadingPhi  = new TH2F("NeutralLeadingPhi","#phi_{#pi^{0}}  vs p_{T trigger}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 543     fhNeutralLeadingPhi->SetYTitle("#phi_{#pi^{0}} (rad)");
 544     fhNeutralLeadingPhi->SetXTitle("p_{T trigger} (GeV/c)");
 545
 546     fhNeutralLeadingEta  = new TH2F("NeutralLeadingEta","#eta_{#pi^{0}}  vs p_{T trigger}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 547     fhNeutralLeadingEta->SetYTitle("#eta_{#pi^{0}} ");
 548     fhNeutralLeadingEta->SetXTitle("p_{T trigger} (GeV/c)");
 549
 550     fhNeutralLeadingDeltaPt  = new TH2F("NeutralLeadingDeltaPt","p_{T trigger} - p_{T #pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 551     fhNeutralLeadingDeltaPt->SetYTitle("#Delta p_{T} (GeV/c)");
 552     fhNeutralLeadingDeltaPt->SetXTitle("p_{T trigger} (GeV/c)");
 553
 554     fhNeutralLeadingDeltaPhi  = new TH2F("NeutralLeadingDeltaPhi","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 555     fhNeutralLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
 556     fhNeutralLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 557
 558     fhNeutralLeadingDeltaEta  = new TH2F("NeutralLeadingDeltaEta","#eta_{trigger} - #eta_{#pi^{0}} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 559     fhNeutralLeadingDeltaEta->SetYTitle("#Delta #eta");
 560     fhNeutralLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 561
 562     fhNeutralLeadingRatioPt  = new TH2F("NeutralLeadingRatioPt","p_{T leading #pi^{0}} /p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 563     fhNeutralLeadingRatioPt->SetYTitle("p_{T lead #pi^{0}} /p_{T trigger}");
 564     fhNeutralLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 565
 566     fhNeutralLeadingXi  = new TH2F("NeutralLeadingXi","ln(p_{T trigger} / p_{T leading #pi^{0}} ) vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,10);
 567     fhNeutralLeadingXi->SetYTitle("#xi");
 568     fhNeutralLeadingXi->SetXTitle("p_{T trigger} (GeV/c)");
 569
 570     fOutCont->Add(fhNeutralLeadingPt) ;
 571     fOutCont->Add(fhNeutralLeadingPhi) ;
 572     fOutCont->Add(fhNeutralLeadingEta) ;
 573     fOutCont->Add(fhNeutralLeadingDeltaPt) ;
 574     fOutCont->Add(fhNeutralLeadingDeltaPhi) ;
 575     fOutCont->Add(fhNeutralLeadingDeltaEta) ;
 576     fOutCont->Add(fhNeutralLeadingRatioPt) ;
 577     fOutCont->Add(fhNeutralLeadingXi) ;
 578
 579         fhNeutralLeadingDeltaPhiRatioPt30  = new TH2F("NeutralLeadingDeltaPhiRatioPt30","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, neutral leading, p_{T trigger} > 30 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1);
 580         fhNeutralLeadingDeltaPhiRatioPt30->SetXTitle("#Delta #phi (rad)");
 581         fhNeutralLeadingDeltaPhiRatioPt30->SetYTitle("p_{T leading} / p_{T trigger}");
 582
 583         fhNeutralLeadingDeltaPhiRatioPt50  = new TH2F("NeutralLeadingDeltaPhiRatioPt50","#phi_{trigger} - #phi_{#pi^{0}} vs p_{T leading}/p_{T trigger}, neutral leading, p_{T trigger} > 50 GeV/c",120,0,TMath::TwoPi(),nptbins,0,1);
 584         fhNeutralLeadingDeltaPhiRatioPt50->SetXTitle("#Delta #phi (rad)");
 585         fhNeutralLeadingDeltaPhiRatioPt50->SetYTitle("p_{T leading} / p_{T trigger}");
 586         fOutCont->Add(fhNeutralLeadingDeltaPhiRatioPt30) ;
 587         fOutCont->Add(fhNeutralLeadingDeltaPhiRatioPt50) ;
 588
 589   }
 590
 591   if(!fSeveralConeAndPtCuts){// not several cones
 592
 593     //Jet Distributions
 594     fhJetPt  = new TH2F("JetPt","p_{T  jet} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 595     fhJetPt->SetYTitle("p_{T  jet}");
 596     fhJetPt->SetXTitle("p_{T trigger} (GeV/c)");
 597
 598     fhJetRatioPt  = new TH2F("JetRatioPt","p_{T  jet}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 599     fhJetRatioPt->SetYTitle("p_{T  jet}/p_{T trigger}");
 600     fhJetRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 601
 602     fhJetDeltaPhi  = new TH2F("JetDeltaPhi","#phi_{jet} - #phi_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 603     fhJetDeltaPhi->SetYTitle("#Delta #phi (rad)");
 604     fhJetDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 605
 606     fhJetDeltaEta  = new TH2F("JetDeltaEta","#eta_{jet} - #eta_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 607     fhJetDeltaEta->SetYTitle("#Delta #eta");
 608     fhJetDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 609
 610     fhJetLeadingRatioPt  = new TH2F("JetLeadingRatioPt","p_{T  jet} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 611     fhJetLeadingRatioPt->SetYTitle("p_{T  leading}/p_{T jet}");
 612     fhJetLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 613
 614     fhJetLeadingDeltaPhi  = new TH2F("JetLeadingDeltaPhi","#phi_{jet} - #phi_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-TMath::Pi(),TMath::Pi());
 615     fhJetLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
 616     fhJetLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 617
 618     fhJetLeadingDeltaEta  = new TH2F("JetLeadingDeltaEta","#eta_{jet} - #eta_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 619     fhJetLeadingDeltaEta->SetYTitle("#Delta #eta");
 620     fhJetLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 621
 622     fhJetFFz  = new TH2F("JetFFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,2);
 623     fhJetFFz->SetYTitle("z");
 624     fhJetFFz->SetXTitle("p_{T trigger}");
 625
 626     fhJetFFxi  = new TH2F("JetFFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}",nptbins,ptmin,ptmax,100,0.,10.);
 627     fhJetFFxi->SetYTitle("#xi");
 628     fhJetFFxi->SetXTitle("p_{T trigger}");
 629
 630     fhJetFFpt  = new TH2F("JetFFpt","#xi = p_{T i charged}) vs p_{T trigger}",nptbins,ptmin,ptmax,200,0.,50.);
 631     fhJetFFpt->SetYTitle("p_{T charged hadron}");
 632     fhJetFFpt->SetXTitle("p_{T trigger}");
 633
 634     fhJetNTracksInCone  = new TH2F("JetNTracksInCone","N particles in cone vs p_{T trigger}",nptbins,ptmin,ptmax,5000,0, 5000);
 635     fhJetNTracksInCone->SetYTitle("N tracks in jet cone");
 636     fhJetNTracksInCone->SetXTitle("p_{T trigger} (GeV/c)");
 637
 638     fOutCont->Add(fhJetPt) ;
 639     fOutCont->Add(fhJetRatioPt) ;
 640     fOutCont->Add(fhJetDeltaPhi) ;
 641     fOutCont->Add(fhJetDeltaEta) ;
 642     fOutCont->Add(fhJetLeadingRatioPt) ;
 643     fOutCont->Add(fhJetLeadingDeltaPhi) ;
 644     fOutCont->Add(fhJetLeadingDeltaEta) ;
 645     fOutCont->Add(fhJetFFz) ;
 646     fOutCont->Add(fhJetFFxi) ;
 647     fOutCont->Add(fhJetFFpt) ;
 648     fOutCont->Add(fhJetNTracksInCone) ;
 649
 650     //Bkg Distributions
 651     fhBkgPt  = new TH2F("BkgPt","p_{T  bkg} vs p_{T trigger}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 652     fhBkgPt->SetYTitle("p_{T  bkg}");
 653     fhBkgPt->SetXTitle("p_{T trigger} (GeV/c)");
 654
 655     fhBkgRatioPt  = new TH2F("BkgRatioPt","p_{T  bkg}/p_{T trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 656     fhBkgRatioPt->SetYTitle("p_{T  bkg}/p_{T trigger}");
 657     fhBkgRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 658
 659     fhBkgDeltaPhi  = new TH2F("BkgDeltaPhi","#phi_{bkg} - #phi_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 660     fhBkgDeltaPhi->SetYTitle("#Delta #phi (rad)");
 661     fhBkgDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 662
 663     fhBkgDeltaEta  = new TH2F("BkgDeltaEta","#eta_{bkg} - #eta_{trigger} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 664     fhBkgDeltaEta->SetYTitle("#Delta #eta");
 665     fhBkgDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 666
 667     fhBkgLeadingRatioPt  = new TH2F("BkgLeadingRatioPt","p_{T  bkg} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,2);
 668     fhBkgLeadingRatioPt->SetYTitle("p_{T  leading}/p_{T bkg}");
 669     fhBkgLeadingRatioPt->SetXTitle("p_{T trigger} (GeV/c)");
 670
 671     fhBkgLeadingDeltaPhi  = new TH2F("BkgLeadingDeltaPhi","#phi_{bkg} - #phi_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 672     fhBkgLeadingDeltaPhi->SetYTitle("#Delta #phi (rad)");
 673     fhBkgLeadingDeltaPhi->SetXTitle("p_{T trigger} (GeV/c)");
 674
 675     fhBkgLeadingDeltaEta  = new TH2F("BkgLeadingDeltaEta","#eta_{bkg} - #eta_{leading} vs p_{T trigger}",nptbins,ptmin,ptmax,120,-2,2);
 676     fhBkgLeadingDeltaEta->SetYTitle("#Delta #eta");
 677     fhBkgLeadingDeltaEta->SetXTitle("p_{T trigger} (GeV/c)");
 678
 679     fhBkgFFz  = new TH2F("BkgFFz","z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", nptbins,ptmin,ptmax,200,0.,2);
 680     fhBkgFFz->SetYTitle("z");
 681     fhBkgFFz->SetXTitle("p_{T trigger}");
 682
 683     fhBkgFFxi  = new TH2F("BkgFFxi","#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", nptbins,ptmin,ptmax,100,0.,10.);
 684     fhBkgFFxi->SetYTitle("#xi");
 685     fhBkgFFxi->SetXTitle("p_{T trigger}");
 686
 687     fhBkgFFpt  = new TH2F("BkgFFpt","p_{T charged hadron } vs p_{T trigger}", nptbins,ptmin,ptmax,200,0.,50.);
 688     fhBkgFFpt->SetYTitle("p_{T charged} hadron");
 689     fhBkgFFpt->SetXTitle("p_{T trigger}");
 690
 691     fhBkgNTracksInCone  = new TH2F("BkgNTracksInCone","N particles in cone vs p_{T trigger}",nptbins,ptmin,ptmax,5000,0, 5000);
 692     fhBkgNTracksInCone->SetYTitle("N tracks in bkg cone");
 693     fhBkgNTracksInCone->SetXTitle("p_{T trigger} (GeV/c)");
 694
 695     fOutCont->Add(fhBkgPt) ;
 696     fOutCont->Add(fhBkgRatioPt) ;
 697     fOutCont->Add(fhBkgDeltaPhi) ;
 698     fOutCont->Add(fhBkgDeltaEta) ;
 699     fOutCont->Add(fhBkgLeadingRatioPt) ;
 700     fOutCont->Add(fhBkgLeadingDeltaPhi) ;
 701     fOutCont->Add(fhBkgLeadingDeltaEta) ;
 702     fOutCont->Add(fhBkgFFz) ;
 703     fOutCont->Add(fhBkgFFxi) ;
 704     fOutCont->Add(fhBkgFFpt) ;
 705     fOutCont->Add(fhBkgNTracksInCone) ;
 706
 707   }//not several cones
 708   else{ //If we want to study the jet for different cones and pt
 709     for(Int_t icone = 0; icone<fJetNCone; icone++){//icone
 710       for(Int_t ipt = 0; ipt<fJetNPt;ipt++){ //ipt
 711
 712         TString lastnamehist ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt];
 713         TString lastnametitle =", cone ="+fJetNameCones[icone]+", pt > " +fJetNamePtThres[ipt]+" GeV/c";
 714
 715         //Jet Distributions
 716         fhJetPts[icone][ipt] = new TH2F(Form("JetPt%s",lastnamehist.Data()),Form("p_{T  jet} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 717         fhJetPts[icone][ipt]->SetYTitle("p_{T  jet}");
 718         fhJetPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 719
 720         fhJetRatioPts[icone][ipt] = new TH2F(Form("JetRatioPt%s",lastnamehist.Data()),Form("p_{T  jet}/p_{T trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2);
 721         fhJetRatioPts[icone][ipt]->SetYTitle("p_{T  jet}/p_{T trigger}");
 722         fhJetRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 723
 724         fhJetDeltaPhis[icone][ipt] = new TH2F(Form("JetDeltaPhi%s",lastnamehist.Data()),Form("#phi_{jet} - #phi_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 725         fhJetDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
 726         fhJetDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 727
 728         fhJetDeltaEtas[icone][ipt] = new TH2F(Form("JetDeltaEta%s",lastnamehist.Data()),Form("#eta_{jet} - #eta_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2);
 729         fhJetDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
 730         fhJetDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 731
 732         fhJetLeadingRatioPts[icone][ipt] = new TH2F(Form("JetLeadingRatioPt%s",lastnamehist.Data()),Form("p_{T  jet} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2);
 733         fhJetLeadingRatioPts[icone][ipt]->SetYTitle("p_{T  leading}/p_{T jet}");
 734         fhJetLeadingRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 735
 736         fhJetLeadingDeltaPhis[icone][ipt] = new TH2F(Form("JetLeadingDeltaPhi%s",lastnamehist.Data()),Form("#phi_{jet} - #phi_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 737         fhJetLeadingDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
 738         fhJetLeadingDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 739
 740         fhJetLeadingDeltaEtas[icone][ipt] = new TH2F(Form("JetLeadingDeltaEta%s",lastnamehist.Data()),Form("#eta_{jet} - #eta_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2);
 741         fhJetLeadingDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
 742         fhJetLeadingDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 743
 744         fhJetFFzs[icone][ipt] = new TH2F(Form("JetFFz%s",lastnamehist.Data()),"z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2);
 745         fhJetFFzs[icone][ipt]->SetYTitle("z");
 746         fhJetFFzs[icone][ipt]->SetXTitle("p_{T trigger}");
 747
 748         fhJetFFxis[icone][ipt] = new TH2F(Form("JetFFxi%s",lastnamehist.Data()),"#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.);
 749         fhJetFFxis[icone][ipt]->SetYTitle("#xi");
 750         fhJetFFxis[icone][ipt]->SetXTitle("p_{T trigger}");
 751
 752         fhJetFFpts[icone][ipt] = new TH2F(Form("JetFFpt%s",lastnamehist.Data()),"p_{T charged hadron } in jet vs p_{T trigger}", 120,0.,120.,200,0.,50.);
 753         fhJetFFpts[icone][ipt]->SetYTitle("p_{T charged hadron}");
 754         fhJetFFpts[icone][ipt]->SetXTitle("p_{T trigger}");
 755
 756         fhJetNTracksInCones[icone][ipt] = new TH2F(Form("JetNTracksInCone%s",lastnamehist.Data()),Form("N particles in cone vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,5000,0, 5000);
 757         fhJetNTracksInCones[icone][ipt]->SetYTitle("N tracks in jet cone");
 758         fhJetNTracksInCones[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 759
 760         fOutCont->Add(fhJetPts[icone][ipt]) ;
 761         fOutCont->Add(fhJetRatioPts[icone][ipt]) ;
 762         fOutCont->Add(fhJetDeltaPhis[icone][ipt]) ;
 763         fOutCont->Add(fhJetDeltaEtas[icone][ipt]) ;
 764         fOutCont->Add(fhJetLeadingRatioPts[icone][ipt]) ;
 765         fOutCont->Add(fhJetLeadingDeltaPhis[icone][ipt]) ;
 766         fOutCont->Add(fhJetLeadingDeltaEtas[icone][ipt]) ;
 767         fOutCont->Add(fhJetFFzs[icone][ipt]) ;
 768         fOutCont->Add(fhJetFFxis[icone][ipt]) ;
 769         fOutCont->Add(fhJetFFpts[icone][ipt]) ;
 770         fOutCont->Add(fhJetNTracksInCones[icone][ipt]) ;
 771
 772         //Bkg Distributions
 773         fhBkgPts[icone][ipt] = new TH2F(Form("BkgPt%s",lastnamehist.Data()),Form("p_{T  bkg} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 774         fhBkgPts[icone][ipt]->SetYTitle("p_{T  bkg}");
 775         fhBkgPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 776
 777         fhBkgRatioPts[icone][ipt] = new TH2F(Form("BkgRatioPt%s",lastnamehist.Data()),Form("p_{T  bkg}/p_{T trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2);
 778         fhBkgRatioPts[icone][ipt]->SetYTitle("p_{T  bkg}/p_{T trigger}");
 779         fhBkgRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 780
 781         fhBkgDeltaPhis[icone][ipt] = new TH2F(Form("BkgDeltaPhi%s",lastnamehist.Data()),Form("#phi_{bkg} - #phi_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 782         fhBkgDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
 783         fhBkgDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 784
 785         fhBkgDeltaEtas[icone][ipt] = new TH2F(Form("BkgDeltaEta%s",lastnamehist.Data()),Form("#eta_{bkg} - #eta_{trigger} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2);
 786         fhBkgDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
 787         fhBkgDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 788
 789         fhBkgLeadingRatioPts[icone][ipt] = new TH2F(Form("BkgLeadingRatioPt%s",lastnamehist.Data()),Form("p_{T  bkg} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,2);
 790         fhBkgLeadingRatioPts[icone][ipt]->SetYTitle("p_{T  leading}/p_{T bkg}");
 791         fhBkgLeadingRatioPts[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 792
 793         fhBkgLeadingDeltaPhis[icone][ipt] = new TH2F(Form("BkgLeadingDeltaPhi%s",lastnamehist.Data()),Form("#phi_{bkg} - #phi_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,0,TMath::TwoPi());
 794         fhBkgLeadingDeltaPhis[icone][ipt]->SetYTitle("#Delta #phi (rad)");
 795         fhBkgLeadingDeltaPhis[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 796
 797         fhBkgLeadingDeltaEtas[icone][ipt] = new TH2F(Form("BkgLeadingDeltaEta%s",lastnamehist.Data()),Form("#eta_{bkg} - #eta_{leading} vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,120,-2,2);
 798         fhBkgLeadingDeltaEtas[icone][ipt]->SetYTitle("#Delta #eta");
 799         fhBkgLeadingDeltaEtas[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 800
 801         fhBkgFFzs[icone][ipt] = new TH2F(Form("BkgFFz%s",lastnamehist.Data()),"z = p_{T i charged}/p_{T trigger} vs p_{T trigger}", 120,0.,120.,200,0.,2);
 802         fhBkgFFzs[icone][ipt]->SetYTitle("z");
 803         fhBkgFFzs[icone][ipt]->SetXTitle("p_{T trigger}");
 804
 805         fhBkgFFxis[icone][ipt] = new TH2F(Form("BkgFFxi%s",lastnamehist.Data()),"#xi = ln(p_{T trigger}/p_{T i charged}) vs p_{T trigger}", 120,0.,120.,100,0.,10.);
 806         fhBkgFFxis[icone][ipt]->SetYTitle("#xi");
 807         fhBkgFFxis[icone][ipt]->SetXTitle("p_{T trigger}");
 808
 809         fhBkgFFpts[icone][ipt] = new TH2F(Form("BkgFFpt%s",lastnamehist.Data()),"p_{T charged hadron} in jet vs p_{T trigger}", 120,0.,120.,200,0.,50.);
 810         fhBkgFFpts[icone][ipt]->SetYTitle("p_{T charged hadron}");
 811         fhBkgFFpts[icone][ipt]->SetXTitle("p_{T trigger}");
 812
 813         fhBkgNTracksInCones[icone][ipt] = new TH2F(Form("BkgNTracksInCone%s",lastnamehist.Data()),Form("N particles in cone vs p_{T trigger} %s",lastnametitle.Data()),nptbins,ptmin,ptmax,5000,0, 5000);
 814         fhBkgNTracksInCones[icone][ipt]->SetYTitle("N tracks in bkg cone");
 815         fhBkgNTracksInCones[icone][ipt]->SetXTitle("p_{T trigger} (GeV/c)");
 816
 817         fOutCont->Add(fhBkgPts[icone][ipt]) ;
 818         fOutCont->Add(fhBkgRatioPts[icone][ipt]) ;
 819         fOutCont->Add(fhBkgDeltaPhis[icone][ipt]) ;
 820         fOutCont->Add(fhBkgDeltaEtas[icone][ipt]) ;
 821         fOutCont->Add(fhBkgLeadingRatioPts[icone][ipt]) ;
 822         fOutCont->Add(fhBkgLeadingDeltaPhis[icone][ipt]) ;
 823         fOutCont->Add(fhBkgLeadingDeltaEtas[icone][ipt]) ;
 824         fOutCont->Add(fhBkgFFzs[icone][ipt]) ;
 825         fOutCont->Add(fhBkgFFxis[icone][ipt]) ;
 826         fOutCont->Add(fhBkgFFpts[icone][ipt]) ;
 827         fOutCont->Add(fhBkgNTracksInCones[icone][ipt]) ;
 828
 829       }//ipt
 830     } //icone
 831   }//If we want to study any cone or pt threshold
 832
 833   //Keep neutral meson selection histograms if requiered
 834   //Setting done in AliNeutralMesonSelection
 835   if(GetNeutralMesonSelection()){
 836     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
 837     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
 838       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) fOutCont->Add(nmsHistos->At(i)) ;
 839         delete nmsHistos;
 840   }
 841
 842
 843   if(GetDebug() > 2){
 844     printf("AliAnaParticleJetLeadingConeCorrelation::GetCreateOutputObjects() - All histograms names : \n");
 845     for(Int_t i  = 0 ;  i<  fOutCont->GetEntries(); i++)
 846       printf("Histo i %d name %s\n",i,((fOutCont->At(i))->GetName()));
 847     //cout<< (fOutCont->At(i))->GetName()<<endl;
 848   }
 849
 850   return fOutCont;
 851
 852 }
 853
 854 //____________________________________________________________________________
 855 Bool_t  AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle(AliAODPWG4ParticleCorrelation *particle, TLorentzVector & pLeading)
 856 {
 857   //Search Charged or Neutral leading particle, select the highest one and fill AOD
 858
 859   TLorentzVector pLeadingCh(0,0,0,0) ;
 860   TLorentzVector pLeadingPi0(0,0,0,0) ;
 861
 862   GetLeadingCharge(particle, pLeadingCh) ;
 863   if(!fJetsOnlyInCTS) GetLeadingPi0(particle, pLeadingPi0) ;
 864
 865   Double_t ptch = pLeadingCh.Pt();
 866   Double_t ptpi = pLeadingPi0.Pt();
 867   if (ptch > 0 || ptpi > 0){
 868     if((ptch >= ptpi)){
 869       if(GetDebug() > 1)printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Leading found in CTS \n");
 870       pLeading = pLeadingCh;
 871       if(GetDebug() > 1) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Found Leading: pt %2.3f, phi %2.3f deg, eta %2.3f\n",
 872                                 pLeading.Pt(),pLeading.Phi()*TMath::RadToDeg(),pLeading.Eta()) ;
 873       //Put leading in AOD
 874       particle->SetLeading(pLeadingCh);
 875       particle->SetLeadingDetector("CTS");
 876       return kTRUE;
 877     }
 878     else{
 879       if(GetDebug() > 1)
 880         printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Leading found in EMCAL \n");
 881       pLeading = pLeadingPi0;
 882       if(GetDebug() > 1) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - Found Leading: pt %2.3f, phi %2.3f, eta %2.3f\n",
 883                                 pLeading.Pt(),pLeading.Phi()*TMath::RadToDeg(),pLeading.Eta()) ;
 884       //Put leading in AOD
 885       particle->SetLeading(pLeadingPi0);
 886       particle->SetLeadingDetector("EMCAL");
 887       return kTRUE;
 888     }
 889   }
 890
 891   if(GetDebug() > 1)printf ("AliAnaParticleJetLeadingConeCorrelation::GetLeadingParticle() - NO LEADING PARTICLE FOUND \n");
 892
 893   return kFALSE;
 894
 895 }
 896
 897 //____________________________________________________________________________
 898 void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge(AliAODPWG4ParticleCorrelation* const particle, TLorentzVector & pLeading) const
 899 {
 900   //Search for the charged particle with highest pt and with
 901   //Phi=Phi_trigger-Pi and pT=0.1E_gamma
 902
 903   if(GetCTSTracks()){
 904     Double_t ptTrig  = particle->Pt();
 905     Double_t phiTrig = particle->Phi();
 906     Double_t rat     = -100 ;
 907     Double_t ptl     = -100 ;
 908     Double_t phil    = -100 ;
 909     Double_t pt      = -100.;
 910     Double_t phi     = -100.;
 911     TVector3 p3;
 912
 913     for(Int_t ipr = 0;ipr < GetCTSTracks()->GetEntriesFast() ; ipr ++ ){
 914       AliVTrack* track = (AliVTrack *)(GetCTSTracks()->At(ipr)) ;
 915       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
 916       pt   = p3.Pt();
 917       phi  = p3.Phi() ;
 918       if(phi < 0) phi+=TMath::TwoPi();
 919       rat  = pt/ptTrig ;
 920       //printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge() - Tracks: pt %2.3f eta %2.3f phi %2.3f pt/ptTrig %2.3f \n",
 921       //           pt, p3.Eta(), phi,pt/ptTrig) ;
 922       Float_t deltaphi = TMath::Abs(phiTrig-phi);
 923       if((deltaphi > fDeltaPhiMinCut) && (deltaphi < fDeltaPhiMaxCut) &&
 924          (rat > fLeadingRatioMinCut) && (rat < fLeadingRatioMaxCut)  && (pt  > ptl)) {
 925         phil = phi ;
 926         ptl  = pt ;
 927         pLeading.SetVect(p3);
 928       }
 929     }// track loop
 930
 931     if(GetDebug() > 1 && ptl > 0 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingCharge() - Leading in CTS: pt %2.3f eta %2.3f phi %2.3f pt/ptTrig %2.3f, |phiTrig-phi| %2.3f \n",
 932                                           ptl, pLeading.Eta(), phil,ptl/ptTrig, TMath::Abs(phiTrig-phil)) ;
 933
 934   }//CTS list exist
 935 }
 936
 937 //____________________________________________________________________________
 938 void  AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0(AliAODPWG4ParticleCorrelation* const particle, TLorentzVector & pLeading)
 939 {
 940   //Search for the neutral pion with highest pt and with
 941   //Phi=Phi_trigger-Pi and pT=0.1E_gamma
 942
 943   if(GetEMCALClusters()){
 944     Double_t ptTrig  = particle->Pt();
 945     Double_t phiTrig = particle->Phi();
 946     Double_t rat     = -100 ;
 947     Double_t ptl     = -100 ;
 948     Double_t phil    = -100 ;
 949     Double_t pt      = -100.;
 950     Double_t phi     = -100.;
 951
 952     TLorentzVector gammai;
 953     TLorentzVector gammaj;
 954
 955     //Get vertex for photon momentum calculation
 956     Double_t vertex [] = {0,0,0} ; //vertex
 957     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 958     {
 959       GetVertex(vertex);
 960     }
 961
 962     //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
 963     for(Int_t iclus = 0;iclus < GetEMCALClusters()->GetEntriesFast() ; iclus ++ ){
 964       AliVCluster * calo = (AliVCluster *)(GetEMCALClusters()->At(iclus)) ;
 965
 966       //Cluster selection, not charged, with photon or pi0 id and in fiducial cut
 967       Int_t pdgi=0;
 968       if(!SelectCluster(calo, vertex,  gammai, pdgi))  continue ;
 969
 970       if(GetDebug() > 2) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral cluster: pt %2.3f, phi %2.3f \n",
 971                                 gammai.Pt(),gammai.Phi());
 972
 973       //2 gamma overlapped, found with PID
 974       if(pdgi == AliCaloPID::kPi0){
 975
 976         if(GetDebug() > 2) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral cluster ID as Pi0 \n");
 977
 978         pt  = gammai.Pt();
 979         rat = pt/ptTrig;
 980         phi = gammai.Phi();
 981         if(phi < 0) phi+=TMath::TwoPi();
 982
 983         //Selection within angular and energy limits
 984         Float_t deltaphi = TMath::Abs(phiTrig-phi);
 985         if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  &&
 986            deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut )
 987         {
 988           phil = phi ;
 989           ptl  = pt ;
 990           pLeading.SetPxPyPzE(gammai.Px(),gammai.Py(),gammai.Pz(),gammai.E());
 991         }// cuts
 992       }// pdg = AliCaloPID::kPi0
 993       //Make invariant mass analysis
 994       else if(pdgi == AliCaloPID::kPhoton){
 995         //Search the photon companion in case it comes from  a Pi0 decay
 996         //Apply several cuts to select the good pair
 997         for(Int_t jclus = iclus+1; jclus < GetEMCALClusters()->GetEntriesFast() ; jclus ++ ){
 998           AliVCluster * calo2 = (AliVCluster *) (GetEMCALClusters()->At(jclus)) ;
 999
1000           //Cluster selection, not charged with photon or pi0 id and in fiducial cut
1001           Int_t pdgj=0;
1002
1003           if     (!SelectCluster(calo2, vertex,  gammaj, pdgj))  continue ;
1004
1005           if(pdgj == AliCaloPID::kPhoton ){
1006
1007             pt  = (gammai+gammaj).Pt();
1008             phi = (gammai+gammaj).Phi();
1009             if(phi < 0) phi+=TMath::TwoPi();
1010             rat = pt/ptTrig;
1011
1012             //Selection within angular and energy limits
1013             Float_t deltaphi = TMath::Abs(phiTrig-phi);
1014             if(GetDebug() > 3 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral Hadron Correlation: gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, |phiTrig-phi| %2.3f, pt/ptTrig %2.3f, M %2.3f\n",
1015                                        pt,phi,(gammai+gammaj).Eta(), deltaphi, rat, (gammai+gammaj).M());
1016
1017             if(pt > ptl  && rat > fLeadingRatioMinCut  && rat < fLeadingRatioMaxCut  &&
1018                deltaphi > fDeltaPhiMinCut && deltaphi < fDeltaPhiMaxCut ){
1019               //Select good pair (aperture and invariant mass)
1020               if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj,"EMCAL")){
1021                 phil = phi ;
1022                 ptl  = pt ;
1023                 pLeading=(gammai+gammaj);
1024
1025                 if(GetDebug() > 3 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Neutral Hadron Correlation: Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f\n",
1026                                            ptl,phil,(gammai+gammaj).Eta(), (gammai+gammaj).M());
1027               }//pi0 selection
1028
1029
1030             }//Pair selected as leading
1031           }//if pair of gammas
1032         }//2nd loop
1033       }// if pdg = 22
1034     }// 1st Loop
1035
1036     if(GetDebug() > 2 && pLeading.Pt() > 0 ) printf("AliAnaParticleJetLeadingConeCorrelation::GetLeadingPi0() - Leading EMCAL: pt %2.3f eta %2.3f phi %2.3f pt/Eg %2.3f \n",
1037                                                     pLeading.Pt(), pLeading.Eta(), phil,  pLeading.Pt()/ptTrig) ;
1038
1039   }//EMCAL list exists
1040 }
1041
1042 //____________________________________________________________________________
1043 void AliAnaParticleJetLeadingConeCorrelation::InitParameters()
1044 {
1045   //Initialize the parameters of the analysis.
1046   SetInputAODName("PWG4Particle");
1047   SetAODObjArrayName("JetLeadingCone");
1048   AddToHistogramsName("AnaJetLCCorr_");
1049
1050   fJetsOnlyInCTS      = kFALSE ;
1051   fPbPb               = kFALSE ;
1052   fReMakeJet          = kFALSE ;
1053
1054   //Leading selection parameters
1055   fDeltaPhiMinCut     = 2.9 ;
1056   fDeltaPhiMaxCut     = 3.4 ;
1057   fLeadingRatioMinCut = 0.1;
1058   fLeadingRatioMaxCut = 1.5;
1059
1060   //Jet selection parameters
1061   //Fixed cut
1062   fJetRatioMaxCut     = 1.2 ;
1063   fJetRatioMinCut     = 0.3 ;
1064   fJetCTSRatioMaxCut  = 1.2 ;
1065   fJetCTSRatioMinCut  = 0.3 ;
1066   fSelect               = 0  ; //0, Accept all jets, 1, selection depends on energy, 2 fixed selection
1067
1068   fSelectIsolated = kFALSE;
1069
1070   //Cut depending on gamma energy
1071   fPtTriggerSelectionCut = 10.; //For Low pt jets+BKG, another limits applied
1072   //Reconstructed jet energy dependence parameters
1073   //e_jet = a1+e_gamma b2.
1074   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1075   fJetE1[0] = -5.75; fJetE1[1] = -4.1;
1076   fJetE2[0] = 1.005; fJetE2[1] = 1.05;
1077
1078   //Reconstructed sigma of jet energy dependence parameters
1079   //s_jet = a1+e_gamma b2.
1080   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1081   fJetSigma1[0] = 2.65;   fJetSigma1[1] = 2.75;
1082   fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
1083
1084   //Background mean energy and RMS
1085   //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1086   //Index 2-> (low pt jets)BKG > 0.5 GeV;
1087   //Index > 2, same for CTS conf
1088   fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
1089   fBkgMean[3] = 0.; fBkgMean[4] = 6.4;  fBkgMean[5] = 48.6;
1090   fBkgRMS[0]  = 0.; fBkgRMS[1]  = 7.5;  fBkgRMS[2]  = 22.0;
1091   fBkgRMS[3]  = 0.; fBkgRMS[4]  = 5.4;  fBkgRMS[5]  = 13.2;
1092
1093   //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
1094   //limits for monoenergetic jets.
1095   //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1096   //Index 2-> (low pt jets) BKG > 0.5 GeV;
1097   //Index > 2, same for CTS conf
1098
1099   fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
1100   fJetXMin1[3] =-2.0  ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1  ;
1101   fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
1102   fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
1103   fJetXMax1[0] =-3.8  ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6  ;
1104   fJetXMax1[3] =-2.7  ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7  ;
1105   fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
1106   fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
1107
1108
1109   //Different cones and pt thresholds to construct the jet
1110
1111   fJetCone        = 0.3  ;
1112   fJetPtThreshold = 0.5   ;
1113   fJetPtThresPbPb = 2.   ;
1114   fJetNCone       = 4    ;
1115   fJetNPt         = 4    ;
1116   fJetCones[0]    = 0.2  ; fJetNameCones[0]   = "02" ;
1117   fJetCones[1]    = 0.3  ; fJetNameCones[1]   = "03" ;
1118   fJetCones[2]    = 0.4  ; fJetNameCones[2]   = "04" ;
1119   fJetCones[2]    = 0.5  ; fJetNameCones[2]   = "05" ;
1120
1121   fJetPtThres[0]  = 0.0  ; fJetNamePtThres[0] = "00" ;
1122   fJetPtThres[1]  = 0.5  ; fJetNamePtThres[1] = "05" ;
1123   fJetPtThres[2]  = 1.0  ; fJetNamePtThres[2] = "10" ;
1124   fJetPtThres[3]  = 2.0  ; fJetNamePtThres[3] = "20" ;
1125 }
1126
1127 //__________________________________________________________________________-
1128 Bool_t AliAnaParticleJetLeadingConeCorrelation::IsJetSelected(const Double_t ptTrig, const Double_t ptjet) const {
1129   //Given the pt of the jet and the trigger particle, select the jet or not
1130   //3 options, fSelect=0 accepts all, fSelect=1 selects jets depending on a
1131   //function energy dependent and fSelect=2 selects on simple fixed cuts
1132
1133   if(ptjet == 0) return kFALSE;
1134
1135   Double_t rat = ptTrig / ptjet ;
1136
1137   //###############################################################
1138   if(fSelect == 0)
1139     return kTRUE; //Accept all jets, no restriction
1140   //###############################################################
1141   else if(fSelect == 1){
1142     //Check if the energy of the reconstructed jet is within an energy window
1143     //WARNING: to be rechecked, don't remember what all the steps mean
1144     Double_t par[6];
1145     Double_t xmax[2];
1146     Double_t xmin[2];
1147
1148     Int_t iCTS = 0;
1149     if(fJetsOnlyInCTS)
1150       iCTS = 3 ;
1151
1152     if(!fPbPb){
1153       //Phythia alone, jets with pt_th > 0.2, r = 0.3
1154       par[0] = fJetE1[0]; par[1] = fJetE2[0];
1155       //Energy of the jet peak
1156       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1157       par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1158       //Sigma  of the jet peak
1159       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1160       par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
1161       //Parameters reserved for PbPb bkg.
1162       xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
1163       xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
1164       //Factor that multiplies sigma to obtain the best limits,
1165       //by observation, of mono jet ratios (ptjet/ptTrig)
1166       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1167
1168     }
1169     else{
1170       if(ptTrig > fPtTriggerSelectionCut){
1171         //Phythia +PbPb with  pt_th > 2 GeV/c, r = 0.3
1172         par[0] = fJetE1[0]; par[1] = fJetE2[0];
1173         //Energy of the jet peak, same as in pp
1174         //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1175         par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1176         //Sigma  of the jet peak, same as in pp
1177         //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1178         par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
1179         //Mean value and RMS of PbPb Bkg
1180         xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
1181         xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
1182         //Factor that multiplies sigma to obtain the best limits,
1183         //by observation, of mono jet ratios (ptjet/ptTrig) mixed with PbPb Bkg,
1184         //pt_th > 2 GeV, r = 0.3
1185         //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1186
1187       }
1188       else{
1189         //Phythia + PbPb with  pt_th > 0.5 GeV/c, r = 0.3
1190         par[0] = fJetE1[1]; par[1] = fJetE2[1];
1191         //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
1192         //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1193         par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
1194         //Sigma  of the jet peak, pt_th > 2 GeV/c, r = 0.3
1195         //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1196         par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
1197         //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
1198         xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
1199         xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
1200         //Factor that multiplies sigma to obtain the best limits,
1201         //by observation, of mono jet ratios (ptjet/ptTrig) mixed with PbPb Bkg,
1202         //pt_th > 2 GeV, r = 0.3
1203         //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1204
1205       }//If low pt jet in bkg
1206     }//if Bkg
1207
1208     //Calculate minimum and maximum limits of the jet ratio.
1209     Double_t min = CalculateJetRatioLimit(ptTrig, par, xmin);
1210     Double_t max = CalculateJetRatioLimit(ptTrig, par, xmax);
1211
1212     if(GetDebug() > 3)printf("AliAnaParticleJetLeadingConeCorrelation::IsJetSelected() - Jet selection?  : Limits min %2.3f, max %2.3f,  pt_jet %2.3f,  pt_gamma %2.3f, pt_jet / pt_gamma %2.3f\n",min,max,ptjet,ptTrig,rat);
1213
1214     if(( min < rat ) && ( max > ptjet/rat))
1215       return kTRUE;
1216     else
1217       return kFALSE;
1218   }//fSelect = 1
1219   //###############################################################
1220   else if(fSelect == 2){
1221     //Simple selection
1222     if(!fJetsOnlyInCTS){
1223       if((rat <  fJetRatioMaxCut) && (rat > fJetRatioMinCut )) return kTRUE;
1224     }
1225     else{
1226       if((rat <  fJetCTSRatioMaxCut) && (rat > fJetCTSRatioMinCut )) return kTRUE;
1227     }
1228   }// fSelect = 2
1229   //###############################################################
1230   else{
1231     printf("AliAnaParticleJetLeadingConeCorrelation::IsJetSelected() - Jet selection option larger than 2, DON'T SELECT JETS\n");
1232     return kFALSE ;
1233   }
1234
1235   return kFALSE;
1236
1237 }
1238
1239 //___________________________________________________________________
1240 Bool_t AliAnaParticleJetLeadingConeCorrelation::IsParticleInJetCone(const Double_t eta, Double_t phi, const Double_t etal,  Double_t phil)
1241   const {
1242   //Check if the particle is inside the cone defined by the leading particle
1243   //WARNING: To be rechecked
1244
1245   if(phi < 0) phi+=TMath::TwoPi();
1246   if(phil < 0) phil+=TMath::TwoPi();
1247   Double_t  rad = 10000 + fJetCone;
1248
1249   if(TMath::Abs(phi-phil) <= (TMath::TwoPi() - fJetCone))
1250     rad = TMath::Sqrt(TMath::Power(eta-etal,2)+TMath::Power(phi-phil,2));
1251   else{
1252     if(phi-phil > TMath::TwoPi() - fJetCone)
1253       rad = TMath::Sqrt(TMath::Power(eta-etal,2)+TMath::Power((phi-TMath::TwoPi())-phil,2));
1254     if(phi-phil < -(TMath::TwoPi() - fJetCone))
1255       rad = TMath::Sqrt(TMath::Power(eta-etal,2)+TMath::Power((phi+TMath::TwoPi())-phil,2));
1256   }
1257
1258   if(rad < fJetCone) return kTRUE ;
1259   else return kFALSE ;
1260
1261 }
1262
1263 //__________________________________________________________________
1264 void  AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD()
1265 {
1266   //Particle-Hadron Correlation Analysis, fill AODs
1267
1268   if(!GetInputAODBranch()){
1269     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - No input particles in AOD with name branch < %s > \n",
1270            GetInputAODName().Data());
1271     abort();
1272   }
1273
1274   if(strcmp(GetInputAODBranch()->GetClass()->GetName(), "AliAODPWG4ParticleCorrelation")){
1275     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - Wrong type of AOD object, change AOD class name in input AOD: It should be <AliAODPWG4ParticleCorrelation> and not <%s> \n",GetInputAODBranch()->GetClass()->GetName());
1276     abort();
1277   }
1278
1279   if(GetDebug() > 1){
1280     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - Begin jet leading cone  correlation analysis, fill AODs \n");
1281     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
1282     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetCTSTracks()->GetEntriesFast());
1283     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
1284   }
1285
1286   TLorentzVector pLeading(0,0,0,0); //It will contain the kinematics of the found leading particle
1287
1288   //Loop on stored AOD particles, trigger
1289   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1290   for(Int_t iaod = 0; iaod < naod ; iaod++){
1291     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1292
1293     //  printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - Trigger : pt %3.2f, phi %2.2f, eta %2.2f\n",particle->Pt(), particle->Phi(), particle->Eta());
1294
1295     //Search leading particles in CTS and EMCAL
1296     if(GetLeadingParticle(particle, pLeading)){
1297
1298       //Construct the jet around the leading, Fill AOD jet particle list, select jet
1299       //and fill AOD with jet and background
1300       MakeAODJet(particle, pLeading);
1301
1302     }//Leading found
1303   }//AOD trigger particle loop
1304
1305   if(GetDebug() >1)printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillAOD() - End of jet leading cone analysis, fill AODs \n");
1306
1307 }
1308
1309 //__________________________________________________________________
1310 void  AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms()
1311 {
1312
1313   //Particle-Hadron Correlation Analysis, fill histograms
1314
1315   if(!GetInputAODBranch()){
1316     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - No input particles in AOD with name branch < %s > \n",
1317            GetInputAODName().Data());
1318     abort();
1319   }
1320   if(GetDebug() > 1){
1321     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - Begin jet leading cone  correlation analysis, fill histograms \n");
1322     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - In particle branch aod entries %d\n", GetInputAODBranch()->GetEntriesFast());
1323     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - In CTS aod entries %d\n", GetCTSTracks()->GetEntriesFast());
1324     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
1325   }
1326
1327   TLorentzVector pLeading(0,0,0,0) ;
1328
1329   //Loop on stored AOD particles, trigger
1330   Int_t naod = GetInputAODBranch()->GetEntriesFast();
1331   for(Int_t iaod = 0; iaod < naod ; iaod++){
1332     AliAODPWG4ParticleCorrelation* particle =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
1333
1334     if(OnlyIsolated() && !particle->IsIsolated()) continue;
1335
1336     Double_t pt = particle->Pt();
1337     Double_t phi = particle->Phi();
1338     Double_t eta = particle->Eta();
1339
1340     //Get leading particle, fill histograms
1341     pLeading = particle->GetLeading();
1342     TString det = particle->GetLeadingDetector();
1343
1344     if(det!="" && pLeading.Pt() > 0){
1345       Double_t ptL = pLeading.Pt();
1346       Double_t phiL = pLeading.Phi();
1347       if(phiL < 0 ) phiL+=TMath::TwoPi();
1348       Double_t etaL = pLeading.Eta();
1349
1350       if(GetDebug() > 1)   printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - Trigger with pt %3.2f, phi %2.2f, eta %2.2f\n", pt, phi, eta);
1351
1352       if(det == "CTS"){
1353         fhChargedLeadingPt->Fill(pt,ptL);
1354         fhChargedLeadingPhi->Fill(pt,phiL);
1355         fhChargedLeadingEta->Fill(pt,etaL);
1356         fhChargedLeadingDeltaPt->Fill(pt,pt-ptL);
1357         fhChargedLeadingDeltaPhi->Fill(pt,TMath::Abs(phi-phiL));
1358         fhChargedLeadingDeltaEta->Fill(pt,eta-etaL);
1359         fhChargedLeadingRatioPt->Fill(pt,ptL/pt);
1360         fhChargedLeadingXi->Fill(pt,TMath::Log(pt/ptL));
1361         if(pt > 30) fhChargedLeadingDeltaPhiRatioPt30->Fill(TMath::Abs(phi-phiL),ptL/pt);
1362         if(pt > 50) fhChargedLeadingDeltaPhiRatioPt50->Fill(TMath::Abs(phi-phiL),ptL/pt);
1363       }
1364       else if(det== "EMCAL"){
1365         fhNeutralLeadingPt->Fill(pt,ptL);
1366         fhNeutralLeadingPhi->Fill(pt,phiL);
1367         fhNeutralLeadingEta->Fill(pt,etaL);
1368         fhNeutralLeadingDeltaPt->Fill(pt,pt-ptL);
1369         fhNeutralLeadingDeltaPhi->Fill(pt,TMath::Abs(phi-phiL));
1370         fhNeutralLeadingDeltaEta->Fill(pt,eta-etaL);
1371         fhNeutralLeadingRatioPt->Fill(pt,ptL/pt);
1372         fhNeutralLeadingXi->Fill(pt,TMath::Log(pt/ptL));
1373         if(pt > 30) fhNeutralLeadingDeltaPhiRatioPt30->Fill(TMath::Abs(phi-phiL),ptL/pt);
1374         if(pt > 50) fhNeutralLeadingDeltaPhiRatioPt50->Fill(TMath::Abs(phi-phiL),ptL/pt);
1375
1376       }
1377
1378       //Fill Jet histograms
1379       TLorentzVector bkg(0,0,0,0);
1380       TLorentzVector jet(0,0,0,0);
1381       if(!fSeveralConeAndPtCuts){//just fill histograms
1382         if(!fReMakeJet){
1383           jet=particle->GetCorrelatedJet();
1384           bkg=particle->GetCorrelatedBackground();
1385         }
1386         else  MakeJetFromAOD(particle, pLeading, jet,bkg);
1387
1388         if(jet.Pt() > 0){//Jet was found
1389           FillJetHistos(particle, pLeading, jet,"Jet","");
1390           FillJetHistos(particle, pLeading, bkg,"Bkg","");
1391         }
1392       }
1393       else if(fSeveralConeAndPtCuts){
1394         for(Int_t icone = 0; icone<fJetNCone; icone++) {
1395           fJetCone=fJetCones[icone];
1396           for(Int_t ipt = 0; ipt<fJetNPt;ipt++) {
1397             TString lastname ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt];
1398             fJetPtThreshold=fJetPtThres[ipt];
1399             MakeJetFromAOD(particle, pLeading, jet,bkg);
1400             if(jet.Pt() > 0) {//Jet was found
1401               FillJetHistos(particle, pLeading, jet,"Jet",lastname);
1402               FillJetHistos(particle, pLeading, bkg,"Bkg",lastname);
1403             }
1404           }//icone
1405         }//ipt
1406       }//fSeveralConeAndPtCuts
1407     }//Leading
1408   }//AOD trigger particle loop
1409
1410   if(GetDebug() >1)printf("AliAnaParticleJetLeadingConeCorrelation::MakeAnalysisFillHistograms() - End of jet leading cone analysis, fill histograms \n");
1411
1412 }
1413
1414 //____________________________________________________________________________
1415 void AliAnaParticleJetLeadingConeCorrelation::MakeAODJet(AliAODPWG4ParticleCorrelation *particle, const TLorentzVector  pLeading)
1416   const {
1417   //Fill the jet with the particles around the leading particle with
1418   //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and
1419   //fill aod with found information
1420
1421   TLorentzVector bkg(0,0,0,0);
1422   TLorentzVector jet(0,0,0,0);
1423   TLorentzVector lv (0,0,0,0); //Temporal container for jet particles kinematics
1424
1425   Double_t ptTrig   = particle->Pt();
1426   Double_t phiTrig  = particle->Phi();
1427   Double_t phil     = pLeading.Phi();
1428   if(phil<0) phil+=TMath::TwoPi();
1429   Double_t etal     = pLeading.Eta();
1430
1431   //Different pt cut for jet particles in different collisions systems
1432   Float_t ptcut = fJetPtThreshold;
1433   if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
1434
1435   //Add charged particles to jet if they are in cone around the leading particle
1436   if(!GetCTSTracks()) {
1437     printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - Cannot construct jets without tracks, STOP analysis");
1438     return;
1439   }
1440
1441   //Fill jet with tracks
1442   TVector3 p3;
1443   //Initialize reference arrays that will contain jet and background tracks
1444   TObjArray * reftracks  = new TObjArray;
1445   TObjArray * reftracksbkg  = new TObjArray;
1446
1447   for(Int_t ipr = 0;ipr < (GetCTSTracks())->GetEntriesFast() ; ipr ++ ){
1448     AliVTrack* track = (AliVTrack *)((GetCTSTracks())->At(ipr)) ;
1449     p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1450
1451     //Particles in jet
1452     if(IsParticleInJetCone(p3.Eta(), p3.Phi(), etal, phil)){
1453
1454       reftracks->Add(track);
1455
1456       if(p3.Pt() > ptcut ){
1457         lv.SetVect(p3);
1458         jet+=lv;
1459       }
1460     }
1461
1462     //Background around (phi_gamma-pi, eta_leading)
1463     else if(IsParticleInJetCone(p3.Eta(),p3.Phi(),etal, phiTrig)) {
1464
1465       reftracksbkg->Add(track);
1466
1467       if(p3.Pt() > ptcut ){
1468         lv.SetVect(p3);
1469         bkg+=lv;
1470       }
1471     }
1472   }//Track loop
1473
1474   //Add referenced tracks to AOD
1475   if(reftracks->GetEntriesFast() > 0 ){
1476     reftracks->SetName(Form("%sTracks",GetAODObjArrayName().Data()));
1477     particle->AddObjArray(reftracks);
1478   }
1479   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No tracks in jet cone\n");
1480   if(reftracksbkg->GetEntriesFast() > 0 ){
1481     reftracksbkg->SetName(Form("%sTracksBkg",GetAODObjArrayName().Data()));
1482     particle->AddObjArray(reftracksbkg);
1483   }
1484   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No background tracks in jet cone\n");
1485
1486   //Add neutral particles to jet
1487   //Initialize reference arrays that will contain jet and background tracks
1488   TObjArray * refclusters  = new TObjArray;
1489   TObjArray * refclustersbkg  = new TObjArray;
1490   if(!fJetsOnlyInCTS && GetEMCALClusters()){
1491
1492         //Get vertex for photon momentum calculation
1493         Double_t vertex[]  = {0,0,0} ; //vertex
1494         if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
1495         {
1496                 GetReader()->GetVertex(vertex);
1497                 //if(GetReader()->GetSecondInputAODTree()) GetReader()->GetSecondInputAODVertex(vertex2);
1498         }
1499
1500     for(Int_t iclus = 0;iclus < (GetEMCALClusters())->GetEntriesFast() ; iclus ++ ){
1501       AliVCluster * calo = (AliVCluster *) (GetEMCALClusters()->At(iclus)) ;
1502
1503       //Cluster selection, not charged
1504       if(IsTrackMatched(calo,GetReader()->GetInputEvent())) continue ;
1505
1506       //Get Momentum vector,
1507       calo->GetMomentum(lv,vertex) ;//Assume that come from vertex in straight line
1508
1509       //Particles in jet
1510       if(IsParticleInJetCone(lv.Eta(),lv.Phi(), etal, phil)){
1511
1512         refclusters->Add(calo);
1513
1514         if(lv.Pt() > ptcut )  jet+=lv;
1515       }
1516       //Background around (phi_gamma-pi, eta_leading)
1517       else if(IsParticleInJetCone(lv.Eta(),lv.Phi(),etal, phiTrig)){
1518
1519
1520         refclustersbkg->Add(calo);
1521
1522         if(lv.Pt() > ptcut )  bkg+=lv;
1523       }
1524     }//cluster loop
1525   }//jets with neutral particles
1526
1527   //Add referenced clusters to AOD
1528   if(refclusters->GetEntriesFast() > 0 ){
1529     refclusters->SetName(Form("%sClusters",GetAODObjArrayName().Data()));
1530     particle->AddObjArray(refclusters);
1531   }
1532   else  if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No clusters in jet cone\n");
1533   if(refclustersbkg->GetEntriesFast() > 0 ){
1534     refclustersbkg->SetName(Form("%sClustersBkg",GetAODObjArrayName().Data()));
1535     particle->AddObjArray(refclustersbkg);
1536   }
1537   else if(GetDebug() > 2 ) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - No background clusters in jet cone\n");
1538
1539   //If there is any jet found, select after some criteria and
1540   //and fill AOD with corresponding TLorentzVector kinematics
1541   if(IsJetSelected(particle->Pt(), jet.Pt())) {
1542     particle->SetCorrelatedJet(jet);
1543     particle->SetCorrelatedBackground(bkg);
1544     if(GetDebug()>1) printf("AliAnaParticleJetLeadingConeCorrelation::MakeAODJet() - Found jet: Trigger  pt %2.3f, Jet pt %2.3f, Bkg pt %2.3f\n",ptTrig,jet.Pt(),bkg.Pt());
1545   }
1546
1547 }
1548
1549 //______________________________________________________________________________________________________
1550 void AliAnaParticleJetLeadingConeCorrelation::MakeJetFromAOD(AliAODPWG4ParticleCorrelation *particle,
1551                                                              const TLorentzVector  pLeading,
1552                                                              TLorentzVector & jet, TLorentzVector & bkg) const
1553 {
1554   //Fill the jet with the particles around the leading particle with
1555   //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and
1556   //fill aod tlorentzvectors with jet and bakcground found
1557
1558   TLorentzVector lv (0,0,0,0); //Temporal container for jet particles kinematics
1559
1560   Double_t ptTrig  = particle->Pt();
1561   Double_t phiTrig = particle->Phi();
1562   Double_t phil = pLeading.Phi();
1563   if(phil < 0) phil+=TMath::TwoPi();
1564   Double_t etal = pLeading.Eta();
1565
1566   TObjArray * refclusters    = particle->GetObjArray(Form("Clusters%s"   ,GetAODObjArrayName().Data()));
1567   TObjArray * reftracks      = particle->GetObjArray(Form("Tracks%s"     ,GetAODObjArrayName().Data()));
1568   TObjArray * refclustersbkg = particle->GetObjArray(Form("ClustersBkg%s",GetAODObjArrayName().Data()));
1569   TObjArray * reftracksbkg   = particle->GetObjArray(Form("TracksBkg%s"  ,GetAODObjArrayName().Data()));
1570
1571   //Different pt cut for jet particles in different collisions systems
1572   Float_t ptcut = fJetPtThreshold;
1573   if(fPbPb && !fSeveralConeAndPtCuts && ptTrig > fPtTriggerSelectionCut)  ptcut = fJetPtThresPbPb ;
1574
1575   //Fill jet with tracks
1576   //Particles in jet
1577   TVector3 p3;
1578   if(reftracks){
1579     for(Int_t ipr = 0;ipr < reftracks->GetEntriesFast() ; ipr ++ ){
1580       AliVTrack* track = (AliVTrack *) reftracks->At(ipr) ;
1581       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1582       Float_t phi = p3.Phi();
1583       if(phi < 0) phi+=TMath::TwoPi();
1584       if(p3.Pt() > ptcut && IsParticleInJetCone(p3.Eta(), phi, etal, phil) ){
1585         lv.SetVect(p3);
1586         jet+=lv;
1587       }
1588     }//jet Track loop
1589   }
1590   //Particles in background
1591   if(reftracksbkg){
1592     for(Int_t ipr = 0;ipr < reftracksbkg->GetEntriesFast() ; ipr ++ ){
1593       AliVTrack* track = (AliVTrack *) reftracksbkg->At(ipr) ;
1594       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
1595       if(p3.Pt() > ptcut && IsParticleInJetCone(p3.Eta(),p3.Phi(),etal, phiTrig) ) {
1596         lv.SetVect(p3);
1597         bkg+=lv;
1598       }
1599     }//background Track loop
1600   }
1601
1602   //Add neutral particles to jet
1603   if(!fJetsOnlyInCTS && refclusters){
1604
1605     //Get vertex for photon momentum calculation
1606     Double_t vertex[]  = {0,0,0} ; //vertex
1607     if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
1608     {
1609       GetReader()->GetVertex(vertex);
1610     }
1611
1612     //Loop on jet particles
1613     if(refclusters){
1614       for(Int_t iclus = 0;iclus < refclusters->GetEntriesFast() ; iclus ++ ){
1615         AliVCluster * calo = (AliVCluster *) refclusters->At(iclus) ;
1616
1617         calo->GetMomentum(lv,vertex) ;//Assume that come from vertex in straight line
1618
1619         if(lv.Pt() > ptcut && IsParticleInJetCone(lv.Eta(),lv.Phi(), etal, phil)) jet+=lv;
1620       }//jet cluster loop
1621     }
1622
1623     //Loop on background particles
1624     if(refclustersbkg){
1625       for(Int_t iclus = 0;iclus < refclustersbkg->GetEntriesFast() ; iclus ++ ){
1626         AliVCluster * calo = (AliVCluster *) refclustersbkg->At(iclus) ;
1627
1628         calo->GetMomentum(lv,vertex) ;//Assume that come from vertex in straight line
1629
1630         if( lv.Pt() > ptcut && IsParticleInJetCone(lv.Eta(),lv.Phi(),etal, phiTrig)) bkg+=lv;
1631       }//background cluster loop
1632     }
1633   }//clusters in jet
1634
1635   //If there is any jet found, leave jet and bkg as they are,
1636   //if not set them to 0.
1637   if(!IsJetSelected(particle->Pt(), jet.Pt())) {
1638     jet.SetPxPyPzE(0.,0.,0.,0.);
1639     bkg.SetPxPyPzE(0.,0.,0.,0.);
1640   }
1641   else
1642     if(GetDebug()>1) printf("AliAnaParticleJetLeadingConeCorrelation::MakeJetFromAOD() - Found jet: Trigger  pt %2.3f, Jet pt %2.3f, Bkg pt %2.3f\n",ptTrig,jet.Pt(),bkg.Pt());
1643
1644 }
1645
1646 //__________________________________________________________________
1647 void AliAnaParticleJetLeadingConeCorrelation::Print(const Option_t * opt) const
1648 {
1649
1650   //Print some relevant parameters set for the analysis
1651   if(! opt)
1652     return;
1653
1654   printf("**** Print  %s %s ****\n", GetName(), GetTitle() ) ;
1655   AliAnaPartCorrBaseClass::Print(" ");
1656
1657   if(fJetsOnlyInCTS)printf("Jets reconstructed in CTS \n");
1658   else printf("Jets reconstructed in CTS+EMCAL \n");
1659
1660   if(fPbPb) printf("PbPb events, pT cut in jet cone energy reconstruction %2.1f \n", fJetPtThreshold);
1661   else printf("pp events, pT cut in jet cone energy reconstruction %2.1f \n", fJetPtThresPbPb);
1662
1663   printf("If pT of trigger < %2.3f, select jets as in pp? \n", fPtTriggerSelectionCut);
1664
1665   printf("Phi gamma-Leading        <     %3.2f\n", fDeltaPhiMaxCut) ;
1666   printf("Phi gamma-Leading        >     %3.2f\n", fDeltaPhiMinCut) ;
1667   printf("pT Leading / pT Trigger  <     %3.2f\n", fLeadingRatioMaxCut) ;
1668   printf("pT Leading / pT Trigger  >     %3.2f\n", fLeadingRatioMinCut) ;
1669
1670   if(fSelect == 2){
1671     printf("pT Jet / pT Gamma                     <    %3.2f\n", fJetRatioMaxCut) ;
1672     printf("pT Jet / pT Gamma                     >    %3.2f\n", fJetRatioMinCut) ;
1673     printf("pT Jet (Only CTS)/ pT Trigger   <    %3.2f\n", fJetCTSRatioMaxCut) ;
1674     printf("pT Jet (Only CTS)/ pT Trigger   >    %3.2f\n", fJetCTSRatioMinCut) ;
1675   }
1676   else if(fSelect == 0)
1677     printf("Accept all reconstructed jets \n") ;
1678   else   if(fSelect == 1)
1679     printf("Accept jets depending on trigger energy \n") ;
1680   else
1681     printf("Wrong jet selection option:   %d \n", fSelect) ;
1682
1683   printf("Isolated Trigger?  %d\n", fSelectIsolated) ;
1684
1685 }