New more general analysis implemention for particle identification and correlation...
[u/mrichter/AliRoot.git] / PWG4 / AliAnaParticleJetLeadingConeCorrelation.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 //_________________________________________________________________________
25 // Class for the analysis of gamma-jet correlations:
26 // 1)Take the trigger particle stored in AliAODParticleCorrelation,
27 // 2) Search for the highest pt leading particle opposite to the photon within a phi, pt window
28 // 3) Take all particles around leading in a cone R with pt larger than threshold and construct the jet
29 //
30 //  Class created from old AliPHOSGammaJet 
31 //  (see AliRoot versions previous Release 4-09)
32 //
33 //*-- Author: Gustavo Conesa (LNF-INFN) 
34 //////////////////////////////////////////////////////////////////////////////
35
36
37 // --- ROOT system ---
38
39 //---- AliRoot system ----
40 #include "AliNeutralMesonSelection.h"
41 #include "AliAnaParticleJetLeadingCone.h"  
42 #include "AliLog.h"
43
44 ClassImp(AliAnaParticleJetLeadingCone)
45
46
47 //____________________________________________________________________________
48   AliAnaParticleJetLeadingCone::AliAnaParticleJetLeadingCone() : 
49     AliAnaBaseClass(),  fPbPb(kFALSE),     
50     fSeveralConeAndPtCuts(0),  
51     fJetCTSRatioMaxCut(0.), fJetCTSRatioMinCut(0.), fJetRatioMaxCut(0.),
52     fJetRatioMinCut(0.), 
53     fJetNCone(0),fJetNPt(0), fJetCone(0), 
54     fJetPtThreshold(0),fJetPtThresPbPb(0),
55     fPtJetSelectionCut(0.0), fSelect(0),
56     fhPhiCharged(0), fhPhiNeutral(0), fhEtaCharged(0), fhEtaNeutral(0), 
57     fhDeltaPhiGammaCharged(0),  fhDeltaPhiGammaNeutral(0), 
58     fhDeltaEtaGammaCharged(0), fhDeltaEtaGammaNeutral(0), 
59     fhAnglePairLeading(), fhInvMassPairLeading(), 
60     fhChargedRatio(0), fhNeutralRatio (0), 
61     fhNBkg (0), fhNLeading(0), fhNJet(0), fhJetRatio(0), fhJetPt (0), 
62     fhBkgRatio (0), fhBkgPt(0),  fhJetFragment(0), fhBkgFragment(0), 
63     fhJetPtDist(0),  fhBkgPtDist(0) 
64 {
65   //Default Ctor
66
67   //Initialize parameters
68
69   SetCorrelationType(kJetLeadCone);
70
71   for(Int_t i = 0; i<10; i++){
72     fJetCones[i]         = 0.0 ;
73     fJetNameCones[i]     = ""  ;
74     fJetPtThres[i]      = 0.0 ;
75     fJetNamePtThres[i]  = ""  ;
76     if( i < 6 ){
77       fJetXMin1[i]     = 0.0 ;
78       fJetXMin2[i]     = 0.0 ;
79       fJetXMax1[i]     = 0.0 ;
80       fJetXMax2[i]     = 0.0 ;
81       fBkgMean[i]      = 0.0 ;
82       fBkgRMS[i]       = 0.0 ;
83       if( i < 2 ){
84         fJetE1[i]        = 0.0 ;
85         fJetE2[i]        = 0.0 ;
86         fJetSigma1[i]    = 0.0 ;
87         fJetSigma2[i]    = 0.0 ;
88       }
89     }
90   }
91
92   InitParameters();
93 }
94
95 //____________________________________________________________________________
96 AliAnaParticleJetLeadingCone::AliAnaParticleJetLeadingCone(const AliAnaParticleJetLeadingCone & g) :   
97   AliAnaBaseClass(g), fPbPb(g.fPbPb), 
98   fSeveralConeAndPtCuts(g.fSeveralConeAndPtCuts), 
99   fJetCTSRatioMaxCut(g.fJetCTSRatioMaxCut),
100   fJetCTSRatioMinCut(g.fJetCTSRatioMinCut), fJetRatioMaxCut(g.fJetRatioMaxCut),
101   fJetRatioMinCut(g.fJetRatioMinCut),  fJetNCone(g.fJetNCone),
102   fJetNPt(g.fJetNPt), fJetCone(g.fJetCone),
103   fJetPtThreshold(g.fJetPtThreshold),fJetPtThresPbPb(g.fJetPtThresPbPb),
104   fPtJetSelectionCut(g.fPtJetSelectionCut), fSelect(g.fSelect),  
105   fhPhiCharged(g.fhPhiCharged), fhPhiNeutral(g.fhPhiNeutral), 
106   fhEtaCharged(g.fhEtaCharged), fhEtaNeutral(g.fhEtaNeutral), 
107   fhDeltaPhiGammaCharged(g.fhDeltaPhiGammaCharged),  
108   fhDeltaPhiGammaNeutral(g.fhDeltaPhiGammaNeutral), 
109   fhDeltaEtaGammaCharged(g.fhDeltaEtaGammaCharged), 
110   fhDeltaEtaGammaNeutral(g.fhDeltaEtaGammaNeutral), 
111   fhAnglePairLeading(g.fhAnglePairLeading), fhInvMassPairLeading(g.fhInvMassPairLeading), 
112   fhChargedRatio(g.fhChargedRatio), fhNeutralRatio(g.fhNeutralRatio), 
113   fhNBkg(g. fhNBkg), fhNLeading(g. fhNLeading), fhNJet(g.fhNJet), fhJetRatio(g.fhJetRatio), fhJetPt(g.fhJetPt), 
114   fhBkgRatio (g.fhBkgRatio), fhBkgPt(g.fhBkgPt),  fhJetFragment(g.fhJetFragment), fhBkgFragment(g.fhBkgFragment), 
115   fhJetPtDist(g.fhJetPtDist),  fhBkgPtDist(g.fhBkgPtDist)   
116 {
117   // cpy ctor
118
119   for(Int_t i = 0; i<10; i++){
120     fJetCones[i]        = g.fJetCones[i] ;
121     fJetNameCones[i]    = g.fJetNameCones[i] ;
122     fJetPtThres[i]      = g.fJetPtThres[i] ;
123     fJetNamePtThres[i]  = g.fJetNamePtThres[i] ;
124     if( i < 6 ){
125       fJetXMin1[i]       = g.fJetXMin1[i] ;
126       fJetXMin2[i]       = g.fJetXMin2[i] ;
127       fJetXMax1[i]       = g.fJetXMax1[i] ;
128       fJetXMax2[i]       = g.fJetXMax2[i] ;
129       fBkgMean[i]        = g.fBkgMean[i] ;
130       fBkgRMS[i]         = g.fBkgRMS[i] ;
131       if( i < 2 ){
132         fJetE1[i]        = g.fJetE1[i] ;
133         fJetE2[i]        = g.fJetE2[i] ;
134         fJetSigma1[i]    = g.fJetSigma1[i] ;
135         fJetSigma2[i]    = g.fJetSigma2[i] ;
136       }
137     }          
138   } 
139
140   
141 }
142
143 //_________________________________________________________________________
144 AliAnaParticleJetLeadingCone & AliAnaParticleJetLeadingCone::operator = (const AliAnaParticleJetLeadingCone & source)
145 {
146   // assignment operator
147
148   if(this == &source)return *this;
149   ((AliAnaBaseClass *)this)->operator=(source);
150
151   fSeveralConeAndPtCuts = source.fSeveralConeAndPtCuts ; 
152   fPbPb = source.fPbPb ;
153   fJetCTSRatioMaxCut = source.fJetCTSRatioMaxCut ;
154   fJetCTSRatioMinCut = source.fJetCTSRatioMinCut ; fJetRatioMaxCut = source.fJetRatioMaxCut ;
155   fJetRatioMinCut = source.fJetRatioMinCut ;  fJetNCone = source.fJetNCone ;
156   fJetNPt = source.fJetNPt ; fJetCone = source.fJetCone ; 
157   fJetPtThreshold = source.fJetPtThreshold ;
158   fJetPtThresPbPb = source.fJetPtThresPbPb ;
159   fPtJetSelectionCut = source.fPtJetSelectionCut ;
160   fSelect = source.fSelect ;  fhChargedRatio = source.fhChargedRatio ; fhNeutralRatio = source.fhNeutralRatio ; 
161
162   fhPhiCharged = source.fhPhiCharged ; fhPhiNeutral = source.fhPhiNeutral ; 
163   fhEtaCharged = source.fhEtaCharged ; fhEtaNeutral = source.fhEtaNeutral ; 
164   fhDeltaPhiGammaCharged = source.fhDeltaPhiGammaCharged ;  
165   fhDeltaPhiGammaNeutral = source.fhDeltaPhiGammaNeutral ; 
166   fhDeltaEtaGammaCharged = source.fhDeltaEtaGammaCharged ; 
167   fhDeltaEtaGammaNeutral = source.fhDeltaEtaGammaNeutral ; 
168   
169   fhAnglePairLeading = source.fhAnglePairLeading ; 
170   fhInvMassPairLeading = source.fhInvMassPairLeading ; 
171   fhNBkg = source. fhNBkg ; fhNLeading = source. fhNLeading ; 
172   fhNJet = source.fhNJet ; fhJetRatio = source.fhJetRatio ; fhJetPt = source.fhJetPt ; 
173   fhBkgRatio  = source.fhBkgRatio ; fhBkgPt = source.fhBkgPt ;  
174   fhJetFragment = source.fhJetFragment ; fhBkgFragment = source.fhBkgFragment ; 
175   fhJetPtDist = source.fhJetPtDist ;  fhBkgPtDist = source.fhBkgPtDist ;
176
177
178   for(Int_t i = 0; i<10; i++){
179     fJetCones[i]        = source.fJetCones[i] ;
180     fJetNameCones[i]    = source.fJetNameCones[i] ;
181     fJetPtThres[i]      = source.fJetPtThres[i] ;
182     fJetNamePtThres[i]  = source.fJetNamePtThres[i] ;
183     if( i < 6 ){
184       fJetXMin1[i]       = source.fJetXMin1[i] ;
185       fJetXMin2[i]       = source.fJetXMin2[i] ;
186       fJetXMax1[i]       = source.fJetXMax1[i] ;
187       fJetXMax2[i]       = source.fJetXMax2[i] ;
188       fBkgMean[i]        = source.fBkgMean[i] ;
189       fBkgRMS[i]         = source.fBkgRMS[i] ;
190       if( i < 2 ){
191         fJetE1[i]        = source.fJetE1[i] ;
192         fJetE2[i]        = source.fJetE2[i] ;
193         fJetSigma1[i]    = source.fJetSigma1[i] ;
194         fJetSigma2[i]    = source.fJetSigma2[i] ;
195
196       }
197     }          
198   } 
199
200   return *this;
201
202 }
203
204 //____________________________________________________________________________
205 AliAnaParticleJetLeadingCone::~AliAnaParticleJetLeadingCone() 
206 {
207    // Remove all pointers except analysis output pointers.
208   delete [] fJetE1;  
209   delete [] fJetE2;    
210   delete [] fJetSigma1;
211   delete [] fJetSigma2;
212   delete [] fBkgMean; 
213   delete [] fBkgRMS;  
214   delete [] fJetXMin1;
215   delete [] fJetXMin2;
216   delete [] fJetXMax1;
217   delete [] fJetXMax2; 
218   delete [] fJetCones;         
219   delete [] fJetNameCones;   
220   delete [] fJetPtThres;       
221   delete [] fJetNamePtThres;  
222 }
223
224
225
226 //________________________________________________________________________
227 TList *  AliAnaParticleJetLeadingCone::GetCreateOutputObjects()
228 {  
229   // Create histograms to be saved in output file and 
230   // store them in fOutputContainer
231  
232   AliDebug(1,"Init jet in leading cone histograms");
233   
234   TList * outputContainer = new TList() ; 
235   outputContainer->SetName("GammaJetCorrelationHistos") ; 
236
237   fhChargedRatio  = new TH2F
238     ("ChargedRatio","p_{T leading charge} /p_{T #gamma} vs p_{T #gamma}",
239      120,0,120,120,0,1); 
240   fhChargedRatio->SetYTitle("p_{T lead charge} /p_{T #gamma}");
241   fhChargedRatio->SetXTitle("p_{T #gamma} (GeV/c)");
242   
243   fhPhiCharged  = new TH2F
244     ("PhiCharged","#phi_{#pi^{#pm}}  vs p_{T #gamma}",
245      120,0,120,120,0,7); 
246   fhPhiCharged->SetYTitle("#phi_{#pi^{#pm}} (rad)");
247   fhPhiCharged->SetXTitle("p_{T #gamma} (GeV/c)");
248   
249   fhEtaCharged  = new TH2F
250     ("EtaCharged","#eta_{#pi^{#pm}}  vs p_{T #gamma}",
251      120,0,120,120,-1,1); 
252   fhEtaCharged->SetYTitle("#eta_{#pi^{#pm}} (rad)");
253   fhEtaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
254   
255   fhDeltaPhiGammaCharged  = new TH2F
256     ("DeltaPhiGammaCharged","#phi_{#gamma} - #phi_{charged #pi} vs p_{T #gamma}",
257      200,0,120,200,0,6.4); 
258   fhDeltaPhiGammaCharged->SetYTitle("#Delta #phi");
259   fhDeltaPhiGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
260   
261   fhDeltaEtaGammaCharged  = new TH2F
262     ("DeltaEtaGammaCharged","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
263      200,0,120,200,-2,2); 
264   fhDeltaEtaGammaCharged->SetYTitle("#Delta #eta");
265   fhDeltaEtaGammaCharged->SetXTitle("p_{T #gamma} (GeV/c)");
266   
267   outputContainer->Add(fhPhiCharged) ;
268   outputContainer->Add(fhEtaCharged) ;
269   outputContainer->Add(fhChargedRatio) ;
270   outputContainer->Add(fhDeltaPhiGammaCharged) ; 
271   outputContainer->Add(fhDeltaEtaGammaCharged) ; 
272
273   if(!AreJetsOnlyInCTS()){
274     
275     fhNeutralRatio  = new TH2F
276       ("NeutralRatio","p_{T leading  #pi^{0}} /p_{T #gamma} vs p_{T #gamma}",
277        120,0,120,120,0,1); 
278     fhNeutralRatio->SetYTitle("p_{T lead  #pi^{0}} /p_{T #gamma}");
279     fhNeutralRatio->SetXTitle("p_{T #gamma} (GeV/c)");
280     
281     //
282     fhAnglePairLeading  = new TH2F
283       ("AnglePairLeading",
284        "Angle between all #gamma pair finally selected vs p_{T  #pi^{0}}",
285        200,0,50,200,0,0.2); 
286     fhAnglePairLeading->SetYTitle("Angle (rad)");
287     fhAnglePairLeading->SetXTitle("E_{ #pi^{0}} (GeV/c)");
288     
289     fhInvMassPairLeading  = new TH2F
290       ("InvMassPairLeading",
291        "Invariant Mass of #gamma pair selected vs p_{T #gamma}",
292        120,0,120,360,0,0.5); 
293     fhInvMassPairLeading->SetYTitle("Invariant Mass (GeV/c^{2})");
294     fhInvMassPairLeading->SetXTitle("p_{T #gamma} (GeV/c)");
295
296     fhPhiNeutral  = new TH2F
297       ("PhiNeutral","#phi_{#pi^{0}}  vs p_{T #gamma}",
298        120,0,120,120,0,7); 
299     fhPhiNeutral->SetYTitle("#phi_{#pi^{0}} (rad)");
300     fhPhiNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
301     
302     fhEtaNeutral  = new TH2F
303       ("EtaNeutral","#eta_{#pi^{0}}  vs p_{T #gamma}",
304        120,0,120,120,-1,1); 
305     fhEtaNeutral->SetYTitle("#eta_{#pi^{0}} (rad)");
306     fhEtaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
307     
308     fhDeltaPhiGammaNeutral  = new TH2F
309       ("DeltaPhiGammaNeutral","#phi_{#gamma} - #phi_{#pi^{0}} vs p_{T #gamma}",
310        200,0,120,200,0,6.4); 
311     fhDeltaPhiGammaNeutral->SetYTitle("#Delta #phi");
312     fhDeltaPhiGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
313     
314     fhDeltaEtaGammaNeutral  = new TH2F
315       ("DeltaEtaGammaNeutral","#eta_{#gamma} - #eta_{#pi^{#pm}} vs p_{T #gamma}",
316        200,0,120,200,-2,2); 
317     fhDeltaEtaGammaNeutral->SetYTitle("#Delta #eta");
318     fhDeltaEtaGammaNeutral->SetXTitle("p_{T #gamma} (GeV/c)");
319
320     outputContainer->Add(fhPhiNeutral) ;  
321     outputContainer->Add(fhEtaNeutral) ;  
322     outputContainer->Add(fhNeutralRatio) ; 
323     outputContainer->Add(fhDeltaPhiGammaNeutral) ; 
324     outputContainer->Add(fhDeltaEtaGammaNeutral) ;
325     
326     outputContainer->Add(fhInvMassPairLeading) ; 
327     outputContainer->Add(fhAnglePairLeading) ; 
328   }
329   
330   if(!fSeveralConeAndPtCuts){// not several cones
331     
332     //Count
333     fhNBkg = new TH1F("NBkg","bkg multiplicity",9000,0,9000); 
334     fhNBkg->SetYTitle("counts");
335     fhNBkg->SetXTitle("N");
336     outputContainer->Add(fhNBkg) ; 
337     
338     fhNLeading  = new TH2F
339       ("NLeading","Accepted Jet Leading", 240,0,120,240,0,120); 
340     fhNLeading->SetYTitle("p_{T charge} (GeV/c)");
341     fhNLeading->SetXTitle("p_{T #gamma}(GeV/c)");
342     outputContainer->Add(fhNLeading) ; 
343     
344     fhNJet  = new TH1F("NJet","Accepted jets",240,0,120); 
345     fhNJet->SetYTitle("N");
346     fhNJet->SetXTitle("p_{T #gamma}(GeV/c)");
347     outputContainer->Add(fhNJet) ; 
348     
349     //Ratios and Pt dist of reconstructed (not selected) jets
350     //Jet
351     fhJetRatio  = new TH2F
352       ("JetRatio","p_{T jet lead}/p_{T #gamma} vs p_{T #gamma}",
353        240,0,120,200,0,10);
354     fhJetRatio->SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
355     fhJetRatio->SetXTitle("p_{T #gamma} (GeV/c)");
356     outputContainer->Add(fhJetRatio) ; 
357     
358     fhJetPt  = new TH2F
359       ("JetPt", "p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
360     fhJetPt->SetYTitle("p_{T jet}");
361     fhJetPt->SetXTitle("p_{T #gamma} (GeV/c)");
362     outputContainer->Add(fhJetPt) ; 
363     
364     //Bkg
365     
366     fhBkgRatio  = new TH2F
367       ("BkgRatio","p_{T bkg lead}/p_{T #gamma} vs p_{T #gamma}",
368        240,0,120,200,0,10);
369     fhBkgRatio->SetYTitle("p_{T bkg lead charge}/p_{T #gamma}");
370     fhBkgRatio->SetXTitle("p_{T #gamma} (GeV/c)");
371     outputContainer->Add(fhBkgRatio) ;
372     
373     fhBkgPt  = new TH2F
374       ("BkgPt","p_{T jet lead} vs p_{T #gamma}",240,0,120,400,0,200);
375     fhBkgPt->SetYTitle("p_{T jet lead charge}/p_{T #gamma}");
376     fhBkgPt->SetXTitle("p_{T #gamma} (GeV/c)");
377     outputContainer->Add(fhBkgPt) ;
378     
379     //Jet Distributions
380     
381     fhJetFragment  = 
382       new TH2F("JetFragment","x = p_{T i charged}/p_{T #gamma}",
383                240,0.,120.,1000,0.,1.2); 
384     fhJetFragment->SetYTitle("x_{T}");
385     fhJetFragment->SetXTitle("p_{T #gamma}");
386     outputContainer->Add(fhJetFragment) ;
387     
388     fhBkgFragment  = new TH2F
389       ("BkgFragment","x = p_{T i charged}/p_{T #gamma}",
390        240,0.,120.,1000,0.,1.2);
391     fhBkgFragment->SetYTitle("x_{T}");
392     fhBkgFragment->SetXTitle("p_{T #gamma}");
393     outputContainer->Add(fhBkgFragment) ;
394     
395     fhJetPtDist  = 
396       new TH2F("JetPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); 
397     fhJetPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
398     outputContainer->Add(fhJetPtDist) ;
399     
400     fhBkgPtDist  = new TH2F
401       ("BkgPtDist","x = p_{T i charged}",240,0.,120.,400,0.,200.); 
402     fhBkgPtDist->SetXTitle("p_{T #gamma} (GeV/c)");
403     outputContainer->Add(fhBkgPtDist) ;
404     
405   }//not several cones
406   else{ //If we want to study the jet for different cones and pt
407     for(Int_t icone = 0; icone<fJetNCone; icone++){//icone
408       for(Int_t ipt = 0; ipt<fJetNPt;ipt++){ //ipt
409         
410         //Jet
411         
412         fhJetRatios[icone][ipt]  = new TH2F
413           ("JetRatioCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], 
414            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
415            +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c",
416            240,0,120,200,0,10);
417         fhJetRatios[icone][ipt]->
418           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
419         fhJetRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
420         outputContainer->Add(fhJetRatios[icone][ipt]) ; 
421         
422         
423         fhJetPts[icone][ipt]  = new TH2F
424           ("JetPtCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], 
425            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
426            +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c",
427            240,0,120,400,0,200);
428         fhJetPts[icone][ipt]->
429           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
430         fhJetPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
431         outputContainer->Add(fhJetPts[icone][ipt]) ; 
432         
433         //Bkg
434         fhBkgRatios[icone][ipt]  = new TH2F
435           ("BkgRatioCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], 
436            "p_{T bkg lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
437            +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c",
438            240,0,120,200,0,10);
439         fhBkgRatios[icone][ipt]->
440           SetYTitle("p_{T bkg lead #pi^{0}}/p_{T #gamma}");
441         fhBkgRatios[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
442         outputContainer->Add(fhBkgRatios[icone][ipt]) ; 
443         
444         fhBkgPts[icone][ipt]  = new TH2F
445           ("BkgPtCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt], 
446            "p_{T jet lead #pi^{0}}/p_{T #gamma} vs p_{T #gamma}, cone ="
447            +fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+" GeV/c",
448            240,0,120,400,0,200);
449         fhBkgPts[icone][ipt]->
450           SetYTitle("p_{T jet lead #pi^{0}}/p_{T #gamma}");
451         fhBkgPts[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
452         outputContainer->Add(fhBkgPts[icone][ipt]) ; 
453         
454         //Counts
455         fhNBkgs[icone][ipt]  = new TH1F
456           ("NBkgCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
457            "bkg multiplicity cone ="+fJetNameCones[icone]+", pt>" 
458            +fJetNamePtThres[ipt]+" GeV/c",9000,0,9000); 
459         fhNBkgs[icone][ipt]->SetYTitle("counts");
460         fhNBkgs[icone][ipt]->SetXTitle("N");
461         outputContainer->Add(fhNBkgs[icone][ipt]) ; 
462         
463         fhNLeadings[icone][ipt]  = new TH2F
464           ("NLeadingCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
465            "p_{T #gamma} vs p_{T #pi^{0}} cone ="+fJetNameCones[icone]+", pt>" 
466            +fJetNamePtThres[ipt]+" GeV/c",120,0,120,120,0,120); 
467         fhNLeadings[icone][ipt]->SetYTitle("p_{T #pi^{0}}(GeV/c)");
468         fhNLeadings[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
469         outputContainer->Add(fhNLeadings[icone][ipt]) ; 
470         
471         fhNJets[icone][ipt]  = new TH1F
472           ("NJetCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
473            "Number of neutral jets, cone ="+fJetNameCones[icone]+", pt>" 
474            +fJetNamePtThres[ipt]+" GeV/c",120,0,120); 
475         fhNJets[icone][ipt]->SetYTitle("N");
476         fhNJets[icone][ipt]->SetXTitle("p_{T #gamma}(GeV/c)");
477         outputContainer->Add(fhNJets[icone][ipt]) ; 
478         
479         //Fragmentation Function
480         fhJetFragments[icone][ipt]  = new TH2F
481           ("JetFragmentCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
482            "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fJetNameCones[icone]+", pt>" 
483            +fJetNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); 
484         fhJetFragments[icone][ipt]->SetYTitle("x_{T}");
485         fhJetFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
486         outputContainer->Add(fhJetFragments[icone][ipt]) ; 
487         
488         fhBkgFragments[icone][ipt]  = new TH2F
489           ("BkgFragmentCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
490            "x_{T} = p_{T i}/p_{T #gamma}, cone ="+fJetNameCones[icone]+", pt>" 
491            +fJetNamePtThres[ipt]+" GeV/c",120,0.,120.,240,0.,1.2); 
492         fhBkgFragments[icone][ipt]->SetYTitle("x_{T}");
493         fhBkgFragments[icone][ipt]->SetXTitle("p_{T #gamma}");
494         outputContainer->Add(fhBkgFragments[icone][ipt]) ; 
495         
496         //Jet particle distribution
497         
498         fhJetPtDists[icone][ipt]  = new TH2F
499           ("JetPtDistCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
500            "p_{T i}, cone ="+fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+
501            " GeV/c",120,0.,120.,120,0.,120.); 
502         fhJetPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
503         outputContainer->Add(fhJetPtDists[icone][ipt]) ; 
504         
505         fhBkgPtDists[icone][ipt]  = new TH2F
506           ("BkgPtDistCone"+fJetNameCones[icone]+"Pt"+fJetNamePtThres[ipt],
507            "p_{T i}, cone ="+fJetNameCones[icone]+", pt>" +fJetNamePtThres[ipt]+
508            " GeV/c",120,0.,120.,120,0.,120.); 
509         fhBkgPtDists[icone][ipt]->SetXTitle("p_{T #gamma} (GeV/c)");
510         outputContainer->Add(fhBkgPtDists[icone][ipt]) ; 
511         
512       }//ipt
513     } //icone
514   }//If we want to study any cone or pt threshold
515
516   return outputContainer;
517 }
518
519   //____________________________________________________________________________
520   void AliAnaParticleJetLeadingCone::InitParameters()
521 {
522  
523   //Initialize the parameters of the analysis.
524   SetJetsOnlyInCTS(kFALSE) ;
525   fPbPb                = kFALSE ;
526
527   SetDeltaPhiCutRange(2.9,3.4) ; 
528   SetRatioCutRange(0.1,1.0) ; 
529
530   //Jet selection parameters
531   //Fixed cut   
532   fJetRatioMaxCut = 1.2 ; 
533   fJetRatioMinCut = 0.3 ; 
534   fJetCTSRatioMaxCut = 1.2 ;
535   fJetCTSRatioMinCut = 0.3 ;
536   fSelect         = 0  ;
537
538   //Cut depending on gamma energy
539
540   fPtJetSelectionCut = 20.; //For Low pt jets+BKG, another limits applied
541   //Reconstructed jet energy dependence parameters 
542   //e_jet = a1+e_gamma b2. 
543   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
544   fJetE1[0] = -5.75; fJetE1[1] = -4.1;
545   fJetE2[0] = 1.005; fJetE2[1] = 1.05;
546
547   //Reconstructed sigma of jet energy dependence parameters 
548   //s_jet = a1+e_gamma b2. 
549   //Index 0-> Pt>2 GeV r = 0.3; Index 1-> Pt>0.5 GeV r = 0.3
550   fJetSigma1[0] = 2.65;   fJetSigma1[1] = 2.75;
551   fJetSigma2[0] = 0.0018; fJetSigma2[1] = 0.033;
552
553   //Background mean energy and RMS
554   //Index 0-> No BKG; Index 1-> BKG > 2 GeV; 
555   //Index 2-> (low pt jets)BKG > 0.5 GeV;
556   //Index > 2, same for CTS conf
557   fBkgMean[0] = 0.; fBkgMean[1] = 8.8 ; fBkgMean[2] = 69.5;
558   fBkgMean[3] = 0.; fBkgMean[4] = 6.4;  fBkgMean[5] = 48.6;
559   fBkgRMS[0]  = 0.; fBkgRMS[1]  = 7.5;  fBkgRMS[2]  = 22.0; 
560   fBkgRMS[3]  = 0.; fBkgRMS[4]  = 5.4;  fBkgRMS[5]  = 13.2; 
561
562   //Factor x of min/max = E -+ x * sigma. Obtained after selecting the
563   //limits for monoenergetic jets.
564   //Index 0-> No BKG; Index 1-> BKG > 2 GeV; 
565   //Index 2-> (low pt jets) BKG > 0.5 GeV;
566   //Index > 2, same for CTS conf
567
568   fJetXMin1[0] =-0.69 ; fJetXMin1[1] = 0.39 ; fJetXMin1[2] =-0.88 ; 
569   fJetXMin1[3] =-2.0  ; fJetXMin1[4] =-0.442 ; fJetXMin1[5] =-1.1  ;
570   fJetXMin2[0] = 0.066; fJetXMin2[1] = 0.038; fJetXMin2[2] = 0.034; 
571   fJetXMin2[3] = 0.25 ; fJetXMin2[4] = 0.113; fJetXMin2[5] = 0.077 ;
572   fJetXMax1[0] =-3.8  ; fJetXMax1[1] =-0.76 ; fJetXMax1[2] =-3.6  ; 
573   fJetXMax1[3] =-2.7  ; fJetXMax1[4] =-1.21 ; fJetXMax1[5] =-3.7  ;
574   fJetXMax2[0] =-0.012; fJetXMax2[1] =-0.022; fJetXMax2[2] = 0.016; 
575   fJetXMax2[3] =-0.024; fJetXMax2[4] =-0.008; fJetXMax2[5] = 0.027;
576
577
578   //Different cones and pt thresholds to construct the jet
579
580   fJetCone        = 0.3  ;
581   fJetPtThreshold = 0.5   ;
582   fJetPtThresPbPb = 2.   ;
583   fJetNCone       = 4    ;
584   fJetNPt         = 4    ;
585   fJetCones[0]    = 0.2  ; fJetNameCones[0]   = "02" ;
586   fJetCones[1]    = 0.3  ; fJetNameCones[1]   = "03" ;
587   fJetCones[2]    = 0.4  ; fJetNameCones[2]   = "04" ;
588   fJetCones[2]    = 0.5  ; fJetNameCones[2]   = "05" ;
589
590   fJetPtThres[0]  = 0.0  ; fJetNamePtThres[0] = "00" ;
591   fJetPtThres[1]  = 0.5  ; fJetNamePtThres[1] = "05" ;
592   fJetPtThres[2]  = 1.0  ; fJetNamePtThres[2] = "10" ;
593   fJetPtThres[3]  = 2.0  ; fJetNamePtThres[3] = "20" ;
594 }
595
596 //__________________________________________________________________
597 void AliAnaParticleJetLeadingCone::Print(const Option_t * opt) const
598 {
599
600   //Print some relevant parameters set for the analysis
601   if(! opt)
602     return;
603   
604   Info("Print", "%s %s", GetName(), GetTitle() ) ;
605   printf("Correlation analysis           =     %d\n",kJetLeadCone) ;
606   
607   
608   printf("Phi gamma-Leading      <     %f\n", GetDeltaPhiMaxCut()) ; 
609   printf("Phi gamma-Leading      >     %f\n", GetDeltaPhiMinCut()) ;
610   printf("pT Leading / pT Gamma             <     %f\n", GetRatioMaxCut()) ; 
611   printf("pT Leading / pT Gamma             >     %f\n", GetRatioMinCut()) ;
612   
613   if(fSelect == 2){
614     printf("pT Jet / pT Gamma                     <    %f\n", fJetRatioMaxCut) ; 
615     printf("pT Jet / pT Gamma                     >    %f\n", fJetRatioMinCut) ;
616     printf("pT Jet (Only CTS)/ pT Gamma   <    %f\n", fJetCTSRatioMaxCut) ; 
617     printf("pT Jet (Only CTS)/ pT Gamma   >    %f\n", fJetCTSRatioMinCut) ;
618   }
619   
620
621
622 //__________________________________________________________________
623 void  AliAnaParticleJetLeadingCone::MakeAnalysisFillAOD() 
624 {
625   
626   //Particle-Hadron Correlation Analysis, fill AODs
627   if(GetDebug() > 1){
628     printf("Begin jet leading cone  correlation analysis, fill AODs \n");
629     printf("In particle branch aod entries %d\n", GetAODBranch()->GetEntries());
630     printf("In CTS aod entries %d\n", GetAODCTS()->GetEntries());
631     printf("In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntries());
632   }
633  
634   TLorentzVector * pLeading = new TLorentVector; //It will contain the kinematics of the found leading particle
635   
636   //Loop on stored AOD particles, trigger
637   Int_t naod = GetAODBranch()->GetEntriesFast();
638   for(Int_t iaod = 0; iaod < naod ; iaod++){
639     AliAODParticleCorrelation* particle =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
640     
641     //Search leading particles in CTS and EMCAL 
642     if(GetLeadingParticle(particle, pLeading)){
643       
644       if(GetDebug() > 1) printf("Leading: pt %f, phi %f, eta %f", pLeading->Pt(),pLeading->Phi(),pLeading->Eta())) ;
645       
646     MakeJet(particle, pLeading,"", kFALSE);
647
648     }//Leading
649   }//AOD trigger particle loop
650   
651   if(GetDebug() >1)printf("End of jet leading cone analysis, fill AODs \n");
652
653
654
655 //__________________________________________________________________
656 void  AliAnaParticleJetLeadingCone::MakeAnalysisFillHistograms() 
657 {
658   
659   //Particle-Hadron Correlation Analysis, fill AODs
660   if(GetDebug() > 1){
661     printf("Begin jet leading cone  correlation analysis, fill histograms \n");
662     printf("In particle branch aod entries %d\n", GetAODBranch()->GetEntries());
663     printf("In CTS aod entries %d\n", GetAODCTS()->GetEntries());
664     printf("In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntries());
665   }
666  
667   TLorentzVector * pLeading = new TLorentzVector;
668   
669   //Loop on stored AOD particles, trigger
670   Int_t naod = GetAODBranch()->GetEntriesFast();
671   for(Int_t iaod = 0; iaod < naod ; iaod++){
672     AliAODParticleCorrelation* particle =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
673     
674     Double_t pt = particle->Pt();
675     Double_t phi = particle->Phi();
676     Double_t eta = particle->Eta();
677     
678     //Get leading particle, fill histograms
679     pLeading = particle->GetLeading();
680     TString det = particle->GetLeadingDetector();      
681
682     if(det!="" && pLeading){
683       Double_t ptL = pLeading->Pt(); 
684       Double_t phiL = pLeading->Phi(); 
685       if(phiL < 0 ) phiL+=TMath::TwoPi();
686       Double_t etaL = pLeading->Eta(); 
687       
688       if(GetDebug() > 1) printf("Leading found in %s, with pt %3.2f, phi %2.2f, eta %2.2f",det.Data(), ptL, phiL, etaL);
689       if(det == "CTS"){
690         fhChargedRatio->Fill(pt,ptL/particle->Pt());
691         fhEtaCharged->Fill(pt,etaL);
692         fhPhiCharged->Fill(pt,phiL);
693         fhDeltaPhiGammaCharged->Fill(pt,phi-phiL);
694         fhDeltaEtaGammaCharged->Fill(pt,eta-etaL);
695       }
696       else if(det== "EMCAL"){
697         fhNeutralRatio->Fill(pt,ptL/particle->Pt());
698         fhEtaNeutral->Fill(pt,etaL);
699         fhPhiNeutral->Fill(pt,phiL);
700         fhDeltaPhiGammaNeutral->Fill(pt,phi-phiL);
701         fhDeltaEtaGammaNeutral->Fill(pt,eta-etaL);
702       }
703     }
704       
705       //Search Jet
706       if(!fSeveralConeAndPtCuts)
707         MakeJet(particle, pLeading, "", kTRUE);
708       else{
709         for(Int_t icone = 0; icone<fJetNCone; icone++) {
710           for(Int_t ipt = 0; ipt<fJetNPt;ipt++) {  
711             TString lastname ="Cone"+ fJetNameCones[icone]+"Pt"+ fJetNamePtThres[ipt];
712             fJetCone=fJetCones[icone];
713             fJetPtThreshold=fJetPtThres[ipt];
714             MakeJet(particle, pLeading, lastname);
715           }//icone
716         }//ipt
717       }//fSeveralConeAndPtCuts
718     }//Leading
719   }//AOD trigger particle loop
720   
721   if(GetDebug() >1)printf("End of jet leading cone analysis, fill histograms \n");
722
723
724
725 //____________________________________________________________________________
726 Bool_t  AliAnaParticleJetLeadingCone::GetLeadingParticle(AliAODParticleCorrelation *particle, TLorentzVector * pLeading) 
727 {
728   //Search Charged or Neutral leading particle, select the highest one.
729   
730   TLorentzVector *pLeadingCh = new TLorentzVector;
731   TLorentzVector *pLeadingPi0 = new TLorentzVector;
732   
733   Double_t pt  =  particle->Pt(); 
734   Double_t phi = particle->Phi(); 
735   Double_t eta = particle->Eta(); 
736
737   GetLeadingCharge(particle, pLeadingCh) ;
738   if(!AreJetsOnlyInCTS()) GetLeadingPi0(particle, pLeadingPi0) ;
739       
740   Double_t ptch = pLeadingCh->Pt(); 
741   Double_t phich = pLeadingCh->Phi(); 
742   if(phich < 0 ) phich+=TMath::TwoPi();
743   Double_t etach = pLeadingCh->Eta(); 
744   Double_t ptpi = pLeadingPi0->Pt(); 
745   Double_t phipi = pLeadingPi0->Phi(); 
746   if(phipi < 0 ) phipi+=TMath::TwoPi();
747   Double_t etapi = pLeadingPi0->Eta(); 
748       
749   if (ptch > 0 || ptpi > 0){
750     if((ptch >= ptpi)){
751       if(GetDebug() > 1)printf("Leading found in CTS \n");
752       pLeading = pLeadingCh;
753       particle->SetLeading(pLeading);
754       return kTRUE;
755     }
756     else{
757       if(GetDebug() > 1)printf("Leading found in EMCAL \n");
758       pLeading = pLeadingPi0;
759       particle->SetLeading(pLeading);    
760       return kTRUE;   
761     }
762     
763     if(GetDebug() > 1)printf ("NO LEADING PARTICLE FOUND \n");
764     return kFALSE; 
765   
766 }
767
768 //____________________________________________________________________________
769 void  AliAnaParticleJetLeadingCone::GetLeadingCharge(AliAODParticleCorrelation * particle, TLorentzVector * pLeading) 
770 {  
771   //Search for the charged particle with highest pt and with 
772   //Phi=Phi_trigger-Pi and pT=0.1E_gamma 
773
774   if(GetAODCTS()){
775     Double_t ptTrig = particle->Pt();
776     Double_t phiTrig = particle->Phi();
777     Double_t rat  = -100 ;
778     Double_t ptl  = -100 ;
779     Double_t phil = -100 ;
780     Double_t pt  = -100.;
781     Double_t phi = -100.;
782     TVector3 p3;
783     
784     for(Int_t ipr = 0;ipr < GetAODCTS()->GetEntries() ; ipr ++ ){
785       AliAODTrack* track = dynamic_cast<AliAODTrack *>(GetAODCTS()->At(ipr)) ;
786       p3.SetXYZ(track->Px(),track->Py(),track->Pz());
787       pt    = p3.Pt();
788       phi  = p3.Phi() ;
789       if(phi<0) phi+=TMath::TwoPi();
790       rat  = ptl/particle->Pt() ;
791       
792       //Selection within angular and energy limits
793       if(((phiTrig-phi) > GetDeltaPhiMinCut()) && ((phiTrig-phi)<GetDeltaPhiMaxCut()) &&
794          (rat > GetRatioMinCut()) && (rat < GetRatioMaxCut())  && (pt  > ptl)) {
795         phil = phi ;
796         ptl  = pt ;
797         pLeading->SetMomentum(p3.Px(),p3.Py(),p3.Pz(),0);
798       }
799     }// track loop
800   }//CTS list exist
801   
802   if(GetDebug() > 1)printf("Leading in CTS: pt %f eta %f phi %f pt/ptTrig %f \n", ptl, pLeading->Eta(), phil,ptl/ptTrig)) ;
803  
804 }
805
806 //____________________________________________________________________________
807 void  AliAnaParticleJetLeadingCone::GetLeadingPi0(AliAODParticleCorrelation * particle, TLorentzVector * pLeading) 
808 {  
809   //Search for the neutral pion with highest pt and with 
810   //Phi=Phi_trigger-Pi and pT=0.1E_gamma 
811   if(GetAODEMCAL()){
812     Double_t ptTrig = particle->Pt();
813     Double_t phiTrig = particle->Phi();
814     Double_t rat  = -100 ;
815     Double_t ptl  = -100 ;
816     Double_t phil = -100 ;
817     Double_t pt  = -100.;
818     Double_t phi = -100.;
819     
820     TLorentzVector gammai;
821     TLorentzVector gammaj;
822     
823     Double_t vertex[] = {0,0,0};
824     
825     if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
826     
827     //Cluster loop, select pairs with good pt, phi and fill AODs or histograms
828     for(Int_t iclus = 0;iclus < GetAODEMCAL()->GetEntries() ; iclus ++ ){
829       AliAODCaloCluster * calo = dynamic_cast< AliAODCaloCluster *>(GetAODEMCAL()->At(iclus)) ;
830       
831       //Cluster selection, not charged, with photon or pi0 id and in fidutial cut
832       Int_t pdg=0;
833       if(!SelectCluster(calo2,vertex, gammaj, pdgj)) continue ;
834       
835       if(GetDebug() > 2) printf("neutral cluster: pt %f, phi %f \n", gammai.Pt(),gammai.Phi());
836         
837       //2 gamma overlapped, found with PID
838       if(pdg == AliCaloPID::kPi0){ 
839         pt = gammai->Pt();
840         rat = pt/ptTrig;
841         phi = gammai->Phi();
842         if(phi<0) phi+=TMath::TwoPi();
843         
844         //Selection within angular and energy limits
845         if(ptl > pt  && rat > GetRatioMinCut()  && rat < GetRatioMaxCut()  && 
846            (phig-phil) > GetDeltaPhiMinCut() && (phig-phil) < GetDeltaPhiMaxCut() )
847           {
848             phi = phil ;
849             pt  = ptl ;
850             pLeading->SetMomentum(gammai.Px(),gammai.Py(),gammai.Pz(),gammai.E());
851           }// cuts
852       }// pdg = 111 
853       //Make invariant mass analysis
854       else if(pdg == AliCaloPID::kPhoton){      
855         //Search the photon companion in case it comes from  a Pi0 decay
856         //Apply several cuts to select the good pair;
857         for(Int_t jclus = iclus+1; jclus < GetAODEMCAL()->GetEntries() ; jclus ++ ){
858           AliAODCaloCluster * calo2 = dynamic_cast< AliAODCaloCluster *>(GetAODEMCAL()->At(jclus)) ;
859           
860           //Cluster selection, not charged with photon or pi0 id and in fidutial cut
861           Int_t pdgj=0;
862           if(!SelectCluster(calo2,vertex, gammaj, pdgj)) continue ;
863           
864           if(pdgj == AliCaloPID::kPhoton ){
865             
866             pt  = (gammai+gammaj).Pt();
867             phi = (gammai+gammaj).Phi();
868             rat = pt/ptTrig;
869             
870             //Selection within angular and energy limits
871             if(ptl > pt  && rat > GetRatioMinCut()  && rat < GetRatioMaxCut()  && 
872                (phig-phil) > GetDeltaPhiMinCut() && (phig-phil) < GetDeltaPhiMaxCut() ){
873               //Select good pair (aperture and invariant mass)
874               if(GetNeutralMesonSelection()->SelectPair(gammai, gammaj)){
875                 phi = phil ;
876                 pt  = ptl ;
877                 pLeading->SetMomentum(gammai.Px(),gammai.Py(),gammai.Pz(),gammai.E());
878               }//pi0 selection
879               
880               if(GetDebug() > 3 ) printf("Neutral Hadron Correlation: Selected gamma pair: pt %2.2f, phi %2.2f, eta %2.2f, M %2.3f",
881                                          (gammai+gammaj).Pt(),(gammai+gammaj).Phi(),(gammai+gammaj).Eta(), (gammai+gammaj).M());
882             }//Pair selected as leading
883           }//if pair of gammas
884         }//2nd loop
885       }// if pdg = 22
886     }// 1st Loop
887   }//EMCAL list exists
888
889   if(GetDebug()>2) printf("Leading EMCAL: pt %f eta %f phi %f pt/Eg %f \n",  pLeading->Pt(), pLeading->Eta(),  pLeading->Phi(),  pLeading->Pt()/ptTrig) ;
890   
891 }
892
893
894 //__________________________________________________________________________-
895 Bool_t AliAnaParticleJetLeadingCone::IsJetSelected(const Double_t ptg, const Double_t ptj){
896   //Check if the energy of the reconstructed jet is within an energy window
897
898   Double_t par[6];
899   Double_t xmax[2];
900   Double_t xmin[2];
901
902   Int_t iCTS = 0;
903   if(AreJetsOnlyInCTS())
904     iCTS = 3 ;
905
906   if(!fPbPb){
907     //Phythia alone, jets with pt_th > 0.2, r = 0.3 
908     par[0] = fJetE1[0]; par[1] = fJetE2[0]; 
909     //Energy of the jet peak
910     //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
911     par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
912     //Sigma  of the jet peak
913     //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
914     par[4] = fBkgMean[0 + iCTS]; par[5] = fBkgRMS[0 + iCTS];
915     //Parameters reserved for PbPb bkg.
916     xmax[0] = fJetXMax1[0 + iCTS]; xmax[1] = fJetXMax2[0 + iCTS];
917     xmin[0] = fJetXMin1[0 + iCTS]; xmin[1] = fJetXMin2[0 + iCTS];
918     //Factor that multiplies sigma to obtain the best limits, 
919     //by observation, of mono jet ratios (ptjet/ptg)
920     //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
921    
922   }
923   else{
924     if(ptg > fPtJetSelectionCut){
925       //Phythia +PbPb with  pt_th > 2 GeV/c, r = 0.3 
926       par[0] = fJetE1[0]; par[1] = fJetE2[0]; 
927       //Energy of the jet peak, same as in pp
928       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
929       par[2] = fJetSigma1[0]; par[3] = fJetSigma2[0];
930       //Sigma  of the jet peak, same as in pp
931       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
932       par[4] = fBkgMean[1 + iCTS]; par[5] = fBkgRMS[1 + iCTS];
933       //Mean value and RMS of PbPb Bkg 
934       xmax[0] = fJetXMax1[1 + iCTS]; xmax[1] = fJetXMax2[1 + iCTS];
935       xmin[0] = fJetXMin1[1 + iCTS]; xmin[1] = fJetXMin2[1 + iCTS];
936       //Factor that multiplies sigma to obtain the best limits, 
937       //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, 
938       //pt_th > 2 GeV, r = 0.3
939       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
940      
941     }
942     else{
943       //Phythia + PbPb with  pt_th > 0.5 GeV/c, r = 0.3
944       par[0] = fJetE1[1]; par[1] = fJetE2[1]; 
945       //Energy of the jet peak, pt_th > 2 GeV/c, r = 0.3 
946       //e_jet = fJetE1[0]+fJetE2[0]*e_gamma, simulation fit
947       par[2] = fJetSigma1[1]; par[3] = fJetSigma2[1];
948       //Sigma  of the jet peak, pt_th > 2 GeV/c, r = 0.3
949       //sigma_jet = fJetSigma1[0]+fJetSigma2[0]*e_gamma, simulation fit
950       par[4] = fBkgMean[2 + iCTS]; par[5] = fBkgRMS[2 + iCTS];
951       //Mean value and RMS of PbPb Bkg in a 0.3 cone, pt > 2 GeV.
952       xmax[0] = fJetXMax1[2 + iCTS]; xmax[1] = fJetXMax2[2 + iCTS];
953       xmin[0] = fJetXMin1[2 + iCTS]; xmin[1] = fJetXMin2[2 + iCTS];
954       //Factor that multiplies sigma to obtain the best limits, 
955       //by observation, of mono jet ratios (ptjet/ptg) mixed with PbPb Bkg, 
956       //pt_th > 2 GeV, r = 0.3
957       //X_jet = fJetX1[0]+fJetX2[0]*e_gamma
958      
959     }//If low pt jet in bkg
960   }//if Bkg
961
962  //Calculate minimum and maximum limits of the jet ratio.
963   Double_t min = CalculateJetRatioLimit(ptg, par, xmin);
964   Double_t max = CalculateJetRatioLimit(ptg, par, xmax);
965   
966   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));
967
968   if(( min < ptj/ptg ) && ( max > ptj/ptg))
969     return kTRUE;
970   else
971     return kFALSE;
972
973 }
974
975 //____________________________________________________________________________
976 void AliAnaParticleJetLeadingCone::MakeJet(AliAODParticleCorrelation *particle, TParticle* pLeading,TString lastname)
977 {
978   //Fill the jet with the particles around the leading particle with 
979   //R=fJetCone and pt_th = fJetPtThres. Calculate the energy of the jet and 
980   //check if we select it. Fill jet histograms
981   
982   TClonesArray * jetList = new TClonesArray("TParticle",1000);
983   TClonesArray * bkgList = new TClonesArray("TParticle",1000);
984
985   TLorentzVector jet   (0,0,0,0);  
986   TLorentzVector bkg(0,0,0,0);
987   TLorentzVector lv (0,0,0,0);
988
989   Double_t ptjet = 0.0;
990   Double_t ptbkg = 0.0;
991   Int_t n0 = 0;
992   Int_t n1 = 0;  
993   Bool_t b1 = kFALSE;
994   Bool_t b0 = kFALSE;
995   
996   Double_t ptg  = pGamma->Pt();
997   Double_t phig = pGamma->Phi();
998   Double_t ptl  = pLeading->Pt();
999   Double_t phil = pLeading->Phi();
1000   Double_t etal = pLeading->Eta();
1001
1002   Float_t ptcut = fJetPtThreshold;
1003   if(fPbPb && !fSeveralConeAndPtCuts && ptg > fPtJetSelectionCut)  ptcut = fJetPtThresPbPb ;
1004  
1005   //Add charged particles to jet
1006   TIter nextch(plCTS) ; 
1007   TParticle * particle = 0 ; 
1008   while ( (particle = dynamic_cast<TParticle*>(nextch())) ) {
1009     
1010     b0 = kFALSE;
1011     b1 = kFALSE;
1012
1013     //Particles in jet 
1014     SetJet(particle, b0, fJetCone, etal, phil) ;  
1015
1016     if(b0){
1017       new((*jetList)[n0++]) TParticle(*particle) ;
1018       particle->Momentum(lv);
1019       if(particle->Pt() > ptcut ){
1020         jet+=lv;
1021         ptjet+=particle->Pt();
1022       }
1023     }
1024
1025     //Background around (phi_gamma-pi, eta_leading)
1026     SetJet(particle, b1, fJetCone,etal, phig) ;
1027
1028     if(b1) { 
1029       new((*bkgList)[n1++]) TParticle(*particle) ;
1030       particle->Momentum(lv);
1031       if(particle->Pt() > ptcut ){
1032         bkg+=lv;
1033         ptbkg+=particle->Pt();    
1034       }  
1035     }
1036   }
1037
1038    //Add neutral particles to jet
1039   TIter nextne(plNe) ; 
1040   particle = 0 ; 
1041   while ( (particle = dynamic_cast<TParticle*>(nextne())) ) {
1042     
1043     b0 = kFALSE;
1044     b1 = kFALSE;
1045
1046     //Particles in jet 
1047     SetJet(particle, b0, fJetCone, etal, phil) ;  
1048
1049     if(b0){
1050       new((*jetList)[n0++]) TParticle(*particle) ;
1051       particle->Momentum(lv);
1052       if(particle->Pt() > ptcut ){
1053         jet+=lv;
1054         ptjet+=particle->Pt();
1055       }
1056     }
1057
1058     //Background around (phi_gamma-pi, eta_leading)
1059     SetJet(particle, b1, fJetCone,etal, phig) ;
1060
1061     if(b1) { 
1062       new((*bkgList)[n1++]) TParticle(*particle) ;
1063       particle->Momentum(lv);
1064       if(particle->Pt() > ptcut ){
1065         bkg+=lv;
1066         ptbkg+=particle->Pt();    
1067       }  
1068     }
1069   }
1070   
1071   ptjet = jet.Pt();
1072   ptbkg = bkg.Pt();
1073
1074   if(ptjet > 0.) {
1075
1076     AliDebug(2,Form("Gamma   pt %f, Jet pt %f, Bkg pt %f",ptg,ptjet,ptbkg));
1077     
1078     //Fill histograms
1079     
1080     Double_t ratjet   = ptjet/ptg ;
1081     Double_t ratbkg  = ptbkg/ptg ;
1082
1083     dynamic_cast<TH2F*>
1084       (GetOutputContainer()->FindObject("JetRatio"+lastname))
1085       ->Fill(ptg,ratjet);        
1086     dynamic_cast<TH2F*>
1087       (GetOutputContainer()->FindObject("JetPt"+lastname))
1088       ->Fill(ptg,ptjet);
1089     
1090     dynamic_cast<TH2F*>
1091       (GetOutputContainer()->FindObject("BkgRatio"+lastname))
1092       ->Fill(ptg,ratbkg);
1093     
1094     dynamic_cast<TH2F*>
1095       (GetOutputContainer()->FindObject("BkgPt"+lastname))
1096       ->Fill(ptg,ptbkg);
1097
1098     //Jet selection
1099     Bool_t kSelect = kFALSE;
1100     if(fSelect == 0)
1101       kSelect = kTRUE; //Accept all jets, no restriction
1102     else if(fSelect == 1){
1103       //Selection with parametrized cuts
1104       if(IsJetSelected(ptg,ptjet))   kSelect = kTRUE;
1105     }
1106     else if(fSelect == 2){
1107       //Simple selection
1108       if(!AreJetsOnlyInCTS()){
1109         if((ratjet <  fJetRatioMaxCut) && (ratjet > fJetRatioMinCut )) kSelect = kTRUE;
1110       }
1111       else{
1112         if((ratjet <  fJetCTSRatioMaxCut) && (ratjet > fJetCTSRatioMinCut )) kSelect = kTRUE;
1113       }
1114     }
1115     else
1116       AliError("Jet selection option larger than 2, DON'T SELECT JETS");
1117     
1118     
1119      if(kSelect){
1120        AliDebug(1,Form("Jet Selected: pt %f ", ptjet)) ;
1121       
1122        FillJetHistos(jetList, ptg, ptl,"Jet",lastname);
1123        FillJetHistos(bkgList, ptg, ptl, "Bkg",lastname);
1124      }
1125   } //ptjet > 0
1126   
1127   jetList ->Delete();
1128   bkgList ->Delete();
1129   
1130 }
1131
1132 //___________________________________________________________________
1133 void AliAnaParticleJetLeadingCone::SetJet(TParticle * part, Bool_t & b, Float_t cone, 
1134                              Double_t eta, Double_t phi)
1135 {
1136
1137   //Check if the particle is inside the cone defined by the leading particle
1138   b = kFALSE;
1139   
1140   if(phi > TMath::TwoPi())
1141     phi-=TMath::TwoPi();
1142   if(phi < 0.)
1143     phi+=TMath::TwoPi();
1144   
1145   Double_t  rad = 10000 + cone;
1146   
1147   if(TMath::Abs(part->Phi()-phi) <= (TMath::TwoPi() - cone))
1148     rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1149                       TMath::Power(part->Phi()-phi,2));
1150   else{
1151     if(part->Phi()-phi > TMath::TwoPi() - cone)
1152       rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1153                         TMath::Power((part->Phi()-TMath::TwoPi())-phi,2));
1154     if(part->Phi()-phi < -(TMath::TwoPi() - cone))
1155       rad = TMath::Sqrt(TMath::Power(part->Eta()-eta,2)+
1156                         TMath::Power((part->Phi()+TMath::TwoPi())-phi,2));
1157   }
1158
1159   if(rad < cone )
1160     b = kTRUE;
1161   
1162 }
1163
1164 //____________________________________________________________________________
1165 void AliAnaParticleJetLeadingCone::FillJetHistos(TClonesArray * pl, Double_t ptg, Double_t ptl, TString type, TString lastname)
1166 {
1167   //Fill histograms wth jet fragmentation 
1168   //and number of selected jets and leading particles
1169   //and the background multiplicity
1170   TParticle * particle = 0 ;
1171   Int_t ipr = 0;
1172   Float_t  charge = 0;
1173   
1174   TIter next(pl) ; 
1175   while ( (particle = dynamic_cast<TParticle*>(next())) ) {
1176     ipr++ ;
1177     Double_t pt = particle->Pt();
1178     
1179     charge = TDatabasePDG::Instance()
1180       ->GetParticle(particle->GetPdgCode())->Charge();
1181     if(charge != 0){//Only jet Charged particles 
1182       dynamic_cast<TH2F*>
1183         (GetOutputContainer()->FindObject(type+"Fragment"+lastname))
1184         ->Fill(ptg,pt/ptg);
1185       dynamic_cast<TH2F*>
1186         (GetOutputContainer()->FindObject(type+"PtDist"+lastname))
1187         ->Fill(ptg,pt);
1188     }//charged
1189     
1190   }//while
1191   
1192   if(type == "Bkg")
1193     dynamic_cast<TH1F*>
1194       (GetOutputContainer()->FindObject("NBkg"+lastname))
1195       ->Fill(ipr);
1196   else{
1197     dynamic_cast<TH1F*>
1198       (GetOutputContainer()->FindObject("NJet"+lastname))->
1199       Fill(ptg);
1200     dynamic_cast<TH2F*>
1201       (GetOutputContainer()->FindObject("NLeading"+lastname))
1202       ->Fill(ptg,ptl);
1203   }
1204   
1205 }
1206
1207
1208 //____________________________________________________________________________
1209 Bool_t  AliAnaParticleJetLeadingConeCorrelation::SelectCluster(AliAODCaloCluster * calo, Double_t *vertex, TLorentzVector & mom, Int_t & pdg){
1210    //Select cluster depending on its pid and acceptance selections
1211    
1212    //Skip matched clusters with tracks
1213   if(calo->GetNTracksMatched() > 0) return kFALSE;
1214    
1215    //Check PID
1216    calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
1217    pdg = AliCaloPID::kPhoton;   
1218    if(IsCaloPIDOn()){
1219      //Get most probable PID, 2 options check PID weights (in MC this option is mandatory)
1220      //or redo PID, recommended option for EMCal.
1221      if(!IsCaloPIDRecalculationOn() || GetReader()->GetDataType() == AliCaloTrackReader::kMC )
1222        pdg = GetCaloPID()->GetPdg("EMCAL",calo->PID(),mom.E());//PID with weights
1223      else
1224        pdg = GetCaloPID()->GetPdg("EMCAL",mom,calo->GetM02(),0,0,0,0);//PID recalculated
1225      
1226      if(GetDebug() > 1) printf("PDG of identified particle %d\n",pdg);
1227      //If it does not pass pid, skip
1228      if(pdg != AliCaloPID::kPhoton || pdg != AliCaloPID::kPi0) {
1229        return kFALSE ;
1230    }
1231    
1232    //Check acceptance selection
1233    if(IsFidutialCutOn()){
1234      Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,"EMCAL") ;
1235      if(! in ) return kFALSE ;
1236    }
1237    
1238    if(GetDebug() > 1) printf("cluster selection cuts passed: pT %3.2f, pdg %d\n",mom.Pt(), pdg);
1239    
1240    return kTRUE;
1241    
1242  }