PHOS/AliPHOSGammaJet.cxx

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