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