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