c418e3d4d5c73d13e1e3fb964d803fa90b04189e
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetCore.cxx
1 #include "TChain.h"
2 #include "TTree.h"
3 #include "TMath.h"
4 #include "TH1F.h"
5 #include "TH2F.h"
6 #include "TH3F.h"
7 #include "THnSparse.h"
8 #include "TCanvas.h"
9
10 #include "AliLog.h"
11
12 #include "AliAnalysisTask.h"
13 #include "AliAnalysisManager.h"
14
15 #include "AliVEvent.h"
16 #include "AliESDEvent.h"
17 #include "AliESDInputHandler.h"
18 #include "AliCentrality.h"
19 #include "AliAnalysisHelperJetTasks.h"
20 #include "AliInputEventHandler.h"
21 #include "AliAODJetEventBackground.h"
22 #include "AliAnalysisTaskFastEmbedding.h"
23
24 #include "AliAODEvent.h"
25 #include "AliAODHandler.h"
26 #include "AliAODJet.h"
27
28 #include "AliAnalysisTaskJetCore.h"
29
30 ClassImp(AliAnalysisTaskJetCore)
31
32 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
33 AliAnalysisTaskSE(),
34 fESD(0x0),
35 fAOD(0x0),
36 fAODExtension(0x0),
37 fBackgroundBranch(""),
38 fNonStdFile(""),
39 fIsPbPb(kTRUE),
40 fOfflineTrgMask(AliVEvent::kAny),
41 fMinContribVtx(1),
42 fVtxZMin(-8.),
43 fVtxZMax(8.),
44 fEvtClassMin(0),
45 fEvtClassMax(4),
46 fRadioFrac(0.2),
47 fMinDist(0.1),
48 fCentMin(0.),
49 fCentMax(100.),
50 fNInputTracksMin(0),
51 fNInputTracksMax(-1),
52 fAngStructCloseTracks(0),
53 fJetEtaMin(-.5),
54 fJetEtaMax(.5),
55 fJetPtMin(20.),
56 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
57 fJetPtFractionMin(0.5),
58 fNMatchJets(4),
59 fMatchMaxDist(0.8),
60 fKeepJets(kFALSE),
61 fkNbranches(2),
62 fkEvtClasses(12),
63 fOutputList(0x0),
64 fbEvent(kTRUE),
65 fHistEvtSelection(0x0),
66 fHistJetSelection(0x0),
67 fh2JetSelection(0x0),
68 fh2JetCoreMethod1C10(0x0),
69 fh2JetCoreMethod2C10(0x0),
70 fh2JetCoreMethod3C10(0x0),
71 fh2JetCoreMethod1C20(0x0),
72 fh2JetCoreMethod2C20(0x0),
73 fh2JetCoreMethod3C20(0x0),
74 fh2JetCoreMethod1C30(0x0),
75 fh2JetCoreMethod2C30(0x0),
76 fh2JetCoreMethod3C30(0x0),
77 fh2JetCoreMethod1C60(0x0),
78 fh2JetCoreMethod2C60(0x0),
79 fh2JetCoreMethod3C60(0x0),
80 fh2SumPtInC10(0x0),
81 fh2SumPtInC20(0x0),
82 fh2SumPtInC30(0x0),
83 fh2SumPtInC60(0x0),
84 fh2SumPtOutC10(0x0),
85 fh2SumPtOutC10b(0x0),
86 fh2SumPtOutC20(0x0),
87 fh2SumPtOutC30(0x0),
88 fh2SumPtOutC60(0x0),
89 fh2SumPtInC10bkg(0x0),
90 fh2SumPtInC20bkg(0x0),
91 fh2SumPtInC30bkg(0x0),
92 fh2SumPtInC60bkg(0x0),
93 fh2SumPtOutC10bkg(0x0),
94 fh2SumPtOutC20bkg(0x0),
95 fh2SumPtOutC30bkg(0x0),
96 fh2SumPtOutC60bkg(0x0),
97 fh2DeltaRC10pt1(0x0),
98 fh2DeltaRC20pt1(0x0),
99 fh2DeltaRC30pt1(0x0),
100 fh2DeltaRC60pt1(0x0),  
101 fh2DeltaRC10pt2(0x0),
102 fh2DeltaRC20pt2(0x0),
103 fh2DeltaRC30pt2(0x0),
104 fh2DeltaRC60pt2(0x0),  
105 fh2DeltaRC10pt3(0x0),
106 fh2DeltaRC20pt3(0x0),
107 fh2DeltaRC30pt3(0x0),
108 fh2DeltaRC60pt3(0x0),  
109 fh2DeltaRC10pt4(0x0),
110 fh2DeltaRC20pt4(0x0),
111 fh2DeltaRC30pt4(0x0),
112 fh2DeltaRC60pt4(0x0),
113 fh2DeltaEtaC10pt1(0x0),
114 fh2DeltaEtaC20pt1(0x0),
115 fh2DeltaEtaC30pt1(0x0),
116 fh2DeltaEtaC60pt1(0x0),  
117 fh2DeltaEtaC10pt2(0x0),
118 fh2DeltaEtaC20pt2(0x0),
119 fh2DeltaEtaC30pt2(0x0),
120 fh2DeltaEtaC60pt2(0x0),  
121 fh2DeltaEtaC10pt3(0x0),
122 fh2DeltaEtaC20pt3(0x0),
123 fh2DeltaEtaC30pt3(0x0),
124 fh2DeltaEtaC60pt3(0x0),  
125 fh2DeltaEtaC10pt4(0x0),
126 fh2DeltaEtaC20pt4(0x0),
127 fh2DeltaEtaC30pt4(0x0),
128 fh2DeltaEtaC60pt4(0x0),
129 fh2DeltaPhiC10pt1(0x0),
130 fh2DeltaPhiC20pt1(0x0),
131 fh2DeltaPhiC30pt1(0x0),
132 fh2DeltaPhiC60pt1(0x0),  
133 fh2DeltaPhiC10pt2(0x0),
134 fh2DeltaPhiC20pt2(0x0),
135 fh2DeltaPhiC30pt2(0x0),
136 fh2DeltaPhiC60pt2(0x0),  
137 fh2DeltaPhiC10pt3(0x0),
138 fh2DeltaPhiC20pt3(0x0),
139 fh2DeltaPhiC30pt3(0x0),
140 fh2DeltaPhiC60pt3(0x0),  
141 fh2DeltaPhiC10pt4(0x0),
142 fh2DeltaPhiC20pt4(0x0),
143 fh2DeltaPhiC30pt4(0x0),
144 fh2DeltaPhiC60pt4(0x0),
145 fh2AngStructpt1C10(0x0),
146 fh2AngStructpt2C10(0x0),
147 fh2AngStructpt3C10(0x0),
148 fh2AngStructpt4C10(0x0),
149 fh2AngStructpt1C20(0x0),
150 fh2AngStructpt2C20(0x0),
151 fh2AngStructpt3C20(0x0),
152 fh2AngStructpt4C20(0x0),    
153 fh2AngStructpt1C30(0x0),
154 fh2AngStructpt2C30(0x0),
155 fh2AngStructpt3C30(0x0),
156 fh2AngStructpt4C30(0x0),   
157 fh2AngStructpt1C60(0x0),
158 fh2AngStructpt2C60(0x0),
159 fh2AngStructpt3C60(0x0),
160 fh2AngStructpt4C60(0x0)
161 {
162    // default Constructor
163
164    fJetBranchName[0] = "";
165    fJetBranchName[1] = "";
166
167    fListJets[0] = new TList;
168    fListJets[1] = new TList;
169 }
170
171 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
172 AliAnalysisTaskSE(name),
173 fESD(0x0),
174 fAOD(0x0),
175 fAODExtension(0x0),
176 fBackgroundBranch(""),
177 fNonStdFile(""),
178 fIsPbPb(kTRUE),
179 fOfflineTrgMask(AliVEvent::kAny),
180 fMinContribVtx(1),
181 fVtxZMin(-8.),
182 fVtxZMax(8.),
183 fEvtClassMin(0),
184 fEvtClassMax(4),
185 fRadioFrac(0.2),
186 fMinDist(0.1),
187 fCentMin(0.),
188 fCentMax(100.),
189 fNInputTracksMin(0),
190 fNInputTracksMax(-1),
191 fAngStructCloseTracks(0),
192 fJetEtaMin(-.5),
193 fJetEtaMax(.5),
194 fJetPtMin(20.),
195 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
196 fJetPtFractionMin(0.5),
197 fNMatchJets(4),
198 fMatchMaxDist(0.8),
199 fKeepJets(kFALSE),
200 fkNbranches(2),
201 fkEvtClasses(12),
202 fOutputList(0x0),
203 fbEvent(kTRUE),
204 fHistEvtSelection(0x0),
205 fHistJetSelection(0x0),
206 fh2JetSelection(0x0),
207 fh2JetCoreMethod1C10(0x0),
208 fh2JetCoreMethod2C10(0x0),
209 fh2JetCoreMethod3C10(0x0),
210 fh2JetCoreMethod1C20(0x0),
211 fh2JetCoreMethod2C20(0x0),
212 fh2JetCoreMethod3C20(0x0),
213 fh2JetCoreMethod1C30(0x0),
214 fh2JetCoreMethod2C30(0x0),
215 fh2JetCoreMethod3C30(0x0),
216 fh2JetCoreMethod1C60(0x0),
217 fh2JetCoreMethod2C60(0x0),
218 fh2JetCoreMethod3C60(0x0),
219 fh2SumPtInC10(0x0),
220 fh2SumPtInC20(0x0),
221 fh2SumPtInC30(0x0),
222 fh2SumPtInC60(0x0),
223 fh2SumPtOutC10(0x0),
224 fh2SumPtOutC10b(0x0),
225 fh2SumPtOutC20(0x0),
226 fh2SumPtOutC30(0x0),
227 fh2SumPtOutC60(0x0),
228 fh2SumPtInC10bkg(0x0),
229 fh2SumPtInC20bkg(0x0),
230 fh2SumPtInC30bkg(0x0),
231 fh2SumPtInC60bkg(0x0),
232 fh2SumPtOutC10bkg(0x0),
233 fh2SumPtOutC20bkg(0x0),
234 fh2SumPtOutC30bkg(0x0),
235 fh2SumPtOutC60bkg(0x0),
236 fh2DeltaRC10pt1(0x0),
237 fh2DeltaRC20pt1(0x0),
238 fh2DeltaRC30pt1(0x0),
239 fh2DeltaRC60pt1(0x0),  
240 fh2DeltaRC10pt2(0x0),
241 fh2DeltaRC20pt2(0x0),
242 fh2DeltaRC30pt2(0x0),
243 fh2DeltaRC60pt2(0x0),  
244 fh2DeltaRC10pt3(0x0),
245 fh2DeltaRC20pt3(0x0),
246 fh2DeltaRC30pt3(0x0),
247 fh2DeltaRC60pt3(0x0),  
248 fh2DeltaRC10pt4(0x0),
249 fh2DeltaRC20pt4(0x0),
250 fh2DeltaRC30pt4(0x0),
251 fh2DeltaRC60pt4(0x0),
252 fh2DeltaEtaC10pt1(0x0),
253 fh2DeltaEtaC20pt1(0x0),
254 fh2DeltaEtaC30pt1(0x0),
255 fh2DeltaEtaC60pt1(0x0),  
256 fh2DeltaEtaC10pt2(0x0),
257 fh2DeltaEtaC20pt2(0x0),
258 fh2DeltaEtaC30pt2(0x0),
259 fh2DeltaEtaC60pt2(0x0),  
260 fh2DeltaEtaC10pt3(0x0),
261 fh2DeltaEtaC20pt3(0x0),
262 fh2DeltaEtaC30pt3(0x0),
263 fh2DeltaEtaC60pt3(0x0),  
264 fh2DeltaEtaC10pt4(0x0),
265 fh2DeltaEtaC20pt4(0x0),
266 fh2DeltaEtaC30pt4(0x0),
267 fh2DeltaEtaC60pt4(0x0),
268 fh2DeltaPhiC10pt1(0x0),
269 fh2DeltaPhiC20pt1(0x0),
270 fh2DeltaPhiC30pt1(0x0),
271 fh2DeltaPhiC60pt1(0x0),  
272 fh2DeltaPhiC10pt2(0x0),
273 fh2DeltaPhiC20pt2(0x0),
274 fh2DeltaPhiC30pt2(0x0),
275 fh2DeltaPhiC60pt2(0x0),  
276 fh2DeltaPhiC10pt3(0x0),
277 fh2DeltaPhiC20pt3(0x0),
278 fh2DeltaPhiC30pt3(0x0),
279 fh2DeltaPhiC60pt3(0x0),  
280 fh2DeltaPhiC10pt4(0x0),
281 fh2DeltaPhiC20pt4(0x0),
282 fh2DeltaPhiC30pt4(0x0),
283 fh2DeltaPhiC60pt4(0x0),
284 fh2AngStructpt1C10(0x0),
285 fh2AngStructpt2C10(0x0),
286 fh2AngStructpt3C10(0x0),
287 fh2AngStructpt4C10(0x0),
288 fh2AngStructpt1C20(0x0),
289 fh2AngStructpt2C20(0x0),
290 fh2AngStructpt3C20(0x0),
291 fh2AngStructpt4C20(0x0),    
292 fh2AngStructpt1C30(0x0),
293 fh2AngStructpt2C30(0x0),
294 fh2AngStructpt3C30(0x0),
295 fh2AngStructpt4C30(0x0),   
296 fh2AngStructpt1C60(0x0),
297 fh2AngStructpt2C60(0x0),
298 fh2AngStructpt3C60(0x0),
299 fh2AngStructpt4C60(0x0)    
300
301  {
302    // Constructor
303
304    fJetBranchName[0] = "";
305    fJetBranchName[1] = "";
306
307    fListJets[0] = new TList;
308    fListJets[1] = new TList;
309
310    DefineOutput(1, TList::Class());
311 }
312
313 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
314 {
315    delete fListJets[0];
316    delete fListJets[1];
317 }
318
319 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
320 {
321    fJetBranchName[0] = branch1;
322    fJetBranchName[1] = branch2;
323 }
324
325 void AliAnalysisTaskJetCore::Init()
326 {
327
328    // check for jet branches
329    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
330       AliError("Jet branch name not set.");
331    }
332
333 }
334
335 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
336 {
337    // Create histograms
338    // Called once
339    OpenFile(1);
340    if(!fOutputList) fOutputList = new TList;
341    fOutputList->SetOwner(kTRUE);
342
343    Bool_t oldStatus = TH1::AddDirectoryStatus();
344    TH1::AddDirectory(kFALSE);
345
346
347    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
348    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
349    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
350    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
351    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
352    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
353    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
354
355    fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
356    fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
357    fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
358    fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
359    fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
360    fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
361    fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
362    fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
363    fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
364
365    fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
366    fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
367    fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
368    fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
369    fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
370    fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
371    fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
372    fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
373    fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
374
375
376    //UInt_t entries = 0; // bit coded, see GetDimParams() below
377    //UInt_t opt = 0;  // bit coded, default (0) or high resolution (1)
378
379    //  Int_t bins[5]={10,10,150,150,15};
380    //Double_t xmin[5]={0.,0.,0.,0.,0.};
381    //Double_t xmax[5]={100.,6.5,150.,1.5,1.5};     
382    //fhnDeltaRjets = new THnSparseF("fhnDeltaRjets", "fhnDeltaRjets",5,bins,xmin,xmax);
383
384
385
386     fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
387     fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
388     fh2JetCoreMethod3C10 = new TH2F("JetCoreMethod3C10","",150, 0., 150.,100, 0., 1.5); 
389     fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
390     fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
391     fh2JetCoreMethod3C20 = new TH2F("JetCoreMethod3C20","",150, 0., 150.,100, 0., 1.5); 
392     fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
393     fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
394     fh2JetCoreMethod3C30 = new TH2F("JetCoreMethod3C30","",150, 0., 150.,100, 0., 1.5); 
395     fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
396     fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);
397     fh2JetCoreMethod3C60 = new TH2F("JetCoreMethod3C60","",150, 0., 150.,100, 0., 1.5); 
398    
399     fh2SumPtInC10 = new TH2F("PtSumInC10","",150, 0., 150.,5000, 0., 50.); 
400     fh2SumPtInC20 = new TH2F("PtSumInC20","",150, 0., 150.,5000, 0., 50.); 
401     fh2SumPtInC30 = new TH2F("PtSumInC30","",150, 0., 150.,5000, 0., 50.); 
402     fh2SumPtInC60 = new TH2F("PtSumInC60","",150, 0., 150.,5000, 0., 50.); 
403     fh2SumPtOutC10 = new TH2F("PtSumOutC10","",150, 0., 150.,5000, 0., 50.);
404     fh2SumPtOutC10b= new TH2F("PtSumOutC10b","",150, 0., 150.,500, 0., 50.);
405     fh2SumPtOutC20 = new TH2F("PtSumOutC20","",150, 0., 150.,5000, 0., 50.); 
406     fh2SumPtOutC30 = new TH2F("PtSumOutC30","",150, 0., 150.,5000, 0., 50.); 
407     fh2SumPtOutC60 = new TH2F("PtSumOutC60","",150, 0., 150.,5000, 0., 50.); 
408      
409     fh2SumPtInC10bkg = new TH2F("PtSumInC10bkg","",150, 0., 150.,5000, 0., 50.); 
410     fh2SumPtInC20bkg = new TH2F("PtSumInC20bkg","",150, 0., 150.,5000, 0., 50.); 
411     fh2SumPtInC30bkg = new TH2F("PtSumInC30bkg","",150, 0., 150.,5000, 0., 50.); 
412     fh2SumPtInC60bkg = new TH2F("PtSumInC60bkg","",150, 0., 150.,5000, 0., 50.); 
413     fh2SumPtOutC10bkg = new TH2F("PtSumOutC10bkg","",150, 0., 150.,5000, 0., 50.); 
414     fh2SumPtOutC20bkg = new TH2F("PtSumOutC20bkg","",150, 0., 150.,5000, 0., 50.); 
415     fh2SumPtOutC30bkg = new TH2F("PtSumOutC30bkg","",150, 0., 150.,5000, 0., 50.); 
416     fh2SumPtOutC60bkg = new TH2F("PtSumOutC60bkg","",150, 0., 150.,5000, 0., 50.); 
417
418
419
420     fh2DeltaRC10pt1 = new TH2F("DeltaRC10pt1","",150, 0., 150.,100,0.,1.5); 
421     fh2DeltaRC20pt1 = new TH2F("DeltaRC20pt1","",150, 0., 150.,100,0.,1.5); 
422     fh2DeltaRC30pt1 = new TH2F("DeltaRC30pt1","",150, 0., 150.,100,0.,1.5); 
423     fh2DeltaRC60pt1 = new TH2F("DeltaRC60pt1","",150, 0., 150.,100,0.,1.5); 
424     fh2DeltaRC10pt2 = new TH2F("DeltaRC10pt2","",150, 0., 150.,100,0.,1.5); 
425     fh2DeltaRC20pt2 = new TH2F("DeltaRC20pt2","",150, 0., 150.,100,0.,1.5); 
426     fh2DeltaRC30pt2 = new TH2F("DeltaRC30pt2","",150, 0., 150.,100,0.,1.5); 
427     fh2DeltaRC60pt2 = new TH2F("DeltaRC60pt2","",150, 0., 150.,100,0.,1.5); 
428     fh2DeltaRC10pt3 = new TH2F("DeltaRC10pt3","",150, 0., 150.,100,0.,1.5); 
429     fh2DeltaRC20pt3 = new TH2F("DeltaRC20pt3","",150, 0., 150.,100,0.,1.5); 
430     fh2DeltaRC30pt3 = new TH2F("DeltaRC30pt3","",150, 0., 150.,100,0.,1.5); 
431     fh2DeltaRC60pt3 = new TH2F("DeltaRC60pt3","",150, 0., 150.,100,0.,1.5); 
432     fh2DeltaRC10pt4 = new TH2F("DeltaRC10pt4","",150, 0., 150.,100,0.,1.5); 
433     fh2DeltaRC20pt4 = new TH2F("DeltaRC20pt4","",150, 0., 150.,100,0.,1.5); 
434     fh2DeltaRC30pt4 = new TH2F("DeltaRC30pt4","",150, 0., 150.,100,0.,1.5); 
435     fh2DeltaRC60pt4 = new TH2F("DeltaRC60pt4","",150, 0., 150.,100,0.,1.5); 
436
437     fh2DeltaEtaC10pt1 = new TH2F("DeltaEtaC10pt1","",150, 0., 150.,100,-1.5,1.5); 
438     fh2DeltaEtaC20pt1 = new TH2F("DeltaEtaC20pt1","",150, 0., 150.,100,-1.5,1.5); 
439     fh2DeltaEtaC30pt1 = new TH2F("DeltaEtaC30pt1","",150, 0., 150.,100,-1.5,1.5); 
440     fh2DeltaEtaC60pt1 = new TH2F("DeltaEtaC60pt1","",150, 0., 150.,100,-1.5,1.5); 
441     fh2DeltaEtaC10pt2 = new TH2F("DeltaEtaC10pt2","",150, 0., 150.,100,-1.5,1.5); 
442     fh2DeltaEtaC20pt2 = new TH2F("DeltaEtaC20pt2","",150, 0., 150.,100,-1.5,1.5); 
443     fh2DeltaEtaC30pt2 = new TH2F("DeltaEtaC30pt2","",150, 0., 150.,100,-1.5,1.5); 
444     fh2DeltaEtaC60pt2 = new TH2F("DeltaEtaC60pt2","",150, 0., 150.,100,-1.5,1.5); 
445     fh2DeltaEtaC10pt3 = new TH2F("DeltaEtaC10pt3","",150, 0., 150.,100,-1.5,1.5); 
446     fh2DeltaEtaC20pt3 = new TH2F("DeltaEtaC20pt3","",150, 0., 150.,100,-1.5,1.5); 
447     fh2DeltaEtaC30pt3 = new TH2F("DeltaEtaC30pt3","",150, 0., 150.,100,-1.5,1.5); 
448     fh2DeltaEtaC60pt3 = new TH2F("DeltaEtaC60pt3","",150, 0., 150.,100,-1.5,1.5); 
449     fh2DeltaEtaC10pt4 = new TH2F("DeltaEtaC10pt4","",150, 0., 150.,100,-1.5,1.5); 
450     fh2DeltaEtaC20pt4 = new TH2F("DeltaEtaC20pt4","",150, 0., 150.,100,-1.5,1.5); 
451     fh2DeltaEtaC30pt4 = new TH2F("DeltaEtaC30pt4","",150, 0., 150.,100,-1.5,1.5); 
452     fh2DeltaEtaC60pt4 = new TH2F("DeltaEtaC60pt4","",150, 0., 150.,100,-1.5,1.5); 
453     fh2DeltaPhiC10pt1 = new TH2F("DeltaPhiC10pt1","",150, 0., 150.,100,-6.5,6.5); 
454     fh2DeltaPhiC20pt1 = new TH2F("DeltaPhiC20pt1","",150, 0., 150.,100,-6.5,6.5); 
455     fh2DeltaPhiC30pt1 = new TH2F("DeltaPhiC30pt1","",150, 0., 150.,100,-6.5,6.5); 
456     fh2DeltaPhiC60pt1 = new TH2F("DeltaPhiC60pt1","",150, 0., 150.,100,-6.5,6.5); 
457     fh2DeltaPhiC10pt2 = new TH2F("DeltaPhiC10pt2","",150, 0., 150.,100,-6.5,6.5); 
458     fh2DeltaPhiC20pt2 = new TH2F("DeltaPhiC20pt2","",150, 0., 150.,100,-6.5,6.5); 
459     fh2DeltaPhiC30pt2 = new TH2F("DeltaPhiC30pt2","",150, 0., 150.,100,-6.5,6.5); 
460     fh2DeltaPhiC60pt2 = new TH2F("DeltaPhiC60pt2","",150, 0., 150.,100,-6.5,6.5); 
461     fh2DeltaPhiC10pt3 = new TH2F("DeltaPhiC10pt3","",150, 0., 150.,100,-6.5,6.5); 
462     fh2DeltaPhiC20pt3 = new TH2F("DeltaPhiC20pt3","",150, 0., 150.,100,-6.5,6.5); 
463     fh2DeltaPhiC30pt3 = new TH2F("DeltaPhiC30pt3","",150, 0., 150.,100,-6.5,6.5); 
464     fh2DeltaPhiC60pt3 = new TH2F("DeltaPhiC60pt3","",150, 0., 150.,100,-6.5,6.5); 
465     fh2DeltaPhiC10pt4 = new TH2F("DeltaPhiC10pt4","",150, 0., 150.,100,-6.5,6.5); 
466     fh2DeltaPhiC20pt4 = new TH2F("DeltaPhiC20pt4","",150, 0., 150.,100,-6.5,6.5); 
467     fh2DeltaPhiC30pt4 = new TH2F("DeltaPhiC30pt4","",150, 0., 150.,100,-6.5,6.5); 
468     fh2DeltaPhiC60pt4 = new TH2F("DeltaPhiC60pt4","",150, 0., 150.,100,-6.5,6.5); 
469
470     fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,5.);
471     fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,5.);
472     fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,5.);
473     fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,5.);
474     fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,5.);
475     fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,5.);
476     fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,5.);
477     fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,5.);
478     fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,5.);
479     fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,5.);
480     fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,5.);
481     fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,5.);
482     fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,5.);
483     fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,5.);
484     fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,5.);
485     fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,5.); 
486
487    fOutputList->Add(fHistEvtSelection);
488    fOutputList->Add(fHistJetSelection);
489    fOutputList->Add(fh2JetSelection);
490       fOutputList->Add(fh2JetCoreMethod1C10);
491       fOutputList->Add(fh2JetCoreMethod2C10);
492       fOutputList->Add(fh2JetCoreMethod3C10);
493       fOutputList->Add(fh2JetCoreMethod1C20);
494       fOutputList->Add(fh2JetCoreMethod2C20);
495       fOutputList->Add(fh2JetCoreMethod3C20);
496       fOutputList->Add(fh2JetCoreMethod1C30);
497       fOutputList->Add(fh2JetCoreMethod2C30);
498       fOutputList->Add(fh2JetCoreMethod3C30);
499       fOutputList->Add(fh2JetCoreMethod1C60);
500       fOutputList->Add(fh2JetCoreMethod2C60);
501       fOutputList->Add(fh2JetCoreMethod3C60);
502       fOutputList->Add(fh2SumPtInC10);
503       fOutputList->Add(fh2SumPtInC20);
504       fOutputList->Add(fh2SumPtInC30);
505       fOutputList->Add(fh2SumPtInC60);
506       fOutputList->Add(fh2SumPtOutC10);
507       fOutputList->Add(fh2SumPtOutC10b);
508       fOutputList->Add(fh2SumPtOutC20);
509       fOutputList->Add(fh2SumPtOutC30);
510       fOutputList->Add(fh2SumPtOutC60);
511       
512       fOutputList->Add(fh2SumPtInC10bkg);
513       fOutputList->Add(fh2SumPtInC20bkg);
514       fOutputList->Add(fh2SumPtInC30bkg);
515       fOutputList->Add(fh2SumPtInC60bkg);
516       fOutputList->Add(fh2SumPtOutC10bkg);
517       fOutputList->Add(fh2SumPtOutC20bkg);
518       fOutputList->Add(fh2SumPtOutC30bkg);
519       fOutputList->Add(fh2SumPtOutC60bkg);     
520       fOutputList->Add(fh2DeltaRC10pt1);
521       fOutputList->Add(fh2DeltaRC20pt1);
522       fOutputList->Add(fh2DeltaRC30pt1);
523       fOutputList->Add(fh2DeltaRC60pt1); 
524       fOutputList->Add(fh2DeltaRC10pt2);
525       fOutputList->Add(fh2DeltaRC20pt2);
526       fOutputList->Add(fh2DeltaRC30pt2);
527       fOutputList->Add(fh2DeltaRC60pt2); 
528       fOutputList->Add(fh2DeltaRC10pt3);
529       fOutputList->Add(fh2DeltaRC20pt3);
530       fOutputList->Add(fh2DeltaRC30pt3);
531       fOutputList->Add(fh2DeltaRC60pt3); 
532       fOutputList->Add(fh2DeltaRC10pt4);
533       fOutputList->Add(fh2DeltaRC20pt4);
534       fOutputList->Add(fh2DeltaRC30pt4);
535       fOutputList->Add(fh2DeltaRC60pt4);
536       
537       fOutputList->Add(fh2DeltaEtaC10pt1);
538       fOutputList->Add(fh2DeltaEtaC20pt1);
539       fOutputList->Add(fh2DeltaEtaC30pt1);
540       fOutputList->Add(fh2DeltaEtaC60pt1); 
541       fOutputList->Add(fh2DeltaEtaC10pt2);
542       fOutputList->Add(fh2DeltaEtaC20pt2);
543       fOutputList->Add(fh2DeltaEtaC30pt2);
544       fOutputList->Add(fh2DeltaEtaC60pt2); 
545       fOutputList->Add(fh2DeltaEtaC10pt3);
546       fOutputList->Add(fh2DeltaEtaC20pt3);
547       fOutputList->Add(fh2DeltaEtaC30pt3);
548       fOutputList->Add(fh2DeltaEtaC60pt3); 
549       fOutputList->Add(fh2DeltaEtaC10pt4);
550       fOutputList->Add(fh2DeltaEtaC20pt4);
551       fOutputList->Add(fh2DeltaEtaC30pt4);
552       fOutputList->Add(fh2DeltaEtaC60pt4); 
553       fOutputList->Add(fh2DeltaPhiC10pt1);
554       fOutputList->Add(fh2DeltaPhiC20pt1);
555       fOutputList->Add(fh2DeltaPhiC30pt1);
556       fOutputList->Add(fh2DeltaPhiC60pt1); 
557       fOutputList->Add(fh2DeltaPhiC10pt2);
558       fOutputList->Add(fh2DeltaPhiC20pt2);
559       fOutputList->Add(fh2DeltaPhiC30pt2);
560       fOutputList->Add(fh2DeltaPhiC60pt2); 
561       fOutputList->Add(fh2DeltaPhiC10pt3);
562       fOutputList->Add(fh2DeltaPhiC20pt3);
563       fOutputList->Add(fh2DeltaPhiC30pt3);
564       fOutputList->Add(fh2DeltaPhiC60pt3); 
565       fOutputList->Add(fh2DeltaPhiC10pt4);
566       fOutputList->Add(fh2DeltaPhiC20pt4);
567       fOutputList->Add(fh2DeltaPhiC30pt4);
568       fOutputList->Add(fh2DeltaPhiC60pt4);      
569      
570        fOutputList->Add(fh2AngStructpt1C10);
571        fOutputList->Add(fh2AngStructpt2C10);
572        fOutputList->Add(fh2AngStructpt3C10);
573        fOutputList->Add(fh2AngStructpt4C10); 
574        fOutputList->Add(fh2AngStructpt1C20);
575        fOutputList->Add(fh2AngStructpt2C20);
576        fOutputList->Add(fh2AngStructpt3C20);
577        fOutputList->Add(fh2AngStructpt4C20); 
578        fOutputList->Add(fh2AngStructpt1C30);
579        fOutputList->Add(fh2AngStructpt2C30);
580        fOutputList->Add(fh2AngStructpt3C30);
581        fOutputList->Add(fh2AngStructpt4C30);
582        fOutputList->Add(fh2AngStructpt1C60);
583        fOutputList->Add(fh2AngStructpt2C60);
584        fOutputList->Add(fh2AngStructpt3C60);
585        fOutputList->Add(fh2AngStructpt4C60);  
586
587
588
589
590    // =========== Switch on Sumw2 for all histos ===========
591    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
592       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
593       if (h1){
594          h1->Sumw2();
595          continue;
596       }
597     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
598       if (hn){
599          hn->Sumw2();
600       }   
601    }
602    TH1::AddDirectory(oldStatus);
603
604    PostData(1, fOutputList);
605 }
606
607 void AliAnalysisTaskJetCore::UserExec(Option_t *)
608 {
609    
610
611    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
612       AliError("Jet branch name not set.");
613       return;
614    }
615
616    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
617    if (!fESD) {
618       AliError("ESD not available");
619       fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
620    } else {
621       fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
622    }
623  
624     if(fNonStdFile.Length()!=0){
625     // case that we have an AOD extension we need can fetch the jets from the extended output
626     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
627     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
628     if(!fAODExtension){
629       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
630     }}
631     
632
633
634
635
636    // -- event selection --
637    fHistEvtSelection->Fill(1); // number of events before event selection
638
639    // physics selection
640    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
641    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
642    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
643       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
644       fHistEvtSelection->Fill(2);
645       PostData(1, fOutputList);
646       return;
647    }
648
649    // vertex selection
650    AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
651    Int_t nTracksPrim = primVtx->GetNContributors();
652    if ((nTracksPrim < fMinContribVtx) ||
653          (primVtx->GetZ() < fVtxZMin) ||
654          (primVtx->GetZ() > fVtxZMax) ){
655       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
656       fHistEvtSelection->Fill(3);
657       PostData(1, fOutputList);
658       return;
659    }
660
661    // event class selection (from jet helper task)
662    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
663    if(fDebug) Printf("Event class %d", eventClass);
664    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
665       fHistEvtSelection->Fill(4);
666       PostData(1, fOutputList);
667       return;
668    }
669
670    // centrality selection
671    AliCentrality *cent = 0x0;
672    Float_t centValue = 0.; 
673    if(fESD) cent = fESD->GetCentrality();
674    if(cent) centValue = cent->GetCentralityPercentile("V0M");
675    if(fDebug) printf("centrality: %f\n", centValue);
676    if (centValue < fCentMin || centValue > fCentMax){
677       fHistEvtSelection->Fill(4);
678       PostData(1, fOutputList);
679       return;
680    }
681
682
683    // multiplicity due to input tracks
684    //Int_t nInputTracks = GetNInputTracks();
685    //if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
686    //   fHistEvtSelection->Fill(5);
687    //   PostData(1, fOutputList);
688    //   return;
689    // }
690
691    
692    fHistEvtSelection->Fill(0); // accepted events  
693    // -- end event selection --
694   
695    // get background
696    AliAODJetEventBackground* externalBackground = 0;
697    if(fAOD&&!externalBackground&&fBackgroundBranch.Length()){
698       externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
699       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
700    }
701    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
702      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
703       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
704    }
705    
706    Float_t rho = 0;
707    if(externalBackground)rho = externalBackground->GetBackground(0);
708
709
710    // fetch jets
711    TClonesArray *aodJets[2];
712    aodJets[0]=0;
713    if(fAOD&&!aodJets[0]){
714    aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // 
715    aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data()));  }
716    if(fAODExtension && !aodJets[0]){ 
717      aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); // 
718      aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
719
720    TList ParticleList;
721    Int_t nT = GetListOfTracks(&ParticleList);
722      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
723       fListJets[iJetType]->Clear();
724       if (!aodJets[iJetType]) continue;
725
726       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
727       
728       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
729          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
730          if (jet) fListJets[iJetType]->Add(jet);
731       }
732       fListJets[iJetType]->Sort();
733    }
734    
735    Double_t etabig=0;
736    Double_t ptbig=0;
737    Double_t areabig=0;
738    Double_t phibig=0.;
739    Double_t etasmall=0;
740    Double_t ptsmall=0;
741    Double_t areasmall=0;
742    Double_t distr=0.;
743    Double_t phismall=0.;
744
745    Double_t up1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
746    Double_t up2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
747    Double_t up3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
748    Double_t up4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
749    Double_t down1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
750    Double_t down2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
751    Double_t down3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
752    Double_t down4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
753
754
755
756    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
757            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
758            etabig  = jetbig->Eta();
759            phibig  = jetbig->Phi();
760            ptbig   = jetbig->Pt();
761            if(ptbig==0) continue; 
762            areabig = jetbig->EffectiveAreaCharged();
763            Double_t ptcorr=ptbig-rho*areabig;
764            if(ptcorr<=0) continue;
765            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
766                    Double_t dismin=100.;
767                    Double_t ptmax=-10.; 
768                    Int_t index1=-1;
769                    Int_t index2=-1;
770                    Double_t fracin=0.;
771                    Double_t sumPtIn=0.;
772                    Double_t sumPtOut=0.; 
773            //compute sum of the pt of the tracks in a concentric cone
774            TRefArray *genTrackList = jetbig->GetRefTracks();
775            Int_t nTracksGenJet = genTrackList->GetEntriesFast();
776            AliAODTrack* genTrack;
777              for(Int_t ir=0; ir<nTracksGenJet; ++ir){
778              genTrack = (AliAODTrack*)(genTrackList->At(ir));
779              Float_t etr=genTrack->Eta();
780              Float_t phir=genTrack->Phi();
781              distr=(etr-etabig)*(etr-etabig)+(phir-phibig)*(phir-phibig);
782              distr=TMath::Sqrt(distr);
783              if(distr<=fRadioFrac){ fracin=fracin+genTrack->Pt();}}
784     
785              if(centValue<10) fh2JetCoreMethod3C10->Fill(ptcorr,fracin/ptbig);
786               
787              if((centValue>20)&&(centValue<40)) fh2JetCoreMethod3C20->Fill(ptcorr,fracin/ptbig);
788                                                 
789              if((centValue>30)&&(centValue<60)) fh2JetCoreMethod3C30->Fill(ptcorr,fracin/ptbig);
790                                                 
791              if(centValue>60) fh2JetCoreMethod3C60->Fill(ptcorr,fracin/ptbig);
792                              
793          
794              
795
796                 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
797                   AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
798                   etasmall  = jetsmall->Eta();
799                   phismall = jetsmall->Phi();
800                   ptsmall   = jetsmall->Pt();
801                   areasmall = jetsmall->EffectiveAreaCharged();
802                   Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
803                   tmpDeltaR=TMath::Sqrt(tmpDeltaR);
804                      //Fraction in the jet core  
805                     if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;  
806                     index2=j;}  
807                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
808                       index1=j;}} //en of loop over R=0.2 jets
809                 //method1:most concentric jet=core 
810                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));       
811                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
812                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
813                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
814                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
815                 //method2:hardest contained jet=core   
816                 if(index2!=-1){ 
817                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
818                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
819                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig); 
820                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
821                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }  
822
823      Double_t R=fRadioFrac*2.;      
824   
825      for(int it = 0;it<nT;++it){
826
827           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
828           Float_t deltaR = jetbig->DeltaR(part);
829           Float_t deltaEta = part->Eta()-etabig;
830           Float_t deltaPhi = part->Phi()-phibig;  
831           if((deltaR>=0.4)&&(deltaR<=0.6))sumPtIn=sumPtIn+part->Pt();                     
832           if((deltaR>=0.8)&&(deltaR<=1.))sumPtOut=sumPtOut+part->Pt();     
833           if(centValue<10.){
834                if((ptcorr>=70.)&&(ptcorr<85.)) {fh2DeltaRC10pt1->Fill(part->Pt(),deltaR);
835                if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R))fh2DeltaEtaC10pt1->Fill(part->Pt(),deltaEta);
836                if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC10pt1->Fill(part->Pt(),deltaPhi);}
837                if((ptcorr>=85.)&&(ptcorr<100.)) {fh2DeltaRC10pt2->Fill(part->Pt(),deltaR);
838                if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC10pt2->Fill(part->Pt(),deltaEta);
839                if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC10pt2->Fill(part->Pt(),deltaPhi); }
840                if((ptcorr>=100.)&&(ptcorr<120.)) {fh2DeltaRC10pt3->Fill(part->Pt(),deltaR);
841                if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC10pt3->Fill(part->Pt(),deltaEta);
842                if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC10pt3->Fill(part->Pt(),deltaPhi);}
843                if((ptcorr>=120.)&&(ptcorr<140.)) {fh2DeltaRC10pt4->Fill(part->Pt(),deltaR); 
844                if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC10pt4->Fill(part->Pt(),deltaEta);
845                if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC10pt4->Fill(part->Pt(),deltaPhi); }}
846
847           if((centValue>20.)&&(centValue<40.)){
848               if((ptcorr>=70.)&&(ptcorr<85.)) {fh2DeltaRC20pt1->Fill(part->Pt(),deltaR);
849               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R))fh2DeltaEtaC20pt1->Fill(part->Pt(),deltaEta);
850               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC20pt1->Fill(part->Pt(),deltaPhi);}
851               if((ptcorr>=85.)&&(ptcorr<100.)) {fh2DeltaRC20pt2->Fill(part->Pt(),deltaR);
852               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC20pt2->Fill(part->Pt(),deltaEta);
853               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC20pt2->Fill(part->Pt(),deltaPhi); }
854               if((ptcorr>=100.)&&(ptcorr<120.)) {fh2DeltaRC20pt3->Fill(part->Pt(),deltaR);
855               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC20pt3->Fill(part->Pt(),deltaEta);
856               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC20pt3->Fill(part->Pt(),deltaPhi);}
857               if((ptcorr>=120.)&&(ptcorr<140.)) {fh2DeltaRC20pt4->Fill(part->Pt(),deltaR); 
858               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC20pt4->Fill(part->Pt(),deltaEta);
859               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC20pt4->Fill(part->Pt(),deltaPhi); }}
860
861           if((centValue>30.)&&(centValue<60.)){
862
863               if((ptcorr>=70.)&&(ptcorr<85.)) {fh2DeltaRC30pt1->Fill(part->Pt(),deltaR);
864               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R))fh2DeltaEtaC30pt1->Fill(part->Pt(),deltaEta);
865               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC30pt1->Fill(part->Pt(),deltaPhi);}
866               if((ptcorr>=85.)&&(ptcorr<100.)) {fh2DeltaRC30pt2->Fill(part->Pt(),deltaR);
867               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC30pt2->Fill(part->Pt(),deltaEta);
868               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC30pt2->Fill(part->Pt(),deltaPhi); }
869               if((ptcorr>=100.)&&(ptcorr<120.)) {fh2DeltaRC30pt3->Fill(part->Pt(),deltaR);
870               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC30pt3->Fill(part->Pt(),deltaEta);
871               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC30pt3->Fill(part->Pt(),deltaPhi);}
872               if((ptcorr>=120.)&&(ptcorr<140.)) {fh2DeltaRC30pt4->Fill(part->Pt(),deltaR); 
873               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC30pt4->Fill(part->Pt(),deltaEta);
874               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC30pt4->Fill(part->Pt(),deltaPhi); }}
875                  
876
877           if(centValue>60.){
878               if((ptcorr>=70.)&&(ptcorr<85.)) {fh2DeltaRC60pt1->Fill(part->Pt(),deltaR);
879               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R))fh2DeltaEtaC60pt1->Fill(part->Pt(),deltaEta);
880               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC60pt1->Fill(part->Pt(),deltaPhi);}
881               if((ptcorr>=85.)&&(ptcorr<100.)) {fh2DeltaRC60pt2->Fill(part->Pt(),deltaR);
882               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC60pt2->Fill(part->Pt(),deltaEta);
883               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC60pt2->Fill(part->Pt(),deltaPhi); }
884               if((ptcorr>=100.)&&(ptcorr<120.)) {fh2DeltaRC60pt3->Fill(part->Pt(),deltaR);
885               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC60pt3->Fill(part->Pt(),deltaEta);
886               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC60pt3->Fill(part->Pt(),deltaPhi);}
887               if((ptcorr>=120.)&&(ptcorr<140.)) {fh2DeltaRC60pt4->Fill(part->Pt(),deltaR); 
888               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC60pt4->Fill(part->Pt(),deltaEta);
889               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC60pt4->Fill(part->Pt(),deltaPhi);}}
890                  
891      } //end of track loop
892      Double_t coronain=rho*TMath::Pi()*(1.-0.8*0.8);
893      Double_t coronaout=rho*TMath::Pi()*(0.6*0.6-0.4*0.4);
894      if(centValue<10.){  
895      fh2SumPtInC10bkg->Fill(ptcorr,coronain/ptbig);       
896      fh2SumPtOutC10bkg->Fill(ptcorr,coronaout/ptbig);           
897      fh2SumPtInC10->Fill(ptcorr,sumPtIn/ptbig);
898      fh2SumPtOutC10->Fill(ptcorr,sumPtOut/ptbig); 
899      fh2SumPtOutC10b->Fill(ptcorr,sumPtOut/ptbig);      
900 }
901      if((centValue>20.)&&(centValue<40.)){             
902      fh2SumPtInC20bkg->Fill(ptcorr,coronain/ptbig);
903      fh2SumPtOutC20bkg->Fill(ptcorr,coronaout/ptbig);    
904      fh2SumPtInC20->Fill(ptcorr,sumPtIn/ptbig);
905      fh2SumPtOutC20->Fill(ptcorr,sumPtOut/ptbig); } 
906      if((centValue>30.)&&(centValue<60.)){
907      fh2SumPtInC30bkg->Fill(ptcorr,coronain/ptbig);
908      fh2SumPtOutC30bkg->Fill(ptcorr,coronaout/ptbig);                    
909      fh2SumPtInC30->Fill(ptcorr,sumPtIn/ptbig);
910      fh2SumPtOutC30->Fill(ptcorr,sumPtOut/ptbig); } 
911      if(centValue>60.){         
912      fh2SumPtInC60bkg->Fill(ptcorr,coronain/ptbig);
913      fh2SumPtOutC60bkg->Fill(ptcorr,coronaout/ptbig);      
914      fh2SumPtInC60->Fill(ptcorr,sumPtIn/ptbig);
915      fh2SumPtOutC60->Fill(ptcorr,sumPtOut/ptbig);}      
916   
917      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
918
919      //tracks up to R=0.8 distant from the jet axis
920      if(fAngStructCloseTracks==1){
921       TList CloseTrackList;
922       Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
923       Double_t difR=0.04;
924       for(Int_t l=0;l<15;l++){
925       Double_t rr=l*0.1+0.1;
926        for(int it = 0;it<nn;++it){
927            AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
928            for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){      
929            AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);  
930            Double_t ptm=part1->Pt();
931            Double_t ptn=part2->Pt();    
932            Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
933                       Rnm=TMath::Sqrt(Rnm);
934                       Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));      
935                       Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
936                       if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
937                                                         down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
938                       if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
939                                                         down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
940                       if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
941                                                          down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
942                       if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
943                         down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
944      }
945     
946      //only jet constituents
947       if(fAngStructCloseTracks==0){
948
949       Double_t difR=0.04;
950       for(Int_t l=0;l<15;l++){
951       Double_t rr=l*0.1+0.1;
952
953     
954       AliAODTrack* part1;
955       AliAODTrack* part2;
956           for(Int_t it=0; it<nTracksGenJet; ++it){
957              part1 = (AliAODTrack*)(genTrackList->At(it));
958            for(Int_t itu=it+1; itu<nTracksGenJet; ++itu){
959              part2 = (AliAODTrack*)(genTrackList->At(itu));
960            Double_t ptm=part1->Pt();
961            Double_t ptn=part2->Pt();    
962            Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
963                       Rnm=TMath::Sqrt(Rnm);
964                       Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
965                       Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR))); 
966                       if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
967                                                         down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
968                       if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
969                                                         down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}  
970                       if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
971                                                          down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
972                       if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
973                         down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
974
975
976    
977
978
979
980
981
982
983
984
985
986
987
988
989    }
990
991
992      //end loop over R=0.4 jets 
993    
994      for(Int_t l=0;l<15;l++){
995      Double_t rr=l*0.1+0.1;
996         if(down1[l]!=0){  
997         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
998         if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
999         if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
1000         if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
1001         if(down2[l]!=0){  
1002         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
1003         if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
1004         if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
1005         if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
1006         if(down3[l]!=0){  
1007         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
1008         if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
1009         if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
1010         if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
1011         if(down4[l]!=0){  
1012         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
1013         if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
1014         if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
1015         if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}
1016
1017     
1018
1019
1020
1021
1022
1023    PostData(1, fOutputList);
1024 }
1025
1026 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
1027 {
1028    // Draw result to the screen
1029    // Called once at the end of the query
1030
1031    if (!GetOutputData(1))
1032    return;
1033 }
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
1046
1047     Int_t iCount = 0;
1048  
1049   Int_t fFilterMask=272;
1050     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1051       AliAODTrack *tr = fAOD->GetTrack(it);
1052       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1053       if(TMath::Abs(tr->Eta())>0.9)continue;
1054       if(tr->Pt()<0.15)continue;
1055       list->Add(tr);
1056       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1057       iCount++;
1058     }
1059   
1060    list->Sort();
1061   return iCount;
1062  
1063 }
1064
1065  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1066
1067     Int_t iCount = 0;
1068  
1069   Int_t fFilterMask=272;
1070     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1071       AliAODTrack *tr = fAOD->GetTrack(it);
1072       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1073       if(TMath::Abs(tr->Eta())>0.9)continue;
1074       if(tr->Pt()<0.15)continue;
1075       Double_t disR=jetbig->DeltaR(tr);
1076       if(disR>0.8)  continue;
1077       list->Add(tr);
1078       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1079       iCount++;
1080     }
1081   
1082    list->Sort();
1083    return iCount;
1084
1085 }
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1098 {
1099
1100    Int_t nInputTracks = 0;
1101
1102    TString jbname(fJetBranchName[1]);
1103    //needs complete event, use jets without background subtraction
1104    for(Int_t i=1; i<=3; ++i){
1105       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1106    }
1107    // use only HI event
1108    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1109    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1110
1111    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1112    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
1113    if(!tmpAODjets){
1114       Printf("Jet branch %s not found", jbname.Data());
1115       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1116       return -1;
1117    }
1118    
1119    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1120       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1121       if(!jet) continue;
1122       TRefArray *trackList = jet->GetRefTracks();
1123       Int_t nTracks = trackList->GetEntriesFast();
1124       nInputTracks += nTracks;
1125       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1126    }
1127    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1128
1129    return nInputTracks;  
1130 }
1131
1132
1133
1134 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1135
1136   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1137   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1138   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1139   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1140   double dphi = mphi-vphi;
1141   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1142   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1143   return dphi;//dphi in [-Pi, Pi]
1144 }
1145
1146
1147