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