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