PHOS/AliPHOSGammaJet.cxx

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