PHOS/AliPHOSGammaJet.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.10  2006/01/23 18:04:08  hristov
  21  * Removing meaningless const
  22  *
  23  * Revision 1.9  2006/01/12 16:23:26  schutz
  24  * ESD is properly read with methods of macros/AliReadESD.C copied in it
  25  *
  26  * Revision 1.8  2005/12/20 07:08:32  schutz
  27  * corrected error in call AliReadESD
  28  *
  29  * Revision 1.6  2005/05/28 14:19:04  schutz
  30  * Compilation warnings fixed by T.P.
  31  *
  32  */
  33
  34 //_________________________________________________________________________
  35 // Class for the analysis of gamma-jet correlations
  36 //  Basically it seaches for a prompt photon in the PHOS acceptance,
  37 //  if so we construct a jet around the highest pt particle in the opposite
  38 //  side in azimuth, inside the TPC and EMCAL acceptances. First the leading
  39 //  particle and then the jet have to fullfill several conditions
  40 //  (energy, direction ..) to be accepted. Then the fragmentation function
  41 //  of this jet is constructed
  42 //
  43 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
  44 //////////////////////////////////////////////////////////////////////////////
  45
  46
  47 // --- ROOT system ---
  48
  49 #include <TFile.h>
  50 #include <TParticle.h>
  51 #include <TCanvas.h>
  52 #include <TPaveLabel.h>
  53 #include <TPad.h>
  54 #include <TH2.h>
  55
  56 #include "AliPHOSGammaJet.h"
  57 #include "AliPHOSGetter.h"
  58 #include "AliPHOSGeometry.h"
  59 #include "AliPHOSFastGlobalReconstruction.h"
  60 #include "AliESD.h"
  61 #include "AliESDtrack.h"
  62 #include "Riostream.h"
  63
  64
  65 ClassImp(AliPHOSGammaJet)
  66
  67 //____________________________________________________________________________
  68 AliPHOSGammaJet::AliPHOSGammaJet() : TTask()
  69 {
  70   // ctor
  71   fAngleMaxParam.Set(4) ;
  72   fAngleMaxParam.Reset(0.);
  73   fAnyConeOrPt      = 0  ;
  74   fEtaCut           = 0. ;
  75   fESDdata          = 0  ;
  76   fFastRec          = 0  ;
  77   fInvMassMaxCut    = 0. ;
  78   fInvMassMinCut    = 0. ;
  79   fJetRatioMaxCut   = 0. ;
  80   fJetRatioMinCut   = 0. ;
  81   fJetTPCRatioMaxCut   = 0. ;
  82   fJetTPCRatioMinCut   = 0. ;
  83   fMinDistance      = 0. ;
  84   fNEvent           = 0  ;
  85   fNCone            = 0  ;
  86   fNPt              = 0  ;
  87   fOnlyCharged      = 0  ;
  88   fOptFast          = 0  ;
  89   fPhiMaxCut        = 0. ;
  90   fPhiMinCut        = 0. ;
  91   fPtCut            = 0. ;
  92   fNeutralPtCut     = 0. ;
  93   fChargedPtCut     = 0. ;
  94   fRatioMaxCut      = 0. ;
  95   fRatioMinCut      = 0. ;
  96   fCone             = 0  ;
  97   fPtThreshold      = 0  ;
  98   fTPCCutsLikeEMCAL = 0  ;
  99   fSelect           = 0  ;
 100
 101   fDirName          = "" ;
 102   fESDTree          = "" ;
 103   fPattern          = "" ;
 104
 105   for(Int_t i = 0; i<10; i++){
 106     fCones[i]         = 0.0 ;
 107     fNameCones[i]     = ""  ;
 108     fPtThres[i]      = 0.0 ;
 109     fPtJetSelectionCut = 0.0 ;
 110     fNamePtThres[i]  = ""  ;
 111     if( i < 6 ){
 112       fJetXMin1[i]     = 0.0 ;
 113       fJetXMin2[i]     = 0.0 ;
 114       fJetXMax1[i]     = 0.0 ;
 115       fJetXMax2[i]     = 0.0 ;
 116       fBkgMean[i]      = 0.0 ;
 117       fBkgRMS[i]       = 0.0 ;
 118       if( i < 2 ){
 119         fJetE1[i]        = 0.0 ;
 120         fJetE2[i]        = 0.0 ;
 121         fJetSigma1[i]    = 0.0 ;
 122         fJetSigma2[i]    = 0.0 ;
 123         fPhiEMCALCut[i]  = 0.0 ;
 124       }
 125     }
 126   }
 127
 128   fOptionGJ       = "" ;
 129   fOutputFile     = new TFile(gDirectory->GetName()) ;
 130   fInputFileName  = gDirectory->GetName() ;
 131   fOutputFileName = gDirectory->GetName() ;
 132   fHIJINGFileName = gDirectory->GetName() ;
 133   fHIJING        = 0 ;
 134   fPosParaA      = 0. ;
 135   fPosParaB      = 0. ;
 136   fRan           = 0 ;
 137   fResPara1      = 0. ;
 138   fResPara2      = 0. ;
 139   fResPara3      = 0. ;
 140
 141   fListHistos     = new TObjArray(100) ;
 142   TList * list = gDirectory->GetListOfKeys() ;
 143   TIter next(list) ;
 144   TH2F * h = 0 ;
 145   Int_t index ;
 146   for (index = 0 ; index < list->GetSize()-1 ; index++) {
 147     //-1 to avoid GammaJet Task
 148     h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ;
 149     fListHistos->Add(h) ;
 150   }
 151   List() ;
 152 }
 153
 154 //____________________________________________________________________________
 155 AliPHOSGammaJet::AliPHOSGammaJet(const TString inputfilename) :
 156   TTask("GammaJet","Analysis of gamma-jet correlations")
 157 {
 158   // ctor
 159   fInputFileName = inputfilename;
 160   fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);
 161   AliPHOSGetter *  gime = AliPHOSGetter::Instance(fInputFileName) ;
 162   fNEvent = gime->MaxEvent();
 163   InitParameters();
 164   fListHistos = 0 ;
 165 }
 166
 167 //____________________________________________________________________________
 168 AliPHOSGammaJet::AliPHOSGammaJet(const AliPHOSGammaJet & gj) : TTask(gj)
 169 {
 170   // cpy ctor
 171   fAngleMaxParam     = gj.fAngleMaxParam;
 172   fAnyConeOrPt       = gj.fAnyConeOrPt;
 173   fESDdata           = gj.fESDdata;
 174   fEtaCut            = gj.fEtaCut ;
 175   fInvMassMaxCut     = gj.fInvMassMaxCut ;
 176   fInvMassMinCut     = gj.fInvMassMinCut ;
 177   fFastRec           = gj.fFastRec ;
 178   fOptionGJ          = gj.fOptionGJ ;
 179   fMinDistance       = gj.fMinDistance ;
 180   fOptFast           = gj.fOptFast ;
 181   fOnlyCharged       = gj.fOnlyCharged ;
 182   fOutputFile        = gj.fOutputFile ;
 183   fInputFileName     = gj.fInputFileName ;
 184   fOutputFileName    = gj.fOutputFileName ;
 185   fHIJINGFileName    = gj.fHIJINGFileName ;
 186   fHIJING            = gj.fHIJING ;
 187   fRatioMaxCut       = gj.fRatioMaxCut ;
 188   fRatioMinCut       = gj.fRatioMinCut ;
 189   fJetRatioMaxCut    = gj.fJetRatioMaxCut ;
 190   fJetRatioMinCut    = gj.fJetRatioMinCut ;
 191   fJetTPCRatioMaxCut = gj.fJetRatioMaxCut ;
 192   fJetTPCRatioMinCut = gj.fJetRatioMinCut ;
 193   fNEvent            = gj.fNEvent ;
 194   fNCone             = gj.fNCone ;
 195   fNPt               = gj.fNPt ;
 196   fResPara1          = gj.fResPara1 ;
 197   fResPara2          = gj.fResPara2 ;
 198   fResPara3          = gj.fResPara3 ;
 199   fPtCut             = gj.fPtCut ;
 200   fNeutralPtCut      = gj.fNeutralPtCut ;
 201   fChargedPtCut      = gj.fChargedPtCut ;
 202   fPtJetSelectionCut = gj.fPtJetSelectionCut ;
 203   fPhiMaxCut         = gj.fPhiMaxCut  ;
 204   fPhiMinCut         = gj.fPhiMinCut ;
 205   fPosParaA          = gj.fPosParaA ;
 206   fPosParaB          = gj.fPosParaB ;
 207   fSelect            = gj.fSelect   ;
 208   fTPCCutsLikeEMCAL  = gj.fTPCCutsLikeEMCAL ;
 209   fCone              = gj.fCone ;
 210   fPtThreshold       = gj.fPtThreshold  ;
 211
 212   fDirName           = gj.fDirName ;
 213   fESDTree           = gj.fESDTree ;
 214   fPattern           = gj.fPattern ;
 215
 216   SetName (gj.GetName()) ;
 217   SetTitle(gj.GetTitle()) ;
 218
 219   for(Int_t i = 0; i<10; i++){
 220     fCones[i]        = gj.fCones[i] ;
 221     fNameCones[i]    = gj.fNameCones[i] ;
 222     fPtThres[i]      = gj.fPtThres[i] ;
 223     fNamePtThres[i]  = gj.fNamePtThres[i] ;
 224     if( i < 6 ){
 225       fJetXMin1[i]       = gj.fJetXMin1[i] ;
 226       fJetXMin2[i]       = gj.fJetXMin2[i] ;
 227       fJetXMax1[i]       = gj.fJetXMax1[i] ;
 228       fJetXMax2[i]       = gj.fJetXMax2[i] ;
 229       fBkgMean[i]        = gj.fBkgMean[i] ;
 230       fBkgRMS[i]         = gj.fBkgRMS[i] ;
 231       if( i < 2 ){
 232         fJetE1[i]        = gj.fJetE1[i] ;
 233         fJetE2[i]        = gj.fJetE2[i] ;
 234         fJetSigma1[i]    = gj.fJetSigma1[i] ;
 235         fJetSigma2[i]    = gj.fJetSigma2[i] ;
 236         fPhiEMCALCut[i]  = gj.fPhiEMCALCut[i] ;
 237       }
 238     }
 239   }
 240 }
 241
 242 //____________________________________________________________________________
 243 AliPHOSGammaJet::~AliPHOSGammaJet()
 244 {
 245   fOutputFile->Close() ;
 246
 247 }
 248
 249 //____________________________________________________________________________
 250 void AliPHOSGammaJet::AddHIJINGToList(Int_t iEvent, TClonesArray * particleList,
 251                                       TClonesArray * plCh,
 252                                       TClonesArray * plNe,
 253                                       TClonesArray * plNePHOS)
 254 {
 255
 256   // List of particles copied to a root file.
 257
 258 //   Char_t sb[] = "bgrd/";
 259 //   //  cout<<sb<<endl;
 260 //   Char_t si[10];
 261 //   Char_t sf[] = "/list.root";
 262 //   //  cout<<sf<<endl;
 263 //   sprintf(si,"%d",iEvent);
 264 //   strcat(sb,si) ;
 265 //   //  cout<<si<<endl;
 266 //   strcat(sb,sf) ;
 267 //   //  cout<<si<<endl;
 268 //   TFile * f = TFile::Open(sb) ;
 269 //   //cout<<f->GetName()<<endl;
 270
 271   Char_t fi[100];
 272   sprintf(fi,"bgrd/%d/list.root",iEvent);
 273   TFile * f = TFile::Open(fi) ;
 274   //cout<<f->GetName()<<endl;
 275
 276   TParticle *particle = new TParticle();
 277   TTree *t = (TTree*) f->Get("T");
 278   TBranch *branch = t->GetBranch("primaries");
 279   branch->SetAddress(&particle);
 280
 281   Int_t index   = particleList->GetEntries() ;
 282   Int_t indexNe = plNe->GetEntries() ;
 283   Int_t indexCh = plCh->GetEntries() ;
 284   Int_t indexNePHOS = plNePHOS->GetEntries() ;
 285   Double_t charge = 0.;
 286   Int_t iParticle = 0 ;
 287   Int_t m = 0;
 288   Double_t x = 0., z = 0.;
 289   //  cout<<"bkg entries "<<t->GetEntries()<<endl;
 290
 291   AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
 292   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 293
 294   if(!fOptFast){
 295     for (iParticle=0 ; iParticle < t->GetEntries() ; iParticle++) {
 296       t->GetEvent(iParticle) ;
 297
 298       m = 0 ;
 299       x = 0. ;
 300       z = 0. ;
 301
 302       charge = TDatabasePDG::Instance()
 303         ->GetParticle(particle->GetPdgCode())->Charge();
 304
 305       if((charge != 0) && (particle->Pt() > fChargedPtCut)){
 306         if(TMath::Abs(particle->Eta())<fEtaCut){
 307           new((*particleList)[index]) TParticle(*particle) ;
 308           (dynamic_cast<TParticle*>(particleList->At(index)))
 309             ->SetStatusCode(0) ;
 310           index++ ;
 311
 312           new((*plCh)[indexCh])       TParticle(*particle) ;
 313           (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
 314           indexCh++ ;
 315
 316           if(strstr(fOptionGJ,"deb all")||strstr(fOptionGJ,"deb")){
 317             dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 318               Fill(particle->Pt());
 319           }
 320         }
 321         else if((charge == 0) && (particle->Pt() > fNeutralPtCut) ){
 322           geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
 323
 324           if(m != 0)
 325             {//Is in PHOS
 326               if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
 327                 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 328                   Fill(particle->Pt());
 329
 330               new((*plNePHOS)[indexNePHOS])       TParticle(*particle) ;
 331               (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
 332               indexNePHOS++ ;
 333             }
 334
 335           if((particle->Phi()>fPhiEMCALCut[0]) &&
 336              (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
 337             {//Is in EMCAL
 338               if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
 339                 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 340                   Fill(particle->Pt());
 341
 342               new((*particleList)[index]) TParticle(*particle) ;
 343               (dynamic_cast<TParticle*>(particleList->At(index)))
 344                 ->SetStatusCode(0) ;
 345               index++ ;
 346               new((*plNe)[indexNe])       TParticle(*particle) ;
 347               (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
 348               indexNe++ ;
 349             }
 350         }
 351       }
 352     }
 353   } //No OptFast
 354   else{
 355
 356     Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass();
 357     TLorentzVector pPi0, pGamma1, pGamma2 ;
 358     Double_t angle = 0, cellDistance = 0.;
 359     Bool_t p1 = kFALSE;
 360
 361     //     fFastRec = new AliPHOSFastGlobalReconstruction(fHIJINGFileName);
 362     //     fFastRec->FastReconstruction(iiEvent);
 363
 364     for (iParticle=0 ; iParticle <   t->GetEntries() ; iParticle++) {
 365       t->GetEvent(iParticle) ;
 366       m = 0 ;
 367       x = 0. ;
 368       z = 0. ;
 369       charge = TDatabasePDG::Instance()
 370         ->GetParticle(particle->GetPdgCode())->Charge();
 371
 372       if((charge != 0) && (particle->Pt() > fChargedPtCut)
 373          && (TMath::Abs(particle->Eta())<fEtaCut)){
 374
 375         new((*particleList)[index]) TParticle(*particle) ;
 376         (dynamic_cast<TParticle*>(particleList->At(index)))
 377           ->SetStatusCode(0) ;
 378         index++ ;
 379
 380         new((*plCh)[indexCh])       TParticle(*particle) ;
 381         (dynamic_cast<TParticle*>(plCh->At(indexCh)))->SetStatusCode(0) ;
 382         indexCh++ ;
 383       }
 384       else if(charge == 0){
 385         geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
 386         if((particle->GetPdgCode() != 111) && (particle->Pt() > fNeutralPtCut) &&
 387            (TMath::Abs(particle->Eta())<fEtaCut) ){
 388           TLorentzVector part(particle->Px(),particle->Py(),
 389                               particle->Pz(),particle->Energy());
 390           MakePhoton(part);
 391           if(part.Pt() > fNeutralPtCut){
 392
 393             if(particle->Phi()>fPhiEMCALCut[0] &&
 394                particle->Phi()<fPhiEMCALCut[1] && m == 0)
 395               {
 396                 new((*particleList)[index]) TParticle(*particle) ;
 397                 (dynamic_cast<TParticle*>(particleList->At(index)))
 398                   ->SetStatusCode(0) ;
 399                 (dynamic_cast<TParticle*>(particleList->At(index)))
 400                   ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 401                 index++ ;
 402
 403
 404                 new((*plNe)[indexNe])       TParticle(*particle) ;
 405                 (dynamic_cast<TParticle*>(plNe->At(indexNe)))
 406                   ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 407                 (dynamic_cast<TParticle*>(plNe->At(indexNe)))->SetStatusCode(0) ;
 408                 indexNe++ ;
 409               }
 410             if(m != 0)
 411               {
 412                 new((*plNePHOS)[indexNePHOS])       TParticle(*particle) ;
 413                 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
 414                   ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 415                 (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
 416                 indexNePHOS++ ;
 417               }
 418           }
 419         }
 420         if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
 421            (TMath::Abs(particle->Eta())<fEtaCut+1))
 422           {
 423
 424             pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),
 425                             particle->Energy());
 426
 427             //Decay
 428
 429             Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle);
 430             //Check if decay photons are too close for PHOS
 431             cellDistance = angle*460; //cm
 432             if (cellDistance < fMinDistance) {
 433               if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
 434                 dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 435                   Fill(particle->Pt());
 436
 437               //Pi0 inside phi EMCAL acceptance
 438
 439               TLorentzVector part(particle->Px(),particle->Py(),
 440                                   particle->Pz(),particle->Energy());
 441               MakePhoton(part);
 442               if(part.Pt() > fNeutralPtCut){
 443                 if(particle->Phi()>fPhiEMCALCut[0] &&
 444                    particle->Phi()<fPhiEMCALCut[1] && m == 0){
 445
 446                   new((*particleList)[index]) TParticle(*particle) ;
 447                   (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
 448                   (dynamic_cast<TParticle*>(particleList->At(index)))
 449                     ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 450                   index++ ;
 451
 452                   new((*plNe)[indexNe])       TParticle(*particle) ;
 453                   (dynamic_cast<TParticle*>(plNe->At(indexNe)))      ->SetStatusCode(0) ;
 454                   (dynamic_cast<TParticle*>(plNe->At(indexNe)))
 455                     ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 456                   indexNe++ ;
 457                 }
 458                 if(m != 0){
 459                   if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
 460                     dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 461                       Fill(particle->Pt());
 462                   new((*plNePHOS)[indexNePHOS])       TParticle(*particle) ;
 463                   (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
 464                   (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
 465                     ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 466                   indexNePHOS++;
 467                 }//In PHOS
 468               }
 469             }// if cell<distance
 470             else {
 471
 472               dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
 473                 ->Fill(pPi0.E(),angle);
 474
 475               p1 = kFALSE;
 476               if(pGamma1.Pt() > 0. && TMath::Abs(pGamma1.Eta())<fEtaCut){
 477
 478                 MakePhoton(pGamma1);
 479
 480                 if(pGamma1.Pt() > fNeutralPtCut ){
 481
 482                   TParticle * photon1 =
 483                     new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
 484                                   pGamma1.Pz(),pGamma1.E(),0,0,0,0);
 485                   geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x);
 486                   if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
 487                       && m == 0){
 488                     if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
 489                       dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 490                         Fill(photon1->Pt());
 491                     new((*particleList)[index]) TParticle(*photon1) ;
 492                     (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
 493                     index++ ;
 494
 495                     new((*plNe)[indexNe])       TParticle(*photon1) ;
 496                     (dynamic_cast<TParticle*>(plNe->At(indexNe)))      ->SetStatusCode(0) ;
 497                     indexNe++ ;
 498                     p1 = kTRUE;
 499                   }
 500                   if(m != 0){
 501                     if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
 502                       dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 503                         Fill(photon1->Pt());
 504                     new((*plNePHOS)[indexNePHOS])       TParticle(*photon1) ;
 505                     (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))->SetStatusCode(0) ;
 506                     indexNePHOS++;
 507                     p1 = kTRUE;
 508                   }//Is in PHOS
 509                 }
 510               }
 511               if(pGamma2.Pt() > 0. && TMath::Abs(pGamma2.Eta())<fEtaCut){
 512
 513                 MakePhoton(pGamma2);
 514                 if(pGamma2.Pt() > fNeutralPtCut){
 515
 516                   TParticle * photon2 =
 517                     new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(),
 518                                   pGamma2.Pz(),pGamma2.E(),0,0,0,0);
 519                   geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x);
 520                   if(photon2->Phi()>fPhiEMCALCut[0] &&
 521                      photon2->Phi()<fPhiEMCALCut[1] && m == 0){
 522                     if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
 523                       dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 524                         Fill(photon2->Pt());
 525                     new((*particleList)[index]) TParticle(*photon2) ;
 526                     (dynamic_cast<TParticle*>(particleList->At(index)))->SetStatusCode(0) ;
 527                     index++ ;
 528
 529                     new((*plNe)[indexNe])       TParticle(*photon2) ;
 530                     (dynamic_cast<TParticle*>(plNe->At(indexNe)))      ->SetStatusCode(0) ;
 531                     indexNe++ ;
 532                   }
 533                   if(m != 0){
 534                     if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
 535                       dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 536                         Fill(photon2->Pt());
 537                     new((*plNePHOS)[indexNePHOS])       TParticle(*photon2) ;
 538                     (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS))) ->SetStatusCode(0) ;
 539                     indexNePHOS++;
 540                   }
 541                   if(p1){
 542                     //              e = (pGamma1+pGamma2).E();
 543                     //              if(IsAngleInWindow(angle,e))
 544                       dynamic_cast<TH2F*>
 545                         (fListHistos->FindObject("AnglePairAccepted"))->
 546                         Fill(pPi0.E(),angle);
 547                   }
 548                 }
 549               }//photon2 in acceptance
 550             }//if angle > mindist
 551           }//if pi0
 552       }
 553     }//for (iParticle<nParticle)
 554   }
 555
 556   //Info("AddHIJINGToList","End HIJING");
 557 }
 558
 559 //____________________________________________________________________________
 560 Double_t AliPHOSGammaJet::CalculateJetRatioLimit(const Double_t ptg,
 561                                                  const Double_t *par,
 562                                                  const Double_t *x) {
 563
 564   //Info("CalculateLimit","x1 %f, x2%f",x[0],x[1]);
 565   Double_t epp = par[0] + par[1] * ptg ;
 566   Double_t spp = par[2] + par[3] * ptg ;
 567   Double_t f   = x[0]   + x[1]   * ptg ;
 568   Double_t epb = epp + par[4] ;
 569   Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
 570   Double_t rat = (epb - spb * f) / ptg ;
 571   //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
 572   //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
 573   return rat ;
 574 }
 575
 576 //____________________________________________________________________________
 577 void AliPHOSGammaJet::CreateParticleList(Int_t iEvent,
 578                                          TClonesArray * particleList,
 579                                          TClonesArray * plCh,
 580                                          TClonesArray * plNe,
 581                                          TClonesArray * plNePHOS)
 582 {
 583   //Info("CreateParticleList","Inside");
 584   AliPHOSGetter * gime = AliPHOSGetter::Instance(fInputFileName) ;
 585   const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
 586   gime->Event(iEvent, "X") ;
 587
 588
 589   Int_t index = particleList->GetEntries() ;
 590   Int_t indexCh     = plCh->GetEntries() ;
 591   Int_t indexNe     = plNe->GetEntries() ;
 592   Int_t indexNePHOS = plNePHOS->GetEntries() ;
 593   Int_t iParticle = 0 ;
 594   Double_t charge = 0.;
 595   Int_t m = 0;
 596   Double_t x = 0., z = 0.;
 597   if(!fOptFast){
 598
 599     for (iParticle=0 ; iParticle <  gime->NPrimaries() ; iParticle++) {
 600       const TParticle * particle = gime->Primary(iParticle) ;
 601
 602       m = 0 ;
 603       x = 0. ;
 604       z = 0. ;
 605
 606       //Keep Stable particles within eta range
 607       if((particle->GetStatusCode() == 1) &&
 608          (particle->Pt() > 0)){
 609         if(TMath::Abs(particle->Eta())<fEtaCut){
 610           //Fill lists
 611
 612           charge = TDatabasePDG::Instance()
 613             ->GetParticle(particle->GetPdgCode())->Charge();
 614           if((charge != 0) && (particle->Pt() > fChargedPtCut)){
 615
 616             if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
 617               dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 618                 Fill(particle->Pt());
 619             new((*plCh)[indexCh++])       TParticle(*particle) ;
 620             new((*particleList)[index++]) TParticle(*particle) ;
 621           }
 622           else if((charge == 0) && (particle->Pt() > fNeutralPtCut)){
 623             geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
 624             if(m != 0)
 625               {//Is in PHOS
 626                 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
 627                   dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 628                     Fill(particle->Pt());
 629
 630                 new((*plNePHOS)[indexNePHOS++])       TParticle(*particle) ;
 631               }
 632             if((particle->Phi()>fPhiEMCALCut[0]) &&
 633                (particle->Phi()<fPhiEMCALCut[1]) && m == 0)
 634               {//Is in EMCAL
 635                 if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
 636                   dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 637                     Fill(particle->Pt());
 638                 new((*plNe)[indexNe++])       TParticle(*particle) ;
 639                 new((*particleList)[index++]) TParticle(*particle) ;
 640               }
 641           }
 642         }
 643       }//final particle, etacut
 644     }//for (iParticle<nParticle)
 645   }// No Op
 646   else
 647     {
 648       Double_t mass = TDatabasePDG::Instance()->GetParticle(111)->Mass();
 649       TLorentzVector pPi0, pGamma1, pGamma2 ;
 650       Double_t angle = 0, cellDistance = 0.;
 651
 652       fFastRec = new AliPHOSFastGlobalReconstruction(fInputFileName);
 653       fFastRec->FastReconstruction(iEvent);
 654
 655       for (iParticle=0 ; iParticle <  gime->NPrimaries() ; iParticle++) {
 656         const TParticle * particle = gime->Primary(iParticle) ;
 657         m = 0 ;
 658         x = 0. ;
 659         z = 0. ;
 660         //Keep Stable particles within eta range
 661         if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){
 662
 663           //Fill lists
 664
 665           charge = TDatabasePDG::Instance()
 666             ->GetParticle(particle->GetPdgCode())->Charge();
 667           if((charge != 0) && (particle->Pt() > fChargedPtCut) && (TMath::Abs(particle->Eta())<fEtaCut)){
 668             new((*plCh)[indexCh++])       TParticle(*particle) ;
 669             new((*particleList)[index++]) TParticle(*particle) ;
 670           }
 671           else if(charge == 0) {
 672             geom->ImpactOnEmc(particle->Theta(),particle->Phi(), m,z,x);
 673             if((particle->GetPdgCode() != 111) && particle->Pt() > 0 &&
 674                (TMath::Abs(particle->Eta())<fEtaCut))
 675             {
 676
 677               TLorentzVector part(particle->Px(),particle->Py(),
 678                                   particle->Pz(),particle->Energy());
 679
 680               MakePhoton(part);
 681
 682               if(part.Pt() > fNeutralPtCut){
 683                 if(particle->Phi()>fPhiEMCALCut[0] &&
 684                    particle->Phi()<fPhiEMCALCut[1] && m == 0)
 685                   {
 686                     new((*particleList)[index]) TParticle(*particle) ;
 687                     (dynamic_cast<TParticle*>(particleList->At(index)))
 688                       ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 689                     index++ ;
 690
 691                     new((*plNe)[indexNe])       TParticle(*particle) ;
 692                     (dynamic_cast<TParticle*>(plNe->At(indexNe)))
 693                       ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 694                     indexNe++ ;
 695                   }
 696                 if(m != 0)
 697                   {
 698                     new((*plNePHOS)[indexNePHOS])       TParticle(*particle) ;
 699                     (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
 700                       ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 701                     indexNePHOS++ ;
 702                   }
 703               }// Small pt
 704             } //No Pi0
 705             if((particle->GetPdgCode() == 111) && (particle->Pt() > fNeutralPtCut) &&
 706                (TMath::Abs(particle->Eta())<fEtaCut+1))
 707               {
 708
 709                 pPi0.SetPxPyPzE(particle->Px(),particle->Py(),particle->Pz(),
 710                                 particle->Energy());
 711
 712                 //Decay
 713
 714                 Pi0Decay(mass,pPi0,pGamma1,pGamma2,angle);
 715                 //Check if decay photons are too close for PHOS
 716                 cellDistance = angle*460; //cm
 717
 718                 if (cellDistance < fMinDistance) {
 719
 720                   //Pi0 inside phi EMCAL acceptance
 721
 722
 723                     TLorentzVector part(particle->Px(),particle->Py(),
 724                                         particle->Pz(),particle->Energy());
 725                     MakePhoton(part);
 726
 727                     if(part.Pt() > fNeutralPtCut){
 728                       if(particle->Phi()>fPhiEMCALCut[0] &&
 729                          particle->Phi()<fPhiEMCALCut[1] && m == 0){
 730                         if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
 731                           dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 732                             Fill(particle->Pt());
 733
 734                         new((*plNe)[indexNe])       TParticle(*particle) ;
 735                         (dynamic_cast<TParticle*>(plNe->At(indexNe)))
 736                           ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 737                         new((*particleList)[index]) TParticle(*particle) ;
 738                         (dynamic_cast<TParticle*>(particleList->At(index)))
 739                           ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 740                         index++;
 741                         indexNe++;
 742                       }//InEMCAL
 743                       if(m != 0){
 744                         if(strstr(fOptionGJ,"deb all")|| strstr(fOptionGJ,"deb"))
 745                           dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 746                             Fill(particle->Pt());
 747                         new((*plNePHOS)[indexNePHOS])       TParticle(*particle) ;
 748                         (dynamic_cast<TParticle*>(plNePHOS->At(indexNePHOS)))
 749                           ->SetMomentum(part.Px(),part.Py(),part.Pz(),part.E());
 750                         indexNePHOS++;
 751                       }//In PHOS
 752                     }//Small Pt
 753                 }// if cell<distance
 754                 else {
 755
 756                   dynamic_cast<TH2F*>(fListHistos->FindObject("AnglePair"))
 757                     ->Fill(pPi0.E(),angle);
 758
 759                   Bool_t p1 = kFALSE;
 760
 761                   if(pGamma1.Pt() > 0 && TMath::Abs(pGamma1.Eta())<fEtaCut){
 762                     MakePhoton(pGamma1);
 763
 764                     if(pGamma1.Pt() > fNeutralPtCut){
 765                       TParticle * photon1 =
 766                         new TParticle(22,1,0,0,0,0,pGamma1.Px(),pGamma1.Py(),
 767                                       pGamma1.Pz(),pGamma1.E(),0,0,0,0);
 768                       geom->ImpactOnEmc(photon1->Theta(),photon1->Phi(), m,z,x);
 769                       if( photon1->Phi()>fPhiEMCALCut[0] && photon1->Phi()<fPhiEMCALCut[1]
 770                           && m == 0){
 771                       if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
 772                         dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 773                           Fill(photon1->Pt());
 774                       new((*plNe)[indexNe++])       TParticle(*photon1) ;
 775                       new((*particleList)[index++]) TParticle(*photon1) ;
 776                       //photon1->Print();
 777                       p1 = kTRUE;
 778                       }
 779                       if(m != 0){
 780                         if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
 781                           dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 782                             Fill(photon1->Pt());
 783                         new((*plNePHOS)[indexNePHOS++])       TParticle(*photon1) ;
 784                       p1 = kTRUE;
 785                       }
 786                     }
 787                   }
 788                   if(pGamma2.Pt() > 0 && TMath::Abs(pGamma2.Eta())<fEtaCut){
 789                     MakePhoton(pGamma2);
 790
 791                     if(pGamma2.Pt() > fNeutralPtCut ){
 792                       TParticle * photon2 =
 793                         new TParticle(22,1,0,0,0,0,pGamma2.Px(), pGamma2.Py(),
 794                                       pGamma2.Pz(),pGamma2.E(),0,0,0,0);
 795                       geom->ImpactOnEmc(photon2->Theta(),photon2->Phi(), m,z,x);
 796                       if(photon2->Phi()>fPhiEMCALCut[0] &&
 797                          photon2->Phi()<fPhiEMCALCut[1] && m == 0){
 798                         if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
 799                           dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 800                             Fill(photon2->Pt());
 801                         new((*plNe)[indexNe++])       TParticle(*photon2) ;
 802                         new((*particleList)[index++]) TParticle(*photon2) ;
 803                       }
 804                       if(m != 0){
 805                         if(strstr(fOptionGJ,"deb all") || strstr(fOptionGJ,"deb"))
 806                           dynamic_cast<TH1F*>(fListHistos->FindObject("PtSpectra"))->
 807                             Fill(photon2->Pt());
 808                         new((*plNePHOS)[indexNePHOS++])       TParticle(*photon2) ;
 809                       }
 810
 811                       if(p1){
 812                         // Float_t e = (pGamma1+pGamma2).E();
 813                         // if(IsAngleInWindow(angle,e))
 814                           dynamic_cast<TH2F*>
 815                             (fListHistos->FindObject("AnglePairAccepted"))->
 816                             Fill(pPi0.E(),angle);
 817                       }
 818                     }
 819                   }//photon2 in acceptance
 820                 }//if angle > mindist
 821               }//if pi0
 822           }//If neutral
 823         }//final particle, etacut
 824       }//for (iParticle<nParticle)
 825     }//OptFast
 826   //gime->Delete() ;
 827 }
 828 //____________________________________________________________________________
 829 void AliPHOSGammaJet::CreateParticleListFromESD(TClonesArray * pl,
 830                                                 TClonesArray * plCh,
 831                                                 TClonesArray * plNe,
 832                                                 TClonesArray * plNePHOS,
 833                                                 const AliESD * esd){
 834
 835   //Create a list of particles from the ESD. These particles have been measured
 836   //by the Central Tracking system (TPC), PHOS and EMCAL
 837   //(EMCAL not available for the moment).
 838   //Info("CreateParticleListFromESD","Inside");
 839
 840   Int_t index = pl->GetEntries() ;
 841   Int_t npar  = 0 ;
 842   Double_t pid[AliPID::kSPECIESN];
 843
 844   //########### PHOS ##############
 845   //Info("CreateParticleListFromESD","Fill ESD PHOS list");
 846   Int_t begphos = esd->GetFirstPHOSParticle();
 847   Int_t endphos = esd->GetFirstPHOSParticle() +
 848     esd->GetNumberOfPHOSParticles() ;
 849   Int_t indexNePHOS = plNePHOS->GetEntries() ;
 850   if(strstr(fOptionGJ,"deb all"))
 851     Info("CreateParticleListFromESD","PHOS: first particle %d, last particle %d",
 852          begphos,endphos);
 853
 854   for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
 855     AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
 856
 857     //Create a TParticle to fill the particle list
 858
 859     Double_t en = track->GetPHOSsignal() ;
 860     Double_t * p = new Double_t();
 861     track->GetPHOSposition(p) ;
 862     TVector3 pos(p[0],p[1],p[2]) ;
 863     Double_t phi  = pos.Phi();
 864     Double_t theta= pos.Theta();
 865     Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
 866     Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
 867     Double_t pz = en*TMath::Cos(theta);
 868
 869     TParticle * particle = new TParticle() ;
 870     particle->SetMomentum(px,py,pz,en) ;
 871
 872     //Select only photons
 873
 874     track->GetPHOSpid(pid);
 875     //cout<<"pid "<<pid[AliPID::kPhoton]<<endl ;
 876     if( pid[AliPID::kPhoton] > 0.75)
 877       new((*plNePHOS)[indexNePHOS++])   TParticle(*particle) ;
 878   }
 879
 880   //########### TPC #####################
 881   //Info("CreateParticleListFromESD","Fill ESD TPC list");
 882   Int_t begtpc = 0 ;
 883   Int_t endtpc = esd->GetNumberOfTracks() ;
 884   Int_t indexCh     = plCh->GetEntries() ;
 885   if(strstr(fOptionGJ,"deb all"))
 886     Info("CreateParticleListFromESD","TPC: first particle %d, last particle %d",
 887          begtpc,endtpc);
 888
 889   for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
 890     AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
 891
 892     Double_t en = track ->GetTPCsignal() ;
 893     TVector3 mom = track->P3() ;
 894     Double_t px = mom.Px();
 895     Double_t py = mom.Py();
 896     Double_t pz = mom.Pz(); //Check with TPC people if this is correct.
 897
 898     //cout<<"TPC signal "<<en<<endl;
 899     //cout<<"px "<<px<<"; py "<<py<<"; pz "<<pz<<endl;
 900     TParticle * particle = new TParticle() ;
 901     particle->SetMomentum(px,py,pz,en) ;
 902
 903     new((*plCh)[indexCh++])       TParticle(*particle) ;
 904     new((*pl)[index++])           TParticle(*particle) ;
 905
 906   }
 907
 908   //################ EMCAL ##############
 909   Double_t v[3] ; //vertex ;
 910   esd->GetVertex()->GetXYZ(v) ;
 911   //##########Uncomment when ESD for EMCAL works ##########
 912   //Info("CreateParticleListFromESD","Fill ESD EMCAL list");
 913
 914   Int_t begem = esd->GetFirstEMCALParticle();
 915   Int_t endem = esd->GetFirstEMCALParticle() +
 916     esd->GetNumberOfEMCALParticles() ;
 917   Int_t indexNe  = plNe->GetEntries() ;
 918   if(strstr(fOptionGJ,"deb all"))
 919     Info("CreateParticleListFromESD","EMCAL: first particle %d, last particle %d",
 920          begem,endem);
 921
 922   for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
 923      AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
 924
 925     Double_t en = track->GetEMCALsignal() ;
 926     Double_t *p = new Double_t();
 927     track->GetEMCALposition(p) ;
 928     TVector3 pos(p[0],p[1],p[2]) ;
 929     Double_t phi  = pos.Phi();
 930     Double_t theta= pos.Theta();
 931     Double_t px = en*TMath::Cos(phi)*TMath::Sin(theta);;
 932     Double_t py = en*TMath::Sin(phi)*TMath::Sin(theta);
 933     Double_t pz = en*TMath::Cos(theta);
 934     //cout<<"EMCAL signal "<<en<<endl;
 935     //cout<<"px "<<px<<"; py "<<py<<"; pz "<<pz<<endl;
 936     //TParticle * particle = new TParticle() ;
 937     //particle->SetMomentum(px,py,pz,en) ;
 938
 939     Int_t pdg = 0;
 940     //  //Uncomment if PID IS WORKING, photon and pi0 idenitification.
 941     //  //if( pid[AliPID::kPhoton] > 0.75) //This has to be fixen.
 942     //  //pdg = 22;
 943     //  //else if( pid[AliPID::kPi0] > 0.75)
 944     //  //pdg = 111;
 945     pdg = 22; //No PID, assume all photons
 946     TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
 947                                          px, py, pz, en, v[0], v[1], v[2], 0);
 948
 949     new((*plNe)[indexNe++])       TParticle(*particle) ;
 950     new((*pl)[index++])           TParticle(*particle) ;
 951   }
 952
 953   //  Info("CreateParticleListFromESD","End Inside");
 954
 955 }
 956
 957
 958
 959 //____________________________________________________________________________
 960 void AliPHOSGammaJet::Exec(Option_t *option)
 961 {
 962   // does the job
 963   fOptionGJ = option;
 964   MakeHistos() ;
 965
 966   AliESD * esd = 0;
 967   TChain * t = 0 ;
 968
 969   if(fESDdata){
 970     // Create chain of esd trees
 971     const UInt_t kNevent = static_cast<UInt_t>(GetNEvent()) ;
 972     t = ReadESD(kNevent, fDirName, fESDTree, fPattern) ;
 973     if(!t) {
 974       AliError("Could not create the TChain") ;
 975       //break ;
 976     }
 977
 978     // ESD
 979     t->SetBranchAddress("ESD",&esd); // point to the container esd where to put the event from the esdTree
 980
 981   }
 982
 983
 984 //   AliGenPythia* pyth = (AliGenPythia*) gAlice->Generator();
 985 //   pyth->Init();
 986
 987   TClonesArray * particleList = new TClonesArray("TParticle",1000);
 988   TClonesArray * plCh         = new TClonesArray("TParticle",1000);
 989   TClonesArray * plNe         = new TClonesArray("TParticle",1000);
 990   TClonesArray * plNePHOS     = new TClonesArray("TParticle",1000);
 991
 992   for (Int_t iEvent = 0 ; iEvent < fNEvent ; iEvent++) {
 993     if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
 994       Info("Exec", "Event %d", iEvent) ;
 995
 996     fRan.SetSeed(0);
 997
 998     Double_t phig = 0., phil = 0., phich = 0 , phipi = 0;
 999     Double_t etag = 0., etal = 0., etach = 0., etapi = 0. ;
1000     Double_t ptg  = 0., ptl  = 0., ptch  = 0., ptpi  = 0.;
1001
1002     TLorentzVector jet   (0,0,0,0);
1003     TLorentzVector jettpc(0,0,0,0);
1004
1005     if(fESDdata){
1006
1007       Int_t iNbytes = t->GetEntry(iEvent); // store event in esd
1008       //cout<<"nbytes "<<iNbytes<<endl;
1009       if ( iNbytes == 0 ) {
1010         AliError("Empty TChain") ;
1011         break ;
1012       }
1013       CreateParticleListFromESD(particleList, plCh,plNe,plNePHOS, esd); //,iEvent);
1014     }
1015     else{
1016       CreateParticleList(iEvent, particleList, plCh,plNe,plNePHOS);
1017
1018       //    TLorentzVector pyjet(0,0,0,0);
1019
1020       //     Int_t nJ, nJT;
1021       //     Float_t jets[4][10];
1022       //     pyth->SetJetReconstructionMode(1);
1023       //     pyth->LoadEvent();
1024       //     pyth->GetJets(nJ, nJT, jets);
1025
1026       //     Float_t pxJ = jets[0][0];
1027       //     Float_t pyJ = jets[1][0];
1028       //     Float_t pzJ = jets[2][0];
1029       //     Float_t eJ  = jets[3][0];
1030       //     pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
1031
1032       //     if(nJT > 1){
1033       //       //Info("Exec",">>>>>>>>>>Number of jets !!!!   %d",nJT);
1034       //       for (Int_t iJ = 1; iJ < nJT; iJ++) {
1035       //        Float_t pxJ = jets[0][iJ];
1036       //        Float_t pyJ = jets[1][iJ];
1037       //        Float_t pzJ = jets[2][iJ];
1038       //        Float_t eJ  = jets[3][iJ];
1039       //        pyjet.SetPxPyPzE(pxJ,pyJ,pzJ,eJ ) ;
1040       //        //Info("Exec",">>>>>Pythia Jet: %d, Phi %f, Eta %f, Pt %f",
1041       //        //           iJ,pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1042       //       }
1043
1044       //     }
1045
1046       if(fHIJING)
1047         AddHIJINGToList(iEvent, particleList, plCh,plNe, plNePHOS);
1048     }
1049
1050     Bool_t iIsInPHOS = kFALSE ;
1051     GetGammaJet(plNePHOS, ptg, phig, etag, iIsInPHOS) ;
1052
1053     if(iIsInPHOS){
1054
1055       //Info("Exec"," In PHOS") ;
1056       dynamic_cast<TH1F*>(fListHistos->FindObject("NGamma"))->Fill(ptg);
1057       dynamic_cast<TH2F*>(fListHistos->FindObject("PhiGamma"))
1058         ->Fill(ptg,phig);
1059       dynamic_cast<TH2F*>(fListHistos->FindObject("EtaGamma"))
1060         ->Fill(ptg,etag);
1061       if(strstr(fOptionGJ,"deb")||strstr(fOptionGJ,"deb all"))
1062         Info("Exec", "Gamma: pt %f, phi %f, eta %f", ptg,
1063              phig,etag) ;
1064
1065 //       cout<<"n charged "<<plCh->GetEntries()<<endl;
1066 //       cout<<"n neutral "<<plNe->GetEntries()<<endl;
1067 //       cout<<"n All     "<<particleList->GetEntries()<<endl;
1068
1069       GetLeadingCharge(plCh, ptg, phig, ptch, etach, phich) ;
1070       GetLeadingPi0   (plNe, ptg, phig, ptpi, etapi, phipi) ;
1071
1072 //       cout<<"n2 charged "<<plCh->GetEntries()<<endl;
1073 //       cout<<"n2 neutral "<<plNe->GetEntries()<<endl;
1074 //       cout<<"n2 All     "<<particleList->GetEntries()<<endl;
1075
1076
1077       //TPC+EMCAL
1078
1079       //Is the leading cone inside EMCAL?
1080       Bool_t insidech = kFALSE ;
1081       if((phich - fCone) >  fPhiEMCALCut[0] &&
1082          (phich + fCone) <  fPhiEMCALCut[1]){
1083         insidech = kTRUE ;
1084       }
1085       Bool_t insidepi = kFALSE ;
1086       if((phipi - fCone) >  fPhiEMCALCut[0] &&
1087          (phipi + fCone) <  fPhiEMCALCut[1]){
1088         insidepi = kTRUE ;
1089       }
1090
1091       if ((ptch > 0 || ptpi > 0)){
1092         if((ptch > ptpi) && insidech){
1093           phil = phich ;
1094           etal = etach ;
1095           ptl  = ptch ;
1096           dynamic_cast<TH2F*>(fListHistos->FindObject("ChargeRatio"))
1097             ->Fill(ptg,ptch/ptg);
1098           dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiCharge"))
1099             ->Fill(ptg,phig-phich);
1100           dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaCharge"))
1101             ->Fill(ptg,etag-etach);
1102           if(strstr(fOptionGJ,"deb"))
1103           Info("Exec"," Charged Leading") ;
1104         }
1105         if((ptpi > ptch) && insidepi){
1106           phil = phipi ;
1107           etal = etapi ;
1108           ptl  = ptpi ;
1109
1110           dynamic_cast<TH2F*>(fListHistos->FindObject("Pi0Ratio"))
1111             ->Fill(ptg,ptpi/ptg);
1112           dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiPi0"))
1113             ->Fill(ptg,phig-phipi);
1114           dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaPi0"))
1115             ->Fill(ptg,etag-etapi);
1116
1117           if(ptpi > 0. && strstr(fOptionGJ,"deb"))
1118             Info("Exec"," Pi0 Leading") ;
1119         }
1120
1121         if(strstr(fOptionGJ,"deb"))
1122           Info("Exec","Leading pt %f, phi %f",ptl,phil);
1123         if(insidech || insidepi){
1124           if(!fAnyConeOrPt){
1125
1126             MakeJet(particleList, ptg, phig, ptl, phil, etal, "", jet);
1127
1128             if(strstr(fOptionGJ,"deb")){
1129 //            Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
1130 //                 pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1131               Info("Exec","TPC+EMCAL Jet: Phi %f, Eta %f, Pt %f",
1132                    jet.Phi(),jet.Eta(),jet.Pt());
1133             }
1134 //          dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiJet"))
1135 //            ->Fill(ptg,pyjet.Phi()-jet.Phi());
1136 //          dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaJet"))
1137 //            ->Fill(ptg,pyjet.Eta()-jet.Eta());
1138 //          dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtJet"))
1139 //            ->Fill(ptg,pyjet.Pt()-jet.Pt());
1140           }
1141           else
1142             MakeJetAnyConeOrPt(particleList, ptg, phig, ptl, phil, etal, "");
1143         }
1144
1145         //TPC
1146         if(fOnlyCharged && ptch > 0.)
1147           {
1148             if(strstr(fOptionGJ,"deb"))
1149               Info("Exec","Leading TPC pt %f, phi %f",ptch,phich);
1150
1151             dynamic_cast<TH2F*>(fListHistos->FindObject("TPCRatio"))
1152               ->Fill(ptg,ptch/ptg);
1153             dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPC"))
1154               ->Fill(ptg,phig-phich);
1155             dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPC"))
1156               ->Fill(ptg,etag-etach);
1157
1158             if(!fAnyConeOrPt){
1159
1160               MakeJet(plCh, ptg, phig, ptch, phich, etach, "TPC",jettpc);
1161
1162               if(strstr(fOptionGJ,"deb")){
1163 //              Info("Exec","Pythia Jet: Phi %f, Eta %f, Pt %f",
1164 //                   pyjet.Phi(),pyjet.Eta(),pyjet.Pt());
1165                 Info("Exec","TPC Jet: Phi %f, Eta %f, Pt %f",
1166                      jettpc.Phi(),jettpc.Eta(),jettpc.Pt());
1167               }
1168 //            dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPhiTPCJet"))
1169 //              ->Fill(ptg,pyjet.Phi()-jettpc.Phi());
1170 //            dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaEtaTPCJet"))
1171 //              ->Fill(ptg,pyjet.Eta()-jettpc.Eta());
1172 //            dynamic_cast<TH2F*>(fListHistos->FindObject("DeltaPtTPCJet"))
1173 //              ->Fill(ptg,pyjet.Pt()-jettpc.Pt());
1174             }
1175             else
1176               MakeJetAnyConeOrPt(plCh, ptg, phig, ptch, phich, etach, "TPC");
1177
1178           }
1179       }
1180     }
1181
1182     particleList->Delete() ;
1183     plCh->Delete() ;
1184     plNe->Delete() ;
1185     plNePHOS->Delete() ;
1186   }//loop: events
1187
1188   delete plNe ;
1189   delete plCh ;
1190   delete particleList ;
1191
1192   fOutputFile->Write() ;
1193   fOutputFile->cd();
1194   this->Write();
1195 }
1196
1197 //____________________________________________________________________________
1198 void AliPHOSGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg,
1199                                     TString conf, TString type)
1200 {
1201   //Fill jet fragmentation histograms if !fAnyCone,
1202   //only for fCone and fPtThres
1203   TParticle * particle = 0 ;
1204   Int_t ipr = -1 ;
1205   Float_t  charge = 0;
1206
1207   TIter next(pl) ;
1208   while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1209     ipr++ ;
1210     Double_t pt = particle->Pt();
1211
1212     charge = TDatabasePDG::Instance()
1213       ->GetParticle(particle->GetPdgCode())->Charge();
1214     if(charge != 0){//Only jet Charged particles
1215       dynamic_cast<TH2F*>
1216         (fListHistos->FindObject(type+conf+"Fragment"))
1217         ->Fill(ptg,pt/ptg);
1218       dynamic_cast<TH2F*>
1219         (fListHistos->FindObject(type+conf+"PtDist"))
1220         ->Fill(ptg,pt);
1221     }
1222   }
1223   if(type == "Bkg"){
1224     dynamic_cast<TH1F*>
1225       (fListHistos->FindObject("NBkg"+conf))->Fill(ipr);
1226   }
1227 }
1228 //____________________________________________________________________________
1229 void AliPHOSGammaJet::FillJetHistosAnyConeOrPt(TClonesArray * pl, Double_t ptg,
1230                                                TString conf, TString type,
1231                                                TString cone, TString ptcut)
1232 {
1233    //Fill jet fragmentation histograms if fAnyCone,
1234    //for several cones and pt thresholds
1235    TParticle *particle = 0;
1236    Int_t ipr=-1;
1237    Float_t  charge = 0;
1238
1239    TIter next(pl) ;
1240    while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1241      ipr++;
1242      Double_t pt = particle->Pt();
1243      charge = TDatabasePDG::Instance()
1244        ->GetParticle(particle->GetPdgCode())->Charge();
1245      if(charge != 0){//Only jet Charged particles
1246        dynamic_cast<TH2F*>
1247          (fListHistos->FindObject(type+conf+"FragmentCone"+cone+"Pt"+ptcut))
1248          ->Fill(ptg,pt/ptg);
1249        dynamic_cast<TH2F*>
1250          (fListHistos->FindObject(type+conf+"PtDistCone"+cone+"Pt"+ptcut))
1251          ->Fill(ptg,pt);
1252      }
1253    }//while
1254
1255    if(type == "Bkg"){
1256      dynamic_cast<TH1F*>
1257        (fListHistos->FindObject("NBkg"+conf+"Cone"+cone+"Pt"+ptcut))
1258        ->Fill(ipr);
1259    }
1260 }
1261
1262 //____________________________________________________________________________
1263 void AliPHOSGammaJet::GetGammaJet(TClonesArray * pl, Double_t &pt,
1264                                   Double_t &phi, Double_t &eta, Bool_t &Is) const
1265 {
1266   //Search for the prompt photon in PHOS with pt > fPtCut
1267   pt  = -10.;
1268   eta = -10.;
1269   phi = -10.;
1270
1271   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
1272     TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
1273
1274     if((particle->Pt() > fPtCut) && (particle->Pt() > pt)){
1275
1276       pt  = particle->Pt();
1277       phi = particle->Phi() ;
1278       eta = particle->Eta() ;
1279       Is  = kTRUE;
1280     }
1281   }
1282 }
1283
1284 //____________________________________________________________________________
1285 void  AliPHOSGammaJet::GetLeadingCharge(TClonesArray * pl,
1286                                         Double_t ptg, Double_t phig,
1287                                         Double_t &pt, Double_t &eta, Double_t &phi) const
1288 {
1289   //Search for the charged particle with highest with
1290   //Phi=Phi_gamma-Pi and pT=0.1E_gamma
1291   pt  = -100.;
1292   eta = -100;
1293   phi = -100;
1294
1295   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
1296
1297     TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
1298
1299     Double_t ptl  = particle->Pt();
1300     Double_t rat  = ptl/ptg ;
1301     Double_t phil = particle->Phi() ;
1302
1303      if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
1304         (rat > fRatioMinCut) && (rat < fRatioMaxCut)  && (ptl  > pt)) {
1305       eta = particle->Eta() ;
1306       phi = phil ;
1307       pt  = ptl ;
1308       //printf("GetLeadingCharge: %f %f %f %f \n", pt, eta, phi,rat) ;
1309     }
1310   }
1311   //printf("GetLeadingCharge: %f %f %f \n", pt, eta, phi) ;
1312
1313 }
1314
1315
1316 //____________________________________________________________________________
1317 void  AliPHOSGammaJet::GetLeadingPi0(TClonesArray * pl,
1318                                      Double_t ptg, Double_t phig,
1319                                      Double_t &pt,  Double_t &eta, Double_t &phi)
1320 {
1321
1322   //Search for the neutral pion with highest with
1323   //Phi=Phi_gamma-Pi and pT=0.1E_gamma
1324   pt  = -100.;
1325   eta = -100.;
1326   phi = -100.;
1327   Double_t ptl = -100.;
1328   Double_t rat = -100.;
1329   Double_t phil = -100. ;
1330
1331   TIter next(pl);
1332   TParticle * particle = 0;
1333   Float_t ef = 0;
1334   if(!fOptFast){
1335     Float_t e = 0;
1336     while ( (particle = (TParticle*)next()) ) {
1337       if( particle->GetPdgCode() == 111){
1338         ptl  = particle->Pt();
1339         rat  = ptl/ptg ;
1340         phil = particle->Phi() ;
1341         e    = particle->Energy();
1342         dynamic_cast<TH2F*>
1343           (fListHistos->FindObject("AnglePairNoCut"))->
1344           Fill(e,0.1);
1345         dynamic_cast<TH2F*>
1346           (fListHistos->FindObject("InvMassPairNoCut"))->
1347           Fill(ptg,1.35);
1348
1349         if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
1350            (rat > fRatioMinCut) && (rat < fRatioMaxCut)) {
1351
1352           dynamic_cast<TH2F*>
1353             (fListHistos->FindObject("AnglePairLeadingCut"))->
1354             Fill(e,0.1);
1355           dynamic_cast<TH2F*>
1356             (fListHistos->FindObject("InvMassPairLeadingCut"))->
1357             Fill(ptg,1.35);
1358
1359           dynamic_cast<TH2F*>
1360             (fListHistos->FindObject("AnglePairAngleCut"))->
1361             Fill(e,0.15);
1362           dynamic_cast<TH2F*>
1363             (fListHistos->FindObject("InvMassPairAngleCut"))->
1364             Fill(ptg,1.36);
1365
1366           dynamic_cast<TH2F*>
1367             (fListHistos->FindObject("InvMassPairAllCut"))->
1368             Fill(ptg,0.27);
1369           dynamic_cast<TH2F*>
1370             (fListHistos->FindObject("AnglePairAllCut"))->
1371             Fill(e,1.34);
1372
1373
1374           if(ptl  > pt){
1375             eta = particle->Eta() ;
1376             phi = phil ;
1377             pt  = ptl ;
1378             ef  = e;
1379             //printf("GetLeadingPi0: %f %f %f %f %f \n", pt, eta, phi, rat, ptg) ;
1380           }
1381         }
1382       }
1383
1384       dynamic_cast<TH2F*>
1385         (fListHistos->FindObject("InvMassPairLeading"))->
1386         Fill(ptg,1.35);
1387       dynamic_cast<TH2F*>
1388         (fListHistos->FindObject("AnglePairLeading"))->
1389         Fill(ef,0.1);
1390     }
1391   }//No fOptfast
1392   else{
1393     Int_t iPrimary = -1;
1394     TLorentzVector gammai,gammaj;
1395     Double_t angle = 0., e = 0., invmass = 0.;
1396     Double_t anglef = 0., ef = 0., invmassf = 0.;
1397     Int_t ksPdg = 0;
1398     Int_t jPrimary=-1;
1399
1400     while ( (particle = (TParticle*)next()) ) {
1401       iPrimary++;
1402
1403       ksPdg = particle->GetPdgCode();
1404       ptl  = particle->Pt();
1405       if(ksPdg == 111){ //2 gamma
1406         rat = ptl/ptg ;
1407         phil = particle->Phi() ;
1408         if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
1409            ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
1410           eta = particle->Eta() ;
1411           phi = phil ;
1412           pt  = ptl ;
1413         }// cuts
1414       }// pdg = 111
1415       if(ksPdg == 22){//1 gamma
1416         particle->Momentum(gammai);
1417         jPrimary=-1;
1418         TIter next2(pl);
1419         while ( (particle = (TParticle*)next2()) ) {
1420           jPrimary++;
1421           if(jPrimary>iPrimary){
1422             ksPdg = particle->GetPdgCode();
1423             if(ksPdg == 22){
1424               particle->Momentum(gammaj);
1425               //Info("GetLeadingPi0","Egammai %f, Egammaj %f",
1426               //gammai.Pt(),gammaj.Pt());
1427
1428               ptl  = (gammai+gammaj).Pt();
1429               phil = (gammai+gammaj).Phi();
1430               if(phil < 0)
1431                 phil+=TMath::TwoPi();
1432               rat = ptl/ptg ;
1433               invmass = (gammai+gammaj).M();
1434               angle = gammaj.Angle(gammai.Vect());
1435               //Info("GetLeadingPi0","Angle %f", angle);
1436               e = (gammai+gammaj).E();
1437
1438               dynamic_cast<TH2F*>
1439                 (fListHistos->FindObject("AnglePairNoCut"))->
1440                 Fill(e,angle);
1441               dynamic_cast<TH2F*>
1442                 (fListHistos->FindObject("InvMassPairNoCut"))->
1443                 Fill(ptg,invmass);
1444
1445               if((rat > fRatioMinCut) && (rat < fRatioMaxCut) &&
1446                  ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
1447
1448                 dynamic_cast<TH2F*>
1449                   (fListHistos->FindObject("AnglePairLeadingCut"))->
1450                   Fill(e,angle);
1451                 dynamic_cast<TH2F*>
1452                   (fListHistos->FindObject("InvMassPairLeadingCut"))->
1453                   Fill(ptg,invmass);
1454
1455                 if(IsAngleInWindow(angle,e)){
1456                   dynamic_cast<TH2F*>
1457                     (fListHistos->FindObject("AnglePairAngleCut"))->
1458                     Fill(e,angle);
1459                   dynamic_cast<TH2F*>
1460                     (fListHistos->FindObject("InvMassPairAngleCut"))->
1461                     Fill(ptg,invmass);
1462
1463                   //Info("GetLeadingPi0","InvMass %f", invmass);
1464                   if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){
1465                     dynamic_cast<TH2F*>
1466                       (fListHistos->FindObject("InvMassPairAllCut"))->
1467                       Fill(ptg,invmass);
1468                     dynamic_cast<TH2F*>
1469                       (fListHistos->FindObject("AnglePairAllCut"))->
1470                       Fill(e,angle);
1471                     if(ptl > pt ){
1472                       pt       = ptl;
1473                       eta      = particle->Eta() ;
1474                       phi      = phil ;
1475                       ef       = e ;
1476                       anglef   = angle ;
1477                       invmassf = invmass ;
1478
1479                     }
1480                   }//cuts
1481                 }//(invmass>0.125) && (invmass<0.145)
1482               }//gammaj.Angle(gammai.Vect())<0.04
1483             }//if pdg = 22
1484           }//iprimary<jprimary
1485         }//while
1486       }// if pdg = 22
1487       //     }
1488     }//while
1489
1490     if(ef > 0.){
1491       dynamic_cast<TH2F*>
1492         (fListHistos->FindObject("InvMassPairLeading"))->
1493         Fill(ptg,invmassf);
1494       dynamic_cast<TH2F*>
1495         (fListHistos->FindObject("AnglePairLeading"))->
1496         Fill(ef,anglef);
1497     }
1498   }//fOptFast
1499   //  printf("GetLeadingPi0: %f %f %f \n", pt, eta, phi) ;
1500 }
1501
1502
1503 //____________________________________________________________________________
1504 void AliPHOSGammaJet::InitParameters()
1505 {
1506
1507   //Initialize the parameters of the analysis.
1508
1509   fAngleMaxParam.Set(4) ;
1510   fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
1511   fAngleMaxParam.AddAt(-0.25,1) ;
1512   fAngleMaxParam.AddAt(0.025,2) ;
1513   fAngleMaxParam.AddAt(-2e-4,3) ;
1514   fAnyConeOrPt    = kFALSE ;
1515   fOutputFileName = "GammaJet.root" ;
1516   fOptionGJ         = "";
1517   fHIJINGFileName = "galice.root" ;
1518   fHIJING         = kFALSE ;
1519   fMinDistance    = 3.6 ;
1520   fEtaCut         = 0.7 ;
1521   fInvMassMaxCut  = 0.15 ;
1522   fInvMassMinCut  = 0.12 ;
1523   fOnlyCharged    = kFALSE ;
1524   fOptFast        = kFALSE ;
1525   fESDdata        = kTRUE ;
1526   fPhiEMCALCut[0] = 60. *TMath::Pi()/180.;
1527   fPhiEMCALCut[1] = 180.*TMath::Pi()/180.;
1528   fPhiMaxCut      = 3.4 ;
1529   fPhiMinCut      = 2.9 ;
1530   fPtCut          = 10. ;
1531   fNeutralPtCut   = 0.5 ;
1532   fChargedPtCut   = 0.5 ;
1533   fTPCCutsLikeEMCAL  = kFALSE ;
1534   //Jet selection parameters
1535   //Fixed cut (old)
1536   fRatioMaxCut    = 1.0 ;
1537   fRatioMinCut    = 0.1 ;
1538   fJetRatioMaxCut = 1.2 ;
1539   fJetRatioMinCut = 0.8 ;
1540   fJetTPCRatioMaxCut = 1.2 ;
1541   fJetTPCRatioMinCut = 0.3 ;
1542   fSelect         = kFALSE  ;
1543
1544   fDirName      = "./" ;
1545   fESDTree      = "esdTree" ;
1546   fPattern      = "." ;
1547
1548   //Cut depending on gamma energy
1549
1550   fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applyed
1551   //Reconstructed jet energy dependence parameters
1552   //e_jet = a1+e_gamma b2.
1553   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1554   fJetE1[0] = -5.75; fJetE1[1] = -4.1;
1555   fJetE2[0] = 1.005; fJetE2[1] = 1.05;
1556
1557   //Reconstructed sigma of jet energy dependence parameters
1558   //s_jet = a1+e_gamma b2.
1559   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
1560   fJetSigma1[0] = 2.65;   fJetSigma1[1] = 2.75;
1561   fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
1562
1563   //Background mean energy and RMS
1564   //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1565   //Index 2-> (low pt jets)BKG > 0.5 GeV;
1566   //Index > 2, same for TPC conf
1567   fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
1568   fBkgMean[3] = 0.; fBkgMean[4] = 6.4;  fBkgMean[5] = 48.6;
1569   fBkgRMS[0]  = 0.; fBkgRMS[1]  = 7.5;  fBkgRMS[2]  = 22.0;
1570   fBkgRMS[3]  = 0.; fBkgRMS[4]  = 5.4;  fBkgRMS[5]  = 13.2;
1571
1572   //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
1573   //limits for monoenergetic jets.
1574   //Index 0-> No BKG; Index 1-> BKG > 2 GeV;
1575   //Index 2-> (low pt jets) BKG > 0.5 GeV;
1576   //Index > 2, same for TPC conf
1577
1578   fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ;
1579   fJetXMin1[3] =-2.0  ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1  ;
1580   fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034;
1581   fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
1582   fJetXMax1[0] =-3.8  ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6  ;
1583   fJetXMax1[3] =-2.7  ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7  ;
1584   fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016;
1585   fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
1586
1587
1588   //Photon fast reconstruction
1589   fResPara1       = 0.0255 ;    // GeV
1590   fResPara2       = 0.0272 ;
1591   fResPara3       = 0.0129 ;
1592
1593   fPosParaA      = 0.096 ;    // cm
1594   fPosParaB      = 0.229 ;
1595
1596   //Different cones and pt thresholds to construct the jet
1597
1598   fCone        = 0.3  ;
1599   fPtThreshold = 0.   ;
1600   fNCone       = 1    ;
1601   fNPt         = 1    ;
1602   fCones[0]    = 0.3  ; fNameCones[0]   = "03" ;
1603   fPtThres[0]  = 0.5  ; fNamePtThres[0] = "05" ;
1604
1605 }
1606
1607 //__________________________________________________________________________-
1608 Bool_t AliPHOSGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
1609   //Check if the opening angle of the candidate pairs is inside
1610   //our selection windowd
1611   Bool_t result = kFALSE;
1612   Double_t mpi0 = 0.1349766;
1613   Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
1614     +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
1615   Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
1616   Double_t min = 100. ;
1617   if(arg>0.)
1618     min = TMath::ACos(arg);
1619
1620   if((angle<max)&&(angle>=min))
1621     result = kTRUE;
1622
1623   return result;
1624 }
1625
1626 //__________________________________________________________________________-
1627 Bool_t AliPHOSGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj,
1628                                       const TString type ){
1629   //Check if the energy of the reconstructed jet is within an energy window
1630
1631   Double_t par[6];
1632   Double_t xmax[2];
1633   Double_t xmin[2];
1634
1635   Int_t iTPC = 0;
1636
1637   if(type == "TPC" && !fTPCCutsLikeEMCAL){
1638     iTPC = 3 ;//If(fTPCCutsLikeEMCAL) take jet energy cuts like EMCAL
1639   }
1640
1641
1642   if(!fHIJING){
1643     //Phythia alone, jets with pt_th > 0.2, r = 0.3
1644     par[0] = fJetE1[0]; par[1] = fJetE2[0];
1645     //Energy of the jet peak
1646     //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1647     par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1648     //Sigma  of the jet peak
1649     //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1650     par[4] = fBkgMean[0 + iTPC]; par[5] = fBkgRMS[0 + iTPC];
1651     //Parameters reserved for HIJING bkg.
1652     xmax[0] = fJetXMax1[0 + iTPC]; xmax[1] = fJetXMax2[0 + iTPC];
1653     xmin[0] = fJetXMin1[0 + iTPC]; xmin[1] = fJetXMin2[0 + iTPC];
1654     //Factor that multiplies sigma to obtain the best limits,
1655     //by observation, of mono jet ratios (ptjet/ptg)
1656     //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1657
1658   }
1659   else{
1660     if(ptg > fPtJetSelectionCut){
1661       //Phythia +HIJING with  pt_th > 2 GeV/c, r = 0.3
1662       par[0] = fJetE1[0]; par[1] = fJetE2[0];
1663       //Energy of the jet peak, same as in pp
1664       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1665       par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
1666       //Sigma  of the jet peak, same as in pp
1667       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1668       par[4] = fBkgMean[1 + iTPC]; par[5] = fBkgRMS[1 + iTPC];
1669       //Mean value and RMS of HIJING Bkg
1670       xmax[0] = fJetXMax1[1 + iTPC]; xmax[1] = fJetXMax2[1 + iTPC];
1671       xmin[0] = fJetXMin1[1 + iTPC]; xmin[1] = fJetXMin2[1 + iTPC];
1672       //Factor that multiplies sigma to obtain the best limits,
1673       //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg,
1674       //pt_th > 2 GeV, r = 0.3
1675       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1676
1677     }
1678     else{
1679       //Phythia + HIJING with  pt_th > 0.5 GeV/c, r = 0.3
1680       par[0] = fJetE1[1]; par[1] = fJetE2[1];
1681       //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3
1682       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
1683       par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
1684       //Sigma  of the jet peak, pt_th > 2 GeV/c, r = 0.3
1685       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
1686       par[4] = fBkgMean[2 + iTPC]; par[5] = fBkgRMS[2 + iTPC];
1687       //Mean value and RMS of HIJING Bkg in a 0.3 cone, pt > 2 GeV.
1688       xmax[0] = fJetXMax1[2 + iTPC]; xmax[1] = fJetXMax2[2 + iTPC];
1689       xmin[0] = fJetXMin1[2 + iTPC]; xmin[1] = fJetXMin2[2 + iTPC];
1690       //Factor that multiplies sigma to obtain the best limits,
1691       //by observation, of mono jet ratios (ptjet/ptg) mixed with HIJING Bkg,
1692       //pt_th > 2 GeV, r = 0.3
1693       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
1694
1695     }//If low pt jet in bkg
1696   }//if Bkg
1697
1698  //Calculate minimum and maximum limits of the jet ratio.
1699   Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
1700   Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
1701
1702   //Info("IsJetSeleted","%s : Limits min %f, max %f, ptg / ptj %f",
1703   //   type.Data(),min,max,ptj/ptg);
1704   if(( min < ptj/ptg ) && ( max > ptj/ptg))
1705     return kTRUE;
1706   else
1707     return kFALSE;
1708
1709 }
1710
1711 //____________________________________________________________________________
1712 void AliPHOSGammaJet::List() const
1713 {
1714   // List the histos
1715
1716   Info("List", "%d histograms found", fListHistos->GetEntries() ) ;
1717   TIter next(fListHistos) ;
1718   TH2F * h = 0 ;
1719   while ( (h = dynamic_cast<TH2F*>(next())) )
1720     Info("List", "%s", h->GetName()) ;
1721 }
1722
1723 //____________________________________________________________________________
1724 Double_t AliPHOSGammaJet::MakeEnergy(const Double_t energy)
1725 {
1726   // Smears the energy according to the energy dependent energy resolution.
1727   // A gaussian distribution is assumed
1728
1729   Double_t sigma  = SigmaE(energy) ;
1730   return  fRan.Gaus(energy, sigma) ;
1731
1732
1733 }
1734 //____________________________________________________________________________
1735 void AliPHOSGammaJet::MakeHistos()
1736 {
1737   // Create histograms to be saved in output file and
1738   // stores them in a TObjectArray
1739
1740   fOutputFile = new TFile(fOutputFileName, "recreate") ;
1741
1742   fListHistos = new TObjArray(10000) ;
1743
1744    // Histos gamma pt vs leading pt
1745
1746   TH1F * hPtSpectra  = new TH1F
1747     ("PtSpectra","p_{T i} vs p_{T #gamma}",200,0,200);
1748   hPtSpectra->SetXTitle("p_{T} (GeV/c)");
1749   fListHistos->Add(hPtSpectra) ;
1750
1751   //Histos ratio charged leading pt / gamma pt vs pt
1752
1753   TH2F * hChargeRatio  = new TH2F
1754     ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
1755      120,0,120,120,0,1);
1756   hChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
1757   hChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1758   fListHistos->Add(hChargeRatio) ;
1759
1760   TH2F * hTPCRatio  = new TH2F
1761     ("TPCRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
1762      120,0,120,120,0,1);
1763   hTPCRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
1764   hTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1765   fListHistos->Add(hTPCRatio) ;
1766
1767
1768   TH2F * hPi0Ratio  = new TH2F
1769     ("Pi0Ratio","p_{T leading  #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
1770      120,0,120,120,0,1);
1771   hPi0Ratio->SetYTitle("p_{T lead  #pi^{0}} /p_{T #gamma}");
1772   hPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
1773   fListHistos->Add(hPi0Ratio) ;
1774
1775   TH2F * hPhiGamma  = new TH2F
1776     ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
1777   hPhiGamma->SetYTitle("#phi");
1778   hPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
1779   fListHistos->Add(hPhiGamma) ;
1780
1781   TH2F * hEtaGamma  = new TH2F
1782     ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
1783   hEtaGamma->SetYTitle("#eta");
1784   hEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
1785   fListHistos->Add(hEtaGamma) ;
1786
1787 //   //Jet reconstruction check
1788 //   TH2F * hDeltaPhiJet  = new TH2F
1789 //     ("DeltaPhiJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1790 //      200,0,120,200,-1.5,1.5);
1791 //   hDeltaPhiJet->SetYTitle("#Delta #phi");
1792 //   hDeltaPhiJet->SetXTitle("p_{T #gamma} (GeV/c)");
1793 //   fListHistos->Add(hDeltaPhiJet) ;
1794
1795 //   TH2F * hDeltaPhiTPCJet  = new TH2F
1796 //     ("DeltaPhiTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1797 //      200,0,120,200,-1.5,1.5);
1798 //   hDeltaPhiTPCJet->SetYTitle("#Delta #phi");
1799 //   hDeltaPhiTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1800 //   fListHistos->Add(hDeltaPhiTPCJet) ;
1801
1802 //   TH2F * hDeltaEtaJet  = new TH2F
1803 //     ("DeltaEtaJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1804 //      200,0,120,200,-1.5,1.5);
1805 //   hDeltaEtaJet->SetYTitle("#Delta #phi");
1806 //   hDeltaEtaJet->SetXTitle("p_{T #gamma} (GeV/c)");
1807 //   fListHistos->Add(hDeltaEtaJet) ;
1808
1809 //   TH2F * hDeltaEtaTPCJet  = new TH2F
1810 //     ("DeltaEtaTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1811 //      200,0,120,200,-1.5,1.5);
1812 //   hDeltaEtaTPCJet->SetYTitle("#Delta #phi");
1813 //   hDeltaEtaTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1814 //   fListHistos->Add(hDeltaEtaTPCJet) ;
1815
1816 //   TH2F * hDeltaPtJet  = new TH2F
1817 //     ("DeltaPtJet","#phi_{jet} - #phi_{pyth jet} vs p_{T #gamma}",
1818 //      200,0,120,200,0.,100.);
1819 //   hDeltaPtJet->SetYTitle("#Delta #phi");
1820 //   hDeltaPtJet->SetXTitle("p_{T #gamma} (GeV/c)");
1821 //   fListHistos->Add(hDeltaPtJet) ;
1822
1823 //   TH2F * hDeltaPtTPCJet  = new TH2F
1824 //     ("DeltaPtTPCJet","#phi_{jet TPC} - #phi_{pyth jet} vs p_{T #gamma}",
1825 //      200,0,120,200,0.,100.);
1826 //   hDeltaPtTPCJet->SetYTitle("#Delta #phi");
1827 //   hDeltaPtTPCJet->SetXTitle("p_{T #gamma} (GeV/c)");
1828 //   fListHistos->Add(hDeltaPtTPCJet) ;
1829
1830   //
1831   TH2F * hDeltaPhiCharge  = new TH2F
1832     ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
1833      200,0,120,200,0,6.4);
1834   hDeltaPhiCharge->SetYTitle("#Delta #phi");
1835   hDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
1836   fListHistos->Add(hDeltaPhiCharge) ;
1837
1838   TH2F * hDeltaPhiTPC  = new TH2F
1839     ("DeltaPhiTPC","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
1840      200,0,120,200,0,6.4);
1841   hDeltaPhiTPC->SetYTitle("#Delta #phi");
1842   hDeltaPhiTPC->SetXTitle("p_{T #gamma} (GeV/c)");
1843   fListHistos->Add(hDeltaPhiTPC) ;
1844   TH2F * hDeltaPhiPi0  = new TH2F
1845     ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
1846      200,0,120,200,0,6.4);
1847   hDeltaPhiPi0->SetYTitle("#Delta #phi");
1848   hDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1849   fListHistos->Add(hDeltaPhiPi0) ;
1850
1851   TH2F * hDeltaEtaCharge  = new TH2F
1852     ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
1853      200,0,120,200,-2,2);
1854   hDeltaEtaCharge->SetYTitle("#Delta #eta");
1855   hDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
1856   fListHistos->Add(hDeltaEtaCharge) ;
1857
1858   TH2F * hDeltaEtaTPC  = new TH2F
1859     ("DeltaEtaTPC","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
1860      200,0,120,200,-2,2);
1861   hDeltaEtaTPC->SetYTitle("#Delta #eta");
1862   hDeltaEtaTPC->SetXTitle("p_{T #gamma} (GeV/c)");
1863   fListHistos->Add(hDeltaEtaTPC) ;
1864
1865   TH2F * hDeltaEtaPi0  = new TH2F
1866     ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
1867      200,0,120,200,-2,2);
1868   hDeltaEtaPi0->SetYTitle("#Delta #eta");
1869   hDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
1870   fListHistos->Add(hDeltaEtaPi0) ;
1871
1872   if(fOptFast){
1873
1874     TH2F * hAnglePair  = new TH2F
1875       ("AnglePair",
1876        "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}",
1877        200,0,50,200,0,0.2);
1878     hAnglePair->SetYTitle("Angle (rad)");
1879     hAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1880     fListHistos->Add(hAnglePair) ;
1881
1882     TH2F * hAnglePairAccepted  = new TH2F
1883       ("AnglePairAccepted",
1884        "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}, both #gamma in eta<0.7, inside window",
1885        200,0,50,200,0,0.2);
1886     hAnglePairAccepted->SetYTitle("Angle (rad)");
1887     hAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1888     fListHistos->Add(hAnglePairAccepted) ;
1889
1890     TH2F * hAnglePairNoCut  = new TH2F
1891       ("AnglePairNoCut",
1892        "Angle between all #gamma pair vs p_{T  #pi^{0}}",200,0,50,200,0,0.2);
1893     hAnglePairNoCut->SetYTitle("Angle (rad)");
1894     hAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1895     fListHistos->Add(hAnglePairNoCut) ;
1896
1897     TH2F * hAnglePairLeadingCut  = new TH2F
1898       ("AnglePairLeadingCut",
1899        "Angle between all #gamma pair that have a good phi and pt vs p_{T  #pi^{0}}",
1900        200,0,50,200,0,0.2);
1901     hAnglePairLeadingCut->SetYTitle("Angle (rad)");
1902     hAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1903     fListHistos->Add(hAnglePairLeadingCut) ;
1904
1905     TH2F * hAnglePairAngleCut  = new TH2F
1906       ("AnglePairAngleCut",
1907        "Angle between all #gamma pair (angle + leading cut) vs p_{T  #pi^{0}}"
1908        ,200,0,50,200,0,0.2);
1909     hAnglePairAngleCut->SetYTitle("Angle (rad)");
1910     hAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1911     fListHistos->Add(hAnglePairAngleCut) ;
1912
1913     TH2F * hAnglePairAllCut  = new TH2F
1914       ("AnglePairAllCut",
1915        "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T  #pi^{0}}"
1916        ,200,0,50,200,0,0.2);
1917     hAnglePairAllCut->SetYTitle("Angle (rad)");
1918     hAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1919     fListHistos->Add(hAnglePairAllCut) ;
1920
1921     TH2F * hAnglePairLeading  = new TH2F
1922       ("AnglePairLeading",
1923        "Angle between all #gamma pair finally selected vs p_{T  #pi^{0}}",
1924        200,0,50,200,0,0.2);
1925     hAnglePairLeading->SetYTitle("Angle (rad)");
1926     hAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
1927     fListHistos->Add(hAnglePairLeading) ;
1928
1929
1930     TH2F * hInvMassPairNoCut  = new TH2F
1931       ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
1932        120,0,120,360,0,0.5);
1933     hInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1934     hInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
1935     fListHistos->Add(hInvMassPairNoCut) ;
1936
1937     TH2F * hInvMassPairLeadingCut  = new TH2F
1938       ("InvMassPairLeadingCut",
1939        "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
1940        120,0,120,360,0,0.5);
1941     hInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1942     hInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
1943     fListHistos->Add(hInvMassPairLeadingCut) ;
1944
1945     TH2F * hInvMassPairAngleCut  = new TH2F
1946       ("InvMassPairAngleCut",
1947        "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
1948        120,0,120,360,0,0.5);
1949     hInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1950     hInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
1951     fListHistos->Add(hInvMassPairAngleCut) ;
1952
1953
1954     TH2F * hInvMassPairAllCut  = new TH2F
1955       ("InvMassPairAllCut",
1956        "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
1957        120,0,120,360,0,0.5);
1958     hInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
1959     hInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
1960     fListHistos->Add(hInvMassPairAllCut) ;
1961
1962     TH2F * hInvMassPairLeading  = new TH2F
1963       ("InvMassPairLeading",
1964        "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
1965        120,0,120,360,0,0.5);
1966     hInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
1967     hInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
1968     fListHistos->Add(hInvMassPairLeading) ;
1969   }
1970
1971   //Count
1972
1973   TH1F * hNGamma  = new TH1F("NGamma","Number of #gamma over PHOS",240,0,120);
1974   hNGamma->SetYTitle("N");
1975   hNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
1976   fListHistos->Add(hNGamma) ;
1977
1978   TH1F * hNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000);
1979   hNBkg->SetYTitle("counts");
1980   hNBkg->SetXTitle("N");
1981   fListHistos->Add(hNBkg) ;
1982
1983   TH2F * hNLeading  = new TH2F
1984     ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120);
1985   hNLeading->SetYTitle("p_{T charge} (GeV/c)");
1986   hNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
1987   fListHistos->Add(hNLeading) ;
1988
1989
1990   TH1F * hN  = new TH1F("NJet","Accepted jets",240,0,120);
1991   hN->SetYTitle("N");
1992   hN->SetXTitle("p_{T #gamma}(GeV/c)");
1993   fListHistos->Add(hN) ;
1994
1995
1996   //Ratios and Pt dist of reconstructed (not selected) jets
1997   //Jet
1998   TH2F * hJetRatio  = new TH2F
1999     ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
2000      240,0,120,200,0,10);
2001   hJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2002   hJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2003   fListHistos->Add(hJetRatio) ;
2004
2005
2006   TH2F * hJetPt  = new TH2F
2007     ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
2008   hJetPt->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2009   hJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
2010   fListHistos->Add(hJetPt) ;
2011
2012
2013   //Bkg
2014
2015   TH2F * hBkgRatio  = new TH2F
2016     ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
2017      240,0,120,200,0,10);
2018   hBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
2019   hBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2020   fListHistos->Add(hBkgRatio) ;
2021
2022
2023   TH2F * hBkgPt  = new TH2F
2024     ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
2025   hBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
2026   hBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
2027   fListHistos->Add(hBkgPt) ;
2028
2029
2030   //Jet Distributions
2031
2032   TH2F * hJetFragment  =
2033     new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
2034              240,0.,120.,1000,0.,1.2);
2035   hJetFragment->SetYTitle("x_{T}");
2036   hJetFragment->SetXTitle("p_{T #gamma}");
2037   fListHistos->Add(hJetFragment) ;
2038
2039   TH2F * hBkgFragment  = new TH2F
2040     ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
2041      240,0.,120.,1000,0.,1.2);
2042   hBkgFragment->SetYTitle("x_{T}");
2043   hBkgFragment->SetXTitle("p_{T #gamma}");
2044   fListHistos->Add(hBkgFragment) ;
2045
2046   TH2F * hJetPtDist  =
2047     new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2048   hJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2049   fListHistos->Add(hJetPtDist) ;
2050
2051   TH2F * hBkgPtDist  = new TH2F
2052     ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2053   hBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2054   fListHistos->Add(hBkgPtDist) ;
2055
2056
2057   if(fOnlyCharged){
2058     //Counts
2059     TH1F * hNBkgTPC  = new TH1F
2060       ("NBkgTPC","TPC bkg multiplicity ",9000,0,9000);
2061     hNBkgTPC->SetYTitle("counts");
2062     hNBkgTPC->SetXTitle("N");
2063     fListHistos->Add(hNBkgTPC) ;
2064
2065     TH2F * hNTPCLeading  = new TH2F
2066       ("NTPCLeading","Accepted TPC jet leading",240,0,120,240,0,120);
2067     hNTPCLeading->SetYTitle("p_{T charge} (GeV/c)");
2068     hNTPCLeading->SetXTitle("p_{T #gamma}(GeV/c)");
2069     fListHistos->Add(hNTPCLeading) ;
2070
2071     TH1F * hNTPC  = new TH1F("NTPCJet","Number of TPC jets",240,0,120);
2072     hNTPC->SetYTitle("N");
2073     hNTPC->SetXTitle("p_{T #gamma}(GeV/c)");
2074     fListHistos->Add(hNTPC) ;
2075
2076     TH2F * hJetTPCRatio  = new TH2F
2077       ("JetTPCRatio", "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}",
2078        240,0,120,200,0,10);
2079     hJetTPCRatio->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2080     hJetTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2081     fListHistos->Add(hJetTPCRatio) ;
2082
2083     TH2F * hBkgTPCRatio  = new TH2F
2084       ("BkgTPCRatio","p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}",
2085        240,0,120,200,0,10);
2086     hBkgTPCRatio->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2087     hBkgTPCRatio->SetXTitle("p_{T #gamma} (GeV/c)");
2088     fListHistos->Add(hBkgTPCRatio) ;
2089
2090     TH2F * hJetTPCPt  = new TH2F
2091       ("JetTPCPt", "p_{T jet lead TPC} vs p_{T #gamma}",240,0,120,400,0,200);
2092     hJetTPCPt->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2093     hJetTPCPt->SetXTitle("p_{T #gamma} (GeV/c)");
2094     fListHistos->Add(hJetTPCPt) ;
2095
2096     TH2F * hBkgTPCPt  = new TH2F
2097       ("BkgTPCPt", "p_{T bkg lead TPC} vs p_{T #gamma}",240,0,120,400,0,200);
2098     hBkgTPCPt->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2099     hBkgTPCPt->SetXTitle("p_{T #gamma} (GeV/c)");
2100     fListHistos->Add(hBkgTPCPt) ;
2101
2102     //JetDistributions
2103
2104     TH2F * hJetTPCFragment  =
2105       new TH2F("JetTPCFragment","x = p_{T i charged}/p_{T #gamma}",
2106                240,0.,120.,1000,0.,1.2);
2107     hJetTPCFragment->SetYTitle("x_{T}");
2108     hJetTPCFragment->SetXTitle("p_{T #gamma}");
2109     fListHistos->Add(hJetTPCFragment) ;
2110
2111     TH2F * hBkgTPCFragment  = new TH2F
2112       ("BkgTPCFragment","x = p_{T i charged}/p_{T #gamma}",
2113        240,0.,120.,1000,0.,1.2);
2114     hBkgTPCFragment->SetYTitle("x_{T}");
2115     hBkgTPCFragment->SetXTitle("p_{T #gamma}");
2116     fListHistos->Add(hBkgTPCFragment) ;
2117
2118
2119     TH2F * hJetTPCPtDist  = new TH2F("JetTPCPtDist",
2120        "x = p_{T i charged}",240,0.,120.,400,0.,200.);
2121     hJetTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2122     fListHistos->Add(hJetTPCPtDist) ;
2123
2124     TH2F * hBkgTPCPtDist  = new TH2F
2125       ("BkgTPCPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.);
2126     hBkgTPCPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
2127     fListHistos->Add(hBkgTPCPtDist) ;
2128
2129   }
2130
2131
2132   if(fAnyConeOrPt){
2133     //If we want to study the jet for different cones and pt. Old version
2134
2135     TH2F * hJetRatios[5][5];
2136     TH2F * hJetTPCRatios[5][5];
2137
2138     TH2F * hJetPts[5][5];
2139     TH2F * hJetTPCPts[5][5];
2140
2141     TH2F * hBkgRatios[5][5];
2142     TH2F * hBkgTPCRatios[5][5];
2143
2144     TH2F * hBkgPts[5][5];
2145     TH2F * hBkgTPCPts[5][5];
2146
2147     TH2F * hNLeadings[5][5];
2148     TH2F * hNTPCLeadings[5][5];
2149
2150     TH1F * hNs[5][5];
2151     TH1F * hNTPCs[5][5];
2152
2153     TH1F * hNBkgs[5][5];
2154     TH1F * hNBkgTPCs[5][5];
2155
2156     TH2F * hJetFragments[5][5];
2157     TH2F * hBkgFragments[5][5];
2158     TH2F * hJetPtDists[5][5];
2159     TH2F * hBkgPtDists[5][5];
2160
2161     TH2F * hJetTPCFragments[5][5];
2162     TH2F * hBkgTPCFragments[5][5];
2163     TH2F * hJetTPCPtDists[5][5];
2164     TH2F * hBkgTPCPtDists[5][5];
2165
2166
2167     for(Int_t icone = 0; icone<fNCone; icone++){
2168       for(Int_t ipt = 0; ipt<fNPt;ipt++){
2169         //Jet Pt / Gamma Pt
2170
2171         //Jet
2172
2173         hJetRatios[icone][ipt]  = new TH2F
2174           ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2175          "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2176            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2177            240,0,120,200,0,10);
2178         hJetRatios[icone][ipt]->
2179           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2180         hJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2181         fListHistos->Add(hJetRatios[icone][ipt]) ;
2182
2183
2184         hJetPts[icone][ipt]  = new TH2F
2185           ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2186            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2187            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2188            240,0,120,400,0,200);
2189         hJetPts[icone][ipt]->
2190           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2191         hJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2192         fListHistos->Add(hJetPts[icone][ipt]) ;
2193
2194
2195         //Bkg
2196         hBkgRatios[icone][ipt]  = new TH2F
2197           ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2198            "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2199            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2200            240,0,120,200,0,10);
2201         hBkgRatios[icone][ipt]->
2202           SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
2203         hBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2204         fListHistos->Add(hBkgRatios[icone][ipt]) ;
2205
2206
2207
2208         hBkgPts[icone][ipt]  = new TH2F
2209           ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2210            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
2211            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2212            240,0,120,400,0,200);
2213         hBkgPts[icone][ipt]->
2214           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
2215         hBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2216         fListHistos->Add(hBkgPts[icone][ipt]) ;
2217
2218
2219         if(fOnlyCharged){
2220
2221           hJetTPCRatios[icone][ipt]  = new TH2F
2222             ("JetTPCRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2223              "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2224              +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2225              240,0,120,200,0,10);
2226           hJetTPCRatios[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2227           hJetTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2228           fListHistos->Add(hJetTPCRatios[icone][ipt]) ;
2229
2230           hBkgTPCRatios[icone][ipt]  = new TH2F
2231             ("BkgTPCRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2232              "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2233              +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2234              240,0,120,200,0,10);
2235           hBkgTPCRatios[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2236           hBkgTPCRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2237           fListHistos->Add(hBkgTPCRatios[icone][ipt]) ;
2238
2239           hJetTPCPts[icone][ipt]  = new TH2F
2240             ("JetTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2241              "p_{T jet lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2242              +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2243              240,0,120,400,0,200);
2244           hJetTPCPts[icone][ipt]->SetYTitle("p_{T jet lead TPC}/p_{T #gamma}");
2245           hJetTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2246           fListHistos->Add(hJetTPCPts[icone][ipt]) ;
2247
2248           hBkgTPCPts[icone][ipt]  = new TH2F
2249             ("BkgTPCPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2250              "p_{T bkg lead TPC}/p_{T #gamma} vs p_{T #gamma}, cone ="
2251              +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
2252              240,0,120,400,0,200);
2253           hBkgTPCPts[icone][ipt]->SetYTitle("p_{T bkg lead TPC}/p_{T #gamma}");
2254           hBkgTPCPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2255           fListHistos->Add(hBkgTPCPts[icone][ipt]) ;
2256
2257         }
2258
2259         //Counts
2260         hNBkgs[icone][ipt]  = new TH1F
2261           ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2262            "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
2263            +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
2264         hNBkgs[icone][ipt]->SetYTitle("counts");
2265         hNBkgs[icone][ipt]->SetXTitle("N");
2266         fListHistos->Add(hNBkgs[icone][ipt]) ;
2267
2268         hNBkgTPCs[icone][ipt]  = new TH1F
2269           ("NBkgTPCCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2270            "bkg multiplicity cone ="+fNameCones[icone]+", pt>"
2271            +fNamePtThres[ipt]+" GeV/c",9000,0,9000);
2272         hNBkgTPCs[icone][ipt]->SetYTitle("counts");
2273         hNBkgTPCs[icone][ipt]->SetXTitle("N");
2274         fListHistos->Add(hNBkgTPCs[icone][ipt]) ;
2275
2276
2277         hNLeadings[icone][ipt]  = new TH2F
2278           ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2279            "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>"
2280            +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
2281         hNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
2282         hNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2283         fListHistos->Add(hNLeadings[icone][ipt]) ;
2284
2285         hNs[icone][ipt]  = new TH1F
2286           ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2287            "Number of neutral jets, cone ="+fNameCones[icone]+", pt>"
2288            +fNamePtThres[ipt]+" GeV/c",120,0,120);
2289         hNs[icone][ipt]->SetYTitle("N");
2290         hNs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2291         fListHistos->Add(hNs[icone][ipt]) ;
2292
2293         //Fragmentation Function
2294         hJetFragments[icone][ipt]  = new TH2F
2295           ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2296            "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2297            +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2298         hJetFragments[icone][ipt]->SetYTitle("x_{T}");
2299         hJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2300         fListHistos->Add(hJetFragments[icone][ipt]) ;
2301
2302         hBkgFragments[icone][ipt]  = new TH2F
2303           ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2304            "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2305            +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2306         hBkgFragments[icone][ipt]->SetYTitle("x_{T}");
2307         hBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2308         fListHistos->Add(hBkgFragments[icone][ipt]) ;
2309
2310
2311         //Jet particle distribution
2312
2313         hJetPtDists[icone][ipt]  = new TH2F
2314           ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2315            "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
2316            " GeV/c",120,0.,120.,120,0.,120.);
2317         hJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2318         fListHistos->Add(hJetPtDists[icone][ipt]) ;
2319
2320         hBkgPtDists[icone][ipt]  = new TH2F
2321           ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2322            "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
2323            " GeV/c",120,0.,120.,120,0.,120.);
2324         hBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2325         fListHistos->Add(hBkgPtDists[icone][ipt]) ;
2326
2327
2328         if(fOnlyCharged){
2329           hNTPCLeadings[icone][ipt]  = new TH2F
2330             ("NTPCLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2331              "p_{T #gamma} vs p_{T charge} cone ="+fNameCones[icone]+", pt>"
2332              +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120);
2333           hNTPCLeadings[icone][ipt]->SetYTitle("p_{T charge} (GeV/c)");
2334           hNTPCLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2335           fListHistos->Add(hNTPCLeadings[icone][ipt]) ;
2336
2337           hNTPCs[icone][ipt]  = new TH1F
2338             ("NTPCJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2339              "Number of charged jets, cone ="+fNameCones[icone]+", pt>"
2340              +fNamePtThres[ipt]+" GeV/c",120,0,120);
2341           hNTPCs[icone][ipt]->SetYTitle("N");
2342           hNTPCs[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
2343           fListHistos->Add(hNTPCs[icone][ipt]) ;
2344
2345           hJetTPCFragments[icone][ipt]  = new TH2F
2346             ("JetTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2347              "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2348              +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2349           hJetTPCFragments[icone][ipt]->SetYTitle("x_{T}");
2350           hJetTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2351           fListHistos->Add(hJetTPCFragments[icone][ipt]) ;
2352
2353           hBkgTPCFragments[icone][ipt]  = new TH2F
2354             ("BkgTPCFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2355              "x = p_{T i charged}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>"
2356              +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2);
2357           hBkgTPCFragments[icone][ipt]->SetYTitle("x_{T}");
2358           hBkgTPCFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
2359           fListHistos->Add(hBkgTPCFragments[icone][ipt]) ;
2360
2361           hJetTPCPtDists[icone][ipt]  = new TH2F
2362             ("JetTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2363              "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>"
2364              +fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.);
2365           hJetTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2366           fListHistos->Add(hJetTPCPtDists[icone][ipt]) ;
2367
2368           hBkgTPCPtDists[icone][ipt]  = new TH2F
2369             ("BkgTPCPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
2370              "x = p_{T i charged}, cone ="+fNameCones[icone]+", pt>" +
2371              fNamePtThres[ipt]+" GeV/c",120,0.,120.,120,0.,120.);
2372           hBkgTPCPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
2373           fListHistos->Add(hBkgTPCPtDists[icone][ipt]) ;
2374         }
2375       }//ipt
2376     } //icone
2377   }//If we want to study any cone or pt threshold
2378 }
2379
2380
2381 //____________________________________________________________________________
2382 void AliPHOSGammaJet::MakeJet(TClonesArray * pl,
2383                               Double_t ptg, Double_t phig,
2384                               Double_t ptl, Double_t  phil, Double_t  etal,
2385                               TString conf, TLorentzVector & jet)
2386 {
2387   //Fill the jet with the particles around the leading particle with
2388   //R=fCone and pt_th = fPtThres. Calculate the energy of the jet and
2389   //check if we select it. Fill jet histograms
2390   Float_t ptcut = 0. ;
2391   if(fHIJING){
2392     if(ptg > fPtJetSelectionCut)  ptcut = 2. ;
2393     else                          ptcut = 0.5;
2394   }
2395
2396   TClonesArray * jetList = new TClonesArray("TParticle",1000);
2397   TClonesArray * bkgList = new TClonesArray("TParticle",1000);
2398
2399   TLorentzVector bkg(0,0,0,0);
2400   TLorentzVector lv (0,0,0,0);
2401
2402   Double_t ptjet = 0.0;
2403   Double_t ptbkg = 0.0;
2404   Int_t n0 = 0;
2405   Int_t n1 = 0;
2406   Bool_t b1 = kFALSE;
2407   Bool_t b0 = kFALSE;
2408
2409   TIter next(pl) ;
2410   TParticle * particle = 0 ;
2411
2412   while ( (particle = dynamic_cast<TParticle*>(next())) ) {
2413
2414     b0 = kFALSE;
2415     b1 = kFALSE;
2416
2417     //Particles in jet
2418     SetJet(particle, b0, fCone, etal, phil) ;
2419
2420     if(b0){
2421       new((*jetList)[n0++]) TParticle(*particle) ;
2422       particle->Momentum(lv);
2423       if(particle->Pt() > ptcut ){
2424         jet+=lv;
2425         ptjet+=particle->Pt();
2426       }
2427     }
2428
2429     //Background around (phi_gamma-pi, eta_leading)
2430     SetJet(particle, b1, fCone,etal, phig) ;
2431
2432     if(b1) {
2433       new((*bkgList)[n1++]) TParticle(*particle) ;
2434       if(particle->Pt() > ptcut ){
2435         bkg+=lv;
2436         ptbkg+=particle->Pt();
2437       }
2438     }
2439   }
2440
2441   ptjet = jet.Pt();
2442   ptbkg = bkg.Pt();
2443
2444   if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
2445     Info("MakeJet","Gamma   pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg);
2446
2447
2448   //Fill histograms
2449
2450   Double_t rat     = ptjet/ptg ;
2451   Double_t ratbkg  = ptbkg/ptg ;
2452
2453   dynamic_cast<TH2F*>
2454     (fListHistos->FindObject("Jet"+conf+"Ratio"))->Fill(ptg,rat);
2455   dynamic_cast<TH2F*>
2456     (fListHistos->FindObject("Jet"+conf+"Pt"))   ->Fill(ptg,ptjet);
2457   dynamic_cast<TH2F*>
2458     (fListHistos->FindObject("Bkg"+conf+"Ratio"))->Fill(ptg,ratbkg);
2459   dynamic_cast<TH2F*>
2460     (fListHistos->FindObject("Bkg"+conf+"Pt"))   ->Fill(ptg,ptbkg);
2461
2462
2463   if(IsJetSelected(ptg,ptjet,conf) || fSelect){
2464     if(strstr(fOptionGJ,"deb") || strstr(fOptionGJ,"deb all"))
2465       Info("MakeJet","JetSelected");
2466     dynamic_cast<TH1F*>(fListHistos->FindObject("N"+conf+"Jet"))->
2467       Fill(ptg);
2468     dynamic_cast<TH2F*>(fListHistos->FindObject("N"+conf+"Leading"))
2469       ->Fill(ptg,ptl);
2470     FillJetHistos(jetList, ptg, conf, "Jet");
2471     FillJetHistos(bkgList, ptg, conf, "Bkg");
2472   }
2473
2474   jetList ->Delete();
2475   bkgList ->Delete();
2476 }
2477
2478 //____________________________________________________________________________
2479 void AliPHOSGammaJet::MakeJetAnyConeOrPt(TClonesArray * pl, Double_t ptg,
2480                                          Double_t  phig, Double_t ptl,
2481                                          Double_t  phil, Double_t  etal,
2482                                          TString conf){
2483
2484   //Fill the jet with the particles around the leading particle with
2485   //R=fCone(i) and pt_th = fPtThres(i). Calculate the energy of the jet and
2486   //check if we select it. Fill jet i histograms
2487
2488   TClonesArray * jetList = new TClonesArray("TParticle",1000);
2489   TClonesArray * bkgList = new TClonesArray("TParticle",1000);
2490
2491   Double_t ptjet  = 0.0;
2492   Double_t ptbkg  = 0.0;
2493
2494   Int_t n1  = 0;
2495   Int_t n0  = 0;
2496   Bool_t b1 = kFALSE;
2497   Bool_t b0 = kFALSE;
2498
2499   //Create as many jets as cones and pt thresholds are defined
2500   Double_t maxcut = fJetRatioMaxCut;
2501   Double_t mincut = fJetRatioMinCut;
2502
2503   if(conf == "TPC"){
2504     maxcut = fJetTPCRatioMaxCut;
2505     mincut = fJetTPCRatioMinCut;
2506   }
2507
2508   Double_t ratjet  = 0;
2509   Double_t ratbkg  = 0;
2510
2511   for(Int_t icone = 0; icone<fNCone; icone++) {
2512     for(Int_t ipt = 0; ipt<fNPt;ipt++) {
2513
2514       TString cone  = fNameCones[icone]  ;
2515       TString ptcut = fNamePtThres[ipt] ;
2516
2517       TIter next(pl) ;
2518       TParticle * particle = 0 ;
2519
2520       ptjet  = 0 ;
2521       ptbkg = 0 ;
2522
2523       while ( (particle = dynamic_cast<TParticle*>(next())) ) {
2524         b1 = kFALSE;
2525         b0 = kFALSE;
2526
2527         SetJet(particle, b0, fCones[icone],etal, phil) ;
2528         SetJet(particle, b1, fCones[icone],etal, phig) ;
2529
2530         if(b0){
2531           new((*jetList)[n0++]) TParticle(*particle) ;
2532           if(particle->Pt() > fPtThres[ipt] )
2533             ptjet+=particle->Pt();
2534         }
2535         if(b1) {
2536           new((*bkgList)[n1++]) TParticle(*particle) ;
2537           if(particle->Pt() > fPtThres[ipt] )
2538             ptbkg+=particle->Pt();
2539         }
2540
2541       }
2542
2543       //Fill histograms
2544       if(ptjet > 0.) {
2545
2546         if(strstr(fOptionGJ,"deb")){
2547           Info("MakeJetAnyPt","cone %f, ptcut %f",fCones[icone],fPtThres[ipt]);
2548           Info("MakeJetAnyPt","pT: Gamma %f, Jet %f, Bkg  %f",ptg,ptjet,ptbkg);
2549         }
2550
2551         ratjet  = ptjet /ptg;
2552         ratbkg  = ptbkg/ptg;
2553
2554         dynamic_cast<TH2F*>
2555           (fListHistos->FindObject("Jet"+conf+"RatioCone"+cone+"Pt"+ptcut))
2556           ->Fill(ptg,ratjet);
2557         dynamic_cast<TH2F*>
2558           (fListHistos->FindObject("Jet"+conf+"PtCone"+cone+"Pt"+ptcut))
2559           ->Fill(ptg,ptjet);
2560
2561         dynamic_cast<TH2F*>
2562           (fListHistos->FindObject("Bkg"+conf+"RatioCone"+cone+"Pt"+ptcut))
2563           ->Fill(ptg,ratbkg);
2564
2565         dynamic_cast<TH2F*>
2566           (fListHistos->FindObject("Bkg"+conf+"PtCone"+cone+"Pt"+ptcut))
2567           ->Fill(ptg,ptbkg);
2568
2569
2570         //Select Jet
2571         if((ratjet < maxcut) && (ratjet > mincut)){
2572
2573           dynamic_cast<TH1F*>
2574             (fListHistos->FindObject("N"+conf+"JetCone"+cone+"Pt"+ptcut))->
2575             Fill(ptg);
2576           dynamic_cast<TH2F*>
2577             (fListHistos->FindObject("N"+conf+"LeadingCone"+cone+"Pt"+ptcut))
2578             ->Fill(ptg,ptl);
2579
2580             FillJetHistosAnyConeOrPt
2581               (jetList,ptg,conf,"Jet",fNameCones[icone],fNamePtThres[ipt]);
2582             FillJetHistosAnyConeOrPt
2583               (bkgList,ptg,conf,"Bkg",fNameCones[icone],fNamePtThres[ipt]);
2584
2585         }
2586       } //ptjet > 0
2587       jetList ->Delete();
2588       bkgList ->Delete();
2589     }//for pt threshold
2590   }// for cone
2591 }
2592
2593 //____________________________________________________________________________
2594 void  AliPHOSGammaJet::MakePhoton(TLorentzVector & particle)
2595 {
2596   //Fast reconstruction for photons
2597   Double_t energy = particle.E()  ;
2598   Double_t modenergy = MakeEnergy(energy) ;
2599   //Info("MakePhoton","Energy %f, Modif %f",energy,modenergy);
2600
2601   // get the detected direction
2602     TVector3 pos = particle.Vect();
2603     pos*=460./energy;
2604     TVector3 modpos = MakePosition(energy, pos) ;
2605     modpos *= modenergy / 460.;
2606
2607     Float_t modtheta = modpos.Theta();
2608     Float_t modphi   = modpos.Phi();
2609
2610     // Set the modified 4-momentum of the reconstructed particle
2611     Float_t py = modenergy*TMath::Sin(modphi)*TMath::Sin(modtheta);
2612     Float_t px = modenergy*TMath::Cos(modphi)*TMath::Sin(modtheta);
2613     Float_t pz = modenergy*TMath::Cos(modtheta);
2614
2615     particle.SetPxPyPzE(px,py,pz,modenergy);
2616
2617 }
2618
2619 //____________________________________________________________________________
2620 TVector3 AliPHOSGammaJet::MakePosition(const Double_t energy, const TVector3 pos)
2621 {
2622   // Smears the impact position according to the energy dependent position resolution
2623   // A gaussian position distribution is assumed
2624
2625   TVector3 newpos ;
2626
2627   Double_t sigma = SigmaP(energy) ;
2628   Double_t x = fRan.Gaus( pos.X(), sigma ) ;
2629   Double_t z = fRan.Gaus( pos.Z(), sigma ) ;
2630   Double_t y = pos.Y() ;
2631
2632   newpos.SetX(x) ;
2633   newpos.SetY(y) ;
2634   newpos.SetZ(z) ;
2635
2636   // Info("MakePosition","Theta dif %f",pos.Theta()-newpos.Theta());
2637 //   Info("MakePosition","Phi   dif %f",pos.Phi()-newpos.Phi());
2638   return newpos ;
2639 }
2640
2641 //____________________________________________________________________________
2642 void AliPHOSGammaJet::Pi0Decay(Double_t mPi0, TLorentzVector &p0,
2643                                TLorentzVector &p1, TLorentzVector &p2, Double_t &angle) {
2644   // Perform isotropic decay pi0 -> 2 photons
2645   // p0 is pi0 4-momentum (inut)
2646   // p1 and p2 are photon 4-momenta (output)
2647   //  cout<<"Boost vector"<<endl;
2648   TVector3 b = p0.BoostVector();
2649   //cout<<"Parameters"<<endl;
2650   //Double_t mPi0   = p0.M();
2651   Double_t phi    = TMath::TwoPi() * gRandom->Rndm();
2652   Double_t cosThe = 2 * gRandom->Rndm() - 1;
2653   Double_t cosPhi = TMath::Cos(phi);
2654   Double_t sinPhi = TMath::Sin(phi);
2655   Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
2656   Double_t ePi0   = mPi0/2.;
2657   //cout<<"ePi0 "<<ePi0<<endl;
2658   //cout<<"Components"<<endl;
2659   p1.SetPx(+ePi0*cosPhi*sinThe);
2660   p1.SetPy(+ePi0*sinPhi*sinThe);
2661   p1.SetPz(+ePi0*cosThe);
2662   p1.SetE(ePi0);
2663   //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
2664   //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
2665   p2.SetPx(-ePi0*cosPhi*sinThe);
2666   p2.SetPy(-ePi0*sinPhi*sinThe);
2667   p2.SetPz(-ePi0*cosThe);
2668   p2.SetE(ePi0);
2669   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
2670   //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
2671   //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
2672   p1.Boost(b);
2673   //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
2674   p2.Boost(b);
2675   //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
2676   //cout<<"angle"<<endl;
2677   angle = p1.Angle(p2.Vect());
2678   //cout<<angle<<endl;
2679 }
2680
2681 //____________________________________________________________________________
2682 void AliPHOSGammaJet::Plot(TString what, Option_t * option) const
2683 {
2684   //Plot some relevant histograms of the analysis
2685   TH2F * h = dynamic_cast<TH2F*>(fOutputFile->Get(what));
2686   if(h){
2687     h->Draw();
2688   }
2689   else if (what == "all") {
2690   TCanvas * can = new TCanvas("GammaJet", "Gamma-Jet Study",10,40,1600,1200);
2691   can->cd() ;
2692   can->Range(0,0,22,20);
2693   TPaveLabel *pl1 = new TPaveLabel(1,18,20,19.5,"Titre","br");
2694   pl1->SetFillColor(18);
2695   pl1->SetTextFont(32);
2696   pl1->SetTextColor(49);
2697   pl1->Draw();
2698   Int_t index ;
2699   TPad * pad = 0;
2700   Float_t begx = -0.29, begy = 0., endx = 0., endy = 0.30 ;
2701   for (index = 0 ; index < fListHistos->GetEntries() ; index++) {
2702     TString name("pad") ;
2703     name += index ;
2704     begx += 0.30 ;
2705     endx += 0.30 ;
2706     if (begx >= 1.0 || endx >= 1.0) {
2707       begx = 0.01 ;
2708       endx = 0.30 ;
2709       begy += 0.30 ;
2710       endy += 0.30 ;
2711     }
2712     printf("%f %f %f %f \n", begx, begy, endx, endy) ;
2713     pad = new TPad(name,"This is a pad",begx,begy,endx,endy,33);
2714     pad->Draw();
2715     pad->cd() ;
2716     fListHistos->At(index)->Draw(option) ;
2717     pad->Modified() ;
2718     pad->Update() ;
2719     can->cd() ;
2720   }
2721
2722   }
2723   else{
2724     Info("Draw", "Histogram %s does not exist or unknown option", what.Data()) ;
2725      fOutputFile->ls();
2726   }
2727 }
2728
2729 //____________________________________________________________________________
2730 void AliPHOSGammaJet::Print(const Option_t * opt) const
2731 {
2732
2733   //Print some relevant parameters set for the analysis
2734   if(! opt)
2735     return;
2736
2737   Info("Print", "%s %s", GetName(), GetTitle() ) ;
2738
2739   printf("Eta cut           : %f\n", fEtaCut) ;
2740   printf("D phi max cut     : %f\n", fPhiMaxCut) ;
2741   printf("D phi min cut     : %f\n", fPhiMinCut) ;
2742   printf("Leading Ratio max cut     : %f\n", fRatioMaxCut) ;
2743   printf("Leading Ratio min cut     : %f\n", fRatioMinCut) ;
2744   printf("Jet Ratio max cut     : %f\n", fJetRatioMaxCut) ;
2745   printf("Jet Ratio min cut     : %f\n", fJetRatioMinCut) ;
2746   printf("Jet TPC Ratio max cut     : %f\n", fJetTPCRatioMaxCut) ;
2747   printf("Jet TPC Ratio min cut     : %f\n", fJetTPCRatioMinCut) ;
2748   printf("Fast recons       : %d\n", fOptFast);
2749   printf("Inv Mass max cut  : %f\n", fInvMassMaxCut) ;
2750   printf("Inv Mass min cut  : %f\n", fInvMassMinCut) ;
2751
2752 }
2753
2754 //__________________________________________________________________________
2755 TChain * AliPHOSGammaJet::ReadESDfromdisk(const UInt_t eventsToRead,
2756                          const TString dirName,
2757                          const TString esdTreeName,
2758                          const char *  pattern)
2759 {
2760   // Reads ESDs from Disk
2761  TChain *  rv = 0;
2762
2763   AliInfo( Form("\nReading files in %s \nESD tree name is %s \nReading %d events",
2764                 dirName.Data(), esdTreeName.Data(), eventsToRead) ) ;
2765
2766   // create a TChain of all the files
2767   TChain * cESDTree = new TChain(esdTreeName) ;
2768
2769   // read from the directory file until the require number of events are collected
2770   void * from = gSystem->OpenDirectory(dirName) ;
2771   if (!from) {
2772     AliError( Form("Directory %s does not exist") ) ;
2773     rv = 0 ;
2774   }
2775   else{ // reading file names from directory
2776     const char * subdir ;
2777     // search all subdirectories witch matching pattern
2778     while( (subdir = gSystem->GetDirEntry(from))  &&
2779            (cESDTree->GetEntries() < eventsToRead)) {
2780       if ( strstr(subdir, pattern) != 0 ) {
2781         char file[200] ;
2782         sprintf(file, "%s%s/AliESDs.root", dirName.Data(), subdir);
2783         AliInfo( Form("Adding %s\n", file) );
2784         cESDTree->Add(file) ;
2785       }
2786     } // while file
2787
2788     AliInfo( Form(" %d events read", cESDTree->GetEntriesFast()) ) ;
2789     rv = cESDTree ;
2790
2791   } // reading file names from directory
2792   return rv ;
2793 }
2794
2795 //__________________________________________________________________________
2796 TChain * AliPHOSGammaJet::ReadESD(const UInt_t eventsToRead,
2797                   const TString dirName,
2798                   const TString esdTreeName,
2799                   const char *  pattern)
2800 {
2801   // Read AliESDs files and return a Chain of events
2802
2803   if ( dirName == "" ) {
2804     AliError("Give the name of the DIR where to find files") ;
2805     return 0 ;
2806   }
2807   if ( esdTreeName == "" )
2808     return ReadESDfromdisk(eventsToRead, dirName) ;
2809   else if ( strcmp(pattern, "") == 0 )
2810     return ReadESDfromdisk(eventsToRead, dirName, esdTreeName) ;
2811   else
2812     return ReadESDfromdisk(eventsToRead, dirName, esdTreeName, pattern) ;
2813 }
2814
2815 //___________________________________________________________________
2816 void AliPHOSGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone,
2817                              Double_t eta, Double_t phi)
2818 {
2819
2820   //Check if the particle is inside the cone defined by the leading particle
2821   b = kFALSE;
2822
2823   if(phi > TMath::TwoPi())
2824     phi-=TMath::TwoPi();
2825   if(phi < 0.)
2826     phi+=TMath::TwoPi();
2827
2828   Double_t  rad = 10000 + cone;
2829
2830   if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
2831     rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2832                       TMath::Power(part->Phi()-phi,2));
2833   else{
2834     if(part->Phi()-phi > TMath::TwoPi() - cone)
2835       rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2836                         TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
2837     if(part->Phi()-phi < -(TMath::TwoPi() - cone))
2838       rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
2839                         TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
2840   }
2841
2842   if(rad < cone )
2843     b = kTRUE;
2844
2845 }
2846
2847 //____________________________________________________________________________
2848 Double_t AliPHOSGammaJet::SigmaE(Double_t energy)
2849 {
2850   // Calculates the energy dependent energy resolution
2851
2852   Double_t rv = -1 ;
2853
2854   rv = TMath::Sqrt( TMath::Power(fResPara1/energy, 2)
2855                + TMath::Power(fResPara2/TMath::Sqrt(energy), 2)
2856                + TMath::Power(fResPara3, 2) ) ;
2857
2858   return rv * energy ;
2859 }
2860
2861 //____________________________________________________________________________
2862 Double_t AliPHOSGammaJet::SigmaP(Double_t energy)
2863 {
2864   // Calculates the energy dependent position resolution
2865
2866   Double_t sigma = TMath::Sqrt(TMath::Power(fPosParaA,2) +
2867                                TMath::Power(fPosParaB,2) / energy) ;
2868
2869
2870   return sigma   ; // in cm
2871 }
2872