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