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