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                   if(indexTriggGena == triggersMCa[ia]) continue;
1332
1333                   AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCa[ia]);
1334                   if(!tGenTz) continue;
1335
1336                   deltaPhia = RelativePhi(tGenTa->Phi(),tGenTz->Phi());
1337                   deltaEtaa = tGenTa->Eta()-tGenTz->Eta();
1338                   deltaRa = sqrt(deltaPhia*deltaPhia + deltaEtaa*deltaEtaa);
1339
1340                   if(deltaRa<0.4) aa++;
1341                }
1342             }
1343             fhNofMultipleTriggersConeGenA->Fill(aa);
1344          }
1345
1346          if(ntriggersMCb>0){
1347             Int_t rndb     = fRandom->Integer(ntriggersMCb); //0 to ntriggers-1
1348             Int_t indexTriggGenb = triggersMCb[rndb];
1349
1350             Double_t deltaPhib, deltaEtab, deltaRb;
1351             Int_t bb = 0;
1352
1353             //Correlation with single inclusive  TRIGGER
1354             AliVParticle* tGenTb = (AliVParticle*) particleListGen.At(indexTriggGenb);
1355             if(tGenTb){
1356                for(Int_t ib=0; ib<ntriggersMCb; ib++){
1357                   if(indexTriggGenb == triggersMCb[ib]) continue;
1358
1359                   AliVParticle* tGenTz = (AliVParticle*) particleListGen.At(triggersMCb[ib]);
1360                   if(!tGenTz) continue;
1361
1362                   deltaPhib = RelativePhi(tGenTb->Phi(),tGenTz->Phi());
1363                   deltaEtab = tGenTb->Eta()-tGenTz->Eta();
1364                   deltaRb = sqrt(deltaPhib*deltaPhib + deltaEtab*deltaEtab);
1365
1366                   if(deltaRb<0.4) bb++;
1367                }
1368             }
1369             fhNofMultipleTriggersConeGenB->Fill(bb);
1370          }
1371       }
1372
1373
1374
1375       //==============  Estimate bg in generated events ==============
1376       Double_t rhoFromCellMedianGen=0.0, rhoConeGen=0.0;
1377
1378       //Estimate rho from cone
1379       EstimateBgCone(fListJetsGen, &particleListGen, rhoConeGen);
1380
1381       //Estimate rho from cell median minus jets
1382       EstimateBgRhoMedian(fListJetsBgGen, &particleListGen, rhoFromCellMedianGen,1);//mc data
1383
1384       //============  Generator trigger+jet ==================
1385       if(fHardest==0){ //single inclusive trigger
1386          if(ntriggersMC>0){ //there is at least one trigger
1387             Int_t rnd     = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
1388             indexTriggGen = triggersMC[rnd];
1389
1390             fhNofMultipleTriggersGen->Fill(ntriggersMC-1);
1391
1392             Double_t deltaPhi, deltaEta, deltaR;
1393             Int_t k = 0;
1394
1395             //Correlation with single inclusive  TRIGGER
1396             AliVParticle* tGenT1 = (AliVParticle*) particleListGen.At(indexTriggGen);
1397             if(tGenT1){
1398                for(Int_t ia=0; ia<ntriggersMC; ia++){
1399                   if(indexTriggGen == triggersMC[ia]) continue;
1400
1401                   AliVParticle* tGenT2 = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1402                   if(!tGenT2) continue;
1403
1404                   deltaPhi = RelativePhi(tGenT1->Phi(),tGenT2->Phi());
1405                   deltaEta = tGenT1->Eta()-tGenT2->Eta();
1406                   deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1407
1408                   if(deltaR<0.4) k++;
1409                }
1410             }
1411             fhNofMultipleTriggersConeGen->Fill(k);
1412
1413             if(ntriggersMC>1){
1414                //Correlation of each trigger with any other trigger
1415                for(Int_t ia=0; ia<ntriggersMC-1; ia++){
1416                   AliVParticle* tGenI = (AliVParticle*) particleListGen.At(triggersMC[ia]);
1417                   if(!tGenI) continue;
1418                   for(Int_t ib=ia+1; ib<ntriggersMC; ib++){
1419                      AliVParticle* tGenII = (AliVParticle*) particleListGen.At(triggersMC[ib]);
1420                      if(!tGenII) continue;
1421
1422                         deltaPhi = RelativePhi(tGenI->Phi(),tGenII->Phi());
1423                         deltaEta = tGenI->Eta()-tGenII->Eta();
1424                         deltaR = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1425                         fhDeltaRMultTriggersGen->Fill(deltaR);
1426                   }
1427                }
1428             }
1429
1430          }else{
1431             indexTriggGen = -1; //trigger not found
1432          }
1433       }
1434
1435       //-----------
1436       Int_t ilowGen  = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
1437       Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
1438       Bool_t fillOnceGen = kTRUE;
1439       //-----------
1440
1441       for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
1442          indexTriggGen = igen; //trigger hadron
1443
1444          if(indexTriggGen == -1) continue;
1445          AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
1446          if(!triggerGen) continue;
1447
1448          if(fHardest >= 2){
1449             if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
1450          }
1451          if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
1452
1453          ptTriggGen = triggerGen->Pt(); //count triggers
1454          fh2NtriggersGen->Fill(centValue, ptTriggGen);
1455
1456          //Count jets and trigger-jet pairs at MC  generator level
1457          for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1458             AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1459             if(!jet) continue;
1460             Double_t etaJetGen = jet->Eta();
1461             Double_t areaJetGen = jet->EffectiveAreaCharged();
1462
1463             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1464
1465                Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1466                Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1467
1468                //Rongrong's analysis
1469                Double_t dPhiGen = jet->Phi() - triggerGen->Phi();
1470                if(dPhiGen>2*TMath::Pi())    dPhiGen -= 2*TMath::Pi();
1471                if(dPhiGen<-2*TMath::Pi())   dPhiGen += 2*TMath::Pi();
1472                if(dPhiGen<-0.5*TMath::Pi()) dPhiGen += 2*TMath::Pi();
1473                if(dPhiGen>1.5*TMath::Pi())  dPhiGen -= 2*TMath::Pi();
1474
1475                //[A,pTjet,pTtrig,dphi]
1476                tmpArrayFour[0] = areaJetGen;
1477                tmpArrayFour[1] = jet->Pt();
1478                tmpArrayFour[2] = ptTriggGen;
1479                tmpArrayFour[3] = dPhiGen;
1480
1481                fHJetPhiCorrGen->Fill(tmpArrayFour); // Dphi distribution jet-triger
1482
1483                //[A,pTjet-pTUe,pTtrig,dphi]
1484                tmpArrayFour[1] = jet->Pt() - ptUeFromCellMedianGen;
1485                fHJetPhiCorrSubUeMedianGen->Fill(tmpArrayFour);
1486
1487                //[A,pTjet-pTUe,pTtrig,dphi]
1488                tmpArrayFour[1] = jet->Pt() - ptUeConeGen;
1489                fHJetPhiCorrSubUeConeGen->Fill(tmpArrayFour);
1490
1491                //Leticia's analysis
1492                Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
1493                if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1494
1495                //Centrality, A, pT, pTtrigg
1496                tmpArrayFive[0] = centValue;
1497                tmpArrayFive[1] = areaJetGen;
1498                tmpArrayFive[2] = jet->Pt();
1499                tmpArrayFive[3] = ptTriggGen;
1500                tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1501                fHJetSpecGen->Fill(tmpArrayFive);
1502
1503                //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1504                tmpArrayFive[0] = centValue;
1505                tmpArrayFive[1] = areaJetGen;
1506                tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1507                tmpArrayFive[3] = ptTriggGen;
1508                tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1509                fHJetSpecSubUeMedianGen->Fill(tmpArrayFive);
1510
1511                //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1512                tmpArrayFive[0] = centValue;
1513                tmpArrayFive[1] = areaJetGen;
1514                tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1515                tmpArrayFive[3] = ptTriggGen;
1516                tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1517                fHJetSpecSubUeConeGen->Fill(tmpArrayFive);
1518
1519                //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1520                tmpArrayFive[0] = areaJetGen;
1521                tmpArrayFive[1] = jet->Pt();
1522                tmpArrayFive[2] = jet->Pt() - ptUeFromCellMedianGen;
1523                tmpArrayFive[3] = ptUeFromCellMedianGen;
1524                tmpArrayFive[4] = rhoFromCellMedianGen;
1525                fHJetUeMedianGen->Fill(tmpArrayFive);
1526
1527                //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", perpendicular Cone
1528                tmpArrayFive[0] = areaJetGen;
1529                tmpArrayFive[1] = jet->Pt();
1530                tmpArrayFive[2] = jet->Pt() - ptUeConeGen;
1531                tmpArrayFive[3] = ptUeConeGen;
1532                tmpArrayFive[4] = rhoConeGen;
1533                fHJetUeConeGen->Fill(tmpArrayFive);
1534
1535                if(fillOnceGen){
1536                   Double_t fillRhoGen[] = { rhoConeGen, rhoFromCellMedianGen};
1537                   fHRhoUeMedianVsConeGen->Fill(fillRhoGen);
1538                   fillOnceGen = kFALSE;
1539                }
1540             }//back to back jet-trigger pair
1541          }//jet loop
1542       }//trigger loop
1543
1544       if(fFillRespMx && !fIsKine){
1545          //================ RESPONSE MATRIX ===============
1546          //Count jets and trigger-jet pairs at MC  generator level
1547          for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
1548             AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
1549             if(!jet) continue;
1550             Double_t etaJetGen = jet->Eta();
1551             Double_t ptJetGen  = jet->Pt();
1552
1553             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
1554                fhJetPtGen->Fill(ptJetGen); // gen level pt spectum of jets response mx normalization
1555
1556                Double_t areaJetGen = jet->EffectiveAreaCharged();
1557                Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1558                Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1559
1560                fhJetPtSubUeMedianGen->Fill(ptJetGen - ptUeFromCellMedianGen);
1561                fhJetPtSubUeConeGen->Fill(ptJetGen   - ptUeConeGen);
1562             }
1563          }
1564          if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
1565
1566             Int_t ng = (Int_t) fListJetsGen->GetEntries();
1567             Int_t nr = (Int_t) fListJets->GetEntries();
1568
1569             //Find closest MC generator - reconstructed jets
1570             if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
1571             if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
1572
1573             if(fDebug){
1574                Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
1575                Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
1576             }
1577             //matching of MC genrator level and reconstructed jets
1578             AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
1579
1580             // Fill response matrix
1581             for(Int_t ir = 0; ir < nr; ir++){
1582                AliAODJet *recJet   = (AliAODJet*) fListJets->At(ir);
1583                Double_t etaJetRec  = recJet->Eta();
1584                Double_t ptJetRec   = recJet->Pt();
1585                Double_t areaJetRec = recJet->EffectiveAreaCharged();
1586                //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1587
1588                if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
1589                   Int_t ig = faGenIndex[ir]; //associated generator level jet
1590                   if(ig >= 0 && ig < ng){
1591                      if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1592                      AliAODJet *genJet  = (AliAODJet*) fListJetsGen->At(ig);
1593                      Double_t ptJetGen  = genJet->Pt();
1594                      Double_t etaJetGen = genJet->Eta();
1595
1596                      //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1597                      if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
1598                         fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
1599
1600                         Double_t areaJetGen  = genJet->EffectiveAreaCharged();
1601                         Double_t ptUeFromCellMedianGen = rhoFromCellMedianGen*areaJetGen;
1602                         Double_t ptUeConeGen         = rhoConeGen*areaJetGen;
1603                         Double_t ptUeFromCellMedianRec = rhoFromCellMedian*areaJetRec;
1604                         Double_t ptUeConeRec         = rhoCone*areaJetRec;
1605                         fhJetPtGenVsJetPtRecSubUeMedian->Fill(ptJetRec-ptUeFromCellMedianRec,
1606                                                            ptJetGen-ptUeFromCellMedianGen);
1607                         fhJetPtGenVsJetPtRecSubUeCone->Fill(ptJetRec-ptUeConeRec, ptJetGen-ptUeConeGen);
1608                      }
1609                   }//ig>=0
1610                }//rec jet in eta acceptance
1611             }//loop over reconstructed jets
1612          }// # of  rec jets >0
1613
1614          //=========================== Full jet vs charged jet matrix  ==========================
1615          if(fIsFullMC){
1616             //Count full jets and charged-jet pairs at MC  generator level
1617             for(Int_t ij=0; ij<fListJetsGenFull->GetEntries(); ij++){
1618                AliAODJet* jetFull = (AliAODJet*)(fListJetsGenFull->At(ij));
1619                if(!jetFull) continue;
1620                Double_t etaJetFull = jetFull->Eta();
1621                Double_t ptJetFull  = jetFull->Pt();
1622
1623                if((fJetEtaMin<=etaJetFull) && (etaJetFull<=fJetEtaMax)){
1624                   fhJetPtGenFull->Fill(ptJetFull); // generator level pt spectum of full jets
1625                }
1626             }
1627             if(fListJetsGen->GetEntries()>0 && fListJetsGenFull->GetEntries()>0){ //at least some reconstructed jets
1628                Int_t nful = (Int_t) fListJetsGenFull->GetEntries();
1629                Int_t nchr = (Int_t) fListJetsGen->GetEntries();
1630
1631                //Find closest MC generator full - charged jet
1632                if(faGenIndex.GetSize()<nchr) faGenIndex.Set(nchr); //idx of gen FULL jet assoc to gen CHARGED jet
1633                if(faRecIndex.GetSize()<nful) faRecIndex.Set(nful); //idx of gen CHARGED jet assoc to gen FULL jet
1634
1635                if(fDebug){
1636                   Printf("New Charg List %d Full index Array %d",nchr,faGenIndex.GetSize());
1637                   Printf("New Full List %d Charg index Array %d",nful,faRecIndex.GetSize());
1638                }
1639                //matching of MC genrator level and reconstructed jets
1640                AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGenFull,nful,fListJetsGen,nchr,faGenIndex,faRecIndex,fDebug);
1641
1642               // Fill response matrix
1643                for(Int_t ichr = 0; ichr < nchr; ichr++){ //charged jet loop
1644                   AliAODJet *chJet  = (AliAODJet*) fListJetsGen->At(ichr);
1645                   Double_t etaJetCh = chJet->Eta();
1646                   Double_t ptJetCh  = chJet->Pt();
1647                   //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1648
1649                   if((fJetEtaMin <= etaJetCh) && (etaJetCh <= fJetEtaMax)){
1650                      Int_t iful = faGenIndex[ichr]; //associated generator level jet
1651                      if(iful >= 0 && iful < nful){
1652                         if(fDebug > 10) Printf("%s:%d iful = %d ichr = %d",(char*)__FILE__,__LINE__,iful,ichr);
1653                         AliAODJet *genJetFull  = (AliAODJet*) fListJetsGenFull->At(iful);
1654                         Double_t ptJetFull  = genJetFull->Pt();
1655                         Double_t etaJetFull = genJetFull->Eta();
1656
1657                         //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
1658                         if((fJetEtaMin <= etaJetFull) && (etaJetFull <= fJetEtaMax)){
1659                            fhJetPtGenChargVsJetPtGenFull->Fill(ptJetFull,ptJetCh);
1660                         }
1661                      }//iful>=0
1662                   }//rec jet in eta acceptance
1663                }//loop over reconstructed jets
1664             }// # of  rec jets >0
1665          }//pointer MC generator jets
1666       } //fill resp mx only for bin
1667    }//analyze generator level MC
1668
1669
1670    if(fIsKine){ //skip reconstructed data analysis in case of kine
1671       PostData(1, fOutputList);
1672       return;
1673    }
1674    //=============  RECONSTRUCTED INCLUSIVE JETS ===============
1675
1676    Double_t etaJet  = 0.0;
1677    Double_t pTJet   = 0.0;
1678    Double_t areaJet = 0.0;
1679    Double_t phiJet  = 0.0;
1680    //Int_t indexLeadingJet     = -1;
1681    //Double_t pTLeadingJet     = -10.0;
1682    //Double_t areaLeadingJet   = -10.0;
1683
1684    for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1685       AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1686       if(!jet) continue;
1687       etaJet  = jet->Eta();
1688       phiJet  = jet->Phi();
1689       pTJet   = jet->Pt();
1690       if(pTJet==0) continue;
1691
1692       if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1693       /*areaJet = jet->EffectiveAreaCharged();*/
1694
1695       //Jet Diagnostics---------------------------------
1696       fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
1697       fhJetEta->Fill(pTJet, etaJet);
1698       //search for leading jet
1699       /*if(pTJet > pTLeadingJet){
1700          indexLeadingJet  = ij;
1701          pTLeadingJet     = pTJet;
1702          areaLeadingJet   = areaJet;
1703       }
1704
1705       // raw spectra of INCLUSIVE jets
1706       //Centrality, pTjet, A
1707       Double_t fillraw[] = { centValue,
1708                              pTJet,
1709                              areaJet
1710                            };
1711       fHJetPtRaw->Fill(fillraw);*/
1712    }
1713    /*
1714    if(indexLeadingJet > -1){
1715       // raw spectra of LEADING jets
1716       //Centrality, pTjet,  A
1717       Double_t fillleading[] = { centValue,
1718                                  pTLeadingJet,
1719                                  areaLeadingJet
1720                                };
1721       fHLeadingJetPtRaw->Fill(fillleading);
1722    }
1723    */
1724
1725    // ===============  RECONSTRUCTED TRIGGER-JET PAIRS ================
1726    if(fIsChargedMC && fFillRespMx){
1727      FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1728    }
1729    Bool_t filledOnce = kTRUE; //fill rho histogram only once per event
1730    //set ranges of the trigger loop
1731    Int_t ilow  = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1732    Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
1733
1734    for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1735       indexTrigg = itrk; //trigger hadron with pT >10 GeV
1736
1737       if(indexTrigg < 0) continue;
1738
1739       AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1740       if(!triggerHadron){
1741          PostData(1, fOutputList);
1742          return;
1743       }
1744       if(fHardest >= 2){
1745          if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
1746       }
1747       if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
1748
1749       //Fill trigger histograms
1750       fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1751       fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1752       fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
1753
1754
1755       //---------- make trigger-jet pairs ---------
1756       //Int_t injet4     = 0;
1757       //Int_t injet      = 0;
1758
1759       for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1760          AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1761          if(!jet) continue;
1762          etaJet  = jet->Eta();
1763          phiJet  = jet->Phi();
1764          pTJet   = jet->Pt();
1765          if(pTJet==0) continue;
1766
1767          if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1768          areaJet = jet->EffectiveAreaCharged();
1769
1770         //subtract bg using estinates base on median of kT jets
1771          Double_t ptUeFromCellMedian = rhoFromCellMedian*areaJet;
1772          Double_t ptUeCone         = rhoCone*areaJet;
1773
1774          //if(areaJet >= 0.07) injet++;
1775          //if(areaJet >= 0.4)  injet4++;
1776
1777          Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1778          fhDphiTriggerJet->Fill(dphi); //Input
1779
1780          //Dphi versus jet pT
1781          //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
1782          /*Double_t filldp[] = { centValue,
1783                                dphi,
1784                                pTJet,
1785                                triggerHadron->Pt()
1786                              };
1787          fHDphiVsJetPtAll->Fill(filldp);
1788          */
1789          //Rongrong's analysis
1790
1791          Double_t dPhi = phiJet - triggerHadron->Phi();
1792          if(dPhi>2*TMath::Pi())    dPhi -= 2*TMath::Pi();
1793          if(dPhi<-2*TMath::Pi())   dPhi += 2*TMath::Pi();
1794          if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
1795          if(dPhi>1.5*TMath::Pi())  dPhi -= 2*TMath::Pi();
1796
1797          //[A,pTjet,pTtrig,dphi]
1798          tmpArrayFour[0] = areaJet;
1799          tmpArrayFour[1] = pTJet;
1800          tmpArrayFour[2] = triggerHadron->Pt();
1801          tmpArrayFour[3] = dPhi;
1802
1803          fHJetPhiCorr->Fill(tmpArrayFour); // Dphi distribution jet-triger
1804
1805          //[A,pTjet-pTUe,pTtrig,dphi]
1806          tmpArrayFour[1] = pTJet - ptUeFromCellMedian;
1807          fHJetPhiCorrSubUeMedian->Fill(tmpArrayFour);
1808
1809          //[A,pTjet-pTUe,pTtrig,dphi]
1810          tmpArrayFour[1] = pTJet - ptUeCone;
1811          fHJetPhiCorrSubUeCone->Fill(tmpArrayFour);
1812
1813          //Leticia's analysis
1814
1815          // Select back to back trigger - jet pairs
1816          if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1817          fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1818
1819          //NO bg subtraction
1820          //Centrality, A, pTjet, pTtrigg, dphi
1821          tmpArrayFive[0] = centValue;
1822          tmpArrayFive[1] = areaJet;
1823          tmpArrayFive[2] = pTJet;
1824          tmpArrayFive[3] = triggerHadron->Pt();
1825          tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1826          fHJetSpec->Fill(tmpArrayFive);
1827
1828           //Centrality, A, pTjet-pTbgCellMedian, pTtrigg, dphi
1829          tmpArrayFive[0] = centValue;
1830          tmpArrayFive[1] = areaJet;
1831          tmpArrayFive[2] = pTJet-ptUeFromCellMedian;
1832          tmpArrayFive[3] = triggerHadron->Pt();
1833          tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1834          fHJetSpecSubUeMedian->Fill(tmpArrayFive);
1835
1836          //Centrality, A, pTjet-pTbgCone, pTtrigg, dphi
1837          tmpArrayFive[0] = centValue;
1838          tmpArrayFive[1] = areaJet;
1839          tmpArrayFive[2] = pTJet - ptUeCone;
1840          tmpArrayFive[3] = triggerHadron->Pt();
1841          tmpArrayFive[4] = TMath::Abs((Double_t) dphi);
1842          fHJetSpecSubUeCone->Fill(tmpArrayFive);
1843
1844          //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", kT median
1845          tmpArrayFive[0] = areaJet;
1846          tmpArrayFive[1] = pTJet;
1847          tmpArrayFive[2] = pTJet - ptUeFromCellMedian;
1848          tmpArrayFive[3] = ptUeFromCellMedian;
1849          tmpArrayFive[4] = rhoFromCellMedian;
1850          fHJetUeMedian->Fill(tmpArrayFive);
1851
1852          //Ue diagnostics  "[A,pTjet,pTjet-pTUe, pTUe, rhoUe]", Cone median
1853          tmpArrayFive[0] = areaJet;
1854          tmpArrayFive[1] = pTJet;
1855          tmpArrayFive[2] = pTJet - ptUeCone;
1856          tmpArrayFive[3] = ptUeCone;
1857          tmpArrayFive[4] = rhoCone;
1858          fHJetUeCone->Fill(tmpArrayFive);
1859
1860          if(filledOnce){ //fill for each event only once
1861             Double_t fillRho[] = { rhoCone,rhoFromCellMedian};
1862             fHRhoUeMedianVsCone->Fill(fillRho);
1863             filledOnce = kFALSE;
1864          }
1865       }//jet loop
1866
1867       //Fill Jet Density In the Event A>0.07
1868       /*if(injet>0){
1869          Double_t filldens[]={ (Double_t) particleList.GetEntries(),
1870                             injet/fkAcceptance,
1871                             triggerHadron->Pt()
1872                           };
1873          fHJetDensity->Fill(filldens);
1874       }
1875
1876       //Fill Jet Density In the Event A>0.4
1877       if(injet4>0){
1878          Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1879                                 injet4/fkAcceptance,
1880                                 triggerHadron->Pt()
1881                               };
1882          fHJetDensityA4->Fill(filldens4);
1883       }*/
1884    }
1885
1886
1887    PostData(1, fOutputList);
1888 }
1889
1890 //----------------------------------------------------------------------------
1891 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1892 {
1893    // Draw result to the screen
1894    // Called once at the end of the query
1895
1896    if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1897    if(!GetOutputData(1)) return;
1898 }
1899
1900
1901 //----------------------------------------------------------------------------
1902 Int_t  AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1903    //Fill the list of accepted tracks (passed track cut)
1904    //return consecutive index of the hardest ch hadron in the list
1905    Int_t iCount        = 0;
1906    AliAODEvent *aodevt = NULL;
1907
1908    if(!fESD) aodevt = fAODIn;
1909    else      aodevt = fAODOut;
1910
1911    if(!aodevt) return -1;
1912
1913    Int_t    index = -1; //index of the highest particle in the list
1914    Double_t ptmax = -10;
1915    Int_t    triggers[200];
1916    Int_t    ntriggers = 0; //index in triggers array
1917
1918
1919    for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
1920       //AliAODTrack *tr = aodevt->GetTrack(it);
1921       AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aodevt->GetTrack(it));
1922       if(!tr) AliFatal("Not a standard AOD");
1923
1924       if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1925       //if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1926       if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1927       if(tr->Pt() < fTrackLowPtCut) continue;
1928       list->Add(tr);
1929
1930       //Search trigger:
1931       if(fHardest>0){ //leading particle
1932          if(tr->Pt()>ptmax){
1933             ptmax = tr->Pt();
1934             index = iCount;
1935          }
1936       }
1937
1938       if(fHardest==0 && ntriggers<200){ //single inclusive
1939          if(fTriggerPtRangeLow <= tr->Pt() &&
1940             tr->Pt() < fTriggerPtRangeHigh){
1941             triggers[ntriggers] = iCount;
1942             ntriggers++;
1943          }
1944       }
1945
1946       iCount++;
1947    }
1948
1949    if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1950       Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1951       index     = triggers[rnd];
1952
1953       fhNofMultipleTriggers->Fill(ntriggers-1);
1954
1955       Double_t deltaPhi, deltaEta, deltaR;
1956       Int_t k=0;
1957       //Correlation with single inclusive trigger
1958       AliVParticle* tGent1 = (AliVParticle*) list->At(index);
1959       if(tGent1){
1960          for(Int_t ia=0; ia<ntriggers; ia++){
1961             if(triggers[ia]==index) continue;
1962             AliVParticle* tGent2 = (AliVParticle*) list->At(triggers[ia]);
1963             if(!tGent2) continue;
1964             deltaPhi = RelativePhi(tGent1->Phi(),tGent2->Phi());
1965             deltaEta = tGent1->Eta()-tGent2->Eta();
1966             deltaR   = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1967             if(deltaR<0.4) k++;
1968          }
1969       }
1970       fhNofMultipleTriggersCone->Fill(k);
1971
1972       if(ntriggers>1){
1973          //Correlation with any other trigger
1974          for(Int_t ia=0; ia<ntriggers-1; ia++){
1975             AliVParticle* tGeni = (AliVParticle*) list->At(triggers[ia]);
1976             if(!tGeni) continue;
1977             for(Int_t ib=ia+1; ib<ntriggers; ib++){
1978                AliVParticle* tGenii = (AliVParticle*) list->At(triggers[ib]);
1979                if(!tGenii) continue;
1980                deltaPhi = RelativePhi(tGeni->Phi(),tGenii->Phi());
1981                deltaEta = tGeni->Eta()-tGenii->Eta();
1982                deltaR   = sqrt(deltaPhi*deltaPhi + deltaEta*deltaEta);
1983                fhDeltaRMultTriggers->Fill(deltaR);
1984             }
1985          }
1986       }
1987    }
1988
1989    return index;
1990 }
1991
1992 //----------------------------------------------------------------------------
1993 /*
1994 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1995    //calculate sum of track pT in the cone perpendicular in phi to the jet
1996    //jetR = cone radius
1997    // jetPhi, jetEta = direction of the jet
1998    Int_t numberOfTrks = trkList->GetEntries();
1999    Double_t pTsum = 0.0;
2000    Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
2001    for(Int_t it=0; it<numberOfTrks; it++){
2002       AliVParticle *trk = (AliVParticle*) trkList->At(it);
2003       Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
2004       Double_t deta = trk->Eta()-jetEta;
2005
2006       if( (dphi*dphi + deta*deta)< (jetR*jetR)){
2007          pTsum += trk->Pt();
2008       }
2009    }
2010
2011    return pTsum;
2012 }
2013 */
2014 //----------------------------------------------------------------------------
2015
2016 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
2017    //Get relative azimuthal angle of two particles -pi to pi
2018    if      (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
2019    else if (vphi > TMath::Pi())  vphi -= TMath::TwoPi();
2020
2021    if      (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
2022    else if (mphi > TMath::Pi())  mphi -= TMath::TwoPi();
2023
2024    Double_t dphi = mphi - vphi;
2025    if      (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
2026    else if (dphi > TMath::Pi())  dphi -= TMath::TwoPi();
2027
2028    return dphi;//dphi in [-Pi, Pi]
2029 }
2030
2031
2032 //----------------------------------------------------------------------------
2033 Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
2034    //fill track efficiency denominator
2035    if(trk->Pt() < fTrackLowPtCut) return kFALSE;
2036    if(trk->Charge()==0)        return kFALSE;
2037    if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
2038    trkList->Add(trk);
2039    if(fFillRespMx) fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
2040
2041    //Search trigger:
2042    if(fHardest>0){ //leading particle
2043       if(ptLeading < trk->Pt()){
2044          index      = counter;
2045          ptLeading  = trk->Pt();
2046       }
2047    }
2048
2049    if(fHardest==0){ //single inclusive
2050       index = -1;
2051       if(fTriggerPtRangeLow <= trk->Pt() &&
2052             trk->Pt() < fTriggerPtRangeHigh){
2053          index      = counter;
2054       }
2055    }
2056
2057    return kTRUE;
2058 }
2059
2060 //----------------------------------------------------------------------------
2061 void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
2062    //fill track efficiency numerator
2063    if(!recList) return;
2064    if(!genList) return;
2065    Int_t nRec = recList->GetEntries();
2066    Int_t nGen = genList->GetEntries();
2067    for(Int_t ir=0; ir<nRec; ir++){
2068       AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
2069       if(!trkRec) continue;
2070       Int_t recLabel = TMath::Abs(trkRec->GetLabel());
2071       Bool_t matched = kFALSE;
2072
2073       for(Int_t ig=0; ig<nGen; ig++){
2074         AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
2075         if(!trkGen) continue;
2076         Int_t genLabel = TMath::Abs(trkGen->GetLabel());
2077         if(recLabel==genLabel){
2078           fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
2079           matched = kTRUE;
2080           break;
2081         }
2082       }
2083       if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
2084    }
2085
2086    return;
2087 }
2088 //________________________________________________________________
2089 void AliAnalysisTaskJetCorePP::EstimateBgRhoMedian(TList *listJet, TList* listPart,  Double_t &rhoMedian, Int_t mode){
2090    //Estimate background rho by means of integrating track pT outside identified jet cones
2091
2092    rhoMedian = 0;
2093
2094    //identify leading jet in the full track acceptance
2095    Int_t idxLeadingJet    = -1;
2096    Double_t pTleading     = 0.0;
2097    AliAODJet* jet = NULL;
2098
2099    for(Int_t ij = 0; ij < listJet->GetEntries(); ij++){ //loop over bg kT jets
2100       jet = (AliAODJet*)(listJet->At(ij));
2101       if(!jet) continue;
2102       if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue; //track acceptance cut
2103       if(jet->Pt() > pTleading){
2104          idxLeadingJet = ij;
2105          pTleading     = jet->Pt();
2106       }
2107    }
2108
2109    Int_t  njets = 0;
2110    static Double_t jphi[100];
2111    static Double_t jeta[100];
2112    static Double_t jRsquared[100];
2113
2114    for(Int_t ij = 0; ij< listJet->GetEntries(); ij++){ //loop over bg kT jets
2115
2116       jet = (AliAODJet*)(listJet->At(ij));
2117       if(!jet) continue;
2118       if(TMath::Abs(jet->Eta()) > fTrackEtaCut) continue;
2119
2120       if((jet->Pt() > fBgMaxJetPt || ij== idxLeadingJet) && njets<100){
2121          jphi[njets]      = RelativePhi(jet->Phi(),0.0); //-pi,pi
2122          jeta[njets]      = jet->Eta();
2123          jRsquared[njets] = fSafetyMargin*jet->EffectiveAreaCharged()/TMath::Pi(); //scale jet radius
2124          njets++;
2125       }
2126    }
2127
2128
2129    static Double_t nOutCone[10][4];
2130    static Double_t sumPtOutOfCone[10][4];
2131
2132    for(Int_t ie=0; ie<fnEta; ie++){
2133       for(Int_t ip=0; ip<fnPhi; ip++){
2134          nOutCone[ip][ie]       = 0.0;     //initialize counter
2135          sumPtOutOfCone[ip][ie] = 0.0;
2136       }
2137    }
2138
2139    Double_t rndphi, rndeta;
2140    Double_t rndphishift, rndetashift;
2141    Double_t dphi, deta;
2142    Bool_t   bIsInCone;
2143
2144    for(Int_t it=0; it<fnTrials; it++){
2145
2146       rndphi = fRandom->Uniform(0, fPhiSize);
2147       rndeta = fRandom->Uniform(0, fEtaSize);
2148
2149       for(Int_t ip=0; ip<fnPhi; ip++){  //move radom position to each cell
2150          rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
2151          for(Int_t ie=0; ie<fnEta; ie++){
2152             rndetashift = rndeta + ie*fEtaSize - fEtaSize;
2153
2154             bIsInCone = 0; //tag if trial is in the jet cone
2155             for(Int_t ij=0; ij<njets; ij++){
2156                deta = jeta[ij] - rndetashift;
2157                dphi = RelativePhi(rndphishift,jphi[ij]);
2158                if((dphi*dphi + deta*deta) < jRsquared[ij]){
2159                   bIsInCone = 1;
2160                   break;
2161                }
2162             }
2163             if(!bIsInCone) nOutCone[ip][ie]++;
2164          }
2165       }
2166    }
2167
2168
2169    //Sum particle pT outside jets in cells
2170    Int_t npart = listPart->GetEntries();
2171    Int_t phicell,etacell;
2172    AliVParticle *part = NULL;
2173    for(Int_t ip=0; ip < npart; ip++){
2174
2175       part = (AliVParticle*) listPart->At(ip);
2176       if(!part) continue;
2177       if( TMath::Abs(part->Eta()) > fEtaSize) continue; //skip part outside +-0.7 in eta
2178
2179       bIsInCone = 0; //init
2180       for(Int_t ij=0; ij<njets; ij++){
2181          dphi = RelativePhi(jphi[ij], part->Phi());
2182          deta = jeta[ij] - part->Eta();
2183          if((dphi*dphi + deta*deta) < jRsquared[ij]){
2184             bIsInCone = 1;
2185             break;
2186          }
2187       }
2188       if(!bIsInCone){
2189          phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
2190          etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
2191          sumPtOutOfCone[phicell][etacell]+= part->Pt();
2192       }
2193    }
2194
2195    static Double_t rhoInCells[20];
2196    Double_t  relativeArea;
2197    Int_t  nCells=0;
2198    Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
2199    for(Int_t ip=0; ip<fnPhi; ip++){
2200       for(Int_t ie=0; ie<fnEta; ie++){
2201          relativeArea = nOutCone[ip][ie]/fnTrials;
2202          //cout<<"RA=  "<< relativeArea<<"  BA="<<bufferArea<<"    BPT="<< bufferPt <<endl;
2203
2204          bufferArea += relativeArea;
2205          bufferPt   += sumPtOutOfCone[ip][ie];
2206          if(bufferArea > fJetFreeAreaFrac){
2207             rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
2208
2209             if(mode==1) fhCellAreaToMedianGen->Fill(bufferArea*fCellArea);
2210             else fhCellAreaToMedian->Fill(bufferArea*fCellArea);
2211
2212             bufferArea = 0.0;
2213             bufferPt   = 0.0;
2214             nCells++;
2215          }
2216       }
2217    }
2218
2219
2220    if(nCells>0){
2221      rhoMedian = TMath::Median(nCells, rhoInCells);
2222      if(mode==1){ //mc data
2223         fhEntriesToMedianGen->Fill(nCells);
2224      }else{ //real data
2225         fhEntriesToMedian->Fill(nCells);
2226      }
2227    }
2228
2229 }
2230
2231 //__________________________________________________________________
2232 void AliAnalysisTaskJetCorePP::EstimateBgCone(TList *listJet, TList* listPart,  Double_t &rhoPerpCone){
2233
2234    rhoPerpCone = 0.0; //init
2235
2236    if(!listJet) return;
2237    Int_t njet  = listJet->GetEntries();
2238    Int_t npart = listPart->GetEntries();
2239    Double_t pTleading  = 0.0;
2240    Double_t phiLeading = 1000.;
2241    Double_t etaLeading = 1000.;
2242
2243    for(Int_t jsig=0; jsig < njet; jsig++){
2244       AliAODJet* signaljet = (AliAODJet*)(listJet->At(jsig));
2245       if(!signaljet) continue;
2246       if((signaljet->Eta()<fJetEtaMin) && (fJetEtaMax<signaljet->Eta())) continue; //acceptance cut
2247       if(signaljet->Pt() >= pTleading){ //replace leading and subleading jet
2248          pTleading  = signaljet->Pt();
2249          phiLeading = signaljet->Phi();
2250          etaLeading = signaljet->Eta();
2251       }
2252    }
2253
2254    Double_t phileftcone  = phiLeading + TMath::Pi()/2;
2255    Double_t phirightcone = phiLeading - TMath::Pi()/2;
2256    Double_t dp, de;
2257
2258    for(Int_t ip=0; ip < npart; ip++){
2259
2260       AliVParticle *part = (AliVParticle*) listPart->At(ip);
2261       if(!part){
2262          continue;
2263       }
2264
2265       dp = RelativePhi(phileftcone, part->Phi());
2266       de = etaLeading - part->Eta();
2267       if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2268
2269       dp = RelativePhi(phirightcone, part->Phi());
2270       if( sqrt(dp*dp + de*de)< fBgConeR) rhoPerpCone += part->Pt();
2271
2272    }
2273
2274    //normalize total pT by two times cone are
2275    rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fBgConeR*fBgConeR);
2276
2277 }
2278 //__________________________________________________________________
2279
2280 void AliAnalysisTaskJetCorePP::ReadTClonesArray(TString bname, TList *list){
2281    //Convert TClones array of jets to TList
2282
2283    if(!strlen(bname.Data())){
2284       AliError(Form("Jet branch %s not set.", bname.Data()));
2285       return;
2286    }
2287
2288    TClonesArray *array=0x0;
2289    if(fUseExchContainer){ //take input from exchange containers
2290       if(fJetBranchNameKine.Length()>0 && bname.CompareTo(fJetBranchNameKine.Data())==0){
2291          array = dynamic_cast<TClonesArray*>(GetInputData(1)); //connect exchange container slot 1
2292       }
2293       if(fJetBranchNameBgKine.Length()>0 && bname.CompareTo(fJetBranchNameBgKine.Data())==0){
2294         array = dynamic_cast<TClonesArray*>(GetInputData(2)); //connect exchange container slot 2
2295       }
2296    }else{ //take input from AOD
2297       if(fAODOut&&!array){
2298          array = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(bname.Data()));
2299       }
2300       if(fAODExtension&&!array){
2301          array = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(bname.Data()));
2302       }
2303       if(fAODIn&&!array){
2304          array = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(bname.Data()));
2305       }
2306    }
2307
2308    list->Clear(); //clear list beforehand
2309
2310    if(array){
2311       if(fDebug){
2312          Printf("########## %s: %d jets",bname.Data(), array->GetEntriesFast());
2313       }
2314
2315       for(Int_t iJet = 0; iJet < array->GetEntriesFast(); iJet++) {
2316          AliAODJet *jet = dynamic_cast<AliAODJet*>((*array)[iJet]);
2317          if (jet) list->Add(jet);
2318       }
2319    }
2320
2321    return;
2322 }