New Gamma package
[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaJet.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 //_________________________________________________________________________
24 // Class for the analysis of gamma-jet correlations 
25 //  Basically it seaches for a prompt photon in the Calorimeters acceptance, 
26 //  if so we construct a jet around the highest pt particle in the opposite 
27 //  side in azimuth, inside the Central Tracking System (ITS+TPC) and 
28 //  EMCAL acceptances (only when PHOS detects the gamma). First the leading 
29 //  particle and then the jet have to fullfill several conditions 
30 //  (energy, direction ..) to be accepted. Then the fragmentation function 
31 //  of this jet is constructed   
32 //  Class created from old AliPHOSGammaJet 
33 //  (see AliRoot versions previous Release 4-09)
34 //
35 //*-- Author: Gustavo Conesa (LNF-INFN) 
36 //////////////////////////////////////////////////////////////////////////////
37
38
39 // --- ROOT system ---
40
41 #include <TFile.h>
42 #include <TParticle.h>
43 #include <TH2.h>
44
45 #include "AliAnaGammaJet.h" 
46 #include "AliESD.h"
47 #include "AliESDtrack.h"
48 #include "AliESDCaloCluster.h"
49 #include "Riostream.h"
50 #include "AliLog.h"
51
52 ClassImp(AliAnaGammaJet)
53
54 //____________________________________________________________________________
55 AliAnaGammaJet::AliAnaGammaJet(const char *name) : 
56   AliAnaGammaDirect(name), 
57   fSeveralConeAndPtCuts(0),
58   fPbPb(kFALSE), 
59   fJetsOnlyInCTS(0),
60   fEtaEMCALCut(0.),fPhiMaxCut(0.),
61   fPhiMinCut(0.), 
62   fInvMassMaxCut(0.), fInvMassMinCut(0.),
63   fJetCTSRatioMaxCut(0.),
64   fJetCTSRatioMinCut(0.), fJetRatioMaxCut(0.),
65   fJetRatioMinCut(0.), fNCone(0),
66   fNPt(0), fCone(0), fPtThreshold(0),
67   fPtJetSelectionCut(0.0),
68   fAngleMaxParam(), fSelect(0)
69 {
70
71   // ctor
72   fAngleMaxParam.Set(4) ;
73   fAngleMaxParam.Reset(0.);
74
75   for(Int_t i = 0; i<10; i++){
76     fCones[i]         = 0.0 ;
77     fNameCones[i]     = ""  ;
78     fPtThres[i]      = 0.0 ;
79     fNamePtThres[i]  = ""  ;
80     if( i < 6 ){
81       fJetXMin1[i]     = 0.0 ;
82       fJetXMin2[i]     = 0.0 ;
83       fJetXMax1[i]     = 0.0 ;
84       fJetXMax2[i]     = 0.0 ;
85       fBkgMean[i]      = 0.0 ;
86       fBkgRMS[i]       = 0.0 ;
87       if( i < 2 ){
88         fJetE1[i]        = 0.0 ;
89         fJetE2[i]        = 0.0 ;
90         fJetSigma1[i]    = 0.0 ;
91         fJetSigma2[i]    = 0.0 ;
92         fPhiEMCALCut[i]  = 0.0 ;
93       }
94     }
95   }
96         
97   TList * list = gDirectory->GetListOfKeys() ; 
98   TIter next(list) ; 
99   TH2F * h = 0 ;
100   Int_t index ; 
101   for (index = 0 ; index < list->GetSize()-1 ; index++) { 
102     //-1 to avoid GammaJet Task
103     h = dynamic_cast<TH2F*>(gDirectory->Get(list->At(index)->GetName())) ; 
104     fOutputContainer->Add(h) ; 
105   }
106
107   // Input slot #0 works with an Ntuple
108   DefineInput(0, TChain::Class());
109   // Output slot #0 writes into a TH1 container
110   DefineOutput(0,  TObjArray::Class()) ; 
111   
112 }
113
114
115 //____________________________________________________________________________
116 AliAnaGammaJet::AliAnaGammaJet(const AliAnaGammaJet & gj) : 
117   AliAnaGammaDirect(gj), 
118   fSeveralConeAndPtCuts(gj.fSeveralConeAndPtCuts), 
119   fPbPb(gj.fPbPb), fJetsOnlyInCTS(gj.fJetsOnlyInCTS),
120   fEtaEMCALCut(gj.fEtaEMCALCut),
121   fPhiMaxCut(gj.fPhiMaxCut), fPhiMinCut(gj.fPhiMinCut), 
122   fInvMassMaxCut(gj.fInvMassMaxCut), fInvMassMinCut(gj.fInvMassMinCut),
123   fRatioMinCut(gj.fRatioMinCut), 
124   fJetCTSRatioMaxCut(gj.fJetCTSRatioMaxCut),
125   fJetCTSRatioMinCut(gj.fJetCTSRatioMinCut), fJetRatioMaxCut(gj.fJetRatioMaxCut),
126   fJetRatioMinCut(gj.fJetRatioMinCut),  fNCone(gj.fNCone),
127   fNPt(gj.fNPt), fCone(gj.fCone), fPtThreshold(gj.fPtThreshold),
128   fPtJetSelectionCut(gj.fPtJetSelectionCut),
129   fOutputContainer(0), fAngleMaxParam(gj.fAngleMaxParam), 
130   fSelect(gj.fSelect)
131 {
132   // cpy ctor
133   SetName (gj.GetName()) ; 
134   SetTitle(gj.GetTitle()) ; 
135
136   for(Int_t i = 0; i<10; i++){
137     fCones[i]        = gj.fCones[i] ;
138     fNameCones[i]    = gj.fNameCones[i] ;
139     fPtThres[i]      = gj.fPtThres[i] ;
140     fNamePtThres[i]  = gj.fNamePtThres[i] ;
141     if( i < 6 ){
142       fJetXMin1[i]       = gj.fJetXMin1[i] ;
143       fJetXMin2[i]       = gj.fJetXMin2[i] ;
144       fJetXMax1[i]       = gj.fJetXMax1[i] ;
145       fJetXMax2[i]       = gj.fJetXMax2[i] ;
146       fBkgMean[i]        = gj.fBkgMean[i] ;
147       fBkgRMS[i]         = gj.fBkgRMS[i] ;
148       if( i < 2 ){
149         fJetE1[i]        = gj.fJetE1[i] ;
150         fJetE2[i]        = gj.fJetE2[i] ;
151         fJetSigma1[i]    = gj.fJetSigma1[i] ;
152         fJetSigma2[i]    = gj.fJetSigma2[i] ;
153         fPhiEMCALCut[i]  = gj.fPhiEMCALCut[i] ;
154       }
155     }          
156   } 
157 }
158
159 //____________________________________________________________________________
160 AliAnaGammaJet::~AliAnaGammaJet() 
161 {
162   // Remove all pointers
163   fOutputContainer->Clear() ; 
164   delete fOutputContainer ;
165   
166   delete fhChargeRatio  ; 
167   delete fhPi0Ratio   ; 
168   delete fhDeltaPhiCharge  ;  
169   delete fhDeltaPhiPi0   ; 
170   delete fhDeltaEtaCharge  ; 
171   delete fhDeltaEtaPi0  ; 
172   delete fhAnglePair  ; 
173   delete fhAnglePairAccepted  ; 
174   delete fhAnglePairNoCut  ; 
175   delete fhAnglePairLeadingCut  ; 
176   delete fhAnglePairAngleCut   ; 
177   delete fhAnglePairAllCut   ; 
178   delete fhAnglePairLeading  ; 
179   delete fhInvMassPairNoCut    ; 
180   delete fhInvMassPairLeadingCut  ; 
181   delete fhInvMassPairAngleCut  ; 
182   delete fhInvMassPairAllCut   ; 
183   delete fhInvMassPairLeading  ; 
184   delete fhNBkg   ; 
185   delete fhNLeading  ; 
186   delete fhNJet  ; 
187   delete fhJetRatio  ; 
188   delete fhJetPt   ; 
189   delete fhBkgRatio   ; 
190   delete fhBkgPt  ; 
191   delete fhJetFragment  ; 
192   delete fhBkgFragment  ; 
193   delete fhJetPtDist  ; 
194   delete fhBkgPtDist  ; 
195   
196   delete [] fhJetRatios;  
197   delete [] fhJetPts;  
198   delete [] fhBkgRatios;
199   delete [] fhBkgPts;  
200
201   delete [] fhNLeadings;
202   delete [] fhNJets;  
203   delete [] fhNBkgs;
204   
205   delete [] fhJetFragments;
206   delete [] fhBkgFragments;
207   delete [] fhJetPtDists;
208   delete [] fhBkgPtDists;
209   
210 }
211
212 //____________________________________________________________________________
213 Double_t AliAnaGammaJet::CalculateJetRatioLimit(const Double_t ptg, 
214                                                  const Double_t *par, 
215                                                  const Double_t *x) {
216
217   //Parametrized cut for the energy of the jet.
218
219   Double_t epp = par[0] + par[1] * ptg ;
220   Double_t spp = par[2] + par[3] * ptg ;
221   Double_t f   = x[0]   + x[1]   * ptg ;
222   Double_t epb = epp + par[4] ;
223   Double_t spb = TMath::Sqrt(spp*spp+ par[5]*par[5]) ;
224   Double_t rat = (epb - spb * f) / ptg ;
225   //Info("CalculateLimit","epp %f, spp %f, f %f", epp, spp, f);
226   //Info("CalculateLimit","epb %f, spb %f, rat %f", epb, spb, rat);
227   return rat ;
228 }
229
230
231 //____________________________________________________________________________
232 void AliAnaGammaJet::Exec(Option_t *) 
233 {
234   
235   // Processing of one event
236   
237   //Get ESDs
238   Long64_t entry = GetChain()->GetReadEntry() ;
239   
240   if (!GetESD()) {
241     AliError("fESD is not connected to the input!") ; 
242     return ; 
243   }
244   
245   if (GetPrintInfo()) 
246     AliInfo(Form("%s ----> Processing event # %lld",  (dynamic_cast<TChain *>(GetChain()))->GetFile()->GetName(), entry)) ; 
247   
248   //CreateTLists with arrays of TParticles. Filled with particles only relevant for the analysis.
249   
250   TClonesArray * particleList = new TClonesArray("TParticle",1000); // All particles refitted in CTS and detected in EMCAL (jet)
251   TClonesArray * plCTS         = new TClonesArray("TParticle",1000); // All particles refitted in Central Tracking System (ITS+TPC)
252   TClonesArray * plNe         = new TClonesArray("TParticle",1000);   // All particles measured in Jet Calorimeter
253   TClonesArray * plCalo     = new TClonesArray("TParticle",1000);  // All particles measured in Prompt Gamma calorimeter 
254   
255   
256   TParticle *pGamma = new TParticle(); //It will contain the kinematics of the found prompt gamma
257   TParticle *pLeading = new TParticle(); //It will contain the kinematics of the found leading particle
258   
259   Bool_t iIsInPHOSorEMCAL = kFALSE ; //To check if Gamma was in any calorimeter
260   
261   //Fill lists with photons, neutral particles and charged particles
262   //look for the highest energy photon in the event inside fCalorimeter
263   AliDebug(2, "Fill particle lists, get prompt gamma");
264   
265   //Fill particle lists 
266   
267   if(GetCalorimeter() == "PHOS")
268     CreateParticleList(particleList, plCTS,plNe,plCalo); 
269   else if(GetCalorimeter() == "EMCAL")
270     CreateParticleList(particleList, plCTS,plCalo,plNe); 
271   else
272     AliError("No calorimeter selected");
273   
274   //Search highest energy prompt gamma in calorimeter
275   GetPromptGamma(plCalo,  plCTS, pGamma, iIsInPHOSorEMCAL) ; 
276   
277   AliDebug(1, Form("Is Gamma in %s? %d",GetCalorimeter().Data(),iIsInPHOSorEMCAL));
278   AliDebug(3,Form("Charged list entries %d, Neutral list entries %d, %s list entries %d",
279                   plCTS->GetEntries(),plNe->GetEntries(), GetCalorimeter().Data(), plCalo->GetEntries()));
280
281   //If there is any prompt photon in fCalorimeter, 
282   //search jet leading particle
283   
284   if(iIsInPHOSorEMCAL){
285     if (GetPrintInfo())
286       AliInfo(Form("Prompt Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
287  
288     AliDebug(2, "Get Leading Particles, Make jet");
289
290     //Search leading particles in CTS and EMCAL 
291     if(GetLeadingParticle(plCTS, plNe, pGamma, pLeading)){
292       if (GetPrintInfo())
293         AliInfo(Form("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ;
294
295       //Search Jet
296       if(!fSeveralConeAndPtCuts)
297         MakeJet(particleList,pGamma,pLeading,"");
298       else{
299         for(Int_t icone = 0; icone<fNCone; icone++) {
300           for(Int_t ipt = 0; ipt<fNPt;ipt++) {  
301             TString lastname ="Cone"+ fNameCones[icone]+"Pt"+ fNamePtThres[ipt];
302             MakeJet(particleList, pGamma, pLeading,lastname);
303           }//icone
304         }//ipt
305       }//fSeveralConeAndPtCuts
306     }//Leading
307   }//Gamma in Calo
308      
309   AliDebug(2, "End of analysis, delete pointers");
310
311   particleList->Delete() ; 
312   plCTS->Delete() ;
313   plNe->Delete() ;
314   plCalo->Delete() ;
315   pLeading->Delete();
316   pGamma->Delete();
317
318   delete plNe ;
319   delete plCalo ;
320   delete plCTS ;
321   delete particleList ;
322   //  delete pLeading;
323   //  delete pGamma;
324
325   PostData(0, fOutputContainer);
326 }    
327
328 //____________________________________________________________________________
329 void AliAnaGammaJet::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname)
330 {
331   //Fill histograms wth jet fragmentation 
332   //and number of selected jets and leading particles
333   //and the background multiplicity
334   TParticle * particle = 0 ;
335   Int_t ipr = 0;
336   Float_t  charge = 0;
337
338   TIter next(pl) ; 
339   while ( (particle = dynamic_cast<TParticle*>(next())) ) {
340     ipr++ ;
341     Double_t pt = particle->Pt();
342
343     charge = TDatabasePDG::Instance()
344       ->GetParticle(particle->GetPdgCode())->Charge();
345     if(charge != 0){//Only jet Charged particles 
346       dynamic_cast<TH2F*>
347         (fOutputContainer->FindObject(type+"Fragment"+lastname))
348         ->Fill(ptg,pt/ptg);
349       dynamic_cast<TH2F*>
350         (fOutputContainer->FindObject(type+"PtDist"+lastname))
351         ->Fill(ptg,pt);
352     }//charged
353
354   }//while
355
356   if(type == "Bkg")
357     dynamic_cast<TH1F*>
358       (fOutputContainer->FindObject("NBkg"+lastname))
359       ->Fill(ipr);
360   else{
361     dynamic_cast<TH1F*>
362       (fOutputContainer->FindObject("NJet"+lastname))->
363       Fill(ptg);
364     dynamic_cast<TH2F*>
365       (fOutputContainer->FindObject("NLeading"+lastname))
366       ->Fill(ptg,ptl);
367   }
368   
369 }
370
371 //____________________________________________________________________________
372 void  AliAnaGammaJet::GetLeadingCharge(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading) const 
373 {  
374   //Search for the charged particle with highest with 
375   //Phi=Phi_gamma-Pi and pT=0.1E_gamma 
376   Double_t pt  = -100.;
377   Double_t phi = -100.;
378
379   for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
380
381     TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
382
383     Double_t ptl  = particle->Pt();
384     Double_t rat  = ptl/pGamma->Pt() ;
385     Double_t phil = particle->Phi() ;
386     Double_t phig = pGamma->Phi();
387
388     //Selection within angular and energy limits
389     if(((phig-phil)> fPhiMinCut) && ((phig-phil)<fPhiMaxCut) &&
390         (rat > fRatioMinCut) && (rat < fRatioMaxCut)  && (ptl  > pt)) {
391        phi = phil ;
392        pt  = ptl ;
393        pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
394        AliDebug(4,Form("Charge in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, particle->Eta(), phi,rat)) ;
395      }
396   }
397   
398   AliDebug(3,Form("Leading in CTS: pt %f eta %f phi %f pt/Eg %f \n", pt, pLeading->Eta(), phi,pt/pGamma->Pt())) ;
399
400 }
401
402
403 //____________________________________________________________________________
404 void  AliAnaGammaJet::GetLeadingPi0(TClonesArray * pl, TParticle * pGamma, TParticle * pLeading)  
405 {  
406
407   //Search for the neutral pion with highest with 
408   //Phi=Phi_gamma-Pi and pT=0.1E_gamma 
409   Double_t pt  = -100.;
410   Double_t phi = -100.;
411   Double_t ptl = -100.;
412   Double_t rat = -100.; 
413   Double_t phil = -100. ;
414   Double_t ptg  = pGamma->Pt();
415   Double_t phig = pGamma->Phi();
416
417   TIter next(pl);
418   TParticle * particle = 0;
419   
420   Int_t iPrimary = -1;
421   TLorentzVector gammai,gammaj;
422   Double_t angle = 0., e = 0., invmass = 0.;
423   Double_t anglef = 0., ef = 0., invmassf = 0.;
424   Int_t ksPdg = 0;
425   Int_t jPrimary=-1;
426
427   while ( (particle = (TParticle*)next()) ) {
428     iPrimary++;   
429     
430     ksPdg = particle->GetPdgCode();
431     ptl  = particle->Pt();
432     if(ksPdg == 111){ //2 gamma overlapped, found with PID
433       rat = ptl/ptg ;
434       phil = particle->Phi() ;
435       //Selection within angular and energy limits
436       if((ptl> pt)&& (rat > fRatioMinCut) && (rat < fRatioMaxCut) && 
437          ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
438         phi = phil ;
439         pt  = ptl ;
440         pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
441         AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
442       }// cuts
443     }// pdg = 111
444     else if(ksPdg == 22){//1 gamma
445       //Search the photon companion in case it comes from  a Pi0 decay
446       //Apply several cuts to select the good pair
447       particle->Momentum(gammai);
448       jPrimary=-1;
449       TIter next2(pl);
450       while ( (particle = (TParticle*)next2()) ) {
451         jPrimary++;
452         if(jPrimary>iPrimary){
453           ksPdg = particle->GetPdgCode();
454
455           if(ksPdg == 22){
456             particle->Momentum(gammaj);
457             //Info("GetLeadingPi0","Egammai %f, Egammaj %f", 
458             //gammai.Pt(),gammaj.Pt());
459             ptl  = (gammai+gammaj).Pt();
460             phil = (gammai+gammaj).Phi();
461             if(phil < 0)
462               phil+=TMath::TwoPi();
463             rat = ptl/ptg ;
464             invmass = (gammai+gammaj).M();
465             angle = gammaj.Angle(gammai.Vect());
466             //Info("GetLeadingPi0","Angle %f", angle);
467             e = (gammai+gammaj).E();
468             //Fill histograms with no cuts applied.
469             fhAnglePairNoCut->Fill(e,angle);
470             fhInvMassPairNoCut->Fill(ptg,invmass);
471             //First cut on the energy and azimuth of the pair
472             if((rat > fRatioMinCut) && (rat < fRatioMaxCut) && 
473                ((phig-phil)>fPhiMinCut)&&((phig-phil)<fPhiMaxCut)){
474               
475               fhAnglePairLeadingCut->Fill(e,angle);
476               fhInvMassPairLeadingCut->Fill(ptg,invmass);
477               //Second cut on the aperture of the pair
478               if(IsAngleInWindow(angle,e)){
479                 fhAnglePairAngleCut->Fill(e,angle);
480                 fhInvMassPairAngleCut->Fill(ptg,invmass);
481                 
482                 //Info("GetLeadingPi0","InvMass %f", invmass);
483                 //Third cut on the invariant mass of the pair
484                 if((invmass>fInvMassMinCut) && (invmass<fInvMassMaxCut)){ 
485                   fhInvMassPairAllCut->Fill(ptg,invmass);
486                   fhAnglePairAllCut->Fill(e,angle);
487                   if(ptl > pt ){
488                     pt       = ptl;
489                     phi      = phil ;
490                     ef       = e ;
491                     anglef   = angle ;
492                     invmassf = invmass ;
493                     pLeading->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
494                     AliDebug(4,Form("Pi0 candidate: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
495                   }
496                 }//cuts
497               }//(invmass>0.125) && (invmass<0.145)
498             }//gammaj.Angle(gammai.Vect())<0.04
499           }//if pdg = 22
500         }//iprimary<jprimary
501       }//while
502     }// if pdg = 22
503     //     }
504   }//while
505   
506   if(ef > 0.){//Final pi0 found, highest pair energy, fill histograms
507     fhInvMassPairLeading->Fill(ptg,invmassf);
508     fhAnglePairLeading->Fill(ef,anglef);
509   }
510   
511   AliDebug(3,Form("Leading EMCAL: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/pGamma->Pt())) ;
512 }
513
514 //____________________________________________________________________________
515 Bool_t  AliAnaGammaJet::GetLeadingParticle(TClonesArray * plCTS, TClonesArray * plNe,  
516                                          TParticle * pGamma, TParticle * pLeading) 
517 {
518   //Search Charged or Neutral leading particle, select the highest one.
519   
520   TParticle * pLeadingCh = new TParticle();
521   TParticle * pLeadingPi0 = new TParticle();
522   
523   Double_t ptg  =  pGamma->Pt(); 
524   Double_t phig = pGamma->Phi(); 
525   Double_t etag = pGamma->Eta(); 
526   
527   if(GetCalorimeter() == "PHOS" && !fJetsOnlyInCTS)
528     {
529       AliDebug(3, "GetLeadingPi0");
530       GetLeadingPi0   (plNe, pGamma, pLeadingPi0) ;
531       AliDebug(3, "GetLeadingCharge");
532       GetLeadingCharge(plCTS, pGamma, pLeadingCh) ;
533       
534       Double_t ptch = pLeadingCh->Pt(); 
535       Double_t phich = pLeadingCh->Phi(); 
536       Double_t etach = pLeadingCh->Eta(); 
537       Double_t ptpi = pLeadingPi0->Pt(); 
538       Double_t phipi = pLeadingPi0->Phi(); 
539       Double_t etapi = pLeadingPi0->Eta(); 
540
541       //Is leading cone inside EMCAL acceptance?
542       
543       Bool_t insidech = kFALSE ;
544       if((phich - fCone) >  fPhiEMCALCut[0] && 
545          (phich + fCone) <  fPhiEMCALCut[1] && 
546         (etach-fCone) < fEtaEMCALCut )
547         insidech = kTRUE ;
548       
549       Bool_t insidepi = kFALSE ;
550       if((phipi - fCone) >  fPhiEMCALCut[0] && 
551          (phipi + fCone) <  fPhiEMCALCut[1] &&
552         (etapi-fCone) < fEtaEMCALCut )
553         insidepi = kTRUE ;
554
555       AliDebug(2,Form("Leading:  charged pt %f, pi0 pt  %f",ptch,ptpi)) ;
556       
557       if (ptch > 0 || ptpi > 0)
558         {
559           if(insidech && (ptch > ptpi))
560             {
561               if (GetPrintInfo())
562                 AliInfo("Leading found in CTS");
563               pLeading->SetMomentum(pLeadingCh->Px(),pLeadingCh->Py(),pLeadingCh->Pz(),pLeadingCh->Energy());
564               AliDebug(3,Form("Final leading found in CTS, pt %f, phi %f, eta %f",ptch,phich,etach)) ;
565               fhChargeRatio->Fill(ptg,ptch/pGamma->Pt());
566               fhDeltaPhiCharge->Fill(ptg,pGamma->Phi()-phich);
567               fhDeltaEtaCharge->Fill(ptg,pGamma->Eta()-etach);
568               return 1;
569             }
570           
571           else if((ptpi > ptch) && insidepi)
572             {
573               if (GetPrintInfo())
574                 AliInfo("Leading found in EMCAL");
575               pLeading->SetMomentum(pLeadingPi0->Px(),pLeadingPi0->Py(),pLeadingPi0->Pz(),pLeadingPi0->Energy());
576               AliDebug(3,Form("Final leading found in EMCAL, pt %f, phi %f, eta %f",ptpi,phipi,etapi)) ;
577               fhPi0Ratio     ->Fill(ptg,ptpi/ptg);
578               fhDeltaPhiPi0->Fill(ptg,phig-phipi);
579               fhDeltaEtaPi0->Fill(ptg,etag-etapi);
580               return 1;     
581             }
582
583           else{
584             AliDebug(3,"NO LEADING PARTICLE FOUND");}
585           return 0; 
586         }
587       else{
588         AliDebug(3,"NO LEADING PARTICLE FOUND");
589         return 0;
590       }
591     }
592
593   else
594     {
595       //No calorimeter present for Leading particle detection
596       GetLeadingCharge(plCTS, pGamma, pLeading) ;
597       Double_t ptch = pLeading->Pt(); 
598       Double_t phich = pLeading->Phi(); 
599       Double_t etach = pLeading->Eta(); 
600       if(ptch > 0){
601         fhChargeRatio->Fill(ptg,ptch/ptg);
602         fhDeltaPhiCharge->Fill(ptg,phig-phich);
603         fhDeltaEtaCharge->Fill(ptg,etag-etach);
604         AliDebug(3,Form("Leading found :  pt %f, phi %f, eta %f",ptch,phich,etach)) ;
605         return 1;
606       }
607       else
608         {
609           AliDebug(3,"NO LEADING PARTICLE FOUND");      
610           return 0;
611         }
612     }
613 }
614
615   //____________________________________________________________________________
616 void AliAnaGammaJet::Init(const Option_t * )
617 {
618 //   // Initialisation of branch container 
619   AliAnaGammaDirect::Init();
620  
621   //Initialize the parameters of the analysis.
622   //fCalorimeter="PHOS";
623   fAngleMaxParam.Set(4) ;
624   fAngleMaxParam.AddAt(0.4,0);//={0.4,-0.25,0.025,-2e-4};
625   fAngleMaxParam.AddAt(-0.25,1) ;
626   fAngleMaxParam.AddAt(0.025,2) ;
627   fAngleMaxParam.AddAt(-2e-4,3) ;
628   fSeveralConeAndPtCuts   = kFALSE ;
629   //fPrintInfo           = kTRUE;
630   fPbPb                = kFALSE ;
631   fInvMassMaxCut  = 0.16 ;
632   fInvMassMinCut  = 0.11 ;
633   //fJetsOnlyInCTS    = kTRUE ;
634   fEtaEMCALCut     = 0.7 ;
635   fPhiEMCALCut[0] = 80. *TMath::Pi()/180.;
636   fPhiEMCALCut[1] = 190.*TMath::Pi()/180.;
637   fPhiMaxCut      = 3.4 ;
638   fPhiMinCut      = 2.9 ;
639
640   //Jet selection parameters
641   //Fixed cut (old)
642   fRatioMaxCut    = 1.0 ;
643   fRatioMinCut    = 0.1 ; 
644   fJetRatioMaxCut = 1.2 ; 
645   fJetRatioMinCut = 0.3 ; 
646   fJetsOnlyInCTS = kFALSE ;
647   fJetCTSRatioMaxCut = 1.2 ;
648   fJetCTSRatioMinCut = 0.3 ;
649   fSelect         = 0  ;
650
651   //Cut depending on gamma energy
652
653   fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applied
654   //Reconstructed jet energy dependence parameters 
655   //e_jet = a1+e_gamma b2. 
656   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
657   fJetE1[0] = -5.75; fJetE1[1] = -4.1;
658   fJetE2[0] = 1.005; fJetE2[1] = 1.05;
659
660   //Reconstructed sigma of jet energy dependence parameters 
661   //s_jet = a1+e_gamma b2. 
662   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
663   fJetSigma1[0] = 2.65;   fJetSigma1[1] = 2.75;
664   fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
665
666   //Background mean energy and RMS
667   //Index 0-> No BKG; Index 1-> BKG > 2 GeV; 
668   //Index 2-> (low pt jets)BKG > 0.5 GeV;
669   //Index > 2, same for CTS conf
670   fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
671   fBkgMean[3] = 0.; fBkgMean[4] = 6.4;  fBkgMean[5] = 48.6;
672   fBkgRMS[0]  = 0.; fBkgRMS[1]  = 7.5;  fBkgRMS[2]  = 22.0; 
673   fBkgRMS[3]  = 0.; fBkgRMS[4]  = 5.4;  fBkgRMS[5]  = 13.2; 
674
675   //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
676   //limits for monoenergetic jets.
677   //Index 0-> No BKG; Index 1-> BKG > 2 GeV; 
678   //Index 2-> (low pt jets) BKG > 0.5 GeV;
679   //Index > 2, same for CTS conf
680
681   fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ; 
682   fJetXMin1[3] =-2.0  ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1  ;
683   fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034; 
684   fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
685   fJetXMax1[0] =-3.8  ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6  ; 
686   fJetXMax1[3] =-2.7  ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7  ;
687   fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016; 
688   fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
689
690
691   //Different cones and pt thresholds to construct the jet
692
693   fCone        = 0.3  ;
694   fPtThreshold = 0.   ;
695   fNCone       = 3    ;
696   fNPt         = 3    ;
697   fCones[1]    = 0.2  ; fNameCones[1]   = "02" ;
698   fPtThres[0]  = 0.0  ; fNamePtThres[0] = "00" ;
699   fCones[0]    = 0.3  ; fNameCones[0]   = "03" ;
700   fPtThres[1]  = 0.5  ; fNamePtThres[1] = "05" ;
701   fCones[2]    = 0.4  ; fNameCones[2]   = "04" ;
702   fPtThres[2]  = 1.0  ; fNamePtThres[2] = "10" ;
703   //Fill particle lists when PID is ok
704   // fEMCALPID = kFALSE;
705   // fPHOSPID = kFALSE;
706
707   //Initialization of histograms 
708
709   MakeHistos() ; 
710
711 }
712
713 //__________________________________________________________________________-
714 Bool_t AliAnaGammaJet::IsAngleInWindow(const Float_t angle,const Float_t e) {
715   //Check if the opening angle of the candidate pairs is inside 
716   //our selection windowd
717
718   Bool_t result = kFALSE;
719   Double_t mpi0 = 0.1349766;
720   Double_t max =  fAngleMaxParam.At(0)*TMath::Exp(fAngleMaxParam.At(1)*e)
721     +fAngleMaxParam.At(2)+fAngleMaxParam.At(3)*e;
722   Double_t arg = (e*e-2*mpi0*mpi0)/(e*e);
723   Double_t min = 100. ;
724   if(arg>0.)
725     min = TMath::ACos(arg);
726
727   if((angle<max)&&(angle>=min))
728     result = kTRUE;
729  
730   return result;
731 }
732
733 //__________________________________________________________________________-
734 Bool_t AliAnaGammaJet::IsJetSelected(const Double_t ptg, const Double_t ptj){
735   //Check if the energy of the reconstructed jet is within an energy window
736
737   Double_t par[6];
738   Double_t xmax[2];
739   Double_t xmin[2];
740
741   Int_t iCTS = 0;
742   if(fJetsOnlyInCTS)
743     iCTS = 3 ;
744
745   if(!fPbPb){
746     //Phythia alone, jets with pt_th > 0.2, r = 0.3 
747     par[0] = fJetE1[0]; par[1] = fJetE2[0]; 
748     //Energy of the jet peak
749     //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
750     par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
751     //Sigma  of the jet peak
752     //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
753     par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
754     //Parameters reserved for PbPb bkg.
755     xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
756     xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
757     //Factor that multiplies sigma to obtain the best limits, 
758     //by observation, of mono jet ratios (ptjet/ptg)
759     //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
760    
761   }
762   else{
763     if(ptg > fPtJetSelectionCut){
764       //Phythia +PbPb with  pt_th > 2 GeV/c, r = 0.3 
765       par[0] = fJetE1[0]; par[1] = fJetE2[0]; 
766       //Energy of the jet peak, same as in pp
767       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
768       par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
769       //Sigma  of the jet peak, same as in pp
770       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
771       par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
772       //Mean value and RMS of PbPb Bkg 
773       xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
774       xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
775       //Factor that multiplies sigma to obtain the best limits, 
776       //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, 
777       //pt_th > 2 GeV, r = 0.3
778       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
779      
780     }
781     else{
782       //Phythia + PbPb with  pt_th > 0.5 GeV/c, r = 0.3
783       par[0] = fJetE1[1]; par[1] = fJetE2[1]; 
784       //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3 
785       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
786       par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
787       //Sigma  of the jet peak, pt_th > 2 GeV/c, r = 0.3
788       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
789       par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
790       //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
791       xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
792       xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
793       //Factor that multiplies sigma to obtain the best limits, 
794       //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, 
795       //pt_th > 2 GeV, r = 0.3
796       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
797      
798     }//If low pt jet in bkg
799   }//if Bkg
800
801  //Calculate minimum and maximum limits of the jet ratio.
802   Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
803   Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
804   
805   AliDebug(3,Form("Jet selection?  : Limits min %f, max %f,  pt_jet %f,  pt_gamma %f, pt_jet / pt_gamma %f",min,max,ptj,ptg,ptj/ptg));
806
807   if(( min < ptj/ptg ) && ( max > ptj/ptg))
808     return kTRUE;
809   else
810     return kFALSE;
811
812 }
813
814 //____________________________________________________________________________
815 void AliAnaGammaJet::MakeHistos()
816 {
817   // Create histograms to be saved in output file and 
818   // stores them in fOutputContainer
819   
820   fOutputContainer = new TObjArray(10000) ;
821  
822   TObjArray  * outputContainer =GetOutputContainer();
823   for(Int_t i = 0; i < outputContainer->GetEntries(); i++ )
824     fOutputContainer->Add(outputContainer->At(i)) ;
825
826   //
827   fhChargeRatio  = new TH2F
828     ("ChargeRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
829      120,0,120,120,0,1); 
830   fhChargeRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
831   fhChargeRatio->SetXTitle("p_{T #gamma} (GeV/c)");
832   fOutputContainer->Add(fhChargeRatio) ;
833   
834   fhDeltaPhiCharge  = new TH2F
835     ("DeltaPhiCharge","#phi_{#gamma} - #phi_{charge} vs p_{T #gamma}",
836      200,0,120,200,0,6.4); 
837   fhDeltaPhiCharge->SetYTitle("#Delta #phi");
838   fhDeltaPhiCharge->SetXTitle("p_{T #gamma} (GeV/c)");
839   fOutputContainer->Add(fhDeltaPhiCharge) ; 
840   
841   fhDeltaEtaCharge  = new TH2F
842     ("DeltaEtaCharge","#eta_{#gamma} - #eta_{charge} vs p_{T #gamma}",
843      200,0,120,200,-2,2); 
844   fhDeltaEtaCharge->SetYTitle("#Delta #eta");
845   fhDeltaEtaCharge->SetXTitle("p_{T #gamma} (GeV/c)");
846   fOutputContainer->Add(fhDeltaEtaCharge) ; 
847   
848   //
849   if(!fJetsOnlyInCTS){
850     fhPi0Ratio  = new TH2F
851       ("Pi0Ratio","p_{T leading  #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
852        120,0,120,120,0,1); 
853     fhPi0Ratio->SetYTitle("p_{T lead  #pi^{0}} /p_{T #gamma}");
854     fhPi0Ratio->SetXTitle("p_{T #gamma} (GeV/c)");
855     fOutputContainer->Add(fhPi0Ratio) ; 
856     
857     fhDeltaPhiPi0  = new TH2F
858       ("DeltaPhiPi0","#phi_{#gamma} - #phi_{ #pi^{0}} vs p_{T #gamma}",
859        200,0,120,200,0,6.4); 
860     fhDeltaPhiPi0->SetYTitle("#Delta #phi");
861     fhDeltaPhiPi0->SetXTitle("p_{T #gamma} (GeV/c)");
862     fOutputContainer->Add(fhDeltaPhiPi0) ; 
863     
864     fhDeltaEtaPi0  = new TH2F
865       ("DeltaEtaPi0","#eta_{#gamma} - #eta_{ #pi^{0}} vs p_{T #gamma}",
866        200,0,120,200,-2,2); 
867     fhDeltaEtaPi0->SetYTitle("#Delta #eta");
868     fhDeltaEtaPi0->SetXTitle("p_{T #gamma} (GeV/c)");
869     fOutputContainer->Add(fhDeltaEtaPi0) ; 
870  
871     //
872     fhAnglePair  = new TH2F
873       ("AnglePair",
874        "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}",
875        200,0,50,200,0,0.2); 
876     fhAnglePair->SetYTitle("Angle (rad)");
877     fhAnglePair->SetXTitle("E_{ #pi^{0}} (GeV/c)");
878     fOutputContainer->Add(fhAnglePair) ; 
879     
880     fhAnglePairAccepted  = new TH2F
881       ("AnglePairAccepted",
882        "Angle between #pi^{0} #gamma pair vs p_{T  #pi^{0}}, both #gamma in eta<0.7, inside window",
883        200,0,50,200,0,0.2); 
884     fhAnglePairAccepted->SetYTitle("Angle (rad)");
885     fhAnglePairAccepted->SetXTitle("E_{ #pi^{0}} (GeV/c)");
886     fOutputContainer->Add(fhAnglePairAccepted) ; 
887     
888     fhAnglePairNoCut  = new TH2F
889       ("AnglePairNoCut",
890        "Angle between all #gamma pair vs p_{T  #pi^{0}}",200,0,50,200,0,0.2); 
891     fhAnglePairNoCut->SetYTitle("Angle (rad)");
892     fhAnglePairNoCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
893     fOutputContainer->Add(fhAnglePairNoCut) ; 
894     
895     fhAnglePairLeadingCut  = new TH2F
896       ("AnglePairLeadingCut",
897        "Angle between all #gamma pair that have a good phi and pt vs p_{T  #pi^{0}}",
898        200,0,50,200,0,0.2); 
899     fhAnglePairLeadingCut->SetYTitle("Angle (rad)");
900     fhAnglePairLeadingCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
901     fOutputContainer->Add(fhAnglePairLeadingCut) ; 
902     
903     fhAnglePairAngleCut  = new TH2F
904       ("AnglePairAngleCut",
905        "Angle between all #gamma pair (angle + leading cut) vs p_{T  #pi^{0}}"
906        ,200,0,50,200,0,0.2); 
907     fhAnglePairAngleCut->SetYTitle("Angle (rad)");
908     fhAnglePairAngleCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
909     fOutputContainer->Add(fhAnglePairAngleCut) ;
910     
911     fhAnglePairAllCut  = new TH2F
912       ("AnglePairAllCut",
913        "Angle between all #gamma pair (angle + inv mass cut+leading) vs p_{T  #pi^{0}}"
914        ,200,0,50,200,0,0.2); 
915     fhAnglePairAllCut->SetYTitle("Angle (rad)");
916     fhAnglePairAllCut->SetXTitle("E_{ #pi^{0}} (GeV/c)");
917     fOutputContainer->Add(fhAnglePairAllCut) ; 
918     
919     fhAnglePairLeading  = new TH2F
920       ("AnglePairLeading",
921        "Angle between all #gamma pair finally selected vs p_{T  #pi^{0}}",
922        200,0,50,200,0,0.2); 
923     fhAnglePairLeading->SetYTitle("Angle (rad)");
924     fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
925     fOutputContainer->Add(fhAnglePairLeading) ; 
926     
927     //
928     fhInvMassPairNoCut  = new TH2F
929       ("InvMassPairNoCut","Invariant Mass of all #gamma pair vs p_{T #gamma}",
930        120,0,120,360,0,0.5); 
931     fhInvMassPairNoCut->SetYTitle("Invariant Mass (GeV/c^{2})");
932     fhInvMassPairNoCut->SetXTitle("p_{T #gamma} (GeV/c)");
933     fOutputContainer->Add(fhInvMassPairNoCut) ; 
934     
935     fhInvMassPairLeadingCut  = new TH2F
936       ("InvMassPairLeadingCut",
937        "Invariant Mass of #gamma pair (leading cuts) vs p_{T #gamma}",
938        120,0,120,360,0,0.5); 
939     fhInvMassPairLeadingCut->SetYTitle("Invariant Mass (GeV/c^{2})");
940     fhInvMassPairLeadingCut->SetXTitle("p_{T #gamma} (GeV/c)");
941     fOutputContainer->Add(fhInvMassPairLeadingCut) ; 
942     
943     fhInvMassPairAngleCut  = new TH2F
944       ("InvMassPairAngleCut",
945        "Invariant Mass of #gamma pair (angle cut) vs p_{T #gamma}",
946        120,0,120,360,0,0.5); 
947     fhInvMassPairAngleCut->SetYTitle("Invariant Mass (GeV/c^{2})");
948     fhInvMassPairAngleCut->SetXTitle("p_{T #gamma} (GeV/c)");
949     fOutputContainer->Add(fhInvMassPairAngleCut) ; 
950     
951     fhInvMassPairAllCut  = new TH2F
952       ("InvMassPairAllCut",
953        "Invariant Mass of #gamma pair (angle+invmass cut+leading) vs p_{T #gamma}",
954        120,0,120,360,0,0.5); 
955     fhInvMassPairAllCut->SetYTitle("Invariant Mass (GeV/c^{2})");
956     fhInvMassPairAllCut->SetXTitle("p_{T #gamma} (GeV/c)");
957     fOutputContainer->Add(fhInvMassPairAllCut) ; 
958     
959     fhInvMassPairLeading  = new TH2F
960       ("InvMassPairLeading",
961        "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
962        120,0,120,360,0,0.5); 
963     fhInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
964     fhInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
965     fOutputContainer->Add(fhInvMassPairLeading) ; 
966   }
967   
968   
969   if(!fSeveralConeAndPtCuts){
970     
971     //Count
972     fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000); 
973     fhNBkg->SetYTitle("counts");
974     fhNBkg->SetXTitle("N");
975     fOutputContainer->Add(fhNBkg) ; 
976     
977     fhNLeading  = new TH2F
978       ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120); 
979     fhNLeading->SetYTitle("p_{T charge} (GeV/c)");
980     fhNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
981     fOutputContainer->Add(fhNLeading) ; 
982     
983     fhNJet  = new TH1F("NJet","Accepted jets",240,0,120); 
984     fhNJet->SetYTitle("N");
985     fhNJet->SetXTitle("p_{T #gamma}(GeV/c)");
986     fOutputContainer->Add(fhNJet) ; 
987     
988     //Ratios and Pt dist of reconstructed (not selected) jets
989     //Jet
990     fhJetRatio  = new TH2F
991       ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
992        240,0,120,200,0,10);
993     fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
994     fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
995     fOutputContainer->Add(fhJetRatio) ; 
996     
997     fhJetPt  = new TH2F
998       ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
999     fhJetPt->SetYTitle("p_{T jet}");
1000     fhJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
1001     fOutputContainer->Add(fhJetPt) ; 
1002     
1003     //Bkg
1004     
1005     fhBkgRatio  = new TH2F
1006       ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
1007        240,0,120,200,0,10);
1008     fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
1009     fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
1010     fOutputContainer->Add(fhBkgRatio) ;
1011     
1012     fhBkgPt  = new TH2F
1013       ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
1014     fhBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
1015     fhBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
1016     fOutputContainer->Add(fhBkgPt) ;
1017     
1018     //Jet Distributions
1019     
1020     fhJetFragment  = 
1021       new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
1022                240,0.,120.,1000,0.,1.2); 
1023     fhJetFragment->SetYTitle("x_{T}");
1024     fhJetFragment->SetXTitle("p_{T #gamma}");
1025     fOutputContainer->Add(fhJetFragment) ;
1026     
1027     fhBkgFragment  = new TH2F
1028       ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
1029        240,0.,120.,1000,0.,1.2);
1030     fhBkgFragment->SetYTitle("x_{T}");
1031     fhBkgFragment->SetXTitle("p_{T #gamma}");
1032     fOutputContainer->Add(fhBkgFragment) ;
1033     
1034     fhJetPtDist  = 
1035       new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); 
1036     fhJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
1037     fOutputContainer->Add(fhJetPtDist) ;
1038     
1039     fhBkgPtDist  = new TH2F
1040       ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); 
1041     fhBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
1042     fOutputContainer->Add(fhBkgPtDist) ;
1043
1044   }
1045   else{
1046     //If we want to study the jet for different cones and pt
1047     
1048     for(Int_t icone = 0; icone<fNCone; icone++){
1049       for(Int_t ipt = 0; ipt<fNPt;ipt++){ 
1050         
1051         //Jet
1052         
1053         fhJetRatios[icone][ipt]  = new TH2F
1054           ("JetRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
1055            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
1056            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
1057            240,0,120,200,0,10);
1058         fhJetRatios[icone][ipt]->
1059           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1060         fhJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1061         fOutputContainer->Add(fhJetRatios[icone][ipt]) ; 
1062         
1063         
1064         fhJetPts[icone][ipt]  = new TH2F
1065           ("JetPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
1066            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
1067            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
1068            240,0,120,400,0,200);
1069         fhJetPts[icone][ipt]->
1070           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1071         fhJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1072         fOutputContainer->Add(fhJetPts[icone][ipt]) ; 
1073         
1074         //Bkg
1075         fhBkgRatios[icone][ipt]  = new TH2F
1076           ("BkgRatioCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
1077            "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
1078            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
1079            240,0,120,200,0,10);
1080         fhBkgRatios[icone][ipt]->
1081           SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
1082         fhBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1083         fOutputContainer->Add(fhBkgRatios[icone][ipt]) ; 
1084         
1085         fhBkgPts[icone][ipt]  = new TH2F
1086           ("BkgPtCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt], 
1087            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
1088            +fNameCones[icone]+", pt>" +fNamePtThres[ipt]+" GeV/c",
1089            240,0,120,400,0,200);
1090         fhBkgPts[icone][ipt]->
1091           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
1092         fhBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1093         fOutputContainer->Add(fhBkgPts[icone][ipt]) ; 
1094         
1095         //Counts
1096         fhNBkgs[icone][ipt]  = new TH1F
1097           ("NBkgCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1098            "bkg multiplicity cone ="+fNameCones[icone]+", pt>" 
1099            +fNamePtThres[ipt]+" GeV/c",9000,0,9000); 
1100         fhNBkgs[icone][ipt]->SetYTitle("counts");
1101         fhNBkgs[icone][ipt]->SetXTitle("N");
1102         fOutputContainer->Add(fhNBkgs[icone][ipt]) ; 
1103         
1104         fhNLeadings[icone][ipt]  = new TH2F
1105           ("NLeadingCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1106            "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fNameCones[icone]+", pt>" 
1107            +fNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120); 
1108         fhNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
1109         fhNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
1110         fOutputContainer->Add(fhNLeadings[icone][ipt]) ; 
1111         
1112         fhNJets[icone][ipt]  = new TH1F
1113           ("NJetCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1114            "Number of neutral jets, cone ="+fNameCones[icone]+", pt>" 
1115            +fNamePtThres[ipt]+" GeV/c",120,0,120); 
1116         fhNJets[icone][ipt]->SetYTitle("N");
1117         fhNJets[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
1118         fOutputContainer->Add(fhNJets[icone][ipt]) ; 
1119         
1120         //Fragmentation Function
1121         fhJetFragments[icone][ipt]  = new TH2F
1122           ("JetFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1123            "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" 
1124            +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); 
1125         fhJetFragments[icone][ipt]->SetYTitle("x_{T}");
1126         fhJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
1127         fOutputContainer->Add(fhJetFragments[icone][ipt]) ; 
1128         
1129         fhBkgFragments[icone][ipt]  = new TH2F
1130           ("BkgFragmentCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1131            "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fNameCones[icone]+", pt>" 
1132            +fNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); 
1133         fhBkgFragments[icone][ipt]->SetYTitle("x_{T}");
1134         fhBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
1135         fOutputContainer->Add(fhBkgFragments[icone][ipt]) ; 
1136         
1137         //Jet particle distribution
1138         
1139         fhJetPtDists[icone][ipt]  = new TH2F
1140           ("JetPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1141            "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
1142            " GeV/c",120,0.,120.,120,0.,120.); 
1143         fhJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1144         fOutputContainer->Add(fhJetPtDists[icone][ipt]) ; 
1145         
1146         fhBkgPtDists[icone][ipt]  = new TH2F
1147           ("BkgPtDistCone"+fNameCones[icone]+"Pt"+fNamePtThres[ipt],
1148            "p_{T i}, cone ="+fNameCones[icone]+", pt>" +fNamePtThres[ipt]+
1149            " GeV/c",120,0.,120.,120,0.,120.); 
1150         fhBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
1151         fOutputContainer->Add(fhBkgPtDists[icone][ipt]) ; 
1152         
1153       }//ipt
1154     } //icone
1155   }//If we want to study any cone or pt threshold
1156 }
1157
1158 //____________________________________________________________________________
1159 void AliAnaGammaJet::MakeJet(TClonesArray * pl, TParticle * pGamma, TParticle* pLeading,TString lastname)
1160 {
1161   //Fill the jet with the particles around the leading particle with 
1162   //R=fCone and pt_th = fPtThres. Caclulate the energy of the jet and 
1163   //check if we select it. Fill jet histograms
1164   
1165   TClonesArray * jetList = new TClonesArray("TParticle",1000);
1166   TClonesArray * bkgList = new TClonesArray("TParticle",1000);
1167
1168   TLorentzVector jet   (0,0,0,0);  
1169   TLorentzVector bkg(0,0,0,0);
1170   TLorentzVector lv (0,0,0,0);
1171
1172   Double_t ptjet = 0.0;
1173   Double_t ptbkg = 0.0;
1174   Int_t n0 = 0;
1175   Int_t n1 = 0;  
1176   Bool_t b1 = kFALSE;
1177   Bool_t b0 = kFALSE;
1178   
1179   Double_t ptg  = pGamma->Pt();
1180   Double_t phig = pGamma->Phi();
1181   Double_t ptl  = pLeading->Pt();
1182   Double_t phil = pLeading->Phi();
1183   Double_t etal = pLeading->Eta();
1184
1185   Float_t ptcut = 0. ;
1186   if(fPbPb){
1187     if(ptg > fPtJetSelectionCut)  ptcut = 2. ;
1188     else                          ptcut = 0.5;
1189   }
1190
1191   TIter next(pl) ; 
1192   TParticle * particle = 0 ; 
1193   while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1194     
1195     b0 = kFALSE;
1196     b1 = kFALSE;
1197
1198     //Particles in jet 
1199     SetJet(particle, b0, fCone, etal, phil) ;  
1200
1201     if(b0){
1202       new((*jetList)[n0++]) TParticle(*particle) ;
1203       particle->Momentum(lv);
1204       if(particle->Pt() > ptcut ){
1205         jet+=lv;
1206         ptjet+=particle->Pt();
1207       }
1208     }
1209
1210     //Background around (phi_gamma-pi, eta_leading)
1211     SetJet(particle, b1, fCone,etal, phig) ;
1212
1213     if(b1) { 
1214       new((*bkgList)[n1++]) TParticle(*particle) ;
1215       particle->Momentum(lv);
1216       if(particle->Pt() > ptcut ){
1217         bkg+=lv;
1218         ptbkg+=particle->Pt();    
1219       }  
1220     }
1221   }
1222   
1223   ptjet = jet.Pt();
1224   ptbkg = bkg.Pt();
1225
1226   if(ptjet > 0.) {
1227
1228     AliDebug(2,Form("Gamma   pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg));
1229     
1230     //Fill histograms
1231     
1232     Double_t ratjet   = ptjet/ptg ;
1233     Double_t ratbkg  = ptbkg/ptg ;
1234     
1235     dynamic_cast<TH2F*>
1236       (fOutputContainer->FindObject("JetRatio"+lastname))
1237       ->Fill(ptg,ratjet);        
1238     dynamic_cast<TH2F*>
1239       (fOutputContainer->FindObject("JetPt"+lastname))
1240       ->Fill(ptg,ptjet);
1241     
1242     dynamic_cast<TH2F*>
1243       (fOutputContainer->FindObject("BkgRatio"+lastname))
1244       ->Fill(ptg,ratbkg);
1245     
1246     dynamic_cast<TH2F*>
1247       (fOutputContainer->FindObject("BkgPt"+lastname))
1248       ->Fill(ptg,ptbkg);
1249
1250
1251     //Jet selection
1252     Bool_t kSelect = kFALSE;
1253     if(fSelect == 0)
1254       kSelect = kTRUE; //Accept all jets, no restriction
1255     else if(fSelect == 1){
1256       //Selection with parametrized cuts
1257       if(IsJetSelected(ptg,ptjet))   kSelect = kTRUE;
1258     }
1259     else if(fSelect == 2){
1260       //Simple selection
1261       if(!fJetsOnlyInCTS){
1262         if((ratjet <  fJetRatioMaxCut) && (ratjet > fJetRatioMinCut )) kSelect = kTRUE;
1263       }
1264       else{
1265         if((ratjet <  fJetCTSRatioMaxCut) && (ratjet > fJetCTSRatioMinCut )) kSelect = kTRUE;
1266       }
1267     }
1268     else
1269       AliError("Jet selection option larger than 2, DONT SELECT JETS");
1270     
1271     
1272     if(kSelect){
1273       if (GetPrintInfo())
1274         AliInfo(Form("Jet Selected: pt %f ", ptjet)) ;
1275       
1276       FillJetHistos(jetList, ptg, ptl,"Jet",lastname);
1277       FillJetHistos(bkgList, ptg, ptl, "Bkg",lastname);
1278     }
1279   } //ptjet > 0
1280   
1281   jetList ->Delete();
1282   bkgList ->Delete();
1283   
1284 }
1285
1286 //____________________________________________________________________________
1287 void AliAnaGammaJet::Print(const Option_t * opt) const
1288 {
1289
1290   //Print some relevant parameters set for the analysis
1291   if(! opt)
1292     return;
1293
1294   Info("Print", "%s %s", GetName(), GetTitle() ) ;
1295   if(!fJetsOnlyInCTS)
1296     printf("EMCAL Acceptance cut  : phi min %f, phi max %f, eta %f \n", fPhiEMCALCut[0],  fPhiEMCALCut[1], fEtaEMCALCut) ;
1297   
1298   printf("Phi_Leading                                <     %f\n", fPhiMaxCut) ; 
1299   printf("Phi_Leading                                >     %f\n", fPhiMinCut) ;
1300   printf("pT Leading / pT Gamma             <     %f\n", fRatioMaxCut) ; 
1301   printf("pT Leading / pT Gamma             >     %f\n", fRatioMinCut) ;
1302   printf("M_pair                                        <     %f\n", fInvMassMaxCut) ; 
1303   printf("M_pair                                        >     %f\n", fInvMassMinCut) ; 
1304   if(fSelect == 2){
1305     printf("pT Jet / pT Gamma                     <    %f\n", fJetRatioMaxCut) ; 
1306     printf("pT Jet / pT Gamma                     >    %f\n", fJetRatioMinCut) ;
1307     printf("pT Jet (Only CTS)/ pT Gamma   <    %f\n", fJetCTSRatioMaxCut) ; 
1308     printf("pT Jet (Only CTS)/ pT Gamma   >    %f\n", fJetCTSRatioMinCut) ;
1309   }
1310
1311  
1312
1313
1314 //___________________________________________________________________
1315 void AliAnaGammaJet::SetJet(TParticle * part, Bool_t & b, Float_t cone, 
1316                              Double_t eta, Double_t phi)
1317 {
1318
1319   //Check if the particle is inside the cone defined by the leading particle
1320   b = kFALSE;
1321   
1322   if(phi > TMath::TwoPi())
1323     phi-=TMath::TwoPi();
1324   if(phi < 0.)
1325     phi+=TMath::TwoPi();
1326   
1327   Double_t  rad = 10000 + cone;
1328   
1329   if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
1330     rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1331                       TMath::Power(part->Phi()-phi,2));
1332   else{
1333     if(part->Phi()-phi > TMath::TwoPi() - cone)
1334       rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1335                         TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
1336     if(part->Phi()-phi < -(TMath::TwoPi() - cone))
1337       rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1338                         TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
1339   }
1340
1341   if(rad < cone )
1342     b = kTRUE;
1343   
1344 }
1345
1346
1347 void AliAnaGammaJet::Terminate(Option_t *)
1348 {
1349    // The Terminate() function is the last function to be called during
1350    // a query. It always runs on the client, it can be used to present
1351    // the results graphically or save the results to file.
1352     
1353
1354 }