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 fhDeltaRMultTriggersLow(0x0),
 146 fhDeltaRMultTriggersHigh(0x0),
 147 fhDeltaPhiMultTriggersLow(0x0),
 148 fhDeltaPhiMultTriggersHigh(0x0),
 149 fhDeltaPhiMultTriggersInclLow(0x0),
 150 fhDeltaPhiMultTriggersInclHigh(0x0),
 151 fhInclTrigCounter(0x0),
 152 //fHJetPtRaw(0x0),
 153 //fHLeadingJetPtRaw(0x0),
 154 //fHDphiVsJetPtAll(0x0),
 155 fhJetPtGenVsJetPtRec(0x0),
 156 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
 157 fhJetPtGenVsJetPtRecSubUeCone(0x0),
 158 fhJetPtGen(0x0),
 159 fhJetPtSubUeMedianGen(0x0),
 160 fhJetPtSubUeConeGen(0x0),
 161 fhJetPtResolutionVsPtGen(0x0),
 162 fhJetPtResolutionVsPtConeGen(0x0),
 163 fhJetPtGenChargVsJetPtGenFull(0x0),
 164 fhJetPtGenFull(0x0),
 165 fh2NtriggersGen(0x0),
 166 fHJetSpecGen(0x0),
 167 fHJetSpecSubUeMedianGen(0x0),
 168 fHJetSpecSubUeConeGen(0x0),
 169 fHJetPhiCorrGen(0x0),
 170 fHJetPhiCorrSubUeMedianGen(0x0),
 171 fHJetPhiCorrSubUeConeGen(0x0),
 172 fHJetUeMedianGen(0x0),
 173 fHJetUeConeGen(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),
 214 fIsChargedMC(0),
 215 fIsKine(0),
 216 fIsFullMC(0),
 217 faGenIndex(0),
 218 faRecIndex(0),
 219 fkAcceptance(2.0*TMath::Pi()*1.8),
 220 fkDeltaPhiCut(TMath::Pi()-0.8),
 221 fh1Xsec(0x0),
 222 fh1Trials(0x0),
 223 fh1AvgTrials(0x0),
 224 fh1PtHard(0x0),
 225 fh1PtHardNoW(0x0),
 226 fh1PtHardTrials(0x0),
 227 fAvgTrials(1),
 228 fHardest(0),
 229 fEventNumberRangeLow(0),
 230 fEventNumberRangeHigh(99),
 231 fTriggerPtRangeLow(0.0),
 232 fTriggerPtRangeHigh(10000.0),
 233 fFillRespMx(0),
 234 fRandom(0x0),
 235 fnTrials(1000),
 236 fJetFreeAreaFrac(0.5),
 237 fnPhi(9),
 238 fnEta(2),
 239 fEtaSize(0.7),
 240 fPhiSize(2*TMath::Pi()/fnPhi),
 241 fCellArea(fPhiSize*fEtaSize),
 242 fSafetyMargin(1.1),
 243 fDoubleBinning(kFALSE)
 244 {
 245    // default Constructor
 246 }
 247
 248 //---------------------------------------------------------------------
 249
 250 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
 251 AliAnalysisTaskSE(name),
 252 fESD(0x0),
 253 fAODIn(0x0),
 254 fAODOut(0x0),
 255 fAODExtension(0x0),
 256 fMcEvent(0x0),
 257 fMcHandler(0x0),
 258 fJetBranchName(""),
 259 fJetBranchNameChargMC(""),
 260 fJetBranchNameKine(""),
 261 fJetBranchNameFullMC(""),
 262 fJetBranchNameBg(""),
 263 fJetBranchNameBgChargMC(""),
 264 fJetBranchNameBgKine(""),
 265 fListJets(0x0),
 266 fListJetsGen(0x0),
 267 fListJetsGenFull(0x0),
 268 fListJetsBg(0x0),
 269 fListJetsBgGen(0x0),
 270 fNonStdFile(""),
 271 fSystem(0),  //pp=0   pPb=1
 272 fJetParamR(0.4),
 273 fBgJetParamR(0.3),
 274 fBgMaxJetPt(8.0),
 275 fBgConeR(0.4),
 276 fOfflineTrgMask(AliVEvent::kAny),
 277 fMinContribVtx(1),
 278 fVtxZMin(-10.0),
 279 fVtxZMax(10.0),
 280 fFilterMask(0),
 281 fCentMin(0.0),
 282 fCentMax(100.0),
 283 fJetEtaMin(-0.5),
 284 fJetEtaMax(0.5),
 285 fTriggerEtaCut(0.9),
 286 fTrackEtaCut(0.9),
 287 fTrackLowPtCut(0.15),
 288 fUseExchContainer(0),
 289 fOutputList(0x0),
 290 fHistEvtSelection(0x0),
 291 fh2Ntriggers(0x0),
 292 fHJetSpec(0x0),
 293 fHJetSpecSubUeMedian(0x0),
 294 fHJetSpecSubUeCone(0x0),
 295 fHJetPhiCorr(0x0),
 296 fHJetPhiCorrSubUeMedian(0x0),
 297 fHJetPhiCorrSubUeCone(0x0),
 298 fHJetUeMedian(0x0),
 299 fHJetUeCone(0x0),
 300 fHRhoUeMedianVsCone(0x0),
 301 //fHJetDensity(0x0),
 302 //fHJetDensityA4(0x0),
 303 fhJetPhi(0x0),
 304 fhTriggerPhi(0x0),
 305 fhJetEta(0x0),
 306 fhTriggerEta(0x0),
 307 fhVertexZ(0x0),
 308 fhVertexZAccept(0x0),
 309 fhContribVtx(0x0),
 310 fhContribVtxAccept(0x0),
 311 fhDphiTriggerJet(0x0),
 312 fhDphiTriggerJetAccept(0x0),
 313 fhCentrality(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),
 326 //fHJetPtRaw(0x0),
 327 //fHLeadingJetPtRaw(0x0),
 328 //fHDphiVsJetPtAll(0x0),
 329 fhJetPtGenVsJetPtRec(0x0),
 330 fhJetPtGenVsJetPtRecSubUeMedian(0x0),
 331 fhJetPtGenVsJetPtRecSubUeCone(0x0),
 332 fhJetPtGen(0x0),
 333 fhJetPtSubUeMedianGen(0x0),
 334 fhJetPtSubUeConeGen(0x0),
 335 fhJetPtResolutionVsPtGen(0x0),
 336 fhJetPtResolutionVsPtConeGen(0x0),
 337 fhJetPtGenChargVsJetPtGenFull(0x0),
 338 fhJetPtGenFull(0x0),
 339 fh2NtriggersGen(0x0),
 340 fHJetSpecGen(0x0),
 341 fHJetSpecSubUeMedianGen(0x0),
 342 fHJetSpecSubUeConeGen(0x0),
 343 fHJetPhiCorrGen(0x0),
 344 fHJetPhiCorrSubUeMedianGen(0x0),
 345 fHJetPhiCorrSubUeConeGen(0x0),
 346 fHJetUeMedianGen(0x0),
 347 fHJetUeConeGen(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),
 388 fIsChargedMC(0),
 389 fIsKine(0),
 390 fIsFullMC(0),
 391 faGenIndex(0),
 392 faRecIndex(0),
 393 fkAcceptance(2.0*TMath::Pi()*1.8),
 394 fkDeltaPhiCut(TMath::Pi()-0.8),
 395 fh1Xsec(0x0),
 396 fh1Trials(0x0),
 397 fh1AvgTrials(0x0),
 398 fh1PtHard(0x0),
 399 fh1PtHardNoW(0x0),
 400 fh1PtHardTrials(0x0),
 401 fAvgTrials(1),
 402 fHardest(0),
 403 fEventNumberRangeLow(0),
 404 fEventNumberRangeHigh(99),
 405 fTriggerPtRangeLow(0.0),
 406 fTriggerPtRangeHigh(10000.0),
 407 fFillRespMx(0),
 408 fRandom(0x0),
 409 fnTrials(1000),
 410 fJetFreeAreaFrac(0.5),
 411 fnPhi(9),
 412 fnEta(2),
 413 fEtaSize(0.7),
 414 fPhiSize(2*TMath::Pi()/fnPhi),
 415 fCellArea(fPhiSize*fEtaSize),
 416 fSafetyMargin(1.1),
 417 fDoubleBinning(kFALSE)
 418 {
 419 // Constructor
 420    DefineOutput(1, TList::Class());
 421
 422    TString dummy(name);
 423    if(dummy.Contains("KINE")){
 424       DefineInput(1, TClonesArray::Class());
 425       DefineInput(2, TClonesArray::Class());
 426    }
 427 }
 428
 429 //--------------------------------------------------------------
 430 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
 431 AliAnalysisTaskSE(a.GetName()),
 432 fESD(a.fESD),
 433 fAODIn(a.fAODIn),
 434 fAODOut(a.fAODOut),
 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),
 451 fSystem(a.fSystem),
 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),
 569 fIsKine(a.fIsKine),
 570 fIsFullMC(a.fIsFullMC),
 571 faGenIndex(a.faGenIndex),
 572 faRecIndex(a.faRecIndex),
 573 fkAcceptance(a.fkAcceptance),
 574 fkDeltaPhiCut(a.fkDeltaPhiCut),
 575 fh1Xsec(a.fh1Xsec),
 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),
 588 fRandom(a.fRandom),
 589 fnTrials(a.fnTrials),
 590 fJetFreeAreaFrac(a.fJetFreeAreaFrac),
 591 fnPhi(a.fnPhi),
 592 fnEta(a.fnEta),
 593 fEtaSize(a.fEtaSize),
 594 fPhiSize(a.fPhiSize),
 595 fCellArea(a.fCellArea),
 596 fSafetyMargin(a.fSafetyMargin),
 597 fDoubleBinning(a.fDoubleBinning)
 598 {
 599    //Copy Constructor
 600    fRandom->SetSeed(0);
 601 }
 602 //--------------------------------------------------------------
 603
 604 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
 605   // assignment operator
 606   this->~AliAnalysisTaskJetCorePP();
 607   new(this) AliAnalysisTaskJetCorePP(a);
 608   return *this;
 609 }
 610 //--------------------------------------------------------------
 611
 612 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
 613 {
 614    //Destructor
 615    delete fListJets;
 616    delete fListJetsGen;
 617    delete fListJetsGenFull;
 618    delete fListJetsBg;
 619    delete fListJetsBgGen;
 620    delete fOutputList; // ????
 621    delete fRandom;
 622 }
 623
 624 //--------------------------------------------------------------
 625
 626
 627 Bool_t AliAnalysisTaskJetCorePP::Notify()
 628 {
 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;
 634    Float_t trials  = 1;
 635    fAvgTrials = 1;
 636
 637    if(fIsChargedMC){
 638       fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 639       if(!fESD){
 640          if(fDebug>1) AliError("ESD not available");
 641          fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 642       }
 643
 644       fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 645
 646
 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;
 651          if(!fAODExtension){
 652             if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
 653          }
 654       }
 655
 656       TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 657
 658       if(tree){
 659          TFile *curfile = tree->GetCurrentFile();
 660          if(!curfile) {
 661             Error("Notify","No current file");
 662             return kFALSE;
 663          }
 664          if(!fh1Xsec || !fh1Trials){
 665             Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 666             return kFALSE;
 667          }
 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);
 674       }
 675
 676       if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
 677    }
 678
 679
 680    return kTRUE;
 681 }
 682 //--------------------------------------------------------------
 683
 684 void AliAnalysisTaskJetCorePP::Init()
 685 {
 686    // check for jet branches
 687    if(fJetBranchNameKine.Length()==0){
 688       if(!strlen(fJetBranchName.Data())){
 689          AliError("Jet branch name not set.");
 690       }
 691    }
 692    fMcHandler = dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 693
 694
 695 }
 696
 697 //--------------------------------------------------------------
 698
 699 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
 700 {
 701    // Create histograms and initilize variables
 702    fSafetyMargin = fBgConeR*fBgConeR /(fBgJetParamR*fBgJetParamR);
 703
 704    // Called once
 705    fListJets   = new TList();  //reconstructed level antikt jets
 706    fListJetsBg = new TList();  //reconstructed jets to be removed from bg
 707
 708    fIsChargedMC = (fJetBranchNameChargMC.Length()>0) ? kTRUE : kFALSE;
 709    fIsKine      = (fJetBranchNameKine.Length()>0)    ? kTRUE : kFALSE;
 710    fIsFullMC    = (fJetBranchNameFullMC.Length()>0)  ? kTRUE : kFALSE;
 711
 712    fRandom = new TRandom3(0);
 713
 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
 717
 718       if(fIsFullMC)
 719          fListJetsGenFull = new TList(); //generator level full jets
 720    }
 721    OpenFile(1);
 722    if(!fOutputList) fOutputList = new TList();
 723    fOutputList->SetOwner(kTRUE);
 724
 725    Bool_t oldStatus = TH1::AddDirectoryStatus();
 726    TH1::AddDirectory(kFALSE);
 727
 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)");
 735
 736    fOutputList->Add(fHistEvtSelection);
 737
 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);
 743
 744    Int_t bw = (fDoubleBinning==0) ? 1 : 2; //make larger bin width
 745
 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);
 755
 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);
 764
 765
 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);
 774
 775    if(!fIsKine) fOutputList->Add(fHJetPhiCorr); // Dphi distribution jet-triger
 776
 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);
 785
 786
 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);
 797
 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);
 802
 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};
 808
 809    fHRhoUeMedianVsCone = new THnSparseF("hRhoUeMedianVsCone","[Rho Cone, Rho Median]",
 810                                       dimRho, nBinsRho, lowBinRho, hiBinRho);
 811    if(!fIsKine) fOutputList->Add(fHRhoUeMedianVsCone);
 812
 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};
 818
 819    fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
 820                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
 821
 822    fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
 823                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
 824
 825    fOutputList->Add(fHJetDensity);
 826    fOutputList->Add(fHJetDensityA4);
 827    */
 828
 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);
 838
 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);
 849
 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");
 854
 855    fhDeltaRMultTriggersLow = new  TH1D("fhDeltaRMultTriggersLow","fhDeltaRMultTriggersLow", 200,0,4);
 856    fhDeltaRMultTriggersHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersHigh");
 857
 858    fhDeltaPhiMultTriggersLow  = new  TH1D("fhDeltaPhiMultTriggersLow","fhDeltaPhiRultTriggers", 320,-3.2,3.2); //single incl trigger
 859    fhDeltaPhiMultTriggersHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersHigh");  //single incl trigger
 860
 861    fhDeltaPhiMultTriggersInclLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersInclLow");
 862    fhDeltaPhiMultTriggersInclHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersInclHigh");
 863    fhInclTrigCounter = new TH1D("fhInclTrigCounter","fhInclTrigCounter",1,0,1);
 864
 865    if(!fIsKine){
 866       fOutputList->Add(fhJetPhi);
 867       fOutputList->Add(fhTriggerPhi);
 868       fOutputList->Add(fhJetEta);
 869       fOutputList->Add(fhTriggerEta);
 870    }
 871    fOutputList->Add(fhVertexZ);
 872    fOutputList->Add(fhVertexZAccept);
 873    if(!fIsKine){
 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);
 893    }
 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);
 904
 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);
 911
 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);
 922    */
 923
 924    //analyze MC generator level
 925    if(fIsChargedMC || fIsKine){
 926       if(fFillRespMx){
 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
 930          //....
 931          fhJetPtGenVsJetPtRecSubUeMedian = new TH2D("fhJetPtGenVsJetPtRecSubUeMedian","fhJetPtGenVsJetPtRecSubUeMedian", bw*110,-20,200, bw*110,-20,200);
 932          fOutputList->Add(fhJetPtGenVsJetPtRecSubUeMedian); // with kT median bg subtr
 933          //....
 934          fhJetPtGenVsJetPtRecSubUeCone=(TH2D*)fhJetPtGenVsJetPtRecSubUeMedian->Clone("fhJetPtGenVsJetPtRecSubUeCone");
 935          fhJetPtGenVsJetPtRecSubUeCone->SetTitle("fhJetPtGenVsJetPtRecSubUeCone");
 936          fOutputList->Add(fhJetPtGenVsJetPtRecSubUeCone); // with weighted kT median bg subtr
 937          //....
 938          fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",bw*100,0,200); //MC generator charged jet pt spectrum
 939          fOutputList->Add(fhJetPtGen);
 940          //....
 941          fhJetPtSubUeMedianGen = new TH1D("fhJetPtSubUeMedianGen","Jet Pt - UE pT (MC Gen)",bw*110,-20,200);
 942          fOutputList->Add(fhJetPtSubUeMedianGen);  // with kT median bg subtr
 943          //....
 944          fhJetPtSubUeConeGen = (TH1D*) fhJetPtSubUeMedianGen->Clone("fhJetPtSubUeConeGen");
 945          fOutputList->Add(fhJetPtSubUeConeGen); // with weighted kT median bg subtr
 946          //....
 947          fhJetPtResolutionVsPtGen = new TH2D("fhJetPtResolutionVsPtGen","fhJetPtResolutionVsPtGen",20,0,100, 35,-1.,0.4);
 948          fOutputList->Add(fhJetPtResolutionVsPtGen); // with weighted kT median bg subtr
 949           //....
 950          fhJetPtResolutionVsPtConeGen = (TH2D*) fhJetPtResolutionVsPtGen->Clone("fhJetPtResolutionVsPtConeGen");
 951          fOutputList->Add(fhJetPtResolutionVsPtConeGen); // with weighted kT median bg subtr
 952
 953          //....
 954          if(fIsFullMC){
 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
 957          //....
 958             fhJetPtGenFull = new TH1D("fhJetPtGenFull","Jet Pt (MC Full jets Gen)",100,0,200); //MC generator full jet pt spectrum
 959             fOutputList->Add(fhJetPtGenFull);
 960          }
 961       }
 962       //....
 963       fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
 964       fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
 965       fOutputList->Add(fh2NtriggersGen);
 966
 967       //Centrality, A, pT, pTtrigg
 968       fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
 969       fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
 970       fOutputList->Add(fHJetSpecGen);
 971
 972       fHJetSpecSubUeMedianGen = (THnSparseF*)  fHJetSpecSubUeMedian->Clone("fHJetSpecSubUeMedianGen");
 973       fHJetSpecSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeMedian->GetTitle()));
 974       fOutputList->Add(fHJetSpecSubUeMedianGen);
 975
 976       fHJetSpecSubUeConeGen =  (THnSparseF*) fHJetSpecSubUeCone->Clone("fHJetSpecSubUeConeGen");
 977       fHJetSpecSubUeConeGen->SetTitle(Form("%s Gen MC",fHJetSpecSubUeCone->GetTitle()));
 978       fOutputList->Add(fHJetSpecSubUeConeGen);
 979       //---
 980       fHJetPhiCorrGen = (THnSparseF*) fHJetPhiCorr->Clone("fHJetPhiCorrGen");
 981       fHJetPhiCorrGen->SetTitle(Form("%s Gen MC",fHJetPhiCorr->GetTitle()));
 982       fOutputList->Add(fHJetPhiCorrGen);
 983
 984       fHJetPhiCorrSubUeMedianGen = (THnSparseF*) fHJetPhiCorrSubUeMedian->Clone("fHJetPhiCorrSubUeMedianGen");
 985       fHJetPhiCorrSubUeMedianGen->SetTitle(Form("%s Gen MC",fHJetPhiCorrSubUeMedian->GetTitle()));
 986       fOutputList->Add(fHJetPhiCorrSubUeMedianGen);
 987
 988       fHJetPhiCorrSubUeConeGen = (THnSparseF*) fHJetPhiCorrSubUeCone->Clone("fHJetPhiCorrSubUeConeGen");
 989       fHJetPhiCorrSubUeConeGen->SetTitle(Form("%s Gen MC", fHJetPhiCorrSubUeCone->GetTitle()));
 990       fOutputList->Add(fHJetPhiCorrSubUeConeGen);
 991
 992       //---
 993       fHJetUeMedianGen  =  (THnSparseF*) fHJetUeMedian->Clone("fHJetUeMedianGen");
 994       fHJetUeMedianGen->SetTitle(Form("%s Gen MC", fHJetUeMedian->GetTitle()));
 995       fOutputList->Add(fHJetUeMedianGen);
 996
 997       fHJetUeConeGen =  (THnSparseF*) fHJetUeCone->Clone("fHJetUeConeGen");
 998       fHJetUeConeGen->SetTitle(Form("%s Gen MC", fHJetUeCone->GetTitle()));
 999       fOutputList->Add(fHJetUeConeGen);
