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