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 // *******************************************
11 /**************************************************************************
12 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
14 * Author: The ALICE Off-line Project. *
15 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
28 #include "TClonesArray.h"
36 #include "THnSparse.h"
46 #include "AliAnalysisTask.h"
47 #include "AliAnalysisManager.h"
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//
71 ClassImp(AliAnalysisTaskJetCorePP)
73 //Filip Krizek 1st March 2013
75 //---------------------------------------------------------------------
76 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
85 fJetBranchNameChargMC(""),
86 fJetBranchNameKine(""),
87 fJetBranchNameFullMC(""),
89 fJetBranchNameBgChargMC(""),
90 fJetBranchNameBgKine(""),
93 fListJetsGenFull(0x0),
97 fSystem(0), //pp=0 pPb=1
102 fOfflineTrgMask(AliVEvent::kAny),
113 fTrackLowPtCut(0.15),
114 fUseExchContainer(0),
116 fHistEvtSelection(0x0),
119 fHJetSpecSubUeMedian(0x0),
120 fHJetSpecSubUeCone(0x0),
122 fHJetPhiCorrSubUeMedian(0x0),
123 fHJetPhiCorrSubUeCone(0x0),
126 fHRhoUeMedianVsCone(0x0),
128 //fHJetDensityA4(0x0),
134 fhVertexZAccept(0x0),
136 fhContribVtxAccept(0x0),
137 fhDphiTriggerJet(0x0),
138 fhDphiTriggerJetAccept(0x0),
140 fhCentralityAccept(0x0),
141 fhNofMultipleTriggers(0x0),
142 fhNofMultipleTriggersCone(0x0),
143 fhNofMultipleTriggersConeLow(0x0),
144 fhNofMultipleTriggersConeHigh(0x0),
145 fhDeltaRMultTriggersLow(0x0),
146 fhDeltaRMultTriggersHigh(0x0),
147 fhDeltaPhiMultTriggersLow(0x0),
148 fhDeltaPhiMultTriggersHigh(0x0),
149 fhDeltaPhiMultTriggersInclLow(0x0),
150 fhDeltaPhiMultTriggersInclHigh(0x0),
151 fhInclTrigCounter(0x0),
153 //fHLeadingJetPtRaw(0x0),
154 //fHDphiVsJetPtAll(0x0),
155 fhJetPtGenVsJetPtRec(0x0),
156 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
157 fhJetPtGenVsJetPtRecSubUeCone(0x0),
159 fhJetPtSubUeMedianGen(0x0),
160 fhJetPtSubUeConeGen(0x0),
161 fhJetPtResolutionVsPtGen(0x0),
162 fhJetPtResolutionVsPtConeGen(0x0),
163 fhJetPtGenChargVsJetPtGenFull(0x0),
165 fh2NtriggersGen(0x0),
167 fHJetSpecSubUeMedianGen(0x0),
168 fHJetSpecSubUeConeGen(0x0),
169 fHJetPhiCorrGen(0x0),
170 fHJetPhiCorrSubUeMedianGen(0x0),
171 fHJetPhiCorrSubUeConeGen(0x0),
172 fHJetUeMedianGen(0x0),
174 fhPtTrkTruePrimRec(0x0),
175 fhPtTrkTruePrimGen(0x0),
176 fhPtTrkSecOrFakeRec(0x0),
177 fHRhoUeMedianVsConeGen(0x0),
178 fhEntriesToMedian(0x0),
179 fhEntriesToMedianGen(0x0),
180 fhCellAreaToMedian(0x0),
181 fhCellAreaToMedianGen(0x0),
182 fhNofMultipleTriggersGen(0x0),
183 fhNofMultipleTriggersConeGen(0x0),
184 fhNofMultipleTriggersConeGenLow(0x0),
185 fhNofMultipleTriggersConeGenHigh(0x0),
186 fhDeltaRMultTriggersGenLow(0x0),
187 fhDeltaRMultTriggersGenHigh(0x0),
188 fhDeltaPhiMultTriggersGenLow(0x0),
189 fhDeltaPhiMultTriggersGenHigh(0x0),
190 fhDeltaPhiMultTriggersInclGenLow(0x0),
191 fhDeltaPhiMultTriggersInclGenHigh(0x0),
192 fhInclTrigCounterGen(0x0),
193 fhNofMultipleTriggersConeGenA(0x0),
194 fhNofMultipleTriggersConeGenALow(0x0),
195 fhNofMultipleTriggersConeGenAHigh(0x0),
196 fhDeltaRMultTriggersGenALow(0x0),
197 fhDeltaPhiMultTriggersGenALow(0x0),
198 fhDeltaRMultTriggersGenAHigh(0x0),
199 fhDeltaPhiMultTriggersGenAHigh(0x0),
200 fhNofTriggersGenA(0x0),
201 fhNofMultipleTriggersConeGenB(0x0),
202 fhNofMultipleTriggersConeGenBLow(0x0),
203 fhNofMultipleTriggersConeGenBHigh(0x0),
204 fhDeltaRMultTriggersGenBLow(0x0),
205 fhDeltaPhiMultTriggersGenBLow(0x0),
206 fhDeltaRMultTriggersGenBHigh(0x0),
207 fhDeltaPhiMultTriggersGenBHigh(0x0),
208 fhNofTriggersGenB(0x0),
209 fhTriggerCounterGenLevel(0x0),
210 fhDeltaRMultTriggersGenLevelLow(0x0),
211 fhDeltaPhiMultTriggersGenLevelLow(0x0),
212 fhDeltaRMultTriggersGenLevelHigh(0x0),
213 fhDeltaPhiMultTriggersGenLevelHigh(0x0),
219 fkAcceptance(2.0*TMath::Pi()*1.8),
220 fkDeltaPhiCut(TMath::Pi()-0.8),
226 fh1PtHardTrials(0x0),
229 fEventNumberRangeLow(0),
230 fEventNumberRangeHigh(99),
231 fTriggerPtRangeLow(0.0),
232 fTriggerPtRangeHigh(10000.0),
236 fJetFreeAreaFrac(0.5),
240 fPhiSize(2*TMath::Pi()/fnPhi),
241 fCellArea(fPhiSize*fEtaSize),
243 fDoubleBinning(kFALSE)
245 // default Constructor
248 //---------------------------------------------------------------------
250 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
251 AliAnalysisTaskSE(name),
259 fJetBranchNameChargMC(""),
260 fJetBranchNameKine(""),
261 fJetBranchNameFullMC(""),
262 fJetBranchNameBg(""),
263 fJetBranchNameBgChargMC(""),
264 fJetBranchNameBgKine(""),
267 fListJetsGenFull(0x0),
271 fSystem(0), //pp=0 pPb=1
276 fOfflineTrgMask(AliVEvent::kAny),
287 fTrackLowPtCut(0.15),
288 fUseExchContainer(0),
290 fHistEvtSelection(0x0),
293 fHJetSpecSubUeMedian(0x0),
294 fHJetSpecSubUeCone(0x0),
296 fHJetPhiCorrSubUeMedian(0x0),
297 fHJetPhiCorrSubUeCone(0x0),
300 fHRhoUeMedianVsCone(0x0),
302 //fHJetDensityA4(0x0),
308 fhVertexZAccept(0x0),
310 fhContribVtxAccept(0x0),
311 fhDphiTriggerJet(0x0),
312 fhDphiTriggerJetAccept(0x0),
314 fhCentralityAccept(0x0),
315 fhNofMultipleTriggers(0x0),
316 fhNofMultipleTriggersCone(0x0),
317 fhNofMultipleTriggersConeLow(0x0),
318 fhNofMultipleTriggersConeHigh(0x0),
319 fhDeltaRMultTriggersLow(0x0),
320 fhDeltaRMultTriggersHigh(0x0),
321 fhDeltaPhiMultTriggersLow(0x0),
322 fhDeltaPhiMultTriggersHigh(0x0),
323 fhDeltaPhiMultTriggersInclLow(0x0),
324 fhDeltaPhiMultTriggersInclHigh(0x0),
325 fhInclTrigCounter(0x0),
327 //fHLeadingJetPtRaw(0x0),
328 //fHDphiVsJetPtAll(0x0),
329 fhJetPtGenVsJetPtRec(0x0),
330 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
331 fhJetPtGenVsJetPtRecSubUeCone(0x0),
333 fhJetPtSubUeMedianGen(0x0),
334 fhJetPtSubUeConeGen(0x0),
335 fhJetPtResolutionVsPtGen(0x0),
336 fhJetPtResolutionVsPtConeGen(0x0),
337 fhJetPtGenChargVsJetPtGenFull(0x0),
339 fh2NtriggersGen(0x0),
341 fHJetSpecSubUeMedianGen(0x0),
342 fHJetSpecSubUeConeGen(0x0),
343 fHJetPhiCorrGen(0x0),
344 fHJetPhiCorrSubUeMedianGen(0x0),
345 fHJetPhiCorrSubUeConeGen(0x0),
346 fHJetUeMedianGen(0x0),
348 fhPtTrkTruePrimRec(0x0),
349 fhPtTrkTruePrimGen(0x0),
350 fhPtTrkSecOrFakeRec(0x0),
351 fHRhoUeMedianVsConeGen(0x0),
352 fhEntriesToMedian(0x0),
353 fhEntriesToMedianGen(0x0),
354 fhCellAreaToMedian(0x0),
355 fhCellAreaToMedianGen(0x0),
356 fhNofMultipleTriggersGen(0x0),
357 fhNofMultipleTriggersConeGen(0x0),
358 fhNofMultipleTriggersConeGenLow(0x0),
359 fhNofMultipleTriggersConeGenHigh(0x0),
360 fhDeltaRMultTriggersGenLow(0x0),
361 fhDeltaRMultTriggersGenHigh(0x0),
362 fhDeltaPhiMultTriggersGenLow(0x0),
363 fhDeltaPhiMultTriggersGenHigh(0x0),
364 fhDeltaPhiMultTriggersInclGenLow(0x0),
365 fhDeltaPhiMultTriggersInclGenHigh(0x0),
366 fhInclTrigCounterGen(0x0),
367 fhNofMultipleTriggersConeGenA(0x0),
368 fhNofMultipleTriggersConeGenALow(0x0),
369 fhNofMultipleTriggersConeGenAHigh(0x0),
370 fhDeltaRMultTriggersGenALow(0x0),
371 fhDeltaPhiMultTriggersGenALow(0x0),
372 fhDeltaRMultTriggersGenAHigh(0x0),
373 fhDeltaPhiMultTriggersGenAHigh(0x0),
374 fhNofTriggersGenA(0x0),
375 fhNofMultipleTriggersConeGenB(0x0),
376 fhNofMultipleTriggersConeGenBLow(0x0),
377 fhNofMultipleTriggersConeGenBHigh(0x0),
378 fhDeltaRMultTriggersGenBLow(0x0),
379 fhDeltaPhiMultTriggersGenBLow(0x0),
380 fhDeltaRMultTriggersGenBHigh(0x0),
381 fhDeltaPhiMultTriggersGenBHigh(0x0),
382 fhNofTriggersGenB(0x0),
383 fhTriggerCounterGenLevel(0x0),
384 fhDeltaRMultTriggersGenLevelLow(0x0),
385 fhDeltaPhiMultTriggersGenLevelLow(0x0),
386 fhDeltaRMultTriggersGenLevelHigh(0x0),
387 fhDeltaPhiMultTriggersGenLevelHigh(0x0),
393 fkAcceptance(2.0*TMath::Pi()*1.8),
394 fkDeltaPhiCut(TMath::Pi()-0.8),
400 fh1PtHardTrials(0x0),
403 fEventNumberRangeLow(0),
404 fEventNumberRangeHigh(99),
405 fTriggerPtRangeLow(0.0),
406 fTriggerPtRangeHigh(10000.0),
410 fJetFreeAreaFrac(0.5),
414 fPhiSize(2*TMath::Pi()/fnPhi),
415 fCellArea(fPhiSize*fEtaSize),
417 fDoubleBinning(kFALSE)
420 DefineOutput(1, TList::Class());
423 if(dummy.Contains("KINE")){
424 DefineInput(1, TClonesArray::Class());
425 DefineInput(2, TClonesArray::Class());
429 //--------------------------------------------------------------
430 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
431 AliAnalysisTaskSE(a.GetName()),
435 fAODExtension(a.fAODExtension),
436 fMcEvent(a.fMcEvent),
437 fMcHandler(a.fMcHandler),
438 fJetBranchName(a.fJetBranchName),
439 fJetBranchNameChargMC(a.fJetBranchNameChargMC),
440 fJetBranchNameKine(a.fJetBranchNameKine),
441 fJetBranchNameFullMC(a.fJetBranchNameFullMC),
442 fJetBranchNameBg(a.fJetBranchNameBg),
443 fJetBranchNameBgChargMC(a.fJetBranchNameBgChargMC),
444 fJetBranchNameBgKine(a.fJetBranchNameBgKine),
445 fListJets(a.fListJets),
446 fListJetsGen(a.fListJetsGen),
447 fListJetsGenFull(a.fListJetsGenFull),
448 fListJetsBg(a.fListJetsBg),
449 fListJetsBgGen(a.fListJetsBgGen),
450 fNonStdFile(a.fNonStdFile),
452 fJetParamR(a.fJetParamR),
453 fBgJetParamR(a.fBgJetParamR),
454 fBgMaxJetPt(a.fBgMaxJetPt),
455 fBgConeR(a.fBgConeR),
456 fOfflineTrgMask(a.fOfflineTrgMask),
457 fMinContribVtx(a.fMinContribVtx),
458 fVtxZMin(a.fVtxZMin),
459 fVtxZMax(a.fVtxZMax),
460 fFilterMask(a.fFilterMask),
461 fCentMin(a.fCentMin),
462 fCentMax(a.fCentMax),
463 fJetEtaMin(a.fJetEtaMin),
464 fJetEtaMax(a.fJetEtaMax),
465 fTriggerEtaCut(a.fTriggerEtaCut),
466 fTrackEtaCut(a.fTrackEtaCut),
467 fTrackLowPtCut(a.fTrackLowPtCut),
468 fUseExchContainer(a.fUseExchContainer),
469 fOutputList(a.fOutputList),
470 fHistEvtSelection(a.fHistEvtSelection),
471 fh2Ntriggers(a.fh2Ntriggers),
472 fHJetSpec(a.fHJetSpec),
473 fHJetSpecSubUeMedian(a.fHJetSpecSubUeMedian),
474 fHJetSpecSubUeCone(a.fHJetSpecSubUeCone),
475 fHJetPhiCorr(a.fHJetPhiCorr),
476 fHJetPhiCorrSubUeMedian(a.fHJetPhiCorrSubUeMedian),
477 fHJetPhiCorrSubUeCone(a.fHJetPhiCorrSubUeCone),
478 fHJetUeMedian(a.fHJetUeMedian),
479 fHJetUeCone(a.fHJetUeCone),
480 fHRhoUeMedianVsCone(a.fHRhoUeMedianVsCone),
481 //fHJetDensity(a.fHJetDensity),
482 //fHJetDensityA4(a.fHJetDensityA4),
483 fhJetPhi(a.fhJetPhi),
484 fhTriggerPhi(a.fhTriggerPhi),
485 fhJetEta(a.fhJetEta),
486 fhTriggerEta(a.fhTriggerEta),
487 fhVertexZ(a.fhVertexZ),
488 fhVertexZAccept(a.fhVertexZAccept),
489 fhContribVtx(a.fhContribVtx),
490 fhContribVtxAccept(a.fhContribVtxAccept),
491 fhDphiTriggerJet(a.fhDphiTriggerJet),
492 fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
493 fhCentrality(a.fhCentrality),
494 fhCentralityAccept(a.fhCentralityAccept),
495 fhNofMultipleTriggers(a.fhNofMultipleTriggers),
496 fhNofMultipleTriggersCone(a.fhNofMultipleTriggersCone),
497 fhNofMultipleTriggersConeLow(a.fhNofMultipleTriggersConeLow),
498 fhNofMultipleTriggersConeHigh(a.fhNofMultipleTriggersConeHigh),
499 fhDeltaRMultTriggersLow(a.fhDeltaRMultTriggersLow),
500 fhDeltaRMultTriggersHigh(a.fhDeltaRMultTriggersHigh),
501 fhDeltaPhiMultTriggersLow(a.fhDeltaPhiMultTriggersLow),
502 fhDeltaPhiMultTriggersHigh(a.fhDeltaPhiMultTriggersHigh),
503 fhDeltaPhiMultTriggersInclLow(a.fhDeltaPhiMultTriggersInclLow),
504 fhDeltaPhiMultTriggersInclHigh(a.fhDeltaPhiMultTriggersInclHigh),
505 fhInclTrigCounter(a.fhInclTrigCounter),
506 //fHJetPtRaw(a.fHJetPtRaw),
507 //fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
508 //fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
509 fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
510 fhJetPtGenVsJetPtRecSubUeMedian(a.fhJetPtGenVsJetPtRecSubUeMedian),
511 fhJetPtGenVsJetPtRecSubUeCone(a.fhJetPtGenVsJetPtRecSubUeCone),
512 fhJetPtGen(a.fhJetPtGen),
513 fhJetPtSubUeMedianGen(a.fhJetPtSubUeMedianGen),
514 fhJetPtSubUeConeGen(a.fhJetPtSubUeConeGen),
515 fhJetPtResolutionVsPtGen(a.fhJetPtResolutionVsPtGen),
516 fhJetPtResolutionVsPtConeGen(a.fhJetPtResolutionVsPtConeGen),
517 fhJetPtGenChargVsJetPtGenFull(a.fhJetPtGenChargVsJetPtGenFull),
518 fhJetPtGenFull(a.fhJetPtGenFull),
519 fh2NtriggersGen(a.fh2NtriggersGen),
520 fHJetSpecGen(a.fHJetSpecGen),
521 fHJetSpecSubUeMedianGen(a.fHJetSpecSubUeMedianGen),
522 fHJetSpecSubUeConeGen(a.fHJetSpecSubUeConeGen),
523 fHJetPhiCorrGen(a.fHJetPhiCorrGen),
524 fHJetPhiCorrSubUeMedianGen(a.fHJetPhiCorrSubUeMedianGen),
525 fHJetPhiCorrSubUeConeGen(a.fHJetPhiCorrSubUeConeGen),
526 fHJetUeMedianGen(a.fHJetUeMedianGen),
527 fHJetUeConeGen(a.fHJetUeConeGen),
528 fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
529 fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
530 fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
531 fHRhoUeMedianVsConeGen(a.fHRhoUeMedianVsConeGen),
532 fhEntriesToMedian(a.fhEntriesToMedian),
533 fhEntriesToMedianGen(a.fhEntriesToMedianGen),
534 fhCellAreaToMedian(a.fhCellAreaToMedian),
535 fhCellAreaToMedianGen(a.fhCellAreaToMedianGen),
536 fhNofMultipleTriggersGen(a.fhNofMultipleTriggersGen),
537 fhNofMultipleTriggersConeGen(a.fhNofMultipleTriggersConeGen),
538 fhNofMultipleTriggersConeGenLow(a.fhNofMultipleTriggersConeGenLow),
539 fhNofMultipleTriggersConeGenHigh(a.fhNofMultipleTriggersConeGenHigh),
540 fhDeltaRMultTriggersGenLow(a.fhDeltaRMultTriggersGenLow),
541 fhDeltaRMultTriggersGenHigh(a.fhDeltaRMultTriggersGenHigh),
542 fhDeltaPhiMultTriggersGenLow(a.fhDeltaPhiMultTriggersGenLow),
543 fhDeltaPhiMultTriggersGenHigh(a.fhDeltaPhiMultTriggersGenHigh),
544 fhDeltaPhiMultTriggersInclGenLow(a.fhDeltaPhiMultTriggersInclGenLow),
545 fhDeltaPhiMultTriggersInclGenHigh(a.fhDeltaPhiMultTriggersInclGenHigh),
546 fhInclTrigCounterGen(a.fhInclTrigCounterGen),
547 fhNofMultipleTriggersConeGenA(a.fhNofMultipleTriggersConeGenA),
548 fhNofMultipleTriggersConeGenALow(a.fhNofMultipleTriggersConeGenALow),
549 fhNofMultipleTriggersConeGenAHigh(a.fhNofMultipleTriggersConeGenAHigh),
550 fhDeltaRMultTriggersGenALow(a.fhDeltaRMultTriggersGenALow),
551 fhDeltaPhiMultTriggersGenALow(a.fhDeltaPhiMultTriggersGenALow),
552 fhDeltaRMultTriggersGenAHigh(a.fhDeltaRMultTriggersGenAHigh),
553 fhDeltaPhiMultTriggersGenAHigh(a.fhDeltaPhiMultTriggersGenAHigh),
554 fhNofTriggersGenA(a.fhNofTriggersGenA),
555 fhNofMultipleTriggersConeGenB(a.fhNofMultipleTriggersConeGenB),
556 fhNofMultipleTriggersConeGenBLow(a.fhNofMultipleTriggersConeGenBLow),
557 fhNofMultipleTriggersConeGenBHigh(a.fhNofMultipleTriggersConeGenBHigh),
558 fhDeltaRMultTriggersGenBLow(a.fhDeltaRMultTriggersGenBLow),
559 fhDeltaPhiMultTriggersGenBLow(a.fhDeltaPhiMultTriggersGenBLow),
560 fhDeltaRMultTriggersGenBHigh(a.fhDeltaRMultTriggersGenBHigh),
561 fhDeltaPhiMultTriggersGenBHigh(a.fhDeltaPhiMultTriggersGenBHigh),
562 fhNofTriggersGenB(a.fhNofTriggersGenB),
563 fhTriggerCounterGenLevel(a.fhTriggerCounterGenLevel),
564 fhDeltaRMultTriggersGenLevelLow(a.fhDeltaRMultTriggersGenLevelLow),
565 fhDeltaPhiMultTriggersGenLevelLow(a.fhDeltaPhiMultTriggersGenLevelLow),
566 fhDeltaRMultTriggersGenLevelHigh(a.fhDeltaRMultTriggersGenLevelHigh),
567 fhDeltaPhiMultTriggersGenLevelHigh(a.fhDeltaPhiMultTriggersGenLevelHigh),
568 fIsChargedMC(a.fIsChargedMC),
570 fIsFullMC(a.fIsFullMC),
571 faGenIndex(a.faGenIndex),
572 faRecIndex(a.faRecIndex),
573 fkAcceptance(a.fkAcceptance),
574 fkDeltaPhiCut(a.fkDeltaPhiCut),
576 fh1Trials(a.fh1Trials),
577 fh1AvgTrials(a.fh1AvgTrials),
578 fh1PtHard(a.fh1PtHard),
579 fh1PtHardNoW(a.fh1PtHardNoW),
580 fh1PtHardTrials(a.fh1PtHardTrials),
581 fAvgTrials(a.fAvgTrials),
582 fHardest(a.fHardest),
583 fEventNumberRangeLow(a.fEventNumberRangeLow),
584 fEventNumberRangeHigh(a.fEventNumberRangeHigh),
585 fTriggerPtRangeLow(a.fTriggerPtRangeLow),
586 fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
587 fFillRespMx(a.fFillRespMx),
589 fnTrials(a.fnTrials),
590 fJetFreeAreaFrac(a.fJetFreeAreaFrac),
593 fEtaSize(a.fEtaSize),
594 fPhiSize(a.fPhiSize),
595 fCellArea(a.fCellArea),
596 fSafetyMargin(a.fSafetyMargin),
597 fDoubleBinning(a.fDoubleBinning)
602 //--------------------------------------------------------------
604 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
605 // assignment operator
606 this->~AliAnalysisTaskJetCorePP();
607 new(this) AliAnalysisTaskJetCorePP(a);
610 //--------------------------------------------------------------
612 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
617 delete fListJetsGenFull;
619 delete fListJetsBgGen;
620 delete fOutputList; // ????
624 //--------------------------------------------------------------
627 Bool_t AliAnalysisTaskJetCorePP::Notify()
629 //Implemented Notify() to read the cross sections
630 //and number of trials from pyxsec.root
631 //inspired by AliAnalysisTaskJetSpectrum2::Notify()
632 if(!(fIsChargedMC || fIsKine)) return kFALSE;
633 Float_t xsection = 0;
638 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
640 if(fDebug>1) AliError("ESD not available");
641 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
644 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
647 if(fNonStdFile.Length()!=0){
648 // case that we have an AOD extension we can fetch the jets from the extended output
649 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
650 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
652 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
656 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
659 TFile *curfile = tree->GetCurrentFile();
661 Error("Notify","No current file");
664 if(!fh1Xsec || !fh1Trials){
665 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
668 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
669 fh1Xsec->Fill("<#sigma>",xsection);
670 // construct a poor man average trials
671 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
672 if(trials>=nEntries && nEntries>0.) fAvgTrials = trials/nEntries;
673 fh1Trials->Fill("#sum{ntrials}",trials);
676 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
682 //--------------------------------------------------------------
684 void AliAnalysisTaskJetCorePP::Init()
686 // check for jet branches
687 if(fJetBranchNameKine.Length()==0){
688 if(!strlen(fJetBranchName.Data())){
689 AliError("Jet branch name not set.");
692 fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
697 //--------------------------------------------------------------
699 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
701 // Create histograms and initilize variables
702 fSafetyMargin = fBgConeR*fBgConeR /(fBgJetParamR*fBgJetParamR);
705 fListJets = new TList(); //reconstructed level antikt jets
706 fListJetsBg = new TList(); //reconstructed jets to be removed from bg
708 fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
709 fIsKine = (fJetBranchNameKine.Length()>0) ? kTRUE : kFALSE;
710 fIsFullMC = (fJetBranchNameFullMC.Length()>0) ? kTRUE : kFALSE;
712 fRandom = new TRandom3(0);
714 if(fIsChargedMC || fIsKine){ //full MC or pure kine
715 fListJetsGen = new TList(); //generator level charged antikt jets
716 fListJetsBgGen = new TList(); //generator level jets to be removed from bg
719 fListJetsGenFull = new TList(); //generator level full jets
722 if(!fOutputList) fOutputList = new TList();
723 fOutputList->SetOwner(kTRUE);
725 Bool_t oldStatus = TH1::AddDirectoryStatus();
726 TH1::AddDirectory(kFALSE);
728 fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
729 fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
730 fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
731 fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
732 fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
733 fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
734 fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
736 fOutputList->Add(fHistEvtSelection);
738 Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
739 //trigger pt spectrum (reconstructed)
740 fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
741 nBinsCentrality,0.0,100.0,50,0.0,50.0);
742 if(!fIsKine) fOutputList->Add(fh2Ntriggers);
744 Int_t bw = (fDoubleBinning==0) ? 1 : 2; //make larger bin width
746 //Centrality, A, pTjet, pTtrigg, dphi
747 const Int_t dimSpec = 5;
748 const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 50, bw*110, 50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
749 const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, -20, 0.0, fkDeltaPhiCut};
750 const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0, TMath::Pi()};
751 fHJetSpec = new THnSparseF("fHJetSpec",
752 "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
753 dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
754 if(!fIsKine) fOutputList->Add(fHJetSpec);
756 //background estimated as median of kT jets
757 fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian");
758 fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
759 if(!fIsKine) fOutputList->Add(fHJetSpecSubUeMedian);
760 //background estimated as weighted median of kT jets ala Cone
761 fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone");
762 fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe,pTtrig,dphi]");
763 if(!fIsKine) fOutputList->Add(fHJetSpecSubUeCone);
766 //A, pTjet, pTtrigg, dphi
767 const Int_t dimCor = 4;
768 const Int_t nBinsCor[dimCor] = {50, 110, 50, 100};
769 const Double_t lowBinCor[dimCor] = {0.0, -20, 0.0, -0.5*TMath::Pi()};
770 const Double_t hiBinCor[dimCor] = {1.0, 200, 50.0, 1.5*TMath::Pi()};
771 fHJetPhiCorr = new THnSparseF("fHJetPhiCorr",
772 "Recoil jet spectrum [A,pTjet,pTtrig,dphi]",
773 dimCor,nBinsCor,lowBinCor,hiBinCor);
775 if(!fIsKine) fOutputList->Add(fHJetPhiCorr); // Dphi distribution jet-triger
777 //background estimated as median of kT jets
778 fHJetPhiCorrSubUeMedian = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeMedian");
779 fHJetPhiCorrSubUeMedian->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
780 if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeMedian);
781 //background estimated as weighted median of kT jets ala Cone
782 fHJetPhiCorrSubUeCone = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrSubUeCone");
783 fHJetPhiCorrSubUeCone->SetTitle("Recoil jet spectrum [A,pTjet-pTUe,pTtrig,dphi]");
784 if(!fIsKine) fOutputList->Add(fHJetPhiCorrSubUeCone);
787 //------------------- HISTOS FOR DIAGNOSTIC ----------------------
788 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median
789 const Int_t dimSpecMed = 5;
790 const Int_t nBinsSpecMed[dimSpecMed] = {25, 50, 60, 20, 20};
791 const Double_t lowBinSpecMed[dimSpecMed] = {0.0, 0.0, -20.0, 0.0, 0.0};
792 const Double_t hiBinSpecMed[dimSpecMed] = {1.0, 100.0, 100.0, 10.0, 20.0};
793 fHJetUeMedian = new THnSparseF("fHJetUeMedian",
794 "Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]",
795 dimSpecMed, nBinsSpecMed, lowBinSpecMed, hiBinSpecMed);
796 if(!fIsKine) fOutputList->Add(fHJetUeMedian);
798 //A, pTjet, pTjet-pTUe, pTUe, rhoUe bg estimate from kT median Cone
799 fHJetUeCone = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeCone");
800 fHJetUeCone->SetTitle("Recoil jet spectrum [A,pTjet,pTjet-pTUe, pTUe, rhoUe]");
801 if(!fIsKine) fOutputList->Add(fHJetUeCone);
803 //rho bacground reconstructed data
804 const Int_t dimRho = 2;
805 const Int_t nBinsRho[dimRho] = {50 , 50};
806 const Double_t lowBinRho[dimRho] = {0.0 , 0.0};
807 const Double_t hiBinRho[dimRho] = {20.0 , 20.0};
809 fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",
810 dimRho, nBinsRho, lowBinRho, hiBinRho);
811 if(!fIsKine) fOutputList->Add(fHRhoUeMedianVsCone);
813 //Jet number density histos [Trk Mult, jet density, pT trigger]
814 /*const Int_t dimJetDens = 3;
815 const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10};
816 const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0};
817 const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0};
819 fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
820 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
822 fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
823 dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
825 fOutputList->Add(fHJetDensity);
826 fOutputList->Add(fHJetDensityA4);
829 //inclusive azimuthal and pseudorapidity histograms
830 fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
831 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
832 fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
833 50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
834 fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
835 50,0, 100, 40,-0.9,0.9);
836 fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
837 50, 0, 50, 40,-0.9,0.9);
839 fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
840 fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
841 fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
842 fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
843 fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
844 fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
845 fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
846 fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
847 fhEntriesToMedian = new TH1D("fhEntriesToMedian","fhEntriesToMedian",30,0,30);
848 fhCellAreaToMedian = new TH1D("fhCellAreaToMedian", "fhCellAreaToMedian", 75,0,1.5);
850 fhNofMultipleTriggers = new TH1D("fhNofMultipleTriggers","fhNofMultipleTriggers",20,0,20);
851 fhNofMultipleTriggersCone = new TH1D("fhNofMultipleTriggersCone","fhNofMultipleTriggersCone R<0.4",20,0,20);
852 fhNofMultipleTriggersConeLow = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConepTa15to20");
853 fhNofMultipleTriggersConeHigh = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConepTa20to50");
855 fhDeltaRMultTriggersLow = new TH1D("fhDeltaRMultTriggersLow","fhDeltaRMultTriggersLow", 200,0,4);
856 fhDeltaRMultTriggersHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersHigh");
858 fhDeltaPhiMultTriggersLow = new TH1D("fhDeltaPhiMultTriggersLow","fhDeltaPhiRultTriggers", 320,-3.2,3.2); //single incl trigger
859 fhDeltaPhiMultTriggersHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersHigh"); //single incl trigger
861 fhDeltaPhiMultTriggersInclLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersInclLow");
862 fhDeltaPhiMultTriggersInclHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersInclHigh");
863 fhInclTrigCounter = new TH1D("fhInclTrigCounter","fhInclTrigCounter",1,0,1);
866 fOutputList->Add(fhJetPhi);
867 fOutputList->Add(fhTriggerPhi);
868 fOutputList->Add(fhJetEta);
869 fOutputList->Add(fhTriggerEta);
871 fOutputList->Add(fhVertexZ);
872 fOutputList->Add(fhVertexZAccept);
874 fOutputList->Add(fhContribVtx);
875 fOutputList->Add(fhContribVtxAccept);
876 fOutputList->Add(fhDphiTriggerJet);
877 fOutputList->Add(fhDphiTriggerJetAccept);
878 fOutputList->Add(fhCentrality);
879 fOutputList->Add(fhCentralityAccept);
880 fOutputList->Add(fhEntriesToMedian);
881 fOutputList->Add(fhCellAreaToMedian);
882 fOutputList->Add(fhNofMultipleTriggers);
883 fOutputList->Add(fhNofMultipleTriggersCone);
884 fOutputList->Add(fhNofMultipleTriggersConeLow);
885 fOutputList->Add(fhNofMultipleTriggersConeHigh);
886 fOutputList->Add(fhDeltaRMultTriggersLow);
887 fOutputList->Add(fhDeltaRMultTriggersHigh);
888 fOutputList->Add(fhDeltaPhiMultTriggersLow);
889 fOutputList->Add(fhDeltaPhiMultTriggersHigh);
890 fOutputList->Add(fhDeltaPhiMultTriggersInclLow);
891 fOutputList->Add(fhDeltaPhiMultTriggersInclHigh);
892 fOutputList->Add(fhInclTrigCounter);
894 // raw spectra of INCLUSIVE jets
895 //Centrality, pTjet, A
896 /*const Int_t dimRaw = 3;
897 const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100};
898 const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0};
899 const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0};
900 fHJetPtRaw = new THnSparseF("fHJetPtRaw",
901 "Incl. jet spectrum [cent,pTjet,A]",
902 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
903 fOutputList->Add(fHJetPtRaw);
905 // raw spectra of LEADING jets
906 //Centrality, pTjet, A
907 fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
908 "Leading jet spectrum [cent,pTjet,A]",
909 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
910 fOutputList->Add(fHLeadingJetPtRaw);
912 // Dphi versus pT jet
913 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
914 const Int_t dimDp = 4;
915 const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50};
916 const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0};
917 const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0};
918 fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
919 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
920 dimDp,nBinsDp,lowBinDp,hiBinDp);
921 fOutputList->Add(fHDphiVsJetPtAll);
924 //analyze MC generator level
925 if(fIsChargedMC || fIsKine){
927 //Fill response matrix only once
928 fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", bw*100,0,200, bw*100,0,200);
929 fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC charg jet pt spectrum versus rec charged jet pt spectrum
931 fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", bw*110,-20,200, bw*110,-20,200);
932 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr
934 fhJetPtGenVsJetPtRecSubUeCone=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCone");
935 fhJetPtGenVsJetPtRecSubUeCone->SetTitle("fhJetPtGenVsJetPtRecSubUeCone");
936 fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCone); // with weighted kT median bg subtr
938 fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",bw*100,0,200); //MC generator charged jet pt spectrum
939 fOutputList->Add(fhJetPtGen);
941 fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",bw*110,-20,200);
942 fOutputList->Add(fhJetPtSubUeMedianGen); // with kT median bg subtr
944 fhJetPtSubUeConeGen = (TH1D*) fhJetPtSubUeMedianGen->Clone("fhJetPtSubUeConeGen");
945 fOutputList->Add(fhJetPtSubUeConeGen); // with weighted kT median bg subtr
947 fhJetPtResolutionVsPtGen = new TH2D("fhJetPtResolutionVsPtGen","fhJetPtResolutionVsPtGen",20,0,100, 35,-1.,0.4);
948 fOutputList->Add(fhJetPtResolutionVsPtGen); // with weighted kT median bg subtr
950 fhJetPtResolutionVsPtConeGen = (TH2D*) fhJetPtResolutionVsPtGen->Clone("fhJetPtResolutionVsPtConeGen");
951 fOutputList->Add(fhJetPtResolutionVsPtConeGen); // with weighted kT median bg subtr
955 fhJetPtGenChargVsJetPtGenFull = new TH2D("fhJetPtGenChargVsJetPtGenFull","fhJetPtGenChargVsJetPtGenFull", 100,0,200, 100,0,200);
956 fOutputList->Add(fhJetPtGenChargVsJetPtGenFull); //gen full MC jet pt versus gen charged jet MC pt
958 fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",100,0,200); //MC generator full jet pt spectrum
959 fOutputList->Add(fhJetPtGenFull);
963 fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
964 fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
965 fOutputList->Add(fh2NtriggersGen);
967 //Centrality, A, pT, pTtrigg
968 fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
969 fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
970 fOutputList->Add(fHJetSpecGen);
972 fHJetSpecSubUeMedianGen = (THnSparseF*) fHJetSpecSubUeMedian->Clone("fHJetSpecSubUeMedianGen");
973 fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
974 fOutputList->Add(fHJetSpecSubUeMedianGen);
976 fHJetSpecSubUeConeGen = (THnSparseF*) fHJetSpecSubUeCone->Clone("fHJetSpecSubUeConeGen");
977 fHJetSpecSubUeConeGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCone->GetTitle()));
978 fOutputList->Add(fHJetSpecSubUeConeGen);
980 fHJetPhiCorrGen = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrGen");
981 fHJetPhiCorrGen->SetTitle(Form("%s Gen MC",fHJetPhiCorr->GetTitle()));
982 fOutputList->Add(fHJetPhiCorrGen);
984 fHJetPhiCorrSubUeMedianGen = (THnSparseF*) fHJetPhiCorrSubUeMedian->Clone("fHJetPhiCorrSubUeMedianGen");
985 fHJetPhiCorrSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetPhiCorrSubUeMedian->GetTitle()));
986 fOutputList->Add(fHJetPhiCorrSubUeMedianGen);
988 fHJetPhiCorrSubUeConeGen = (THnSparseF*) fHJetPhiCorrSubUeCone->Clone("fHJetPhiCorrSubUeConeGen");
989 fHJetPhiCorrSubUeConeGen->SetTitle(Form("%s Gen MC", fHJetPhiCorrSubUeCone->GetTitle()));
990 fOutputList->Add(fHJetPhiCorrSubUeConeGen);
993 fHJetUeMedianGen = (THnSparseF*) fHJetUeMedian->Clone("fHJetUeMedianGen");
994 fHJetUeMedianGen->SetTitle(Form("%s Gen MC", fHJetUeMedian->GetTitle()));
995 fOutputList->Add(fHJetUeMedianGen);
997 fHJetUeConeGen = (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen");
998 fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle()));
999 fOutputList->Add(fHJetUeConeGen);
1002 //track efficiency/contamination histograms eta versus pT
1003 Double_t bins [] = {0, 0.2,0.4,0.6, 0.8, 1., 1.2, 1.4, 1.6, 1.8, 2., 2.5, 3., 3.5, 4., 5., 6., 8., 10., 20., 50.};
1004 Int_t nbins = sizeof(bins)/sizeof(Double_t)-1;
1005 fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",nbins, bins, 18,-0.9,0.9);
1006 fOutputList->Add(fhPtTrkTruePrimRec);
1008 fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
1009 fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
1010 fOutputList->Add(fhPtTrkTruePrimGen);
1012 fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
1013 fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
1014 fOutputList->Add(fhPtTrkSecOrFakeRec);
1017 fHRhoUeMedianVsConeGen = (THnSparseF*) fHRhoUeMedianVsCone->Clone("hRhoUeMedianVsConeGen");
1018 fHRhoUeMedianVsConeGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCone->GetTitle()));
1019 fOutputList->Add(fHRhoUeMedianVsConeGen);
1021 fhEntriesToMedianGen = (TH1D*) fhEntriesToMedian->Clone("fhEntriesToMedianGen");
1022 fhEntriesToMedianGen->SetTitle(Form("%s Gen MC", fhEntriesToMedian->GetTitle()));
1023 fOutputList->Add(fhEntriesToMedianGen);
1025 fhCellAreaToMedianGen = (TH1D*) fhCellAreaToMedian->Clone("fhCellAreaToMedianGen");
1026 fhCellAreaToMedianGen->SetTitle(Form("%s Gen MC", fhCellAreaToMedian->GetTitle()));
1027 fOutputList->Add(fhCellAreaToMedianGen);
1029 fhNofMultipleTriggersGen = (TH1D*) fhNofMultipleTriggers->Clone("fhNofMultipleTriggersGen");
1030 fOutputList->Add(fhNofMultipleTriggersGen);
1032 fhNofMultipleTriggersConeGen = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGen");
1033 fOutputList->Add(fhNofMultipleTriggersConeGen);
1035 fhNofMultipleTriggersConeGenLow = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGenpta15to20");
1036 fOutputList->Add(fhNofMultipleTriggersConeGenLow);
1038 fhNofMultipleTriggersConeGenHigh = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGenpta20to50");
1039 fOutputList->Add(fhNofMultipleTriggersConeGenHigh);
1041 fhDeltaRMultTriggersGenLow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenLow");
1042 fOutputList->Add(fhDeltaRMultTriggersGenLow);
1044 fhDeltaRMultTriggersGenHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenHigh");
1045 fOutputList->Add(fhDeltaRMultTriggersGenHigh);
1047 fhDeltaPhiMultTriggersGenLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenLow");
1048 fOutputList->Add(fhDeltaPhiMultTriggersGenLow);
1050 fhDeltaPhiMultTriggersGenHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenHigh");
1051 fOutputList->Add(fhDeltaPhiMultTriggersGenHigh);
1053 fhDeltaPhiMultTriggersInclGenLow = (TH1D*) fhDeltaPhiMultTriggersInclLow->Clone("fhDeltaPhiMultTriggersInclGenLow");
1054 fOutputList->Add(fhDeltaPhiMultTriggersInclGenLow);
1056 fhDeltaPhiMultTriggersInclGenHigh = (TH1D*) fhDeltaPhiMultTriggersInclHigh->Clone("fhDeltaPhiMultTriggersInclGenHigh");
1057 fOutputList->Add(fhDeltaPhiMultTriggersInclGenHigh);
1059 fhInclTrigCounterGen= (TH1D*) fhInclTrigCounter->Clone("fhInclTrigCounterGen");
1060 fOutputList->Add(fhInclTrigCounterGen);
1062 fhNofMultipleTriggersConeGenA = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10");
1063 fOutputList->Add(fhNofMultipleTriggersConeGenA);
1065 fhNofMultipleTriggersConeGenALow = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10pTa15to20");
1066 fOutputList->Add(fhNofMultipleTriggersConeGenALow);
1068 fhNofMultipleTriggersConeGenAHigh = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10pTa20to50");
1069 fOutputList->Add(fhNofMultipleTriggersConeGenAHigh);
1071 fhDeltaRMultTriggersGenALow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen10pTa15to20");
1072 fOutputList->Add(fhDeltaRMultTriggersGenALow);
1074 fhDeltaPhiMultTriggersGenALow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen10pTa15to20");
1075 fOutputList->Add(fhDeltaPhiMultTriggersGenALow);
1077 fhDeltaRMultTriggersGenAHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen10pTa20to50");
1078 fOutputList->Add(fhDeltaRMultTriggersGenAHigh);
1080 fhDeltaPhiMultTriggersGenAHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen10pTa20to50");
1081 fOutputList->Add(fhDeltaPhiMultTriggersGenAHigh);
1083 fhNofTriggersGenA = (TH1D*) fhInclTrigCounter->Clone("fhNofTriggersGen10");
1084 fOutputList->Add(fhNofTriggersGenA);
1087 fhNofMultipleTriggersConeGenB = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5");
1088 fOutputList->Add(fhNofMultipleTriggersConeGenB);
1090 fhNofMultipleTriggersConeGenBLow = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5pTa15to20");
1091 fOutputList->Add(fhNofMultipleTriggersConeGenBLow);
1093 fhNofMultipleTriggersConeGenBHigh = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5pTa20to50");
1094 fOutputList->Add(fhNofMultipleTriggersConeGenBHigh);
1096 fhDeltaRMultTriggersGenBLow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen5pTa15to20");
1097 fOutputList->Add(fhDeltaRMultTriggersGenBLow);
1099 fhDeltaPhiMultTriggersGenBLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen5pTa15to20");
1100 fOutputList->Add(fhDeltaPhiMultTriggersGenBLow);
1102 fhDeltaRMultTriggersGenBHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen5pTa20to50");
1103 fOutputList->Add(fhDeltaRMultTriggersGenBHigh);
1105 fhDeltaPhiMultTriggersGenBHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen5pTa20to50");
1106 fOutputList->Add(fhDeltaPhiMultTriggersGenBHigh);
1108 fhNofTriggersGenB = (TH1D*) fhInclTrigCounter->Clone("fhNofTriggersGen5");
1109 fOutputList->Add(fhNofTriggersGenB);
1114 if(fIsChargedMC){ // Filled from gen level if well reconstructed trigger exists
1115 fhTriggerCounterGenLevel = (TH1D*) fhInclTrigCounter->Clone("fhTriggerCounterGenLevel");
1116 fOutputList->Add(fhTriggerCounterGenLevel);
1118 fhDeltaRMultTriggersGenLevelLow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenLevelLow");
1119 fOutputList->Add(fhDeltaRMultTriggersGenLevelLow);
1121 fhDeltaPhiMultTriggersGenLevelLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenLevelLow");
1122 fOutputList->Add(fhDeltaPhiMultTriggersGenLevelLow);
1124 fhDeltaRMultTriggersGenLevelHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenLevelHigh");
1125 fOutputList->Add(fhDeltaRMultTriggersGenLevelHigh);
1127 fhDeltaPhiMultTriggersGenLevelHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenLevelHigh");
1128 fOutputList->Add(fhDeltaPhiMultTriggersGenLevelHigh);
1131 //-------------------------------------
1132 // pythia histograms
1133 const Int_t nBinPt = 150;
1134 Double_t binLimitsPt[nBinPt+1];
1135 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
1137 binLimitsPt[iPt] = -50.0;
1139 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
1143 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
1144 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1145 fOutputList->Add(fh1Xsec);
1146 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
1147 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1148 fOutputList->Add(fh1Trials);
1149 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
1150 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
1151 fOutputList->Add(fh1AvgTrials);
1152 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
1153 fOutputList->Add(fh1PtHard);
1154 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
1155 fOutputList->Add(fh1PtHardNoW);
1156 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
1157 fOutputList->Add(fh1PtHardTrials);
1160 // =========== Switch on Sumw2 for all histos ===========
1161 for(Int_t i=0; i<fOutputList->GetEntries(); i++){
1162 TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
1167 THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
1172 TH1::AddDirectory(oldStatus);
1174 PostData(1, fOutputList);
1177 //--------------------------------------------------------------------
1179 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
1185 Double_t eventW = 1.0;
1186 Double_t ptHard = 0.0;
1187 Double_t nTrials = 1.0; // Trials for MC trigger
1188 if(fIsChargedMC || fIsKine) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
1190 if(TMath::Abs((Float_t) fJetParamR)<0.00001){
1191 AliError("ANTIKT Cone radius is set to zero.");
1195 if(TMath::Abs((Float_t) fBgJetParamR)<0.00001){
1196 AliError("ANTIKT Cone radius is set to zero.");
1200 if(!fIsKine){ //real data or full MC
1201 if(!strlen(fJetBranchName.Data())){
1202 AliError("Name of jet branch not set.");
1205 if(!strlen(fJetBranchNameBg.Data())){
1206 AliError("Name of jet bg branch not set.");
1210 if(!strlen(fJetBranchNameBgKine.Data())){
1211 AliError("Name of jet bg branch for kine not set.");
1217 fMcEvent = fMcHandler->MCEvent();
1219 if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcHandler=NULL\n");
1220 PostData(1, fOutputList);
1224 if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcEvent=NULL \n");
1225 PostData(1, fOutputList);
1229 Float_t xsection = 0;
1232 AliGenPythiaEventHeader *genPH =
1233 dynamic_cast<AliGenPythiaEventHeader*> (fMcEvent->GenEventHeader());
1235 xsection = genPH->GetXsection();
1236 trials = genPH->Trials();
1237 ptHard = genPH->GetPtHard();
1239 fh1Xsec->Fill("<#sigma>",xsection);
1240 fh1Trials->Fill("#sum{ntrials}",trials);
1241 fh1PtHard->Fill(ptHard,eventW);
1242 fh1PtHardNoW->Fill(ptHard,1);
1243 fh1PtHardTrials->Fill(ptHard,trials);
1247 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
1249 if(fDebug>1) AliError("ESD not available");
1250 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1253 fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
1254 AliAODEvent* aod = NULL;
1255 // take all other information from the aod we take the tracks from
1256 if(!aod && !fIsKine){
1257 if(!fESD) aod = fAODIn; //not ESD and not kine => input AOD
1258 else aod = fAODOut;// ESD or kine
1262 if(fNonStdFile.Length()!=0){
1263 // case that we have an AOD extension we can fetch the jets from the extended output
1264 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1265 fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
1267 if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
1271 // ----------------- EVENT SELECTION --------------------------
1272 fHistEvtSelection->Fill(1); // number of events before event selection
1275 // physics selection
1276 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
1277 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1279 if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
1280 if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
1281 fHistEvtSelection->Fill(2);
1282 PostData(1, fOutputList);
1288 if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1289 fHistEvtSelection->Fill(3);
1290 PostData(1, fOutputList);
1294 // vertex selection for reconstructed data
1295 AliAODVertex* primVtx = aod->GetPrimaryVertex();
1298 if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
1299 fHistEvtSelection->Fill(3);
1300 PostData(1, fOutputList);
1304 Int_t nTracksPrim = primVtx->GetNContributors();
1305 Float_t vtxz = primVtx->GetZ();
1307 fhContribVtx->Fill(nTracksPrim);
1308 if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
1310 if((nTracksPrim < fMinContribVtx) ||
1311 (vtxz < fVtxZMin) ||
1313 if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
1314 (char*)__FILE__,__LINE__,vtxz);
1315 fHistEvtSelection->Fill(3);
1316 PostData(1, fOutputList);
1320 fhContribVtxAccept->Fill(nTracksPrim);
1321 fhVertexZAccept->Fill(vtxz);
1323 }else{ //KINE cut on vertex
1324 const AliVVertex *vtxMC = fMcEvent->GetPrimaryVertex();
1325 Float_t zVtx = vtxMC->GetZ();
1326 fhVertexZ->Fill(zVtx);
1327 if(zVtx < fVtxZMin || zVtx > fVtxZMax){ //vertex cut
1328 fHistEvtSelection->Fill(3);
1329 PostData(1, fOutputList);
1332 fhVertexZAccept->Fill(zVtx);
1335 //FK// No event class selection imposed
1336 // event class selection (from jet helper task)
1337 //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
1338 //if(fDebug) Printf("Event class %d", eventClass);
1339 //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
1340 // fHistEvtSelection->Fill(4);
1341 // PostData(1, fOutputList);
1345 //------------------ CENTRALITY SELECTION ---------------
1346 AliCentrality *cent = 0x0;
1347 Double_t centValue = 0.0;
1348 if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb
1350 cent = fESD->GetCentrality();
1351 if(cent) centValue = cent->GetCentralityPercentile("V0M");
1353 //centValue = aod->GetHeader()->GetCentrality();
1354 centValue = ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
1356 if(fDebug) printf("centrality: %f\n", centValue);
1358 fhCentrality->Fill(centValue);
1360 if(centValue < fCentMin || centValue > fCentMax){
1361 fHistEvtSelection->Fill(4);
1362 PostData(1, fOutputList);
1366 fhCentralityAccept->Fill( centValue );
1370 //-----------------select disjunct event subsamples ----------------
1371 if(!fIsKine){ //reconstructed data
1372 //Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile();
1373 AliAODHeader * header = dynamic_cast<AliAODHeader*>(aod->GetHeader());
1374 if(!header) AliFatal("Not a standard AOD");
1376 Int_t eventnum = header->GetEventNumberESDFile();
1377 Int_t lastdigit = eventnum % 10;
1378 if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
1379 fHistEvtSelection->Fill(5);
1380 PostData(1, fOutputList);
1385 if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
1386 fHistEvtSelection->Fill(0); // accepted events
1387 // ==================== end event selection ============================
1389 Double_t tmpArrayFive[5];
1390 Double_t tmpArrayFour[4];
1393 TList particleList; //list of tracks
1394 Int_t indexTrigg = -1;
1395 Double_t rhoFromCellMedian=0.0, rhoCone=0.0;
1398 //=============== Reconstructed antikt jets ===============
1399 ReadTClonesArray(fJetBranchName.Data() , fListJets);
1400 ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg);
1402 //============ Estimate background in reconstructed events ===========
1404 //Find Hadron trigger and Estimate rho from cone
1405 indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1406 EstimateBgCone(fListJets, &particleList, rhoCone);
1408 //Estimate rho from cell median minus jets
1409 EstimateBgRhoMedian(fListJetsBg, &particleList, rhoFromCellMedian,0); //real data
1411 //============== analyze generator level MC ================
1412 TList particleListGen; //list of tracks in MC
1414 if(fIsChargedMC || fIsKine){
1416 if(fIsKine){ //pure kine
1418 //================= generated charged antikt jets ================
1419 ReadTClonesArray(fJetBranchNameKine.Data(), fListJetsGen);
1420 ReadTClonesArray(fJetBranchNameBgKine.Data(), fListJetsBgGen);
1422 //================= generated charged antikt jets ================
1423 ReadTClonesArray(fJetBranchNameChargMC.Data(), fListJetsGen);
1424 ReadTClonesArray(fJetBranchNameBgChargMC.Data(), fListJetsBgGen);
1426 if(fIsFullMC){ //generated full jets
1427 ReadTClonesArray(fJetBranchNameFullMC.Data(), fListJetsGenFull);
1430 //========================================================
1431 //serarch for charged trigger at the MC generator level
1433 Int_t indexTriggGen = -1;
1434 Double_t ptTriggGen = -1;
1435 Int_t iCounterGen = 0; //number of entries in particleListGen array
1436 Int_t triggersMC[200];//list of trigger candidates
1437 Int_t ntriggersMC = 0; //index in triggers array
1438 Int_t triggersMCa[200]; //list of trigger candidates 10%eloss
1439 Int_t ntriggersMCa = 0; //index in triggers array 10%eloss
1440 Int_t triggersMCb[200]; //list of trigger candidates 5%eloss
1441 Int_t ntriggersMCb = 0; //index in triggers array 5%eloss
1445 if(fESD){//ESD input
1447 AliMCEvent* mcEvent = MCEvent();
1449 PostData(1, fOutputList);
1453 AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
1454 if(pythiaGenHeader){
1455 nTrials = pythiaGenHeader->Trials();
1456 ptHard = pythiaGenHeader->GetPtHard();
1457 fh1PtHard->Fill(ptHard,eventW);
1458 fh1PtHardNoW->Fill(ptHard,1);
1459 fh1PtHardTrials->Fill(ptHard,nTrials);
1462 for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
1463 if(!mcEvent->IsPhysicalPrimary(it)) continue;
1464 AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
1466 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1468 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1469 if(indexTriggGen > -1){//trigger candidate was found
1470 triggersMC[ntriggersMC] = indexTriggGen;
1475 iCounterGen++;//index in particleListGen array
1480 PostData(1, fOutputList);
1483 TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
1485 PostData(1, fOutputList);
1488 for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
1489 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1491 if(!part->IsPhysicalPrimary()) continue;
1492 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1494 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1495 if(indexTriggGen > -1){ //trigger candidater was found
1496 triggersMC[ntriggersMC] = indexTriggGen;
1501 iCounterGen++;//number of entries in particleListGen array
1505 }else{ //analyze kine
1507 for(Int_t it = 0; it < fMcEvent->GetNumberOfTracks(); it++){
1508 if(!fMcEvent->IsPhysicalPrimary(it)) continue;
1509 AliMCParticle* part = (AliMCParticle*) fMcEvent->GetTrack(it);
1511 if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1513 if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1514 if(indexTriggGen > -1){//trigger candidate was found
1515 triggersMC[ntriggersMC] = indexTriggGen;
1520 iCounterGen++;//index in particleListGen array
1526 Int_t npar = particleListGen.GetEntries();
1527 for(Int_t ip=0; ip < npar; ip++){
1528 AliVParticle *part = (AliVParticle*) particleListGen.At(ip);
1531 Double_t pta = 0.9 * part->Pt(); //10% energy loss
1532 Double_t ptb = 0.95 * part->Pt(); //5% energy loss
1533 if(fTriggerPtRangeLow <= pta && pta < fTriggerPtRangeHigh &&
1534 TMath::Abs(part->Eta()) < fTriggerEtaCut && ntriggersMCa<200){
1535 triggersMCa[ntriggersMCa] = ip;
1539 if(fTriggerPtRangeLow <= ptb && ptb < fTriggerPtRangeHigh &&
1540 TMath::Abs(part->Eta()) < fTriggerEtaCut && ntriggersMCb<200){
1541 triggersMCb[ntriggersMCb] = ip;
1547 Int_t rnda = fRandom->Integer(ntriggersMCa); //0 to ntriggers-1
1548 Int_t indexTriggGena = triggersMCa[rnda];
1550 Double_t deltaPhia, deltaEtaa, deltaRa;
1554 //Correlation with single inclusive TRIGGER
1555 AliVParticle* tGenTa = (AliVParticle*) particleListGen.At(indexTriggGena);
1557 fhNofTriggersGenA->Fill(0.5); // 15-50
1558 //for(Int_t ia=0; ia<ntriggersMCa; ia++)
1559 for(Int_t ia=0; ia< npar; ia++){
1560 //if(indexTriggGena == triggersMCa[ia]) continue;
1561 if(indexTriggGena == ia) continue;
1563 //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCa[ia]);
1564 AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ia);
1565 if(!tGenTz) continue;
1566 if(tGenTz->Pt()*0.9<15.0) continue;
1567 if(tGenTz->Pt()*0.9>50.0) continue;
1569 deltaPhia = RelativePhi(tGenTa->Phi(),tGenTz->Phi());
1570 deltaEtaa = tGenTa->Eta()-tGenTz->Eta();
1571 deltaRa = sqrt(deltaPhia*deltaPhia + deltaEtaa*deltaEtaa);
1573 if(tGenTz->Pt()*0.9<20.0){
1574 fhDeltaRMultTriggersGenALow->Fill(deltaRa);
1575 fhDeltaPhiMultTriggersGenALow->Fill(deltaPhia);
1577 fhDeltaRMultTriggersGenAHigh->Fill(deltaRa);
1578 fhDeltaPhiMultTriggersGenAHigh->Fill(deltaPhia);
1582 if(tGenTz->Pt()*0.9<20.0) aaLow++; //15-20
1583 else aaHigh++; //20-50
1587 fhNofMultipleTriggersConeGenA->Fill(aaLow+aaHigh); // 15-50
1588 fhNofMultipleTriggersConeGenALow->Fill(aaLow); //15-20
1589 fhNofMultipleTriggersConeGenAHigh->Fill(aaHigh);//20-50
1593 Int_t rndb = fRandom->Integer(ntriggersMCb); //0 to ntriggers-1
1594 Int_t indexTriggGenb = triggersMCb[rndb];
1596 Double_t deltaPhib, deltaEtab, deltaRb;
1600 //Correlation with single inclusive TRIGGER
1601 AliVParticle* tGenTb = (AliVParticle*) particleListGen.At(indexTriggGenb);
1603 fhNofTriggersGenB->Fill(0.5); // 20-50
1604 //for(Int_t ib=0; ib<ntriggersMCb; ib++)
1605 for(Int_t ib=0; ib<npar; ib++){
1606 //if(indexTriggGenb == triggersMCb[ib]) continue;
1607 if(indexTriggGenb == ib) continue;
1609 //AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCb[ib]);
1610 AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(ib);
1611 if(!tGenTz) continue;
1612 if(tGenTz->Pt()*0.95<15.0) continue;
1613 if(tGenTz->Pt()*0.95>50.0) continue;
1615 deltaPhib = RelativePhi(tGenTb->Phi(),tGenTz->Phi());
1616 deltaEtab = tGenTb->Eta()-tGenTz->Eta();
1617 deltaRb = sqrt(deltaPhib*deltaPhib + deltaEtab*deltaEtab);
1619 if(tGenTz->Pt()*0.95<20.0){
1620 fhDeltaRMultTriggersGenBLow->Fill(deltaRb);
1621 fhDeltaPhiMultTriggersGenBLow->Fill(deltaPhib);
1623 fhDeltaRMultTriggersGenBHigh->Fill(deltaRb);
1624 fhDeltaPhiMultTriggersGenBHigh->Fill(deltaPhib);
1628 if(tGenTz->Pt()*0.95<20.0) bbLow++; //15-20
1633 fhNofMultipleTriggersConeGenB->Fill(bbLow+bbHigh); //15-50
1634 fhNofMultipleTriggersConeGenBLow->Fill(bbLow);//15-20
1635 fhNofMultipleTriggersConeGenBHigh->Fill(bbHigh);//20-50
1641 //============== Estimate bg in generated events ==============
1642 Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
1644 //Estimate rho from cone
1645 EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
1647 //Estimate rho from cell median minus jets
1648 EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
1650 if(ntriggersMC > 0){ //select random inclusive trigger
1651 fhInclTrigCounterGen->Fill(0.5,ntriggersMC); //count inclusive triggers
1653 //Lower pTa bin 15-20
1654 for(Int_t it1=0; it1<ntriggersMC; it1++){
1655 AliVParticle* tGent1 = (AliVParticle*) particleListGen.At(triggersMC[it1]);
1656 if(!tGent1) continue;
1657 for(Int_t ia=0; ia<particleListGen.GetEntries(); ia++){
1658 if(ia == triggersMC[it1]) continue;
1659 AliVParticle* tGent2 = (AliVParticle*) particleListGen.At(ia);
1660 if(!tGent2) continue;
1661 if(tGent2->Pt()<15.0) continue;
1662 if(tGent2->Pt()>20.0) continue;
1663 fhDeltaPhiMultTriggersInclGenLow->Fill(RelativePhi(tGent1->Phi(),tGent2->Phi()));
1666 //Higher pTa bin 20-50
1667 for(Int_t it1=0; it1<ntriggersMC-1; it1++){
1668 AliVParticle* tGent1 = (AliVParticle*) particleListGen.At(triggersMC[it1]);
1669 if(!tGent1) continue;
1671 for(Int_t it2=it1+1; it2<ntriggersMC; it2++){
1672 AliVParticle* tGent2 = (AliVParticle*) particleListGen.At(triggersMC[it2]);
1673 if(!tGent2) continue;
1674 fhDeltaPhiMultTriggersInclGenHigh->Fill(RelativePhi(tGent1->Phi(),tGent2->Phi()));
1680 //============ Generator trigger+jet ==================
1681 if(fHardest==0){ //single inclusive trigger
1682 if(ntriggersMC>0){ //there is at least one trigger
1683 Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
1684 indexTriggGen = triggersMC[rnd];
1686 fhNofMultipleTriggersGen->Fill(ntriggersMC-1);
1688 Double_t deltaPhi, deltaEta, deltaR;
1692 //Correlation with single inclusive TRIGGER
1693 AliVParticle* tGenT1 = (AliVParticle*) particleListGen.At(indexTriggGen);
1695 //for(Int_t ia=0; ia<ntriggersMC; ia++)
1696 for(Int_t ia=0; ia< particleListGen.GetEntries(); ia++){
1697 //if(indexTriggGen == triggersMC[ia]) continue;
1698 if(indexTriggGen == ia) continue;
1700 //AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1701 AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(ia);
1702 if(!tGenT2) continue;
1703 if(tGenT2->Pt()<15.0) continue;
1704 if(tGenT2->Pt()>50.0) continue;
1705 deltaPhi = RelativePhi(tGenT1->Phi(),tGenT2->Phi());
1706 deltaEta = tGenT1->Eta()-tGenT2->Eta();
1707 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1709 if(tGenT2->Pt()<20.0){ //single inclusive trigger + assoc 15-20
1710 fhDeltaRMultTriggersGenLow->Fill(deltaR);
1711 fhDeltaPhiMultTriggersGenLow->Fill(deltaPhi);
1712 }else{ //single inclusive trigger + assoc 20-50
1713 fhDeltaRMultTriggersGenHigh->Fill(deltaR);
1714 fhDeltaPhiMultTriggersGenHigh->Fill(deltaPhi);
1718 if(tGenT2->Pt()<20.0) iLow++;
1723 fhNofMultipleTriggersConeGen->Fill(iLow+iHigh);
1724 fhNofMultipleTriggersConeGenLow->Fill(iLow);
1725 fhNofMultipleTriggersConeGenHigh->Fill(iHigh);
1727 //Single inclusive trigger within |Delta eta
1731 indexTriggGen = -1; //trigger not found
1736 Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1737 Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
1738 Bool_t fillOnceGen = kTRUE;
1741 for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1742 indexTriggGen = igen; //trigger hadron
1744 if(indexTriggGen == -1) continue;
1745 AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
1746 if(!triggerGen) continue;
1749 if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
1751 if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1753 ptTriggGen = triggerGen->Pt(); //count triggers
1754 fh2NtriggersGen->Fill(centValue, ptTriggGen);
1756 //Count jets and trigger-jet pairs at MC generator level
1757 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1758 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1760 Double_t etaJetGen = jet->Eta();
1761 Double_t areaJetGen = jet->EffectiveAreaCharged();
1763 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1765 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1766 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1768 //Rongrong's analysis
1769 Double_t dPhiGen = jet->Phi() - triggerGen->Phi();
1770 if(dPhiGen>2*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1771 if(dPhiGen<-2*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1772 if(dPhiGen<-0.5*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1773 if(dPhiGen>1.5*TMath::Pi()) dPhiGen -= 2*TMath::Pi();
1775 //[A,pTjet,pTtrig,dphi]
1776 tmpArrayFour[0] = areaJetGen;
1777 tmpArrayFour[1] = jet->Pt();
1778 tmpArrayFour[2] = ptTriggGen;
1779 tmpArrayFour[3] = dPhiGen;
1781 fHJetPhiCorrGen->Fill(tmpArrayFour); // Dphi distribution jet-triger
1783 //[A,pTjet-pTUe,pTtrig,dphi]
1784 tmpArrayFour[1] = jet->Pt() - ptUeFromCellMedianGen;
1785 fHJetPhiCorrSubUeMedianGen->Fill(tmpArrayFour);
1787 //[A,pTjet-pTUe,pTtrig,dphi]
1788 tmpArrayFour[1] = jet->Pt() - ptUeConeGen;
1789 fHJetPhiCorrSubUeConeGen->Fill(tmpArrayFour);
1791 //Leticia's analysis
1792 Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
1793 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1795 //Centrality, A, pT, pTtrigg
1796 tmpArrayFive[0] = centValue;
1797 tmpArrayFive[1] = areaJetGen;
1798 tmpArrayFive[2] = jet->Pt();
1799 tmpArrayFive[3] = ptTriggGen;
1800 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1801 fHJetSpecGen->Fill(tmpArrayFive);
1803 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1804 tmpArrayFive[0] = centValue;
1805 tmpArrayFive[1] = areaJetGen;
1806 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1807 tmpArrayFive[3] = ptTriggGen;
1808 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1809 fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1811 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1812 tmpArrayFive[0] = centValue;
1813 tmpArrayFive[1] = areaJetGen;
1814 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1815 tmpArrayFive[3] = ptTriggGen;
1816 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1817 fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
1819 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1820 tmpArrayFive[0] = areaJetGen;
1821 tmpArrayFive[1] = jet->Pt();
1822 tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1823 tmpArrayFive[3] = ptUeFromCellMedianGen;
1824 tmpArrayFive[4] = rhoFromCellMedianGen;
1825 fHJetUeMedianGen->Fill(tmpArrayFive);
1827 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
1828 tmpArrayFive[0] = areaJetGen;
1829 tmpArrayFive[1] = jet->Pt();
1830 tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1831 tmpArrayFive[3] = ptUeConeGen;
1832 tmpArrayFive[4] = rhoConeGen;
1833 fHJetUeConeGen->Fill(tmpArrayFive);
1836 Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1837 fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
1838 fillOnceGen = kFALSE;
1840 }//back to back jet-trigger pair
1844 if(fFillRespMx && !fIsKine){
1845 //================ RESPONSE MATRIX ===============
1846 //Count jets and trigger-jet pairs at MC generator level
1847 for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1848 AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1850 Double_t etaJetGen = jet->Eta();
1851 Double_t ptJetGen = jet->Pt();
1853 if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1854 fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
1856 Double_t areaJetGen = jet->EffectiveAreaCharged();
1857 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1858 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1860 fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
1861 fhJetPtSubUeConeGen->Fill(ptJetGen - ptUeConeGen);
1864 if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1866 Int_t ng = (Int_t) fListJetsGen->GetEntries();
1867 Int_t nr = (Int_t) fListJets->GetEntries();
1869 //Find closest MC generator - reconstructed jets
1870 if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
1871 if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
1874 Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
1875 Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
1877 //matching of MC genrator level and reconstructed jets
1878 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
1880 // Fill response matrix
1881 for(Int_t ir = 0; ir < nr; ir++){
1882 AliAODJet *recJet = (AliAODJet*) fListJets->At(ir);
1883 Double_t etaJetRec = recJet->Eta();
1884 Double_t ptJetRec = recJet->Pt();
1885 Double_t areaJetRec = recJet->EffectiveAreaCharged();
1886 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1888 if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
1889 Int_t ig = faGenIndex[ir]; //associated generator level jet
1890 if(ig >= 0 && ig < ng){
1891 if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1892 AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig);
1893 Double_t ptJetGen = genJet->Pt();
1894 Double_t etaJetGen = genJet->Eta();
1896 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1897 if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1898 fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
1900 fhJetPtResolutionVsPtGen->Fill(ptJetGen,(ptJetRec-ptJetGen)/ptJetGen);
1902 Double_t areaJetGen = genJet->EffectiveAreaCharged();
1903 Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1904 Double_t ptUeConeGen = rhoConeGen*areaJetGen;
1905 Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
1906 Double_t ptUeConeRec = rhoCone*areaJetRec;
1907 fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec,
1908 ptJetGen-ptUeFromCellMedianGen);
1909 fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
1911 if((ptJetGen-ptUeConeGen)>0){
1912 fhJetPtResolutionVsPtConeGen->Fill(ptJetGen-ptUeConeGen,((ptJetRec-ptUeConeRec)- (ptJetGen-ptUeConeGen))/ (ptJetGen-ptUeConeGen));
1916 }//rec jet in eta acceptance
1917 }//loop over reconstructed jets
1918 }// # of rec jets >0
1920 //=========================== Full jet vs charged jet matrix ==========================
1922 //Count full jets and charged-jet pairs at MC generator level
1923 for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1924 AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1925 if(!jetFull) continue;
1926 Double_t etaJetFull = jetFull->Eta();
1927 Double_t ptJetFull = jetFull->Pt();
1929 if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1930 fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
1933 if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1934 Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1935 Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1937 //Find closest MC generator full - charged jet
1938 if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1939 if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1942 Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
1943 Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
1945 //matching of MC genrator level and reconstructed jets
1946 AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1948 // Fill response matrix
1949 for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1950 AliAODJet *chJet = (AliAODJet*) fListJetsGen->At(ichr);
1951 Double_t etaJetCh = chJet->Eta();
1952 Double_t ptJetCh = chJet->Pt();
1953 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1955 if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1956 Int_t iful = faGenIndex[ichr]; //associated generator level jet
1957 if(iful >= 0 && iful < nful){
1958 if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1959 AliAODJet *genJetFull = (AliAODJet*) fListJetsGenFull->At(iful);
1960 Double_t ptJetFull = genJetFull->Pt();
1961 Double_t etaJetFull = genJetFull->Eta();
1963 //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1964 if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1965 fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1968 }//rec jet in eta acceptance
1969 }//loop over reconstructed jets
1970 }// # of rec jets >0
1971 }//pointer MC generator jets
1972 } //fill resp mx only for bin
1973 }//analyze generator level MC
1976 if(fIsKine){ //skip reconstructed data analysis in case of kine
1977 PostData(1, fOutputList);
1980 //============= RECONSTRUCTED INCLUSIVE JETS ===============
1982 Double_t etaJet = 0.0;
1983 Double_t pTJet = 0.0;
1984 Double_t areaJet = 0.0;
1985 Double_t phiJet = 0.0;
1986 //Int_t indexLeadingJet = -1;
1987 //Double_t pTLeadingJet = -10.0;
1988 //Double_t areaLeadingJet = -10.0;
1990 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1991 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1993 etaJet = jet->Eta();
1994 phiJet = jet->Phi();
1996 if(pTJet==0) continue;
1998 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1999 /*areaJet = jet->EffectiveAreaCharged();*/
2001 //Jet Diagnostics---------------------------------
2002 fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
2003 fhJetEta->Fill(pTJet, etaJet);
2004 //search for leading jet
2005 /*if(pTJet > pTLeadingJet){
2006 indexLeadingJet = ij;
2007 pTLeadingJet = pTJet;
2008 areaLeadingJet = areaJet;
2011 // raw spectra of INCLUSIVE jets
2012 //Centrality, pTjet, A
2013 Double_t fillraw[] = { centValue,
2017 fHJetPtRaw->Fill(fillraw);*/
2020 if(indexLeadingJet > -1){
2021 // raw spectra of LEADING jets
2022 //Centrality, pTjet, A
2023 Double_t fillleading[] = { centValue,
2027 fHLeadingJetPtRaw->Fill(fillleading);
2030 //================ TWO PARTICLE CORRELATION CORRECTION ===============
2031 if(fIsChargedMC && fHardest==0 ){
2032 if(indexTrigg>=0){ // Reconstructed trigger
2033 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
2034 Int_t genTriggIndex = -1;
2036 Int_t trigLabel = TMath::Abs(triggerHadron->GetLabel());
2038 Int_t nGen = particleListGen.GetEntries();
2039 for(Int_t ig=0; ig<nGen; ig++){
2040 AliVParticle *trkGen = (AliVParticle*) particleListGen.At(ig);
2041 if(!trkGen) continue;
2043 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
2044 if(trigLabel==genLabel){
2045 if(fTriggerPtRangeLow <= trkGen->Pt() &&
2046 trkGen->Pt() < fTriggerPtRangeHigh &&
2047 TMath::Abs((Float_t) trkGen->Eta()) < fTriggerEtaCut){
2048 genTriggIndex = ig; //There is a gen level particle track matching to rec trigger
2053 if(genTriggIndex>-1){ // make generator level TT x Assoc
2054 AliVParticle *trkGenTT = (AliVParticle*) particleListGen.At(genTriggIndex);
2056 fhTriggerCounterGenLevel->Fill(0.5); // count gen triggers that have reconstructed analog
2057 Double_t deltaPhi,deltaEta,deltaR;
2059 for(Int_t ig=0; ig<nGen; ig++){
2060 if(ig==genTriggIndex) continue; //skip trigger
2061 AliVParticle *trkGenA = (AliVParticle*) particleListGen.At(ig);
2062 if(!trkGenA) continue;
2063 if(trkGenA->Pt()<15.0) continue;
2064 if(trkGenA->Pt()>50.0) continue;
2065 deltaPhi = RelativePhi(trkGenTT->Phi(), trkGenA->Phi());
2066 deltaEta = trkGenTT->Eta() - trkGenA->Eta();
2067 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
2069 if(trkGenA->Pt()<20.0){ // gen level assoc particles
2070 fhDeltaRMultTriggersGenLevelLow->Fill(deltaR);
2071 fhDeltaPhiMultTriggersGenLevelLow->Fill(deltaPhi);
2073 fhDeltaRMultTriggersGenLevelHigh->Fill(deltaR);
2074 fhDeltaPhiMultTriggersGenLevelHigh->Fill(deltaPhi);
2083 // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================
2084 if(fIsChargedMC && fFillRespMx){
2085 FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
2087 Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
2088 //set ranges of the trigger loop
2089 Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
2090 Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
2092 for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
2093 indexTrigg = itrk; //trigger hadron with pT >10 GeV
2095 if(indexTrigg < 0) continue;
2097 AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
2099 PostData(1, fOutputList);
2103 if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
2105 if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
2107 //Fill trigger histograms
2108 fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
2109 fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
2110 fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
2113 //---------- make trigger-jet pairs ---------
2117 for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
2118 AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
2120 etaJet = jet->Eta();
2121 phiJet = jet->Phi();
2123 if(pTJet==0) continue;
2125 if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
2126 areaJet = jet->EffectiveAreaCharged();
2128 //subtract bg using estinates base on median of kT jets
2129 Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
2130 Double_t ptUeCone = rhoCone*areaJet;
2132 //if(areaJet >= 0.07) injet++;
2133 //if(areaJet >= 0.4) injet4++;
2135 Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
2136 fhDphiTriggerJet->Fill(dphi); //Input
2138 //Dphi versus jet pT
2139 //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
2140 /*Double_t filldp[] = { centValue,
2145 fHDphiVsJetPtAll->Fill(filldp);
2147 //Rongrong's analysis
2149 Double_t dPhi = phiJet - triggerHadron->Phi();
2150 if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
2151 if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
2152 if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
2153 if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
2155 //[A,pTjet,pTtrig,dphi]
2156 tmpArrayFour[0] = areaJet;
2157 tmpArrayFour[1] = pTJet;
2158 tmpArrayFour[2] = triggerHadron->Pt();
2159 tmpArrayFour[3] = dPhi;
2161 fHJetPhiCorr->Fill(tmpArrayFour); // Dphi distribution jet-triger
2163 //[A,pTjet-pTUe,pTtrig,dphi]
2164 tmpArrayFour[1] = pTJet - ptUeFromCellMedian;
2165 fHJetPhiCorrSubUeMedian->Fill(tmpArrayFour);
2167 //[A,pTjet-pTUe,pTtrig,dphi]
2168 tmpArrayFour[1] = pTJet - ptUeCone;
2169 fHJetPhiCorrSubUeCone->Fill(tmpArrayFour);
2171 //Leticia's analysis
2173 // Select back to back trigger - jet pairs
2174 if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
2175 fhDphiTriggerJetAccept->Fill(dphi); //Accepted
2178 //Centrality, A, pTjet, pTtrigg, dphi
2179 tmpArrayFive[0] = centValue;
2180 tmpArrayFive[1] = areaJet;
2181 tmpArrayFive[2] = pTJet;
2182 tmpArrayFive[3] = triggerHadron->Pt();
2183 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
2184 fHJetSpec->Fill(tmpArrayFive);
2186 //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
2187 tmpArrayFive[0] = centValue;
2188 tmpArrayFive[1] = areaJet;
2189 tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
2190 tmpArrayFive[3] = triggerHadron->Pt();
2191 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
2192 fHJetSpecSubUeMedian->Fill(tmpArrayFive);
2194 //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
2195 tmpArrayFive[0] = centValue;
2196 tmpArrayFive[1] = areaJet;
2197 tmpArrayFive[2] = pTJet - ptUeCone;
2198 tmpArrayFive[3] = triggerHadron->Pt();
2199 tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
2200 fHJetSpecSubUeCone->Fill(tmpArrayFive);
2202 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
2203 tmpArrayFive[0] = areaJet;
2204 tmpArrayFive[1] = pTJet;
2205 tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
2206 tmpArrayFive[3] = ptUeFromCellMedian;
2207 tmpArrayFive[4] = rhoFromCellMedian;
2208 fHJetUeMedian->Fill(tmpArrayFive);
2210 //Ue diagnostics "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
2211 tmpArrayFive[0] = areaJet;
2212 tmpArrayFive[1] = pTJet;
2213 tmpArrayFive[2] = pTJet - ptUeCone;
2214 tmpArrayFive[3] = ptUeCone;
2215 tmpArrayFive[4] = rhoCone;
2216 fHJetUeCone->Fill(tmpArrayFive);
2218 if(filledOnce){ //fill for each event only once
2219 Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
2220 fHRhoUeMedianVsCone->Fill(fillRho);
2221 filledOnce = kFALSE;
2225 //Fill Jet Density In the Event A>0.07
2227 Double_t filldens[]={ (Double_t) particleList.GetEntries(),
2231 fHJetDensity->Fill(filldens);
2234 //Fill Jet Density In the Event A>0.4
2236 Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
2237 injet4/fkAcceptance,
2240 fHJetDensityA4->Fill(filldens4);
2245 PostData(1, fOutputList);
2248 //----------------------------------------------------------------------------
2249 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
2251 // Draw result to the screen
2252 // Called once at the end of the query
2254 if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
2255 if(!GetOutputData(1)) return;
2259 //----------------------------------------------------------------------------
2260 Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
2261 //Fill the list of accepted tracks (passed track cut)
2262 //return consecutive index of the hardest ch hadron in the list
2264 AliAODEvent *aodevt = NULL;
2266 if(!fESD) aodevt = fAODIn;
2267 else aodevt = fAODOut;
2269 if(!aodevt) return -1;
2271 Int_t index = -1; //index of the highest particle in the list
2272 Double_t ptmax = -10;
2273 Int_t triggers[200];
2274 Int_t ntriggers = 0; //index in triggers array
2277 for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
2278 //AliAODTrack *tr = aodevt->GetTrack(it);
2279 AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodevt->GetTrack(it));
2280 if(!tr) AliFatal("Not a standard AOD");
2282 if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
2283 //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
2284 if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
2285 if(tr->Pt() < fTrackLowPtCut) continue;
2289 if(fHardest>0){ //leading particle
2290 if(tr->Pt()>ptmax &&
2291 TMath::Abs((Float_t) tr->Eta()) < fTriggerEtaCut){
2297 if(fHardest==0 && ntriggers<200){ //single inclusive
2298 if(fTriggerPtRangeLow <= tr->Pt() &&
2299 tr->Pt() < fTriggerPtRangeHigh &&
2300 TMath::Abs((Float_t) tr->Eta()) < fTriggerEtaCut){
2301 triggers[ntriggers] = iCount;
2309 if(ntriggers>0){ //select random inclusive trigger
2310 fhInclTrigCounter->Fill(0.5,ntriggers); //count inclusive triggers
2313 for(Int_t it1=0; it1<ntriggers; it1++){
2314 AliVParticle* tGent1 = (AliVParticle*) list->At(triggers[it1]);
2315 if(!tGent1) continue;
2316 for(Int_t ia=0; ia<list->GetEntries(); ia++){
2317 if(ia == triggers[it1]) continue;
2318 AliVParticle* tGent2 = (AliVParticle*) list->At(ia);
2319 if(!tGent2) continue;
2320 if(tGent2->Pt()<15.0) continue;
2321 if(tGent2->Pt()>20.0) continue;
2322 fhDeltaPhiMultTriggersInclLow->Fill(RelativePhi(tGent1->Phi(),tGent2->Phi()));
2326 for(Int_t it1=0; it1<ntriggers-1; it1++){
2327 AliVParticle* tGent1 = (AliVParticle*) list->At(triggers[it1]);
2328 if(!tGent1) continue;
2330 for(Int_t it2=it1+1; it2<ntriggers; it2++){
2331 AliVParticle* tGent2 = (AliVParticle*) list->At(triggers[it2]);
2332 if(!tGent2) continue;
2333 fhDeltaPhiMultTriggersInclHigh->Fill(RelativePhi(tGent1->Phi(),tGent2->Phi()));
2338 //single inclusive trigger
2339 if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
2340 Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
2341 index = triggers[rnd];
2343 fhNofMultipleTriggers->Fill(ntriggers-1);
2345 Double_t deltaPhi, deltaEta, deltaR;
2348 //Correlation with single inclusive trigger
2349 AliVParticle* tGent1 = (AliVParticle*) list->At(index);
2351 //for(Int_t ia=0; ia<ntriggers; ia++)
2352 for(Int_t ia=0; ia<list->GetEntries(); ia++){
2353 //if(triggers[ia]==index) continue;
2354 if(ia==index) continue;
2355 //AliVParticle* tGent2 = (AliVParticle*) list->At(triggers[ia]);
2356 AliVParticle* tGent2 = (AliVParticle*) list->At(ia);
2357 if(!tGent2) continue;
2358 if(tGent2->Pt()<15.0) continue;
2359 if(tGent2->Pt()>50.0) continue;
2360 deltaPhi = RelativePhi(tGent1->Phi(),tGent2->Phi());
2361 deltaEta = tGent1->Eta()-tGent2->Eta();
2362 deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
2364 if(tGent2->Pt()<20.0){
2365 fhDeltaRMultTriggersLow->Fill(deltaR);
2366 fhDeltaPhiMultTriggersLow->Fill(deltaPhi);
2368 fhDeltaRMultTriggersHigh->Fill(deltaR);
2369 fhDeltaPhiMultTriggersHigh->Fill(deltaPhi);
2373 if(tGent2->Pt()<20.0) iLow++;
2378 fhNofMultipleTriggersCone->Fill(iLow+iHigh);
2379 fhNofMultipleTriggersConeLow->Fill(iLow);
2380 fhNofMultipleTriggersConeHigh->Fill(iHigh);
2387 //----------------------------------------------------------------------------
2389 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
2390 //calculate sum of track pT in the cone perpendicular in phi to the jet
2391 //jetR = cone radius
2392 // jetPhi, jetEta = direction of the jet
2393 Int_t numberOfTrks = trkList->GetEntries();
2394 Double_t pTsum = 0.0;
2395 Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
2396 for(Int_t it=0; it<numberOfTrks; it++){
2397 AliVParticle *trk = (AliVParticle*) trkList->At(it);
2398 Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
2399 Double_t deta = trk->Eta()-jetEta;
2401 if( (dphi*dphi + deta*deta)< (jetR*jetR)){
2409 //----------------------------------------------------------------------------
2411 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
2412 //Get relative azimuthal angle of two particles -pi to pi
2413 if (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
2414 else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi();
2416 if (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
2417 else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi();
2419 Double_t dphi = mphi - vphi;
2420 if (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
2421 else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi();
2423 return dphi;//dphi in [-Pi, Pi]
2427 //----------------------------------------------------------------------------
2428 Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
2429 //fill track efficiency denominator
2430 if(trk->Pt() < fTrackLowPtCut) return kFALSE;
2431 if(trk->Charge()==0) return kFALSE;
2432 if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
2434 if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
2437 if(fHardest>0){ //leading particle
2438 if(ptLeading < trk->Pt() &&
2439 TMath::Abs((Float_t) trk->Eta()) < fTriggerEtaCut){
2441 ptLeading = trk->Pt();
2445 if(fHardest==0){ //single inclusive
2447 if(fTriggerPtRangeLow <= trk->Pt() &&
2448 trk->Pt() < fTriggerPtRangeHigh &&
2449 TMath::Abs((Float_t) trk->Eta()) < fTriggerEtaCut){
2457 //----------------------------------------------------------------------------
2458 void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
2459 //fill track efficiency numerator
2460 if(!recList) return;
2461 if(!genList) return;
2462 Int_t nRec = recList->GetEntries();
2463 Int_t nGen = genList->GetEntries();
2464 for(Int_t ir=0; ir<nRec; ir++){
2465 AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
2466 if(!trkRec) continue;
2467 Int_t recLabel = TMath::Abs(trkRec->GetLabel());
2468 Bool_t matched = kFALSE;
2470 for(Int_t ig=0; ig<nGen; ig++){
2471 AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
2472 if(!trkGen) continue;
2473 Int_t genLabel = TMath::Abs(trkGen->GetLabel());
2474 if(recLabel==genLabel){
2475 fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
2480 if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
2485 //________________________________________________________________
2486 void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart, Double_t &rhoMedian, Int_t mode){
2487 //Estimate background rho by means of integrating track pT outside identified jet cones
2491 //identify leading jet in the full track acceptance
2492 Int_t idxLeadingJet = -1;
2493 Double_t pTleading = 0.0;
2494 AliAODJet* jet = NULL;
2496 for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
2497 jet = (AliAODJet*)(listJet->At(ij));
2499 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
2500 if(jet->Pt() > pTleading){
2502 pTleading = jet->Pt();
2507 static Double_t jphi[100];
2508 static Double_t jeta[100];
2509 static Double_t jRsquared[100];
2511 for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
2513 jet = (AliAODJet*)(listJet->At(ij));
2515 if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
2517 if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
2518 jphi[njets] = RelativePhi(jet->Phi(),0.0); //-pi,pi
2519 jeta[njets] = jet->Eta();
2520 jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
2526 static Double_t nOutCone[10][4];
2527 static Double_t sumPtOutOfCone[10][4];
2529 for(Int_t ie=0; ie<fnEta; ie++){
2530 for(Int_t ip=0; ip<fnPhi; ip++){
2531 nOutCone[ip][ie] = 0.0; //initialize counter
2532 sumPtOutOfCone[ip][ie] = 0.0;
2536 Double_t rndphi, rndeta;
2537 Double_t rndphishift, rndetashift;
2538 Double_t dphi, deta;
2541 for(Int_t it=0; it<fnTrials; it++){
2543 rndphi = fRandom->Uniform(0, fPhiSize);
2544 rndeta = fRandom->Uniform(0, fEtaSize);
2546 for(Int_t ip=0; ip<fnPhi; ip++){ //move radom position to each cell
2547 rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
2548 for(Int_t ie=0; ie<fnEta; ie++){
2549 rndetashift = rndeta + ie*fEtaSize - fEtaSize;
2551 bIsInCone = 0; //tag if trial is in the jet cone
2552 for(Int_t ij=0; ij<njets; ij++){
2553 deta = jeta[ij] - rndetashift;
2554 dphi = RelativePhi(rndphishift,jphi[ij]);
2555 if((dphi*dphi + deta*deta) < jRsquared[ij]){
2560 if(!bIsInCone) nOutCone[ip][ie]++;
2566 //Sum particle pT outside jets in cells
2567 Int_t npart = listPart->GetEntries();
2568 Int_t phicell,etacell;
2569 AliVParticle *part = NULL;
2570 for(Int_t ip=0; ip < npart; ip++){
2572 part = (AliVParticle*) listPart->At(ip);
2574 if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
2576 bIsInCone = 0; //init
2577 for(Int_t ij=0; ij<njets; ij++){
2578 dphi = RelativePhi(jphi[ij], part->Phi());
2579 deta = jeta[ij] - part->Eta();
2580 if((dphi*dphi + deta*deta) < jRsquared[ij]){
2586 phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
2587 etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
2588 sumPtOutOfCone[phicell][etacell]+= part->Pt();
2592 static Double_t rhoInCells[20];
2593 Double_t relativeArea;
2595 Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
2596 for(Int_t ip=0; ip<fnPhi; ip++){
2597 for(Int_t ie=0; ie<fnEta; ie++){
2598 relativeArea = nOutCone[ip][ie]/fnTrials;
2599 //cout<<"RA= "<< relativeArea<<" BA="<<bufferArea<<" BPT="<< bufferPt <<endl;
2601 bufferArea += relativeArea;
2602 bufferPt += sumPtOutOfCone[ip][ie];
2603 if(bufferArea > fJetFreeAreaFrac){
2604 rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
2606 if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
2607 else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
2618 rhoMedian = TMath::Median(nCells, rhoInCells);
2619 if(mode==1){ //mc data
2620 fhEntriesToMedianGen->Fill(nCells);
2622 fhEntriesToMedian->Fill(nCells);
2628 //__________________________________________________________________
2629 void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart, Double_t &rhoPerpCone){
2631 rhoPerpCone = 0.0; //init
2633 if(!listJet) return;
2634 Int_t njet = listJet->GetEntries();
2635 Int_t npart = listPart->GetEntries();
2636 Double_t pTleading = 0.0;
2637 Double_t phiLeading = 1000.;
2638 Double_t etaLeading = 1000.;
2640 for(Int_t jsig=0; jsig < njet; jsig++){
2641 AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
2642 if(!signaljet) continue;
2643 if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
2644 if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
2645 pTleading = signaljet->Pt();
2646 phiLeading = signaljet->Phi();
2647 etaLeading = signaljet->Eta();
2651 Double_t phileftcone = phiLeading + TMath::Pi()/2;
2652 Double_t phirightcone = phiLeading - TMath::Pi()/2;
2655 for(Int_t ip=0; ip < npart; ip++){
2657 AliVParticle *part = (AliVParticle*) listPart->At(ip);
2662 dp = RelativePhi(phileftcone, part->Phi());
2663 de = etaLeading - part->Eta();
2664 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2666 dp = RelativePhi(phirightcone, part->Phi());
2667 if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2671 //normalize total pT by two times cone are
2672 rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
2675 //__________________________________________________________________
2677 void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
2678 //Convert TClones array of jets to TList
2680 if(!strlen(bname.Data())){
2681 AliError(Form("Jet branch %s not set.", bname.Data()));
2685 TClonesArray *array=0x0;
2686 if(fUseExchContainer){ //take input from exchange containers
2687 if(fJetBranchNameKine.Length()>0 && bname.CompareTo(fJetBranchNameKine.Data())==0){
2688 array = dynamic_cast<TClonesArray*>(GetInputData(1)); //connect exchange container slot 1
2690 if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
2691 array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
2693 }else{ //take input from AOD
2694 if(fAODOut&&!array){
2695 array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
2697 if(fAODExtension&&!array){
2698 array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
2701 array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
2705 list->Clear(); //clear list beforehand
2709 Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
2712 for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
2713 AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
2714 if (jet) list->Add(jet);