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