PHOS/AliPHOSGammaJet.cxx

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