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