PWG4/AliAnaGammaJetLeadCone.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 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.3  2007/10/29 13:48:42  gustavo
  21  * Corrected coding violations
  22  *
  23  * Revision 1.2  2007/08/17 12:40:04  schutz
  24  * New analysis classes by Gustavo Conesa
  25  *
  26  * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
  27  * new analysis classes in the the new analysis framework
  28  *
  29  *
  30  */
  31
  32 //_________________________________________________________________________
  33 // Class for the analysis of gamma-jet correlations:
  34 // 1)Take the prompt photon found with AliAnaGammaDirect,
  35 // 2) Search for the highest pt leading particle opposite to the photon within a phi, pt window
  36 // 3) Take all particles around leading in a cone R with pt larger than threshold and construct the jet
  37 //
  38 //  Class created from old AliPHOSGammaJet
  39 //  (see AliRoot versions previous Release 4-09)
  40 //
  41 //*-- Author: Gustavo Conesa (LNF-INFN)
  42 //////////////////////////////////////////////////////////////////////////////
  43
  44
  45 // --- ROOT system ---
  46
  47 //---- AliRoot system ----
  48 #include "AliNeutralMesonSelection.h"
  49 #include "AliAnaGammaJetLeadCone.h"
  50 #include "AliLog.h"
  51
  52 ClassImp(AliAnaGammaJetLeadCone)
  53
  54
  55 //____________________________________________________________________________
  56   AliAnaGammaJetLeadCone::AliAnaGammaJetLeadCone() :
  57     AliAnaGammaCorrelation(),  fPbPb(kFALSE),
  58     fSeveralConeAndPtCuts(0),
  59     fEtaEMCALCut(0.),
  60     fJetCTSRatioMaxCut(0.), fJetCTSRatioMinCut(0.), fJetRatioMaxCut(0.),
  61     fJetRatioMinCut(0.),
  62     fJetNCone(0),fJetNPt(0), fJetCone(0),
  63     fJetPtThreshold(0),fJetPtThresPbPb(0),
  64     fPtJetSelectionCut(0.0), fSelect(0),
  65     fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0),
  66     fhDeltaPhiGammaCharged(0),  fhDeltaPhiGammaNeutral(0),
  67     fhDeltaEtaGammaCharged(0), fhDeltaEtaGammaNeutral(0),
  68     fhAnglePairLeading(), fhInvMassPairLeading(),
  69     fhChargedRatio(0), fhNeutralRatio (0),
  70     fhNBkg (0), fhNLeading(0), fhNJet(0), fhJetRatio(0), fhJetPt (0),
  71     fhBkgRatio (0), fhBkgPt(0),  fhJetFragment(0), fhBkgFragment(0),
  72     fhJetPtDist(0),  fhBkgPtDist(0)
  73 {
  74   //Default Ctor
  75
  76   //Initialize parameters
  77
  78   SetCorrelationType(kJetLeadCone);
  79
  80   for(Int_t i = 0; i<10; i++){
  81     fJetCones[i]         = 0.0 ;
  82     fJetNameCones[i]     = ""  ;
  83     fJetPtThres[i]      = 0.0 ;
  84     fJetNamePtThres[i]  = ""  ;
  85     if( i < 6 ){
  86       fJetXMin1[i]     = 0.0 ;
  87       fJetXMin2[i]     = 0.0 ;
  88       fJetXMax1[i]     = 0.0 ;
  89       fJetXMax2[i]     = 0.0 ;
  90       fBkgMean[i]      = 0.0 ;
  91       fBkgRMS[i]       = 0.0 ;
  92       if( i < 2 ){
  93         fJetE1[i]        = 0.0 ;
  94         fJetE2[i]        = 0.0 ;
  95         fJetSigma1[i]    = 0.0 ;
  96         fJetSigma2[i]    = 0.0 ;
  97         fPhiEMCALCut[i]  = 0.0 ;
  98       }
  99     }
 100   }
 101
 102   InitParameters();
 103 }
 104
 105 //____________________________________________________________________________
 106 AliAnaGammaJetLeadCone::AliAnaGammaJetLeadCone(const AliAnaGammaJetLeadCone & g) :
 107   AliAnaGammaCorrelation(g), fPbPb(g.fPbPb),
 108   fSeveralConeAndPtCuts(g.fSeveralConeAndPtCuts),
 109   fEtaEMCALCut(g.fEtaEMCALCut),
 110   fJetCTSRatioMaxCut(g.fJetCTSRatioMaxCut),
 111   fJetCTSRatioMinCut(g.fJetCTSRatioMinCut), fJetRatioMaxCut(g.fJetRatioMaxCut),
 112   fJetRatioMinCut(g.fJetRatioMinCut),  fJetNCone(g.fJetNCone),
 113   fJetNPt(g.fJetNPt), fJetCone(g.fJetCone),
 114   fJetPtThreshold(g.fJetPtThreshold),fJetPtThresPbPb(g.fJetPtThresPbPb),
 115   fPtJetSelectionCut(g.fPtJetSelectionCut), fSelect(g.fSelect),
 116   fhPhiCharged(g.fhPhiCharged), fhPhiNeutral(g.fhPhiNeutral),
 117   fhEtaCharged(g.fhEtaCharged), fhEtaNeutral(g.fhEtaNeutral),
 118   fhDeltaPhiGammaCharged(g.fhDeltaPhiGammaCharged),
 119   fhDeltaPhiGammaNeutral(g.fhDeltaPhiGammaNeutral),
 120   fhDeltaEtaGammaCharged(g.fhDeltaEtaGammaCharged),
 121   fhDeltaEtaGammaNeutral(g.fhDeltaEtaGammaNeutral),
 122   fhAnglePairLeading(g.fhAnglePairLeading), fhInvMassPairLeading(g.fhInvMassPairLeading),
 123   fhChargedRatio(g.fhChargedRatio), fhNeutralRatio(g.fhNeutralRatio),
 124   fhNBkg(g. fhNBkg), fhNLeading(g. fhNLeading), fhNJet(g.fhNJet), fhJetRatio(g.fhJetRatio), fhJetPt(g.fhJetPt),
 125   fhBkgRatio (g.fhBkgRatio), fhBkgPt(g.fhBkgPt),  fhJetFragment(g.fhJetFragment), fhBkgFragment(g.fhBkgFragment),
 126   fhJetPtDist(g.fhJetPtDist),  fhBkgPtDist(g.fhBkgPtDist)
 127 {
 128   // cpy ctor
 129
 130   for(Int_t i = 0; i<10; i++){
 131     fJetCones[i]        = g.fJetCones[i] ;
 132     fJetNameCones[i]    = g.fJetNameCones[i] ;
 133     fJetPtThres[i]      = g.fJetPtThres[i] ;
 134     fJetNamePtThres[i]  = g.fJetNamePtThres[i] ;
 135     if( i < 6 ){
 136       fJetXMin1[i]       = g.fJetXMin1[i] ;
 137       fJetXMin2[i]       = g.fJetXMin2[i] ;
 138       fJetXMax1[i]       = g.fJetXMax1[i] ;
 139       fJetXMax2[i]       = g.fJetXMax2[i] ;
 140       fBkgMean[i]        = g.fBkgMean[i] ;
 141       fBkgRMS[i]         = g.fBkgRMS[i] ;
 142       if( i < 2 ){
 143         fJetE1[i]        = g.fJetE1[i] ;
 144         fJetE2[i]        = g.fJetE2[i] ;
 145         fJetSigma1[i]    = g.fJetSigma1[i] ;
 146         fJetSigma2[i]    = g.fJetSigma2[i] ;
 147         fPhiEMCALCut[i]  = g.fPhiEMCALCut[i] ;
 148       }
 149     }
 150   }
 151
 152
 153 }
 154
 155 //_________________________________________________________________________
 156 AliAnaGammaJetLeadCone & AliAnaGammaJetLeadCone::operator = (const AliAnaGammaJetLeadCone & source)
 157 {
 158   // assignment operator
 159
 160   if(this == &source)return *this;
 161   ((AliAnaGammaCorrelation *)this)->operator=(source);
 162
 163   fSeveralConeAndPtCuts = source.fSeveralConeAndPtCuts ;
 164   fPbPb = source.fPbPb ;
 165   fEtaEMCALCut = source.fEtaEMCALCut ;
 166   fJetCTSRatioMaxCut = source.fJetCTSRatioMaxCut ;
 167   fJetCTSRatioMinCut = source.fJetCTSRatioMinCut ; fJetRatioMaxCut = source.fJetRatioMaxCut ;
 168   fJetRatioMinCut = source.fJetRatioMinCut ;  fJetNCone = source.fJetNCone ;
 169   fJetNPt = source.fJetNPt ; fJetCone = source.fJetCone ;
 170   fJetPtThreshold = source.fJetPtThreshold ;
 171   fJetPtThresPbPb = source.fJetPtThresPbPb ;
 172   fPtJetSelectionCut = source.fPtJetSelectionCut ;
 173   fSelect = source.fSelect ;  fhChargedRatio = source.fhChargedRatio ; fhNeutralRatio = source.fhNeutralRatio ;
 174
 175   fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ;
 176   fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ;
 177   fhDeltaPhiGammaCharged = source.fhDeltaPhiGammaCharged ;
 178   fhDeltaPhiGammaNeutral = source.fhDeltaPhiGammaNeutral ;
 179   fhDeltaEtaGammaCharged = source.fhDeltaEtaGammaCharged ;
 180   fhDeltaEtaGammaNeutral = source.fhDeltaEtaGammaNeutral ;
 181
 182   fhAnglePairLeading = source.fhAnglePairLeading ;
 183   fhInvMassPairLeading = source.fhInvMassPairLeading ;
 184   fhNBkg = source. fhNBkg ; fhNLeading = source. fhNLeading ;
 185   fhNJet = source.fhNJet ; fhJetRatio = source.fhJetRatio ; fhJetPt = source.fhJetPt ;
 186   fhBkgRatio  = source.fhBkgRatio ; fhBkgPt = source.fhBkgPt ;
 187   fhJetFragment = source.fhJetFragment ; fhBkgFragment = source.fhBkgFragment ;
 188   fhJetPtDist = source.fhJetPtDist ;  fhBkgPtDist = source.fhBkgPtDist ;
 189
 190
 191   for(Int_t i = 0; i<10; i++){
 192     fJetCones[i]        = source.fJetCones[i] ;
 193     fJetNameCones[i]    = source.fJetNameCones[i] ;
 194     fJetPtThres[i]      = source.fJetPtThres[i] ;
 195     fJetNamePtThres[i]  = source.fJetNamePtThres[i] ;
 196     if( i < 6 ){
 197       fJetXMin1[i]       = source.fJetXMin1[i] ;
 198       fJetXMin2[i]       = source.fJetXMin2[i] ;
 199       fJetXMax1[i]       = source.fJetXMax1[i] ;
 200       fJetXMax2[i]       = source.fJetXMax2[i] ;
 201       fBkgMean[i]        = source.fBkgMean[i] ;
 202       fBkgRMS[i]         = source.fBkgRMS[i] ;
 203       if( i < 2 ){
 204         fJetE1[i]        = source.fJetE1[i] ;
 205         fJetE2[i]        = source.fJetE2[i] ;
 206         fJetSigma1[i]    = source.fJetSigma1[i] ;
 207         fJetSigma2[i]    = source.fJetSigma2[i] ;
 208         fPhiEMCALCut[i]  = source.fPhiEMCALCut[i] ;
 209       }
 210     }
 211   }
 212
 213   return *this;
 214
 215 }
 216
 217 //____________________________________________________________________________
 218 AliAnaGammaJetLeadCone::~AliAnaGammaJetLeadCone()
 219 {
 220    // Remove all pointers except analysis output pointers.
 221
 222 }
 223
 224
 225
 226 //________________________________________________________________________
 227 TList *  AliAnaGammaJetLeadCone::GetCreateOutputObjects()
 228 {
 229   // Create histograms to be saved in output file and
 230   // store them in fOutputContainer
 231
 232   AliDebug(1,"Init jet in leading cone histograms");
 233
 234   TList * outputContainer = new TList() ;
 235   outputContainer->SetName("GammaJetCorrelationHistos") ;
 236
 237   fhChargedRatio  = new TH2F
 238     ("ChargedRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
 239      120,0,120,120,0,1);
 240   fhChargedRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
 241   fhChargedRatio->SetXTitle("p_{T #gamma} (GeV/c)");
 242
 243   fhPhiCharged  = new TH2F
 244     ("PhiCharged","#phi_{#pi^{#pm}}  vs p_{T #gamma}",
 245      120,0,120,120,0,7);
 246   fhPhiCharged->SetYTitle("#phi_{#pi^{#pm}} (rad)");
 247   fhPhiCharged->SetXTitle("p_{T #gamma} (GeV/c)");
 248
 249   fhEtaCharged  = new TH2F
 250     ("EtaCharged","#eta_{#pi^{#pm}}  vs p_{T #gamma}",
 251      120,0,120,120,-1,1);
 252   fhEtaCharged->SetYTitle("#eta_{#pi^{#pm}} (rad)");
 253   fhEtaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
 254
 255   fhDeltaPhiGammaCharged  = new TH2F
 256     ("DeltaPhiGammaCharged","#phi_{#gamma} - #phi_{charged #pi} vs p_{T #gamma}",
 257      200,0,120,200,0,6.4);
 258   fhDeltaPhiGammaCharged->SetYTitle("#Delta #phi");
 259   fhDeltaPhiGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
 260
 261   fhDeltaEtaGammaCharged  = new TH2F
 262     ("DeltaEtaGammaCharged","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
 263      200,0,120,200,-2,2);
 264   fhDeltaEtaGammaCharged->SetYTitle("#Delta #eta");
 265   fhDeltaEtaGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
 266
 267   outputContainer->Add(fhPhiCharged) ;
 268   outputContainer->Add(fhEtaCharged) ;
 269   outputContainer->Add(fhChargedRatio) ;
 270   outputContainer->Add(fhDeltaPhiGammaCharged) ;
 271   outputContainer->Add(fhDeltaEtaGammaCharged) ;
 272
 273   if(!AreJetsOnlyInCTS()){
 274
 275     fhNeutralRatio  = new TH2F
 276       ("NeutralRatio","p_{T leading  #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
 277        120,0,120,120,0,1);
 278     fhNeutralRatio->SetYTitle("p_{T lead  #pi^{0}} /p_{T #gamma}");
 279     fhNeutralRatio->SetXTitle("p_{T #gamma} (GeV/c)");
 280
 281     //
 282     fhAnglePairLeading  = new TH2F
 283       ("AnglePairLeading",
 284        "Angle between all #gamma pair finally selected vs p_{T  #pi^{0}}",
 285        200,0,50,200,0,0.2);
 286     fhAnglePairLeading->SetYTitle("Angle (rad)");
 287     fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
 288
 289     fhInvMassPairLeading  = new TH2F
 290       ("InvMassPairLeading",
 291        "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
 292        120,0,120,360,0,0.5);
 293     fhInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
 294     fhInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
 295
 296     fhPhiNeutral  = new TH2F
 297       ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T #gamma}",
 298        120,0,120,120,0,7);
 299     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
 300     fhPhiNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
 301
 302     fhEtaNeutral  = new TH2F
 303       ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T #gamma}",
 304        120,0,120,120,-1,1);
 305     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
 306     fhEtaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
 307
 308     fhDeltaPhiGammaNeutral  = new TH2F
 309       ("DeltaPhiGammaNeutral","#phi_{#gamma} - #phi_{#pi^{0}} vs p_{T #gamma}",
 310        200,0,120,200,0,6.4);
 311     fhDeltaPhiGammaNeutral->SetYTitle("#Delta #phi");
 312     fhDeltaPhiGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
 313
 314     fhDeltaEtaGammaNeutral  = new TH2F
 315       ("DeltaEtaGammaNeutral","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
 316        200,0,120,200,-2,2);
 317     fhDeltaEtaGammaNeutral->SetYTitle("#Delta #eta");
 318     fhDeltaEtaGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
 319
 320     outputContainer->Add(fhPhiNeutral) ;
 321     outputContainer->Add(fhEtaNeutral) ;
 322     outputContainer->Add(fhNeutralRatio) ;
 323     outputContainer->Add(fhDeltaPhiGammaNeutral) ;
 324     outputContainer->Add(fhDeltaEtaGammaNeutral) ;
 325
 326     outputContainer->Add(fhInvMassPairLeading) ;
 327     outputContainer->Add(fhAnglePairLeading) ;
 328   }
 329
 330   if(!fSeveralConeAndPtCuts){// not several cones
 331
 332     //Count
 333     fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000);
 334     fhNBkg->SetYTitle("counts");
 335     fhNBkg->SetXTitle("N");
 336     outputContainer->Add(fhNBkg) ;
 337
 338     fhNLeading  = new TH2F
 339       ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120);
 340     fhNLeading->SetYTitle("p_{T charge} (GeV/c)");
 341     fhNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
 342     outputContainer->Add(fhNLeading) ;
 343
 344     fhNJet  = new TH1F("NJet","Accepted jets",240,0,120);
 345     fhNJet->SetYTitle("N");
 346     fhNJet->SetXTitle("p_{T #gamma}(GeV/c)");
 347     outputContainer->Add(fhNJet) ;
 348
 349     //Ratios and Pt dist of reconstructed (not selected) jets
 350     //Jet
 351     fhJetRatio  = new TH2F
 352       ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
 353        240,0,120,200,0,10);
 354     fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
 355     fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
 356     outputContainer->Add(fhJetRatio) ;
 357
 358     fhJetPt  = new TH2F
 359       ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
 360     fhJetPt->SetYTitle("p_{T jet}");
 361     fhJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
 362     outputContainer->Add(fhJetPt) ;
 363
 364     //Bkg
 365
 366     fhBkgRatio  = new TH2F
 367       ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
 368        240,0,120,200,0,10);
 369     fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
 370     fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
 371     outputContainer->Add(fhBkgRatio) ;
 372
 373     fhBkgPt  = new TH2F
 374       ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
 375     fhBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
 376     fhBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
 377     outputContainer->Add(fhBkgPt) ;
 378
 379     //Jet Distributions
 380
 381     fhJetFragment  =
 382       new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
 383                240,0.,120.,1000,0.,1.2);
 384     fhJetFragment->SetYTitle("x_{T}");
 385     fhJetFragment->SetXTitle("p_{T #gamma}");
 386     outputContainer->Add(fhJetFragment) ;
 387
 388     fhBkgFragment  = new TH2F
 389       ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
 390        240,0.,120.,1000,0.,1.2);
 391     fhBkgFragment->SetYTitle("x_{T}");
 392     fhBkgFragment->SetXTitle("p_{T #gamma}");
 393     outputContainer->Add(fhBkgFragment) ;
 394
 395     fhJetPtDist  =
 396       new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
 397     fhJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
 398     outputContainer->Add(fhJetPtDist) ;
 399
 400     fhBkgPtDist  = new TH2F
 401       ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
 402     fhBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
 403     outputContainer->Add(fhBkgPtDist) ;
 404
 405   }//not several cones
 406   else{ //If we want to study the jet for different cones and pt
 407     for(Int_t icone = 0; icone<fJetNCone; icone++){//icone
 408       for(Int_t ipt = 0; ipt<fJetNPt;ipt++){ //ipt
 409
 410         //Jet
 411
 412         fhJetRatios[icone][ipt]  = new TH2F
 413           ("JetRatioCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 414            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
 415            +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c",
 416            240,0,120,200,0,10);
 417         fhJetRatios[icone][ipt]->
 418           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
 419         fhJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 420         outputContainer->Add(fhJetRatios[icone][ipt]) ;
 421
 422
 423         fhJetPts[icone][ipt]  = new TH2F
 424           ("JetPtCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 425            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
 426            +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c",
 427            240,0,120,400,0,200);
 428         fhJetPts[icone][ipt]->
 429           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
 430         fhJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 431         outputContainer->Add(fhJetPts[icone][ipt]) ;
 432
 433         //Bkg
 434         fhBkgRatios[icone][ipt]  = new TH2F
 435           ("BkgRatioCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 436            "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
 437            +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c",
 438            240,0,120,200,0,10);
 439         fhBkgRatios[icone][ipt]->
 440           SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
 441         fhBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 442         outputContainer->Add(fhBkgRatios[icone][ipt]) ;
 443
 444         fhBkgPts[icone][ipt]  = new TH2F
 445           ("BkgPtCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 446            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
 447            +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c",
 448            240,0,120,400,0,200);
 449         fhBkgPts[icone][ipt]->
 450           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
 451         fhBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 452         outputContainer->Add(fhBkgPts[icone][ipt]) ;
 453
 454         //Counts
 455         fhNBkgs[icone][ipt]  = new TH1F
 456           ("NBkgCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 457            "bkg multiplicity cone ="+fJetNameCones[icone]+", pt>"
 458            +fJetNamePtThres[ipt]+" GeV/c",9000,0,9000);
 459         fhNBkgs[icone][ipt]->SetYTitle("counts");
 460         fhNBkgs[icone][ipt]->SetXTitle("N");
 461         outputContainer->Add(fhNBkgs[icone][ipt]) ;
 462
 463         fhNLeadings[icone][ipt]  = new TH2F
 464           ("NLeadingCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 465            "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fJetNameCones[icone]+", pt>"
 466            +fJetNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
 467         fhNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
 468         fhNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
 469         outputContainer->Add(fhNLeadings[icone][ipt]) ;
 470
 471         fhNJets[icone][ipt]  = new TH1F
 472           ("NJetCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 473            "Number of neutral jets, cone ="+fJetNameCones[icone]+", pt>"
 474            +fJetNamePtThres[ipt]+" GeV/c",120,0,120);
 475         fhNJets[icone][ipt]->SetYTitle("N");
 476         fhNJets[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
 477         outputContainer->Add(fhNJets[icone][ipt]) ;
 478
 479         //Fragmentation Function
 480         fhJetFragments[icone][ipt]  = new TH2F
 481           ("JetFragmentCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 482            "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fJetNameCones[icone]+", pt>"
 483            +fJetNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
 484         fhJetFragments[icone][ipt]->SetYTitle("x_{T}");
 485         fhJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
 486         outputContainer->Add(fhJetFragments[icone][ipt]) ;
 487
 488         fhBkgFragments[icone][ipt]  = new TH2F
 489           ("BkgFragmentCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 490            "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fJetNameCones[icone]+", pt>"
 491            +fJetNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
 492         fhBkgFragments[icone][ipt]->SetYTitle("x_{T}");
 493         fhBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
 494         outputContainer->Add(fhBkgFragments[icone][ipt]) ;
 495
 496         //Jet particle distribution
 497
 498         fhJetPtDists[icone][ipt]  = new TH2F
 499           ("JetPtDistCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 500            "p_{T i}, cone ="+fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+
 501            " GeV/c",120,0.,120.,120,0.,120.);
 502         fhJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 503         outputContainer->Add(fhJetPtDists[icone][ipt]) ;
 504
 505         fhBkgPtDists[icone][ipt]  = new TH2F
 506           ("BkgPtDistCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
 507            "p_{T i}, cone ="+fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+
 508            " GeV/c",120,0.,120.,120,0.,120.);
 509         fhBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 510         outputContainer->Add(fhBkgPtDists[icone][ipt]) ;
 511
 512       }//ipt
 513     } //icone
 514   }//If we want to study any cone or pt threshold
 515
 516   SetOutputContainer(outputContainer);
 517
 518   return outputContainer;
 519 }
 520
 521   //____________________________________________________________________________
 522   void AliAnaGammaJetLeadCone::InitParameters()
 523 {
 524
 525   //Initialize the parameters of the analysis.
 526   SetJetsOnlyInCTS(kFALSE) ;
 527   fPbPb                = kFALSE ;
 528   fEtaEMCALCut     = 0.7 ;
 529   fPhiEMCALCut[0] = 80. *TMath::Pi()/180.;
 530   fPhiEMCALCut[1] = 190.*TMath::Pi()/180.;
 531   SetDeltaPhiCutRange(2.9,3.4) ;
 532   SetRatioCutRange(0.1,1.0) ;
 533
 534   //Jet selection parameters
 535   //Fixed cut
 536   fJetRatioMaxCut = 1.2 ;
 537   fJetRatioMinCut = 0.3 ;
 538   fJetCTSRatioMaxCut = 1.2 ;
 539   fJetCTSRatioMinCut = 0.3 ;
 540   fSelect         = 0  ;
 541
 542   //Cut depending on gamma energy
 543
 544   fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applied
 545   //Reconstructed jet energy dependence parameters
 546   //e_jet = a1+e_gamma b2.
 547   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
 548   fJetE1[0] = -5.75; fJetE1[1] = -4.1;
 549   fJetE2[0] = 1.005; fJetE2[1] = 1.05;
 550
 551   //Reconstructed sigma of jet energy dependence parameters
 552   //s_jet = a1+e_gamma b2.
 553   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
 554   fJetSigma1[0] = 2.65;   fJetSigma1[1] = 2.75;
 555   fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
 556
 557   //Background mean energy and RMS
 558   //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
 559   //Index 2-> (low pt jets)BKG > 0.5 GeV;
 560   //Index > 2, same for CTS conf
 561   fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
 562   fBkgMean[3] = 0.; fBkgMean[4] = 6.4;  fBkgMean[5] = 48.6;
 563   fBkgRMS[0]  = 0.; fBkgRMS[1]  = 7.5;  fBkgRMS[2]  = 22.0;
 564   fBkgRMS[3]  = 0.; fBkgRMS[4]  = 5.4;  fBkgRMS[5]  = 13.2;
 565
 566   //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
 567   //limits for monoenergetic jets.
 568   //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
 569   //Index 2-> (low pt jets) BKG > 0.5 GeV;
 570   //Index > 2, same for CTS conf
 571
 572   fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
 573   fJetXMin1[3] =-2.0  ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1  ;
 574   fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
 575   fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
 576   fJetXMax1[0] =-3.8  ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6  ;
 577   fJetXMax1[3] =-2.7  ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7  ;
 578   fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
 579   fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
 580
 581
 582   //Different cones and pt thresholds to construct the jet
 583
 584   fJetCone        = 0.3  ;
 585   fJetPtThreshold = 0.5   ;
 586   fJetPtThresPbPb = 2.   ;
 587   fJetNCone       = 4    ;
 588   fJetNPt         = 4    ;
 589   fJetCones[0]    = 0.2  ; fJetNameCones[0]   = "02" ;
 590   fJetCones[1]    = 0.3  ; fJetNameCones[1]   = "03" ;
 591   fJetCones[2]    = 0.4  ; fJetNameCones[2]   = "04" ;
 592   fJetCones[2]    = 0.5  ; fJetNameCones[2]   = "05" ;
 593
 594   fJetPtThres[0]  = 0.0  ; fJetNamePtThres[0] = "00" ;
 595   fJetPtThres[1]  = 0.5  ; fJetNamePtThres[1] = "05" ;
 596   fJetPtThres[2]  = 1.0  ; fJetNamePtThres[2] = "10" ;
 597   fJetPtThres[3]  = 2.0  ; fJetNamePtThres[3] = "20" ;
 598 }
 599
 600 //__________________________________________________________________
 601 void AliAnaGammaJetLeadCone::Print(const Option_t * opt) const
 602 {
 603
 604   //Print some relevant parameters set for the analysis
 605   if(! opt)
 606     return;
 607
 608   Info("Print", "%s %s", GetName(), GetTitle() ) ;
 609   printf("Correlation analysis           =     %d\n",kJetLeadCone) ;
 610
 611
 612   printf("Phi gamma-Leading      <     %f\n", GetDeltaPhiMaxCut()) ;
 613   printf("Phi gamma-Leading      >     %f\n", GetDeltaPhiMinCut()) ;
 614   printf("pT Leading / pT Gamma             <     %f\n", GetRatioMaxCut()) ;
 615   printf("pT Leading / pT Gamma             >     %f\n", GetRatioMinCut()) ;
 616
 617   if(fSelect == 2){
 618     printf("pT Jet / pT Gamma                     <    %f\n", fJetRatioMaxCut) ;
 619     printf("pT Jet / pT Gamma                     >    %f\n", fJetRatioMinCut) ;
 620     printf("pT Jet (Only CTS)/ pT Gamma   <    %f\n", fJetCTSRatioMaxCut) ;
 621     printf("pT Jet (Only CTS)/ pT Gamma   >    %f\n", fJetCTSRatioMinCut) ;
 622   }
 623
 624 }
 625
 626 //__________________________________________________________________
 627 void  AliAnaGammaJetLeadCone::MakeGammaCorrelation(TParticle * pGamma, TClonesArray * plCTS,  TClonesArray * plNe)
 628 {
 629   //Gamma Jet Correlation Analysis
 630   AliDebug(2, "Make correlation");
 631
 632   TParticle *pLeading = new TParticle; //It will contain the kinematics of the found leading particle
 633
 634   //Search leading particles in CTS and EMCAL
 635   if(GetLeadingParticle(plCTS, plNe, pGamma, pLeading)){
 636
 637     AliDebug(1,Form("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ;
 638
 639     //Search Jet
 640     if(!fSeveralConeAndPtCuts)
 641       MakeJet(plCTS, plNe, pGamma, pLeading, "");
 642     else{
 643       for(Int_t icone = 0; icone<fJetNCone; icone++) {
 644         for(Int_t ipt = 0; ipt<fJetNPt;ipt++) {
 645           TString lastname ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt];
 646           fJetCone=fJetCones[icone];
 647           fJetPtThreshold=fJetPtThres[ipt];
 648           MakeJet(plCTS, plNe, pGamma, pLeading, lastname);
 649         }//icone
 650       }//ipt
 651     }//fSeveralConeAndPtCuts
 652   }//Leading
 653
 654   AliDebug(2, "End of analysis, delete pointers");
 655
 656   pLeading->Delete();
 657
 658 }
 659
 660 //____________________________________________________________________________
 661 Bool_t  AliAnaGammaJetLeadCone::GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe,
 662                                             TParticle * pGamma, TParticle * pLeading)
 663 {
 664   //Search Charged or Neutral leading particle, select the highest one.
 665
 666   TParticle * pLeadingCh = new TParticle();
 667   TParticle * pLeadingPi0 = new TParticle();
 668
 669   Double_t ptg  =  pGamma->Pt();
 670   Double_t phig = pGamma->Phi();
 671   Double_t etag = pGamma->Eta();
 672
 673   if(!AreJetsOnlyInCTS())
 674     {
 675       AliDebug(3, "GetLeadingPi0");
 676       GetLeadingPi0(plNe, pGamma, pLeadingPi0) ;
 677       AliDebug(3, "GetLeadingCharge");
 678       GetLeadingCharge(plCTS, pGamma, pLeadingCh) ;
 679
 680       Double_t ptch = pLeadingCh->Pt();
 681       Double_t phich = pLeadingCh->Phi();
 682       Double_t etach = pLeadingCh->Eta();
 683       Double_t ptpi = pLeadingPi0->Pt();
 684       Double_t phipi = pLeadingPi0->Phi();
 685       Double_t etapi = pLeadingPi0->Eta();
 686
 687       //Is leading cone inside EMCAL acceptance?
 688
 689       Bool_t insidech = kFALSE ;
 690       if((phich - fJetCone) >  fPhiEMCALCut[0] &&
 691          (phich + fJetCone) <  fPhiEMCALCut[1] &&
 692         (etach-fJetCone) < fEtaEMCALCut )
 693         insidech = kTRUE ;
 694
 695       Bool_t insidepi = kFALSE ;
 696       if((phipi - fJetCone) >  fPhiEMCALCut[0] &&
 697          (phipi + fJetCone) <  fPhiEMCALCut[1] &&
 698         (etapi-fJetCone) < fEtaEMCALCut )
 699         insidepi = kTRUE ;
 700
 701       AliDebug(2,Form("Leading:  charged pt %f, pi0 pt  %f",ptch,ptpi)) ;
 702
 703       if (ptch > 0 || ptpi > 0)
 704         {
 705           if(insidech && (ptch > ptpi))
 706             {
 707
 708               AliDebug(1,"Leading found in CTS");
 709               pLeading->SetMomentum(pLeadingCh->Px(),pLeadingCh->Py(),pLeadingCh->Pz(),pLeadingCh->Energy());
 710               AliDebug(3,Form("Final leading found in CTS, pt %f, phi %f, eta %f",ptch,phich,etach)) ;
 711               fhChargedRatio->Fill(ptg,ptch/pGamma->Pt());
 712               fhEtaCharged->Fill(ptg,etach);
 713               fhPhiCharged->Fill(ptg,phich);
 714               fhDeltaPhiGammaCharged->Fill(ptg,pGamma->Phi()-phich);
 715               fhDeltaEtaGammaCharged->Fill(ptg,pGamma->Eta()-etach);
 716               return 1;
 717             }
 718
 719           else if((ptpi > ptch) && insidepi)
 720             {
 721               AliDebug(1,"Leading found in EMCAL");
 722               pLeading->SetMomentum(pLeadingPi0->Px(),pLeadingPi0->Py(),pLeadingPi0->Pz(),pLeadingPi0->Energy());
 723               AliDebug(3,Form("Final leading found in EMCAL, pt %f, phi %f, eta %f",ptpi,phipi,etapi)) ;
 724               fhNeutralRatio ->Fill(ptg,ptpi/ptg);
 725               fhEtaNeutral->Fill(ptg,etapi);
 726               fhPhiNeutral->Fill(ptg,phipi);
 727               fhDeltaPhiGammaNeutral->Fill(ptg,phig-phipi);
 728               fhDeltaEtaGammaNeutral->Fill(ptg,etag-etapi);
 729               return 1;
 730             }
 731
 732           else{
 733             AliDebug(3,"NO LEADING PARTICLE FOUND");}
 734           return 0;
 735         }
 736       else{
 737         AliDebug(3,"NO LEADING PARTICLE FOUND");
 738         return 0;
 739       }
 740     }
 741   else
 742     {
 743       //No calorimeter present for Leading particle detection
 744       GetLeadingCharge(plCTS, pGamma, pLeading) ;
 745       Double_t ptch = pLeading->Pt();
 746       Double_t phich = pLeading->Phi();
 747       Double_t etach = pLeading->Eta();
 748       if(ptch > 0){
 749         fhChargedRatio->Fill(ptg,ptch/ptg);
 750         fhEtaCharged->Fill(ptg,etach);
 751         fhPhiCharged->Fill(ptg,phich);
 752         fhDeltaPhiGammaCharged->Fill(ptg,phig-phich);
 753         fhDeltaEtaGammaCharged->Fill(ptg,etag-etach);
 754         AliDebug(3,Form("Leading found :  pt %f, phi %f, eta %f",ptch,phich,etach)) ;
 755         return 1;
 756       }
 757       else
 758         {
 759           AliDebug(3,"NO LEADING PARTICLE FOUND");
 760           return 0;
 761         }
 762     }
 763 }
 764
 765 //____________________________________________________________________________
 766 void  AliAnaGammaJetLeadCone::GetLeadingCharge(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) const
 767 {
 768   //Search for the charged particle with highest with
 769   //Phi=Phi_gamma-Pi and pT=0.1E_gamma
 770   Double_t pt  = -100.;
 771   Double_t phi = -100.;
 772
 773   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
 774
 775     TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
 776
 777     Double_t ptl  = particle->Pt();
 778     Double_t rat  = ptl/pGamma->Pt() ;
 779     Double_t phil = particle->Phi() ;
 780     Double_t phig = pGamma->Phi();
 781
 782     //Selection within angular and energy limits
 783     if(((phig-phil)> GetDeltaPhiMinCut()) && ((phig-phil)<GetDeltaPhiMaxCut()) &&
 784         (rat > GetRatioMinCut()) && (rat < GetRatioMaxCut())  && (ptl  > pt)) {
 785        phi = phil ;
 786        pt  = ptl ;
 787        pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
 788        AliDebug(4,Form("Charge in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, particle->Eta(), phi,rat)) ;
 789      }
 790   }
 791
 792   AliDebug(1,Form("Leading in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, pLeading->Eta(), phi,pt/pGamma->Pt())) ;
 793
 794 }
 795
 796
 797 //____________________________________________________________________________
 798 void  AliAnaGammaJetLeadCone::GetLeadingPi0(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading)
 799 {
 800
 801   //Search for the neutral pion with highest with
 802   //Phi=Phi_gamma-Pi and pT=0.1E_gamma
 803   Double_t pt  = -100.;
 804   Double_t phi = -100.;
 805   Double_t ptl = -100.;
 806   Double_t rat = -100.;
 807   Double_t phil = -100. ;
 808   Double_t ptg  = pGamma->Pt();
 809   Double_t phig = pGamma->Phi();
 810
 811   TIter next(pl);
 812   TParticle * particlei = 0;
 813   TParticle * particlej = 0;
 814   TLorentzVector gammai;
 815   TLorentzVector gammaj;
 816
 817   Int_t iPrimary = -1;
 818   Double_t angle = 0., e = 0., invmass = 0.;
 819   Int_t ksPdg = 0;
 820   Int_t jPrimary=-1;
 821
 822   while ( (particlei = (TParticle*)next()) ) {
 823     iPrimary++;
 824     ksPdg = particlei->GetPdgCode();
 825
 826     ptl  = particlei->Pt();
 827     //2 gamma overlapped, found with PID
 828      if(ksPdg == 111){
 829        rat = ptl/ptg ;
 830        phil = particlei->Phi() ;
 831        //Selection within angular and energy limits
 832        if(ptl > pt  && rat > GetRatioMinCut()  && rat < GetRatioMaxCut()  &&
 833           (phig-phil) > GetDeltaPhiMinCut() && (phig-phil) < GetDeltaPhiMaxCut() )
 834          {
 835            phi = phil ;
 836            pt  = ptl ;
 837            pLeading->SetMomentum(particlei->Px(),particlei->Py(),particlei->Pz(),particlei->Energy());
 838            AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
 839          }// cuts
 840      }// pdg = 111
 841
 842     //Make invariant mass analysis
 843     else if(ksPdg == 22){//  gamma i
 844
 845       //Search the photon companion in case it comes from  a Pi0 decay
 846       //Apply several cuts to select the good pair
 847       particlei->Momentum(gammai);
 848       jPrimary=-1;
 849       TIter next2(pl);
 850       while ( (particlej = (TParticle*)next2()) ) {
 851         jPrimary++;
 852         if(jPrimary>iPrimary){
 853           ksPdg = particlej->GetPdgCode();
 854           particlej->Momentum(gammaj);
 855           if(ksPdg == 22 && GetNeutralMesonSelection()->SelectPair(pGamma, gammai, gammaj))
 856           {
 857             //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
 858             //gammai.Pt(),gammaj.Pt());
 859             ptl  = (gammai+gammaj).Pt();
 860             phil = (gammai+gammaj).Phi();
 861             if(phil < 0)
 862               phil+=TMath::TwoPi();
 863             rat = ptl/ptg ;
 864
 865             if(ptl > pt ){
 866               pt       = ptl;
 867               phi      = phil ;
 868               e       = (gammai+gammaj).E();
 869               angle   = gammaj.Angle(gammai.Vect());
 870               invmass = (gammai+gammaj).M();
 871               pLeading->SetMomentum( (gammai+gammaj).Px(), (gammai+gammaj).Py(), (gammai+gammaj).Pz(), (gammai+gammaj).E());
 872               AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
 873             }
 874           }//if pdg = 22, pair selected
 875         }//iprimary<jprimary
 876       }//while
 877     }// if pdg = 22
 878   }//while
 879
 880   //Final pi0 found, highest pair energy, fill histograms
 881   if(e > 0.){
 882     fhInvMassPairLeading->Fill(e,invmass);
 883     fhAnglePairLeading->Fill(e,angle);
 884   }
 885
 886   AliDebug(3,Form("Leading EMCAL: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
 887 }
 888
 889
 890 //__________________________________________________________________________-
 891 Bool_t AliAnaGammaJetLeadCone::IsJetSelected(const Double_t ptg, const Double_t ptj){
 892   //Check if the energy of the reconstructed jet is within an energy window
 893
 894   Double_t par[6];
 895   Double_t xmax[2];
 896   Double_t xmin[2];
 897
 898   Int_t iCTS = 0;
 899   if(AreJetsOnlyInCTS())
 900     iCTS = 3 ;
 901
 902   if(!fPbPb){
 903     //Phythia alone, jets with pt_th > 0.2, r = 0.3
 904     par[0] = fJetE1[0]; par[1] = fJetE2[0];
 905     //Energy of the jet peak
 906     //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
 907     par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
 908     //Sigma  of the jet peak
 909     //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
 910     par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
 911     //Parameters reserved for PbPb bkg.
 912     xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
 913     xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
 914     //Factor that multiplies sigma to obtain the best limits,
 915     //by observation, of mono jet ratios (ptjet/ptg)
 916     //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
 917
 918   }
 919   else{
 920     if(ptg > fPtJetSelectionCut){
 921       //Phythia +PbPb with  pt_th > 2 GeV/c, r = 0.3
 922       par[0] = fJetE1[0]; par[1] = fJetE2[0];
 923       //Energy of the jet peak, same as in pp
 924       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
 925       par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
 926       //Sigma  of the jet peak, same as in pp
 927       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
 928       par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
 929       //Mean value and RMS of PbPb Bkg
 930       xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
 931       xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
 932       //Factor that multiplies sigma to obtain the best limits,
 933       //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg,
 934       //pt_th > 2 GeV, r = 0.3
 935       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
 936
 937     }
 938     else{
 939       //Phythia + PbPb with  pt_th > 0.5 GeV/c, r = 0.3
 940       par[0] = fJetE1[1]; par[1] = fJetE2[1];
 941       //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
 942       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
 943       par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
 944       //Sigma  of the jet peak, pt_th > 2 GeV/c, r = 0.3
 945       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
 946       par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
 947       //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
 948       xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
 949       xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
 950       //Factor that multiplies sigma to obtain the best limits,
 951       //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg,
 952       //pt_th > 2 GeV, r = 0.3
 953       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
 954
 955     }//If low pt jet in bkg
 956   }//if Bkg
 957
 958  //Calculate minimum and maximum limits of the jet ratio.
 959   Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
 960   Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
 961
 962   AliDebug(3,Form("Jet selection?  : Limits min %f, max %f,  pt_jet %f,  pt_gamma %f, pt_jet / pt_gamma %f",min,max,ptj,ptg,ptj/ptg));
 963
 964   if(( min < ptj/ptg ) && ( max > ptj/ptg))
 965     return kTRUE;
 966   else
 967     return kFALSE;
 968
 969 }
 970
 971 //____________________________________________________________________________
 972 void AliAnaGammaJetLeadCone::MakeJet(TClonesArray * plCTS, TClonesArray * plNe,
 973                                      TParticle * pGamma, TParticle* pLeading,TString lastname)
 974 {
 975   //Fill the jet with the particles around the leading particle with
 976   //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and
 977   //check if we select it. Fill jet histograms
 978
 979   TClonesArray * jetList = new TClonesArray("TParticle",1000);
 980   TClonesArray * bkgList = new TClonesArray("TParticle",1000);
 981
 982   TLorentzVector jet   (0,0,0,0);
 983   TLorentzVector bkg(0,0,0,0);
 984   TLorentzVector lv (0,0,0,0);
 985
 986   Double_t ptjet = 0.0;
 987   Double_t ptbkg = 0.0;
 988   Int_t n0 = 0;
 989   Int_t n1 = 0;
 990   Bool_t b1 = kFALSE;
 991   Bool_t b0 = kFALSE;
 992
 993   Double_t ptg  = pGamma->Pt();
 994   Double_t phig = pGamma->Phi();
 995   Double_t ptl  = pLeading->Pt();
 996   Double_t phil = pLeading->Phi();
 997   Double_t etal = pLeading->Eta();
 998
 999   Float_t ptcut = fJetPtThreshold;
1000   if(fPbPb && !fSeveralConeAndPtCuts && ptg > fPtJetSelectionCut)  ptcut = fJetPtThresPbPb ;
1001
1002   //Add charged particles to jet
1003   TIter nextch(plCTS) ;
1004   TParticle * particle = 0 ;
1005   while ( (particle = dynamic_cast<TParticle*>(nextch())) ) {
1006
1007     b0 = kFALSE;
1008     b1 = kFALSE;
1009
1010     //Particles in jet
1011     SetJet(particle, b0, fJetCone, etal, phil) ;
1012
1013     if(b0){
1014       new((*jetList)[n0++]) TParticle(*particle) ;
1015       particle->Momentum(lv);
1016       if(particle->Pt() > ptcut ){
1017         jet+=lv;
1018         ptjet+=particle->Pt();
1019       }
1020     }
1021
1022     //Background around (phi_gamma-pi, eta_leading)
1023     SetJet(particle, b1, fJetCone,etal, phig) ;
1024
1025     if(b1) {
1026       new((*bkgList)[n1++]) TParticle(*particle) ;
1027       particle->Momentum(lv);
1028       if(particle->Pt() > ptcut ){
1029         bkg+=lv;
1030         ptbkg+=particle->Pt();
1031       }
1032     }
1033   }
1034
1035    //Add neutral particles to jet
1036   TIter nextne(plNe) ;
1037   particle = 0 ;
1038   while ( (particle = dynamic_cast<TParticle*>(nextne())) ) {
1039
1040     b0 = kFALSE;
1041     b1 = kFALSE;
1042
1043     //Particles in jet
1044     SetJet(particle, b0, fJetCone, etal, phil) ;
1045
1046     if(b0){
1047       new((*jetList)[n0++]) TParticle(*particle) ;
1048       particle->Momentum(lv);
1049       if(particle->Pt() > ptcut ){
1050         jet+=lv;
1051         ptjet+=particle->Pt();
1052       }
1053     }
1054
1055     //Background around (phi_gamma-pi, eta_leading)
1056     SetJet(particle, b1, fJetCone,etal, phig) ;
1057
1058     if(b1) {
1059       new((*bkgList)[n1++]) TParticle(*particle) ;
1060       particle->Momentum(lv);
1061       if(particle->Pt() > ptcut ){
1062         bkg+=lv;
1063         ptbkg+=particle->Pt();
1064       }
1065     }
1066   }
1067
1068   ptjet = jet.Pt();
1069   ptbkg = bkg.Pt();
1070
1071   if(ptjet > 0.) {
1072
1073     AliDebug(2,Form("Gamma   pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg));
1074
1075     //Fill histograms
1076
1077     Double_t ratjet   = ptjet/ptg ;
1078     Double_t ratbkg  = ptbkg/ptg ;
1079
1080     dynamic_cast<TH2F*>
1081       (GetOutputContainer()->FindObject("JetRatio"+lastname))
1082       ->Fill(ptg,ratjet);
1083     dynamic_cast<TH2F*>
1084       (GetOutputContainer()->FindObject("JetPt"+lastname))
1085       ->Fill(ptg,ptjet);
1086
1087     dynamic_cast<TH2F*>
1088       (GetOutputContainer()->FindObject("BkgRatio"+lastname))
1089       ->Fill(ptg,ratbkg);
1090
1091     dynamic_cast<TH2F*>
1092       (GetOutputContainer()->FindObject("BkgPt"+lastname))
1093       ->Fill(ptg,ptbkg);
1094
1095     //Jet selection
1096     Bool_t kSelect = kFALSE;
1097     if(fSelect == 0)
1098       kSelect = kTRUE; //Accept all jets, no restriction
1099     else if(fSelect == 1){
1100       //Selection with parametrized cuts
1101       if(IsJetSelected(ptg,ptjet))   kSelect = kTRUE;
1102     }
1103     else if(fSelect == 2){
1104       //Simple selection
1105       if(!AreJetsOnlyInCTS()){
1106         if((ratjet <  fJetRatioMaxCut) && (ratjet > fJetRatioMinCut )) kSelect = kTRUE;
1107       }
1108       else{
1109         if((ratjet <  fJetCTSRatioMaxCut) && (ratjet > fJetCTSRatioMinCut )) kSelect = kTRUE;
1110       }
1111     }
1112     else
1113       AliError("Jet selection option larger than 2, DON'T SELECT JETS");
1114
1115
1116      if(kSelect){
1117        AliDebug(1,Form("Jet Selected: pt %f ", ptjet)) ;
1118
1119        FillJetHistos(jetList, ptg, ptl,"Jet",lastname);
1120        FillJetHistos(bkgList, ptg, ptl, "Bkg",lastname);
1121      }
1122   } //ptjet > 0
1123
1124   jetList ->Delete();
1125   bkgList ->Delete();
1126
1127 }
1128
1129 //___________________________________________________________________
1130 void AliAnaGammaJetLeadCone::SetJet(TParticle * part, Bool_t & b, Float_t cone,
1131                              Double_t eta, Double_t phi)
1132 {
1133
1134   //Check if the particle is inside the cone defined by the leading particle
1135   b = kFALSE;
1136
1137   if(phi > TMath::TwoPi())
1138     phi-=TMath::TwoPi();
1139   if(phi < 0.)
1140     phi+=TMath::TwoPi();
1141
1142   Double_t  rad = 10000 + cone;
1143
1144   if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
1145     rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1146                       TMath::Power(part->Phi()-phi,2));
1147   else{
1148     if(part->Phi()-phi > TMath::TwoPi() - cone)
1149       rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1150                         TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
1151     if(part->Phi()-phi < -(TMath::TwoPi() - cone))
1152       rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1153                         TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
1154   }
1155
1156   if(rad < cone )
1157     b = kTRUE;
1158
1159 }
1160
1161 //____________________________________________________________________________
1162 void AliAnaGammaJetLeadCone::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname)
1163 {
1164   //Fill histograms wth jet fragmentation
1165   //and number of selected jets and leading particles
1166   //and the background multiplicity
1167   TParticle * particle = 0 ;
1168   Int_t ipr = 0;
1169   Float_t  charge = 0;
1170
1171   TIter next(pl) ;
1172   while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1173     ipr++ ;
1174     Double_t pt = particle->Pt();
1175
1176     charge = TDatabasePDG::Instance()
1177       ->GetParticle(particle->GetPdgCode())->Charge();
1178     if(charge != 0){//Only jet Charged particles
1179       dynamic_cast<TH2F*>
1180         (GetOutputContainer()->FindObject(type+"Fragment"+lastname))
1181         ->Fill(ptg,pt/ptg);
1182       dynamic_cast<TH2F*>
1183         (GetOutputContainer()->FindObject(type+"PtDist"+lastname))
1184         ->Fill(ptg,pt);
1185     }//charged
1186
1187   }//while
1188
1189   if(type == "Bkg")
1190     dynamic_cast<TH1F*>
1191       (GetOutputContainer()->FindObject("NBkg"+lastname))
1192       ->Fill(ipr);
1193   else{
1194     dynamic_cast<TH1F*>
1195       (GetOutputContainer()->FindObject("NJet"+lastname))->
1196       Fill(ptg);
1197     dynamic_cast<TH2F*>
1198       (GetOutputContainer()->FindObject("NLeading"+lastname))
1199       ->Fill(ptg,ptl);
1200   }
1201
1202 }