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