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