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