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