]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWG4/JetTasks/AliAnalysisTaskFragmentationFunction.cxx
97ea0bdaa9b4197b9c06a4d8d2678782e3512003
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskFragmentationFunction.cxx
1 /*************************************************************************
2  *                                                                       *
3  * Task for Fragmentation Function Analysis in PWG4 Jet Task Force Train *
4  *                                                                       *
5  *************************************************************************/
6
7
8 /**************************************************************************
9  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10  *                                                                        *
11  * Author: The ALICE Off-line Project.                                    *
12  * Contributors are mentioned in the code where appropriate.              *
13  *                                                                        *
14  * Permission to use, copy, modify and distribute this software and its   *
15  * documentation strictly for non-commercial purposes is hereby granted   *
16  * without fee, provided that the above copyright notice appears in all   *
17  * copies and that both the copyright notice and this permission notice   *
18  * appear in the supporting documentation. The authors make no claims     *
19  * about the suitability of this software for any purpose. It is          *
20  * provided "as is" without express or implied warranty.                  *
21  **************************************************************************/
22
23 /* $Id: */
24
25 #include "TList.h"
26 #include "TH1F.h"
27 #include "TH2F.h"
28 #include "TString.h"
29 #include "THnSparse.h"
30 #include "TProfile.h"
31 #include "TFile.h"
32 #include "TKey.h"
33
34 #include "AliAODInputHandler.h" 
35 #include "AliAODHandler.h" 
36 #include "AliESDEvent.h"
37 #include "AliAODMCParticle.h"
38 #include "AliAODJet.h"
39 #include "AliGenPythiaEventHeader.h"
40 #include "AliInputEventHandler.h"
41
42 #include "AliAnalysisHelperJetTasks.h"
43 #include "AliAnalysisManager.h"
44 #include "AliAnalysisTaskSE.h"
45
46 #include "AliAnalysisTaskFragmentationFunction.h"
47
48
49 ClassImp(AliAnalysisTaskFragmentationFunction)
50
51 //____________________________________________________________________________
52 AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction()
53    : AliAnalysisTaskSE()
54    ,fESD(0)
55    ,fAOD(0)
56    ,fMCEvent(0)
57    ,fBranchRecJets("jets")
58    ,fBranchGenJets("")
59    ,fTrackTypeGen(0)
60    ,fJetTypeGen(0)
61    ,fJetTypeRecEff(0)
62    ,fFilterMask(0)
63    ,fUsePhysicsSelection(kTRUE)
64    ,fTrackPtCut(0)
65    ,fTrackEtaMin(0)
66    ,fTrackEtaMax(0)
67    ,fTrackPhiMin(0)
68    ,fTrackPhiMax(0)
69    ,fJetPtCut(0)
70    ,fJetEtaMin(0)
71    ,fJetEtaMax(0)
72    ,fJetPhiMin(0)
73    ,fJetPhiMax(0)
74    ,fDiJetCut(0)
75    ,fDiJetDeltaPhiCut(0)
76    ,fDiJetPtFractionCut(0)
77    ,fDiJetCDFCut(0)
78    ,fDiJetKindBins(0)
79    ,fFFRadius(0)
80    ,fTracksRec(0)
81    ,fTracksRecCuts(0)
82    ,fTracksGen(0)
83    ,fTracksAODMCCharged(0)
84    ,fTracksRecQualityCuts(0)
85    ,fJetsRec(0)
86    ,fJetsRecCuts(0)
87    ,fJetsGen(0)
88    ,fJetsRecEff(0)
89    ,fQATrackHistosRec(0)
90    ,fQATrackHistosRecCuts(0)
91    ,fQATrackHistosGen(0)
92    ,fQAJetHistosRec(0)
93    ,fQAJetHistosRecCuts(0)
94    ,fQAJetHistosRecCutsLeading(0)
95    ,fQAJetHistosGen(0)
96    ,fQAJetHistosGenLeading(0)
97    ,fQAJetHistosRecEffLeading(0)
98    ,fFFHistosRecCuts(0)
99    ,fFFHistosRecLeading(0)
100    ,fFFHistosRecLeadingTrack(0)
101    ,fFFHistosGen(0)
102    ,fFFHistosGenLeading(0)
103    ,fFFHistosGenLeadingTrack(0)
104    ,fIJHistosRecCuts(0)
105    ,fIJHistosRecLeading(0)
106    ,fIJHistosRecLeadingTrack(0)
107    ,fIJHistosGen(0)
108    ,fIJHistosGenLeading(0)
109    ,fIJHistosGenLeadingTrack(0)
110    ,fFFDiJetHistosRecCuts(0)
111    ,fFFDiJetHistosRecLeading(0)
112    ,fFFDiJetHistosRecLeadingTrack(0)
113    ,fFFDiJetHistosGen(0)
114    ,fFFDiJetHistosGenLeading(0)
115    ,fFFDiJetHistosGenLeadingTrack(0)
116    ,fQADiJetHistosRecCuts(0)
117    ,fQADiJetHistosGen(0)
118    ,fQATrackHighPtThreshold(0)
119    ,fFFNBinsJetPt(0)    
120    ,fFFJetPtMin(0) 
121    ,fFFJetPtMax(0)
122    ,fFFNBinsPt(0)      
123    ,fFFPtMin(0)        
124    ,fFFPtMax(0)        
125    ,fFFNBinsXi(0)      
126    ,fFFXiMin(0)        
127    ,fFFXiMax(0)        
128    ,fFFNBinsZ(0)       
129    ,fFFZMin(0)         
130    ,fFFZMax(0)         
131    ,fQAJetNBinsPt(0)   
132    ,fQAJetPtMin(0)     
133    ,fQAJetPtMax(0)     
134    ,fQAJetNBinsEta(0)  
135    ,fQAJetEtaMin(0)    
136    ,fQAJetEtaMax(0)    
137    ,fQAJetNBinsPhi(0)  
138    ,fQAJetPhiMin(0)    
139    ,fQAJetPhiMax(0)    
140    ,fQATrackNBinsPt(0) 
141    ,fQATrackPtMin(0)   
142    ,fQATrackPtMax(0)   
143    ,fQATrackNBinsEta(0)
144    ,fQATrackEtaMin(0)  
145    ,fQATrackEtaMax(0)  
146    ,fQATrackNBinsPhi(0)
147    ,fQATrackPhiMin(0)  
148    ,fQATrackPhiMax(0)
149    ,fIJNBinsJetPt(0)
150    ,fIJJetPtMin(0)
151    ,fIJJetPtMax(0)
152    ,fIJNBinsPt(0)
153    ,fIJPtMin(0)
154    ,fIJPtMax(0)
155    ,fIJNBinsZ(0)
156    ,fIJZMin(0)
157    ,fIJZMax(0)
158    ,fIJNBinsCosTheta(0)
159    ,fIJCosThetaMin(0)
160    ,fIJCosThetaMax(0)
161    ,fIJNBinsTheta(0)
162    ,fIJThetaMin(0)
163    ,fIJThetaMax(0)
164    ,fIJNBinsJt(0)
165    ,fIJJtMin(0)
166    ,fIJJtMax(0)
167    ,fDiJetNBinsJetInvMass(0)
168    ,fDiJetJetInvMassMin(0)
169    ,fDiJetJetInvMassMax(0)
170    ,fDiJetNBinsJetPt(0)
171    ,fDiJetJetPtMin(0)
172    ,fDiJetJetPtMax(0)
173    ,fDiJetNBinsPt(0)
174    ,fDiJetPtMin(0)
175    ,fDiJetPtMax(0)
176    ,fDiJetNBinsXi(0)
177    ,fDiJetXiMin(0)
178    ,fDiJetXiMax(0)
179    ,fDiJetNBinsZ(0)
180    ,fDiJetZMin(0)
181    ,fDiJetZMax(0)
182    ,fQADiJetNBinsInvMass(0)
183    ,fQADiJetInvMassMin(0)
184    ,fQADiJetInvMassMax(0)
185    ,fQADiJetNBinsJetPt(0)
186    ,fQADiJetJetPtMin(0)
187    ,fQADiJetJetPtMax(0)
188    ,fQADiJetNBinsDeltaPhi(0)
189    ,fQADiJetDeltaPhiMin(0)
190    ,fQADiJetDeltaPhiMax(0)
191    ,fQADiJetNBinsDeltaEta(0)
192    ,fQADiJetDeltaEtaMin(0)
193    ,fQADiJetDeltaEtaMax(0)
194    ,fQADiJetNBinsDeltaPt(0)
195    ,fQADiJetDeltaPtMin(0)
196    ,fQADiJetDeltaPtMax(0)
197    ,fCommonHistList(0)
198    ,fh1EvtSelection(0)
199    ,fh1VertexNContributors(0)
200    ,fh1VertexZ(0)
201    ,fh1EvtMult(0)
202    ,fh1Xsec(0)
203    ,fh1Trials(0)
204    ,fh1PtHard(0)
205    ,fh1PtHardTrials(0)
206    ,fh1nRecJetsCuts(0)
207    ,fh1nGenJets(0)
208    ,fh1nRecEffJets(0)
209    ,fhnSingleTrackRecEffHisto(0)
210    ,fhnJetTrackRecEffHisto(0)
211 {
212    // default constructor
213 }
214
215 //__________________________________________________________________________________________
216 AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const char *name) 
217   : AliAnalysisTaskSE(name)
218   ,fESD(0)
219   ,fAOD(0)
220   ,fMCEvent(0)
221   ,fBranchRecJets("jets")
222   ,fBranchGenJets("")
223   ,fTrackTypeGen(0)
224   ,fJetTypeGen(0)
225   ,fJetTypeRecEff(0)
226   ,fFilterMask(0)
227   ,fUsePhysicsSelection(kTRUE)
228   ,fTrackPtCut(0)
229   ,fTrackEtaMin(0)
230   ,fTrackEtaMax(0)
231   ,fTrackPhiMin(0)
232   ,fTrackPhiMax(0)
233   ,fJetPtCut(0)
234   ,fJetEtaMin(0)
235   ,fJetEtaMax(0)
236   ,fJetPhiMin(0)
237   ,fJetPhiMax(0)
238   ,fDiJetCut(0)
239   ,fDiJetDeltaPhiCut(0)
240   ,fDiJetPtFractionCut(0)
241   ,fDiJetCDFCut(0)
242   ,fDiJetKindBins(0)
243   ,fFFRadius(0)
244   ,fTracksRec(0)
245   ,fTracksRecCuts(0)
246   ,fTracksGen(0)
247   ,fTracksAODMCCharged(0)
248   ,fTracksRecQualityCuts(0)
249   ,fJetsRec(0)
250   ,fJetsRecCuts(0)
251   ,fJetsGen(0)
252   ,fJetsRecEff(0)
253   ,fQATrackHistosRec(0)
254   ,fQATrackHistosRecCuts(0)
255   ,fQATrackHistosGen(0)
256   ,fQAJetHistosRec(0)
257   ,fQAJetHistosRecCuts(0)
258   ,fQAJetHistosRecCutsLeading(0)
259   ,fQAJetHistosGen(0)
260   ,fQAJetHistosGenLeading(0)
261   ,fQAJetHistosRecEffLeading(0)
262   ,fFFHistosRecCuts(0)
263   ,fFFHistosRecLeading(0)
264   ,fFFHistosRecLeadingTrack(0)
265   ,fFFHistosGen(0)
266   ,fFFHistosGenLeading(0)
267   ,fFFHistosGenLeadingTrack(0)
268   ,fIJHistosRecCuts(0)
269   ,fIJHistosRecLeading(0)
270   ,fIJHistosRecLeadingTrack(0)
271   ,fIJHistosGen(0)
272   ,fIJHistosGenLeading(0)
273   ,fIJHistosGenLeadingTrack(0)
274   ,fFFDiJetHistosRecCuts(0)
275   ,fFFDiJetHistosRecLeading(0)
276   ,fFFDiJetHistosRecLeadingTrack(0)
277   ,fFFDiJetHistosGen(0)
278   ,fFFDiJetHistosGenLeading(0)
279   ,fFFDiJetHistosGenLeadingTrack(0)
280   ,fQADiJetHistosRecCuts(0)
281   ,fQADiJetHistosGen(0)
282   ,fQATrackHighPtThreshold(0) 
283   ,fFFNBinsJetPt(0)    
284   ,fFFJetPtMin(0) 
285   ,fFFJetPtMax(0)
286   ,fFFNBinsPt(0)      
287   ,fFFPtMin(0)        
288   ,fFFPtMax(0)        
289   ,fFFNBinsXi(0)      
290   ,fFFXiMin(0)        
291   ,fFFXiMax(0)        
292   ,fFFNBinsZ(0)       
293   ,fFFZMin(0)         
294   ,fFFZMax(0)         
295   ,fQAJetNBinsPt(0)   
296   ,fQAJetPtMin(0)     
297   ,fQAJetPtMax(0)     
298   ,fQAJetNBinsEta(0)  
299   ,fQAJetEtaMin(0)    
300   ,fQAJetEtaMax(0)    
301   ,fQAJetNBinsPhi(0)  
302   ,fQAJetPhiMin(0)    
303   ,fQAJetPhiMax(0)    
304   ,fQATrackNBinsPt(0) 
305   ,fQATrackPtMin(0)   
306   ,fQATrackPtMax(0)   
307   ,fQATrackNBinsEta(0)
308   ,fQATrackEtaMin(0)  
309   ,fQATrackEtaMax(0)  
310   ,fQATrackNBinsPhi(0)
311   ,fQATrackPhiMin(0)  
312   ,fQATrackPhiMax(0)  
313   ,fIJNBinsJetPt(0)
314   ,fIJJetPtMin(0)
315   ,fIJJetPtMax(0)
316   ,fIJNBinsPt(0)
317   ,fIJPtMin(0)
318   ,fIJPtMax(0)
319   ,fIJNBinsZ(0)
320   ,fIJZMin(0)
321   ,fIJZMax(0)
322   ,fIJNBinsCosTheta(0)
323   ,fIJCosThetaMin(0)
324   ,fIJCosThetaMax(0)
325   ,fIJNBinsTheta(0)
326   ,fIJThetaMin(0)
327   ,fIJThetaMax(0)
328   ,fIJNBinsJt(0)
329   ,fIJJtMin(0)
330   ,fIJJtMax(0)
331   ,fDiJetNBinsJetInvMass(0)
332   ,fDiJetJetInvMassMin(0)
333   ,fDiJetJetInvMassMax(0)
334   ,fDiJetNBinsJetPt(0)
335   ,fDiJetJetPtMin(0)
336   ,fDiJetJetPtMax(0)
337   ,fDiJetNBinsPt(0)
338   ,fDiJetPtMin(0)
339   ,fDiJetPtMax(0)
340   ,fDiJetNBinsXi(0)
341   ,fDiJetXiMin(0)
342   ,fDiJetXiMax(0)
343   ,fDiJetNBinsZ(0)
344   ,fDiJetZMin(0)
345   ,fDiJetZMax(0)
346   ,fQADiJetNBinsInvMass(0)
347   ,fQADiJetInvMassMin(0)
348   ,fQADiJetInvMassMax(0)
349   ,fQADiJetNBinsJetPt(0)
350   ,fQADiJetJetPtMin(0)
351   ,fQADiJetJetPtMax(0)
352   ,fQADiJetNBinsDeltaPhi(0)
353   ,fQADiJetDeltaPhiMin(0)
354   ,fQADiJetDeltaPhiMax(0)
355   ,fQADiJetNBinsDeltaEta(0)
356   ,fQADiJetDeltaEtaMin(0)
357   ,fQADiJetDeltaEtaMax(0)
358   ,fQADiJetNBinsDeltaPt(0)
359   ,fQADiJetDeltaPtMin(0)
360   ,fQADiJetDeltaPtMax(0)
361   ,fCommonHistList(0)
362   ,fh1EvtSelection(0)
363   ,fh1VertexNContributors(0)
364   ,fh1VertexZ(0)
365   ,fh1EvtMult(0)
366   ,fh1Xsec(0)
367   ,fh1Trials(0)
368   ,fh1PtHard(0)
369   ,fh1PtHardTrials(0)
370   ,fh1nRecJetsCuts(0)
371   ,fh1nGenJets(0)
372   ,fh1nRecEffJets(0)
373   ,fhnSingleTrackRecEffHisto(0)
374   ,fhnJetTrackRecEffHisto(0)
375 {
376   // constructor
377   
378   DefineOutput(1,TList::Class());
379   
380
381 }
382
383 //__________________________________________________________________________________________________________________________
384 AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const  AliAnalysisTaskFragmentationFunction &copy)
385   : AliAnalysisTaskSE()
386   ,fESD(copy.fESD)
387   ,fAOD(copy.fAOD)
388   ,fMCEvent(copy.fMCEvent)
389   ,fBranchRecJets(copy.fBranchRecJets)
390   ,fBranchGenJets(copy.fBranchGenJets)
391   ,fTrackTypeGen(copy.fTrackTypeGen)
392   ,fJetTypeGen(copy.fJetTypeGen)
393   ,fJetTypeRecEff(copy.fJetTypeRecEff)
394   ,fFilterMask(copy.fFilterMask)
395   ,fUsePhysicsSelection(copy.fUsePhysicsSelection)
396   ,fTrackPtCut(copy.fTrackPtCut)
397   ,fTrackEtaMin(copy.fTrackEtaMin)
398   ,fTrackEtaMax(copy.fTrackEtaMax)
399   ,fTrackPhiMin(copy.fTrackPhiMin)
400   ,fTrackPhiMax(copy.fTrackPhiMax)
401   ,fJetPtCut(copy.fJetPtCut)
402   ,fJetEtaMin(copy.fJetEtaMin)
403   ,fJetEtaMax(copy.fJetEtaMax)
404   ,fJetPhiMin(copy.fJetPhiMin)
405   ,fJetPhiMax(copy.fJetPhiMax)
406   ,fDiJetCut(copy.fDiJetCut)
407   ,fDiJetDeltaPhiCut(copy.fDiJetDeltaPhiCut)
408   ,fDiJetPtFractionCut(copy.fDiJetPtFractionCut)
409   ,fDiJetCDFCut(copy.fDiJetCDFCut)
410   ,fDiJetKindBins(copy.fDiJetKindBins)
411   ,fFFRadius(copy.fFFRadius)
412   ,fTracksRec(copy.fTracksRec)
413   ,fTracksRecCuts(copy.fTracksRecCuts)
414   ,fTracksGen(copy.fTracksGen)
415   ,fTracksAODMCCharged(copy.fTracksAODMCCharged)
416   ,fTracksRecQualityCuts(copy.fTracksRecQualityCuts)
417   ,fJetsRec(copy.fJetsRec)
418   ,fJetsRecCuts(copy.fJetsRecCuts)
419   ,fJetsGen(copy.fJetsGen)
420   ,fJetsRecEff(copy.fJetsRecEff)
421   ,fQATrackHistosRec(copy.fQATrackHistosRec)
422   ,fQATrackHistosRecCuts(copy.fQATrackHistosRecCuts)
423   ,fQATrackHistosGen(copy.fQATrackHistosGen)
424   ,fQAJetHistosRec(copy.fQAJetHistosRec)
425   ,fQAJetHistosRecCuts(copy.fQAJetHistosRecCuts)
426   ,fQAJetHistosRecCutsLeading(copy.fQAJetHistosRecCutsLeading)
427   ,fQAJetHistosGen(copy.fQAJetHistosGen)
428   ,fQAJetHistosGenLeading(copy.fQAJetHistosGenLeading)
429   ,fQAJetHistosRecEffLeading(copy.fQAJetHistosRecEffLeading)
430   ,fFFHistosRecCuts(copy.fFFHistosRecCuts)
431   ,fFFHistosRecLeading(copy.fFFHistosRecLeading)
432   ,fFFHistosRecLeadingTrack(copy.fFFHistosRecLeadingTrack)
433   ,fFFHistosGen(copy.fFFHistosGen)
434   ,fFFHistosGenLeading(copy.fFFHistosGenLeading)
435   ,fFFHistosGenLeadingTrack(copy.fFFHistosGenLeadingTrack)
436   ,fIJHistosRecCuts(copy.fIJHistosRecCuts)
437   ,fIJHistosRecLeading(copy.fIJHistosRecLeading)
438   ,fIJHistosRecLeadingTrack(copy.fIJHistosRecLeadingTrack)
439   ,fIJHistosGen(copy.fIJHistosGen)
440   ,fIJHistosGenLeading(copy.fIJHistosGenLeading)
441   ,fIJHistosGenLeadingTrack(copy.fIJHistosGenLeadingTrack)
442   ,fFFDiJetHistosRecCuts(copy.fFFDiJetHistosRecCuts)
443   ,fFFDiJetHistosRecLeading(copy.fFFDiJetHistosRecLeading)
444   ,fFFDiJetHistosRecLeadingTrack(copy.fFFDiJetHistosRecLeadingTrack)
445   ,fFFDiJetHistosGen(copy.fFFDiJetHistosGen)
446   ,fFFDiJetHistosGenLeading(copy.fFFDiJetHistosGenLeading)
447   ,fFFDiJetHistosGenLeadingTrack(copy.fFFDiJetHistosGenLeadingTrack)
448   ,fQADiJetHistosRecCuts(copy.fQADiJetHistosRecCuts)
449   ,fQADiJetHistosGen(copy.fQADiJetHistosGen)
450   ,fQATrackHighPtThreshold(copy.fQATrackHighPtThreshold) 
451   ,fFFNBinsJetPt(copy.fFFNBinsJetPt)    
452   ,fFFJetPtMin(copy.fFFJetPtMin) 
453   ,fFFJetPtMax(copy.fFFJetPtMax)
454   ,fFFNBinsPt(copy.fFFNBinsPt)      
455   ,fFFPtMin(copy.fFFPtMin)        
456   ,fFFPtMax(copy.fFFPtMax)        
457   ,fFFNBinsXi(copy.fFFNBinsXi)      
458   ,fFFXiMin(copy.fFFXiMin)        
459   ,fFFXiMax(copy.fFFXiMax)        
460   ,fFFNBinsZ(copy.fFFNBinsZ)       
461   ,fFFZMin(copy.fFFZMin)         
462   ,fFFZMax(copy.fFFZMax)         
463   ,fQAJetNBinsPt(copy.fQAJetNBinsPt)   
464   ,fQAJetPtMin(copy.fQAJetPtMin)     
465   ,fQAJetPtMax(copy.fQAJetPtMax)     
466   ,fQAJetNBinsEta(copy.fQAJetNBinsEta)  
467   ,fQAJetEtaMin(copy.fQAJetEtaMin)    
468   ,fQAJetEtaMax(copy.fQAJetEtaMax)    
469   ,fQAJetNBinsPhi(copy.fQAJetNBinsPhi)  
470   ,fQAJetPhiMin(copy.fQAJetPhiMin)    
471   ,fQAJetPhiMax(copy.fQAJetPhiMax)    
472   ,fQATrackNBinsPt(copy.fQATrackNBinsPt) 
473   ,fQATrackPtMin(copy.fQATrackPtMin)   
474   ,fQATrackPtMax(copy.fQATrackPtMax)   
475   ,fQATrackNBinsEta(copy.fQATrackNBinsEta)
476   ,fQATrackEtaMin(copy.fQATrackEtaMin)  
477   ,fQATrackEtaMax(copy.fQATrackEtaMax)  
478   ,fQATrackNBinsPhi(copy.fQATrackNBinsPhi)
479   ,fQATrackPhiMin(copy.fQATrackPhiMin)  
480   ,fQATrackPhiMax(copy.fQATrackPhiMax)
481   ,fIJNBinsJetPt(copy.fIJNBinsJetPt)
482   ,fIJJetPtMin(copy.fIJJetPtMin)
483   ,fIJJetPtMax(copy.fIJJetPtMax)
484   ,fIJNBinsPt(copy.fIJNBinsPt)
485   ,fIJPtMin(copy.fIJPtMin)
486   ,fIJPtMax(copy.fIJPtMax)
487   ,fIJNBinsZ(copy.fIJNBinsZ)
488   ,fIJZMin(copy.fIJZMin)
489   ,fIJZMax(copy.fIJZMax)
490   ,fIJNBinsCosTheta(copy.fIJNBinsCosTheta)
491   ,fIJCosThetaMin(copy.fIJCosThetaMin)
492   ,fIJCosThetaMax(copy.fIJCosThetaMax)
493   ,fIJNBinsTheta(copy.fIJNBinsTheta)
494   ,fIJThetaMin(copy.fIJThetaMin)
495   ,fIJThetaMax(copy.fIJThetaMax)
496   ,fIJNBinsJt(copy.fIJNBinsJt)
497   ,fIJJtMin(copy.fIJJtMin)
498   ,fIJJtMax(copy.fIJJtMax)
499   ,fDiJetNBinsJetInvMass(copy.fDiJetNBinsJetInvMass)
500   ,fDiJetJetInvMassMin(copy.fDiJetJetInvMassMin)
501   ,fDiJetJetInvMassMax(copy.fDiJetJetInvMassMax)
502   ,fDiJetNBinsJetPt(copy.fDiJetNBinsJetPt)
503   ,fDiJetJetPtMin(copy.fDiJetJetPtMin)
504   ,fDiJetJetPtMax(copy.fDiJetJetPtMax)
505   ,fDiJetNBinsPt(copy.fDiJetNBinsPt)
506   ,fDiJetPtMin(copy.fDiJetPtMin)
507   ,fDiJetPtMax(copy.fDiJetPtMax)
508   ,fDiJetNBinsXi(copy.fDiJetNBinsXi)
509   ,fDiJetXiMin(copy.fDiJetXiMin)
510   ,fDiJetXiMax(copy.fDiJetXiMax)
511   ,fDiJetNBinsZ(copy.fDiJetNBinsZ)
512   ,fDiJetZMin(copy.fDiJetZMin)
513   ,fDiJetZMax(copy.fDiJetZMax)
514   ,fQADiJetNBinsInvMass(copy.fQADiJetNBinsInvMass)
515   ,fQADiJetInvMassMin(copy.fQADiJetInvMassMin)
516   ,fQADiJetInvMassMax(copy.fQADiJetInvMassMax)
517   ,fQADiJetNBinsJetPt(copy.fQADiJetNBinsJetPt)
518   ,fQADiJetJetPtMin(copy.fQADiJetJetPtMin)
519   ,fQADiJetJetPtMax(copy.fQADiJetJetPtMax)
520   ,fQADiJetNBinsDeltaPhi(copy.fQADiJetNBinsDeltaPhi)
521   ,fQADiJetDeltaPhiMin(copy.fQADiJetDeltaPhiMin)
522   ,fQADiJetDeltaPhiMax(copy.fQADiJetDeltaPhiMax)
523   ,fQADiJetNBinsDeltaEta(copy.fQADiJetNBinsDeltaEta)
524   ,fQADiJetDeltaEtaMin(copy.fQADiJetDeltaEtaMin)
525   ,fQADiJetDeltaEtaMax(copy.fQADiJetDeltaEtaMax)
526   ,fQADiJetNBinsDeltaPt(copy.fQADiJetNBinsDeltaPt)
527   ,fQADiJetDeltaPtMin(copy.fQADiJetDeltaPtMin)
528   ,fQADiJetDeltaPtMax(copy.fQADiJetDeltaPtMax)
529   ,fCommonHistList(copy.fCommonHistList)
530   ,fh1EvtSelection(copy.fh1EvtSelection)
531   ,fh1VertexNContributors(copy.fh1VertexNContributors)
532   ,fh1VertexZ(copy.fh1VertexZ)
533   ,fh1EvtMult(copy.fh1EvtMult)
534   ,fh1Xsec(copy.fh1Xsec)
535   ,fh1Trials(copy.fh1Trials)
536   ,fh1PtHard(copy.fh1PtHard)  
537   ,fh1PtHardTrials(copy.fh1PtHardTrials)  
538   ,fh1nRecJetsCuts(copy.fh1nRecJetsCuts)
539   ,fh1nGenJets(copy.fh1nGenJets)
540   ,fh1nRecEffJets(copy.fh1nRecEffJets)
541   ,fhnSingleTrackRecEffHisto(copy.fhnSingleTrackRecEffHisto)
542   ,fhnJetTrackRecEffHisto(copy.fhnJetTrackRecEffHisto)
543 {
544   // copy constructor
545
546 }
547
548 // _________________________________________________________________________________________________________________________________
549 AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::operator=(const AliAnalysisTaskFragmentationFunction& o)
550 {
551   // assignment
552   
553   if(this!=&o){
554
555     AliAnalysisTaskSE::operator=(o);
556     fESD                          = o.fESD;
557     fAOD                          = o.fAOD;
558     fMCEvent                      = o.fMCEvent;
559     fBranchRecJets                = o.fBranchRecJets;
560     fBranchGenJets                = o.fBranchGenJets;
561     fTrackTypeGen                 = o.fTrackTypeGen;
562     fJetTypeGen                   = o.fJetTypeGen;
563     fJetTypeRecEff                = o.fJetTypeRecEff;
564     fFilterMask                   = o.fFilterMask;
565     fUsePhysicsSelection          = o.fUsePhysicsSelection;
566     fTrackPtCut                   = o.fTrackPtCut;
567     fTrackEtaMin                  = o.fTrackEtaMin;
568     fTrackEtaMax                  = o.fTrackEtaMax;
569     fTrackPhiMin                  = o.fTrackPhiMin;
570     fTrackPhiMax                  = o.fTrackPhiMax;
571     fJetPtCut                     = o.fJetPtCut;
572     fJetEtaMin                    = o.fJetEtaMin;
573     fJetEtaMax                    = o.fJetEtaMax;
574     fJetPhiMin                    = o.fJetPhiMin;
575     fJetPhiMax                    = o.fJetPhiMin;
576     fDiJetCut                     = o.fDiJetCut;
577     fDiJetDeltaPhiCut             = o.fDiJetDeltaPhiCut;
578     fDiJetPtFractionCut           = o.fDiJetPtFractionCut;
579     fDiJetCDFCut                  = o.fDiJetCDFCut;
580     fDiJetKindBins                = o.fDiJetKindBins;
581     fFFRadius                     = o.fFFRadius;
582     fTracksRec                    = o.fTracksRec;
583     fTracksRecCuts                = o.fTracksRecCuts;
584     fTracksGen                    = o.fTracksGen;
585     fTracksAODMCCharged           = o.fTracksAODMCCharged;
586     fTracksRecQualityCuts         = o.fTracksRecQualityCuts;
587     fJetsRec                      = o.fJetsRec;
588     fJetsRecCuts                  = o.fJetsRecCuts;
589     fJetsGen                      = o.fJetsGen;
590     fJetsRecEff                   = o.fJetsRecEff;
591     fQATrackHistosRec             = o.fQATrackHistosRec;
592     fQATrackHistosRecCuts         = o.fQATrackHistosRecCuts;
593     fQATrackHistosGen             = o.fQATrackHistosGen;
594     fQAJetHistosRec               = o.fQAJetHistosRec;
595     fQAJetHistosRecCuts           = o.fQAJetHistosRecCuts;
596     fQAJetHistosRecCutsLeading    = o.fQAJetHistosRecCutsLeading;
597     fQAJetHistosGen               = o.fQAJetHistosGen;
598     fQAJetHistosGenLeading        = o.fQAJetHistosGenLeading;
599     fQAJetHistosRecEffLeading     = o.fQAJetHistosRecEffLeading;
600     fFFHistosRecCuts              = o.fFFHistosRecCuts;
601     fFFHistosRecLeading           = o.fFFHistosRecLeading;
602     fFFHistosRecLeadingTrack      = o.fFFHistosRecLeadingTrack;
603     fFFHistosGen                  = o.fFFHistosGen;
604     fFFHistosGenLeading           = o.fFFHistosGenLeading;
605     fFFHistosGenLeadingTrack      = o.fFFHistosGenLeadingTrack;
606     fIJHistosRecCuts              = o.fIJHistosRecCuts;
607     fIJHistosRecLeading           = o.fIJHistosRecLeading;
608     fIJHistosRecLeadingTrack      = o.fIJHistosRecLeadingTrack;
609     fIJHistosGen                  = o.fIJHistosGen;
610     fIJHistosGenLeading           = o.fIJHistosGenLeading;
611     fIJHistosGenLeadingTrack      = o.fIJHistosGenLeadingTrack;
612     fFFDiJetHistosRecCuts         = o.fFFDiJetHistosRecCuts;
613     fFFDiJetHistosRecLeading      = o.fFFDiJetHistosRecLeading;
614     fFFDiJetHistosRecLeadingTrack = o.fFFDiJetHistosRecLeadingTrack;
615     fFFDiJetHistosGen             = o.fFFDiJetHistosGen;
616     fFFDiJetHistosGenLeading      = o.fFFDiJetHistosGenLeading;
617     fFFDiJetHistosGenLeadingTrack = o.fFFDiJetHistosGenLeadingTrack;
618     fQADiJetHistosRecCuts         = o.fQADiJetHistosRecCuts;
619     fQADiJetHistosGen             = o.fQADiJetHistosGen;
620     fQATrackHighPtThreshold       = o.fQATrackHighPtThreshold; 
621     fFFNBinsJetPt                 = o.fFFNBinsJetPt;    
622     fFFJetPtMin                   = o.fFFJetPtMin; 
623     fFFJetPtMax                   = o.fFFJetPtMax;
624     fFFNBinsPt                    = o.fFFNBinsPt;      
625     fFFPtMin                      = o.fFFPtMin;        
626     fFFPtMax                      = o.fFFPtMax;        
627     fFFNBinsXi                    = o.fFFNBinsXi;      
628     fFFXiMin                      = o.fFFXiMin;        
629     fFFXiMax                      = o.fFFXiMax;        
630     fFFNBinsZ                     = o.fFFNBinsZ;       
631     fFFZMin                       = o.fFFZMin;         
632     fFFZMax                       = o.fFFZMax;         
633     fQAJetNBinsPt                 = o.fQAJetNBinsPt;   
634     fQAJetPtMin                   = o.fQAJetPtMin;     
635     fQAJetPtMax                   = o.fQAJetPtMax;     
636     fQAJetNBinsEta                = o.fQAJetNBinsEta;  
637     fQAJetEtaMin                  = o.fQAJetEtaMin;    
638     fQAJetEtaMax                  = o.fQAJetEtaMax;    
639     fQAJetNBinsPhi                = o.fQAJetNBinsPhi;  
640     fQAJetPhiMin                  = o.fQAJetPhiMin;    
641     fQAJetPhiMax                  = o.fQAJetPhiMax;    
642     fQATrackNBinsPt               = o.fQATrackNBinsPt; 
643     fQATrackPtMin                 = o.fQATrackPtMin;   
644     fQATrackPtMax                 = o.fQATrackPtMax;   
645     fQATrackNBinsEta              = o.fQATrackNBinsEta;
646     fQATrackEtaMin                = o.fQATrackEtaMin;  
647     fQATrackEtaMax                = o.fQATrackEtaMax;  
648     fQATrackNBinsPhi              = o.fQATrackNBinsPhi;
649     fQATrackPhiMin                = o.fQATrackPhiMin;  
650     fQATrackPhiMax                = o.fQATrackPhiMax;  
651     fIJNBinsJetPt                 = o.fIJNBinsJetPt;
652     fIJJetPtMin                   = o.fIJJetPtMin;
653     fIJJetPtMax                   = o.fIJJetPtMax;
654     fIJNBinsPt                    = o.fIJNBinsPt;
655     fIJPtMin                      = o.fIJPtMin;
656     fIJPtMax                      = o.fIJPtMax;
657     fIJNBinsZ                     = o.fIJNBinsZ;
658     fIJZMin                       = o.fIJZMin;
659     fIJZMax                       = o.fIJZMax;
660     fIJNBinsCosTheta              = o.fIJNBinsCosTheta;
661     fIJCosThetaMin                = o.fIJCosThetaMin;
662     fIJCosThetaMax                = o.fIJCosThetaMax;
663     fIJNBinsTheta                 = o.fIJNBinsTheta;
664     fIJThetaMin                   = o.fIJThetaMin;
665     fIJThetaMax                   = o.fIJThetaMax;
666     fIJNBinsJt                    = o.fIJNBinsJt;
667     fIJJtMin                      = o.fIJJtMin;
668     fIJJtMax                      = o.fIJJtMax;
669     fDiJetNBinsJetInvMass         = o.fDiJetNBinsJetInvMass;
670     fDiJetJetInvMassMin           = o.fDiJetJetInvMassMin;
671     fDiJetJetInvMassMax           = o.fDiJetJetInvMassMax;
672     fDiJetNBinsJetPt              = o.fDiJetNBinsJetPt;
673     fDiJetJetPtMin                = o.fDiJetJetPtMin;
674     fDiJetJetPtMax                = o.fDiJetJetPtMax;
675     fDiJetNBinsPt                 = o.fDiJetNBinsPt;
676     fDiJetPtMin                   = o.fDiJetPtMin;
677     fDiJetPtMax                   = o.fDiJetPtMax;
678     fDiJetNBinsXi                 = o.fDiJetNBinsXi;
679     fDiJetXiMin                   = o.fDiJetXiMin;
680     fDiJetXiMax                   = o.fDiJetXiMax;
681     fDiJetNBinsZ                  = o.fDiJetNBinsZ;
682     fDiJetZMin                    = o.fDiJetZMin;
683     fDiJetZMax                    = o.fDiJetZMax;
684     fQADiJetNBinsInvMass          = o.fQADiJetNBinsInvMass;
685     fQADiJetInvMassMin            = o.fQADiJetInvMassMin;
686     fQADiJetInvMassMax            = o.fQADiJetInvMassMax;
687     fQADiJetNBinsJetPt            = o.fQADiJetNBinsJetPt;
688     fQADiJetJetPtMin              = o.fQADiJetJetPtMin;
689     fQADiJetJetPtMax              = o.fQADiJetJetPtMax;
690     fQADiJetNBinsDeltaPhi         = o.fQADiJetNBinsDeltaPhi;
691     fQADiJetDeltaPhiMin           = o.fQADiJetDeltaPhiMin;
692     fQADiJetDeltaPhiMax           = o.fQADiJetDeltaPhiMax;
693     fQADiJetNBinsDeltaEta         = o.fQADiJetNBinsDeltaEta;
694     fQADiJetDeltaEtaMin           = o.fQADiJetDeltaEtaMin;
695     fQADiJetDeltaEtaMax           = o.fQADiJetDeltaEtaMax;
696     fQADiJetNBinsDeltaPt          = o.fQADiJetNBinsDeltaPt;
697     fQADiJetDeltaPtMin            = o.fQADiJetDeltaPtMin;
698     fQADiJetDeltaPtMax            = o.fQADiJetDeltaPtMax;
699     fCommonHistList               = o.fCommonHistList;
700     fh1EvtSelection               = o.fh1EvtSelection;
701     fh1VertexNContributors        = o.fh1VertexNContributors;
702     fh1VertexZ                    = o.fh1VertexZ;
703     fh1EvtMult                    = o.fh1EvtMult;
704     fh1Xsec                       = o.fh1Xsec;
705     fh1Trials                     = o.fh1Trials;
706     fh1PtHard                     = o.fh1PtHard;
707     fh1PtHardTrials               = o.fh1PtHardTrials;
708     fh1nRecJetsCuts               = o.fh1nRecJetsCuts;
709     fh1nGenJets                   = o.fh1nGenJets; 
710     fh1nRecEffJets                = o.fh1nRecEffJets;
711     fhnSingleTrackRecEffHisto     = o.fhnSingleTrackRecEffHisto;
712     fhnJetTrackRecEffHisto        = o.fhnJetTrackRecEffHisto;
713   }
714     
715   return *this;
716 }
717
718 //___________________________________________________________________________
719 AliAnalysisTaskFragmentationFunction::~AliAnalysisTaskFragmentationFunction()
720 {
721   // destructor
722   
723   if(fTracksRec)            delete fTracksRec;
724   if(fTracksRecCuts)        delete fTracksRecCuts;
725   if(fTracksGen)            delete fTracksGen;
726   if(fTracksAODMCCharged)   delete fTracksAODMCCharged;  
727   if(fTracksRecQualityCuts) delete fTracksRecQualityCuts; 
728   if(fJetsRec)              delete fJetsRec;
729   if(fJetsRecCuts)          delete fJetsRecCuts;
730   if(fJetsGen)              delete fJetsGen;
731   if(fJetsRecEff)           delete fJetsRecEff;
732
733   //  if(fDiJetBins)     delete fDiJetBins;
734
735 }
736
737 //______________________________________________________________________________________________________
738 AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AliFragFuncHistos(const char* name, 
739                                                          Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
740                                                          Int_t nPt, Float_t ptMin, Float_t ptMax,
741                                                          Int_t nXi, Float_t xiMin, Float_t xiMax,
742                                                          Int_t nZ , Float_t zMin , Float_t zMax )
743   : TObject()
744   ,fNBinsJetPt(nJetPt)
745   ,fJetPtMin(jetPtMin)
746   ,fJetPtMax(jetPtMax)
747   ,fNBinsPt(nPt) 
748   ,fPtMin(ptMin)   
749   ,fPtMax(ptMax)   
750   ,fNBinsXi(nXi) 
751   ,fXiMin(xiMin)   
752   ,fXiMax(xiMax)   
753   ,fNBinsZ(nZ)  
754   ,fZMin(zMin)    
755   ,fZMax(zMax)    
756   ,fh2TrackPt(0)
757   ,fh2Xi(0)
758   ,fh2Z(0)
759   ,fh1JetPt(0)
760   ,fName(name)
761 {
762   // default constructor
763
764 }
765
766 //___________________________________________________________________________
767 AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AliFragFuncHistos(const AliFragFuncHistos& copy)
768   : TObject()
769   ,fNBinsJetPt(copy.fNBinsJetPt)
770   ,fJetPtMin(copy.fJetPtMin)
771   ,fJetPtMax(copy.fJetPtMax)
772   ,fNBinsPt(copy.fNBinsPt) 
773   ,fPtMin(copy.fPtMin)   
774   ,fPtMax(copy.fPtMax)   
775   ,fNBinsXi(copy.fNBinsXi) 
776   ,fXiMin(copy.fXiMin)   
777   ,fXiMax(copy.fXiMax)   
778   ,fNBinsZ(copy.fNBinsZ)  
779   ,fZMin(copy.fZMin)    
780   ,fZMax(copy.fZMax)    
781   ,fh2TrackPt(copy.fh2TrackPt)
782   ,fh2Xi(copy.fh2Xi)
783   ,fh2Z(copy.fh2Z)
784   ,fh1JetPt(copy.fh1JetPt)
785   ,fName(copy.fName)
786 {
787   // copy constructor
788 }
789
790 //_______________________________________________________________________________________________________________________________________________________________
791 AliAnalysisTaskFragmentationFunction::AliFragFuncHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncHistos& o)
792 {
793   // assignment
794   
795   if(this!=&o){
796     TObject::operator=(o);
797     fNBinsJetPt = o.fNBinsJetPt;
798     fJetPtMin   = o.fJetPtMin;
799     fJetPtMax   = o.fJetPtMax;
800     fNBinsPt    = o.fNBinsPt; 
801     fPtMin      = o.fPtMin;   
802     fPtMax      = o.fPtMax;   
803     fNBinsXi    = o.fNBinsXi; 
804     fXiMin      = o.fXiMin;   
805     fXiMax      = o.fXiMax;   
806     fNBinsZ     = o.fNBinsZ;  
807     fZMin       = o.fZMin;    
808     fZMax       = o.fZMax;    
809     fh2TrackPt  = o.fh2TrackPt;
810     fh2Xi       = o.fh2Xi;
811     fh2Z        = o.fh2Z;
812     fh1JetPt    = o.fh1JetPt;
813     fName       = o.fName;
814   }
815     
816   return *this;
817 }
818
819 //_________________________________________________________
820 AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::~AliFragFuncHistos()
821 {
822   // destructor 
823
824   if(fh1JetPt)   delete fh1JetPt;
825   if(fh2TrackPt) delete fh2TrackPt;
826   if(fh2Xi)      delete fh2Xi;
827   if(fh2Z)       delete fh2Z;
828 }
829
830 //_________________________________________________________________
831 void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::DefineHistos()
832 {
833   // book FF histos
834
835   fh1JetPt   = new TH1F(Form("fh1FFJetPt%s", fName.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
836   fh2TrackPt = new TH2F(Form("fh2FFTrackPt%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax,fNBinsPt, fPtMin, fPtMax);
837   fh2Xi      = new TH2F(Form("fh2FFXi%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsXi, fXiMin, fXiMax);
838   fh2Z       = new TH2F(Form("fh2FFZ%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsZ, fZMin, fZMax);
839
840   AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{T} [GeV/c]", "entries"); 
841   AliAnalysisTaskFragmentationFunction::SetProperties(fh2TrackPt,"jet p_{T} [GeV/c]","p_{T} [GeV/c]","entries");
842   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Xi,"jet p_{T} [GeV/c]","#xi", "entries");
843   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Z,"jet p_{T} [GeV/c]","z","entries");
844 }
845
846 //_______________________________________________________________________________________________________________
847 void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::FillFF(Float_t trackPt, Float_t jetPt, Bool_t incrementJetPt)
848 {
849   // fill FF
850  
851   if(incrementJetPt) fh1JetPt->Fill(jetPt);    
852   fh2TrackPt->Fill(jetPt,trackPt);
853   
854   Double_t z = 0.;
855   if(jetPt>0) z = trackPt / jetPt;
856   Double_t xi = 0;
857   if(z>0) xi = TMath::Log(1/z);
858   
859   fh2Xi->Fill(jetPt,xi);
860   fh2Z->Fill(jetPt,z);
861 }
862
863 //_________________________________________________________________________________
864 void AliAnalysisTaskFragmentationFunction::AliFragFuncHistos::AddToOutput(TList* list) const
865 {
866   // add histos to list
867
868   list->Add(fh1JetPt);
869   
870   list->Add(fh2TrackPt);
871   list->Add(fh2Xi);
872   list->Add(fh2Z);
873 }
874
875 //_________________________________________________________________________________________________________
876 AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::AliFragFuncQAJetHistos(const char* name,
877                                                                Int_t nPt,   Float_t ptMin,   Float_t ptMax,
878                                                                Int_t nEta,  Float_t etaMin,  Float_t etaMax,
879                                                                Int_t nPhi,  Float_t phiMin,  Float_t phiMax)
880   : TObject()
881   ,fNBinsPt(nPt)
882   ,fPtMin(ptMin)
883   ,fPtMax(ptMax)
884   ,fNBinsEta(nEta)
885   ,fEtaMin(etaMin)
886   ,fEtaMax(etaMax)
887   ,fNBinsPhi(nPhi)
888   ,fPhiMin(phiMin)
889   ,fPhiMax(phiMax)
890   ,fh2EtaPhi(0)
891   ,fh1Pt(0)
892   ,fName(name)
893 {
894   // default constructor
895 }
896
897 //____________________________________________________________________________________
898 AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::AliFragFuncQAJetHistos(const AliFragFuncQAJetHistos& copy)
899   : TObject()
900   ,fNBinsPt(copy.fNBinsPt)
901   ,fPtMin(copy.fPtMin)
902   ,fPtMax(copy.fPtMax)
903   ,fNBinsEta(copy.fNBinsEta)
904   ,fEtaMin(copy.fEtaMin)
905   ,fEtaMax(copy.fEtaMax)
906   ,fNBinsPhi(copy.fNBinsPhi)
907   ,fPhiMin(copy.fPhiMin)
908   ,fPhiMax(copy.fPhiMax)
909   ,fh2EtaPhi(copy.fh2EtaPhi)
910   ,fh1Pt(copy.fh1Pt)
911   ,fName(copy.fName)
912 {
913   // copy constructor
914 }
915
916 //________________________________________________________________________________________________________________________________________________________________________
917 AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos& o)
918 {
919   // assignment
920   
921   if(this!=&o){
922     TObject::operator=(o);
923     fNBinsPt  = o.fNBinsPt;
924     fPtMin    = o.fPtMin;
925     fPtMax    = o.fPtMax;
926     fNBinsEta = o.fNBinsEta;
927     fEtaMin   = o.fEtaMin;
928     fEtaMax   = o.fEtaMax;
929     fNBinsPhi = o.fNBinsPhi;
930     fPhiMin   = o.fPhiMin;
931     fPhiMax   = o.fPhiMax;
932     fh2EtaPhi = o.fh2EtaPhi;
933     fh1Pt     = o.fh1Pt;
934     fName     = o.fName;
935   }
936   
937   return *this;
938 }
939
940 //______________________________________________________________
941 AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::~AliFragFuncQAJetHistos()
942 {
943   // destructor 
944   
945   if(fh2EtaPhi) delete fh2EtaPhi;
946   if(fh1Pt)     delete fh1Pt;
947 }
948
949 //____________________________________________________________________
950 void AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::DefineHistos()
951 {
952   // book jet QA histos
953
954   fh2EtaPhi  = new TH2F(Form("fh2JetQAEtaPhi%s", fName.Data()), Form("%s: #eta - #phi distribution", fName.Data()), fNBinsEta, fEtaMin, fEtaMax, fNBinsPhi, fPhiMin, fPhiMax);
955   fh1Pt      = new TH1F(Form("fh1JetQAPt%s", fName.Data()), Form("%s: p_{T} distribution", fName.Data()), fNBinsPt, fPtMin, fPtMax);
956         
957   AliAnalysisTaskFragmentationFunction::SetProperties(fh2EtaPhi, "#eta", "#phi"); 
958   AliAnalysisTaskFragmentationFunction::SetProperties(fh1Pt, "p_{T} [GeV/c]", "entries");
959 }
960
961 //____________________________________________________________________________________________________
962 void AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::FillJetQA(Float_t eta, Float_t phi, Float_t pt)
963 {
964   // fill jet QA histos 
965
966   fh2EtaPhi->Fill( eta, phi);
967   fh1Pt->Fill( pt );
968 }
969
970 //____________________________________________________________________________________
971 void AliAnalysisTaskFragmentationFunction::AliFragFuncQAJetHistos::AddToOutput(TList* list) const 
972 {
973   // add histos to list
974
975   list->Add(fh2EtaPhi);
976   list->Add(fh1Pt);
977 }
978
979 //___________________________________________________________________________________________________________
980 AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::AliFragFuncQATrackHistos(const char* name,
981                                                                    Int_t nPt, Float_t ptMin, Float_t ptMax,
982                                                                    Int_t nEta, Float_t etaMin, Float_t etaMax,
983                                                                    Int_t nPhi, Float_t phiMin, Float_t phiMax,
984                                                                    Float_t ptThresh) 
985   : TObject()
986   ,fNBinsPt(nPt)
987   ,fPtMin(ptMin)
988   ,fPtMax(ptMax)
989   ,fNBinsEta(nEta)
990   ,fEtaMin(etaMin)
991   ,fEtaMax(etaMax)
992   ,fNBinsPhi(nPhi)
993   ,fPhiMin(phiMin)
994   ,fPhiMax(phiMax)
995   ,fHighPtThreshold(ptThresh)
996   ,fh2EtaPhi(0)
997   ,fh1Pt(0)
998   ,fh2HighPtEtaPhi(0)
999   ,fName(name)
1000 {
1001   // default constructor
1002 }
1003
1004 //__________________________________________________________________________________________
1005 AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::AliFragFuncQATrackHistos(const AliFragFuncQATrackHistos& copy)
1006   : TObject()
1007   ,fNBinsPt(copy.fNBinsPt)
1008   ,fPtMin(copy.fPtMin)
1009   ,fPtMax(copy.fPtMax)
1010   ,fNBinsEta(copy.fNBinsEta)
1011   ,fEtaMin(copy.fEtaMin)
1012   ,fEtaMax(copy.fEtaMax)
1013   ,fNBinsPhi(copy.fNBinsPhi)
1014   ,fPhiMin(copy.fPhiMin)
1015   ,fPhiMax(copy.fPhiMax)
1016   ,fHighPtThreshold(copy.fHighPtThreshold)
1017   ,fh2EtaPhi(copy.fh2EtaPhi)
1018   ,fh1Pt(copy.fh1Pt)
1019   ,fh2HighPtEtaPhi(copy.fh2HighPtEtaPhi)
1020   ,fName(copy.fName)
1021 {
1022   // copy constructor
1023 }
1024
1025 // _____________________________________________________________________________________________________________________________________________________________________________
1026 AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos& o)
1027 {
1028   // assignment
1029   
1030   if(this!=&o){
1031     TObject::operator=(o);
1032     fNBinsPt         = o.fNBinsPt;
1033     fPtMin           = o.fPtMin;
1034     fPtMax           = o.fPtMax;
1035     fNBinsEta        = o.fNBinsEta;
1036     fEtaMin          = o.fEtaMin;
1037     fEtaMax          = o.fEtaMax;
1038     fNBinsPhi        = o.fNBinsPhi;
1039     fPhiMin          = o.fPhiMin;
1040     fPhiMax          = o.fPhiMax;
1041     fHighPtThreshold = o.fHighPtThreshold;
1042     fh2EtaPhi        = o.fh2EtaPhi;
1043     fh1Pt            = o.fh1Pt;
1044     fh2HighPtEtaPhi  = o.fh2HighPtEtaPhi;
1045     fName            = o.fName;
1046   }
1047   
1048   return *this;
1049 }
1050
1051 //___________________________________________________________________
1052 AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::~AliFragFuncQATrackHistos()
1053 {
1054   // destructor 
1055   
1056   if(fh2EtaPhi)       delete fh2EtaPhi;
1057   if(fh2HighPtEtaPhi) delete fh2HighPtEtaPhi;
1058   if(fh1Pt)           delete fh1Pt;
1059 }
1060
1061 //______________________________________________________________________
1062 void AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::DefineHistos()
1063 {
1064   // book track QA histos
1065
1066   fh2EtaPhi       = new TH2F(Form("fh2TrackQAEtaPhi%s", fName.Data()), Form("%s: #eta - #phi distribution", fName.Data()), fNBinsEta, fEtaMin, fEtaMax, fNBinsPhi, fPhiMin, fPhiMax);
1067   fh2HighPtEtaPhi = new TH2F(Form("fh2TrackQAHighPtEtaPhi%s", fName.Data()), Form("%s: #eta - #phi distribution for high-p_{T}", fName.Data()), fNBinsEta, fEtaMin, fEtaMax, fNBinsPhi, fPhiMin, fPhiMax);
1068   fh1Pt           = new TH1F(Form("fh1TrackQAPt%s", fName.Data()), Form("%s: p_{T} distribution", fName.Data()), fNBinsPt, fPtMin, fPtMax);
1069   
1070   AliAnalysisTaskFragmentationFunction::SetProperties(fh2EtaPhi, "#eta", "#phi"); 
1071   AliAnalysisTaskFragmentationFunction::SetProperties(fh2HighPtEtaPhi, "#eta", "#phi");
1072   AliAnalysisTaskFragmentationFunction::SetProperties(fh1Pt, "p_{T} [GeV/c]", "entries");
1073 }
1074
1075 //________________________________________________________________________________________________________
1076 void AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::FillTrackQA(Float_t eta, Float_t phi, Float_t pt)
1077 {
1078   // fill track QA histos
1079     
1080   fh2EtaPhi->Fill( eta, phi);
1081   if(pt > fHighPtThreshold) fh2HighPtEtaPhi->Fill( eta, phi);
1082   fh1Pt->Fill( pt );    
1083 }
1084
1085 //______________________________________________________________________________________
1086 void AliAnalysisTaskFragmentationFunction::AliFragFuncQATrackHistos::AddToOutput(TList* list) const
1087 {
1088   // add histos to list
1089
1090   list->Add(fh2EtaPhi);
1091   list->Add(fh2HighPtEtaPhi);
1092   list->Add(fh1Pt);
1093 }
1094
1095 //______________________________________________________________________________________________________
1096 AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::AliFragFuncIntraJetHistos(const char* name, 
1097                                                          Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
1098                                                          Int_t nPt, Float_t ptMin, Float_t ptMax,
1099                                                          Int_t nZ , Float_t zMin , Float_t zMax,
1100                                                          Int_t nCosTheta , Float_t costhetaMin , Float_t costhetaMax,
1101                                                          Int_t nTheta , Float_t thetaMin , Float_t thetaMax,
1102                                                          Int_t nJt , Float_t jtMin , Float_t jtMax)
1103   : TObject()
1104   ,fNBinsJetPt(nJetPt)
1105   ,fJetPtMin(jetPtMin)
1106   ,fJetPtMax(jetPtMax)
1107   ,fNBinsPt(nPt) 
1108   ,fPtMin(ptMin)   
1109   ,fPtMax(ptMax)   
1110   ,fNBinsZ(nZ) 
1111   ,fZMin(zMin)   
1112   ,fZMax(zMax)   
1113   ,fNBinsJt(nJt)
1114   ,fJtMin(jtMin)
1115   ,fJtMax(jtMax)
1116   ,fNBinsTheta(nTheta)
1117   ,fThetaMin(thetaMin)
1118   ,fThetaMax(thetaMax)
1119   ,fNBinsCosTheta(nCosTheta)
1120   ,fCosThetaMin(costhetaMin)
1121   ,fCosThetaMax(costhetaMax)
1122   ,fh2Theta(0)
1123   ,fh2CosTheta(0)
1124   ,fh2Jt(0)
1125   ,fh2PtvsZ(0)
1126   ,fhnIntraJet(0)
1127   ,fnDim(6)
1128   ,fName(name)
1129 {
1130   // default constructor
1131
1132 }
1133
1134 //___________________________________________________________________________
1135 AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::AliFragFuncIntraJetHistos(const AliFragFuncIntraJetHistos& copy)
1136   : TObject()
1137   ,fNBinsJetPt(copy.fNBinsJetPt)
1138   ,fJetPtMin(copy.fJetPtMin)
1139   ,fJetPtMax(copy.fJetPtMax)
1140   ,fNBinsPt(copy.fNBinsPt) 
1141   ,fPtMin(copy.fPtMin)   
1142   ,fPtMax(copy.fPtMax)   
1143   ,fNBinsZ(copy.fNBinsZ) 
1144   ,fZMin(copy.fZMin)   
1145   ,fZMax(copy.fZMax)   
1146   ,fNBinsJt(copy.fNBinsJt)
1147   ,fJtMin(copy.fJtMin)
1148   ,fJtMax(copy.fJtMax)
1149   ,fNBinsTheta(copy.fNBinsTheta)
1150   ,fThetaMin(copy.fThetaMin)
1151   ,fThetaMax(copy.fThetaMax)
1152   ,fNBinsCosTheta(copy.fNBinsCosTheta)
1153   ,fCosThetaMin(copy.fCosThetaMin)
1154   ,fCosThetaMax(copy.fCosThetaMax)
1155   ,fh2Theta(copy.fh2Theta)
1156   ,fh2CosTheta(copy.fh2CosTheta)
1157   ,fh2Jt(copy.fh2Jt)
1158   ,fh2PtvsZ(copy.fh2PtvsZ)
1159   ,fhnIntraJet(copy.fhnIntraJet)
1160   ,fnDim(copy.fnDim)
1161   ,fName(copy.fName)
1162 {
1163   // copy constructor
1164 }
1165
1166 //_______________________________________________________________________________________________________________________________________________________________
1167 AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos& o)
1168 {
1169   // assignment
1170   
1171   if(this!=&o){
1172     TObject::operator=(o);
1173     fNBinsJetPt       = o.fNBinsJetPt;
1174     fJetPtMin         = o.fJetPtMin;
1175     fJetPtMax         = o.fJetPtMax;
1176     fNBinsPt          = o.fNBinsPt; 
1177     fPtMin            = o.fPtMin;   
1178     fPtMax            = o.fPtMax;   
1179     fNBinsZ           = o.fNBinsZ; 
1180     fZMin             = o.fZMin;   
1181     fZMax             = o.fZMax;   
1182     fNBinsJt          = o.fNBinsJt;
1183     fJtMin            = o.fJtMin;
1184     fJtMax            = o.fJtMax;
1185     fNBinsTheta       = o.fNBinsTheta;
1186     fThetaMin         = o.fThetaMin;
1187     fThetaMax         = o.fThetaMax;
1188     fNBinsCosTheta    = o.fNBinsCosTheta;
1189     fCosThetaMin      = o.fCosThetaMin;
1190     fCosThetaMax      = o.fCosThetaMax;
1191     fh2Theta          = o.fh2Theta;
1192     fh2CosTheta       = o.fh2CosTheta;
1193     fh2Jt             = o.fh2Jt;
1194     fh2PtvsZ          = o.fh2PtvsZ;
1195     fhnIntraJet       = o.fhnIntraJet;
1196     fnDim             = o.fnDim;
1197     fName             = o.fName;
1198   }
1199     
1200   return *this;
1201 }
1202
1203 //_________________________________________________________
1204 AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::~AliFragFuncIntraJetHistos()
1205 {
1206   // destructor 
1207
1208
1209   if(fh2Theta)          delete fh2Theta;
1210   if(fh2CosTheta)       delete fh2CosTheta;
1211   if(fh2Jt)             delete fh2Jt;
1212   if(fh2PtvsZ)          delete fh2PtvsZ;
1213   if(fhnIntraJet)       delete fhnIntraJet;
1214
1215 }
1216
1217 //_________________________________________________________________
1218 void AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::DefineHistos()
1219 {
1220   // book FF histos
1221
1222   fh2Theta    = new TH2F(Form("fh2IJTheta%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax,fNBinsTheta, fThetaMin, fThetaMax);
1223   fh2CosTheta = new TH2F(Form("fh2IJcosTheta%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax,fNBinsCosTheta, fCosThetaMin, fCosThetaMax);
1224   fh2Jt       = new TH2F(Form("fh2IJJt%s",fName.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsJt, fJtMin, fJtMax);
1225
1226   // Create 3D histograms
1227   Int_t    *iBin = new Int_t[fnDim];
1228   Double_t *min  = new Double_t[fnDim];
1229   Double_t *max  = new Double_t[fnDim];
1230
1231   iBin[0] = fNBinsJetPt; iBin[1] = fNBinsTheta; iBin[2] = fNBinsCosTheta; iBin[3] = fNBinsJt; iBin[4] = fNBinsZ; iBin[5] = fNBinsPt;
1232   min[0]  = fJetPtMin; min[1] = fThetaMin; min[2] = fCosThetaMin; min[3] = fJtMin; min[4] = fZMin; min[5] = fPtMin; 
1233   max[0]  = fJetPtMax; max[1] = fThetaMax; max[2] = fCosThetaMax; max[3] = fJtMax; max[4] = fZMax; max[5] = fPtMax;
1234
1235   const char* title = Form("fhnIntraJetPart%s",fName.Data());
1236   const char* comment = "THnSparseF p_{T} jet [GeV/c] : #Theta : cos(#Theta) : j_{T} : Z : p_{T} part [GeV/c]";
1237   fhnIntraJet = new THnSparseF(title,comment,fnDim,iBin,min,max);
1238
1239   const char** axisTitle = new const char*[fnDim];
1240   axisTitle[0] = "p_{T}^{jet} [GeV/c]";
1241   axisTitle[1] = "#Theta";
1242   axisTitle[2] = "Cos(#Theta)";
1243   axisTitle[3] = "j_{T} [GeV]";
1244   axisTitle[4] = "z = p_{T}^{had}/p_{T}^{jet}";
1245   axisTitle[5] = "p_{T}^{had} [GeV/c]";
1246   
1247   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Theta,"jet p_{T} [GeV/c]","#Theta","entries");
1248   AliAnalysisTaskFragmentationFunction::SetProperties(fh2CosTheta,"jet p_{T} [GeV/c]","cos(#Theta)", "entries");
1249   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Jt,"jet p_{T} [GeV/c]","j_{T}","entries");
1250   AliAnalysisTaskFragmentationFunction::SetProperties(fhnIntraJet,fnDim,axisTitle);
1251   delete[] iBin;
1252   delete[] min;
1253   delete[] max;
1254   delete[] axisTitle;
1255 }
1256
1257 //_______________________________________________________________________________________________________________
1258 void AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::FillIntraJet(TLorentzVector* trackV, TLorentzVector* jetV)
1259 {
1260   // fill IntraJet histos
1261  
1262   Float_t cosTheta = 0.; Float_t theta = 0.; 
1263   Float_t jt = 0.; Float_t z = 0.; 
1264   // For Theta distribution
1265   Float_t pxT  = trackV->Px();
1266   Float_t pyT  = trackV->Py();
1267   Float_t pzT  = trackV->Pz();
1268   Float_t ptT  = trackV->Pt();
1269   Float_t pT   = trackV->P();
1270   Float_t etaT = trackV->Eta();
1271   Float_t phiT = trackV->Phi(); // Check the value returned
1272   Float_t pxJ = jetV->Px();
1273   Float_t pyJ = jetV->Py();
1274   Float_t pzJ = jetV->Pz();
1275   Float_t ptJ = jetV->Pt();
1276   Float_t pJ  = jetV->P();
1277
1278   // Compute z
1279   z = (Float_t)(ptT/ptJ);
1280
1281   // Compute theta
1282   cosTheta = (pxT*pxJ+pyT*pyJ+pzT*pzJ)/(pT*pJ);
1283   theta = TMath::ACos(cosTheta);
1284
1285   // Compute jt
1286   TVector3 trackP; TVector3 jetP;
1287   jetP[0] = pxJ;
1288   jetP[1] = pyJ;
1289   jetP[2] = pzJ;
1290   trackP.SetPtEtaPhi(ptT,etaT,phiT);
1291   jt = TMath::Sin(trackP.Angle(jetP))*trackP.Mag();
1292
1293   // Fill histos and THnSparse
1294   fh2CosTheta->Fill(ptJ,cosTheta);
1295   fh2Theta->Fill(ptJ,theta);
1296   fh2Jt->Fill(ptJ,jt);
1297
1298   // Fill THnSparse
1299   Double_t *content = new Double_t[fnDim];
1300   content[0]= ptJ; content[1] = theta; content[2] = cosTheta; content[3] = jt; content[4] = z; content[5] = ptT; 
1301
1302   fhnIntraJet->Fill(content);
1303
1304   delete content;
1305
1306 }
1307
1308 //______________________________________________________________________________________________________
1309 AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::AliFragFuncDiJetHistos(const char* name, Int_t kindSlices,
1310                                                          Int_t nJetInvMass, Float_t jetInvMassMin, Float_t jetInvMassMax,  
1311                                                          Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
1312                                                          Int_t nPt, Float_t ptMin, Float_t ptMax,
1313                                                          Int_t nXi, Float_t xiMin, Float_t xiMax,
1314                                                          Int_t nZ , Float_t zMin , Float_t zMax)
1315   : TObject()
1316   ,fKindSlices(kindSlices)
1317   ,fNBinsJetInvMass(nJetInvMass)
1318   ,fJetInvMassMin(jetInvMassMin)
1319   ,fJetInvMassMax(jetInvMassMax)
1320   ,fNBinsJetPt(nJetPt)
1321   ,fJetPtMin(jetPtMin)
1322   ,fJetPtMax(jetPtMax)
1323   ,fNBinsPt(nPt) 
1324   ,fPtMin(ptMin)   
1325   ,fPtMax(ptMax)   
1326   ,fNBinsXi(nXi) 
1327   ,fXiMin(xiMin)   
1328   ,fXiMax(xiMax)   
1329   ,fNBinsZ(nZ)  
1330   ,fZMin(zMin)    
1331   ,fZMax(zMax)
1332   ,fh2TrackPtJet1(0)
1333   ,fh2TrackPtJet2(0)
1334   ,fh2TrackPtJet(0)
1335   ,fh1Jet1Pt(0)
1336   ,fh1Jet2Pt(0)
1337   ,fh1JetPt(0)
1338   ,fh2Xi1(0)
1339   ,fh2Xi2(0)
1340   ,fh2Xi(0)
1341   ,fh2Z1(0)
1342   ,fh2Z2(0)
1343   ,fh2Z(0)
1344   ,fh2Pt1(0)
1345   ,fh2Pt2(0)
1346   ,fh2Pt(0)
1347   ,fName(name)
1348 {
1349   // default constructor
1350
1351 }
1352
1353 //______________________________________________________________________________________________________
1354 AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::AliFragFuncDiJetHistos(const AliFragFuncDiJetHistos& copy)
1355   : TObject()
1356   ,fKindSlices(copy.fKindSlices)
1357   ,fNBinsJetInvMass(copy.fNBinsJetInvMass)
1358   ,fJetInvMassMin(copy.fJetInvMassMin)
1359   ,fJetInvMassMax(copy.fJetInvMassMax)
1360   ,fNBinsJetPt(copy.fNBinsJetPt)
1361   ,fJetPtMin(copy.fJetPtMin)
1362   ,fJetPtMax(copy.fJetPtMax)
1363   ,fNBinsPt(copy.fNBinsPt) 
1364   ,fPtMin(copy.fPtMin)   
1365   ,fPtMax(copy.fPtMax)   
1366   ,fNBinsXi(copy.fNBinsXi) 
1367   ,fXiMin(copy.fXiMin)   
1368   ,fXiMax(copy.fXiMax)   
1369   ,fNBinsZ(copy.fNBinsZ)  
1370   ,fZMin(copy.fZMin)    
1371   ,fZMax(copy.fZMax)
1372   ,fh2TrackPtJet1(copy.fh2TrackPtJet1)
1373   ,fh2TrackPtJet2(copy.fh2TrackPtJet2)
1374   ,fh2TrackPtJet(copy.fh2TrackPtJet)
1375   ,fh1Jet1Pt(copy.fh1Jet1Pt)
1376   ,fh1Jet2Pt(copy.fh1Jet2Pt)
1377   ,fh1JetPt(copy.fh1JetPt)
1378   ,fh2Xi1(copy.fh2Xi1)
1379   ,fh2Xi2(copy.fh2Xi2)
1380   ,fh2Xi(copy.fh2Xi2)
1381   ,fh2Z1(copy.fh2Z1)
1382   ,fh2Z2(copy.fh2Z2)
1383   ,fh2Z(copy.fh2Z)
1384   ,fh2Pt1(copy.fh2Pt1)
1385   ,fh2Pt2(copy.fh2Pt2)
1386   ,fh2Pt(copy.fh2Pt)
1387   ,fName(copy.fName)
1388 {
1389   // default constructor
1390
1391 }
1392
1393 //_______________________________________________________________________________________________________________________________________________________________
1394 AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos& o)
1395 {
1396   // assignment
1397   
1398   if(this!=&o){
1399     TObject::operator=(o);
1400     fKindSlices      = o.fKindSlices;
1401     fNBinsJetInvMass = o.fNBinsJetInvMass;
1402     fJetInvMassMin   = o.fJetInvMassMin;
1403     fJetInvMassMax   = o.fJetInvMassMax;
1404     fNBinsJetPt      = o.fNBinsJetPt;
1405     fJetPtMin        = o.fJetPtMin;
1406     fJetPtMax        = o.fJetPtMax;
1407     fNBinsPt         = o.fNBinsPt; 
1408     fPtMin           = o.fPtMin;   
1409     fPtMax           = o.fPtMax;   
1410     fNBinsXi         = o.fNBinsXi; 
1411     fXiMin           = o.fXiMin;   
1412     fXiMax           = o.fXiMax;   
1413     fNBinsZ          = o.fNBinsZ;  
1414     fZMin            = o.fZMin;    
1415     fZMax            = o.fZMax;   
1416     fh2TrackPtJet1   = o.fh2TrackPtJet1;
1417     fh2TrackPtJet2   = o.fh2TrackPtJet2;
1418     fh2TrackPtJet    = o.fh2TrackPtJet;
1419     fh1Jet1Pt        = o.fh1Jet1Pt;
1420     fh1Jet2Pt        = o.fh1Jet2Pt;
1421     fh1JetPt         = o.fh1JetPt;
1422     fh2Xi1           = o.fh2Xi1;
1423     fh2Xi2           = o.fh2Xi2;
1424     fh2Xi            = o.fh2Xi;
1425     fh2Z1            = o.fh2Z1;
1426     fh2Z2            = o.fh2Z2;
1427     fh2Z             = o.fh2Z;
1428     fh2Pt1           = o.fh2Pt1;
1429     fh2Pt2           = o.fh2Pt2;
1430     fh2Pt            = o.fh2Pt;
1431     fName            = o.fName;
1432   }
1433     
1434   return *this;
1435 }
1436
1437 //_________________________________________________________
1438 AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::~AliFragFuncDiJetHistos()
1439 {
1440   // destructor 
1441
1442   if(fh2TrackPtJet1) delete fh2TrackPtJet1;
1443   if(fh2TrackPtJet2) delete fh2TrackPtJet2;
1444   if(fh2TrackPtJet ) delete fh2TrackPtJet;
1445   if(fh1Jet1Pt)      delete fh1Jet1Pt;
1446   if(fh1Jet2Pt)      delete fh1Jet2Pt;
1447   if(fh1JetPt)       delete fh1JetPt;
1448   if(fh2Xi1)         delete fh2Xi1;
1449   if(fh2Xi2)         delete fh2Xi2;
1450   if(fh2Xi)          delete fh2Xi;
1451   if(fh2Z1)          delete fh2Z1;
1452   if(fh2Z2)          delete fh2Z2;
1453   if(fh2Z)           delete fh2Z;
1454   if(fh2Pt1)         delete fh2Pt1;
1455   if(fh2Pt2)         delete fh2Pt2;
1456   if(fh2Pt)          delete fh2Pt;
1457 }
1458
1459 //________________________________________________________________________
1460 void AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::DefineDiJetHistos()
1461 {
1462
1463   Int_t nBins = 0;
1464   Double_t min = 0.;
1465   Double_t max = 0.;
1466   const char *xaxis = "";
1467   if(fKindSlices == 1)
1468     {
1469       nBins = fNBinsJetInvMass;
1470       min   = fJetInvMassMin;
1471       max   = fJetInvMassMax;
1472       xaxis = "M_{JJ} [GeV]";
1473     }
1474   if(fKindSlices == 2 || fKindSlices == 3)
1475     {
1476       nBins = fNBinsJetPt;
1477       min   = fJetPtMin;
1478       max   = fJetPtMax;
1479       if(fKindSlices == 2) xaxis = "E_{Tmean} [GeV]";
1480       if(fKindSlices == 3) xaxis ="leading jet p_{T} [GeV/c]";
1481     }
1482   
1483   fh1Jet1Pt      = new TH1F(Form("fh1DJJet1Pt%s", fName.Data()), "", fNBinsJetPt, fJetPtMin, fJetPtMax);
1484   fh1Jet2Pt      = new TH1F(Form("fh1DJJet2Pt%s", fName.Data()), "", fNBinsJetPt, fJetPtMin, fJetPtMax);
1485   fh1JetPt       = new TH1F(Form("fh1DJJetPt%s",  fName.Data()), "", fNBinsJetPt, fJetPtMin, fJetPtMax);
1486   
1487   fh2TrackPtJet1 = new TH2F(Form("fh2DJTrackPtJet1%s", fName.Data()), "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
1488   fh2TrackPtJet2 = new TH2F(Form("fh2DJTrackPtJet2%s", fName.Data()), "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
1489   fh2TrackPtJet  = new TH2F(Form("fh2DJTrackPtJet%s", fName.Data()),  "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
1490   
1491   fh2Xi1         = new TH2F(Form("fh2DJXi1%s", fName.Data()), "",nBins, min, max, fNBinsXi, fXiMin, fXiMax);
1492   fh2Xi2         = new TH2F(Form("fh2DJXi2%s", fName.Data()), "",nBins, min, max, fNBinsXi, fXiMin, fXiMax);
1493   fh2Xi          = new TH2F(Form("fh2DJXi%s", fName.Data()),  "",nBins, min, max, fNBinsXi, fXiMin, fXiMax);
1494   
1495   fh2Z1          = new TH2F(Form("fh2DJZ1%s", fName.Data()), "",nBins, min, max, fNBinsZ, fZMin, fZMax);
1496   fh2Z2          = new TH2F(Form("fh2DJZ2%s", fName.Data()), "",nBins, min, max, fNBinsZ, fZMin, fZMax);
1497   fh2Z           = new TH2F(Form("fh2DJZ%s", fName.Data()),  "",nBins, min, max, fNBinsZ, fZMin, fZMax);
1498   
1499   fh2Pt1         = new TH2F(Form("fh2DJPt1%s", fName.Data()), "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
1500   fh2Pt2         = new TH2F(Form("fh2DJPt2%s", fName.Data()), "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
1501   fh2Pt          = new TH2F(Form("fh2DJPtZ%s", fName.Data()),  "",nBins, min, max, fNBinsPt, fPtMin, fPtMax);
1502       
1503   AliAnalysisTaskFragmentationFunction::SetProperties(fh1Jet1Pt, "p_{T} [GeV/c]", "entries");
1504   AliAnalysisTaskFragmentationFunction::SetProperties(fh1Jet2Pt, "p_{T} [GeV/c]", "entries");
1505   AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{T} [GeV/c]", "entries");
1506
1507   AliAnalysisTaskFragmentationFunction::SetProperties(fh2TrackPtJet1, xaxis, "p_{T} [GeV/c]", "Entries");
1508   AliAnalysisTaskFragmentationFunction::SetProperties(fh2TrackPtJet2, xaxis, "p_{T} [GeV/c]", "Entries");
1509   AliAnalysisTaskFragmentationFunction::SetProperties(fh2TrackPtJet, xaxis, "p_{T} [GeV/c]", "Entries");
1510   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Xi1, xaxis, "#xi", "Entries");
1511   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Xi2, xaxis, "#xi", "Entries");
1512   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Xi, xaxis, "#xi", "Entries");
1513   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Z1, xaxis, "z", "Entries");
1514   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Z2, xaxis, "z", "Entries");
1515   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Z, xaxis, "z", "Entries");
1516   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Pt1, xaxis, "p_{T} [GeV/c]", "Entries");
1517   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Pt2, xaxis, "p_{T} [GeV/c]", "Entries");
1518   AliAnalysisTaskFragmentationFunction::SetProperties(fh2Pt, xaxis, "p_{T} [GeV/c]", "Entries");
1519   
1520 }
1521
1522 //________________________________________________________________________
1523 void AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::FillDiJetFF(Int_t jetType, Float_t trackPt, Float_t jetPt, Double_t jetBin, Bool_t incrementJetPt)
1524 {
1525   if(jetType == 0)
1526     {
1527       if(incrementJetPt) fh1JetPt->Fill(jetPt);  
1528       
1529       fh2TrackPtJet->Fill(jetBin, trackPt);
1530       
1531       Double_t z = trackPt / jetPt;
1532       Double_t xi = 0;
1533       if(z!=0) xi = TMath::Log(1/z);
1534       
1535       fh2Xi->Fill(jetBin, xi);
1536       fh2Z->Fill(jetBin, z);
1537     }
1538   if(jetType == 1)
1539     {
1540       if(incrementJetPt) fh1Jet1Pt->Fill(jetPt);
1541       
1542       fh2TrackPtJet1->Fill(jetBin, trackPt);
1543       
1544       Double_t z = trackPt / jetPt;
1545       Double_t xi = 0;
1546       if(z!=0) xi = TMath::Log(1/z);
1547       
1548       fh2Xi1->Fill(jetBin, xi);
1549       fh2Z1->Fill(jetBin, z);
1550     }
1551   if(jetType == 2)
1552     {
1553       if(incrementJetPt) fh1Jet2Pt->Fill(jetPt);
1554       
1555       fh2TrackPtJet2->Fill(jetBin, trackPt);
1556       
1557       Double_t z = trackPt / jetPt;
1558       Double_t xi = 0;
1559       if(z!=0) xi = TMath::Log(1/z);
1560       
1561       fh2Xi2->Fill(jetBin, xi);
1562       fh2Z2->Fill(jetBin, z);
1563     }
1564
1565
1566 }
1567
1568 //________________________________________________________________________
1569 void AliAnalysisTaskFragmentationFunction::AliFragFuncDiJetHistos::AddToOutput(TList* list)const
1570 {
1571   list->Add(fh1Jet1Pt);
1572   list->Add(fh1Jet2Pt);
1573   list->Add(fh1JetPt);
1574   list->Add(fh2TrackPtJet1);
1575   list->Add(fh2TrackPtJet2);
1576   list->Add(fh2TrackPtJet);
1577   list->Add(fh2Xi1);
1578   list->Add(fh2Xi2);
1579   list->Add(fh2Xi);
1580   list->Add(fh2Z1);
1581   list->Add(fh2Z2);
1582   list->Add(fh2Z);
1583 }
1584
1585 //______________________________________________________________________________________________________
1586 AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::AliFragFuncQADiJetHistos(const char* name, Int_t kindSlices,
1587                                             Int_t nInvMass, Float_t invMassMin, Float_t invMassMax, 
1588                                              Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
1589                                             Int_t nDeltaPhi, Float_t deltaPhiMin, Float_t deltaPhiMax, 
1590                                            Int_t nDeltaEta, Float_t deltaEtaMin, Float_t deltaEtaMax, 
1591                                           Int_t nDeltaPt, Float_t deltaPtMin, Float_t deltaPtMax)
1592   : TObject()
1593   ,fKindSlices(kindSlices)
1594   ,fNBinsJetInvMass(nInvMass)
1595   ,fJetInvMassMin(invMassMin)
1596   ,fJetInvMassMax(invMassMax)
1597   ,fNBinsJetPt(nJetPt)
1598   ,fJetPtMin(jetPtMin)
1599   ,fJetPtMax(jetPtMax)
1600   ,fNBinsDeltaPhi(nDeltaPhi)
1601   ,fDeltaPhiMin(deltaPhiMin)
1602   ,fDeltaPhiMax(deltaPhiMax)
1603   ,fNBinsDeltaEta(nDeltaEta)
1604   ,fDeltaEtaMin(deltaEtaMin)
1605   ,fDeltaEtaMax(deltaEtaMax)
1606   ,fNBinsDeltaPt(nDeltaPt)
1607   ,fDeltaPtMin(deltaPtMin)
1608   ,fDeltaPtMax(deltaPtMax)
1609   ,fh2InvMass(0)
1610   ,fh2DeltaPhi(0)
1611   ,fh2DeltaEta(0)
1612   ,fh2DeltaPt(0)
1613   ,fName(name)
1614 {
1615   // default constructor
1616
1617 }
1618
1619 //______________________________________________________________________________________________________
1620 AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::AliFragFuncQADiJetHistos(const AliFragFuncQADiJetHistos& copy)
1621   : TObject()
1622   ,fKindSlices(copy.fKindSlices)
1623   ,fNBinsJetInvMass(copy.fNBinsJetInvMass)
1624   ,fJetInvMassMin(copy.fJetInvMassMin)
1625   ,fJetInvMassMax(copy.fJetInvMassMax)
1626   ,fNBinsJetPt(copy.fNBinsJetPt)
1627   ,fJetPtMin(copy.fJetPtMin)
1628   ,fJetPtMax(copy.fJetPtMax)
1629   ,fNBinsDeltaPhi(copy.fNBinsDeltaPhi)
1630   ,fDeltaPhiMin(copy.fDeltaPhiMin)
1631   ,fDeltaPhiMax(copy.fDeltaPhiMax)
1632   ,fNBinsDeltaEta(copy.fNBinsDeltaEta)
1633   ,fDeltaEtaMin(copy.fDeltaEtaMin)
1634   ,fDeltaEtaMax(copy.fDeltaEtaMax)
1635   ,fNBinsDeltaPt(copy.fNBinsDeltaPt)
1636   ,fDeltaPtMin(copy.fDeltaPtMin)
1637   ,fDeltaPtMax(copy.fDeltaPtMax)
1638   ,fh2InvMass(copy.fh2InvMass)
1639   ,fh2DeltaPhi(copy.fh2DeltaPhi)
1640   ,fh2DeltaEta(copy.fh2DeltaEta)
1641   ,fh2DeltaPt(copy.fh2DeltaPt)
1642   ,fName(copy.fName)
1643 {
1644   // default constructor
1645
1646 }
1647
1648 //_______________________________________________________________________________________________________________________________________________________________
1649 AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos& AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::operator=(const AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos& o)
1650 {
1651   // assignment
1652   
1653   if(this!=&o){
1654     TObject::operator=(o);
1655     fKindSlices       = o.fKindSlices;
1656     fNBinsJetInvMass  = o.fNBinsJetInvMass;
1657     fJetInvMassMin    = o.fJetInvMassMin;
1658     fJetInvMassMax    = o.fJetInvMassMax;
1659     fNBinsJetPt       = o.fNBinsJetPt;
1660     fJetPtMin         = o.fJetPtMin;
1661     fJetPtMax         = o.fJetPtMax;
1662     fNBinsDeltaPhi    = o.fNBinsDeltaPhi;
1663     fDeltaPhiMin      = o.fDeltaPhiMin;
1664     fDeltaPhiMax      = o.fDeltaPhiMax;
1665     fNBinsDeltaEta    = o.fNBinsDeltaEta;
1666     fDeltaEtaMin      = o.fDeltaEtaMin;
1667     fDeltaEtaMax      = o.fDeltaEtaMax;
1668     fNBinsDeltaPt     = o.fNBinsDeltaPt;
1669     fDeltaPtMin       = o.fDeltaPtMin;
1670     fDeltaPtMax       = o.fDeltaPtMax;
1671     fh2InvMass        = o.fh2InvMass;
1672     fh2DeltaPhi       = o.fh2DeltaPhi;
1673     fh2DeltaEta       = o.fh2DeltaEta;
1674     fh2DeltaPt        = o.fh2DeltaPt;
1675     fName             = o.fName;
1676   }
1677     
1678   return *this;
1679 }
1680
1681 //_________________________________________________________
1682 AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::~AliFragFuncQADiJetHistos()
1683 {
1684   // destructor 
1685
1686   if(fh2InvMass)  delete fh2InvMass;
1687   if(fh2DeltaPhi) delete fh2DeltaPhi;
1688   if(fh2DeltaEta) delete fh2DeltaEta;
1689   if(fh2DeltaPt)  delete fh2DeltaPt;
1690 }
1691
1692 //________________________________________________________________________
1693 void AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::DefineQADiJetHistos()
1694 {
1695
1696   Int_t nBins = 0;
1697   Double_t min = 0.;
1698   Double_t max = 0.;
1699   const char *xaxis = "";
1700   if(fKindSlices == 1)
1701     {
1702       nBins = fNBinsJetInvMass;
1703       min   = fJetInvMassMin;
1704       max   = fJetInvMassMax;
1705       xaxis = "M_{JJ} [GeV]";
1706     }
1707   if(fKindSlices == 2 || fKindSlices == 3)
1708     {
1709       nBins = fNBinsJetPt;
1710       min   = fJetPtMin;
1711       max   = fJetPtMax;
1712       if(fKindSlices == 2) xaxis = "E_{Tmean} [GeV]";
1713       if(fKindSlices == 3) xaxis ="leading jet p_{T} [GeV/c]";
1714     }
1715   
1716   
1717   fh2InvMass  = new TH2F(Form("fh2DJInvMassPositionCut%s",  fName.Data()), "",nBins, min, max, fNBinsJetInvMass, fJetInvMassMin, fJetInvMassMax);
1718   fh2DeltaPhi = new TH2F(Form("fh2DJDeltaPhiPositionCut%s", fName.Data()), "",nBins, min, max, fNBinsDeltaPhi, fDeltaPhiMin, fDeltaPhiMax);
1719   fh2DeltaEta = new TH2F(Form("fh2DJDeltaEtaPositionCut%s", fName.Data()), "",nBins, min, max, fNBinsDeltaEta, fDeltaEtaMin, fDeltaEtaMax);
1720   fh2DeltaPt  = new TH2F(Form("fh2DJDeltaPtPositionCut%s",  fName.Data()), "",nBins, min, max, fNBinsDeltaPt, fDeltaPtMin, fDeltaPtMax);
1721   
1722   AliAnalysisTaskFragmentationFunction::SetProperties(fh2InvMass, xaxis, "Invariant Mass", "Entries");
1723   AliAnalysisTaskFragmentationFunction::SetProperties(fh2DeltaPhi, xaxis, "#Delta #phi", "Entries");
1724   AliAnalysisTaskFragmentationFunction::SetProperties(fh2DeltaEta, xaxis, "#Delta #eta", "Entries");
1725   AliAnalysisTaskFragmentationFunction::SetProperties(fh2DeltaPt, xaxis, "#Delta p_{T}", "Entries");
1726
1727 }
1728
1729 //________________________________________________________________________
1730 void AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::FillDiJetQA(Double_t invMass, Double_t deltaPhi, Double_t deltaEta,Double_t deltaPt, Double_t jetBin)
1731 {
1732   fh2InvMass->Fill(jetBin, invMass);
1733   fh2DeltaPhi->Fill(jetBin, deltaPhi);
1734   fh2DeltaEta->Fill(jetBin, deltaEta);
1735   fh2DeltaPt->Fill(jetBin, deltaPt);
1736 }
1737
1738 //________________________________________________________________________
1739 void AliAnalysisTaskFragmentationFunction::AliFragFuncQADiJetHistos::AddToOutput(TList* list)const
1740 {
1741   list->Add(fh2InvMass);
1742   list->Add(fh2DeltaPhi);
1743   list->Add(fh2DeltaEta);
1744   list->Add(fh2DeltaPt);
1745 }
1746
1747 //_________________________________________________________________________________
1748 void AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::AddToOutput(TList* list) const
1749 {
1750   // add histos to list
1751
1752   list->Add(fh2CosTheta);
1753   list->Add(fh2Theta);
1754   list->Add(fh2Jt);
1755
1756   list->Add(fhnIntraJet);
1757
1758 }
1759
1760
1761 Bool_t AliAnalysisTaskFragmentationFunction::Notify()
1762 {
1763   //
1764   // Implemented Notify() to read the cross sections
1765   // and number of trials from pyxsec.root
1766   // (taken from AliAnalysisTaskJetSpectrum)
1767   // 
1768   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
1769   Double_t xsection = 0;
1770   UInt_t   ntrials  = 0;
1771   Float_t   ftrials  = 0;
1772   if(tree){
1773     TFile *curfile = tree->GetCurrentFile();
1774     if (!curfile) {
1775       Error("Notify","No current file");
1776       return kFALSE;
1777     }
1778     if(!fh1Xsec||!fh1Trials){
1779       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
1780       return kFALSE;
1781     }
1782
1783     TString fileName(curfile->GetName());
1784     if(fileName.Contains("AliESDs.root")){
1785         fileName.ReplaceAll("AliESDs.root", "");
1786     }
1787     else if(fileName.Contains("AliAOD.root")){
1788         fileName.ReplaceAll("AliAOD.root", "");
1789     }
1790     else if(fileName.Contains("AliAODs.root")){
1791         fileName.ReplaceAll("AliAODs.root", "");
1792     }
1793     else if(fileName.Contains("galice.root")){
1794         // for running with galice and kinematics alone...                      
1795         fileName.ReplaceAll("galice.root", "");
1796     }
1797     TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"));
1798     if(!fxsec){
1799       if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec.root"));
1800       // next trial fetch the histgram file
1801       fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
1802       if(!fxsec){
1803         // not a severe condition
1804         if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec_hists.root"));        
1805         return kTRUE;
1806       }
1807       else{
1808         // find the tlist we want to be independtent of the name so use the Tkey
1809         TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
1810         if(!key){
1811           if(fDebug>0)Printf("%s:%d key not found in the file",(char*)__FILE__,__LINE__);       
1812           return kTRUE;
1813         }
1814         TList *list = dynamic_cast<TList*>(key->ReadObj());
1815         if(!list){
1816           if(fDebug>0)Printf("%s:%d key is not a tlist",(char*)__FILE__,__LINE__);      
1817           return kTRUE;
1818         }
1819         xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
1820         ftrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
1821       }
1822     }
1823     else{
1824       TTree *xtree = (TTree*)fxsec->Get("Xsection");
1825       if(!xtree){
1826         Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
1827       }
1828       xtree->SetBranchAddress("xsection",&xsection);
1829       xtree->SetBranchAddress("ntrials",&ntrials);
1830       ftrials = ntrials;
1831       xtree->GetEntry(0);
1832     }
1833     fh1Xsec->Fill("<#sigma>",xsection);
1834     fh1Trials->Fill("#sum{ntrials}",ftrials);
1835   }
1836   
1837   return kTRUE;
1838 }
1839
1840
1841
1842 //__________________________________________________________________
1843 void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
1844 {
1845   // create output objects
1846
1847   if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()");
1848  
1849   // create list of tracks and jets 
1850   
1851   fTracksRec = new TList();
1852   fTracksRec->SetOwner(kFALSE);  
1853
1854   fTracksRecCuts = new TList();
1855   fTracksRecCuts->SetOwner(kFALSE);  
1856
1857   fTracksGen = new TList();
1858   fTracksGen->SetOwner(kFALSE);
1859
1860   fTracksAODMCCharged = new TList();
1861   fTracksAODMCCharged->SetOwner(kFALSE);
1862     
1863   fTracksRecQualityCuts = new TList(); 
1864   fTracksRecQualityCuts->SetOwner(kFALSE);
1865
1866   fJetsRec = new TList();
1867   fJetsRec->SetOwner(kFALSE);
1868
1869   fJetsRecCuts = new TList();
1870   fJetsRecCuts->SetOwner(kFALSE);
1871
1872   fJetsGen = new TList();
1873   fJetsGen->SetOwner(kFALSE);
1874
1875   fJetsRecEff = new TList();
1876   fJetsRecEff->SetOwner(kFALSE);
1877
1878   // fJetsKine = new TList();
1879   // fJetsKine->SetOwner(kTRUE); // delete AOD jets using mom from Kine Tree via TList::Clear()
1880
1881
1882   //
1883   // Create histograms / output container
1884   //
1885
1886   OpenFile(1);
1887   fCommonHistList = new TList();
1888   
1889   Bool_t oldStatus = TH1::AddDirectoryStatus();
1890   TH1::AddDirectory(kFALSE);
1891   
1892   
1893   // Histograms 
1894   fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
1895   fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
1896   fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
1897   fh1EvtMult                 = new TH1F("fh1EvtMult","Event multiplicity, track pT cut > 150 MeV/c, |#eta| < 0.9",120,0.,120.);
1898
1899   fh1Xsec                    = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
1900   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1901   fh1Trials                  = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
1902   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1903   fh1PtHard                  = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
1904   fh1PtHardTrials            = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
1905
1906   fh1nRecJetsCuts            = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",10,-0.5,9.5);
1907   fh1nGenJets                = new TH1F("fh1nGenJets","generated jets per event",10,-0.5,9.5);
1908   fh1nRecEffJets             = new TH1F("fh1nRecEffJets","reconstruction effiency: jets per event",10,-0.5,9.5);
1909
1910   // 5D single track eff histo: phi:eta:gen pt:rec pt:isReconstructed -  use binning as for track QA
1911   Int_t   nBinsSingleTrackEffHisto[5]      = { fQATrackNBinsPhi, fQATrackNBinsEta, fQATrackNBinsPt, fQATrackNBinsPt, 2 };
1912   Double_t binMinSingleTrackEffHisto[5]    = { fQATrackPhiMin, fQATrackEtaMin, fQATrackPtMin, fQATrackPtMin, 0 };
1913   Double_t binMaxSingleTrackEffHisto[5]    = { fQATrackPhiMax, fQATrackEtaMax, fQATrackPtMax, fQATrackPtMax, 2 };
1914   const char* labelsSingleTrackEffHisto[5] = {"#phi","#eta","gen p_{T} [GeV/c]", "rec p_{T} [GeV/c]", "isRec"};
1915
1916   fhnSingleTrackRecEffHisto = new THnSparseF("fhnSingleTrackRecEffHisto","generated tracks phi:eta:pt:isReconstructed",5,
1917                                              nBinsSingleTrackEffHisto,binMinSingleTrackEffHisto,binMaxSingleTrackEffHisto);
1918   
1919   AliAnalysisTaskFragmentationFunction::SetProperties(fhnSingleTrackRecEffHisto,5,labelsSingleTrackEffHisto); 
1920   
1921   
1922   // 7D jets track eff histo: jet phi:eta:pt:track pt:z:xi:isReconstructed - use binning as for track/jet QA
1923   Int_t    nBinsJetTrackEffHisto[7]     = { fQAJetNBinsPhi, fQAJetNBinsEta, fFFNBinsJetPt, fFFNBinsPt, fFFNBinsZ, fFFNBinsXi, 2};
1924   Double_t binMinJetTrackEffHisto[7]    = { fQAJetPhiMin, fQAJetEtaMin, fFFJetPtMin , fFFPtMin, fFFZMin ,  fFFXiMin, 0 };
1925   Double_t binMaxJetTrackEffHisto[7]    = { fQAJetPhiMax, fQAJetEtaMax, fFFJetPtMax , fFFPtMax, fFFZMax ,  fFFXiMax, 2 };
1926   const char* labelsJetTrackEffHisto[7] = {"jet #phi","jet #eta","jet p_{T} [GeV/c]","track p_{T} [GeV/c]","z","#xi","isRec"};
1927
1928   fhnJetTrackRecEffHisto       = new THnSparseF("fhnJetTrackRecEffHisto","generated tracks - jet phi:jet eta:jet pt:track pt:z:xi:isReconstructed",7,
1929                                                 nBinsJetTrackEffHisto,binMinJetTrackEffHisto,binMaxJetTrackEffHisto);
1930
1931   AliAnalysisTaskFragmentationFunction::SetProperties(fhnJetTrackRecEffHisto,7,labelsJetTrackEffHisto);
1932
1933
1934   fQATrackHistosRec          = new AliFragFuncQATrackHistos("Rec", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax, 
1935                                                             fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1936                                                             fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax, 
1937                                                             fQATrackHighPtThreshold);
1938   fQATrackHistosRecCuts      = new AliFragFuncQATrackHistos("RecCuts", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax, 
1939                                                             fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1940                                                             fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax, 
1941                                                             fQATrackHighPtThreshold);
1942   fQATrackHistosGen          = new AliFragFuncQATrackHistos("Gen", fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax, 
1943                                                             fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1944                                                             fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax, 
1945                                                             fQATrackHighPtThreshold);
1946   
1947
1948   fQAJetHistosRec            = new AliFragFuncQAJetHistos("Rec", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
1949                                                           fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
1950                                                           fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1951   fQAJetHistosRecCuts        = new AliFragFuncQAJetHistos("RecCuts", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
1952                                                           fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
1953                                                           fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1954   fQAJetHistosRecCutsLeading = new AliFragFuncQAJetHistos("RecCutsLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
1955                                                           fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
1956                                                           fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1957   fQAJetHistosGen            = new AliFragFuncQAJetHistos("Gen", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
1958                                                           fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
1959                                                           fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1960   fQAJetHistosGenLeading     = new AliFragFuncQAJetHistos("GenLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
1961                                                           fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,
1962                                                           fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);  
1963   fQAJetHistosRecEffLeading  = new AliFragFuncQAJetHistos("RecEffLeading", fQAJetNBinsPt, fQAJetPtMin, fQAJetPtMax, 
1964                                                           fQAJetNBinsEta, fQAJetEtaMin, fQAJetEtaMax,fQAJetNBinsPhi, fQAJetPhiMin, fQAJetPhiMax);
1965
1966
1967   fFFHistosRecCuts           = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
1968                                                      fFFNBinsPt, fFFPtMin, fFFPtMax, 
1969                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,  
1970                                                      fFFNBinsZ , fFFZMin , fFFZMax);
1971   fFFHistosRecLeading        = new AliFragFuncHistos("RecLeading", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
1972                                                      fFFNBinsPt, fFFPtMin, fFFPtMax, 
1973                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,  
1974                                                      fFFNBinsZ , fFFZMin , fFFZMax);
1975   fFFHistosRecLeadingTrack   = new AliFragFuncHistos("RecLeadingTrack", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
1976                                                      fFFNBinsPt, fFFPtMin, fFFPtMax, 
1977                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,  
1978                                                      fFFNBinsZ , fFFZMin , fFFZMax);
1979   fFFHistosGen               = new AliFragFuncHistos("Gen", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
1980                                                      fFFNBinsPt, fFFPtMin, fFFPtMax, 
1981                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,  
1982                                                      fFFNBinsZ , fFFZMin , fFFZMax);
1983   fFFHistosGenLeading        = new AliFragFuncHistos("GenLeading", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
1984                                                      fFFNBinsPt, fFFPtMin, fFFPtMax, 
1985                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,  
1986                                                      fFFNBinsZ , fFFZMin , fFFZMax);
1987   fFFHistosGenLeadingTrack   = new AliFragFuncHistos("GenLeadingTrack", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
1988                                                      fFFNBinsPt, fFFPtMin, fFFPtMax, 
1989                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,  
1990                                                      fFFNBinsZ , fFFZMin , fFFZMax);
1991
1992
1993   fIJHistosRecCuts           = new AliFragFuncIntraJetHistos("RecCuts", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
1994                                                              fIJNBinsPt, fIJPtMin, fIJPtMax, 
1995                                                              fIJNBinsZ, fIJZMin, fIJZMax,  
1996                                                              fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
1997                                                              fIJNBinsTheta , fIJThetaMin , fIJThetaMax,
1998                                                              fIJNBinsJt , fIJJtMin , fIJJtMax);
1999   fIJHistosRecLeading        = new AliFragFuncIntraJetHistos("RecLeading", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
2000                                                              fIJNBinsPt, fIJPtMin, fIJPtMax, 
2001                                                              fIJNBinsZ, fIJZMin, fIJZMax,  
2002                                                              fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
2003                                                              fIJNBinsTheta , fIJThetaMin , fIJThetaMax, 
2004                                                              fIJNBinsJt , fIJJtMin , fIJJtMax);
2005   fIJHistosRecLeadingTrack   = new AliFragFuncIntraJetHistos("RecLeadingTrack", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
2006                                                              fIJNBinsPt, fIJPtMin, fIJPtMax, 
2007                                                              fIJNBinsZ, fIJZMin, fIJZMax,  
2008                                                              fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
2009                                                              fIJNBinsTheta , fIJThetaMin , fIJThetaMax, 
2010                                                              fIJNBinsJt , fIJJtMin , fIJJtMax);
2011   fIJHistosGen               = new AliFragFuncIntraJetHistos("Gen", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
2012                                                              fIJNBinsPt, fIJPtMin, fIJPtMax, 
2013                                                              fIJNBinsZ, fIJZMin, fIJZMax,  
2014                                                              fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
2015                                                              fIJNBinsTheta , fIJThetaMin , fIJThetaMax, 
2016                                                              fIJNBinsJt , fIJJtMin , fIJJtMax);
2017   fIJHistosGenLeading        = new AliFragFuncIntraJetHistos("GenLeading", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
2018                                                              fIJNBinsPt, fIJPtMin, fIJPtMax, 
2019                                                              fIJNBinsZ, fIJZMin, fIJZMax,  
2020                                                              fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
2021                                                              fIJNBinsTheta , fIJThetaMin , fIJThetaMax, 
2022                                                              fIJNBinsJt , fIJJtMin , fIJJtMax);
2023   fIJHistosGenLeadingTrack   = new AliFragFuncIntraJetHistos("GenLeadingTrack", fIJNBinsJetPt, fIJJetPtMin, fIJJetPtMax, 
2024                                                              fIJNBinsPt, fIJPtMin, fIJPtMax, 
2025                                                              fIJNBinsZ, fIJZMin, fIJZMax,  
2026                                                              fIJNBinsCosTheta , fIJCosThetaMin , fIJCosThetaMax, 
2027                                                              fIJNBinsTheta , fIJThetaMin , fIJThetaMax, 
2028                                                              fIJNBinsJt , fIJJtMin , fIJJtMax);
2029
2030
2031   fFFDiJetHistosRecCuts         = new AliFragFuncDiJetHistos("RecCuts", fDiJetKindBins, 
2032                                                              fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
2033                                                              fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax,
2034                                                              fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax, 
2035                                                              fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax, 
2036                                                              fDiJetNBinsZ, fDiJetZMin, fDiJetZMax);
2037   fFFDiJetHistosRecLeading      = new AliFragFuncDiJetHistos("RecLeading", fDiJetKindBins, 
2038                                                              fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
2039                                                              fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax, 
2040                                                              fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax, 
2041                                                              fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax, 
2042                                                              fDiJetNBinsZ, fDiJetZMin, fDiJetZMax); 
2043   fFFDiJetHistosRecLeadingTrack = new AliFragFuncDiJetHistos("RecLeadingTrack", fDiJetKindBins, 
2044                                                              fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
2045                                                              fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax, 
2046                                                              fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax, 
2047                                                              fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax, 
2048                                                              fDiJetNBinsZ, fDiJetZMin, fDiJetZMax);
2049
2050   fFFDiJetHistosGen             = new AliFragFuncDiJetHistos("Gen", fDiJetKindBins, 
2051                                                              fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
2052                                                              fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax,
2053                                                              fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax, 
2054                                                              fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax,
2055                                                              fDiJetNBinsZ, fDiJetZMin, fDiJetZMax);
2056   fFFDiJetHistosGenLeading      = new AliFragFuncDiJetHistos("GenLeading", fDiJetKindBins, 
2057                                                              fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
2058                                                              fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax,
2059                                                              fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax,
2060                                                              fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax,
2061                                                              fDiJetNBinsZ, fDiJetZMin, fDiJetZMax);
2062   fFFDiJetHistosGenLeadingTrack = new AliFragFuncDiJetHistos("GenLeadingTrack", fDiJetKindBins, 
2063                                                              fDiJetNBinsJetInvMass, fDiJetJetInvMassMin, fDiJetJetInvMassMax,
2064                                                              fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax, 
2065                                                              fDiJetNBinsPt, fDiJetPtMin, fDiJetPtMax, 
2066                                                              fDiJetNBinsXi, fDiJetXiMin, fDiJetXiMax, 
2067                                                              fDiJetNBinsZ, fDiJetZMin, fDiJetZMax);
2068
2069   fQADiJetHistosRecCuts = new AliFragFuncQADiJetHistos("RecCuts", fDiJetKindBins, 
2070                                                        fQADiJetNBinsInvMass, fQADiJetInvMassMin, fQADiJetInvMassMax,
2071                                                        fQADiJetNBinsJetPt, fQADiJetJetPtMin, fQADiJetJetPtMax,
2072                                                        fQADiJetNBinsDeltaPhi, fQADiJetDeltaPhiMin, fQADiJetDeltaPhiMax , 
2073                                                        fQADiJetNBinsDeltaEta, fQADiJetDeltaEtaMin, fQADiJetDeltaEtaMax , 
2074                                                        fQADiJetNBinsDeltaPt, fQADiJetDeltaPtMin, fQADiJetDeltaPtMax);
2075   fQADiJetHistosGen     = new AliFragFuncQADiJetHistos("Gen", fDiJetKindBins, 
2076                                                        fQADiJetNBinsInvMass, fQADiJetInvMassMin,  fQADiJetInvMassMax, 
2077                                                        fDiJetNBinsJetPt, fDiJetJetPtMin, fDiJetJetPtMax,
2078                                                        fQADiJetNBinsDeltaPhi, fQADiJetDeltaPhiMin, fQADiJetDeltaPhiMax,
2079                                                        fQADiJetNBinsDeltaEta, fQADiJetDeltaEtaMin, fQADiJetDeltaEtaMax,
2080                                                        fQADiJetNBinsDeltaPt, fQADiJetDeltaPtMin, fQADiJetDeltaPtMax);
2081
2082
2083   fQATrackHistosRec->DefineHistos();
2084   fQATrackHistosRecCuts->DefineHistos();
2085   fQATrackHistosGen->DefineHistos();
2086
2087   fQAJetHistosRec->DefineHistos();
2088   fQAJetHistosRecCuts->DefineHistos();
2089   fQAJetHistosRecCutsLeading->DefineHistos();
2090   fQAJetHistosGen->DefineHistos();
2091   fQAJetHistosGenLeading->DefineHistos();
2092   fQAJetHistosRecEffLeading->DefineHistos();
2093
2094   fFFHistosRecCuts->DefineHistos();
2095   fFFHistosRecLeading->DefineHistos();
2096   fFFHistosRecLeadingTrack->DefineHistos();
2097   fFFHistosGen->DefineHistos();
2098   fFFHistosGenLeading->DefineHistos();
2099   fFFHistosGenLeadingTrack->DefineHistos();
2100
2101   fIJHistosRecCuts->DefineHistos();
2102   fIJHistosRecLeading->DefineHistos();
2103   fIJHistosRecLeadingTrack->DefineHistos();
2104   fIJHistosGen->DefineHistos();
2105   fIJHistosGenLeading->DefineHistos();
2106   fIJHistosGenLeadingTrack->DefineHistos();
2107   
2108   fFFDiJetHistosRecCuts->DefineDiJetHistos();
2109   fFFDiJetHistosRecLeading->DefineDiJetHistos();
2110   fFFDiJetHistosRecLeadingTrack->DefineDiJetHistos();
2111   fFFDiJetHistosGen->DefineDiJetHistos();
2112   fFFDiJetHistosGenLeading->DefineDiJetHistos();
2113   fFFDiJetHistosGenLeadingTrack->DefineDiJetHistos();
2114   fQADiJetHistosRecCuts->DefineQADiJetHistos();
2115   fQADiJetHistosGen->DefineQADiJetHistos();
2116
2117   Bool_t genJets    = (fJetTypeGen != kJetsUndef) ? kTRUE : kFALSE;
2118   Bool_t genTracks  = (fTrackTypeGen != kTrackUndef) ? kTRUE : kFALSE;
2119   Bool_t recJetsEff = (fJetTypeRecEff != kJetsUndef) ? kTRUE : kFALSE;
2120
2121
2122   const Int_t saveLevel = 5;
2123   if(saveLevel>0){
2124     fCommonHistList->Add(fh1EvtSelection);
2125     fFFHistosRecCuts->AddToOutput(fCommonHistList);
2126     fFFHistosRecLeading->AddToOutput(fCommonHistList);
2127     fFFHistosRecLeadingTrack->AddToOutput(fCommonHistList);
2128
2129     if(genJets && genTracks){
2130        fFFHistosGen->AddToOutput(fCommonHistList);
2131        fFFHistosGenLeading->AddToOutput(fCommonHistList);
2132        fFFHistosGenLeadingTrack->AddToOutput(fCommonHistList);
2133
2134        fCommonHistList->Add(fh1Xsec);
2135        fCommonHistList->Add(fh1Trials);
2136        fCommonHistList->Add(fh1PtHard);
2137        fCommonHistList->Add(fh1PtHardTrials);
2138     }
2139   }
2140   if(saveLevel>1){
2141     fQATrackHistosRec->AddToOutput(fCommonHistList);
2142     fQATrackHistosRecCuts->AddToOutput(fCommonHistList);
2143     if(genTracks) fQATrackHistosGen->AddToOutput(fCommonHistList);
2144     
2145     fQAJetHistosRec->AddToOutput(fCommonHistList);
2146     fQAJetHistosRecCuts->AddToOutput(fCommonHistList);
2147     fQAJetHistosRecCutsLeading->AddToOutput(fCommonHistList);
2148     if(recJetsEff) fQAJetHistosRecEffLeading->AddToOutput(fCommonHistList); 
2149
2150     if(genJets){
2151        fQAJetHistosGen->AddToOutput(fCommonHistList);
2152        fQAJetHistosGenLeading->AddToOutput(fCommonHistList);
2153     }
2154
2155     fCommonHistList->Add(fh1EvtMult);
2156     fCommonHistList->Add(fh1nRecJetsCuts);
2157     if(genJets) fCommonHistList->Add(fh1nGenJets);
2158   }
2159   if(saveLevel>2){
2160     fCommonHistList->Add(fh1VertexNContributors);
2161     fCommonHistList->Add(fh1VertexZ);    
2162   }
2163   if(saveLevel>3){
2164     fIJHistosRecCuts->AddToOutput(fCommonHistList);
2165     fIJHistosRecLeading->AddToOutput(fCommonHistList);
2166     fIJHistosRecLeadingTrack->AddToOutput(fCommonHistList);
2167
2168     if(genJets && genTracks){
2169        fIJHistosGen->AddToOutput(fCommonHistList);
2170        fIJHistosGenLeading->AddToOutput(fCommonHistList);
2171        fIJHistosGenLeadingTrack->AddToOutput(fCommonHistList);
2172     }
2173   }
2174   if(saveLevel>4){
2175     fFFDiJetHistosRecCuts->AddToOutput(fCommonHistList);
2176     fFFDiJetHistosRecLeading->AddToOutput(fCommonHistList);
2177     fFFDiJetHistosRecLeadingTrack->AddToOutput(fCommonHistList);
2178     fQADiJetHistosRecCuts->AddToOutput(fCommonHistList);
2179
2180     if(genJets && genTracks){
2181         fFFDiJetHistosGen->AddToOutput(fCommonHistList);
2182         fFFDiJetHistosGenLeading->AddToOutput(fCommonHistList);
2183         fFFDiJetHistosGenLeadingTrack->AddToOutput(fCommonHistList);
2184         fQADiJetHistosGen->AddToOutput(fCommonHistList);
2185     }
2186
2187     if(recJetsEff && genTracks){
2188        fCommonHistList->Add(fhnSingleTrackRecEffHisto);
2189        fCommonHistList->Add(fhnJetTrackRecEffHisto);
2190        fCommonHistList->Add(fh1nRecEffJets);
2191     }
2192   }
2193
2194   // =========== Switch on Sumw2 for all histos ===========
2195   for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
2196     TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
2197     if (h1) h1->Sumw2();
2198     else{
2199       THnSparse *hnSparse = dynamic_cast<THnSparse*>(fCommonHistList->At(i));
2200       if(hnSparse) hnSparse->Sumw2();
2201     }
2202   }
2203   
2204   TH1::AddDirectory(oldStatus);
2205 }
2206
2207 //_______________________________________________
2208 void AliAnalysisTaskFragmentationFunction::Init()
2209 {
2210   // Initialization
2211   if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::Init()");
2212
2213 }
2214
2215 //_____________________________________________________________
2216 void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *) 
2217 {
2218   // Main loop
2219   // Called for each event
2220   if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::UserExec()");
2221         
2222
2223   if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
2224   // Trigger selection
2225   
2226   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
2227     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
2228   if(inputHandler->IsEventSelected() & AliVEvent::kMB){
2229     if(fDebug > 1)  Printf(" Trigger Selection: event ACCEPTED ... ");
2230     fh1EvtSelection->Fill(1.);
2231   } else {
2232     fh1EvtSelection->Fill(0.);
2233     if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
2234       if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
2235       PostData(1, fCommonHistList);
2236       return;
2237     }
2238   }
2239   
2240   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
2241   if(!fESD){
2242     if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
2243   }
2244   
2245   fMCEvent = MCEvent();
2246   if(!fMCEvent){
2247     if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
2248   }
2249   
2250   // get AOD event from input/ouput
2251   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
2252   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
2253     fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
2254     if (fDebug > 1)  Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
2255   }
2256   else {
2257     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
2258     if( handler && handler->InheritsFrom("AliAODHandler") ) {
2259       fAOD  = ((AliAODHandler*)handler)->GetAOD();
2260       if (fDebug > 1)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
2261     }
2262   }
2263   
2264   if(!fAOD){
2265     Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
2266     return;
2267   }
2268   
2269
2270   // event selection (vertex) *****************************************
2271   
2272   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
2273   Int_t nTracksPrim = primVtx->GetNContributors();
2274   fh1VertexNContributors->Fill(nTracksPrim);
2275   
2276   if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
2277   if(!nTracksPrim){
2278     if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__); 
2279     fh1EvtSelection->Fill(2.);
2280     PostData(1, fCommonHistList);
2281     return;
2282   }
2283
2284   fh1VertexZ->Fill(primVtx->GetZ());
2285   
2286   if(TMath::Abs(primVtx->GetZ())>10){
2287     if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ()); 
2288     fh1EvtSelection->Fill(3.);
2289     PostData(1, fCommonHistList);
2290     return; 
2291   }
2292
2293   TString primVtxName(primVtx->GetName());
2294
2295   if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
2296     if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
2297     fh1EvtSelection->Fill(4.);
2298     PostData(1, fCommonHistList);
2299     return;
2300   }
2301   if (fDebug > 1) Printf("%s:%d primary vertex selection: event ACCEPTED ...",(char*)__FILE__,__LINE__); 
2302   fh1EvtSelection->Fill(5.);
2303
2304
2305   //___ get MC information __________________________________________________________________
2306
2307   Double_t ptHard = 0.;
2308   Double_t nTrials = 1; // trials for MC trigger weight for real data
2309  
2310   AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(fMCEvent);
2311   if(!pythiaGenHeader){
2312      if(fJetTypeGen != kJetsUndef && fTrackTypeGen != kTrackUndef){
2313         Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
2314         return;
2315      }
2316   } else {
2317      nTrials = pythiaGenHeader->Trials();
2318      ptHard  = pythiaGenHeader->GetPtHard();
2319
2320      fh1PtHard->Fill(ptHard);
2321      fh1PtHardTrials->Fill(ptHard,nTrials);
2322   }
2323   
2324   
2325   //___ fetch jets __________________________________________________________________________
2326   
2327   Int_t nJ = GetListOfJets(fJetsRec, kJetsRec);
2328   Int_t nRecJets = 0;
2329   if(nJ>=0) nRecJets = fJetsRec->GetEntries();
2330   if(fDebug>2)Printf("%s:%d Selected Rec jets: %d %d",(char*)__FILE__,__LINE__,nJ,nRecJets);
2331   if(nJ != nRecJets) Printf("%s:%d Mismatch Selected Rec Jets: %d %d",(char*)__FILE__,__LINE__,nJ,nRecJets);
2332
2333   Int_t nJCuts = GetListOfJets(fJetsRecCuts, kJetsRecAcceptance);
2334   Int_t nRecJetsCuts = 0;
2335   if(nJCuts>=0) nRecJetsCuts = fJetsRecCuts->GetEntries();
2336   if(fDebug>2)Printf("%s:%d Selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
2337   if(nRecJetsCuts != nJCuts) Printf("%s:%d Mismatch selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
2338   fh1nRecJetsCuts->Fill(nRecJetsCuts);
2339
2340   
2341   if(fJetTypeGen==kJetsKine || fJetTypeGen == kJetsKineAcceptance) fJetsGen->SetOwner(kTRUE); // kine aod jets allocated on heap, delete them with TList::Clear() 
2342   Int_t nJGen  = GetListOfJets(fJetsGen, fJetTypeGen);
2343   Int_t nGenJets = 0;
2344   if(nJGen>=0) nGenJets = fJetsGen->GetEntries();
2345   if(fDebug>2)Printf("%s:%d Selected Gen jets: %d %d",(char*)__FILE__,__LINE__,nJGen,nGenJets);
2346   if(nJGen != nGenJets) Printf("%s:%d Mismatch selected Gen jets: %d %d",(char*)__FILE__,__LINE__,nJGen,nGenJets);
2347   fh1nGenJets->Fill(nGenJets);
2348
2349
2350   if(fJetTypeRecEff==kJetsKine || fJetTypeRecEff == kJetsKineAcceptance) fJetsRecEff->SetOwner(kTRUE); // kine aod jets allocated on heap, delete them with TList::Clear() 
2351   Int_t nJRecEff  = GetListOfJets(fJetsRecEff, fJetTypeRecEff);
2352   Int_t nRecEffJets = 0;
2353   if(nJRecEff>=0) nRecEffJets = fJetsRecEff->GetEntries();
2354   if(fDebug>2)Printf("%s:%d Selected RecEff jets: %d %d",(char*)__FILE__,__LINE__,nJRecEff,nRecEffJets);
2355   if(nJRecEff != nRecEffJets) Printf("%s:%d Mismatch selected RecEff jets: %d %d",(char*)__FILE__,__LINE__,nJRecEff,nRecEffJets);
2356   fh1nRecEffJets->Fill(nRecEffJets);
2357
2358
2359   //____ fetch particles __________________________________________________________
2360   
2361   Int_t nT = GetListOfTracks(fTracksRec, kTrackAOD);
2362   Int_t nRecPart = 0;
2363   if(nT>=0) nRecPart = fTracksRec->GetEntries();
2364   if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,nRecPart);
2365   if(nRecPart != nT) Printf("%s:%d Mismatch selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,nRecPart);
2366   
2367
2368   Int_t nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts);
2369   Int_t nRecPartCuts = 0;
2370   if(nTCuts>=0) nRecPartCuts = fTracksRecCuts->GetEntries();
2371   if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,nRecPartCuts);
2372   if(nRecPartCuts != nTCuts) Printf("%s:%d Mismatch selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,nRecPartCuts);
2373   fh1EvtMult->Fill(nRecPartCuts);
2374
2375
2376   Int_t nTGen = GetListOfTracks(fTracksGen,fTrackTypeGen);
2377   Int_t nGenPart = 0;
2378   if(nTGen>=0) nGenPart = fTracksGen->GetEntries();
2379   if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nTGen,nGenPart);
2380   if(nGenPart != nTGen) Printf("%s:%d Mismatch selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nTGen,nGenPart);
2381   
2382   
2383   //____ analysis, fill histos ___________________________________________________
2384   
2385   // loop over tracks
2386
2387   for(Int_t it=0; it<nRecPart; ++it){
2388     AliVParticle *part = dynamic_cast<AliVParticle*>(fTracksRec->At(it));
2389     fQATrackHistosRec->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt());
2390   }
2391   for(Int_t it=0; it<nRecPartCuts; ++it){
2392     AliVParticle *part = dynamic_cast<AliVParticle*>(fTracksRecCuts->At(it));
2393     fQATrackHistosRecCuts->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt());
2394   }
2395   for(Int_t it=0; it<nGenPart; ++it){
2396     AliVParticle *part = dynamic_cast<AliVParticle*>(fTracksGen->At(it));
2397     fQATrackHistosGen->FillTrackQA( part->Eta(), TVector2::Phi_0_2pi(part->Phi()), part->Pt());
2398   }
2399   
2400   // loop over jets
2401
2402   for(Int_t ij=0; ij<nRecJets; ++ij){
2403
2404     AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRec->At(ij));
2405     fQAJetHistosRec->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
2406   }
2407   
2408
2409   for(Int_t ij=0; ij<nRecJetsCuts; ++ij){
2410
2411     AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(ij));
2412     fQAJetHistosRecCuts->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
2413
2414     if(ij==0){ // leading jet
2415       
2416       fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
2417       
2418       TList* jettracklist = new TList();
2419       Double_t sumPt = 0.;
2420       Float_t leadTrackPx = 0.;
2421       Float_t leadTrackPy = 0.;
2422       Float_t leadTrackPz = 0.;
2423       Float_t leadTrackP  = 0.;
2424       Float_t leadTrackPt = 0.;
2425       TLorentzVector* leadTrackV = new TLorentzVector();
2426       
2427       if(GetFFRadius()<=0){
2428         GetJetTracksTrackrefs(jettracklist, jet);
2429        } else {
2430         GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt);
2431       }
2432       
2433       for(Int_t it=0; it<jettracklist->GetSize(); ++it){
2434         Float_t trackPx = (dynamic_cast<AliVParticle*> (jettracklist->At(it)))->Px();
2435         Float_t trackPy = (dynamic_cast<AliVParticle*> (jettracklist->At(it)))->Py();
2436         Float_t trackPz = (dynamic_cast<AliVParticle*> (jettracklist->At(it)))->Pz();
2437         Float_t trackP = (dynamic_cast<AliVParticle*> (jettracklist->At(it)))->P();
2438         Float_t trackPt = (dynamic_cast<AliVParticle*> (jettracklist->At(it)))->Pt();
2439         Float_t jetPx = jet->Px();
2440         Float_t jetPy = jet->Py();
2441         Float_t jetPz = jet->Pz();
2442         Float_t jetP  = jet->P();
2443         Float_t jetPt = jet->Pt();
2444         TLorentzVector* trackV = new TLorentzVector();
2445         TLorentzVector *jetV = new TLorentzVector();
2446         trackV->SetPxPyPzE(trackPx,trackPy,trackPz,trackP);
2447         jetV->SetPxPyPzE(jetPx,jetPy,jetPz,jetP);
2448
2449         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2450         
2451         fFFHistosRecCuts->FillFF( trackPt, jetPt, incrementJetPt);
2452         fIJHistosRecCuts->FillIntraJet( trackV, jetV );
2453         
2454         if(it==0){ 
2455           leadTrackPx = trackPx;
2456           leadTrackPy = trackPy;
2457           leadTrackPz = trackPz;
2458           leadTrackP  = trackP;
2459           leadTrackPt = trackPt;
2460           fFFHistosRecLeadingTrack->FillFF( leadTrackPt, jetPt, kTRUE);
2461
2462           leadTrackV->SetPxPyPzE(leadTrackPx,leadTrackPy,leadTrackPz,leadTrackP);
2463           fIJHistosRecLeadingTrack->FillIntraJet( leadTrackV, jetV );
2464         }
2465         fFFHistosRecLeading->FillFF( trackPt, leadTrackPt , incrementJetPt);
2466         fIJHistosRecLeading->FillIntraJet( trackV, leadTrackV );
2467
2468         delete trackV;
2469         delete jetV;
2470       }
2471       
2472       delete leadTrackV;
2473       delete jettracklist;
2474     }
2475   }
2476         
2477
2478   for(Int_t ij=0; ij<nGenJets; ++ij){
2479
2480     AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsGen->At(ij));
2481     fQAJetHistosGen->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
2482     
2483     if(ij==0){ // leading jet
2484     
2485       fQAJetHistosGenLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt());
2486       
2487       TList* jettracklist = new TList();
2488       Double_t sumPt = 0.;
2489       Float_t leadTrackPx = 0.;
2490       Float_t leadTrackPy = 0.;
2491       Float_t leadTrackPz = 0.;
2492       Float_t leadTrackP  = 0.;
2493       Float_t leadTrackPt = 0.;
2494       TLorentzVector* leadTrackV = new TLorentzVector();
2495
2496       if(GetFFRadius()<=0){
2497         GetJetTracksTrackrefs(jettracklist, jet);
2498       } else {
2499         GetJetTracksPointing(fTracksGen, jettracklist, jet, GetFFRadius(), sumPt);
2500       }
2501       
2502       for(Int_t it=0; it<jettracklist->GetSize(); ++it){
2503         Float_t trackPx = (dynamic_cast<AliVParticle*>(jettracklist->At(it)))->Px();
2504         Float_t trackPy = (dynamic_cast<AliVParticle*>(jettracklist->At(it)))->Py();
2505         Float_t trackPz = (dynamic_cast<AliVParticle*>(jettracklist->At(it)))->Pz();
2506         Float_t trackP  = (dynamic_cast<AliVParticle*>(jettracklist->At(it)))->P();
2507         Float_t trackPt = (dynamic_cast<AliVParticle*>(jettracklist->At(it)))->Pt();
2508         Float_t jetPx = jet->Px();
2509         Float_t jetPy = jet->Py();
2510         Float_t jetPz = jet->Pz();
2511         Float_t jetP  = jet->P();
2512         Float_t jetPt = jet->Pt();
2513         TLorentzVector* trackV = new TLorentzVector();
2514         TLorentzVector  *jetV = new TLorentzVector();
2515         trackV->SetPxPyPzE(trackPx,trackPy,trackPz,trackP);
2516         jetV->SetPxPyPzE(jetPx,jetPy,jetPz,jetP);
2517
2518         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2519
2520         fFFHistosGen->FillFF( trackPt, jetPt, incrementJetPt);
2521         fIJHistosGen->FillIntraJet( trackV, jetV );
2522         
2523         if(it==0){ 
2524           leadTrackPx = trackPx;
2525           leadTrackPy = trackPy;
2526           leadTrackPz = trackPz;
2527           leadTrackP  = trackP;
2528           leadTrackPt = trackPt;
2529           fFFHistosGenLeadingTrack->FillFF( leadTrackPt, jetPt, kTRUE);
2530
2531           leadTrackV->SetPxPyPzE(leadTrackPx,leadTrackPy,leadTrackPz,leadTrackP);
2532           fIJHistosGenLeadingTrack->FillIntraJet( leadTrackV, jetV );
2533         }
2534         fFFHistosGenLeading->FillFF( trackPt, leadTrackPt, incrementJetPt);
2535         fIJHistosGenLeading->FillIntraJet( trackV, leadTrackV );
2536
2537         delete trackV;
2538         delete jetV;
2539       }
2540       
2541       delete leadTrackV;
2542       delete jettracklist;
2543     }
2544   }
2545
2546   //_______ DiJet part _____________________________________________________
2547
2548   if (nRecJetsCuts > 1)  // OB: debugged this
2549   {
2550
2551     AliAODJet* jet1 = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(0));
2552     AliAODJet* jet2 = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(1));
2553     
2554     // DiJet deltaphi calculation
2555     Double_t phi1 = TVector2::Phi_0_2pi(jet1->Phi());
2556     Double_t phi2 = TVector2::Phi_0_2pi(jet2->Phi());
2557     Double_t deltaPhi = TMath::Abs(phi1-phi2); 
2558     if (deltaPhi > TMath::Pi() && deltaPhi < 2*TMath::Pi()) deltaPhi = 2*TMath::Pi() - deltaPhi;
2559     
2560     // DiJet CDF cut calculation
2561     Double_t Et1     = TMath::Abs(jet1->E()*TMath::Sin(jet1->Theta()));
2562     Double_t Et2     = TMath::Abs(jet2->E()*TMath::Sin(jet2->Theta()));
2563     Double_t sumEt   = Et1 + Et2;
2564     Double_t normEt1PlusEt2   = TMath::Sqrt(Et1*Et1+Et2*Et2+2*Et1*Et2*TMath::Cos(deltaPhi));
2565     Double_t ratio = (Double_t)(normEt1PlusEt2/sumEt);
2566     
2567     // DiJet events selection
2568     Bool_t positionCut       = 0;
2569     Bool_t positionEnergyCut = 0;
2570     Bool_t cdfCut            = 0; 
2571
2572     // Position cut :
2573     if (deltaPhi > fDiJetDeltaPhiCut) positionCut = 1;
2574     // Position-Energy cut :
2575     if ((deltaPhi > fDiJetDeltaPhiCut) && ((jet2->Pt()) >= fDiJetPtFractionCut*(jet1->Pt()))) positionEnergyCut = 1;
2576     // CDF cut :
2577     if (ratio < fDiJetCDFCut) cdfCut = 1;
2578     
2579     Int_t go = 0;
2580     
2581     if (fDiJetCut == 1 && positionCut == 1) go = 1;
2582     if (fDiJetCut == 2 && positionEnergyCut == 1) go = 1;
2583     if (fDiJetCut == 3 && cdfCut == 1) go = 1;
2584
2585     if (go)
2586       {
2587         Double_t deltaEta      = TMath::Abs(jet1->Eta()-jet2->Eta());
2588         Double_t deltaPt       = TMath::Abs(jet1->Pt()-jet2->Pt());
2589         Double_t meanEt        = (Double_t)((Et1+Et2)/2.);
2590         Double_t invariantMass = (Double_t)InvMass(jet1,jet2);
2591         
2592         Double_t  jetBin = GetDiJetBin(invariantMass, jet1->Pt(), meanEt, fDiJetKindBins);
2593
2594         if (jetBin > 0)
2595           {
2596             fQADiJetHistosRecCuts->FillDiJetQA(invariantMass, deltaPhi, deltaEta, deltaPt, jetBin);
2597             
2598             TList* jettracklist1 = new TList();
2599             Double_t sumPt1      = 0.;
2600             Float_t leadTrackPt1 = 0;
2601             
2602             TList* jettracklist2 = new TList();
2603             Double_t sumPt2      = 0.;
2604             Float_t leadTrackPt2 = 0;
2605             
2606             if(GetFFRadius()<=0)
2607               {
2608                 GetJetTracksTrackrefs(jettracklist1, jet1);
2609                 GetJetTracksTrackrefs(jettracklist2, jet2);
2610               }
2611             else
2612               {
2613                 GetJetTracksPointing(fTracksRecCuts, jettracklist1, jet1, GetFFRadius(), sumPt1);
2614                 GetJetTracksPointing(fTracksRecCuts, jettracklist2, jet2, GetFFRadius(), sumPt2);
2615               }
2616             
2617             Int_t nTracks = jettracklist1->GetSize(); 
2618             if (jettracklist1->GetSize() < jettracklist2->GetSize()) nTracks = jettracklist2->GetSize();
2619             
2620             for(Int_t it=0; it<nTracks; ++it)
2621               {
2622                 if (it < jettracklist1->GetSize())
2623                   { 
2624                     Float_t trackPt1 = (dynamic_cast<AliVParticle*> (jettracklist1->At(it)))->Pt();
2625                     Float_t jetPt1   = jet1->Pt();
2626                     
2627                     Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2628                     
2629                     fFFDiJetHistosRecCuts->FillDiJetFF(1, trackPt1, jetPt1, jetBin, incrementJetPt);
2630                     fFFDiJetHistosRecCuts->FillDiJetFF(0, trackPt1, jetPt1, jetBin, incrementJetPt);
2631                     
2632                     if (it == 0)
2633                       {
2634                         leadTrackPt1 = trackPt1;
2635                         
2636                         fFFDiJetHistosRecLeadingTrack->FillDiJetFF(1, leadTrackPt1, jetPt1, jetBin, kTRUE); 
2637                         fFFDiJetHistosRecLeadingTrack->FillDiJetFF(0, leadTrackPt1, jetPt1, jetBin, kTRUE); 
2638                       }
2639                     
2640                     fFFDiJetHistosRecLeading->FillDiJetFF(1, trackPt1, leadTrackPt1, jetBin, incrementJetPt); 
2641                     fFFDiJetHistosRecLeading->FillDiJetFF(0, trackPt1, leadTrackPt1, jetBin, incrementJetPt); 
2642                   }
2643                 
2644                 if (it < jettracklist2->GetSize())
2645                   { 
2646                     Float_t trackPt2   = (dynamic_cast<AliVParticle*>(jettracklist2->At(it)))->Pt();
2647                     Float_t jetPt2     = jet2->Pt();
2648                     
2649                     Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2650                     
2651                     fFFDiJetHistosRecCuts->FillDiJetFF(2, trackPt2, jetPt2, jetBin, incrementJetPt);
2652                     fFFDiJetHistosRecCuts->FillDiJetFF(0, trackPt2, jetPt2, jetBin, incrementJetPt);
2653                     
2654                     if (it == 0)
2655                       {
2656                         leadTrackPt2 = trackPt2;
2657                         
2658                         fFFDiJetHistosRecLeadingTrack->FillDiJetFF(2, leadTrackPt2, jetPt2, jetBin, kTRUE);
2659                         fFFDiJetHistosRecLeadingTrack->FillDiJetFF(0, leadTrackPt2, jetPt2, jetBin, kTRUE);
2660                       }
2661                     
2662                     fFFDiJetHistosRecLeading->FillDiJetFF(2, trackPt2, leadTrackPt2, jetBin, incrementJetPt);
2663                     fFFDiJetHistosRecLeading->FillDiJetFF(0, trackPt2, leadTrackPt2, jetBin, incrementJetPt);
2664                   }
2665               } // End loop on tracks
2666
2667             delete jettracklist1;
2668             delete jettracklist2;
2669
2670           } // End if(jetBin > 0)
2671         else { Printf("Jet bins for di-jet studies not set !");}
2672       } // End if(go)
2673   } // End if(nRecJets > 1)
2674
2675   if (nGenJets > 1)
2676   {
2677     AliAODJet* jet1 = dynamic_cast<AliAODJet*>(fJetsGen->At(0));
2678     AliAODJet* jet2 = dynamic_cast<AliAODJet*>(fJetsGen->At(1));
2679
2680     Double_t deltaPhi = 0;
2681     Double_t phi1 = TVector2::Phi_0_2pi(jet1->Phi());
2682     Double_t phi2 = TVector2::Phi_0_2pi(jet2->Phi());
2683     deltaPhi      = TMath::Abs(phi1-phi2); 
2684     if (deltaPhi > TMath::Pi() && deltaPhi < 2*TMath::Pi()) deltaPhi = 2*TMath::Pi() - deltaPhi;
2685
2686     Double_t Et1            = TMath::Abs(jet1->E()*TMath::Sin(jet1->Theta()));
2687     Double_t Et2            = TMath::Abs(jet2->E()*TMath::Sin(jet2->Theta()));
2688     Double_t sumEt          = Et1 + Et2;
2689     Double_t normEt1PlusEt2 = TMath::Sqrt(Et1*Et1+Et2*Et2+2*Et1*Et2*TMath::Cos(deltaPhi));
2690     Double_t ratio          = (Double_t)(normEt1PlusEt2/sumEt);
2691
2692     // DiJet events selection
2693     Bool_t positionCut       = 0;
2694     Bool_t positionEnergyCut = 0;
2695     Bool_t cdfCut            = 0; 
2696
2697     // Position cut :
2698     if (deltaPhi > fDiJetDeltaPhiCut) positionCut = 1;
2699     // Position-Energy cut :
2700     if ((deltaPhi > fDiJetDeltaPhiCut) && ((jet2->Pt()) >= fDiJetPtFractionCut*(jet1->Pt()))) positionEnergyCut = 1;
2701     // CDF cut :
2702     if (ratio < fDiJetCDFCut) cdfCut = 1;    
2703
2704     Int_t go = 0;
2705
2706     if (fDiJetCut == 1 && positionCut == 1) go = 1;
2707     if (fDiJetCut == 2 && positionEnergyCut == 1) go = 1;
2708     if (fDiJetCut == 3 && cdfCut == 1) go = 1;
2709
2710     if (go)
2711     {
2712       Double_t deltaEta      = TMath::Abs(jet1->Eta()-jet2->Eta());
2713       Double_t deltaPt       = TMath::Abs(jet1->Pt()-jet2->Pt());
2714       Double_t meanEt        = (Double_t)((Et1+Et2)/2.);
2715       Double_t invariantMass = (Double_t)InvMass(jet1,jet2);
2716
2717       Double_t jetBin = GetDiJetBin(invariantMass, jet1->Pt(), meanEt, fDiJetKindBins);
2718
2719       if(jetBin > 0)
2720       {
2721         fQADiJetHistosGen->FillDiJetQA(invariantMass, deltaPhi, deltaEta, deltaPt, jetBin);
2722
2723         TList* jettracklist1 = new TList();
2724         Double_t sumPt1 = 0.;
2725         Float_t leadTrackPt1 = 0.;
2726
2727         TList* jettracklist2 = new TList();
2728         Double_t sumPt2 = 0.;
2729         Float_t leadTrackPt2 = 0.;
2730       
2731         if(GetFFRadius()<=0)
2732         {
2733           GetJetTracksTrackrefs(jettracklist1, jet1);
2734           GetJetTracksTrackrefs(jettracklist2, jet2);
2735         }
2736         else
2737         {
2738           GetJetTracksPointing(fTracksGen, jettracklist1, jet1, GetFFRadius(), sumPt1);
2739           GetJetTracksPointing(fTracksGen, jettracklist2, jet2, GetFFRadius(), sumPt2);
2740         }
2741       
2742         Int_t nTracks = jettracklist1->GetSize(); 
2743         if (jettracklist1->GetSize() < jettracklist2->GetSize()) nTracks = jettracklist2->GetSize();
2744
2745         for(Int_t it=0; it<nTracks; ++it)
2746         {
2747           if (it < jettracklist1->GetSize())
2748           { 
2749             Float_t trackPt1 = (dynamic_cast<AliAODTrack*>(jettracklist1->At(it)))->Pt();
2750             Float_t jetPt1 = jet1->Pt();
2751
2752             Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2753
2754             fFFDiJetHistosGen->FillDiJetFF( 1, trackPt1, jetPt1, jetBin, incrementJetPt);
2755             fFFDiJetHistosGen->FillDiJetFF( 0, trackPt1, jetPt1, jetBin, incrementJetPt);
2756
2757             if(it==0)
2758             { 
2759               leadTrackPt1 = trackPt1;
2760
2761               fFFDiJetHistosGenLeadingTrack->FillDiJetFF( 1, leadTrackPt1, jetPt1, jetBin, kTRUE);
2762               fFFDiJetHistosGenLeadingTrack->FillDiJetFF( 0, leadTrackPt1, jetPt1, jetBin, kTRUE);
2763             }
2764
2765             fFFDiJetHistosGenLeading->FillDiJetFF( 1, trackPt1, leadTrackPt1, jetBin, incrementJetPt);
2766             fFFDiJetHistosGenLeading->FillDiJetFF( 0, trackPt1, leadTrackPt1, jetBin, incrementJetPt);
2767           }
2768       
2769           if (it < jettracklist2->GetSize())
2770           { 
2771             Float_t trackPt2 = (dynamic_cast<AliAODTrack*>(jettracklist2->At(it)))->Pt();
2772             Float_t jetPt2 = jet2->Pt();
2773
2774             Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2775
2776             fFFDiJetHistosGen->FillDiJetFF( 2, trackPt2, jetPt2, jetBin, incrementJetPt);
2777             fFFDiJetHistosGen->FillDiJetFF( 0, trackPt2, jetPt2, jetBin, incrementJetPt);
2778         
2779             if (it==0)
2780             { 
2781               leadTrackPt2 = trackPt2;
2782
2783               fFFDiJetHistosGenLeadingTrack->FillDiJetFF( 2, leadTrackPt2, jetPt2, jetBin, kTRUE);
2784               fFFDiJetHistosGenLeadingTrack->FillDiJetFF( 0, leadTrackPt2, jetPt2, jetBin, kTRUE);
2785             }
2786
2787             fFFDiJetHistosGenLeading->FillDiJetFF( 2, trackPt2, leadTrackPt2, jetBin, incrementJetPt);
2788             fFFDiJetHistosGenLeading->FillDiJetFF( 0, trackPt2, leadTrackPt2, jetBin, incrementJetPt);
2789           }
2790         } // End loop on tracks
2791
2792         delete jettracklist1;
2793         delete jettracklist2;
2794
2795       } // End if(jetBin > 0)
2796       else { Printf("Jet bins for di-jet studies not set !");}
2797     } // End if (go)
2798   } // End if(nGenJets > 1)
2799
2800   
2801   // ____ efficiency _______________________________
2802
2803   // arrays for generated particles: reconstructed AOD track index, isPrimary flag
2804   TArrayI indexAODTr; 
2805   TArrayS isGenPrim; 
2806   
2807   // array for reconcstructed AOD tracks: generated particle index
2808   TArrayI indexMCTr; 
2809   
2810   Int_t  nTracksAODMCCharged = GetListOfTracks(fTracksAODMCCharged, kTrackAODMCCharged);
2811   if(fDebug>2)Printf("%s:%d selected AODMC tracks: %d ",(char*)__FILE__,__LINE__,nTracksAODMCCharged);
2812   
2813   Int_t  nTracksRecQualityCuts = GetListOfTracks(fTracksRecQualityCuts, kTrackAODQualityCuts);
2814   if(fDebug>2)Printf("%s:%d selected rec tracks quality after cuts, full acceptance/pt : %d ",(char*)__FILE__,__LINE__,nTracksRecQualityCuts);
2815   
2816   // associate gen and rec tracks, store indices in TArrays 
2817   AssociateGenRec(fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,indexMCTr,isGenPrim);
2818   
2819   // single track eff
2820   FillSingleTrackRecEffHisto(fhnSingleTrackRecEffHisto,fTracksAODMCCharged,fTracksRecQualityCuts,indexAODTr,isGenPrim);
2821
2822   // jet track eff
2823   for(Int_t ij=0; ij<nRecEffJets; ++ij){ 
2824     
2825     AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRecEff->At(ij));
2826     
2827     if(ij==0){ // leading jet
2828       
2829       TList* jettracklist = new TList();
2830       Double_t sumPt = 0.;
2831       
2832       GetJetTracksPointing(fTracksGen, jettracklist, jet, GetFFRadius(), sumPt); // for efficiency: gen tracks from pointing with gen/rec jet
2833       
2834       Double_t jetEta = jet->Eta();
2835       Double_t jetPhi = TVector2::Phi_0_2pi(jet->Phi());
2836       Double_t jetPt  = sumPt;
2837       
2838       fQAJetHistosRecEffLeading->FillJetQA( jetEta, jetPhi, jetPt );
2839       FillJetTrackRecEffHisto(fhnJetTrackRecEffHisto,jetPhi,jetEta,jetPt,jettracklist,fTracksAODMCCharged,indexAODTr,isGenPrim);
2840       
2841       delete jettracklist;
2842     }
2843   }
2844    
2845   //___________________
2846   
2847   fTracksRec->Clear();
2848   fTracksRecCuts->Clear();
2849   fTracksGen->Clear();
2850   fTracksAODMCCharged->Clear();
2851   fTracksRecQualityCuts->Clear();
2852
2853   fJetsRec->Clear();
2854   fJetsRecCuts->Clear();
2855   fJetsGen->Clear();
2856   fJetsRecEff->Clear();
2857
2858   //Post output data.
2859   PostData(1, fCommonHistList);
2860   
2861 }
2862
2863 //________________________________________________________________________________________
2864 Double_t AliAnalysisTaskFragmentationFunction::InvMass(AliAODJet* jet1, AliAODJet* jet2)
2865 {
2866
2867   Double_t invMass = 0.;
2868   invMass = TMath::Sqrt(pow(jet1->E()+jet2->E(),2) - pow(jet1->Px()+jet2->Px(),2) - 
2869                         pow(jet1->Py()+jet2->Py(),2) - pow(jet1->Pz()+jet2->Pz(),2));
2870
2871   return invMass;
2872
2873 }
2874
2875 //________________________________________________________________________________________
2876 Double_t AliAnalysisTaskFragmentationFunction::GetDiJetBin(Double_t invMass, Double_t leadingJetPt, Double_t EtMean, Int_t kindBins)
2877 {
2878   Double_t jetBinOk = 0.;
2879   Double_t jetBin = 0.;
2880
2881   Float_t stepInvMass = (fDiJetJetInvMassMax - fDiJetJetInvMassMin)/fDiJetNBinsJetInvMass;
2882   Float_t stepPt = (fDiJetJetPtMax - fDiJetJetPtMin)/fDiJetNBinsJetPt;
2883
2884   if (kindBins == 1)
2885     {
2886       for(Int_t i=0; i<fDiJetNBinsJetInvMass; ++i)
2887         {
2888           jetBin = fDiJetJetInvMassMin + i*stepInvMass/2.;
2889           if(((fDiJetJetInvMassMin+i*stepInvMass) <= invMass) &&
2890              (fDiJetJetInvMassMin + (i+1)*stepInvMass) > invMass) jetBinOk = jetBin;
2891         }
2892       jetBinOk = -1.;
2893     }
2894   else if (kindBins == 3)
2895     {
2896       for(Int_t i=0; i<fDiJetNBinsJetPt; ++i)
2897         {
2898           jetBin = fDiJetJetPtMin + i*stepPt/2.;
2899           if(((fDiJetJetPtMin+i*stepPt) <= EtMean) &&
2900              (fDiJetJetPtMin + (i+1)*stepPt) > EtMean) jetBinOk = jetBin;
2901         }
2902       jetBinOk = -1.;
2903     }
2904   else if (kindBins == 2)
2905     {
2906       for(Int_t i=0; i<fDiJetNBinsJetPt; ++i)
2907         {
2908           jetBin = fDiJetJetPtMin + i*stepPt/2.;
2909           if(((fDiJetJetPtMin+i*stepPt) <= leadingJetPt) &&
2910              (fDiJetJetPtMin + (i+1)*stepPt) > leadingJetPt) jetBinOk = jetBin;
2911         }
2912       jetBinOk = -1.;
2913     }
2914   else {Printf("WARNING: kindBins wrongly set ! Please make sure to call SetKindSlices() and set the kind parameter to 1, 2 or 3.\n");}
2915
2916   return jetBinOk;
2917
2918 }
2919
2920
2921 //______________________________________________________________
2922 void AliAnalysisTaskFragmentationFunction::Terminate(Option_t *) 
2923 {
2924   // terminated
2925
2926   if(fDebug > 1) printf("AliAnalysisTaskFragmentationFunction::Terminate() \n");
2927 }  
2928
2929 //_________________________________________________________________________________
2930 Int_t AliAnalysisTaskFragmentationFunction::GetListOfTracks(TList *list, Int_t type)
2931 {
2932   // fill list of tracks selected according to type
2933
2934   if(fDebug > 2) Printf("%s:%d Selecting tracks with %d", (char*)__FILE__,__LINE__,type);
2935   
2936   if(!list){
2937     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
2938     return -1;
2939   }
2940
2941   if(type==kTrackUndef) return 0;
2942   
2943   Int_t iCount = 0;
2944   if(type==kTrackAODCuts || type==kTrackAOD){
2945
2946     // all rec. tracks, esd filter mask, eta range
2947     if(!fAOD) return -1;
2948     
2949     for(Int_t it=0; it<fAOD->GetNumberOfTracks(); ++it){
2950       AliAODTrack *tr = fAOD->GetTrack(it);
2951       
2952       if(type == kTrackAODCuts){
2953         if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))   continue;
2954         if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
2955         if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue;
2956         if(tr->Pt()  < fTrackPtCut) continue;
2957       }
2958       list->Add(tr);
2959       iCount++;
2960     }
2961   }
2962   else if (type==kTrackKineAll || type==kTrackKineCharged || type==kTrackKineChargedAcceptance){
2963     // kine particles, all or rather charged
2964     if(!fMCEvent) return iCount;
2965     
2966     for(Int_t it=0; it<fMCEvent->GetNumberOfTracks(); ++it){
2967       AliMCParticle* part = (AliMCParticle*) fMCEvent->GetTrack(it);
2968       
2969       if(type == kTrackKineCharged || type == kTrackKineChargedAcceptance){
2970         if(part->Charge()==0) continue;
2971         
2972         if(type == kTrackKineChargedAcceptance && 
2973            (       part->Eta() < fTrackEtaMin
2974                 || part->Eta() > fTrackEtaMax
2975                 || part->Phi() < fTrackPhiMin
2976                 || part->Phi() > fTrackPhiMax 
2977                 || part->Pt()  < fTrackPtCut)) continue;
2978       }
2979       
2980       list->Add(part);
2981       iCount++;
2982     }
2983   }
2984   else if (type==kTrackAODMCCharged || type==kTrackAODMCAll || type==kTrackAODMCChargedAcceptance) {
2985     // MC particles (from AOD), physical primaries, all or rather charged or rather charged within acceptance
2986     if(!fAOD) return -1;
2987     
2988     TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
2989     if(!tca)return iCount;
2990     
2991     for(int it=0; it<tca->GetEntriesFast(); ++it){
2992       AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
2993       if(!part->IsPhysicalPrimary())continue;
2994       
2995       if (type==kTrackAODMCCharged || type==kTrackAODMCChargedAcceptance){
2996         if(part->Charge()==0) continue;
2997         if(type==kTrackAODMCChargedAcceptance && 
2998            (     part->Eta() > fTrackEtaMax
2999               || part->Eta() < fTrackEtaMin
3000               || part->Phi() > fTrackPhiMax
3001               || part->Phi() < fTrackPhiMin
3002               || part->Pt()  < fTrackPtCut)) continue;
3003       }
3004       
3005       list->Add(part);
3006       iCount++;
3007     }
3008   }
3009   
3010   list->Sort();
3011   return iCount;
3012   
3013 }
3014 // _______________________________________________________________________________
3015 Int_t AliAnalysisTaskFragmentationFunction::GetListOfJets(TList *list, Int_t type)
3016 {
3017   // fill list of jets selected according to type
3018   
3019   if(!list){
3020     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
3021     return -1;
3022   }
3023
3024   if(type == kJetsRec || type == kJetsRecAcceptance){ // reconstructed jets
3025
3026     if(fBranchRecJets.Length()==0){
3027       Printf("%s:%d no rec jet branch specified", (char*)__FILE__,__LINE__);
3028       if(fDebug>1)fAOD->Print();
3029       return 0;
3030     }
3031
3032     TClonesArray *aodRecJets = new TClonesArray();
3033     if(fBranchRecJets.Length()) aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRecJets.Data()));\r
3034     if(!aodRecJets)             aodRecJets = dynamic_cast<TClonesArray*>(fAOD->GetList()->FindObject(fBranchRecJets.Data()));\r
3035
3036     if(!aodRecJets){
3037       if(fBranchRecJets.Length()) Printf("%s:%d no reconstructed jet array with name %s in AOD", (char*)__FILE__,__LINE__,fBranchRecJets.Data());
3038
3039       if(fDebug>1)fAOD->Print();
3040       return 0;
3041     }
3042
3043     Int_t nRecJets = 0;
3044     
3045     for(Int_t ij=0; ij<aodRecJets->GetEntries(); ++ij){
3046
3047       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ij));
3048       if(!tmp) continue;
3049         
3050       if( tmp->Pt() < fJetPtCut ) continue;
3051       if( type == kJetsRecAcceptance &&
3052           (    tmp->Eta() < fJetEtaMin
3053             || tmp->Eta() > fJetEtaMax
3054             || tmp->Phi() < fJetPhiMin
3055             || tmp->Phi() > fJetPhiMax )) continue;
3056       
3057       list->Add(tmp);
3058           
3059       nRecJets++;
3060     }
3061
3062     list->Sort();
3063     return nRecJets;
3064     delete aodRecJets;
3065   }
3066   else if(type == kJetsKine || type == kJetsKineAcceptance){
3067     
3068     // generated jets
3069     Int_t nGenJets = 0;
3070     
3071     if(!fMCEvent){
3072       if(fDebug>1) Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
3073       return 0;
3074     }
3075     
3076     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(fMCEvent);
3077     if(!pythiaGenHeader){
3078       Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
3079       return 0;
3080     }
3081     
3082     // fetch the pythia generated jets
3083     for(int ip=0; ip<pythiaGenHeader->NTriggerJets(); ++ip){
3084       
3085       Float_t p[4];
3086       AliAODJet *jet = new AliAODJet();
3087       pythiaGenHeader->TriggerJet(ip, p);
3088       jet->SetPxPyPzE(p[0], p[1], p[2], p[3]);
3089
3090       if( type == kJetsKineAcceptance &&
3091           (    jet->Eta() < fJetEtaMin
3092             || jet->Eta() > fJetEtaMax
3093             || jet->Phi() < fJetPhiMin
3094             || jet->Phi() > fJetPhiMax )) continue;
3095       
3096       list->Add(jet);
3097       nGenJets++;
3098     }
3099     list->Sort();
3100     return nGenJets;
3101   }
3102   else if(type == kJetsGen || type == kJetsGenAcceptance ){
3103
3104     if(fBranchGenJets.Length()==0){
3105       if(fDebug>1) Printf("%s:%d no gen jet branch specified", (char*)__FILE__,__LINE__);
3106       return 0;
3107     }
3108     
3109     TClonesArray *aodGenJets = new TClonesArray();
3110     if(fBranchGenJets.Length()) aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGenJets.Data()));\r
3111     if(!aodGenJets)             aodGenJets = dynamic_cast<TClonesArray*>(fAOD->GetList()->FindObject(fBranchGenJets.Data()));\r
3112
3113     if(!aodGenJets){
3114       if(fDebug>0){
3115         if(fBranchGenJets.Length())         Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGenJets.Data());
3116       }
3117       if(fDebug>1)fAOD->Print();
3118       return 0;
3119     }
3120
3121     Int_t nGenJets = 0;
3122     
3123     for(Int_t ig=0; ig<aodGenJets->GetEntries(); ++ig){
3124           
3125       AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
3126       if(!tmp) continue;
3127           
3128       if( tmp->Pt() < fJetPtCut ) continue;
3129       if( type == kJetsGenAcceptance &&
3130           (    tmp->Eta() < fJetEtaMin
3131             || tmp->Eta() > fJetEtaMax
3132             || tmp->Phi() < fJetPhiMin
3133             || tmp->Phi() > fJetPhiMax )) continue;
3134       
3135       list->Add(tmp);
3136       
3137       nGenJets++;
3138     }
3139     list->Sort();
3140     return nGenJets;
3141     delete aodGenJets;
3142   } 
3143   else{
3144     if(fDebug>0)Printf("%s:%d no such type %d",(char*)__FILE__,__LINE__,type);
3145     return 0;
3146   }
3147 }
3148
3149 // _________________________________________________________________________________________________________
3150 void AliAnalysisTaskFragmentationFunction::SetProperties(THnSparse* h,const Int_t dim, const char** labels)
3151 {
3152   //Set properties of THnSparse 
3153
3154   for(Int_t i=0; i<dim; i++){
3155
3156     h->GetAxis(i)->SetTitle(labels[i]);
3157     h->GetAxis(i)->SetTitleColor(1);
3158   }
3159 }
3160
3161 // __________________________________________________________________________________________
3162 void AliAnalysisTaskFragmentationFunction::SetProperties(TH1* h,const char* x, const char* y)
3163 {
3164   //Set properties of histos (x and y title)
3165
3166   h->SetXTitle(x);
3167   h->SetYTitle(y);
3168   h->GetXaxis()->SetTitleColor(1);
3169   h->GetYaxis()->SetTitleColor(1);
3170 }
3171
3172 // _________________________________________________________________________________________________________
3173 void AliAnalysisTaskFragmentationFunction::SetProperties(TH2* h,const char* x, const char* y, const char* z)
3174 {
3175   //Set properties of histos (x,y and z title)
3176
3177   h->SetXTitle(x);
3178   h->SetYTitle(y);
3179   h->SetZTitle(z);
3180   h->GetXaxis()->SetTitleColor(1);
3181   h->GetYaxis()->SetTitleColor(1);
3182   h->GetZaxis()->SetTitleColor(1);
3183 }
3184
3185 // ________________________________________________________________________________________________________________________________________________________
3186 void AliAnalysisTaskFragmentationFunction::GetJetTracksPointing(TList* inputlist, TList* outputlist, AliAODJet* jet, const Double_t radius,Double_t& sumPt)
3187 {
3188   // fill list of tracks in cone around jet axis  
3189
3190   sumPt = 0;
3191
3192   Double_t jetMom[3];
3193   jet->PxPyPz(jetMom);
3194   TVector3 jet3mom(jetMom);
3195
3196   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
3197
3198     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
3199
3200     Double_t trackMom[3];
3201     track->PxPyPz(trackMom);
3202     TVector3 track3mom(trackMom);
3203
3204     Double_t dR = jet3mom.DeltaR(track3mom);
3205
3206     if(dR<radius){
3207
3208       outputlist->Add(track);
3209       
3210       sumPt += track->Pt();
3211     }
3212   }
3213   
3214   outputlist->Sort();
3215 }
3216
3217 // ___________________________________________________________________________________________
3218 void AliAnalysisTaskFragmentationFunction::GetJetTracksTrackrefs(TList* list, AliAODJet* jet)
3219 {
3220   // list of jet tracks from trackrefs
3221   
3222   Int_t nTracks = jet->GetRefTracks()->GetEntriesFast();
3223
3224   for (Int_t itrack=0; itrack<nTracks; itrack++) {
3225     
3226     AliVParticle* track = dynamic_cast<AliVParticle*>(jet->GetRefTracks()->At(itrack));
3227     if(!track){
3228       AliError("expected ref track not found ");
3229       continue;
3230     }
3231         
3232     list->Add(track);
3233   }
3234   
3235   list->Sort();
3236 }
3237
3238 // _ ________________________________________________________________________________________________________________________________
3239 void  AliAnalysisTaskFragmentationFunction::AssociateGenRec(TList* tracksAODMCCharged,TList* tracksRec, TArrayI& indexAODTr,TArrayI& indexMCTr,TArrayS& isGenPrim)
3240 {
3241   // associate generated and reconstructed tracks, fill TArrays of list indices
3242
3243
3244   Int_t nTracksRec  = tracksRec->GetSize();
3245   Int_t nTracksGen  = tracksAODMCCharged->GetSize();
3246   TClonesArray *tca = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
3247
3248   if(!nTracksGen) return;
3249   if(!tca)        return;
3250   
3251   // set size
3252   indexAODTr.Set(nTracksGen);
3253   indexMCTr.Set(nTracksRec);
3254   isGenPrim.Set(nTracksGen);
3255
3256   indexAODTr.Reset(-1);
3257   indexMCTr.Reset(-1);
3258   isGenPrim.Reset(0);
3259
3260   // loop over reconstructed tracks, get generated track 
3261
3262   for(Int_t iRec=0; iRec<nTracksRec; iRec++){ 
3263       
3264     AliAODTrack* rectrack = dynamic_cast<AliAODTrack*>(tracksRec->At(iRec)); 
3265
3266     Int_t label = rectrack->GetLabel();
3267
3268     // find MC track in our list
3269     AliAODMCParticle* gentrack = 0x0;
3270     if(label>=0) gentrack = dynamic_cast<AliAODMCParticle*> (tca->At(label));
3271
3272     Int_t listIndex = -1;
3273     if(gentrack) listIndex = tracksAODMCCharged->IndexOf(gentrack);
3274
3275     if(listIndex>=0){
3276
3277       indexAODTr[listIndex] = iRec;
3278       indexMCTr[iRec]       = listIndex;
3279     }
3280   } 
3281
3282
3283   // define primary sample for reconstruction efficiency
3284
3285   for(Int_t iGen=0; iGen<nTracksGen; iGen++){
3286
3287     AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (tracksAODMCCharged->At(iGen));
3288
3289     Int_t pdg = gentrack->GetPdgCode();    
3290
3291     // 211 - pi, 2212 - proton, 321 - Kaon, 11 - electron, 13 - muon
3292     if(TMath::Abs(pdg) == 211 || TMath::Abs(pdg) == 2212 || TMath::Abs(pdg) == 321 || 
3293        TMath::Abs(pdg) == 11 || TMath::Abs(pdg) == 13){
3294       
3295       isGenPrim[iGen] = kTRUE;
3296     }
3297   }
3298 }
3299
3300 // _____________________________________________________________________________________________________________________________________________
3301 void AliAnalysisTaskFragmentationFunction::FillSingleTrackRecEffHisto(THnSparse* histo, TList* tracksGen, TList* tracksRec, TArrayI& indexAODTr, TArrayS& isGenPrim){
3302
3303   // fill THnSparse for single track reconstruction efficiency
3304
3305   Int_t nTracksGen  = tracksGen->GetSize();
3306
3307   if(!nTracksGen) return;
3308
3309   for(Int_t iGen=0; iGen<nTracksGen; iGen++){
3310
3311     if(isGenPrim[iGen] != 1) continue; // select primaries
3312
3313     AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (tracksGen->At(iGen));
3314     
3315     Double_t ptGen  = gentrack->Pt();
3316     Double_t etaGen = gentrack->Eta();
3317     Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi());
3318
3319     // apply same acc & pt cuts as for FF 
3320     // could in principle also be done setting THNsparse axis limits before projecting, 
3321     // but then the binning needs to be fine grained enough 
3322
3323     if(etaGen < fTrackEtaMin || etaGen > fTrackEtaMax) continue;
3324     if(phiGen < fTrackPhiMin || phiGen > fTrackPhiMax) continue;
3325     if(ptGen  < fTrackPtCut) continue;
3326
3327     Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track 
3328     Double_t isRec =  0;
3329     Double_t ptRec = -1;
3330
3331     if(iRec>=0){
3332
3333       AliAODTrack* rectrack = dynamic_cast<AliAODTrack*>(tracksRec->At(iRec)); 
3334       ptRec = rectrack->Pt();
3335       isRec = 1;
3336     }
3337
3338     Double_t entries[5] = {phiGen,etaGen,ptGen,ptRec,isRec};
3339     histo->Fill(entries);
3340   }
3341 }
3342
3343 // ______________________________________________________________________________________________________________________________________________________
3344 void AliAnalysisTaskFragmentationFunction::FillJetTrackRecEffHisto(THnSparse* histo,Double_t jetPhi, Double_t jetEta, Double_t jetPt, TList* jetTrackList, 
3345                                                                    TList* tracksGen, TArrayI& indexAODTr, TArrayS& isGenPrim)
3346 {
3347   // fill THnSparse for jet track reconstruction efficiency
3348
3349   Int_t nTracksJet = jetTrackList->GetSize(); // list with AODMC tracks
3350
3351   if(!nTracksJet) return; 
3352
3353   for(Int_t iTr=0; iTr<nTracksJet; iTr++){
3354
3355     AliAODMCParticle* gentrack =  dynamic_cast<AliAODMCParticle*> (jetTrackList->At(iTr));
3356
3357     // find jet track in gen tracks list
3358     Int_t iGen = tracksGen->IndexOf(gentrack); 
3359
3360     if(iGen<0){
3361       if(fDebug>0) Printf("%s:%d gen jet track not found ",(char*)__FILE__,__LINE__);
3362       continue;
3363     }
3364
3365     if(isGenPrim[iGen] != 1) continue; // select primaries
3366     
3367     Double_t ptGen  = gentrack->Pt();
3368     Double_t etaGen = gentrack->Eta();
3369     Double_t phiGen = TVector2::Phi_0_2pi(gentrack->Phi());
3370
3371     // apply same acc & pt cuts as for FF 
3372     // could in principle also be done setting THNsparse axis limits before projecting, 
3373     // but then the binning needs to be fine grained enough 
3374
3375     if(etaGen < fTrackEtaMin || etaGen > fTrackEtaMax) continue;
3376     if(phiGen < fTrackPhiMin || phiGen > fTrackPhiMax) continue;
3377     if(ptGen  < fTrackPtCut) continue;
3378
3379     Double_t z = ptGen / jetPt;
3380     Double_t xi = 0;
3381     if(z>0) xi = TMath::Log(1/z);
3382
3383     Double_t isRec =  0;
3384     Int_t iRec = indexAODTr[iGen]; // can be -1 if no good reconstructed track 
3385     if(iRec>=0) isRec = 1;
3386
3387     Double_t entries[7] = {jetPhi,jetEta,jetPt,ptGen,z,xi,isRec};
3388     histo->Fill(entries);
3389   }
3390 }