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