PWG4/AliAnaGammaJet.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id$ */
  16
  17 /* History of cvs commits:
  18  *
  19  * $Log$
  20  * Revision 1.3  2007/03/08 10:24:32  schutz
  21  * Coding convention
  22  *
  23  * Revision 1.2  2007/02/09 18:40:40  schutz
  24  * B\bNew version from Gustavo
  25  *
  26  * Revision 1.1  2007/01/23 17:17:29  schutz
  27  * New Gamma package
  28  *
  29  *
  30  */
  31
  32 //_________________________________________________________________________
  33 // Class for the analysis of gamma-jet correlations
  34 //  Basically it seaches for a prompt photon in the Calorimeters acceptance,
  35 //  if so we construct a jet around the highest pt particle in the opposite
  36 //  side in azimuth, inside the Central Tracking System (ITS+TPC) and
  37 //  EMCAL acceptances (only when PHOS detects the gamma). First the leading
  38 //  particle and then the jet have to fullfill several conditions
  39 //  (energy, direction ..) to be accepted. Then the fragmentation function
  40 //  of this jet is constructed
  41 //  Class created from old AliPHOSGammaJet
  42 //  (see AliRoot versions previous Release 4-09)
  43 //
  44 //*-- Author: Gustavo Conesa (LNF-INFN)
  45 //////////////////////////////////////////////////////////////////////////////
  46
  47
  48 // --- ROOT system ---
  49
  50 #include <TFile.h>
  51 #include <TParticle.h>
  52 #include <TH2.h>
  53
  54 #include "AliAnaGammaJet.h"
  55 #include "AliESD.h"
  56 #include "AliESDtrack.h"
  57 #include "AliESDCaloCluster.h"
  58 #include "Riostream.h"
  59 #include "AliLog.h"
  60
  61 ClassImp(AliAnaGammaJet)
  62
  63 //____________________________________________________________________________
  64 AliAnaGammaJet::AliAnaGammaJet(const char *name) :
  65   AliAnaGammaDirect(name),
  66   fSeveralConeAndPtCuts(0),
  67   fPbPb(kFALSE),
  68   fJetsOnlyInCTS(0),
  69   fEtaEMCALCut(0.),fPhiMaxCut(0.),
  70   fPhiMinCut(0.),
  71   fInvMassMaxCut(0.), fInvMassMinCut(0.), fRatioMaxCut(0), fRatioMinCut(0),
  72   fJetCTSRatioMaxCut(0.), fJetCTSRatioMinCut(0.), fJetRatioMaxCut(0.),
  73   fJetRatioMinCut(0.), fNCone(0),
  74   fNPt(0), fCone(0), fPtThreshold(0),
  75   fPtJetSelectionCut(0.0),
  76   fOutputContainer(new TObjArray(100)),  fAngleMaxParam(),  fSelect(0),
  77   fhChargeRatio(0), fhPi0Ratio (0),
  78   fhDeltaPhiCharge(0),  fhDeltaPhiPi0 (0), fhDeltaEtaCharge(0), fhDeltaEtaPi0(0),
  79   fhAnglePair(0), fhAnglePairAccepted(0), fhAnglePairNoCut(0), fhAnglePairLeadingCut(0),
  80   fhAnglePairAngleCut (0), fhAnglePairAllCut (0), fhAnglePairLeading(0),
  81   fhInvMassPairNoCut  (0), fhInvMassPairLeadingCut(0), fhInvMassPairAngleCut(0),
  82   fhInvMassPairAllCut (0), fhInvMassPairLeading(0),
  83   fhNBkg (0), fhNLeading(0), fhNJet(0), fhJetRatio(0), fhJetPt (0),
  84   fhBkgRatio (0), fhBkgPt(0),  fhJetFragment(0), fhBkgFragment(0),
  85   fhJetPtDist(0),  fhBkgPtDist(0)
  86
  87 {
  88
  89   // ctor
  90   fAngleMaxParam.Set(4) ;
  91   fAngleMaxParam.Reset(0.);
  92
  93   for(Int_t i = 0; i<10; i++){
  94     fCones[i]         = 0.0 ;
  95     fNameCones[i]     = ""  ;
  96     fPtThres[i]      = 0.0 ;
  97     fNamePtThres[i]  = ""  ;
  98     if( i < 6 ){
  99       fJetXMin1[i]     = 0.0 ;
 100       fJetXMin2[i]     = 0.0 ;
 101       fJetXMax1[i]     = 0.0 ;
 102       fJetXMax2[i]     = 0.0 ;
 103       fBkgMean[i]      = 0.0 ;
 104       fBkgRMS[i]       = 0.0 ;
 105       if( i < 2 ){
 106         fJetE1[i]        = 0.0 ;
 107         fJetE2[i]        = 0.0 ;
 108         fJetSigma1[i]    = 0.0 ;
 109         fJetSigma2[i]    = 0.0 ;
 110         fPhiEMCALCut[i]  = 0.0 ;
 111       }
 112     }
 113   }
 114
 115   //Init parameters
 116   InitParameters();
 117
 118   // Input slot #0 works with an Ntuple
 119   DefineInput(0, TChain::Class());
 120   // Output slot #0 writes into a TH1 container
 121   DefineOutput(0,  TObjArray::Class()) ;
 122
 123 }
 124
 125
 126 //____________________________________________________________________________
 127 AliAnaGammaJet::AliAnaGammaJet(const AliAnaGammaJet & gj) :
 128   AliAnaGammaDirect(gj),
 129   fSeveralConeAndPtCuts(gj.fSeveralConeAndPtCuts),
 130   fPbPb(gj.fPbPb), fJetsOnlyInCTS(gj.fJetsOnlyInCTS),
 131   fEtaEMCALCut(gj.fEtaEMCALCut),
 132   fPhiMaxCut(gj.fPhiMaxCut), fPhiMinCut(gj.fPhiMinCut),
 133   fInvMassMaxCut(gj.fInvMassMaxCut), fInvMassMinCut(gj.fInvMassMinCut),
 134   fRatioMaxCut(gj.fRatioMaxCut), fRatioMinCut(gj.fRatioMinCut),
 135   fJetCTSRatioMaxCut(gj.fJetCTSRatioMaxCut),
 136   fJetCTSRatioMinCut(gj.fJetCTSRatioMinCut), fJetRatioMaxCut(gj.fJetRatioMaxCut),
 137   fJetRatioMinCut(gj.fJetRatioMinCut),  fNCone(gj.fNCone),
 138   fNPt(gj.fNPt), fCone(gj.fCone), fPtThreshold(gj.fPtThreshold),
 139   fPtJetSelectionCut(gj.fPtJetSelectionCut),
 140   fOutputContainer(gj.fOutputContainer), fAngleMaxParam(gj.fAngleMaxParam),
 141   fSelect(gj.fSelect),  fhChargeRatio(gj.fhChargeRatio), fhPi0Ratio(gj.fhPi0Ratio),
 142   fhDeltaPhiCharge(gj.fhDeltaPhiCharge),  fhDeltaPhiPi0(gj.fhDeltaPhiPi0),
 143   fhDeltaEtaCharge(gj.fhDeltaEtaCharge), fhDeltaEtaPi0(gj.fhDeltaEtaPi0),
 144   fhAnglePair(gj.fhAnglePair), fhAnglePairAccepted(gj.fhAnglePairAccepted),
 145   fhAnglePairNoCut(gj.fhAnglePairNoCut), fhAnglePairLeadingCut(gj.fhAnglePairLeadingCut),
 146   fhAnglePairAngleCut(gj.fhAnglePairAngleCut), fhAnglePairAllCut(gj.fhAnglePairAllCut),
 147   fhAnglePairLeading(gj.fhAnglePairLeading),
 148   fhInvMassPairNoCut(gj.fhInvMassPairNoCut), fhInvMassPairLeadingCut(gj.fhInvMassPairLeadingCut),
 149   fhInvMassPairAngleCut(gj.fhInvMassPairAngleCut),
 150   fhInvMassPairAllCut(gj. fhInvMassPairAllCut), fhInvMassPairLeading(gj.fhInvMassPairLeading),
 151   fhNBkg(gj. fhNBkg), fhNLeading(gj. fhNLeading), fhNJet(gj.fhNJet), fhJetRatio(gj.fhJetRatio), fhJetPt(gj.fhJetPt),
 152   fhBkgRatio (gj.fhBkgRatio), fhBkgPt(gj.fhBkgPt),  fhJetFragment(gj.fhJetFragment), fhBkgFragment(gj.fhBkgFragment),
 153   fhJetPtDist(gj.fhJetPtDist),  fhBkgPtDist(gj.fhBkgPtDist)
 154 {
 155   // cpy ctor
 156   SetName (gj.GetName()) ;
 157   SetTitle(gj.GetTitle()) ;
 158
 159   for(Int_t i = 0; i<10; i++){
 160     fCones[i]        = gj.fCones[i] ;
 161     fNameCones[i]    = gj.fNameCones[i] ;
 162     fPtThres[i]      = gj.fPtThres[i] ;
 163     fNamePtThres[i]  = gj.fNamePtThres[i] ;
 164     if( i < 6 ){
 165       fJetXMin1[i]       = gj.fJetXMin1[i] ;
 166       fJetXMin2[i]       = gj.fJetXMin2[i] ;
 167       fJetXMax1[i]       = gj.fJetXMax1[i] ;
 168       fJetXMax2[i]       = gj.fJetXMax2[i] ;
 169       fBkgMean[i]        = gj.fBkgMean[i] ;
 170       fBkgRMS[i]         = gj.fBkgRMS[i] ;
 171       if( i < 2 ){
 172         fJetE1[i]        = gj.fJetE1[i] ;
 173         fJetE2[i]        = gj.fJetE2[i] ;
 174         fJetSigma1[i]    = gj.fJetSigma1[i] ;
 175         fJetSigma2[i]    = gj.fJetSigma2[i] ;
 176         fPhiEMCALCut[i]  = gj.fPhiEMCALCut[i] ;
 177       }
 178     }
 179   }
 180 }
 181
 182 //_________________________________________________________________________
 183 AliAnaGammaJet & AliAnaGammaJet::operator = (const AliAnaGammaJet & source)
 184 {
 185   //assignment operator
 186   if(&source == this) return *this;
 187
 188   fOutputContainer = source.fOutputContainer ;
 189   fSeveralConeAndPtCuts = source.fSeveralConeAndPtCuts ;
 190   fPbPb = source.fPbPb ; fJetsOnlyInCTS = source.fJetsOnlyInCTS ;
 191   fEtaEMCALCut = source.fEtaEMCALCut ;
 192   fPhiMaxCut = source.fPhiMaxCut ; fPhiMinCut = source.fPhiMinCut ;
 193   fInvMassMaxCut = source.fInvMassMaxCut ; fInvMassMinCut = source.fInvMassMinCut ;
 194   fRatioMaxCut = source.fRatioMaxCut ; fRatioMinCut = source.fRatioMinCut ;
 195   fJetCTSRatioMaxCut = source.fJetCTSRatioMaxCut ;
 196   fJetCTSRatioMinCut = source.fJetCTSRatioMinCut ; fJetRatioMaxCut = source.fJetRatioMaxCut ;
 197   fJetRatioMinCut = source.fJetRatioMinCut ;  fNCone = source.fNCone ;
 198   fNPt = source.fNPt ; fCone = source.fCone ; fPtThreshold = source.fPtThreshold ;
 199   fPtJetSelectionCut = source.fPtJetSelectionCut ;
 200   fAngleMaxParam = source.fAngleMaxParam ;
 201   fSelect = source.fSelect ;  fhChargeRatio = source.fhChargeRatio ; fhPi0Ratio = source.fhPi0Ratio ;
 202   fhDeltaPhiCharge = source.fhDeltaPhiCharge ;  fhDeltaPhiPi0 = source.fhDeltaPhiPi0 ;
 203   fhDeltaEtaCharge = source.fhDeltaEtaCharge ; fhDeltaEtaPi0 = source.fhDeltaEtaPi0 ;
 204   fhAnglePair = source.fhAnglePair ; fhAnglePairAccepted = source.fhAnglePairAccepted ;
 205   fhAnglePairNoCut = source.fhAnglePairNoCut ; fhAnglePairLeadingCut = source.fhAnglePairLeadingCut ;
 206   fhAnglePairAngleCut = source.fhAnglePairAngleCut ; fhAnglePairAllCut = source.fhAnglePairAllCut ;
 207   fhAnglePairLeading = source.fhAnglePairLeading ;
 208   fhInvMassPairNoCut = source.fhInvMassPairNoCut ; fhInvMassPairLeadingCut = source.fhInvMassPairLeadingCut ;
 209   fhInvMassPairAngleCut = source.fhInvMassPairAngleCut ;
 210   fhInvMassPairAllCut = source. fhInvMassPairAllCut ; fhInvMassPairLeading = source.fhInvMassPairLeading ;
 211   fhNBkg = source. fhNBkg ; fhNLeading = source. fhNLeading ; fhNJet = source.fhNJet ; fhJetRatio = source.fhJetRatio ; fhJetPt = source.fhJetPt ;
 212   fhBkgRatio  = source.fhBkgRatio ; fhBkgPt = source.fhBkgPt ;  fhJetFragment = source.fhJetFragment ; fhBkgFragment = source.fhBkgFragment ;
 213   fhJetPtDist = source.fhJetPtDist ;  fhBkgPtDist = source.fhBkgPtDist ;
 214
 215   for(Int_t i = 0; i<10; i++){
 216     fCones[i]        = source.fCones[i] ;
 217     fNameCones[i]    = source.fNameCones[i] ;
 218     fPtThres[i]      = source.fPtThres[i] ;
 219     fNamePtThres[i]  = source.fNamePtThres[i] ;
 220     if( i < 6 ){
 221       fJetXMin1[i]       = source.fJetXMin1[i] ;
 222       fJetXMin2[i]       = source.fJetXMin2[i] ;
 223       fJetXMax1[i]       = source.fJetXMax1[i] ;
 224       fJetXMax2[i]       = source.fJetXMax2[i] ;
 225       fBkgMean[i]        = source.fBkgMean[i] ;
 226       fBkgRMS[i]         = source.fBkgRMS[i] ;
 227       if( i < 2 ){
 228         fJetE1[i]        = source.fJetE1[i] ;
 229         fJetE2[i]        = source.fJetE2[i] ;
 230         fJetSigma1[i]    = source.fJetSigma1[i] ;
 231         fJetSigma2[i]    = source.fJetSigma2[i] ;
 232         fPhiEMCALCut[i]  = source.fPhiEMCALCut[i] ;
 233       }
 234     }
 235   }
 236   return *this;
 237 }
 238
 239 //____________________________________________________________________________
 240 AliAnaGammaJet::~AliAnaGammaJet()
 241 {
 242   // Remove all pointers
 243   fOutputContainer->Clear() ;
 244   delete fOutputContainer ;
 245
 246   delete fhChargeRatio  ;
 247   delete fhPi0Ratio   ;
 248   delete fhDeltaPhiCharge  ;
 249   delete fhDeltaPhiPi0   ;
 250   delete fhDeltaEtaCharge  ;
 251   delete fhDeltaEtaPi0  ;
 252   delete fhAnglePair  ;
 253   delete fhAnglePairAccepted  ;
 254   delete fhAnglePairNoCut  ;
 255   delete fhAnglePairLeadingCut  ;
 256   delete fhAnglePairAngleCut   ;
 257   delete fhAnglePairAllCut   ;
 258   delete fhAnglePairLeading  ;
 259   delete fhInvMassPairNoCut    ;
 260   delete fhInvMassPairLeadingCut  ;
 261   delete fhInvMassPairAngleCut  ;
 262   delete fhInvMassPairAllCut   ;
 263   delete fhInvMassPairLeading  ;
 264   delete fhNBkg   ;
 265   delete fhNLeading  ;
 266   delete fhNJet  ;
 267   delete fhJetRatio  ;
 268   delete fhJetPt   ;
 269   delete fhBkgRatio   ;
 270   delete fhBkgPt  ;
 271   delete fhJetFragment  ;
 272   delete fhBkgFragment  ;
 273   delete fhJetPtDist  ;
 274   delete fhBkgPtDist  ;
 275
 276   delete [] fhJetRatios;
 277   delete [] fhJetPts;
 278   delete [] fhBkgRatios;
 279   delete [] fhBkgPts;
 280
 281   delete [] fhNLeadings;
 282   delete [] fhNJets;
 283   delete [] fhNBkgs;
 284
 285   delete [] fhJetFragments;
 286   delete [] fhBkgFragments;
 287   delete [] fhJetPtDists;
 288   delete [] fhBkgPtDists;
 289
 290 }
 291
 292 //____________________________________________________________________________
 293 Double_t AliAnaGammaJet::CalculateJetRatioLimit(const Double_t ptg,
 294                                                  const Double_t *par,
 295                                                  const Double_t *x) {
 296
 297   //Parametrized cut for the energy of the jet.
 298
 299   Double_t epp = par[0] + par[1] * ptg ;
 300   Double_t spp = par[2] + par[3] * ptg ;
 301   Double_t f   = x[0]   + x[1]   * ptg ;
 302   Double_t epb = epp + par[4] ;
 303   Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
 304   Double_t rat = (epb - spb * f) / ptg ;
 305   //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
 306   //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
 307   return rat ;
 308 }
 309
 310 //______________________________________________________________________________
 311 void AliAnaGammaJet::ConnectInputData(const Option_t*)
 312 {
 313   // Initialisation of branch container and histograms
 314   AliAnaGammaDirect::ConnectInputData("");
 315
 316 }
 317
 318 //________________________________________________________________________
 319 void AliAnaGammaJet::CreateOutputObjects()
 320 {
 321
 322   // Create histograms to be saved in output file and
 323   // stores them in fOutputContainer
 324
 325   AliAnaGammaDirect::CreateOutputObjects();
 326
 327   fOutputContainer = new TObjArray(100) ;
 328
 329   TObjArray  * outputContainer =GetOutputContainer();
 330   for(Int_t i = 0; i < outputContainer->GetEntries(); i++ )
 331     fOutputContainer->Add(outputContainer->At(i)) ;
 332
 333   //
 334   fhChargeRatio  = new TH2F
 335     ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
 336      120,0,120,120,0,1);
 337   fhChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
 338   fhChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
 339   fOutputContainer->Add(fhChargeRatio) ;
 340
 341   fhDeltaPhiCharge  = new TH2F
 342     ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
 343      200,0,120,200,0,6.4);
 344   fhDeltaPhiCharge->SetYTitle("#Delta #phi");
 345   fhDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
 346   fOutputContainer->Add(fhDeltaPhiCharge) ;
 347
 348   fhDeltaEtaCharge  = new TH2F
 349     ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
 350      200,0,120,200,-2,2);
 351   fhDeltaEtaCharge->SetYTitle("#Delta #eta");
 352   fhDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
 353   fOutputContainer->Add(fhDeltaEtaCharge) ;
 354
 355   //
 356   if(!fJetsOnlyInCTS){
 357     fhPi0Ratio  = new TH2F
 358       ("Pi0Ratio","p_{T leading  #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
 359        120,0,120,120,0,1);
 360     fhPi0Ratio->SetYTitle("p_{T lead  #pi^{0}} /p_{T #gamma}");
 361     fhPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
 362     fOutputContainer->Add(fhPi0Ratio) ;
 363
 364     fhDeltaPhiPi0  = new TH2F
 365       ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
 366        200,0,120,200,0,6.4);
 367     fhDeltaPhiPi0->SetYTitle("#Delta #phi");
 368     fhDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
 369     fOutputContainer->Add(fhDeltaPhiPi0) ;
 370
 371     fhDeltaEtaPi0  = new TH2F
 372       ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
 373        200,0,120,200,-2,2);
 374     fhDeltaEtaPi0->SetYTitle("#Delta #eta");
 375     fhDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
 376     fOutputContainer->Add(fhDeltaEtaPi0) ;
 377
 378     //
 379     fhAnglePair  = new TH2F
 380       ("AnglePair",
 381        "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}",
 382        200,0,50,200,0,0.2);
 383     fhAnglePair->SetYTitle("Angle (rad)");
 384     fhAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
 385     fOutputContainer->Add(fhAnglePair) ;
 386
 387     fhAnglePairAccepted  = new TH2F
 388       ("AnglePairAccepted",
 389        "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}, both #gamma in eta<0.7, inside window",
 390        200,0,50,200,0,0.2);
 391     fhAnglePairAccepted->SetYTitle("Angle (rad)");
 392     fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
 393     fOutputContainer->Add(fhAnglePairAccepted) ;
 394
 395     fhAnglePairNoCut  = new TH2F
 396       ("AnglePairNoCut",
 397        "Angle between all #gamma pair vs p_{T  #pi^{0}}",200,0,50,200,0,0.2);
 398     fhAnglePairNoCut->SetYTitle("Angle (rad)");
 399     fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
 400     fOutputContainer->Add(fhAnglePairNoCut) ;
 401
 402     fhAnglePairLeadingCut  = new TH2F
 403       ("AnglePairLeadingCut",
 404        "Angle between all #gamma pair that have a good phi and pt vs p_{T  #pi^{0}}",
 405        200,0,50,200,0,0.2);
 406     fhAnglePairLeadingCut->SetYTitle("Angle (rad)");
 407     fhAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
 408     fOutputContainer->Add(fhAnglePairLeadingCut) ;
 409
 410     fhAnglePairAngleCut  = new TH2F
 411       ("AnglePairAngleCut",
 412        "Angle between all #gamma pair (angle + leading cut) vs p_{T  #pi^{0}}"
 413        ,200,0,50,200,0,0.2);
 414     fhAnglePairAngleCut->SetYTitle("Angle (rad)");
 415     fhAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
 416     fOutputContainer->Add(fhAnglePairAngleCut) ;
 417
 418     fhAnglePairAllCut  = new TH2F
 419       ("AnglePairAllCut",
 420        "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T  #pi^{0}}"
 421        ,200,0,50,200,0,0.2);
 422     fhAnglePairAllCut->SetYTitle("Angle (rad)");
 423     fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
 424     fOutputContainer->Add(fhAnglePairAllCut) ;
 425
 426     fhAnglePairLeading  = new TH2F
 427       ("AnglePairLeading",
 428        "Angle between all #gamma pair finally selected vs p_{T  #pi^{0}}",
 429        200,0,50,200,0,0.2);
 430     fhAnglePairLeading->SetYTitle("Angle (rad)");
 431     fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
 432     fOutputContainer->Add(fhAnglePairLeading) ;
 433
 434     //
 435     fhInvMassPairNoCut  = new TH2F
 436       ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
 437        120,0,120,360,0,0.5);
 438     fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 439     fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
 440     fOutputContainer->Add(fhInvMassPairNoCut) ;
 441
 442     fhInvMassPairLeadingCut  = new TH2F
 443       ("InvMassPairLeadingCut",
 444        "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
 445        120,0,120,360,0,0.5);
 446     fhInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 447     fhInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
 448     fOutputContainer->Add(fhInvMassPairLeadingCut) ;
 449
 450     fhInvMassPairAngleCut  = new TH2F
 451       ("InvMassPairAngleCut",
 452        "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
 453        120,0,120,360,0,0.5);
 454     fhInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 455     fhInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
 456     fOutputContainer->Add(fhInvMassPairAngleCut) ;
 457
 458     fhInvMassPairAllCut  = new TH2F
 459       ("InvMassPairAllCut",
 460        "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
 461        120,0,120,360,0,0.5);
 462     fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
 463     fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
 464     fOutputContainer->Add(fhInvMassPairAllCut) ;
 465
 466     fhInvMassPairLeading  = new TH2F
 467       ("InvMassPairLeading",
 468        "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
 469        120,0,120,360,0,0.5);
 470     fhInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
 471     fhInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
 472     fOutputContainer->Add(fhInvMassPairLeading) ;
 473   }
 474
 475
 476   if(!fSeveralConeAndPtCuts){
 477
 478     //Count
 479     fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000);
 480     fhNBkg->SetYTitle("counts");
 481     fhNBkg->SetXTitle("N");
 482     fOutputContainer->Add(fhNBkg) ;
 483
 484     fhNLeading  = new TH2F
 485       ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120);
 486     fhNLeading->SetYTitle("p_{T charge} (GeV/c)");
 487     fhNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
 488     fOutputContainer->Add(fhNLeading) ;
 489
 490     fhNJet  = new TH1F("NJet","Accepted jets",240,0,120);
 491     fhNJet->SetYTitle("N");
 492     fhNJet->SetXTitle("p_{T #gamma}(GeV/c)");
 493     fOutputContainer->Add(fhNJet) ;
 494
 495     //Ratios and Pt dist of reconstructed (not selected) jets
 496     //Jet
 497     fhJetRatio  = new TH2F
 498       ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
 499        240,0,120,200,0,10);
 500     fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
 501     fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
 502     fOutputContainer->Add(fhJetRatio) ;
 503
 504     fhJetPt  = new TH2F
 505       ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
 506     fhJetPt->SetYTitle("p_{T jet}");
 507     fhJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
 508     fOutputContainer->Add(fhJetPt) ;
 509
 510     //Bkg
 511
 512     fhBkgRatio  = new TH2F
 513       ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
 514        240,0,120,200,0,10);
 515     fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
 516     fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
 517     fOutputContainer->Add(fhBkgRatio) ;
 518
 519     fhBkgPt  = new TH2F
 520       ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
 521     fhBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
 522     fhBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
 523     fOutputContainer->Add(fhBkgPt) ;
 524
 525     //Jet Distributions
 526
 527     fhJetFragment  =
 528       new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
 529                240,0.,120.,1000,0.,1.2);
 530     fhJetFragment->SetYTitle("x_{T}");
 531     fhJetFragment->SetXTitle("p_{T #gamma}");
 532     fOutputContainer->Add(fhJetFragment) ;
 533
 534     fhBkgFragment  = new TH2F
 535       ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
 536        240,0.,120.,1000,0.,1.2);
 537     fhBkgFragment->SetYTitle("x_{T}");
 538     fhBkgFragment->SetXTitle("p_{T #gamma}");
 539     fOutputContainer->Add(fhBkgFragment) ;
 540
 541     fhJetPtDist  =
 542       new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
 543     fhJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
 544     fOutputContainer->Add(fhJetPtDist) ;
 545
 546     fhBkgPtDist  = new TH2F
 547       ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
 548     fhBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
 549     fOutputContainer->Add(fhBkgPtDist) ;
 550
 551   }
 552   else{
 553     //If we want to study the jet for different cones and pt
 554
 555     for(Int_t icone = 0; icone<fNCone; icone++){
 556       for(Int_t ipt = 0; ipt<fNPt;ipt++){
 557
 558         //Jet
 559
 560         fhJetRatios[icone][ipt]  = new TH2F
 561           ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 562            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
 563            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
 564            240,0,120,200,0,10);
 565         fhJetRatios[icone][ipt]->
 566           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
 567         fhJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 568         fOutputContainer->Add(fhJetRatios[icone][ipt]) ;
 569
 570
 571         fhJetPts[icone][ipt]  = new TH2F
 572           ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 573            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
 574            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
 575            240,0,120,400,0,200);
 576         fhJetPts[icone][ipt]->
 577           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
 578         fhJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 579         fOutputContainer->Add(fhJetPts[icone][ipt]) ;
 580
 581         //Bkg
 582         fhBkgRatios[icone][ipt]  = new TH2F
 583           ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 584            "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
 585            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
 586            240,0,120,200,0,10);
 587         fhBkgRatios[icone][ipt]->
 588           SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
 589         fhBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 590         fOutputContainer->Add(fhBkgRatios[icone][ipt]) ;
 591
 592         fhBkgPts[icone][ipt]  = new TH2F
 593           ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 594            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
 595            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
 596            240,0,120,400,0,200);
 597         fhBkgPts[icone][ipt]->
 598           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
 599         fhBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 600         fOutputContainer->Add(fhBkgPts[icone][ipt]) ;
 601
 602         //Counts
 603         fhNBkgs[icone][ipt]  = new TH1F
 604           ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 605            "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
 606            +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
 607         fhNBkgs[icone][ipt]->SetYTitle("counts");
 608         fhNBkgs[icone][ipt]->SetXTitle("N");
 609         fOutputContainer->Add(fhNBkgs[icone][ipt]) ;
 610
 611         fhNLeadings[icone][ipt]  = new TH2F
 612           ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 613            "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>"
 614            +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
 615         fhNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
 616         fhNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
 617         fOutputContainer->Add(fhNLeadings[icone][ipt]) ;
 618
 619         fhNJets[icone][ipt]  = new TH1F
 620           ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 621            "Number of neutral jets, cone ="+fNameCones[icone]+", pt>"
 622            +fNamePtThres[ipt]+" GeV/c",120,0,120);
 623         fhNJets[icone][ipt]->SetYTitle("N");
 624         fhNJets[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
 625         fOutputContainer->Add(fhNJets[icone][ipt]) ;
 626
 627         //Fragmentation Function
 628         fhJetFragments[icone][ipt]  = new TH2F
 629           ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 630            "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
 631            +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
 632         fhJetFragments[icone][ipt]->SetYTitle("x_{T}");
 633         fhJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
 634         fOutputContainer->Add(fhJetFragments[icone][ipt]) ;
 635
 636         fhBkgFragments[icone][ipt]  = new TH2F
 637           ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 638            "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
 639            +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
 640         fhBkgFragments[icone][ipt]->SetYTitle("x_{T}");
 641         fhBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
 642         fOutputContainer->Add(fhBkgFragments[icone][ipt]) ;
 643
 644         //Jet particle distribution
 645
 646         fhJetPtDists[icone][ipt]  = new TH2F
 647           ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 648            "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
 649            " GeV/c",120,0.,120.,120,0.,120.);
 650         fhJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 651         fOutputContainer->Add(fhJetPtDists[icone][ipt]) ;
 652
 653         fhBkgPtDists[icone][ipt]  = new TH2F
 654           ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
 655            "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
 656            " GeV/c",120,0.,120.,120,0.,120.);
 657         fhBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
 658         fOutputContainer->Add(fhBkgPtDists[icone][ipt]) ;
 659
 660       }//ipt
 661     } //icone
 662   }//If we want to study any cone or pt threshold
 663 }
 664
 665
 666 //____________________________________________________________________________
 667 void AliAnaGammaJet::Exec(Option_t *)
 668 {
 669
 670   // Processing of one event
 671
 672   //Get ESDs
 673   Long64_t entry = GetChain()->GetReadEntry() ;
 674
 675   if (!GetESD()) {
 676     AliError("fESD is not connected to the input!") ;
 677     return ;
 678   }
 679
 680   if (GetPrintInfo())
 681     AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(GetChain()))->GetFile()->GetName(), entry)) ;
 682
 683   //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
 684
 685   TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
 686   TClonesArray * plCTS         = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
 687   TClonesArray * plNe         = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter
 688   TClonesArray * plCalo     = new TClonesArray("TParticle",1000);  // All particles measured in Prompt Gamma calorimeter
 689
 690
 691   TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
 692   TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
 693
 694   Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
 695
 696   //Fill lists with photons, neutral particles and charged particles
 697   //look for the highest energy photon in the event inside fCalorimeter
 698   AliDebug(2, "Fill particle lists, get prompt gamma");
 699
 700   //Fill particle lists
 701
 702   if(GetCalorimeter() == "PHOS")
 703     CreateParticleList(particleList, plCTS,plNe,plCalo);
 704   else if(GetCalorimeter() == "EMCAL")
 705     CreateParticleList(particleList, plCTS,plCalo,plNe);
 706   else
 707     AliError("No calorimeter selected");
 708
 709   //Search highest energy prompt gamma in calorimeter
 710   GetPromptGamma(plCalo,  plCTS, pGamma, iIsInPHOSorEMCAL) ;
 711
 712   AliDebug(1, Form("Is Gamma in %s? %d",GetCalorimeter().Data(),iIsInPHOSorEMCAL));
 713   AliDebug(3,Form("Charged list entries %d, Neutral list entries %d, %s list entries %d",
 714                   plCTS->GetEntries(),plNe->GetEntries(), GetCalorimeter().Data(), plCalo->GetEntries()));
 715
 716   //If there is any prompt photon in fCalorimeter,
 717   //search jet leading particle
 718
 719   if(iIsInPHOSorEMCAL){
 720     if (GetPrintInfo())
 721       AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
 722
 723     AliDebug(2, "Get Leading Particles, Make jet");
 724
 725     //Search leading particles in CTS and EMCAL
 726     if(GetLeadingParticle(plCTS, plNe, pGamma, pLeading)){
 727       if (GetPrintInfo())
 728         AliInfo(Form("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ;
 729
 730       //Search Jet
 731       if(!fSeveralConeAndPtCuts)
 732         MakeJet(particleList,pGamma,pLeading,"");
 733       else{
 734         for(Int_t icone = 0; icone<fNCone; icone++) {
 735           for(Int_t ipt = 0; ipt<fNPt;ipt++) {
 736             TString lastname ="Cone"+ fNameCones[icone]+"Pt"+ fNamePtThres[ipt];
 737             MakeJet(particleList, pGamma, pLeading,lastname);
 738           }//icone
 739         }//ipt
 740       }//fSeveralConeAndPtCuts
 741     }//Leading
 742   }//Gamma in Calo
 743
 744   AliDebug(2, "End of analysis, delete pointers");
 745
 746   particleList->Delete() ;
 747   plCTS->Delete() ;
 748   plNe->Delete() ;
 749   plCalo->Delete() ;
 750   pLeading->Delete();
 751   pGamma->Delete();
 752
 753   delete plNe ;
 754   delete plCalo ;
 755   delete plCTS ;
 756   delete particleList ;
 757   //  delete pLeading;
 758   //  delete pGamma;
 759
 760   PostData(0, fOutputContainer);
 761 }
 762
 763 //____________________________________________________________________________
 764 void AliAnaGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname)
 765 {
 766   //Fill histograms wth jet fragmentation
 767   //and number of selected jets and leading particles
 768   //and the background multiplicity
 769   TParticle * particle = 0 ;
 770   Int_t ipr = 0;
 771   Float_t  charge = 0;
 772
 773   TIter next(pl) ;
 774   while ( (particle = dynamic_cast<TParticle*>(next())) ) {
 775     ipr++ ;
 776     Double_t pt = particle->Pt();
 777
 778     charge = TDatabasePDG::Instance()
 779       ->GetParticle(particle->GetPdgCode())->Charge();
 780     if(charge != 0){//Only jet Charged particles
 781       dynamic_cast<TH2F*>
 782         (fOutputContainer->FindObject(type+"Fragment"+lastname))
 783         ->Fill(ptg,pt/ptg);
 784       dynamic_cast<TH2F*>
 785         (fOutputContainer->FindObject(type+"PtDist"+lastname))
 786         ->Fill(ptg,pt);
 787     }//charged
 788
 789   }//while
 790
 791   if(type == "Bkg")
 792     dynamic_cast<TH1F*>
 793       (fOutputContainer->FindObject("NBkg"+lastname))
 794       ->Fill(ipr);
 795   else{
 796     dynamic_cast<TH1F*>
 797       (fOutputContainer->FindObject("NJet"+lastname))->
 798       Fill(ptg);
 799     dynamic_cast<TH2F*>
 800       (fOutputContainer->FindObject("NLeading"+lastname))
 801       ->Fill(ptg,ptl);
 802   }
 803
 804 }
 805
 806 //____________________________________________________________________________
 807 void  AliAnaGammaJet::GetLeadingCharge(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) const
 808 {
 809   //Search for the charged particle with highest with
 810   //Phi=Phi_gamma-Pi and pT=0.1E_gamma
 811   Double_t pt  = -100.;
 812   Double_t phi = -100.;
 813
 814   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
 815
 816     TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
 817
 818     Double_t ptl  = particle->Pt();
 819     Double_t rat  = ptl/pGamma->Pt() ;
 820     Double_t phil = particle->Phi() ;
 821     Double_t phig = pGamma->Phi();
 822
 823     //Selection within angular and energy limits
 824     if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
 825         (rat > fRatioMinCut) && (rat < fRatioMaxCut)  && (ptl  > pt)) {
 826        phi = phil ;
 827        pt  = ptl ;
 828        pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
 829        AliDebug(4,Form("Charge in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, particle->Eta(), phi,rat)) ;
 830      }
 831   }
 832
 833   AliDebug(3,Form("Leading in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, pLeading->Eta(), phi,pt/pGamma->Pt())) ;
 834
 835 }
 836
 837
 838 //____________________________________________________________________________
 839 void  AliAnaGammaJet::GetLeadingPi0(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading)
 840 {
 841
 842   //Search for the neutral pion with highest with
 843   //Phi=Phi_gamma-Pi and pT=0.1E_gamma
 844   Double_t pt  = -100.;
 845   Double_t phi = -100.;
 846   Double_t ptl = -100.;
 847   Double_t rat = -100.;
 848   Double_t phil = -100. ;
 849   Double_t ptg  = pGamma->Pt();
 850   Double_t phig = pGamma->Phi();
 851
 852   TIter next(pl);
 853   TParticle * particle = 0;
 854
 855   Int_t iPrimary = -1;
 856   TLorentzVector gammai,gammaj;
 857   Double_t angle = 0., e = 0., invmass = 0.;
 858   Double_t anglef = 0., ef = 0., invmassf = 0.;
 859   Int_t ksPdg = 0;
 860   Int_t jPrimary=-1;
 861
 862   while ( (particle = (TParticle*)next()) ) {
 863     iPrimary++;
 864
 865     ksPdg = particle->GetPdgCode();
 866     ptl  = particle->Pt();
 867     if(ksPdg == 111){ //2 gamma overlapped, found with PID
 868       rat = ptl/ptg ;
 869       phil = particle->Phi() ;
 870       //Selection within angular and energy limits
 871       if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
 872          ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
 873         phi = phil ;
 874         pt  = ptl ;
 875         pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
 876         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())) ;
 877       }// cuts
 878     }// pdg = 111
 879     else if(ksPdg == 22){//1 gamma
 880       //Search the photon companion in case it comes from  a Pi0 decay
 881       //Apply several cuts to select the good pair
 882       particle->Momentum(gammai);
 883       jPrimary=-1;
 884       TIter next2(pl);
 885       while ( (particle = (TParticle*)next2()) ) {
 886         jPrimary++;
 887         if(jPrimary>iPrimary){
 888           ksPdg = particle->GetPdgCode();
 889
 890           if(ksPdg == 22){
 891             particle->Momentum(gammaj);
 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             invmass = (gammai+gammaj).M();
 900             angle = gammaj.Angle(gammai.Vect());
 901             //Info("GetLeadingPi0","Angle %f", angle);
 902             e = (gammai+gammaj).E();
 903             //Fill histograms with no cuts applied.
 904             fhAnglePairNoCut->Fill(e,angle);
 905             fhInvMassPairNoCut->Fill(ptg,invmass);
 906             //First cut on the energy and azimuth of the pair
 907             if((rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
 908                ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
 909
 910               fhAnglePairLeadingCut->Fill(e,angle);
 911               fhInvMassPairLeadingCut->Fill(ptg,invmass);
 912               //Second cut on the aperture of the pair
 913               if(IsAngleInWindow(angle,e)){
 914                 fhAnglePairAngleCut->Fill(e,angle);
 915                 fhInvMassPairAngleCut->Fill(ptg,invmass);
 916
 917                 //Info("GetLeadingPi0","InvMass %f", invmass);
 918                 //Third cut on the invariant mass of the pair
 919                 if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
 920                   fhInvMassPairAllCut->Fill(ptg,invmass);
 921                   fhAnglePairAllCut->Fill(e,angle);
 922                   if(ptl > pt ){
 923                     pt       = ptl;
 924                     phi      = phil ;
 925                     ef       = e ;
 926                     anglef   = angle ;
 927                     invmassf = invmass ;
 928                     pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
 929                     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())) ;
 930                   }
 931                 }//cuts
 932               }//(invmass>0.125) && (invmass<0.145)
 933             }//gammaj.Angle(gammai.Vect())<0.04
 934           }//if pdg = 22
 935         }//iprimary<jprimary
 936       }//while
 937     }// if pdg = 22
 938     //     }
 939   }//while
 940
 941   if(ef > 0.){//Final pi0 found, highest pair energy, fill histograms
 942     fhInvMassPairLeading->Fill(ptg,invmassf);
 943     fhAnglePairLeading->Fill(ef,anglef);
 944   }
 945
 946   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())) ;
 947 }
 948
 949 //____________________________________________________________________________
 950 Bool_t  AliAnaGammaJet::GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe,
 951                                          TParticle * pGamma, TParticle * pLeading)
 952 {
 953   //Search Charged or Neutral leading particle, select the highest one.
 954
 955   TParticle * pLeadingCh = new TParticle();
 956   TParticle * pLeadingPi0 = new TParticle();
 957
 958   Double_t ptg  =  pGamma->Pt();
 959   Double_t phig = pGamma->Phi();
 960   Double_t etag = pGamma->Eta();
 961
 962   if(GetCalorimeter() == "PHOS" && !fJetsOnlyInCTS)
 963     {
 964       AliDebug(3, "GetLeadingPi0");
 965       GetLeadingPi0   (plNe, pGamma, pLeadingPi0) ;
 966       AliDebug(3, "GetLeadingCharge");
 967       GetLeadingCharge(plCTS, pGamma, pLeadingCh) ;
 968
 969       Double_t ptch = pLeadingCh->Pt();
 970       Double_t phich = pLeadingCh->Phi();
 971       Double_t etach = pLeadingCh->Eta();
 972       Double_t ptpi = pLeadingPi0->Pt();
 973       Double_t phipi = pLeadingPi0->Phi();
 974       Double_t etapi = pLeadingPi0->Eta();
 975
 976       //Is leading cone inside EMCAL acceptance?
 977
 978       Bool_t insidech = kFALSE ;
 979       if((phich - fCone) >  fPhiEMCALCut[0] &&
 980          (phich + fCone) <  fPhiEMCALCut[1] &&
 981         (etach-fCone) < fEtaEMCALCut )
 982         insidech = kTRUE ;
 983
 984       Bool_t insidepi = kFALSE ;
 985       if((phipi - fCone) >  fPhiEMCALCut[0] &&
 986          (phipi + fCone) <  fPhiEMCALCut[1] &&
 987         (etapi-fCone) < fEtaEMCALCut )
 988         insidepi = kTRUE ;
 989
 990       AliDebug(2,Form("Leading:  charged pt %f, pi0 pt  %f",ptch,ptpi)) ;
 991
 992       if (ptch > 0 || ptpi > 0)
 993         {
 994           if(insidech && (ptch > ptpi))
 995             {
 996               if (GetPrintInfo())
 997                 AliInfo("Leading found in CTS");
 998               pLeading->SetMomentum(pLeadingCh->Px(),pLeadingCh->Py(),pLeadingCh->Pz(),pLeadingCh->Energy());
 999               AliDebug(3,Form("Final leading found in CTS, pt %f, phi %f, eta %f",ptch,phich,etach)) ;
1000               fhChargeRatio->Fill(ptg,ptch/pGamma->Pt());
1001               fhDeltaPhiCharge->Fill(ptg,pGamma->Phi()-phich);
1002               fhDeltaEtaCharge->Fill(ptg,pGamma->Eta()-etach);
1003               return 1;
1004             }
1005
1006           else if((ptpi > ptch) && insidepi)
1007             {
1008               if (GetPrintInfo())
1009                 AliInfo("Leading found in EMCAL");
1010               pLeading->SetMomentum(pLeadingPi0->Px(),pLeadingPi0->Py(),pLeadingPi0->Pz(),pLeadingPi0->Energy());
1011               AliDebug(3,Form("Final leading found in EMCAL, pt %f, phi %f, eta %f",ptpi,phipi,etapi)) ;
1012               fhPi0Ratio     ->Fill(ptg,ptpi/ptg);
1013               fhDeltaPhiPi0->Fill(ptg,phig-phipi);
1014               fhDeltaEtaPi0->Fill(ptg,etag-etapi);
1015               return 1;
1016             }
1017
1018           else{
1019             AliDebug(3,"NO LEADING PARTICLE FOUND");}
1020           return 0;
1021         }
1022       else{
1023         AliDebug(3,"NO LEADING PARTICLE FOUND");
1024         return 0;
1025       }
1026     }
1027
1028   else
1029     {
1030       //No calorimeter present for Leading particle detection
1031       GetLeadingCharge(plCTS, pGamma, pLeading) ;
1032       Double_t ptch = pLeading->Pt();
1033       Double_t phich = pLeading->Phi();
1034       Double_t etach = pLeading->Eta();
1035       if(ptch > 0){
1036         fhChargeRatio->Fill(ptg,ptch/ptg);
1037         fhDeltaPhiCharge->Fill(ptg,phig-phich);
1038         fhDeltaEtaCharge->Fill(ptg,etag-etach);
1039         AliDebug(3,Form("Leading found :  pt %f, phi %f, eta %f",ptch,phich,etach)) ;
1040         return 1;
1041       }
1042       else
1043         {
1044           AliDebug(3,"NO LEADING PARTICLE FOUND");
1045           return 0;
1046         }
1047     }
1048 }
1049
1050   //____________________________________________________________________________
1051 void AliAnaGammaJet::InitParameters()
1052 {
1053 //   // Initialisation of branch container
1054   //AliAnaGammaDirect::InitParameters();
1055
1056   //Initialize the parameters of the analysis.
1057   //fCalorimeter="PHOS";
1058   fAngleMaxParam.Set(4) ;
1059   fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
1060   fAngleMaxParam.AddAt(-0.25,1) ;
1061   fAngleMaxParam.AddAt(0.025,2) ;
1062   fAngleMaxParam.AddAt(-2e-4,3) ;
1063   fSeveralConeAndPtCuts   = kFALSE ;
1064   //fPrintInfo           = kTRUE;
1065   fPbPb                = kFALSE ;
1066   fInvMassMaxCut  = 0.16 ;
1067   fInvMassMinCut  = 0.11 ;
1068   //fJetsOnlyInCTS    = kTRUE ;
1069   fEtaEMCALCut     = 0.7 ;
1070   fPhiEMCALCut[0] = 80. *TMath::Pi()/180.;
1071   fPhiEMCALCut[1] = 190.*TMath::Pi()/180.;
1072   fPhiMaxCut      = 3.4 ;
1073   fPhiMinCut      = 2.9 ;
1074
1075   //Jet selection parameters
1076   //Fixed cut (old)
1077   fRatioMaxCut    = 1.0 ;
1078   fRatioMinCut    = 0.1 ;
1079   fJetRatioMaxCut = 1.2 ;
1080   fJetRatioMinCut = 0.3 ;
1081   fJetsOnlyInCTS = kFALSE ;
1082   fJetCTSRatioMaxCut = 1.2 ;
1083   fJetCTSRatioMinCut = 0.3 ;
1084   fSelect         = 0  ;
1085
1086   //Cut depending on gamma energy
1087
1088   fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applied
1089   //Reconstructed jet energy dependence parameters
1090   //e_jet = a1+e_gamma b2.
1091   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1092   fJetE1[0] = -5.75; fJetE1[1] = -4.1;
1093   fJetE2[0] = 1.005; fJetE2[1] = 1.05;
1094
1095   //Reconstructed sigma of jet energy dependence parameters
1096   //s_jet = a1+e_gamma b2.
1097   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1098   fJetSigma1[0] = 2.65;   fJetSigma1[1] = 2.75;
1099   fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
1100
1101   //Background mean energy and RMS
1102   //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1103   //Index 2-> (low pt jets)BKG > 0.5 GeV;
1104   //Index > 2, same for CTS conf
1105   fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
1106   fBkgMean[3] = 0.; fBkgMean[4] = 6.4;  fBkgMean[5] = 48.6;
1107   fBkgRMS[0]  = 0.; fBkgRMS[1]  = 7.5;  fBkgRMS[2]  = 22.0;
1108   fBkgRMS[3]  = 0.; fBkgRMS[4]  = 5.4;  fBkgRMS[5]  = 13.2;
1109
1110   //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
1111   //limits for monoenergetic jets.
1112   //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1113   //Index 2-> (low pt jets) BKG > 0.5 GeV;
1114   //Index > 2, same for CTS conf
1115
1116   fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
1117   fJetXMin1[3] =-2.0  ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1  ;
1118   fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
1119   fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
1120   fJetXMax1[0] =-3.8  ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6  ;
1121   fJetXMax1[3] =-2.7  ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7  ;
1122   fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
1123   fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
1124
1125
1126   //Different cones and pt thresholds to construct the jet
1127
1128   fCone        = 0.3  ;
1129   fPtThreshold = 0.   ;
1130   fNCone       = 3    ;
1131   fNPt         = 3    ;
1132   fCones[1]    = 0.2  ; fNameCones[1]   = "02" ;
1133   fPtThres[0]  = 0.0  ; fNamePtThres[0] = "00" ;
1134   fCones[0]    = 0.3  ; fNameCones[0]   = "03" ;
1135   fPtThres[1]  = 0.5  ; fNamePtThres[1] = "05" ;
1136   fCones[2]    = 0.4  ; fNameCones[2]   = "04" ;
1137   fPtThres[2]  = 1.0  ; fNamePtThres[2] = "10" ;
1138   //Fill particle lists when PID is ok
1139   // fEMCALPID = kFALSE;
1140   // fPHOSPID = kFALSE;
1141
1142 }
1143
1144 //__________________________________________________________________________-
1145 Bool_t AliAnaGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
1146   //Check if the opening angle of the candidate pairs is inside
1147   //our selection windowd
1148
1149   Bool_t result = kFALSE;
1150   Double_t mpi0 = 0.1349766;
1151   Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
1152     +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
1153   Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
1154   Double_t min = 100. ;
1155   if(arg>0.)
1156     min = TMath::ACos(arg);
1157
1158   if((angle<max)&&(angle>=min))
1159     result = kTRUE;
1160
1161   return result;
1162 }
1163
1164 //__________________________________________________________________________-
1165 Bool_t AliAnaGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj){
1166   //Check if the energy of the reconstructed jet is within an energy window
1167
1168   Double_t par[6];
1169   Double_t xmax[2];
1170   Double_t xmin[2];
1171
1172   Int_t iCTS = 0;
1173   if(fJetsOnlyInCTS)
1174     iCTS = 3 ;
1175
1176   if(!fPbPb){
1177     //Phythia alone, jets with pt_th > 0.2, r = 0.3
1178     par[0] = fJetE1[0]; par[1] = fJetE2[0];
1179     //Energy of the jet peak
1180     //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1181     par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1182     //Sigma  of the jet peak
1183     //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1184     par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
1185     //Parameters reserved for PbPb bkg.
1186     xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
1187     xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
1188     //Factor that multiplies sigma to obtain the best limits,
1189     //by observation, of mono jet ratios (ptjet/ptg)
1190     //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1191
1192   }
1193   else{
1194     if(ptg > fPtJetSelectionCut){
1195       //Phythia +PbPb with  pt_th > 2 GeV/c, r = 0.3
1196       par[0] = fJetE1[0]; par[1] = fJetE2[0];
1197       //Energy of the jet peak, same as in pp
1198       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1199       par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1200       //Sigma  of the jet peak, same as in pp
1201       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1202       par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
1203       //Mean value and RMS of PbPb Bkg
1204       xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
1205       xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
1206       //Factor that multiplies sigma to obtain the best limits,
1207       //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg,
1208       //pt_th > 2 GeV, r = 0.3
1209       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1210
1211     }
1212     else{
1213       //Phythia + PbPb with  pt_th > 0.5 GeV/c, r = 0.3
1214       par[0] = fJetE1[1]; par[1] = fJetE2[1];
1215       //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
1216       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1217       par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
1218       //Sigma  of the jet peak, pt_th > 2 GeV/c, r = 0.3
1219       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1220       par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
1221       //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
1222       xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
1223       xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
1224       //Factor that multiplies sigma to obtain the best limits,
1225       //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg,
1226       //pt_th > 2 GeV, r = 0.3
1227       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1228
1229     }//If low pt jet in bkg
1230   }//if Bkg
1231
1232  //Calculate minimum and maximum limits of the jet ratio.
1233   Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
1234   Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
1235
1236   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));
1237
1238   if(( min < ptj/ptg ) && ( max > ptj/ptg))
1239     return kTRUE;
1240   else
1241     return kFALSE;
1242
1243 }
1244
1245 //____________________________________________________________________________
1246 void AliAnaGammaJet::MakeJet(TClonesArray * pl, TParticle * pGamma, TParticle* pLeading,TString lastname)
1247 {
1248   //Fill the jet with the particles around the leading particle with
1249   //R=fCone and pt_th = fPtThres. Caclulate the energy of the jet and
1250   //check if we select it. Fill jet histograms
1251
1252   TClonesArray * jetList = new TClonesArray("TParticle",1000);
1253   TClonesArray * bkgList = new TClonesArray("TParticle",1000);
1254
1255   TLorentzVector jet   (0,0,0,0);
1256   TLorentzVector bkg(0,0,0,0);
1257   TLorentzVector lv (0,0,0,0);
1258
1259   Double_t ptjet = 0.0;
1260   Double_t ptbkg = 0.0;
1261   Int_t n0 = 0;
1262   Int_t n1 = 0;
1263   Bool_t b1 = kFALSE;
1264   Bool_t b0 = kFALSE;
1265
1266   Double_t ptg  = pGamma->Pt();
1267   Double_t phig = pGamma->Phi();
1268   Double_t ptl  = pLeading->Pt();
1269   Double_t phil = pLeading->Phi();
1270   Double_t etal = pLeading->Eta();
1271
1272   Float_t ptcut = 0. ;
1273   if(fPbPb){
1274     if(ptg > fPtJetSelectionCut)  ptcut = 2. ;
1275     else                          ptcut = 0.5;
1276   }
1277
1278   TIter next(pl) ;
1279   TParticle * particle = 0 ;
1280   while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1281
1282     b0 = kFALSE;
1283     b1 = kFALSE;
1284
1285     //Particles in jet
1286     SetJet(particle, b0, fCone, etal, phil) ;
1287
1288     if(b0){
1289       new((*jetList)[n0++]) TParticle(*particle) ;
1290       particle->Momentum(lv);
1291       if(particle->Pt() > ptcut ){
1292         jet+=lv;
1293         ptjet+=particle->Pt();
1294       }
1295     }
1296
1297     //Background around (phi_gamma-pi, eta_leading)
1298     SetJet(particle, b1, fCone,etal, phig) ;
1299
1300     if(b1) {
1301       new((*bkgList)[n1++]) TParticle(*particle) ;
1302       particle->Momentum(lv);
1303       if(particle->Pt() > ptcut ){
1304         bkg+=lv;
1305         ptbkg+=particle->Pt();
1306       }
1307     }
1308   }
1309
1310   ptjet = jet.Pt();
1311   ptbkg = bkg.Pt();
1312
1313   if(ptjet > 0.) {
1314
1315     AliDebug(2,Form("Gamma   pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg));
1316
1317     //Fill histograms
1318
1319     Double_t ratjet   = ptjet/ptg ;
1320     Double_t ratbkg  = ptbkg/ptg ;
1321
1322     dynamic_cast<TH2F*>
1323       (fOutputContainer->FindObject("JetRatio"+lastname))
1324       ->Fill(ptg,ratjet);
1325     dynamic_cast<TH2F*>
1326       (fOutputContainer->FindObject("JetPt"+lastname))
1327       ->Fill(ptg,ptjet);
1328
1329     dynamic_cast<TH2F*>
1330       (fOutputContainer->FindObject("BkgRatio"+lastname))
1331       ->Fill(ptg,ratbkg);
1332
1333     dynamic_cast<TH2F*>
1334       (fOutputContainer->FindObject("BkgPt"+lastname))
1335       ->Fill(ptg,ptbkg);
1336
1337
1338     //Jet selection
1339     Bool_t kSelect = kFALSE;
1340     if(fSelect == 0)
1341       kSelect = kTRUE; //Accept all jets, no restriction
1342     else if(fSelect == 1){
1343       //Selection with parametrized cuts
1344       if(IsJetSelected(ptg,ptjet))   kSelect = kTRUE;
1345     }
1346     else if(fSelect == 2){
1347       //Simple selection
1348       if(!fJetsOnlyInCTS){
1349         if((ratjet <  fJetRatioMaxCut) && (ratjet > fJetRatioMinCut )) kSelect = kTRUE;
1350       }
1351       else{
1352         if((ratjet <  fJetCTSRatioMaxCut) && (ratjet > fJetCTSRatioMinCut )) kSelect = kTRUE;
1353       }
1354     }
1355     else
1356       AliError("Jet selection option larger than 2, DONT SELECT JETS");
1357
1358
1359     if(kSelect){
1360       if (GetPrintInfo())
1361         AliInfo(Form("Jet Selected: pt %f ", ptjet)) ;
1362
1363       FillJetHistos(jetList, ptg, ptl,"Jet",lastname);
1364       FillJetHistos(bkgList, ptg, ptl, "Bkg",lastname);
1365     }
1366   } //ptjet > 0
1367
1368   jetList ->Delete();
1369   bkgList ->Delete();
1370
1371 }
1372
1373 //____________________________________________________________________________
1374 void AliAnaGammaJet::Print(const Option_t * opt) const
1375 {
1376
1377   //Print some relevant parameters set for the analysis
1378   if(! opt)
1379     return;
1380
1381   Info("Print", "%s %s", GetName(), GetTitle() ) ;
1382   if(!fJetsOnlyInCTS)
1383     printf("EMCAL Acceptance cut  : phi min %f, phi max %f, eta %f \n", fPhiEMCALCut[0],  fPhiEMCALCut[1], fEtaEMCALCut) ;
1384
1385   printf("Phi_Leading                                <     %f\n", fPhiMaxCut) ;
1386   printf("Phi_Leading                                >     %f\n", fPhiMinCut) ;
1387   printf("pT Leading / pT Gamma             <     %f\n", fRatioMaxCut) ;
1388   printf("pT Leading / pT Gamma             >     %f\n", fRatioMinCut) ;
1389   printf("M_pair                                        <     %f\n", fInvMassMaxCut) ;
1390   printf("M_pair                                        >     %f\n", fInvMassMinCut) ;
1391   if(fSelect == 2){
1392     printf("pT Jet / pT Gamma                     <    %f\n", fJetRatioMaxCut) ;
1393     printf("pT Jet / pT Gamma                     >    %f\n", fJetRatioMinCut) ;
1394     printf("pT Jet (Only CTS)/ pT Gamma   <    %f\n", fJetCTSRatioMaxCut) ;
1395     printf("pT Jet (Only CTS)/ pT Gamma   >    %f\n", fJetCTSRatioMinCut) ;
1396   }
1397
1398
1399 }
1400
1401 //___________________________________________________________________
1402 void AliAnaGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone,
1403                              Double_t eta, Double_t phi)
1404 {
1405
1406   //Check if the particle is inside the cone defined by the leading particle
1407   b = kFALSE;
1408
1409   if(phi > TMath::TwoPi())
1410     phi-=TMath::TwoPi();
1411   if(phi < 0.)
1412     phi+=TMath::TwoPi();
1413
1414   Double_t  rad = 10000 + cone;
1415
1416   if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
1417     rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1418                       TMath::Power(part->Phi()-phi,2));
1419   else{
1420     if(part->Phi()-phi > TMath::TwoPi() - cone)
1421       rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1422                         TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
1423     if(part->Phi()-phi < -(TMath::TwoPi() - cone))
1424       rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1425                         TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
1426   }
1427
1428   if(rad < cone )
1429     b = kTRUE;
1430
1431 }
1432
1433
1434 void AliAnaGammaJet::Terminate(Option_t *)
1435 {
1436    // The Terminate() function is the last function to be called during
1437    // a query. It always runs on the client, it can be used to present
1438    // the results graphically or save the results to file.
1439
1440
1441 }