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