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