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