add flexibility for main trigger
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetCorePP.cxx
1
2 // ******************************************
3 // This task computes several jet observables like 
4 // the fraction of energy in inner and outer coronnas,
5 // jet-track correlations,triggered jet shapes and 
6 // correlation strength distribution of particles inside jets.    
7 // Author: lcunquei@cern.ch
8 // *******************************************
9
10
11 /**************************************************************************
12  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
13  *                                                                        *
14  * Author: The ALICE Off-line Project.                                    *
15  * Contributors are mentioned in the code where appropriate.              *
16  *                                                                        *
17  * Permission to use, copy, modify and distribute this software and its   *
18  * documentation strictly for non-commercial purposes is hereby granted   *
19  * without fee, provided that the above copyright notice appears in all   *
20  * copies and that both the copyright notice and this permission notice   *
21  * appear in the supporting documentation. The authors make no claims     *
22  * about the suitability of this software for any purpose. It is          *
23  * provided "as is" without express or implied warranty.                  *
24  **************************************************************************/
25
26
27 #include "TChain.h"
28 #include "TClonesArray.h"
29 #include "TTree.h"
30 #include "TList.h"
31 #include "TMath.h"
32 #include "TH1F.h"
33 #include "TH1D.h"
34 #include "TH2F.h"
35 #include "TH3F.h"
36 #include "THnSparse.h"
37 #include "TCanvas.h"
38 #include "TArrayI.h" 
39 #include "TProfile.h"
40 #include "TFile.h"
41 #include "TKey.h"
42 #include "TRandom3.h"
43
44 #include "AliLog.h"
45
46 #include "AliAnalysisTask.h"
47 #include "AliAnalysisManager.h"
48
49 #include "AliVEvent.h"
50 #include "AliESDEvent.h"
51 #include "AliMCEvent.h"
52 #include "AliESDInputHandler.h"
53 #include "AliCentrality.h"
54 #include "AliAnalysisHelperJetTasks.h"
55 #include "AliInputEventHandler.h"
56 #include "AliAODJetEventBackground.h"
57 #include "AliGenPythiaEventHeader.h"
58 #include "AliAODMCParticle.h"
59 #include "AliMCParticle.h"
60 #include "AliAODEvent.h"
61 #include "AliAODHandler.h"
62 #include "AliAODJet.h"
63 #include "AliVVertex.h"
64 #include "AliVTrack.h"
65 #include "AliAnalysisTaskJetCorePP.h"
66 #include "AliHeader.h" //KF//
67
68 using std::cout;
69 using std::endl;
70
71 ClassImp(AliAnalysisTaskJetCorePP)
72
73 //Filip Krizek 1st March 2013
74
75 //---------------------------------------------------------------------
76 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
77 AliAnalysisTaskSE(),
78 fESD(0x0),
79 fAODIn(0x0),
80 fAODOut(0x0),
81 fAODExtension(0x0),
82 fMcEvent(0x0),
83 fMcHandler(0x0),
84 fJetBranchName(""),
85 fJetBranchNameChargMC(""),
86 fJetBranchNameKine(""),
87 fJetBranchNameFullMC(""),
88 fJetBranchNameBg(""),
89 fJetBranchNameBgChargMC(""),
90 fJetBranchNameBgKine(""),
91 fListJets(0x0),
92 fListJetsGen(0x0),
93 fListJetsGenFull(0x0),
94 fListJetsBg(0x0),
95 fListJetsBgGen(0x0),
96 fNonStdFile(""),
97 fSystem(0), //pp=0  pPb=1
98 fJetParamR(0.4),
99 fBgJetParamR(0.3),
100 fBgMaxJetPt(8.0),
101 fBgConeR(0.4),
102 fOfflineTrgMask(AliVEvent::kAny),
103 fMinContribVtx(1),
104 fVtxZMin(-10.0),
105 fVtxZMax(10.0),
106 fFilterMask(0),
107 fCentMin(0.0),
108 fCentMax(100.0),
109 fJetEtaMin(-0.5),
110 fJetEtaMax(0.5),
111 fTriggerEtaCut(0.9),
112 fTrackEtaCut(0.9),
113 fTrackLowPtCut(0.15),
114 fUseExchContainer(0),
115 fOutputList(0x0),
116 fHistEvtSelection(0x0),
117 fh2Ntriggers(0x0),
118 fHJetSpec(0x0),
119 fHJetSpecSubUeMedian(0x0),
120 fHJetSpecSubUeCone(0x0),
121 fHJetPhiCorr(0x0),
122 fHJetPhiCorrSubUeMedian(0x0),
123 fHJetPhiCorrSubUeCone(0x0),
124 fHJetUeMedian(0x0),
125 fHJetUeCone(0x0),
126 fHRhoUeMedianVsCone(0x0),
127 //fHJetDensity(0x0),
128 //fHJetDensityA4(0x0),
129 fhJetPhi(0x0),
130 fhTriggerPhi(0x0),
131 fhJetEta(0x0),
132 fhTriggerEta(0x0),
133 fhVertexZ(0x0),
134 fhVertexZAccept(0x0),
135 fhContribVtx(0x0),
136 fhContribVtxAccept(0x0),
137 fhDphiTriggerJet(0x0),
138 fhDphiTriggerJetAccept(0x0),
139 fhCentrality(0x0),
140 fhCentralityAccept(0x0),
141 fhNofMultipleTriggers(0x0),
142 fhNofMultipleTriggersCone(0x0),
143 fhNofMultipleTriggersConeLow(0x0),
144 fhNofMultipleTriggersConeHigh(0x0),
145 fhDeltaRMultTriggers(0x0),
146 fhDeltaPhiMultTriggers(0x0),
147 //fHJetPtRaw(0x0),
148 //fHLeadingJetPtRaw(0x0), 
149 //fHDphiVsJetPtAll(0x0), 
150 fhJetPtGenVsJetPtRec(0x0),
151 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
152 fhJetPtGenVsJetPtRecSubUeCone(0x0),
153 fhJetPtGen(0x0), 
154 fhJetPtSubUeMedianGen(0x0), 
155 fhJetPtSubUeConeGen(0x0), 
156 fhJetPtGenChargVsJetPtGenFull(0x0),
157 fhJetPtGenFull(0x0),
158 fh2NtriggersGen(0x0),
159 fHJetSpecGen(0x0),
160 fHJetSpecSubUeMedianGen(0x0),
161 fHJetSpecSubUeConeGen(0x0),
162 fHJetPhiCorrGen(0x0),
163 fHJetPhiCorrSubUeMedianGen(0x0),
164 fHJetPhiCorrSubUeConeGen(0x0),
165 fHJetUeMedianGen(0x0),
166 fHJetUeConeGen(0x0),
167 fhPtTrkTruePrimRec(0x0),
168 fhPtTrkTruePrimGen(0x0),
169 fhPtTrkSecOrFakeRec(0x0),
170 fHRhoUeMedianVsConeGen(0x0),
171 fhEntriesToMedian(0x0),
172 fhEntriesToMedianGen(0x0),
173 fhCellAreaToMedian(0x0),
174 fhCellAreaToMedianGen(0x0),
175 fhNofMultipleTriggersGen(0x0),
176 fhNofMultipleTriggersConeGen(0x0),
177 fhNofMultipleTriggersConeGenLow(0x0),
178 fhNofMultipleTriggersConeGenHigh(0x0),
179 fhDeltaRMultTriggersGen(0x0),
180 fhDeltaPhiMultTriggersGen(0x0),
181 fhNofMultipleTriggersConeGenA(0x0),
182 fhNofMultipleTriggersConeGenALow(0x0),
183 fhNofMultipleTriggersConeGenAHigh(0x0),
184 fhNofMultipleTriggersConeGenB(0x0),
185 fhNofMultipleTriggersConeGenBLow(0x0),
186 fhNofMultipleTriggersConeGenBHigh(0x0),
187 fIsChargedMC(0),
188 fIsKine(0),
189 fIsFullMC(0),
190 faGenIndex(0),
191 faRecIndex(0),
192 fkAcceptance(2.0*TMath::Pi()*1.8),
193 fkDeltaPhiCut(TMath::Pi()-0.8),
194 fh1Xsec(0x0),
195 fh1Trials(0x0),
196 fh1AvgTrials(0x0),
197 fh1PtHard(0x0),
198 fh1PtHardNoW(0x0),  
199 fh1PtHardTrials(0x0),
200 fAvgTrials(1),
201 fHardest(0),
202 fEventNumberRangeLow(0),
203 fEventNumberRangeHigh(99),
204 fTriggerPtRangeLow(0.0),
205 fTriggerPtRangeHigh(10000.0),
206 fFillRespMx(0),
207 fRandom(0x0),
208 fnTrials(1000),
209 fJetFreeAreaFrac(0.5),
210 fnPhi(9),
211 fnEta(2),
212 fEtaSize(0.7),
213 fPhiSize(2*TMath::Pi()/fnPhi),
214 fCellArea(fPhiSize*fEtaSize),
215 fSafetyMargin(1.1),
216 fDoubleBinning(kFALSE)
217 {
218    // default Constructor
219 }
220
221 //---------------------------------------------------------------------
222
223 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
224 AliAnalysisTaskSE(name),
225 fESD(0x0),
226 fAODIn(0x0),
227 fAODOut(0x0),
228 fAODExtension(0x0),
229 fMcEvent(0x0),
230 fMcHandler(0x0),
231 fJetBranchName(""),
232 fJetBranchNameChargMC(""),
233 fJetBranchNameKine(""),
234 fJetBranchNameFullMC(""),
235 fJetBranchNameBg(""),
236 fJetBranchNameBgChargMC(""),
237 fJetBranchNameBgKine(""),
238 fListJets(0x0),
239 fListJetsGen(0x0),
240 fListJetsGenFull(0x0),
241 fListJetsBg(0x0),
242 fListJetsBgGen(0x0),
243 fNonStdFile(""),
244 fSystem(0),  //pp=0   pPb=1
245 fJetParamR(0.4),
246 fBgJetParamR(0.3),
247 fBgMaxJetPt(8.0),
248 fBgConeR(0.4),
249 fOfflineTrgMask(AliVEvent::kAny),
250 fMinContribVtx(1),
251 fVtxZMin(-10.0),
252 fVtxZMax(10.0),
253 fFilterMask(0),
254 fCentMin(0.0),
255 fCentMax(100.0),
256 fJetEtaMin(-0.5),
257 fJetEtaMax(0.5),
258 fTriggerEtaCut(0.9),
259 fTrackEtaCut(0.9),
260 fTrackLowPtCut(0.15),
261 fUseExchContainer(0),
262 fOutputList(0x0),
263 fHistEvtSelection(0x0),
264 fh2Ntriggers(0x0),
265 fHJetSpec(0x0),
266 fHJetSpecSubUeMedian(0x0),
267 fHJetSpecSubUeCone(0x0),
268 fHJetPhiCorr(0x0),
269 fHJetPhiCorrSubUeMedian(0x0),
270 fHJetPhiCorrSubUeCone(0x0),
271 fHJetUeMedian(0x0),
272 fHJetUeCone(0x0),
273 fHRhoUeMedianVsCone(0x0),
274 //fHJetDensity(0x0),
275 //fHJetDensityA4(0x0),
276 fhJetPhi(0x0),
277 fhTriggerPhi(0x0),
278 fhJetEta(0x0),
279 fhTriggerEta(0x0),
280 fhVertexZ(0x0),
281 fhVertexZAccept(0x0),
282 fhContribVtx(0x0),
283 fhContribVtxAccept(0x0),
284 fhDphiTriggerJet(0x0),
285 fhDphiTriggerJetAccept(0x0),
286 fhCentrality(0x0),
287 fhCentralityAccept(0x0),
288 fhNofMultipleTriggers(0x0),
289 fhNofMultipleTriggersCone(0x0),
290 fhNofMultipleTriggersConeLow(0x0),
291 fhNofMultipleTriggersConeHigh(0x0),
292 fhDeltaRMultTriggers(0x0),
293 fhDeltaPhiMultTriggers(0x0),
294 //fHJetPtRaw(0x0),
295 //fHLeadingJetPtRaw(0x0), 
296 //fHDphiVsJetPtAll(0x0), 
297 fhJetPtGenVsJetPtRec(0x0),
298 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
299 fhJetPtGenVsJetPtRecSubUeCone(0x0),
300 fhJetPtGen(0x0),
301 fhJetPtSubUeMedianGen(0x0), 
302 fhJetPtSubUeConeGen(0x0), 
303 fhJetPtGenChargVsJetPtGenFull(0x0),
304 fhJetPtGenFull(0x0),
305 fh2NtriggersGen(0x0),
306 fHJetSpecGen(0x0),
307 fHJetSpecSubUeMedianGen(0x0),
308 fHJetSpecSubUeConeGen(0x0),
309 fHJetPhiCorrGen(0x0),
310 fHJetPhiCorrSubUeMedianGen(0x0),
311 fHJetPhiCorrSubUeConeGen(0x0),
312 fHJetUeMedianGen(0x0),
313 fHJetUeConeGen(0x0),
314 fhPtTrkTruePrimRec(0x0),
315 fhPtTrkTruePrimGen(0x0),
316 fhPtTrkSecOrFakeRec(0x0),
317 fHRhoUeMedianVsConeGen(0x0),
318 fhEntriesToMedian(0x0),
319 fhEntriesToMedianGen(0x0),
320 fhCellAreaToMedian(0x0),
321 fhCellAreaToMedianGen(0x0),
322 fhNofMultipleTriggersGen(0x0),
323 fhNofMultipleTriggersConeGen(0x0),
324 fhNofMultipleTriggersConeGenLow(0x0),
325 fhNofMultipleTriggersConeGenHigh(0x0),
326 fhDeltaRMultTriggersGen(0x0),
327 fhDeltaPhiMultTriggersGen(0x0),
328 fhNofMultipleTriggersConeGenA(0x0),
329 fhNofMultipleTriggersConeGenALow(0x0),
330 fhNofMultipleTriggersConeGenAHigh(0x0),
331 fhNofMultipleTriggersConeGenB(0x0),
332 fhNofMultipleTriggersConeGenBLow(0x0),
333 fhNofMultipleTriggersConeGenBHigh(0x0),
334 fIsChargedMC(0),
335 fIsKine(0),
336 fIsFullMC(0),
337 faGenIndex(0),
338 faRecIndex(0),
339 fkAcceptance(2.0*TMath::Pi()*1.8),
340 fkDeltaPhiCut(TMath::Pi()-0.8),
341 fh1Xsec(0x0),
342 fh1Trials(0x0), 
343 fh1AvgTrials(0x0),
344 fh1PtHard(0x0),
345 fh1PtHardNoW(0x0),  
346 fh1PtHardTrials(0x0),
347 fAvgTrials(1),
348 fHardest(0),
349 fEventNumberRangeLow(0),
350 fEventNumberRangeHigh(99),
351 fTriggerPtRangeLow(0.0),
352 fTriggerPtRangeHigh(10000.0),
353 fFillRespMx(0),
354 fRandom(0x0),
355 fnTrials(1000),
356 fJetFreeAreaFrac(0.5),
357 fnPhi(9),
358 fnEta(2),
359 fEtaSize(0.7),
360 fPhiSize(2*TMath::Pi()/fnPhi),
361 fCellArea(fPhiSize*fEtaSize),
362 fSafetyMargin(1.1),
363 fDoubleBinning(kFALSE)
364 {
365 // Constructor
366    DefineOutput(1, TList::Class());
367
368    TString dummy(name);
369    if(dummy.Contains("KINE")){
370       DefineInput(1, TClonesArray::Class());
371       DefineInput(2, TClonesArray::Class());
372    }
373 }
374
375 //--------------------------------------------------------------
376 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
377 AliAnalysisTaskSE(a.GetName()),
378 fESD(a.fESD),
379 fAODIn(a.fAODIn),
380 fAODOut(a.fAODOut),
381 fAODExtension(a.fAODExtension),
382 fMcEvent(a.fMcEvent),
383 fMcHandler(a.fMcHandler),
384 fJetBranchName(a.fJetBranchName),
385 fJetBranchNameChargMC(a.fJetBranchNameChargMC),
386 fJetBranchNameKine(a.fJetBranchNameKine),
387 fJetBranchNameFullMC(a.fJetBranchNameFullMC),
388 fJetBranchNameBg(a.fJetBranchNameBg),
389 fJetBranchNameBgChargMC(a.fJetBranchNameBgChargMC),
390 fJetBranchNameBgKine(a.fJetBranchNameBgKine),
391 fListJets(a.fListJets),
392 fListJetsGen(a.fListJetsGen),
393 fListJetsGenFull(a.fListJetsGenFull),
394 fListJetsBg(a.fListJetsBg),
395 fListJetsBgGen(a.fListJetsBgGen),
396 fNonStdFile(a.fNonStdFile),
397 fSystem(a.fSystem),  
398 fJetParamR(a.fJetParamR),
399 fBgJetParamR(a.fBgJetParamR),
400 fBgMaxJetPt(a.fBgMaxJetPt),
401 fBgConeR(a.fBgConeR),
402 fOfflineTrgMask(a.fOfflineTrgMask),
403 fMinContribVtx(a.fMinContribVtx),
404 fVtxZMin(a.fVtxZMin),
405 fVtxZMax(a.fVtxZMax),
406 fFilterMask(a.fFilterMask),
407 fCentMin(a.fCentMin),
408 fCentMax(a.fCentMax),
409 fJetEtaMin(a.fJetEtaMin),
410 fJetEtaMax(a.fJetEtaMax),
411 fTriggerEtaCut(a.fTriggerEtaCut),
412 fTrackEtaCut(a.fTrackEtaCut),
413 fTrackLowPtCut(a.fTrackLowPtCut),
414 fUseExchContainer(a.fUseExchContainer),
415 fOutputList(a.fOutputList),
416 fHistEvtSelection(a.fHistEvtSelection),
417 fh2Ntriggers(a.fh2Ntriggers),
418 fHJetSpec(a.fHJetSpec),
419 fHJetSpecSubUeMedian(a.fHJetSpecSubUeMedian),
420 fHJetSpecSubUeCone(a.fHJetSpecSubUeCone),
421 fHJetPhiCorr(a.fHJetPhiCorr),
422 fHJetPhiCorrSubUeMedian(a.fHJetPhiCorrSubUeMedian),
423 fHJetPhiCorrSubUeCone(a.fHJetPhiCorrSubUeCone),
424 fHJetUeMedian(a.fHJetUeMedian),
425 fHJetUeCone(a.fHJetUeCone),
426 fHRhoUeMedianVsCone(a.fHRhoUeMedianVsCone), 
427 //fHJetDensity(a.fHJetDensity),
428 //fHJetDensityA4(a.fHJetDensityA4),
429 fhJetPhi(a.fhJetPhi),
430 fhTriggerPhi(a.fhTriggerPhi),
431 fhJetEta(a.fhJetEta),
432 fhTriggerEta(a.fhTriggerEta),
433 fhVertexZ(a.fhVertexZ),
434 fhVertexZAccept(a.fhVertexZAccept),
435 fhContribVtx(a.fhContribVtx),
436 fhContribVtxAccept(a.fhContribVtxAccept),
437 fhDphiTriggerJet(a.fhDphiTriggerJet),
438 fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
439 fhCentrality(a.fhCentrality),
440 fhCentralityAccept(a.fhCentralityAccept),
441 fhNofMultipleTriggers(a.fhNofMultipleTriggers),
442 fhNofMultipleTriggersCone(a.fhNofMultipleTriggersCone),
443 fhNofMultipleTriggersConeLow(a.fhNofMultipleTriggersConeLow),
444 fhNofMultipleTriggersConeHigh(a.fhNofMultipleTriggersConeHigh),
445 fhDeltaRMultTriggers(a.fhDeltaRMultTriggers),
446 fhDeltaPhiMultTriggers(a.fhDeltaPhiMultTriggers),
447 //fHJetPtRaw(a.fHJetPtRaw),
448 //fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
449 //fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
450 fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
451 fhJetPtGenVsJetPtRecSubUeMedian(a.fhJetPtGenVsJetPtRecSubUeMedian),
452 fhJetPtGenVsJetPtRecSubUeCone(a.fhJetPtGenVsJetPtRecSubUeCone),
453 fhJetPtGen(a.fhJetPtGen),
454 fhJetPtSubUeMedianGen(a.fhJetPtSubUeMedianGen), 
455 fhJetPtSubUeConeGen(a.fhJetPtSubUeConeGen), 
456 fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
457 fhJetPtGenFull(a.fhJetPtGenFull),
458 fh2NtriggersGen(a.fh2NtriggersGen),
459 fHJetSpecGen(a.fHJetSpecGen),
460 fHJetSpecSubUeMedianGen(a.fHJetSpecSubUeMedianGen),
461 fHJetSpecSubUeConeGen(a.fHJetSpecSubUeConeGen),
462 fHJetPhiCorrGen(a.fHJetPhiCorrGen),
463 fHJetPhiCorrSubUeMedianGen(a.fHJetPhiCorrSubUeMedianGen),
464 fHJetPhiCorrSubUeConeGen(a.fHJetPhiCorrSubUeConeGen),
465 fHJetUeMedianGen(a.fHJetUeMedianGen),
466 fHJetUeConeGen(a.fHJetUeConeGen),
467 fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
468 fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
469 fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
470 fHRhoUeMedianVsConeGen(a.fHRhoUeMedianVsConeGen),
471 fhEntriesToMedian(a.fhEntriesToMedian),
472 fhEntriesToMedianGen(a.fhEntriesToMedianGen),
473 fhCellAreaToMedian(a.fhCellAreaToMedian),
474 fhCellAreaToMedianGen(a.fhCellAreaToMedianGen),
475 fhNofMultipleTriggersGen(a.fhNofMultipleTriggersGen),
476 fhNofMultipleTriggersConeGen(a.fhNofMultipleTriggersConeGen),
477 fhNofMultipleTriggersConeGenLow(a.fhNofMultipleTriggersConeGenLow),
478 fhNofMultipleTriggersConeGenHigh(a.fhNofMultipleTriggersConeGenHigh),
479 fhDeltaRMultTriggersGen(a.fhDeltaRMultTriggersGen),
480 fhDeltaPhiMultTriggersGen(a.fhDeltaPhiMultTriggersGen),
481 fhNofMultipleTriggersConeGenA(a.fhNofMultipleTriggersConeGenA),
482 fhNofMultipleTriggersConeGenALow(a.fhNofMultipleTriggersConeGenALow),
483 fhNofMultipleTriggersConeGenAHigh(a.fhNofMultipleTriggersConeGenAHigh),
484 fhNofMultipleTriggersConeGenB(a.fhNofMultipleTriggersConeGenB),
485 fhNofMultipleTriggersConeGenBLow(a.fhNofMultipleTriggersConeGenBLow),
486 fhNofMultipleTriggersConeGenBHigh(a.fhNofMultipleTriggersConeGenBHigh),
487 fIsChargedMC(a.fIsChargedMC),
488 fIsKine(a.fIsKine),
489 fIsFullMC(a.fIsFullMC),
490 faGenIndex(a.faGenIndex),
491 faRecIndex(a.faRecIndex),
492 fkAcceptance(a.fkAcceptance),
493 fkDeltaPhiCut(a.fkDeltaPhiCut),
494 fh1Xsec(a.fh1Xsec),
495 fh1Trials(a.fh1Trials),
496 fh1AvgTrials(a.fh1AvgTrials),
497 fh1PtHard(a.fh1PtHard),
498 fh1PtHardNoW(a.fh1PtHardNoW),  
499 fh1PtHardTrials(a.fh1PtHardTrials),
500 fAvgTrials(a.fAvgTrials),
501 fHardest(a.fHardest),
502 fEventNumberRangeLow(a.fEventNumberRangeLow),
503 fEventNumberRangeHigh(a.fEventNumberRangeHigh),
504 fTriggerPtRangeLow(a.fTriggerPtRangeLow),
505 fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
506 fFillRespMx(a.fFillRespMx),
507 fRandom(a.fRandom),
508 fnTrials(a.fnTrials),
509 fJetFreeAreaFrac(a.fJetFreeAreaFrac),
510 fnPhi(a.fnPhi),
511 fnEta(a.fnEta),
512 fEtaSize(a.fEtaSize),
513 fPhiSize(a.fPhiSize),
514 fCellArea(a.fCellArea),
515 fSafetyMargin(a.fSafetyMargin),
516 fDoubleBinning(a.fDoubleBinning)
517 {
518    //Copy Constructor
519    fRandom->SetSeed(0);
520 }
521 //--------------------------------------------------------------
522
523 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
524   // assignment operator
525   this->~AliAnalysisTaskJetCorePP();
526   new(this) AliAnalysisTaskJetCorePP(a);
527   return *this;
528 }
529 //--------------------------------------------------------------
530
531 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
532 {
533    //Destructor 
534    delete fListJets;
535    delete fListJetsGen;
536    delete fListJetsGenFull;
537    delete fListJetsBg;
538    delete fListJetsBgGen;
539    delete fOutputList; // ????
540    delete fRandom;
541 }
542
543 //--------------------------------------------------------------
544
545
546 Bool_t AliAnalysisTaskJetCorePP::Notify()
547 {
548    //Implemented Notify() to read the cross sections
549    //and number of trials from pyxsec.root
550    //inspired by AliAnalysisTaskJetSpectrum2::Notify()
551    if(!(fIsChargedMC || fIsKine)) return kFALSE; 
552    Float_t xsection = 0;
553    Float_t trials  = 1;
554    fAvgTrials = 1;
555
556    if(fIsChargedMC){ 
557       fESD = dynamic_cast<AliESDEvent*>(InputEvent());
558       if(!fESD){
559          if(fDebug>1) AliError("ESD not available");
560          fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
561       } 
562  
563       fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
564
565
566       if(fNonStdFile.Length()!=0){
567          // case that we have an AOD extension we can fetch the jets from the extended output
568          AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
569          fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
570          if(!fAODExtension){
571             if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
572          } 
573       }
574  
575       TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
576
577       if(tree){
578          TFile *curfile = tree->GetCurrentFile();
579          if(!curfile) {
580             Error("Notify","No current file");
581             return kFALSE;
582          }
583          if(!fh1Xsec || !fh1Trials){
584             Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
585             return kFALSE;
586          }
587          AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
588          fh1Xsec->Fill("<#sigma>",xsection);
589          // construct a poor man average trials
590          Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
591          if(trials>=nEntries && nEntries>0.) fAvgTrials = trials/nEntries;
592          fh1Trials->Fill("#sum{ntrials}",trials);
593       }  
594
595       if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
596    }
597
598
599    return kTRUE;
600 }
601 //--------------------------------------------------------------
602
603 void AliAnalysisTaskJetCorePP::Init()
604 {
605    // check for jet branches
606    if(fJetBranchNameKine.Length()==0){
607       if(!strlen(fJetBranchName.Data())){
608          AliError("Jet branch name not set.");
609       }
610    }
611    fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); 
612
613
614 }
615
616 //--------------------------------------------------------------
617
618 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
619 {
620    // Create histograms and initilize variables
621    fSafetyMargin = fBgConeR*fBgConeR /(fBgJetParamR*fBgJetParamR);
622
623    // Called once
624    fListJets   = new TList();  //reconstructed level antikt jets
625    fListJetsBg = new TList();  //reconstructed jets to be removed from bg
626
627    fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
628    fIsKine      = (fJetBranchNameKine.Length()>0)    ? kTRUE : kFALSE;
629    fIsFullMC    = (fJetBranchNameFullMC.Length()>0)  ? kTRUE : kFALSE;
630
631    fRandom = new TRandom3(0);
632
633    if(fIsChargedMC || fIsKine){   //full MC or pure kine
634       fListJetsGen   = new TList(); //generator level charged antikt jets
635       fListJetsBgGen = new TList(); //generator level jets to be removed from bg 
636
637       if(fIsFullMC)
638          fListJetsGenFull = new TList(); //generator level full jets
639    }
640    OpenFile(1);
641    if(!fOutputList) fOutputList = new TList();
642    fOutputList->SetOwner(kTRUE);
643
644    Bool_t oldStatus = TH1::AddDirectoryStatus();
645    TH1::AddDirectory(kFALSE);
646
647    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
648    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
649    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
650    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
651    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
652    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
653    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
654    
655    fOutputList->Add(fHistEvtSelection);
656
657    Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
658    //trigger pt spectrum (reconstructed) 
659    fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
660                              nBinsCentrality,0.0,100.0,50,0.0,50.0);
661    if(!fIsKine) fOutputList->Add(fh2Ntriggers);
662
663    Int_t bw = (fDoubleBinning==0) ? 1 : 2; //make larger bin width
664
665    //Centrality, A, pTjet, pTtrigg, dphi
666    const Int_t dimSpec   = 5;
667    const Int_t nBinsSpec[dimSpec]     = {nBinsCentrality,  50, bw*110,  50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
668    const Double_t lowBinSpec[dimSpec] = {0.0,             0.0,    -20, 0.0, fkDeltaPhiCut};
669    const Double_t hiBinSpec[dimSpec]  = {100.0,           1.0,  200.0,50.0, TMath::Pi()};
670    fHJetSpec = new THnSparseF("fHJetSpec",
671                    "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
672                    dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
673    if(!fIsKine) fOutputList->Add(fHJetSpec);  
674
675    //background estimated as  median of kT jets 
676    fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian");
677    fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
678    if(!fIsKine) fOutputList->Add(fHJetSpecSubUeMedian); 
679    //background estimated as weighted  median of kT jets  ala Cone
680    fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone");
681    fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
682    if(!fIsKine) fOutputList->Add(fHJetSpecSubUeCone); 
683
684
685    //A, pTjet, pTtrigg, dphi
686    const Int_t dimCor = 4;
687    const Int_t nBinsCor[dimCor]     = {50,  110,   50,      100};
688    const Double_t lowBinCor[dimCor] = {0.0, -20,  0.0, -0.5*TMath::Pi()};
689    const Double_t hiBinCor[dimCor]  = {1.0, 200, 50.0, 1.5*TMath::Pi()};
690    fHJetPhiCorr = new THnSparseF("fHJetPhiCorr",
691                                  "Recoil jet spectrum [A,pTjet,pTtrig,dphi]",
692                                  dimCor,nBinsCor,lowBinCor,hiBinCor);
693
694    if(!fIsKine) fOutputList->Add(fHJetPhiCorr); // Dphi distribution jet-triger
695
696    //background estimated as  median of kT jets 
697    fHJetPhiCorrSubUeMedian = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeMedian");
698    fHJetPhiCorrSubUeMedian->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
699    if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeMedian); 
700    //background estimated as weighted  median of kT jets  ala Cone
701    fHJetPhiCorrSubUeCone = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeCone");
702    fHJetPhiCorrSubUeCone->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
703    if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeCone); 
704
705
706    //------------------- HISTOS FOR DIAGNOSTIC ----------------------
707    //A, pTjet, pTjet-pTUe, pTUe, rhoUe     bg estimate from kT median
708    const Int_t    dimSpecMed   = 5;
709    const Int_t    nBinsSpecMed[dimSpecMed]  = {25,     50,    60,    20,   20};
710    const Double_t lowBinSpecMed[dimSpecMed] = {0.0,   0.0, -20.0,   0.0,  0.0};
711    const Double_t hiBinSpecMed[dimSpecMed]  = {1.0, 100.0, 100.0,  10.0, 20.0};
712    fHJetUeMedian = new THnSparseF("fHJetUeMedian",
713                    "Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]",
714                    dimSpecMed, nBinsSpecMed, lowBinSpecMed, hiBinSpecMed);
715    if(!fIsKine) fOutputList->Add(fHJetUeMedian);  
716    
717    //A, pTjet, pTjet-pTUe, pTUe, rhoUe     bg estimate from kT median Cone
718    fHJetUeCone = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeCone");
719    fHJetUeCone->SetTitle("Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]");
720    if(!fIsKine) fOutputList->Add(fHJetUeCone); 
721
722    //rho bacground reconstructed data
723    const Int_t    dimRho   = 2;
724    const Int_t    nBinsRho[dimRho]  = {50  ,   50};
725    const Double_t lowBinRho[dimRho] = {0.0  , 0.0};
726    const Double_t hiBinRho[dimRho]  = {20.0 , 20.0};
727
728    fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",  
729                                       dimRho, nBinsRho, lowBinRho, hiBinRho);
730    if(!fIsKine) fOutputList->Add(fHRhoUeMedianVsCone);
731
732    //Jet number density histos [Trk Mult, jet density, pT trigger]
733    /*const Int_t    dimJetDens   = 3;
734    const Int_t    nBinsJetDens[dimJetDens]  = {100,   100,   10};
735    const Double_t lowBinJetDens[dimJetDens] = {0.0,   0.0,  0.0};
736    const Double_t hiBinJetDens[dimJetDens]  = {500.0, 5.0, 50.0};
737
738    fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
739                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
740
741    fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
742                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
743
744    fOutputList->Add(fHJetDensity);
745    fOutputList->Add(fHJetDensityA4);
746    */      
747
748    //inclusive azimuthal and pseudorapidity histograms
749    fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
750                         50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
751    fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
752                         50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
753    fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
754                         50,0, 100, 40,-0.9,0.9);
755    fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
756                         50, 0, 50, 40,-0.9,0.9);
757
758    fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);  
759    fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);  
760    fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);   
761    fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);   
762    fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi()); 
763    fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi()); 
764    fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
765    fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
766    fhEntriesToMedian = new TH1D("fhEntriesToMedian","fhEntriesToMedian",30,0,30);
767    fhCellAreaToMedian =  new TH1D("fhCellAreaToMedian", "fhCellAreaToMedian", 75,0,1.5);
768
769    fhNofMultipleTriggers = new TH1D("fhNofMultipleTriggers","fhNofMultipleTriggers",20,0,20);
770    fhNofMultipleTriggersCone = new TH1D("fhNofMultipleTriggersCone","fhNofMultipleTriggersCone R<0.4",20,0,20);
771    fhNofMultipleTriggersConeLow  = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConepTa15to20");
772    fhNofMultipleTriggersConeHigh = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConepTa20to50");
773    fhDeltaRMultTriggers = new  TH1D("fhDeltaRMultTriggers","fhDeltaRMultTriggers", 100,0,4);  
774    fhDeltaPhiMultTriggers = new  TH1D("fhDeltaPhiMultTriggers","fhDeltaPhiRultTriggers", 100,-TMath::Pi(),TMath::Pi());  
775
776    if(!fIsKine){
777       fOutputList->Add(fhJetPhi);
778       fOutputList->Add(fhTriggerPhi);
779       fOutputList->Add(fhJetEta);
780       fOutputList->Add(fhTriggerEta);
781    }
782    fOutputList->Add(fhVertexZ);    
783    fOutputList->Add(fhVertexZAccept);   
784    if(!fIsKine){
785       fOutputList->Add(fhContribVtx); 
786       fOutputList->Add(fhContribVtxAccept); 
787       fOutputList->Add(fhDphiTriggerJet);
788       fOutputList->Add(fhDphiTriggerJetAccept);
789       fOutputList->Add(fhCentrality); 
790       fOutputList->Add(fhCentralityAccept);
791       fOutputList->Add(fhEntriesToMedian);
792       fOutputList->Add(fhCellAreaToMedian);
793       fOutputList->Add(fhNofMultipleTriggers);
794       fOutputList->Add(fhNofMultipleTriggersCone);
795       fOutputList->Add(fhNofMultipleTriggersConeLow);
796       fOutputList->Add(fhNofMultipleTriggersConeHigh);
797       fOutputList->Add(fhDeltaRMultTriggers);
798       fOutputList->Add(fhDeltaPhiMultTriggers);
799    }
800    // raw spectra of INCLUSIVE jets  
801    //Centrality, pTjet, A
802    /*const Int_t dimRaw   = 3;
803    const Int_t nBinsRaw[dimRaw]     = {nBinsCentrality,  50,   100};
804    const Double_t lowBinRaw[dimRaw] = {0.0,             0.0,   0.0};
805    const Double_t hiBinRaw[dimRaw]  = {100.0,           100,   1.0};
806    fHJetPtRaw = new THnSparseF("fHJetPtRaw",
807                                 "Incl. jet spectrum [cent,pTjet,A]",
808                                 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
809    fOutputList->Add(fHJetPtRaw);  
810
811    // raw spectra of LEADING jets  
812    //Centrality, pTjet, A
813    fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
814                                 "Leading jet spectrum [cent,pTjet,A]",
815                                 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
816    fOutputList->Add(fHLeadingJetPtRaw);  
817
818    // Dphi versus pT jet 
819    //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg 
820    const Int_t dimDp   = 4;
821    const Int_t nBinsDp[dimDp]     = {nBinsCentrality,  50,     50,    50};
822    const Double_t lowBinDp[dimDp] = {0.0,       -TMath::Pi(),   0.0,   0.0};
823    const Double_t hiBinDp[dimDp]  = {100.0,      TMath::Pi(), 100.0, 100.0};
824    fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
825                                 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
826                                 dimDp,nBinsDp,lowBinDp,hiBinDp);
827    fOutputList->Add(fHDphiVsJetPtAll);  
828    */
829
830    //analyze MC generator level 
831    if(fIsChargedMC || fIsKine){   
832       if(fFillRespMx){
833          //Fill response matrix only once 
834          fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", bw*100,0,200, bw*100,0,200); 
835          fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC charg jet pt spectrum versus rec charged jet pt spectrum
836          //....
837          fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", bw*110,-20,200, bw*110,-20,200); 
838          fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr 
839          //....
840          fhJetPtGenVsJetPtRecSubUeCone=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCone");
841          fhJetPtGenVsJetPtRecSubUeCone->SetTitle("fhJetPtGenVsJetPtRecSubUeCone");
842          fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCone); // with weighted kT median bg subtr 
843          //....
844          fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",bw*100,0,200); //MC generator charged jet pt spectrum
845          fOutputList->Add(fhJetPtGen);  
846          //....
847          fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",bw*110,-20,200); 
848          fOutputList->Add(fhJetPtSubUeMedianGen);  // with kT median bg subtr
849          //....
850          fhJetPtSubUeConeGen = (TH1D*) fhJetPtSubUeMedianGen->Clone("fhJetPtSubUeConeGen");
851          fOutputList->Add(fhJetPtSubUeConeGen); // with weighted kT median bg subtr
852       
853          //....
854          if(fIsFullMC){
855             fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 100,0,200, 100,0,200);
856             fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen full MC jet pt versus gen charged jet MC pt
857          //....
858             fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",100,0,200); //MC generator full jet pt spectrum
859             fOutputList->Add(fhJetPtGenFull); 
860          }
861       }
862       //....
863       fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
864       fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
865       fOutputList->Add(fh2NtriggersGen);
866
867       //Centrality, A, pT, pTtrigg
868       fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
869       fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
870       fOutputList->Add(fHJetSpecGen); 
871
872       fHJetSpecSubUeMedianGen = (THnSparseF*)  fHJetSpecSubUeMedian->Clone("fHJetSpecSubUeMedianGen");
873       fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
874       fOutputList->Add(fHJetSpecSubUeMedianGen);  
875
876       fHJetSpecSubUeConeGen =  (THnSparseF*) fHJetSpecSubUeCone->Clone("fHJetSpecSubUeConeGen"); 
877       fHJetSpecSubUeConeGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCone->GetTitle()));
878       fOutputList->Add(fHJetSpecSubUeConeGen);
879       //---
880       fHJetPhiCorrGen = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrGen");
881       fHJetPhiCorrGen->SetTitle(Form("%s Gen MC",fHJetPhiCorr->GetTitle()));
882       fOutputList->Add(fHJetPhiCorrGen); 
883
884       fHJetPhiCorrSubUeMedianGen = (THnSparseF*) fHJetPhiCorrSubUeMedian->Clone("fHJetPhiCorrSubUeMedianGen");
885       fHJetPhiCorrSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetPhiCorrSubUeMedian->GetTitle()));
886       fOutputList->Add(fHJetPhiCorrSubUeMedianGen); 
887  
888       fHJetPhiCorrSubUeConeGen = (THnSparseF*) fHJetPhiCorrSubUeCone->Clone("fHJetPhiCorrSubUeConeGen");
889       fHJetPhiCorrSubUeConeGen->SetTitle(Form("%s Gen MC", fHJetPhiCorrSubUeCone->GetTitle()));
890       fOutputList->Add(fHJetPhiCorrSubUeConeGen);
891
892       //---
893       fHJetUeMedianGen  =  (THnSparseF*) fHJetUeMedian->Clone("fHJetUeMedianGen");  
894       fHJetUeMedianGen->SetTitle(Form("%s Gen MC", fHJetUeMedian->GetTitle()));
895       fOutputList->Add(fHJetUeMedianGen);
896
897       fHJetUeConeGen =  (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen"); 
898       fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle())); 
899       fOutputList->Add(fHJetUeConeGen);
900
901       if(fFillRespMx){
902          //track efficiency/contamination histograms  eta versus pT
903          fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
904          fOutputList->Add(fhPtTrkTruePrimRec); 
905       
906          fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
907          fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");    
908          fOutputList->Add(fhPtTrkTruePrimGen);
909       
910          fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");    
911          fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");    
912          fOutputList->Add(fhPtTrkSecOrFakeRec);
913       }
914
915       fHRhoUeMedianVsConeGen = (THnSparseF*) fHRhoUeMedianVsCone->Clone("hRhoUeMedianVsConeGen");
916       fHRhoUeMedianVsConeGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCone->GetTitle())); 
917       fOutputList->Add(fHRhoUeMedianVsConeGen);
918
919       fhEntriesToMedianGen = (TH1D*) fhEntriesToMedian->Clone("fhEntriesToMedianGen");
920       fhEntriesToMedianGen->SetTitle(Form("%s Gen MC", fhEntriesToMedian->GetTitle())); 
921       fOutputList->Add(fhEntriesToMedianGen);
922
923       fhCellAreaToMedianGen  = (TH1D*) fhCellAreaToMedian->Clone("fhCellAreaToMedianGen");
924       fhCellAreaToMedianGen->SetTitle(Form("%s Gen MC", fhCellAreaToMedian->GetTitle())); 
925       fOutputList->Add(fhCellAreaToMedianGen);
926
927       fhNofMultipleTriggersGen = (TH1D*) fhNofMultipleTriggers->Clone("fhNofMultipleTriggersGen"); 
928       fOutputList->Add(fhNofMultipleTriggersGen);
929
930       fhNofMultipleTriggersConeGen = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGen"); 
931       fOutputList->Add(fhNofMultipleTriggersConeGen);
932
933       fhNofMultipleTriggersConeGenLow = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGenpta15to20"); 
934       fOutputList->Add(fhNofMultipleTriggersConeGenLow);
935
936       fhNofMultipleTriggersConeGenHigh = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGenpta20to50"); 
937       fOutputList->Add(fhNofMultipleTriggersConeGenHigh);
938
939
940
941       fhDeltaRMultTriggersGen  = (TH1D*) fhDeltaRMultTriggers->Clone("fhDeltaRMultTriggersGen");
942       fOutputList->Add(fhDeltaRMultTriggersGen);
943
944       fhDeltaPhiMultTriggersGen  = (TH1D*) fhDeltaPhiMultTriggers->Clone("fhDeltaPhiMultTriggersGen");
945       fOutputList->Add(fhDeltaPhiMultTriggersGen);
946
947
948       fhNofMultipleTriggersConeGenA = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10"); 
949       fOutputList->Add(fhNofMultipleTriggersConeGenA);
950
951       fhNofMultipleTriggersConeGenALow = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10pTa15to20"); 
952       fOutputList->Add(fhNofMultipleTriggersConeGenALow);
953
954       fhNofMultipleTriggersConeGenAHigh = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10pTa20to50"); 
955       fOutputList->Add(fhNofMultipleTriggersConeGenAHigh);
956
957       fhNofMultipleTriggersConeGenB = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5"); 
958       fOutputList->Add(fhNofMultipleTriggersConeGenB);
959
960       fhNofMultipleTriggersConeGenBLow = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5pTa15to20"); 
961       fOutputList->Add(fhNofMultipleTriggersConeGenBLow);
962
963       fhNofMultipleTriggersConeGenBHigh = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5pTa20to50"); 
964       fOutputList->Add(fhNofMultipleTriggersConeGenBHigh);
965
966
967
968    }
969    //-------------------------------------
970    //     pythia histograms
971    const Int_t nBinPt = 150;
972    Double_t binLimitsPt[nBinPt+1];
973    for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
974       if(iPt == 0){
975          binLimitsPt[iPt] = -50.0;
976       }else{// 1.0
977          binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
978       }
979    }
980    
981    fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
982    fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
983    fOutputList->Add(fh1Xsec);
984    fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
985    fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
986    fOutputList->Add(fh1Trials);
987    fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
988    fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
989    fOutputList->Add(fh1AvgTrials);
990    fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
991    fOutputList->Add(fh1PtHard);
992    fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
993    fOutputList->Add(fh1PtHardNoW);
994    fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
995    fOutputList->Add(fh1PtHardTrials);      
996    
997
998    // =========== Switch on Sumw2 for all histos ===========
999    for(Int_t i=0; i<fOutputList->GetEntries(); i++){
1000       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
1001       if(h1){
1002          h1->Sumw2();
1003          continue;
1004       }
1005       THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
1006       if(hn){
1007          hn->Sumw2();
1008       }   
1009    }
1010    TH1::AddDirectory(oldStatus);
1011
1012    PostData(1, fOutputList);
1013 }
1014
1015 //--------------------------------------------------------------------
1016
1017 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
1018 {
1019    //User Exec
1020    
1021
1022    //Event loop
1023    Double_t eventW  = 1.0;
1024    Double_t ptHard  = 0.0;
1025    Double_t nTrials = 1.0; // Trials for MC trigger
1026    if(fIsChargedMC || fIsKine) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials); 
1027
1028    if(TMath::Abs((Float_t) fJetParamR)<0.00001){
1029       AliError("ANTIKT Cone radius is set to zero.");  
1030       return;
1031    }
1032
1033    if(TMath::Abs((Float_t) fBgJetParamR)<0.00001){
1034       AliError("ANTIKT Cone radius is set to zero.");  
1035       return;
1036    }
1037
1038    if(!fIsKine){ //real data or full MC
1039       if(!strlen(fJetBranchName.Data())){
1040          AliError("Name of jet branch not set.");
1041          return;
1042       }
1043       if(!strlen(fJetBranchNameBg.Data())){
1044          AliError("Name of jet bg branch not set.");
1045          return;
1046       }
1047    }else{  //Kine
1048        if(!strlen(fJetBranchNameBgKine.Data())){
1049          AliError("Name of jet bg branch for kine not set.");
1050          return;
1051        }
1052
1053       Init(); 
1054       if(fMcHandler){
1055          fMcEvent = fMcHandler->MCEvent(); 
1056       }else{
1057          if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcHandler=NULL\n");
1058          PostData(1, fOutputList);
1059          return;
1060       } 
1061       if(!fMcEvent){
1062          if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcEvent=NULL \n");
1063          PostData(1, fOutputList);
1064          return;
1065       }
1066  
1067       Float_t xsection = 0;
1068       Float_t trials  = 0;
1069
1070       AliGenPythiaEventHeader *genPH =
1071          dynamic_cast<AliGenPythiaEventHeader*> (fMcEvent->GenEventHeader()); 
1072       if(genPH){
1073          xsection = genPH->GetXsection();
1074          trials   = genPH->Trials();
1075          ptHard   = genPH->GetPtHard();
1076       }
1077       fh1Xsec->Fill("<#sigma>",xsection);
1078       fh1Trials->Fill("#sum{ntrials}",trials);
1079       fh1PtHard->Fill(ptHard,eventW);
1080       fh1PtHardNoW->Fill(ptHard,1);
1081       fh1PtHardTrials->Fill(ptHard,trials);
1082    }
1083
1084
1085    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
1086    if(!fESD){
1087       if(fDebug>1) AliError("ESD not available");
1088       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1089    } 
1090  
1091    fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1092    AliAODEvent* aod = NULL;
1093    // take all other information from the aod we take the tracks from
1094    if(!aod && !fIsKine){
1095       if(!fESD)         aod = fAODIn; //not ESD and not kine  => input AOD
1096       else              aod = fAODOut;// ESD or kine
1097    }
1098
1099
1100    if(fNonStdFile.Length()!=0){
1101       // case that we have an AOD extension we can fetch the jets from the extended output
1102       AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1103       fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
1104       if(!fAODExtension){
1105          if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
1106       } 
1107    }
1108     
1109    // ----------------- EVENT SELECTION  --------------------------
1110    fHistEvtSelection->Fill(1); // number of events before event selection
1111
1112    if(!fIsKine){ 
1113       // physics selection
1114       AliInputEventHandler* inputHandler = (AliInputEventHandler*)
1115            ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1116
1117       if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
1118          if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
1119          fHistEvtSelection->Fill(2);
1120          PostData(1, fOutputList);
1121          return;
1122       }
1123
1124    //check AOD pointer
1125       if(!aod){
1126          if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1127          fHistEvtSelection->Fill(3);
1128          PostData(1, fOutputList);
1129          return;
1130       }
1131
1132       // vertex selection for reconstructed data
1133       AliAODVertex* primVtx = aod->GetPrimaryVertex();
1134
1135       if(!primVtx){
1136          if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
1137          fHistEvtSelection->Fill(3);
1138          PostData(1, fOutputList);
1139          return;
1140       }  
1141
1142       Int_t nTracksPrim = primVtx->GetNContributors();
1143       Float_t vtxz = primVtx->GetZ();
1144       //Input events
1145       fhContribVtx->Fill(nTracksPrim);
1146       if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
1147
1148       if((nTracksPrim < fMinContribVtx) ||
1149          (vtxz < fVtxZMin) ||
1150          (vtxz > fVtxZMax)){
1151          if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
1152                          (char*)__FILE__,__LINE__,vtxz);
1153          fHistEvtSelection->Fill(3);
1154          PostData(1, fOutputList);
1155          return;
1156       }else{
1157       //Accepted events
1158          fhContribVtxAccept->Fill(nTracksPrim);
1159          fhVertexZAccept->Fill(vtxz);
1160       }
1161    }else{ //KINE cut on vertex
1162       const AliVVertex *vtxMC = fMcEvent->GetPrimaryVertex();
1163       Float_t zVtx = vtxMC->GetZ();
1164       fhVertexZ->Fill(zVtx);
1165       if(zVtx < fVtxZMin  || zVtx > fVtxZMax){ //vertex cut
1166          fHistEvtSelection->Fill(3);
1167          PostData(1, fOutputList);
1168          return;
1169       }else{
1170          fhVertexZAccept->Fill(zVtx);
1171       }
1172    }
1173    //FK// No event class selection imposed
1174    // event class selection (from jet helper task)
1175    //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
1176    //if(fDebug) Printf("Event class %d", eventClass);
1177    //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
1178    //   fHistEvtSelection->Fill(4);
1179    //   PostData(1, fOutputList);
1180    //   return;
1181    //}
1182
1183    //------------------ CENTRALITY SELECTION ---------------
1184    AliCentrality *cent = 0x0;
1185    Double_t centValue  = 0.0; 
1186    if(fSystem){  //fSystem=0 for pp,   fSystem=1 for pPb
1187       if(fESD){
1188          cent = fESD->GetCentrality();
1189          if(cent) centValue = cent->GetCentralityPercentile("V0M");
1190       }else{
1191          //centValue = aod->GetHeader()->GetCentrality();
1192          centValue = ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
1193       }   
1194       if(fDebug) printf("centrality: %f\n", centValue);
1195       //Input events
1196       fhCentrality->Fill(centValue); 
1197
1198       if(centValue < fCentMin || centValue > fCentMax){
1199          fHistEvtSelection->Fill(4);
1200          PostData(1, fOutputList);
1201          return;
1202       }else{
1203          //Accepted events
1204          fhCentralityAccept->Fill( centValue );
1205       }
1206    }
1207  
1208    //-----------------select disjunct event subsamples ----------------
1209    if(!fIsKine){ //reconstructed data
1210       //Int_t eventnum  = aod->GetHeader()->GetEventNumberESDFile();
1211       AliAODHeader * header = dynamic_cast<AliAODHeader*>(aod->GetHeader());
1212       if(!header) AliFatal("Not a standard AOD");
1213
1214       Int_t eventnum  = header->GetEventNumberESDFile();
1215       Int_t lastdigit = eventnum % 10;
1216       if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
1217          fHistEvtSelection->Fill(5);
1218          PostData(1, fOutputList);
1219          return;
1220        } 
1221    }  
1222
1223    if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
1224    fHistEvtSelection->Fill(0); // accepted events 
1225    // ==================== end event selection ============================ 
1226  
1227    Double_t tmpArrayFive[5];
1228    Double_t tmpArrayFour[4];
1229
1230
1231    TList particleList; //list of tracks
1232    Int_t indexTrigg = -1; 
1233    Double_t rhoFromCellMedian=0.0,    rhoCone=0.0;
1234
1235    if(!fIsKine){ 
1236       //=============== Reconstructed antikt jets =============== 
1237       ReadTClonesArray(fJetBranchName.Data()  , fListJets); 
1238       ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg); 
1239
1240       //============ Estimate background in reconstructed events ===========
1241
1242       //Find Hadron trigger and Estimate rho from cone
1243       indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1244       EstimateBgCone(fListJets, &particleList, rhoCone);
1245
1246       //Estimate rho from cell median minus jets
1247       EstimateBgRhoMedian(fListJetsBg, &particleList, rhoFromCellMedian,0); //real data
1248    }
1249    //==============  analyze generator level MC  ================ 
1250    TList particleListGen; //list of tracks in MC
1251
1252    if(fIsChargedMC || fIsKine){
1253
1254       if(fIsKine){  //pure kine
1255
1256          //================= generated charged antikt jets ================
1257          ReadTClonesArray(fJetBranchNameKine.Data(),   fListJetsGen); 
1258          ReadTClonesArray(fJetBranchNameBgKine.Data(), fListJetsBgGen); 
1259       }else{  
1260          //================= generated charged antikt jets ================
1261          ReadTClonesArray(fJetBranchNameChargMC.Data(),   fListJetsGen); 
1262          ReadTClonesArray(fJetBranchNameBgChargMC.Data(), fListJetsBgGen); 
1263
1264          if(fIsFullMC){ //generated full jets
1265             ReadTClonesArray(fJetBranchNameFullMC.Data(),  fListJetsGenFull); 
1266          }
1267       }
1268       //========================================================
1269       //serarch for charged trigger at the MC generator level
1270
1271       Int_t    indexTriggGen = -1;
1272       Double_t ptTriggGen    = -1;
1273       Int_t    iCounterGen   =  0; //number of entries in particleListGen array
1274       Int_t    triggersMC[200];//list of trigger candidates
1275       Int_t    ntriggersMC   = 0; //index in triggers array
1276       Int_t    triggersMCa[200];   //list of trigger candidates  10%eloss
1277       Int_t    ntriggersMCa   = 0; //index in triggers array   10%eloss
1278       Int_t    triggersMCb[200];   //list of trigger candidates  5%eloss
1279       Int_t    ntriggersMCb   = 0; //index in triggers array     5%eloss
1280
1281
1282       if(!fIsKine){  
1283          if(fESD){//ESD input
1284
1285             AliMCEvent* mcEvent = MCEvent();
1286             if(!mcEvent){
1287                PostData(1, fOutputList);
1288                return;
1289             }
1290
1291             AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
1292             if(pythiaGenHeader){
1293                nTrials = pythiaGenHeader->Trials();
1294                ptHard  = pythiaGenHeader->GetPtHard();
1295                fh1PtHard->Fill(ptHard,eventW);
1296                fh1PtHardNoW->Fill(ptHard,1);
1297                fh1PtHardTrials->Fill(ptHard,nTrials);
1298             }
1299          
1300             for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
1301                if(!mcEvent->IsPhysicalPrimary(it)) continue;
1302                AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
1303                if(!part) continue;  
1304                if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){ 
1305  
1306                   if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1307                      if(indexTriggGen > -1){//trigger candidate was found
1308                         triggersMC[ntriggersMC] = indexTriggGen;
1309                         ntriggersMC++; 
1310                      }
1311                   }
1312  
1313                   iCounterGen++;//index in particleListGen array 
1314                }
1315             }
1316          }else{  //AOD input
1317             if(!fAODIn){
1318                PostData(1, fOutputList);
1319                return;
1320             }
1321             TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
1322             if(!tca){ 
1323                PostData(1, fOutputList);
1324                return;
1325             }
1326             for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
1327                AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1328                if(!part) continue;  
1329                if(!part->IsPhysicalPrimary()) continue;
1330                if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1331  
1332                   if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1333                      if(indexTriggGen > -1){ //trigger candidater was found
1334                         triggersMC[ntriggersMC] = indexTriggGen;
1335                         ntriggersMC++; 
1336                      }
1337                   }
1338  
1339                   iCounterGen++;//number of entries in particleListGen array
1340                }
1341             }
1342          }
1343       }else{ //analyze kine
1344
1345          for(Int_t it = 0; it < fMcEvent->GetNumberOfTracks(); it++){
1346             if(!fMcEvent->IsPhysicalPrimary(it)) continue;
1347             AliMCParticle* part = (AliMCParticle*) fMcEvent->GetTrack(it);
1348             if(!part) continue;  
1349             if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){ 
1350  
1351                if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1352                   if(indexTriggGen > -1){//trigger candidate was found
1353                      triggersMC[ntriggersMC] = indexTriggGen;
1354                      ntriggersMC++; 
1355                   }
1356                }
1357  
1358                iCounterGen++;//index in particleListGen array 
1359             }
1360          }
1361       }
1362  
1363       if(fHardest==0){
1364          Int_t npar = particleListGen.GetEntries();
1365          for(Int_t ip=0; ip < npar; ip++){
1366             AliVParticle *part = (AliVParticle*) particleListGen.At(ip);
1367             if(!part) continue;
1368             
1369             Double_t pta = 0.9 * part->Pt(); //10% energy loss 
1370             Double_t ptb = 0.95 * part->Pt(); //5% energy loss  
1371             if(fTriggerPtRangeLow <= pta && pta < fTriggerPtRangeHigh  &&  ntriggersMCa<200){
1372                triggersMCa[ntriggersMCa] = ip;
1373                ntriggersMCa++;
1374             }
1375
1376             if(fTriggerPtRangeLow <= ptb && ptb < fTriggerPtRangeHigh  &&  ntriggersMCb<200){
1377                triggersMCb[ntriggersMCb] = ip;
1378                ntriggersMCb++;
1379             }
1380          }
1381
1382          if(ntriggersMCa>0){
1383             Int_t rnda     = fRandom->Integer(ntriggersMCa); //0 to ntriggers-1
1384             Int_t indexTriggGena = triggersMCa[rnda];
1385
1386             Double_t deltaPhia, deltaEtaa, deltaRa;
1387             Int_t aaLow = 0; 
1388             Int_t aaHigh = 0; 
1389
1390             //Correlation with single inclusive  TRIGGER
1391             AliVParticle* tGenTa = (AliVParticle*) particleListGen.At(indexTriggGena);  
1392             if(tGenTa){
1393                //for(Int_t ia=0; ia<ntriggersMCa; ia++)
1394                for(Int_t ia=0; ia< npar; ia++){
1395                   //if(indexTriggGena == triggersMCa[ia]) continue;
1396                   if(indexTriggGena == ia) continue;
1397                
1398                   //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCa[ia]);  
1399                   AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ia);  
1400                   if(!tGenTz) continue;
1401                   if(tGenTz->Pt()*0.9<15.0) continue;
1402                   if(tGenTz->Pt()*0.9>50.0) continue;
1403                
1404                   deltaPhia = RelativePhi(tGenTa->Phi(),tGenTz->Phi());
1405                   deltaEtaa = tGenTa->Eta()-tGenTz->Eta(); 
1406                   deltaRa = sqrt(deltaPhia*deltaPhia + deltaEtaa*deltaEtaa);
1407                   
1408                   if(deltaRa<0.4){ 
1409                      if(tGenTz->Pt()*0.9<20.0) aaLow++; //15-20
1410                      else aaHigh++;  //20-50
1411                   }
1412                }
1413             }
1414             fhNofMultipleTriggersConeGenA->Fill(aaLow+aaHigh); // 15-50
1415             fhNofMultipleTriggersConeGenALow->Fill(aaLow); //15-20
1416             fhNofMultipleTriggersConeGenAHigh->Fill(aaHigh);//20-50
1417          }
1418
1419          if(ntriggersMCb>0){
1420             Int_t rndb     = fRandom->Integer(ntriggersMCb); //0 to ntriggers-1
1421             Int_t indexTriggGenb = triggersMCb[rndb];
1422
1423             Double_t deltaPhib, deltaEtab, deltaRb;
1424             Int_t bbLow = 0; 
1425             Int_t bbHigh = 0; 
1426
1427             //Correlation with single inclusive  TRIGGER
1428             AliVParticle* tGenTb = (AliVParticle*) particleListGen.At(indexTriggGenb);  
1429             if(tGenTb){
1430                //for(Int_t ib=0; ib<ntriggersMCb; ib++)
1431                for(Int_t ib=0; ib<npar; ib++){
1432                   //if(indexTriggGenb == triggersMCb[ib]) continue;
1433                   if(indexTriggGenb == ib) continue;
1434                
1435                   //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCb[ib]);  
1436                   AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ib);  
1437                   if(!tGenTz) continue;
1438                   if(tGenTz->Pt()*0.95<15.0) continue;
1439                   if(tGenTz->Pt()*0.95>50.0) continue;
1440                
1441                   deltaPhib = RelativePhi(tGenTb->Phi(),tGenTz->Phi());
1442                   deltaEtab = tGenTb->Eta()-tGenTz->Eta(); 
1443                   deltaRb = sqrt(deltaPhib*deltaPhib + deltaEtab*deltaEtab);
1444                   
1445                   if(deltaRb<0.4){
1446                      if(tGenTz->Pt()*0.95<20.0) bbLow++; //15-20
1447                      else bbHigh++;
1448                   }
1449                }
1450             }
1451             fhNofMultipleTriggersConeGenB->Fill(bbLow+bbHigh); //15-50
1452             fhNofMultipleTriggersConeGenBLow->Fill(bbLow);//15-20
1453             fhNofMultipleTriggersConeGenBHigh->Fill(bbHigh);//20-50
1454          }
1455       }
1456
1457
1458  
1459       //==============  Estimate bg in generated events ==============
1460       Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
1461
1462       //Estimate rho from cone
1463       EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
1464
1465       //Estimate rho from cell median minus jets
1466       EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
1467
1468       //============  Generator trigger+jet ==================
1469       if(fHardest==0){ //single inclusive trigger
1470          if(ntriggersMC>0){ //there is at least one trigger 
1471             Int_t rnd     = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
1472             indexTriggGen = triggersMC[rnd];
1473
1474             fhNofMultipleTriggersGen->Fill(ntriggersMC-1);
1475
1476             Double_t deltaPhi, deltaEta, deltaR;
1477             Int_t iLow = 0; 
1478             Int_t iHigh = 0; 
1479
1480             //Correlation with single inclusive  TRIGGER
1481             AliVParticle* tGenT1 = (AliVParticle*) particleListGen.At(indexTriggGen);  
1482             if(tGenT1){
1483                //for(Int_t ia=0; ia<ntriggersMC; ia++)
1484                for(Int_t ia=0; ia< particleListGen.GetEntries(); ia++){
1485                   //if(indexTriggGen == triggersMC[ia]) continue;
1486                   if(indexTriggGen == ia) continue;
1487                
1488                   //AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(triggersMC[ia]);  
1489                   AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(ia);  
1490                   if(!tGenT2) continue;
1491                   if(tGenT2->Pt()<15.0) continue; 
1492                   if(tGenT2->Pt()>50.0) continue; 
1493                   deltaPhi = RelativePhi(tGenT1->Phi(),tGenT2->Phi());
1494                   deltaEta = tGenT1->Eta()-tGenT2->Eta(); 
1495                   deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1496                   
1497                   fhDeltaRMultTriggersGen->Fill(deltaR);
1498                   fhDeltaPhiMultTriggersGen->Fill(deltaPhi);
1499
1500                   if(deltaR<0.4){
1501                      if(tGenT2->Pt()<20.0) iLow++;
1502                      else iHigh++;
1503                   }
1504                }
1505             }
1506             fhNofMultipleTriggersConeGen->Fill(iLow+iHigh);
1507             fhNofMultipleTriggersConeGenLow->Fill(iLow);
1508             fhNofMultipleTriggersConeGenHigh->Fill(iHigh);
1509
1510          }else{
1511             indexTriggGen = -1; //trigger not found
1512          }
1513       }
1514
1515       //-----------
1516       Int_t ilowGen  = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1517       Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
1518       Bool_t fillOnceGen = kTRUE;
1519       //-----------
1520       
1521       for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1522          indexTriggGen = igen; //trigger hadron 
1523
1524          if(indexTriggGen == -1) continue;  
1525          AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);  
1526          if(!triggerGen) continue;
1527  
1528          if(fHardest >= 2){ 
1529             if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6  
1530          }
1531          if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1532
1533          ptTriggGen = triggerGen->Pt(); //count triggers
1534          fh2NtriggersGen->Fill(centValue, ptTriggGen);
1535
1536          //Count jets and trigger-jet pairs at MC  generator level
1537          for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1538             AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1539             if(!jet) continue;
1540             Double_t etaJetGen = jet->Eta();
1541             Double_t areaJetGen = jet->EffectiveAreaCharged();
1542  
1543             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){ 
1544
1545                Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1546                Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1547
1548                //Rongrong's analysis
1549                Double_t dPhiGen = jet->Phi() - triggerGen->Phi();
1550                if(dPhiGen>2*TMath::Pi())    dPhiGen -= 2*TMath::Pi();
1551                if(dPhiGen<-2*TMath::Pi())   dPhiGen += 2*TMath::Pi();
1552                if(dPhiGen<-0.5*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1553                if(dPhiGen>1.5*TMath::Pi())  dPhiGen -= 2*TMath::Pi();
1554
1555                //[A,pTjet,pTtrig,dphi]
1556                tmpArrayFour[0] = areaJetGen;
1557                tmpArrayFour[1] = jet->Pt();
1558                tmpArrayFour[2] = ptTriggGen;
1559                tmpArrayFour[3] = dPhiGen;
1560
1561                fHJetPhiCorrGen->Fill(tmpArrayFour); // Dphi distribution jet-triger
1562
1563                //[A,pTjet-pTUe,pTtrig,dphi]
1564                tmpArrayFour[1] = jet->Pt() - ptUeFromCellMedianGen;
1565                fHJetPhiCorrSubUeMedianGen->Fill(tmpArrayFour); 
1566
1567                //[A,pTjet-pTUe,pTtrig,dphi]
1568                tmpArrayFour[1] = jet->Pt() - ptUeConeGen;
1569                fHJetPhiCorrSubUeConeGen->Fill(tmpArrayFour);
1570
1571                //Leticia's analysis
1572                Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi()); 
1573                if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1574
1575                //Centrality, A, pT, pTtrigg
1576                tmpArrayFive[0] = centValue;
1577                tmpArrayFive[1] = areaJetGen;
1578                tmpArrayFive[2] = jet->Pt();
1579                tmpArrayFive[3] = ptTriggGen;
1580                tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1581                fHJetSpecGen->Fill(tmpArrayFive);
1582
1583                //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1584                tmpArrayFive[0] = centValue;
1585                tmpArrayFive[1] = areaJetGen;
1586                tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1587                tmpArrayFive[3] = ptTriggGen;
1588                tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1589                fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1590
1591                //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1592                tmpArrayFive[0] = centValue;
1593                tmpArrayFive[1] = areaJetGen;
1594                tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1595                tmpArrayFive[3] = ptTriggGen;
1596                tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1597                fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
1598
1599                //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1600                tmpArrayFive[0] = areaJetGen;
1601                tmpArrayFive[1] = jet->Pt();
1602                tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1603                tmpArrayFive[3] = ptUeFromCellMedianGen;
1604                tmpArrayFive[4] = rhoFromCellMedianGen;
1605                fHJetUeMedianGen->Fill(tmpArrayFive); 
1606
1607                //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
1608                tmpArrayFive[0] = areaJetGen;
1609                tmpArrayFive[1] = jet->Pt();
1610                tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1611                tmpArrayFive[3] = ptUeConeGen;
1612                tmpArrayFive[4] = rhoConeGen;
1613                fHJetUeConeGen->Fill(tmpArrayFive);
1614
1615                if(fillOnceGen){
1616                   Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1617                   fHRhoUeMedianVsConeGen->Fill(fillRhoGen); 
1618                   fillOnceGen = kFALSE;
1619                }
1620             }//back to back jet-trigger pair
1621          }//jet loop
1622       }//trigger loop
1623
1624       if(fFillRespMx && !fIsKine){
1625          //================ RESPONSE MATRIX ===============
1626          //Count jets and trigger-jet pairs at MC  generator level
1627          for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1628             AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1629             if(!jet) continue;
1630             Double_t etaJetGen = jet->Eta();
1631             Double_t ptJetGen  = jet->Pt();
1632        
1633             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){ 
1634                fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
1635
1636                Double_t areaJetGen = jet->EffectiveAreaCharged();
1637                Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1638                Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1639
1640                fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen); 
1641                fhJetPtSubUeConeGen->Fill(ptJetGen   - ptUeConeGen);    
1642             }
1643          }
1644          if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1645          
1646             Int_t ng = (Int_t) fListJetsGen->GetEntries();
1647             Int_t nr = (Int_t) fListJets->GetEntries();
1648
1649             //Find closest MC generator - reconstructed jets
1650             if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
1651             if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
1652
1653             if(fDebug){
1654                Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
1655                Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
1656             }
1657             //matching of MC genrator level and reconstructed jets
1658             AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug); 
1659
1660             // Fill response matrix
1661             for(Int_t ir = 0; ir < nr; ir++){
1662                AliAODJet *recJet   = (AliAODJet*) fListJets->At(ir);
1663                Double_t etaJetRec  = recJet->Eta();
1664                Double_t ptJetRec   = recJet->Pt();
1665                Double_t areaJetRec = recJet->EffectiveAreaCharged();
1666                //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1667
1668                if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){ 
1669                   Int_t ig = faGenIndex[ir]; //associated generator level jet
1670                   if(ig >= 0 && ig < ng){
1671                      if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1672                      AliAODJet *genJet  = (AliAODJet*) fListJetsGen->At(ig);
1673                      Double_t ptJetGen  = genJet->Pt();
1674                      Double_t etaJetGen = genJet->Eta();
1675
1676                      //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1677                      if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1678                         fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
1679
1680                         Double_t areaJetGen  = genJet->EffectiveAreaCharged();
1681                         Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1682                         Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1683                         Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
1684                         Double_t ptUeConeRec         = rhoCone*areaJetRec;
1685                         fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec, 
1686                                                            ptJetGen-ptUeFromCellMedianGen);
1687                         fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
1688                      }
1689                   }//ig>=0
1690                }//rec jet in eta acceptance
1691             }//loop over reconstructed jets
1692          }// # of  rec jets >0
1693       
1694          //=========================== Full jet vs charged jet matrix  ==========================
1695          if(fIsFullMC){
1696             //Count full jets and charged-jet pairs at MC  generator level
1697             for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1698                AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1699                if(!jetFull) continue;
1700                Double_t etaJetFull = jetFull->Eta();
1701                Double_t ptJetFull  = jetFull->Pt();
1702
1703                if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1704                   fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets 
1705                }
1706             }
1707             if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1708                Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1709                Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1710  
1711                //Find closest MC generator full - charged jet
1712                if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1713                if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1714  
1715                if(fDebug){
1716                   Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
1717                   Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
1718                }
1719                //matching of MC genrator level and reconstructed jets
1720                AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1721
1722               // Fill response matrix
1723                for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1724                   AliAODJet *chJet  = (AliAODJet*) fListJetsGen->At(ichr);
1725                   Double_t etaJetCh = chJet->Eta();
1726                   Double_t ptJetCh  = chJet->Pt();
1727                   //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1728     
1729                   if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1730                      Int_t iful = faGenIndex[ichr]; //associated generator level jet
1731                      if(iful >= 0 && iful < nful){
1732                         if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1733                         AliAODJet *genJetFull  = (AliAODJet*) fListJetsGenFull->At(iful);
1734                         Double_t ptJetFull  = genJetFull->Pt();
1735                         Double_t etaJetFull = genJetFull->Eta();
1736
1737                         //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1738                         if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1739                            fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1740                         }
1741                      }//iful>=0
1742                   }//rec jet in eta acceptance
1743                }//loop over reconstructed jets
1744             }// # of  rec jets >0
1745          }//pointer MC generator jets
1746       } //fill resp mx only for bin 
1747    }//analyze generator level MC
1748     
1749     
1750    if(fIsKine){ //skip reconstructed data analysis in case of kine 
1751       PostData(1, fOutputList);
1752       return;
1753    } 
1754    //=============  RECONSTRUCTED INCLUSIVE JETS ===============
1755
1756    Double_t etaJet  = 0.0;
1757    Double_t pTJet   = 0.0;
1758    Double_t areaJet = 0.0;
1759    Double_t phiJet  = 0.0;
1760    //Int_t indexLeadingJet     = -1;
1761    //Double_t pTLeadingJet     = -10.0; 
1762    //Double_t areaLeadingJet   = -10.0;
1763   
1764    for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1765       AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1766       if(!jet) continue;
1767       etaJet  = jet->Eta();
1768       phiJet  = jet->Phi();
1769       pTJet   = jet->Pt();
1770       if(pTJet==0) continue; 
1771      
1772       if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1773       /*areaJet = jet->EffectiveAreaCharged();*/
1774
1775       //Jet Diagnostics---------------------------------
1776       fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
1777       fhJetEta->Fill(pTJet, etaJet);
1778       //search for leading jet
1779       /*if(pTJet > pTLeadingJet){
1780          indexLeadingJet  = ij; 
1781          pTLeadingJet     = pTJet; 
1782          areaLeadingJet   = areaJet; 
1783       } 
1784  
1785       // raw spectra of INCLUSIVE jets  
1786       //Centrality, pTjet, A
1787       Double_t fillraw[] = { centValue,
1788                              pTJet,
1789                              areaJet
1790                            };
1791       fHJetPtRaw->Fill(fillraw);*/
1792    }
1793    /*
1794    if(indexLeadingJet > -1){ 
1795       // raw spectra of LEADING jets  
1796       //Centrality, pTjet,  A
1797       Double_t fillleading[] = { centValue,
1798                                  pTLeadingJet,
1799                                  areaLeadingJet
1800                                };
1801       fHLeadingJetPtRaw->Fill(fillleading);
1802    } 
1803    */
1804  
1805    // ===============  RECONSTRUCTED TRIGGER-JET PAIRS ================
1806    if(fIsChargedMC && fFillRespMx){
1807      FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1808    }
1809    Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
1810    //set ranges of the trigger loop
1811    Int_t ilow  = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1812    Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
1813
1814    for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1815       indexTrigg = itrk; //trigger hadron with pT >10 GeV 
1816  
1817       if(indexTrigg < 0) continue;
1818
1819       AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);     
1820       if(!triggerHadron){  
1821          PostData(1, fOutputList);
1822          return;
1823       }
1824       if(fHardest >= 2){ 
1825          if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6  
1826       }
1827       if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
1828  
1829       //Fill trigger histograms
1830       fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1831       fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1832       fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
1833
1834   
1835       //---------- make trigger-jet pairs ---------
1836       //Int_t injet4     = 0;
1837       //Int_t injet      = 0; 
1838
1839       for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1840          AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1841          if(!jet) continue;
1842          etaJet  = jet->Eta();
1843          phiJet  = jet->Phi();
1844          pTJet   = jet->Pt();
1845          if(pTJet==0) continue; 
1846      
1847          if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1848          areaJet = jet->EffectiveAreaCharged();
1849
1850         //subtract bg using estinates base on median of kT jets
1851          Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
1852          Double_t ptUeCone         = rhoCone*areaJet;
1853
1854          //if(areaJet >= 0.07) injet++; 
1855          //if(areaJet >= 0.4)  injet4++;
1856
1857          Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet); 
1858          fhDphiTriggerJet->Fill(dphi); //Input
1859
1860          //Dphi versus jet pT   
1861          //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg 
1862          /*Double_t filldp[] = { centValue,
1863                                dphi,
1864                                pTJet,
1865                                triggerHadron->Pt()
1866                              };
1867          fHDphiVsJetPtAll->Fill(filldp);
1868          */ 
1869          //Rongrong's analysis
1870
1871          Double_t dPhi = phiJet - triggerHadron->Phi();
1872          if(dPhi>2*TMath::Pi())    dPhi -= 2*TMath::Pi();
1873          if(dPhi<-2*TMath::Pi())   dPhi += 2*TMath::Pi();
1874          if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
1875          if(dPhi>1.5*TMath::Pi())  dPhi -= 2*TMath::Pi();
1876
1877          //[A,pTjet,pTtrig,dphi]
1878          tmpArrayFour[0] = areaJet;
1879          tmpArrayFour[1] = pTJet;
1880          tmpArrayFour[2] = triggerHadron->Pt();
1881          tmpArrayFour[3] = dPhi;
1882
1883          fHJetPhiCorr->Fill(tmpArrayFour); // Dphi distribution jet-triger
1884
1885          //[A,pTjet-pTUe,pTtrig,dphi]
1886          tmpArrayFour[1] = pTJet - ptUeFromCellMedian;
1887          fHJetPhiCorrSubUeMedian->Fill(tmpArrayFour); 
1888
1889          //[A,pTjet-pTUe,pTtrig,dphi]
1890          tmpArrayFour[1] = pTJet - ptUeCone;
1891          fHJetPhiCorrSubUeCone->Fill(tmpArrayFour);
1892
1893          //Leticia's analysis
1894
1895          // Select back to back trigger - jet pairs
1896          if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1897          fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1898           
1899          //NO bg subtraction
1900          //Centrality, A, pTjet, pTtrigg, dphi
1901          tmpArrayFive[0] = centValue;
1902          tmpArrayFive[1] = areaJet;
1903          tmpArrayFive[2] = pTJet;
1904          tmpArrayFive[3] = triggerHadron->Pt();
1905          tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1906          fHJetSpec->Fill(tmpArrayFive);
1907
1908           //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1909          tmpArrayFive[0] = centValue;
1910          tmpArrayFive[1] = areaJet;
1911          tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
1912          tmpArrayFive[3] = triggerHadron->Pt();
1913          tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1914          fHJetSpecSubUeMedian->Fill(tmpArrayFive);
1915
1916          //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1917          tmpArrayFive[0] = centValue;
1918          tmpArrayFive[1] = areaJet;
1919          tmpArrayFive[2] = pTJet - ptUeCone;
1920          tmpArrayFive[3] = triggerHadron->Pt();
1921          tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1922          fHJetSpecSubUeCone->Fill(tmpArrayFive);
1923
1924          //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1925          tmpArrayFive[0] = areaJet;
1926          tmpArrayFive[1] = pTJet;
1927          tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
1928          tmpArrayFive[3] = ptUeFromCellMedian;
1929          tmpArrayFive[4] = rhoFromCellMedian;
1930          fHJetUeMedian->Fill(tmpArrayFive);
1931  
1932          //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
1933          tmpArrayFive[0] = areaJet;
1934          tmpArrayFive[1] = pTJet;
1935          tmpArrayFive[2] = pTJet - ptUeCone;
1936          tmpArrayFive[3] = ptUeCone;
1937          tmpArrayFive[4] = rhoCone;
1938          fHJetUeCone->Fill(tmpArrayFive);
1939
1940          if(filledOnce){ //fill for each event only once
1941             Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
1942             fHRhoUeMedianVsCone->Fill(fillRho);
1943             filledOnce = kFALSE;
1944          }
1945       }//jet loop
1946  
1947       //Fill Jet Density In the Event A>0.07
1948       /*if(injet>0){
1949          Double_t filldens[]={ (Double_t) particleList.GetEntries(),
1950                             injet/fkAcceptance,
1951                             triggerHadron->Pt()
1952                           };
1953          fHJetDensity->Fill(filldens);
1954       } 
1955
1956       //Fill Jet Density In the Event A>0.4
1957       if(injet4>0){ 
1958          Double_t filldens4[]={ (Double_t) particleList.GetEntries(), 
1959                                 injet4/fkAcceptance,
1960                                 triggerHadron->Pt()
1961                               };
1962          fHJetDensityA4->Fill(filldens4);
1963       }*/
1964    }
1965
1966
1967    PostData(1, fOutputList);
1968 }
1969
1970 //----------------------------------------------------------------------------
1971 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1972 {
1973    // Draw result to the screen
1974    // Called once at the end of the query
1975
1976    if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1977    if(!GetOutputData(1)) return;
1978 }
1979
1980
1981 //----------------------------------------------------------------------------
1982 Int_t  AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1983    //Fill the list of accepted tracks (passed track cut)
1984    //return consecutive index of the hardest ch hadron in the list
1985    Int_t iCount        = 0;
1986    AliAODEvent *aodevt = NULL;
1987
1988    if(!fESD) aodevt = fAODIn;
1989    else      aodevt = fAODOut;   
1990
1991    if(!aodevt) return -1;
1992
1993    Int_t    index = -1; //index of the highest particle in the list
1994    Double_t ptmax = -10;
1995    Int_t    triggers[200];
1996    Int_t    ntriggers = 0; //index in triggers array
1997  
1998
1999    for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
2000       //AliAODTrack *tr = aodevt->GetTrack(it);
2001       AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodevt->GetTrack(it));
2002       if(!tr) AliFatal("Not a standard AOD");
2003
2004       if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
2005       //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
2006       if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
2007       if(tr->Pt() < fTrackLowPtCut) continue;
2008       list->Add(tr);
2009
2010       //Search trigger:
2011       if(fHardest>0){ //leading particle 
2012          if(tr->Pt()>ptmax){ 
2013             ptmax = tr->Pt();   
2014             index = iCount;
2015          }
2016       }
2017     
2018       if(fHardest==0 && ntriggers<200){ //single inclusive 
2019          if(fTriggerPtRangeLow <= tr->Pt() && 
2020             tr->Pt() < fTriggerPtRangeHigh){ 
2021             triggers[ntriggers] = iCount;
2022             ntriggers++;
2023          }
2024       }
2025
2026       iCount++;
2027    }
2028
2029    if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
2030       Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
2031       index     = triggers[rnd];
2032
2033       fhNofMultipleTriggers->Fill(ntriggers-1);
2034      
2035       Double_t deltaPhi, deltaEta, deltaR;
2036       Int_t iLow=0;
2037       Int_t iHigh=0;
2038       //Correlation with single inclusive trigger
2039       AliVParticle* tGent1 = (AliVParticle*) list->At(index);  
2040       if(tGent1){
2041          //for(Int_t ia=0; ia<ntriggers; ia++)
2042          for(Int_t ia=0; ia<list->GetEntries(); ia++){
2043             //if(triggers[ia]==index) continue;
2044             if(ia==index) continue;
2045             //AliVParticle* tGent2 = (AliVParticle*) list->At(triggers[ia]);  
2046             AliVParticle* tGent2 = (AliVParticle*) list->At(ia);  
2047             if(!tGent2) continue;
2048             if(tGent2->Pt()<15.0) continue;
2049             if(tGent2->Pt()>50.0) continue;
2050             deltaPhi = RelativePhi(tGent1->Phi(),tGent2->Phi());
2051             deltaEta = tGent1->Eta()-tGent2->Eta(); 
2052             deltaR   = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
2053
2054             fhDeltaRMultTriggers->Fill(deltaR);
2055             fhDeltaPhiMultTriggers->Fill(deltaPhi);
2056
2057             if(deltaR<0.4){
2058                if(tGent2->Pt()<20.0) iLow++;
2059                else iHigh++;
2060             }
2061          } 
2062       }
2063       fhNofMultipleTriggersCone->Fill(iLow+iHigh);
2064       fhNofMultipleTriggersConeLow->Fill(iLow);
2065       fhNofMultipleTriggersConeHigh->Fill(iHigh);
2066
2067    }
2068
2069    return index;
2070 }
2071
2072 //----------------------------------------------------------------------------
2073 /*
2074 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
2075    //calculate sum of track pT in the cone perpendicular in phi to the jet 
2076    //jetR = cone radius
2077    // jetPhi, jetEta = direction of the jet 
2078    Int_t numberOfTrks = trkList->GetEntries();
2079    Double_t pTsum = 0.0;
2080    Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
2081    for(Int_t it=0; it<numberOfTrks; it++){
2082       AliVParticle *trk = (AliVParticle*) trkList->At(it); 
2083       Double_t dphi = RelativePhi(perpConePhi,trk->Phi());     
2084       Double_t deta = trk->Eta()-jetEta;     
2085
2086       if( (dphi*dphi + deta*deta)< (jetR*jetR)){
2087          pTsum += trk->Pt();
2088       } 
2089    }
2090
2091    return pTsum;
2092 }
2093 */
2094 //----------------------------------------------------------------------------
2095
2096 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
2097    //Get relative azimuthal angle of two particles -pi to pi
2098    if      (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
2099    else if (vphi > TMath::Pi())  vphi -= TMath::TwoPi();
2100
2101    if      (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
2102    else if (mphi > TMath::Pi())  mphi -= TMath::TwoPi();
2103
2104    Double_t dphi = mphi - vphi;
2105    if      (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
2106    else if (dphi > TMath::Pi())  dphi -= TMath::TwoPi();
2107
2108    return dphi;//dphi in [-Pi, Pi]
2109 }
2110
2111
2112 //----------------------------------------------------------------------------
2113 Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
2114    //fill track efficiency denominator
2115    if(trk->Pt() < fTrackLowPtCut) return kFALSE;
2116    if(trk->Charge()==0)        return kFALSE;
2117    if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
2118    trkList->Add(trk);
2119    if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
2120
2121    //Search trigger:
2122    if(fHardest>0){ //leading particle 
2123       if(ptLeading < trk->Pt()){
2124          index      = counter;
2125          ptLeading  = trk->Pt();
2126       }
2127    }
2128
2129    if(fHardest==0){ //single inclusive 
2130       index = -1;  
2131       if(fTriggerPtRangeLow <= trk->Pt() && 
2132             trk->Pt() < fTriggerPtRangeHigh){
2133          index      = counter;
2134       }
2135    }
2136
2137    return kTRUE;
2138
2139
2140 //----------------------------------------------------------------------------
2141 void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
2142    //fill track efficiency numerator 
2143    if(!recList) return;
2144    if(!genList) return;
2145    Int_t nRec = recList->GetEntries();
2146    Int_t nGen = genList->GetEntries();
2147    for(Int_t ir=0; ir<nRec; ir++){ 
2148       AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
2149       if(!trkRec) continue; 
2150       Int_t recLabel = TMath::Abs(trkRec->GetLabel());
2151       Bool_t matched = kFALSE;
2152  
2153       for(Int_t ig=0; ig<nGen; ig++){
2154         AliVParticle *trkGen = (AliVParticle*) genList->At(ig);  
2155         if(!trkGen) continue;
2156         Int_t genLabel = TMath::Abs(trkGen->GetLabel()); 
2157         if(recLabel==genLabel){
2158           fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
2159           matched = kTRUE;
2160           break;
2161         }
2162       }
2163       if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta()); 
2164    }
2165
2166    return;
2167 }
2168 //________________________________________________________________
2169 void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart,  Double_t &rhoMedian, Int_t mode){
2170    //Estimate background rho by means of integrating track pT outside identified jet cones
2171
2172    rhoMedian = 0;
2173
2174    //identify leading jet in the full track acceptance
2175    Int_t idxLeadingJet    = -1;
2176    Double_t pTleading     = 0.0; 
2177    AliAODJet* jet = NULL;
2178  
2179    for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
2180       jet = (AliAODJet*)(listJet->At(ij));
2181       if(!jet) continue;
2182       if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
2183       if(jet->Pt() > pTleading){
2184          idxLeadingJet = ij; 
2185          pTleading     = jet->Pt();
2186       }
2187    } 
2188
2189    Int_t  njets = 0;
2190    static Double_t jphi[100];
2191    static Double_t jeta[100];
2192    static Double_t jRsquared[100];
2193
2194    for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
2195
2196       jet = (AliAODJet*)(listJet->At(ij));
2197       if(!jet) continue;
2198       if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; 
2199
2200       if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
2201          jphi[njets]      = RelativePhi(jet->Phi(),0.0); //-pi,pi
2202          jeta[njets]      = jet->Eta();
2203          jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
2204          njets++;
2205       }
2206    }
2207
2208    
2209    static Double_t nOutCone[10][4];
2210    static Double_t sumPtOutOfCone[10][4];
2211
2212    for(Int_t ie=0; ie<fnEta; ie++){
2213       for(Int_t ip=0; ip<fnPhi; ip++){
2214          nOutCone[ip][ie]       = 0.0;     //initialize counter
2215          sumPtOutOfCone[ip][ie] = 0.0;
2216       }
2217    }
2218
2219    Double_t rndphi, rndeta;
2220    Double_t rndphishift, rndetashift;
2221    Double_t dphi, deta; 
2222    Bool_t   bIsInCone;
2223
2224    for(Int_t it=0; it<fnTrials; it++){
2225
2226       rndphi = fRandom->Uniform(0, fPhiSize);
2227       rndeta = fRandom->Uniform(0, fEtaSize);
2228                               
2229       for(Int_t ip=0; ip<fnPhi; ip++){  //move radom position to each cell
2230          rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
2231          for(Int_t ie=0; ie<fnEta; ie++){
2232             rndetashift = rndeta + ie*fEtaSize - fEtaSize;
2233
2234             bIsInCone = 0; //tag if trial is in the jet cone
2235             for(Int_t ij=0; ij<njets; ij++){
2236                deta = jeta[ij] - rndetashift;
2237                dphi = RelativePhi(rndphishift,jphi[ij]);
2238                if((dphi*dphi + deta*deta) < jRsquared[ij]){
2239                   bIsInCone = 1;
2240                   break;
2241                }
2242             }
2243             if(!bIsInCone) nOutCone[ip][ie]++;
2244          }
2245       }
2246    }
2247   
2248
2249    //Sum particle pT outside jets in cells
2250    Int_t npart = listPart->GetEntries();
2251    Int_t phicell,etacell; 
2252    AliVParticle *part = NULL; 
2253    for(Int_t ip=0; ip < npart; ip++){ 
2254          
2255       part = (AliVParticle*) listPart->At(ip);
2256       if(!part) continue;
2257       if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
2258  
2259       bIsInCone = 0; //init
2260       for(Int_t ij=0; ij<njets; ij++){
2261          dphi = RelativePhi(jphi[ij], part->Phi());
2262          deta = jeta[ij] - part->Eta();
2263          if((dphi*dphi + deta*deta) < jRsquared[ij]){
2264             bIsInCone = 1;
2265             break;
2266          }
2267       } 
2268       if(!bIsInCone){
2269          phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
2270          etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
2271          sumPtOutOfCone[phicell][etacell]+= part->Pt();
2272       }
2273    }
2274
2275    static Double_t rhoInCells[20];
2276    Double_t  relativeArea;
2277    Int_t  nCells=0;
2278    Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
2279    for(Int_t ip=0; ip<fnPhi; ip++){
2280       for(Int_t ie=0; ie<fnEta; ie++){
2281          relativeArea = nOutCone[ip][ie]/fnTrials;
2282          //cout<<"RA=  "<< relativeArea<<"  BA="<<bufferArea<<"    BPT="<< bufferPt <<endl;
2283
2284          bufferArea += relativeArea;
2285          bufferPt   += sumPtOutOfCone[ip][ie];
2286          if(bufferArea > fJetFreeAreaFrac){ 
2287             rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
2288
2289             if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);  
2290             else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
2291
2292             bufferArea = 0.0; 
2293             bufferPt   = 0.0;
2294             nCells++;
2295          }
2296       }
2297    }
2298
2299
2300    if(nCells>0){
2301      rhoMedian = TMath::Median(nCells, rhoInCells);
2302      if(mode==1){ //mc data 
2303         fhEntriesToMedianGen->Fill(nCells);
2304      }else{ //real data
2305         fhEntriesToMedian->Fill(nCells);
2306      } 
2307    } 
2308
2309 }
2310
2311 //__________________________________________________________________
2312 void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart,  Double_t &rhoPerpCone){
2313
2314    rhoPerpCone = 0.0; //init
2315
2316    if(!listJet) return;
2317    Int_t njet  = listJet->GetEntries();
2318    Int_t npart = listPart->GetEntries();
2319    Double_t pTleading  = 0.0;
2320    Double_t phiLeading = 1000.;
2321    Double_t etaLeading = 1000.;
2322    
2323    for(Int_t jsig=0; jsig < njet; jsig++){ 
2324       AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
2325       if(!signaljet) continue;
2326       if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
2327       if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
2328          pTleading  = signaljet->Pt();
2329          phiLeading = signaljet->Phi();
2330          etaLeading = signaljet->Eta();
2331       }
2332    }   
2333   
2334    Double_t phileftcone  = phiLeading + TMath::Pi()/2; 
2335    Double_t phirightcone = phiLeading - TMath::Pi()/2; 
2336    Double_t dp, de;
2337
2338    for(Int_t ip=0; ip < npart; ip++){ 
2339          
2340       AliVParticle *part = (AliVParticle*) listPart->At(ip);
2341       if(!part){
2342          continue;
2343       }
2344       
2345       dp = RelativePhi(phileftcone, part->Phi());
2346       de = etaLeading - part->Eta();
2347       if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2348
2349       dp = RelativePhi(phirightcone, part->Phi());
2350       if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2351  
2352    }
2353
2354    //normalize total pT by two times cone are 
2355    rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
2356
2357 }
2358 //__________________________________________________________________
2359
2360 void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
2361    //Convert TClones array of jets to TList 
2362
2363    if(!strlen(bname.Data())){
2364       AliError(Form("Jet branch %s not set.", bname.Data()));
2365       return;
2366    }
2367
2368    TClonesArray *array=0x0;
2369    if(fUseExchContainer){ //take input from exchange containers
2370       if(fJetBranchNameKine.Length()>0 && bname.CompareTo(fJetBranchNameKine.Data())==0){
2371          array = dynamic_cast<TClonesArray*>(GetInputData(1)); //connect exchange container slot 1
2372       }
2373       if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
2374         array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
2375       }
2376    }else{ //take input from AOD
2377       if(fAODOut&&!array){
2378          array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
2379       }
2380       if(fAODExtension&&!array){
2381          array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
2382       }
2383       if(fAODIn&&!array){
2384          array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
2385       }
2386    }
2387
2388    list->Clear(); //clear list beforehand
2389
2390    if(array){
2391       if(fDebug){
2392          Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
2393       }
2394
2395       for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
2396          AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
2397          if (jet) list->Add(jet);
2398       }
2399    }
2400
2401    return;
2402 }