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