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