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