1000
1001       if(fFillRespMx){
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);
1007
1008          fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
1009          fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
1010          fOutputList->Add(fhPtTrkTruePrimGen);
1011
1012          fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
1013          fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
1014          fOutputList->Add(fhPtTrkSecOrFakeRec);
1015       }
1016
1017       fHRhoUeMedianVsConeGen = (THnSparseF*) fHRhoUeMedianVsCone->Clone("hRhoUeMedianVsConeGen");
1018       fHRhoUeMedianVsConeGen->SetTitle(Form("%s Gen MC", fHRhoUeMedianVsCone->GetTitle()));
1019       fOutputList->Add(fHRhoUeMedianVsConeGen);
1020
1021       fhEntriesToMedianGen = (TH1D*) fhEntriesToMedian->Clone("fhEntriesToMedianGen");
1022       fhEntriesToMedianGen->SetTitle(Form("%s Gen MC", fhEntriesToMedian->GetTitle()));
1023       fOutputList->Add(fhEntriesToMedianGen);
1024
1025       fhCellAreaToMedianGen  = (TH1D*) fhCellAreaToMedian->Clone("fhCellAreaToMedianGen");
1026       fhCellAreaToMedianGen->SetTitle(Form("%s Gen MC", fhCellAreaToMedian->GetTitle()));
1027       fOutputList->Add(fhCellAreaToMedianGen);
1028
1029       fhNofMultipleTriggersGen = (TH1D*) fhNofMultipleTriggers->Clone("fhNofMultipleTriggersGen");
1030       fOutputList->Add(fhNofMultipleTriggersGen);
1031
1032       fhNofMultipleTriggersConeGen = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGen");
1033       fOutputList->Add(fhNofMultipleTriggersConeGen);
1034
1035       fhNofMultipleTriggersConeGenLow = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGenpta15to20");
1036       fOutputList->Add(fhNofMultipleTriggersConeGenLow);
1037
1038       fhNofMultipleTriggersConeGenHigh = (TH1D*) fhNofMultipleTriggersCone->Clone("fhNofMultipleTriggersConeGenpta20to50");
1039       fOutputList->Add(fhNofMultipleTriggersConeGenHigh);
1040
1041       fhDeltaRMultTriggersGenLow  = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenLow");
1042       fOutputList->Add(fhDeltaRMultTriggersGenLow);
1043
1044       fhDeltaRMultTriggersGenHigh  = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenHigh");
1045       fOutputList->Add(fhDeltaRMultTriggersGenHigh);
1046
1047       fhDeltaPhiMultTriggersGenLow  = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenLow");
1048       fOutputList->Add(fhDeltaPhiMultTriggersGenLow);
1049
1050       fhDeltaPhiMultTriggersGenHigh  = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenHigh");
1051       fOutputList->Add(fhDeltaPhiMultTriggersGenHigh);
1052
1053       fhDeltaPhiMultTriggersInclGenLow  = (TH1D*) fhDeltaPhiMultTriggersInclLow->Clone("fhDeltaPhiMultTriggersInclGenLow");
1054       fOutputList->Add(fhDeltaPhiMultTriggersInclGenLow);
1055
1056       fhDeltaPhiMultTriggersInclGenHigh = (TH1D*) fhDeltaPhiMultTriggersInclHigh->Clone("fhDeltaPhiMultTriggersInclGenHigh");
1057       fOutputList->Add(fhDeltaPhiMultTriggersInclGenHigh);
1058
1059       fhInclTrigCounterGen= (TH1D*) fhInclTrigCounter->Clone("fhInclTrigCounterGen");
1060       fOutputList->Add(fhInclTrigCounterGen);
1061
1062       fhNofMultipleTriggersConeGenA = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10");
1063       fOutputList->Add(fhNofMultipleTriggersConeGenA);
1064
1065       fhNofMultipleTriggersConeGenALow = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10pTa15to20");
1066       fOutputList->Add(fhNofMultipleTriggersConeGenALow);
1067
1068       fhNofMultipleTriggersConeGenAHigh = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen10pTa20to50");
1069       fOutputList->Add(fhNofMultipleTriggersConeGenAHigh);
1070
1071       fhDeltaRMultTriggersGenALow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen10pTa15to20");
1072       fOutputList->Add(fhDeltaRMultTriggersGenALow);
1073
1074       fhDeltaPhiMultTriggersGenALow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen10pTa15to20");
1075       fOutputList->Add(fhDeltaPhiMultTriggersGenALow);
1076
1077       fhDeltaRMultTriggersGenAHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen10pTa20to50");
1078       fOutputList->Add(fhDeltaRMultTriggersGenAHigh);
1079
1080       fhDeltaPhiMultTriggersGenAHigh =  (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen10pTa20to50");
1081       fOutputList->Add(fhDeltaPhiMultTriggersGenAHigh);
1082
1083       fhNofTriggersGenA = (TH1D*) fhInclTrigCounter->Clone("fhNofTriggersGen10");
1084       fOutputList->Add(fhNofTriggersGenA);
1085
1086
1087       fhNofMultipleTriggersConeGenB = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5");
1088       fOutputList->Add(fhNofMultipleTriggersConeGenB);
1089
1090       fhNofMultipleTriggersConeGenBLow = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5pTa15to20");
1091       fOutputList->Add(fhNofMultipleTriggersConeGenBLow);
1092
1093       fhNofMultipleTriggersConeGenBHigh = (TH1D*) fhNofMultipleTriggersConeGen->Clone("fhNofMultipleTriggersConeGen5pTa20to50");
1094       fOutputList->Add(fhNofMultipleTriggersConeGenBHigh);
1095
1096       fhDeltaRMultTriggersGenBLow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen5pTa15to20");
1097       fOutputList->Add(fhDeltaRMultTriggersGenBLow);
1098
1099       fhDeltaPhiMultTriggersGenBLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen5pTa15to20");
1100       fOutputList->Add(fhDeltaPhiMultTriggersGenBLow);
1101
1102       fhDeltaRMultTriggersGenBHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGen5pTa20to50");
1103       fOutputList->Add(fhDeltaRMultTriggersGenBHigh);
1104
1105       fhDeltaPhiMultTriggersGenBHigh =  (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGen5pTa20to50");
1106       fOutputList->Add(fhDeltaPhiMultTriggersGenBHigh);
1107
1108       fhNofTriggersGenB = (TH1D*) fhInclTrigCounter->Clone("fhNofTriggersGen5");
1109       fOutputList->Add(fhNofTriggersGenB);
1110
1111
1112
1113
1114       if(fIsChargedMC){  // Filled from gen level if  well reconstructed trigger exists
1115          fhTriggerCounterGenLevel = (TH1D*) fhInclTrigCounter->Clone("fhTriggerCounterGenLevel");
1116          fOutputList->Add(fhTriggerCounterGenLevel);
1117
1118          fhDeltaRMultTriggersGenLevelLow = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenLevelLow");
1119          fOutputList->Add(fhDeltaRMultTriggersGenLevelLow);
1120
1121          fhDeltaPhiMultTriggersGenLevelLow = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenLevelLow");
1122          fOutputList->Add(fhDeltaPhiMultTriggersGenLevelLow);
1123
1124          fhDeltaRMultTriggersGenLevelHigh = (TH1D*) fhDeltaRMultTriggersLow->Clone("fhDeltaRMultTriggersGenLevelHigh");
1125          fOutputList->Add(fhDeltaRMultTriggersGenLevelHigh);
1126
1127          fhDeltaPhiMultTriggersGenLevelHigh = (TH1D*) fhDeltaPhiMultTriggersLow->Clone("fhDeltaPhiMultTriggersGenLevelHigh");
1128          fOutputList->Add(fhDeltaPhiMultTriggersGenLevelHigh);
1129       }
1130    }
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++){
1136       if(iPt == 0){
1137          binLimitsPt[iPt] = -50.0;
1138       }else{// 1.0
1139          binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
1140       }
1141    }
1142
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);
1158
1159
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));
1163       if(h1){
1164          h1->Sumw2();
1165          continue;
1166       }
1167       THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
1168       if(hn){
1169          hn->Sumw2();
1170       }
1171    }
1172    TH1::AddDirectory(oldStatus);
1173
1174    PostData(1, fOutputList);
1175 }
1176
1177 //--------------------------------------------------------------------
1178
1179 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
1180 {
1181    //User Exec
1182
1183
1184    //Event loop
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);
1189
1190    if(TMath::Abs((Float_t) fJetParamR)<0.00001){
1191       AliError("ANTIKT Cone radius is set to zero.");
1192       return;
1193    }
1194
1195    if(TMath::Abs((Float_t) fBgJetParamR)<0.00001){
1196       AliError("ANTIKT Cone radius is set to zero.");
1197       return;
1198    }
1199
1200    if(!fIsKine){ //real data or full MC
1201       if(!strlen(fJetBranchName.Data())){
1202          AliError("Name of jet branch not set.");
1203          return;
1204       }
1205       if(!strlen(fJetBranchNameBg.Data())){
1206          AliError("Name of jet bg branch not set.");
1207          return;
1208       }
1209    }else{  //Kine
1210        if(!strlen(fJetBranchNameBgKine.Data())){
1211          AliError("Name of jet bg branch for kine not set.");
1212          return;
1213        }
1214
1215       Init();
1216       if(fMcHandler){
1217          fMcEvent = fMcHandler->MCEvent();
1218       }else{
1219          if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcHandler=NULL\n");
1220          PostData(1, fOutputList);
1221          return;
1222       }
1223       if(!fMcEvent){
1224          if(fDebug > 1) printf("AliAnalysisTaskJetCorePP::Exec() fMcEvent=NULL \n");
1225          PostData(1, fOutputList);
1226          return;
1227       }
1228
1229       Float_t xsection = 0;
1230       Float_t trials  = 0;
1231
1232       AliGenPythiaEventHeader *genPH =
1233          dynamic_cast<AliGenPythiaEventHeader*> (fMcEvent->GenEventHeader());
1234       if(genPH){
1235          xsection = genPH->GetXsection();
1236          trials   = genPH->Trials();
1237          ptHard   = genPH->GetPtHard();
1238       }
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);
1244    }
1245
1246
1247    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
1248    if(!fESD){
1249       if(fDebug>1) AliError("ESD not available");
1250       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
1251    }
1252
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
1259    }
1260
1261
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;
1266       if(!fAODExtension){
1267          if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
1268       }
1269    }
1270
1271    // ----------------- EVENT SELECTION  --------------------------
1272    fHistEvtSelection->Fill(1); // number of events before event selection
1273
1274    if(!fIsKine){
1275       // physics selection
1276       AliInputEventHandler* inputHandler = (AliInputEventHandler*)
1277            ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1278
1279       if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
1280          if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
1281          fHistEvtSelection->Fill(2);
1282          PostData(1, fOutputList);
1283          return;
1284       }
1285
1286    //check AOD pointer
1287       if(!aod){
1288          if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
1289          fHistEvtSelection->Fill(3);
1290          PostData(1, fOutputList);
1291          return;
1292       }
1293
1294       // vertex selection for reconstructed data
1295       AliAODVertex* primVtx = aod->GetPrimaryVertex();
1296
1297       if(!primVtx){
1298          if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
1299          fHistEvtSelection->Fill(3);
1300          PostData(1, fOutputList);
1301          return;
1302       }
1303
1304       Int_t nTracksPrim = primVtx->GetNContributors();
1305       Float_t vtxz = primVtx->GetZ();
1306       //Input events
1307       fhContribVtx->Fill(nTracksPrim);
1308       if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
1309
1310       if((nTracksPrim < fMinContribVtx) ||
1311          (vtxz < fVtxZMin) ||
1312          (vtxz > fVtxZMax)){
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);
1317          return;
1318       }else{
1319       //Accepted events
1320          fhContribVtxAccept->Fill(nTracksPrim);
1321          fhVertexZAccept->Fill(vtxz);
1322       }
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);
1330          return;
1331       }else{
1332          fhVertexZAccept->Fill(zVtx);
1333       }
1334    }
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);
1342    //   return;
1343    //}
1344
1345    //------------------ CENTRALITY SELECTION ---------------
1346    AliCentrality *cent = 0x0;
1347    Double_t centValue  = 0.0;
1348    if(fSystem){  //fSystem=0 for pp,   fSystem=1 for pPb
1349       if(fESD){
1350          cent = fESD->GetCentrality();
1351          if(cent) centValue = cent->GetCentralityPercentile("V0M");
1352       }else{
1353          //centValue = aod->GetHeader()->GetCentrality();
1354          centValue = ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
1355       }
1356       if(fDebug) printf("centrality: %f\n", centValue);
1357       //Input events
1358       fhCentrality->Fill(centValue);
1359
1360       if(centValue < fCentMin || centValue > fCentMax){
1361          fHistEvtSelection->Fill(4);
1362          PostData(1, fOutputList);
1363          return;
1364       }else{
1365          //Accepted events
1366          fhCentralityAccept->Fill( centValue );
1367       }
1368    }
1369
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");
1375
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);
1381          return;
1382        }
1383    }
1384
1385    if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
1386    fHistEvtSelection->Fill(0); // accepted events
1387    // ==================== end event selection ============================
1388
1389    Double_t tmpArrayFive[5];
1390    Double_t tmpArrayFour[4];
1391
1392
1393    TList particleList; //list of tracks
1394    Int_t indexTrigg = -1;
1395    Double_t rhoFromCellMedian=0.0,    rhoCone=0.0;
1396
1397    if(!fIsKine){
1398       //=============== Reconstructed antikt jets ===============
1399       ReadTClonesArray(fJetBranchName.Data()  , fListJets);
1400       ReadTClonesArray(fJetBranchNameBg.Data(), fListJetsBg);
1401
1402       //============ Estimate background in reconstructed events ===========
1403
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);
1407
1408       //Estimate rho from cell median minus jets
1409       EstimateBgRhoMedian(fListJetsBg, &particleList, rhoFromCellMedian,0); //real data
1410    }
1411    //==============  analyze generator level MC  ================
1412    TList particleListGen; //list of tracks in MC
1413
1414    if(fIsChargedMC || fIsKine){
1415
1416       if(fIsKine){  //pure kine
1417
1418          //================= generated charged antikt jets ================
1419          ReadTClonesArray(fJetBranchNameKine.Data(),   fListJetsGen);
1420          ReadTClonesArray(fJetBranchNameBgKine.Data(), fListJetsBgGen);
1421       }else{
1422          //================= generated charged antikt jets ================
1423          ReadTClonesArray(fJetBranchNameChargMC.Data(),   fListJetsGen);
1424          ReadTClonesArray(fJetBranchNameBgChargMC.Data(), fListJetsBgGen);
1425
1426          if(fIsFullMC){ //generated full jets
1427             ReadTClonesArray(fJetBranchNameFullMC.Data(),  fListJetsGenFull);
1428          }
1429       }
1430       //========================================================
1431       //serarch for charged trigger at the MC generator level
1432
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
1442
1443
1444       if(!fIsKine){
1445          if(fESD){//ESD input
1446
1447             AliMCEvent* mcEvent = MCEvent();
1448             if(!mcEvent){
1449                PostData(1, fOutputList);
1450                return;
1451             }
1452
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);
1460             }
1461
1462             for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
1463                if(!mcEvent->IsPhysicalPrimary(it)) continue;
1464                AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
1465                if(!part) continue;
1466                if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1467
1468                   if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1469                      if(indexTriggGen > -1){//trigger candidate was found
1470                         triggersMC[ntriggersMC] = indexTriggGen;
1471                         ntriggersMC++;
1472                      }
1473                   }
1474
1475                   iCounterGen++;//index in particleListGen array
1476                }
1477             }
1478          }else{  //AOD input
1479             if(!fAODIn){
1480                PostData(1, fOutputList);
1481                return;
1482             }
1483             TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
1484             if(!tca){
1485                PostData(1, fOutputList);
1486                return;
1487             }
1488             for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
1489                AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1490                if(!part) continue;
1491                if(!part->IsPhysicalPrimary()) continue;
1492                if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1493
1494                   if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1495                      if(indexTriggGen > -1){ //trigger candidater was found
1496                         triggersMC[ntriggersMC] = indexTriggGen;
1497                         ntriggersMC++;
1498                      }
1499                   }
1500
1501                   iCounterGen++;//number of entries in particleListGen array
1502                }
1503             }
1504          }
1505       }else{ //analyze kine
1506
1507          for(Int_t it = 0; it < fMcEvent->GetNumberOfTracks(); it++){
1508             if(!fMcEvent->IsPhysicalPrimary(it)) continue;
1509             AliMCParticle* part = (AliMCParticle*) fMcEvent->GetTrack(it);
1510             if(!part) continue;
1511             if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
1512
1513                if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
1514                   if(indexTriggGen > -1){//trigger candidate was found
1515                      triggersMC[ntriggersMC] = indexTriggGen;
1516                      ntriggersMC++;
1517                   }
1518                }
1519
1520                iCounterGen++;//index in particleListGen array
1521             }
1522          }
1523       }
1524
1525       if(fHardest==0){
1526          Int_t npar = particleListGen.GetEntries();
1527          for(Int_t ip=0; ip < npar; ip++){
1528             AliVParticle *part = (AliVParticle*) particleListGen.At(ip);
1529             if(!part) continue;
1530
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;
1536                ntriggersMCa++;
1537             }
1538
1539             if(fTriggerPtRangeLow <= ptb && ptb < fTriggerPtRangeHigh  &&
1540                TMath::Abs(part->Eta()) < fTriggerEtaCut  &&  ntriggersMCb<200){
1541                triggersMCb[ntriggersMCb] = ip;
1542                ntriggersMCb++;
1543             }
1544          }
1545
1546          if(ntriggersMCa>0){
1547             Int_t rnda     = fRandom->Integer(ntriggersMCa); //0 to ntriggers-1
1548             Int_t indexTriggGena = triggersMCa[rnda];
1549
1550             Double_t deltaPhia, deltaEtaa, deltaRa;
1551             Int_t aaLow = 0;
1552             Int_t aaHigh = 0;
1553
1554             //Correlation with single inclusive  TRIGGER
1555             AliVParticle* tGenTa = (AliVParticle*) particleListGen.At(indexTriggGena);
1556             if(tGenTa){
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;
1562
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;
1568
1569                   deltaPhia = RelativePhi(tGenTa->Phi(),tGenTz->Phi());
1570                   deltaEtaa = tGenTa->Eta()-tGenTz->Eta();
1571                   deltaRa = sqrt(deltaPhia*deltaPhia + deltaEtaa*deltaEtaa);
1572
1573                   if(tGenTz->Pt()*0.9<20.0){
1574                      fhDeltaRMultTriggersGenALow->Fill(deltaRa);
1575                      fhDeltaPhiMultTriggersGenALow->Fill(deltaPhia);
1576                   }else{
1577                      fhDeltaRMultTriggersGenAHigh->Fill(deltaRa);
1578                      fhDeltaPhiMultTriggersGenAHigh->Fill(deltaPhia);
1579                   }
1580
1581                   if(deltaRa<0.4){
1582                      if(tGenTz->Pt()*0.9<20.0) aaLow++; //15-20
1583                      else aaHigh++;  //20-50
1584                   }
1585                }
1586             }
1587             fhNofMultipleTriggersConeGenA->Fill(aaLow+aaHigh); // 15-50
1588             fhNofMultipleTriggersConeGenALow->Fill(aaLow); //15-20
1589             fhNofMultipleTriggersConeGenAHigh->Fill(aaHigh);//20-50
1590          }
1591
1592          if(ntriggersMCb>0){
1593             Int_t rndb     = fRandom->Integer(ntriggersMCb); //0 to ntriggers-1
1594             Int_t indexTriggGenb = triggersMCb[rndb];
1595
1596             Double_t deltaPhib, deltaEtab, deltaRb;
1597             Int_t bbLow = 0;
1598             Int_t bbHigh = 0;
1599
1600             //Correlation with single inclusive  TRIGGER
1601             AliVParticle* tGenTb = (AliVParticle*) particleListGen.At(indexTriggGenb);
1602             if(tGenTb){
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;
1608
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;
1614
1615                   deltaPhib = RelativePhi(tGenTb->Phi(),tGenTz->Phi());
1616                   deltaEtab = tGenTb->Eta()-tGenTz->Eta();
1617                   deltaRb = sqrt(deltaPhib*deltaPhib + deltaEtab*deltaEtab);
1618
1619                   if(tGenTz->Pt()*0.95<20.0){
1620                      fhDeltaRMultTriggersGenBLow->Fill(deltaRb);
1621                      fhDeltaPhiMultTriggersGenBLow->Fill(deltaPhib);
1622                   }else{
1623                      fhDeltaRMultTriggersGenBHigh->Fill(deltaRb);
1624                      fhDeltaPhiMultTriggersGenBHigh->Fill(deltaPhib);
1625                   }
1626
1627                   if(deltaRb<0.4){
1628                      if(tGenTz->Pt()*0.95<20.0) bbLow++; //15-20
1629                      else bbHigh++;
1630                   }
1631                }
1632             }
1633             fhNofMultipleTriggersConeGenB->Fill(bbLow+bbHigh); //15-50
1634             fhNofMultipleTriggersConeGenBLow->Fill(bbLow);//15-20
1635             fhNofMultipleTriggersConeGenBHigh->Fill(bbHigh);//20-50
1636          }
1637       }
1638
1639
1640
1641       //==============  Estimate bg in generated events ==============
1642       Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
1643
1644       //Estimate rho from cone
1645       EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
1646
1647       //Estimate rho from cell median minus jets
1648       EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
1649
1650       if(ntriggersMC > 0){ //select random inclusive trigger
1651          fhInclTrigCounterGen->Fill(0.5,ntriggersMC); //count inclusive triggers
1652
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()));
1664             }
1665          }
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;
1670
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()));
1675             }
1676          }
1677       }
1678
1679
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];
1685
1686             fhNofMultipleTriggersGen->Fill(ntriggersMC-1);
1687
1688             Double_t deltaPhi, deltaEta, deltaR;
1689             Int_t iLow = 0;
1690             Int_t iHigh = 0;
1691
1692             //Correlation with single inclusive  TRIGGER
1693             AliVParticle* tGenT1 = (AliVParticle*) particleListGen.At(indexTriggGen);
1694             if(tGenT1){
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;
1699
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);
1708
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);
1715                   }
1716
1717                   if(deltaR<0.4){
1718                      if(tGenT2->Pt()<20.0) iLow++;
1719                      else iHigh++;
1720                   }
1721                }
1722             }
1723             fhNofMultipleTriggersConeGen->Fill(iLow+iHigh);
1724             fhNofMultipleTriggersConeGenLow->Fill(iLow);
1725             fhNofMultipleTriggersConeGenHigh->Fill(iHigh);
1726
1727             //Single inclusive trigger within |Delta eta
1728
1729
1730          }else{
1731             indexTriggGen = -1; //trigger not found
1732          }
1733       }
1734
1735       //-----------
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;
1739       //-----------
1740
1741       for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1742          indexTriggGen = igen; //trigger hadron
1743
1744          if(indexTriggGen == -1) continue;
1745          AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
1746          if(!triggerGen) continue;
1747
1748          if(fHardest >= 2){
1749             if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
1750          }
1751          if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1752
1753          ptTriggGen = triggerGen->Pt(); //count triggers
1754          fh2NtriggersGen->Fill(centValue, ptTriggGen);
1755
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));
1759             if(!jet) continue;
1760             Double_t etaJetGen = jet->Eta();
1761             Double_t areaJetGen = jet->EffectiveAreaCharged();
1762
1763             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1764
1765                Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1766                Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1767
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();
1774
1775                //[A,pTjet,pTtrig,dphi]
1776                tmpArrayFour[0] = areaJetGen;
1777                tmpArrayFour[1] = jet->Pt();
1778                tmpArrayFour[2] = ptTriggGen;
1779                tmpArrayFour[3] = dPhiGen;
1780
1781                fHJetPhiCorrGen->Fill(tmpArrayFour); // Dphi distribution jet-triger
1782
1783                //[A,pTjet-pTUe,pTtrig,dphi]
1784                tmpArrayFour[1] = jet->Pt() - ptUeFromCellMedianGen;
1785                fHJetPhiCorrSubUeMedianGen->Fill(tmpArrayFour);
1786
1787                //[A,pTjet-pTUe,pTtrig,dphi]
1788                tmpArrayFour[1] = jet->Pt() - ptUeConeGen;
1789                fHJetPhiCorrSubUeConeGen->Fill(tmpArrayFour);
1790
1791                //Leticia's analysis
1792                Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
1793                if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1794
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);
1802
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);
1810
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);
1818
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);
1826
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);
1834
1835                if(fillOnceGen){
1836                   Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1837                   fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
1838                   fillOnceGen = kFALSE;
1839                }
1840             }//back to back jet-trigger pair
1841          }//jet loop
1842       }//trigger loop
1843
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));
1849             if(!jet) continue;
1850             Double_t etaJetGen = jet->Eta();
1851             Double_t ptJetGen  = jet->Pt();
1852
1853             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1854                fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
1855
1856                Double_t areaJetGen = jet->EffectiveAreaCharged();
1857                Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1858                Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1859
1860                fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
1861                fhJetPtSubUeConeGen->Fill(ptJetGen   - ptUeConeGen);
1862             }
1863          }
1864          if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1865
1866             Int_t ng = (Int_t) fListJetsGen->GetEntries();
1867             Int_t nr = (Int_t) fListJets->GetEntries();
1868
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
1872
1873             if(fDebug){
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());
1876             }
1877             //matching of MC genrator level and reconstructed jets
1878             AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
1879
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
1887
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();
1895
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);
1899                         if(ptJetGen>0){
1900                            fhJetPtResolutionVsPtGen->Fill(ptJetGen,(ptJetRec-ptJetGen)/ptJetGen);
1901                         }
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);
1910
1911                         if((ptJetGen-ptUeConeGen)>0){
1912                            fhJetPtResolutionVsPtConeGen->Fill(ptJetGen-ptUeConeGen,((ptJetRec-ptUeConeRec)- (ptJetGen-ptUeConeGen))/ (ptJetGen-ptUeConeGen));
1913                         }
1914                      }
1915                   }//ig>=0
1916                }//rec jet in eta acceptance
1917             }//loop over reconstructed jets
1918          }// # of  rec jets >0
1919
1920          //=========================== Full jet vs charged jet matrix  ==========================
1921          if(fIsFullMC){
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();
1928
1929                if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1930                   fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
1931                }
1932             }
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();
1936
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
1940
1941                if(fDebug){
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());
1944                }
1945                //matching of MC genrator level and reconstructed jets
1946                AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1947
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
1954
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();
1962
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);
1966                         }
1967                      }//iful>=0
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
1974
1975
1976    if(fIsKine){ //skip reconstructed data analysis in case of kine
1977       PostData(1, fOutputList);
1978       return;
1979    }
1980    //=============  RECONSTRUCTED INCLUSIVE JETS ===============
1981
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;
1989
1990    for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1991       AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1992       if(!jet) continue;
1993       etaJet  = jet->Eta();
1994       phiJet  = jet->Phi();
1995       pTJet   = jet->Pt();
1996       if(pTJet==0) continue;
1997
1998       if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1999       /*areaJet = jet->EffectiveAreaCharged();*/
2000
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;
2009       }
2010
2011       // raw spectra of INCLUSIVE jets
2012       //Centrality, pTjet, A
2013       Double_t fillraw[] = { centValue,
2014                              pTJet,
2015                              areaJet
2016                            };
2017       fHJetPtRaw->Fill(fillraw);*/
2018    }
2019    /*
2020    if(indexLeadingJet > -1){
2021       // raw spectra of LEADING jets
2022       //Centrality, pTjet,  A
2023       Double_t fillleading[] = { centValue,
2024                                  pTLeadingJet,
2025                                  areaLeadingJet
2026                                };
2027       fHLeadingJetPtRaw->Fill(fillleading);
2028    }
2029    */
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;
2035          if(triggerHadron){
2036             Int_t trigLabel = TMath::Abs(triggerHadron->GetLabel());
2037
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;
2042
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
2049                      break;
2050                   }
2051                }
2052             }
2053             if(genTriggIndex>-1){  // make generator level TT x Assoc
2054                AliVParticle *trkGenTT = (AliVParticle*) particleListGen.At(genTriggIndex);
2055                if(trkGenTT){
2056                   fhTriggerCounterGenLevel->Fill(0.5); // count gen triggers that have reconstructed analog
2057                   Double_t deltaPhi,deltaEta,deltaR;
2058
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);
2068
2069                      if(trkGenA->Pt()<20.0){ // gen level assoc particles
2070                         fhDeltaRMultTriggersGenLevelLow->Fill(deltaR);
2071                         fhDeltaPhiMultTriggersGenLevelLow->Fill(deltaPhi);
2072                      }else{
2073                         fhDeltaRMultTriggersGenLevelHigh->Fill(deltaR);
2074                         fhDeltaPhiMultTriggersGenLevelHigh->Fill(deltaPhi);
2075                      }
2076                   }
2077                }
2078             }
2079          }
2080       }
2081    }
2082
2083    // ===============  RECONSTRUCTED TRIGGER-JET PAIRS ================
2084    if(fIsChargedMC && fFillRespMx){
2085      FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
2086    }
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();
2091
2092    for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
2093       indexTrigg = itrk; //trigger hadron with pT >10 GeV
2094
2095       if(indexTrigg < 0) continue;
2096
2097       AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
2098       if(!triggerHadron){
2099          PostData(1, fOutputList);
2100          return;
2101       }
2102       if(fHardest >= 2){
2103          if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
2104       }
2105       if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
2106
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());
2111
2112
2113       //---------- make trigger-jet pairs ---------
2114       //Int_t injet4     = 0;
2115       //Int_t injet      = 0;
2116
2117       for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
2118          AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
2119          if(!jet) continue;
2120          etaJet  = jet->Eta();
2121          phiJet  = jet->Phi();
2122          pTJet   = jet->Pt();
2123          if(pTJet==0) continue;
2124
2125          if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
2126          areaJet = jet->EffectiveAreaCharged();
2127
2128         //subtract bg using estinates base on median of kT jets
2129          Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
2130          Double_t ptUeCone         = rhoCone*areaJet;
2131
2132          //if(areaJet >= 0.07) injet++;
2133          //if(areaJet >= 0.4)  injet4++;
2134
2135          Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
2136          fhDphiTriggerJet->Fill(dphi); //Input
2137
2138          //Dphi versus jet pT
2139          //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
2140          /*Double_t filldp[] = { centValue,
2141                                dphi,
2142                                pTJet,
2143                                triggerHadron->Pt()
2144                              };
2145          fHDphiVsJetPtAll->Fill(filldp);
2146          */
2147          //Rongrong's analysis
2148
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();
2154
2155          //[A,pTjet,pTtrig,dphi]
2156          tmpArrayFour[0] = areaJet;
2157          tmpArrayFour[1] = pTJet;
2158          tmpArrayFour[2] = triggerHadron->Pt();
2159          tmpArrayFour[3] = dPhi;
2160
2161          fHJetPhiCorr->Fill(tmpArrayFour); // Dphi distribution jet-triger
2162
2163          //[A,pTjet-pTUe,pTtrig,dphi]
2164          tmpArrayFour[1] = pTJet - ptUeFromCellMedian;
2165          fHJetPhiCorrSubUeMedian->Fill(tmpArrayFour);
2166
2167          //[A,pTjet-pTUe,pTtrig,dphi]
2168          tmpArrayFour[1] = pTJet - ptUeCone;
2169          fHJetPhiCorrSubUeCone->Fill(tmpArrayFour);
2170
2171          //Leticia's analysis
2172
2173          // Select back to back trigger - jet pairs
2174          if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
2175          fhDphiTriggerJetAccept->Fill(dphi); //Accepted
2176
2177          //NO bg subtraction
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);
2185
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);
2193
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);
2201
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);
2209
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);
2217
2218          if(filledOnce){ //fill for each event only once
2219             Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
2220             fHRhoUeMedianVsCone->Fill(fillRho);
2221             filledOnce = kFALSE;
2222          }
2223       }//jet loop
2224
2225       //Fill Jet Density In the Event A>0.07
2226       /*if(injet>0){
2227          Double_t filldens[]={ (Double_t) particleList.GetEntries(),
2228                             injet/fkAcceptance,
2229                             triggerHadron->Pt()
2230                           };
2231          fHJetDensity->Fill(filldens);
2232       }
2233
2234       //Fill Jet Density In the Event A>0.4
2235       if(injet4>0){
2236          Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
2237                                 injet4/fkAcceptance,
2238                                 triggerHadron->Pt()
2239                               };
2240          fHJetDensityA4->Fill(filldens4);
2241       }*/
2242    }
2243
2244
2245    PostData(1, fOutputList);
2246 }
2247
2248 //----------------------------------------------------------------------------
2249 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
2250 {
2251    // Draw result to the screen
2252    // Called once at the end of the query
2253
2254    if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
2255    if(!GetOutputData(1)) return;
2256 }
2257
2258
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
2263    Int_t iCount        = 0;
2264    AliAODEvent *aodevt = NULL;
2265
2266    if(!fESD) aodevt = fAODIn;
2267    else      aodevt = fAODOut;
2268
2269    if(!aodevt) return -1;
2270
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
2275
2276
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");
2281
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;
2286       list->Add(tr);
2287
2288       //Search trigger:
2289       if(fHardest>0){ //leading particle
2290          if(tr->Pt()>ptmax &&
2291             TMath::Abs((Float_t) tr->Eta()) < fTriggerEtaCut){
2292             ptmax = tr->Pt();
2293             index = iCount;
2294          }
2295       }
2296
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;
2302             ntriggers++;
2303          }
2304       }
2305
2306       iCount++;
2307    }
2308    //inclusive trigger
2309    if(ntriggers>0){ //select random inclusive trigger
2310       fhInclTrigCounter->Fill(0.5,ntriggers); //count inclusive triggers
2311
2312       //Lower pTa bin
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()));
2323          }
2324       }
2325       //Higher pTa bin
2326       for(Int_t it1=0; it1<ntriggers-1; it1++){
2327          AliVParticle* tGent1 = (AliVParticle*) list->At(triggers[it1]);
2328          if(!tGent1) continue;
2329
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()));
2334          }
2335       }
2336    }
2337
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];
2342
2343       fhNofMultipleTriggers->Fill(ntriggers-1);
2344
2345       Double_t deltaPhi, deltaEta, deltaR;
2346       Int_t iLow=0;
2347       Int_t iHigh=0;
2348       //Correlation with single inclusive trigger
2349       AliVParticle* tGent1 = (AliVParticle*) list->At(index);
2350       if(tGent1){
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);
2363
2364             if(tGent2->Pt()<20.0){
2365                fhDeltaRMultTriggersLow->Fill(deltaR);
2366                fhDeltaPhiMultTriggersLow->Fill(deltaPhi);
2367             }else{
2368                fhDeltaRMultTriggersHigh->Fill(deltaR);
2369                fhDeltaPhiMultTriggersHigh->Fill(deltaPhi);
2370             }
2371
2372             if(deltaR<0.4){
2373                if(tGent2->Pt()<20.0) iLow++;
2374                else iHigh++;
2375             }
2376          }
2377       }
2378       fhNofMultipleTriggersCone->Fill(iLow+iHigh);
2379       fhNofMultipleTriggersConeLow->Fill(iLow);
2380       fhNofMultipleTriggersConeHigh->Fill(iHigh);
2381
2382    }
2383
2384    return index;
2385 }
2386
2387 //----------------------------------------------------------------------------
2388 /*
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;
2400
2401       if( (dphi*dphi + deta*deta)< (jetR*jetR)){
2402          pTsum += trk->Pt();
2403       }
2404    }
2405
2406    return pTsum;
2407 }
2408 */
2409 //----------------------------------------------------------------------------
2410
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();
2415
2416    if      (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
2417    else if (mphi > TMath::Pi())  mphi -= TMath::TwoPi();
2418
2419    Double_t dphi = mphi - vphi;
2420    if      (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
2421    else if (dphi > TMath::Pi())  dphi -= TMath::TwoPi();
2422
2423    return dphi;//dphi in [-Pi, Pi]
2424 }
2425
2426
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;
2433    trkList->Add(trk);
2434    if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
2435
2436    //Search trigger:
2437    if(fHardest>0){ //leading particle
2438       if(ptLeading < trk->Pt() &&
2439          TMath::Abs((Float_t) trk->Eta()) < fTriggerEtaCut){
2440          index      = counter;
2441          ptLeading  = trk->Pt();
2442       }
2443    }
2444
2445    if(fHardest==0){ //single inclusive
2446       index = -1;
2447       if(fTriggerPtRangeLow <= trk->Pt() &&
2448             trk->Pt() < fTriggerPtRangeHigh &&
2449              TMath::Abs((Float_t) trk->Eta()) < fTriggerEtaCut){
2450          index      = counter;
2451       }
2452    }
2453
2454    return kTRUE;
2455 }
2456
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;
2469
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());
2476           matched = kTRUE;
2477           break;
2478         }
2479       }
2480       if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
2481    }
2482
2483    return;
2484 }
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
2488
2489    rhoMedian = 0;
2490
2491    //identify leading jet in the full track acceptance
2492    Int_t idxLeadingJet    = -1;
2493    Double_t pTleading     = 0.0;
2494    AliAODJet* jet = NULL;
2495
2496    for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
2497       jet = (AliAODJet*)(listJet->At(ij));
2498       if(!jet) continue;
2499       if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
2500       if(jet->Pt() > pTleading){
2501          idxLeadingJet = ij;
2502          pTleading     = jet->Pt();
2503       }
2504    }
2505
2506    Int_t  njets = 0;
2507    static Double_t jphi[100];
2508    static Double_t jeta[100];
2509    static Double_t jRsquared[100];
2510
2511    for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
2512
2513       jet = (AliAODJet*)(listJet->At(ij));
2514       if(!jet) continue;
2515       if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
2516
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
2521          njets++;
2522       }
2523    }
2524
2525
2526    static Double_t nOutCone[10][4];
2527    static Double_t sumPtOutOfCone[10][4];
2528
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;
2533       }
2534    }
2535
2536    Double_t rndphi, rndeta;
2537    Double_t rndphishift, rndetashift;
2538    Double_t dphi, deta;
2539    Bool_t   bIsInCone;
2540
2541    for(Int_t it=0; it<fnTrials; it++){
2542
2543       rndphi = fRandom->Uniform(0, fPhiSize);
2544       rndeta = fRandom->Uniform(0, fEtaSize);
2545
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;
2550
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]){
2556                   bIsInCone = 1;
2557                   break;
2558                }
2559             }
2560             if(!bIsInCone) nOutCone[ip][ie]++;
2561          }
2562       }
2563    }
2564
2565
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++){
2571
2572       part = (AliVParticle*) listPart->At(ip);
2573       if(!part) continue;
2574       if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
2575
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]){
2581             bIsInCone = 1;
2582             break;
2583          }
2584       }
2585       if(!bIsInCone){
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();
2589       }
2590    }
2591
2592    static Double_t rhoInCells[20];
2593    Double_t  relativeArea;
2594    Int_t  nCells=0;
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;
2600
2601          bufferArea += relativeArea;
2602          bufferPt   += sumPtOutOfCone[ip][ie];
2603          if(bufferArea > fJetFreeAreaFrac){
2604             rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
2605
2606             if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
2607             else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
2608
2609             bufferArea = 0.0;
2610             bufferPt   = 0.0;
2611             nCells++;
2612          }
2613       }
2614    }
2615
2616
2617    if(nCells>0){
2618      rhoMedian = TMath::Median(nCells, rhoInCells);
2619      if(mode==1){ //mc data
2620         fhEntriesToMedianGen->Fill(nCells);
2621      }else{ //real data
2622         fhEntriesToMedian->Fill(nCells);
2623      }
2624    }
2625
2626 }
2627
2628 //__________________________________________________________________
2629 void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart,  Double_t &rhoPerpCone){
2630
2631    rhoPerpCone = 0.0; //init
2632
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.;
2639
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();
2648       }
2649    }
2650
2651    Double_t phileftcone  = phiLeading + TMath::Pi()/2;
2652    Double_t phirightcone = phiLeading - TMath::Pi()/2;
2653    Double_t dp, de;
2654
2655    for(Int_t ip=0; ip < npart; ip++){
2656
2657       AliVParticle *part = (AliVParticle*) listPart->At(ip);
2658       if(!part){
2659          continue;
2660       }
2661
2662       dp = RelativePhi(phileftcone, part->Phi());
2663       de = etaLeading - part->Eta();
2664       if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2665
2666       dp = RelativePhi(phirightcone, part->Phi());
2667       if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2668
2669    }
2670
2671    //normalize total pT by two times cone are
2672    rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
2673
2674 }
2675 //__________________________________________________________________
2676
2677 void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
2678    //Convert TClones array of jets to TList
2679
2680    if(!strlen(bname.Data())){
2681       AliError(Form("Jet branch %s not set.", bname.Data()));
2682       return;
2683    }
2684
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
2689       }
2690       if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
2691         array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
2692       }
2693    }else{ //take input from AOD
2694       if(fAODOut&&!array){
2695          array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
2696       }
2697       if(fAODExtension&&!array){
2698          array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
2699       }
2700       if(fAODIn&&!array){
2701          array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
2702       }
2703    }
2704
2705    list->Clear(); //clear list beforehand
2706
2707    if(array){
2708       if(fDebug){
2709          Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
2710       }
2711
2712       for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
2713          AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
2714          if (jet) list->Add(jet);
2715       }
2716    }
2717
2718    return;
2719 }