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