PWGJE/AliAnalysisTaskJetCorePP.cxx

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