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