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 fJetBranchNameMC(""),
  81 fListJets(0x0),
  82 fListJetsGen(0x0),
  83 fNonStdFile(""),
  84 fSystem(0), //pp=0  pPb=1
  85 fJetParamR(0.4),
  86 fOfflineTrgMask(AliVEvent::kAny),
  87 fMinContribVtx(1),
  88 fVtxZMin(-10.0),
  89 fVtxZMax(10.0),
  90 fFilterMask(0),
  91 fCentMin(0.0),
  92 fCentMax(100.0),
  93 fJetEtaMin(-0.5),
  94 fJetEtaMax(0.5),
  95 fTriggerEtaCut(0.9),
  96 fTrackEtaCut(0.9),
  97 fTrackLowPtCut(0.15),
  98 fOutputList(0x0),
  99 fHistEvtSelection(0x0),
 100 fh2Ntriggers(0x0),
 101 fHJetSpec(0x0),
 102 fHJetDensity(0x0),
 103 fHJetDensityA4(0x0),
 104 fhJetPhi(0x0),
 105 fhTriggerPhi(0x0),
 106 fhJetEta(0x0),
 107 fhTriggerEta(0x0),
 108 fhVertexZ(0x0),
 109 fhVertexZAccept(0x0),
 110 fhContribVtx(0x0),
 111 fhContribVtxAccept(0x0),
 112 fhDphiTriggerJet(0x0),
 113 fhDphiTriggerJetAccept(0x0),
 114 fhCentrality(0x0),
 115 fhCentralityAccept(0x0),
 116 fHJetPtRaw(0x0),
 117 fHLeadingJetPtRaw(0x0),
 118 fHDphiVsJetPtAll(0x0),
 119 fhJetPtGenVsJetPtRec(0x0),
 120 fhJetPtGen(0x0),
 121 fh2NtriggersGen(0x0),
 122 fHJetSpecGen(0x0),
 123 fhPtTrkTruePrimRec(0x0),
 124 fhPtTrkTruePrimGen(0x0),
 125 fhPtTrkSecOrFakeRec(0x0),
 126 fIsMC(0),
 127 faGenIndex(0),
 128 faRecIndex(0),
 129 fkAcceptance(2.0*TMath::Pi()*1.8),
 130 fkDeltaPhiCut(TMath::Pi()-0.8),
 131 fh1Xsec(0x0),
 132 fh1Trials(0x0),
 133 fh1AvgTrials(0x0),
 134 fh1PtHard(0x0),
 135 fh1PtHardNoW(0x0),
 136 fh1PtHardTrials(0x0),
 137 fAvgTrials(1),
 138 fHardest(0),
 139 fEventNumberRangeLow(0),
 140 fEventNumberRangeHigh(99),
 141 fTriggerPtRangeLow(0.0),
 142 fTriggerPtRangeHigh(10000.0),
 143 fRandom(0x0)
 144 {
 145    // default Constructor
 146 }
 147
 148 //---------------------------------------------------------------------
 149
 150 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
 151 AliAnalysisTaskSE(name),
 152 fESD(0x0),
 153 fAODIn(0x0),
 154 fAODOut(0x0),
 155 fAODExtension(0x0),
 156 fJetBranchName(""),
 157 fJetBranchNameMC(""),
 158 fListJets(0x0),
 159 fListJetsGen(0x0),
 160 fNonStdFile(""),
 161 fSystem(0),  //pp=0   pPb=1
 162 fJetParamR(0.4),
 163 fOfflineTrgMask(AliVEvent::kAny),
 164 fMinContribVtx(1),
 165 fVtxZMin(-10.0),
 166 fVtxZMax(10.0),
 167 fFilterMask(0),
 168 fCentMin(0.0),
 169 fCentMax(100.0),
 170 fJetEtaMin(-0.5),
 171 fJetEtaMax(0.5),
 172 fTriggerEtaCut(0.9),
 173 fTrackEtaCut(0.9),
 174 fTrackLowPtCut(0.15),
 175 fOutputList(0x0),
 176 fHistEvtSelection(0x0),
 177 fh2Ntriggers(0x0),
 178 fHJetSpec(0x0),
 179 fHJetDensity(0x0),
 180 fHJetDensityA4(0x0),
 181 fhJetPhi(0x0),
 182 fhTriggerPhi(0x0),
 183 fhJetEta(0x0),
 184 fhTriggerEta(0x0),
 185 fhVertexZ(0x0),
 186 fhVertexZAccept(0x0),
 187 fhContribVtx(0x0),
 188 fhContribVtxAccept(0x0),
 189 fhDphiTriggerJet(0x0),
 190 fhDphiTriggerJetAccept(0x0),
 191 fhCentrality(0x0),
 192 fhCentralityAccept(0x0),
 193 fHJetPtRaw(0x0),
 194 fHLeadingJetPtRaw(0x0),
 195 fHDphiVsJetPtAll(0x0),
 196 fhJetPtGenVsJetPtRec(0x0),
 197 fhJetPtGen(0x0),
 198 fh2NtriggersGen(0x0),
 199 fHJetSpecGen(0x0),
 200 fhPtTrkTruePrimRec(0x0),
 201 fhPtTrkTruePrimGen(0x0),
 202 fhPtTrkSecOrFakeRec(0x0),
 203 fIsMC(0),
 204 faGenIndex(0),
 205 faRecIndex(0),
 206 fkAcceptance(2.0*TMath::Pi()*1.8),
 207 fkDeltaPhiCut(TMath::Pi()-0.8),
 208 fh1Xsec(0x0),
 209 fh1Trials(0x0),
 210 fh1AvgTrials(0x0),
 211 fh1PtHard(0x0),
 212 fh1PtHardNoW(0x0),
 213 fh1PtHardTrials(0x0),
 214 fAvgTrials(1),
 215 fHardest(0),
 216 fEventNumberRangeLow(0),
 217 fEventNumberRangeHigh(99),
 218 fTriggerPtRangeLow(0.0),
 219 fTriggerPtRangeHigh(10000.0),
 220 fRandom(0x0)
 221 {
 222 // Constructor
 223
 224    DefineOutput(1, TList::Class());
 225 }
 226
 227 //--------------------------------------------------------------
 228 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
 229 AliAnalysisTaskSE(a.GetName()),
 230 fESD(a.fESD),
 231 fAODIn(a.fAODIn),
 232 fAODOut(a.fAODOut),
 233 fAODExtension(a.fAODExtension),
 234 fJetBranchName(a.fJetBranchName),
 235 fJetBranchNameMC(a.fJetBranchNameMC),
 236 fListJets(a.fListJets),
 237 fListJetsGen(a.fListJetsGen),
 238 fNonStdFile(a.fNonStdFile),
 239 fSystem(a.fSystem),
 240 fJetParamR(a.fJetParamR),
 241 fOfflineTrgMask(a.fOfflineTrgMask),
 242 fMinContribVtx(a.fMinContribVtx),
 243 fVtxZMin(a.fVtxZMin),
 244 fVtxZMax(a.fVtxZMax),
 245 fFilterMask(a.fFilterMask),
 246 fCentMin(a.fCentMin),
 247 fCentMax(a.fCentMax),
 248 fJetEtaMin(a.fJetEtaMin),
 249 fJetEtaMax(a.fJetEtaMax),
 250 fTriggerEtaCut(a.fTriggerEtaCut),
 251 fTrackEtaCut(a.fTrackEtaCut),
 252 fTrackLowPtCut(a.fTrackLowPtCut),
 253 fOutputList(a.fOutputList),
 254 fHistEvtSelection(a.fHistEvtSelection),
 255 fh2Ntriggers(a.fh2Ntriggers),
 256 fHJetSpec(a.fHJetSpec),
 257 fHJetDensity(a.fHJetDensity),
 258 fHJetDensityA4(a.fHJetDensityA4),
 259 fhJetPhi(a.fhJetPhi),
 260 fhTriggerPhi(a.fhTriggerPhi),
 261 fhJetEta(a.fhJetEta),
 262 fhTriggerEta(a.fhTriggerEta),
 263 fhVertexZ(a.fhVertexZ),
 264 fhVertexZAccept(a.fhVertexZAccept),
 265 fhContribVtx(a.fhContribVtx),
 266 fhContribVtxAccept(a.fhContribVtxAccept),
 267 fhDphiTriggerJet(a.fhDphiTriggerJet),
 268 fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
 269 fhCentrality(a.fhCentrality),
 270 fhCentralityAccept(a.fhCentralityAccept),
 271 fHJetPtRaw(a.fHJetPtRaw),
 272 fHLeadingJetPtRaw(a.fHLeadingJetPtRaw),
 273 fHDphiVsJetPtAll(a.fHDphiVsJetPtAll),
 274 fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec),
 275 fhJetPtGen(a.fhJetPtGen),
 276 fh2NtriggersGen(a.fh2NtriggersGen),
 277 fHJetSpecGen(a.fHJetSpecGen),
 278 fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec),
 279 fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen),
 280 fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec),
 281 fIsMC(a.fIsMC),
 282 faGenIndex(a.faGenIndex),
 283 faRecIndex(a.faRecIndex),
 284 fkAcceptance(a.fkAcceptance),
 285 fkDeltaPhiCut(a.fkDeltaPhiCut),
 286 fh1Xsec(a.fh1Xsec),
 287 fh1Trials(a.fh1Trials),
 288 fh1AvgTrials(a.fh1AvgTrials),
 289 fh1PtHard(a.fh1PtHard),
 290 fh1PtHardNoW(a.fh1PtHardNoW),
 291 fh1PtHardTrials(a.fh1PtHardTrials),
 292 fAvgTrials(a.fAvgTrials),
 293 fHardest(a.fHardest),
 294 fEventNumberRangeLow(a.fEventNumberRangeLow),
 295 fEventNumberRangeHigh(a.fEventNumberRangeHigh),
 296 fTriggerPtRangeLow(a.fTriggerPtRangeLow),
 297 fTriggerPtRangeHigh(a.fTriggerPtRangeHigh),
 298 fRandom(a.fRandom)
 299 {
 300    //Copy Constructor
 301    fRandom->SetSeed(0);
 302 }
 303 //--------------------------------------------------------------
 304
 305 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
 306   // assignment operator
 307   this->~AliAnalysisTaskJetCorePP();
 308   new(this) AliAnalysisTaskJetCorePP(a);
 309   return *this;
 310 }
 311 //--------------------------------------------------------------
 312
 313 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
 314 {
 315    //Destructor
 316    delete fListJets;
 317    delete fListJetsGen;
 318    delete fOutputList; // ????
 319    delete fRandom;
 320 }
 321
 322 //--------------------------------------------------------------
 323
 324
 325 Bool_t AliAnalysisTaskJetCorePP::Notify()
 326 {
 327    //Implemented Notify() to read the cross sections
 328    //and number of trials from pyxsec.root
 329    //inspired by AliAnalysisTaskJetSpectrum2::Notify()
 330    if(!fIsMC) return kFALSE;
 331
 332    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 333    if(!fESD){
 334       if(fDebug>1) AliError("ESD not available");
 335       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 336    }
 337
 338    fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 339
 340
 341    if(fNonStdFile.Length()!=0){
 342       // case that we have an AOD extension we can fetch the jets from the extended output
 343       AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 344       fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
 345       if(!fAODExtension){
 346          if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
 347       }
 348    }
 349
 350    TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 351    Float_t xsection = 0;
 352    Float_t ftrials  = 1;
 353
 354    fAvgTrials = 1;
 355    if(tree){
 356       TFile *curfile = tree->GetCurrentFile();
 357       if(!curfile) {
 358          Error("Notify","No current file");
 359          return kFALSE;
 360       }
 361       if(!fh1Xsec || !fh1Trials){
 362          Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
 363          return kFALSE;
 364       }
 365       AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
 366       fh1Xsec->Fill("<#sigma>",xsection);
 367       // construct a poor man average trials
 368       Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
 369       if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries;
 370       fh1Trials->Fill("#sum{ntrials}",ftrials);
 371    }
 372
 373    if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
 374
 375    return kTRUE;
 376 }
 377 //--------------------------------------------------------------
 378
 379 void AliAnalysisTaskJetCorePP::Init()
 380 {
 381    // check for jet branches
 382    if(!strlen(fJetBranchName.Data())){
 383       AliError("Jet branch name not set.");
 384    }
 385 }
 386
 387 //--------------------------------------------------------------
 388
 389 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
 390 {
 391   // Create histograms
 392    // Called once
 393    fListJets = new TList();  //reconstructed level
 394
 395    fIsMC = (fJetBranchNameMC.Length()>0) ? kTRUE : kFALSE;
 396
 397    fRandom = new TRandom3(0);
 398
 399    if(fIsMC) fListJetsGen = new TList(); //generator level
 400    OpenFile(1);
 401    if(!fOutputList) fOutputList = new TList();
 402    fOutputList->SetOwner(kTRUE);
 403
 404    Bool_t oldStatus = TH1::AddDirectoryStatus();
 405    TH1::AddDirectory(kFALSE);
 406
 407    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 408    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 409    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 410    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 411    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 412    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 413    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
 414
 415    fOutputList->Add(fHistEvtSelection);
 416
 417    Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
 418    //trigger pt spectrum (reconstructed)
 419    fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
 420                              nBinsCentrality,0.0,100.0,50,0.0,50.0);
 421    fOutputList->Add(fh2Ntriggers);
 422
 423    //Centrality, A, pTjet, pTtrigg, dphi
 424    const Int_t dimSpec   = 5;
 425    const Int_t nBinsSpec[dimSpec]     = {nBinsCentrality, 100,   100,  50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))};
 426    const Double_t lowBinSpec[dimSpec] = {0.0,             0.0,     0, 0.0, fkDeltaPhiCut};
 427    const Double_t hiBinSpec[dimSpec]  = {100.0,           1.0, 200.0,50.0, TMath::Pi()};
 428    fHJetSpec = new THnSparseF("fHJetSpec",
 429                    "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]",
 430                    dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
 431    fOutputList->Add(fHJetSpec);
 432
 433    //------------------- HISTOS FOR DIAGNOSTIC ----------------------
 434    //Jet number density histos [Trk Mult, jet density, pT trigger]
 435    const Int_t    dimJetDens   = 3;
 436    const Int_t    nBinsJetDens[dimJetDens]  = {100,   100,   10};
 437    const Double_t lowBinJetDens[dimJetDens] = {0.0,   0.0,  0.0};
 438    const Double_t hiBinJetDens[dimJetDens]  = {500.0, 5.0, 50.0};
 439
 440    fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
 441                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
 442
 443    fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
 444                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
 445
 446    fOutputList->Add(fHJetDensity);
 447    fOutputList->Add(fHJetDensityA4);
 448
 449
 450    //inclusive azimuthal and pseudorapidity histograms
 451    fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
 452                         50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
 453    fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
 454                         50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
 455    fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
 456                         50,0, 100, 40,-0.9,0.9);
 457    fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
 458                         50, 0, 50, 40,-0.9,0.9);
 459
 460    fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
 461    fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
 462    fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
 463    fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
 464    fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
 465    fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
 466    fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
 467    fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
 468
 469    fOutputList->Add(fhJetPhi);
 470    fOutputList->Add(fhTriggerPhi);
 471    fOutputList->Add(fhJetEta);
 472    fOutputList->Add(fhTriggerEta);
 473    fOutputList->Add(fhVertexZ);
 474    fOutputList->Add(fhVertexZAccept);
 475    fOutputList->Add(fhContribVtx);
 476    fOutputList->Add(fhContribVtxAccept);
 477    fOutputList->Add(fhDphiTriggerJet);
 478    fOutputList->Add(fhDphiTriggerJetAccept);
 479    fOutputList->Add(fhCentrality);
 480    fOutputList->Add(fhCentralityAccept);
 481
 482    // raw spectra of INCLUSIVE jets
 483    //Centrality, pTjet, A
 484    const Int_t dimRaw   = 3;
 485    const Int_t nBinsRaw[dimRaw]     = {nBinsCentrality,  50,   100};
 486    const Double_t lowBinRaw[dimRaw] = {0.0,             0.0,   0.0};
 487    const Double_t hiBinRaw[dimRaw]  = {100.0,           100,   1.0};
 488    fHJetPtRaw = new THnSparseF("fHJetPtRaw",
 489                                 "Incl. jet spectrum [cent,pTjet,A]",
 490                                 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
 491    fOutputList->Add(fHJetPtRaw);
 492
 493    // raw spectra of LEADING jets
 494    //Centrality, pTjet, A
 495    fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
 496                                 "Leading jet spectrum [cent,pTjet,A]",
 497                                 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
 498    fOutputList->Add(fHLeadingJetPtRaw);
 499
 500    // Dphi versus pT jet
 501    //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
 502    const Int_t dimDp   = 4;
 503    const Int_t nBinsDp[dimDp]     = {nBinsCentrality,  50,     50,    50};
 504    const Double_t lowBinDp[dimDp] = {0.0,       -TMath::Pi(),   0.0,   0.0};
 505    const Double_t hiBinDp[dimDp]  = {100.0,      TMath::Pi(), 100.0, 100.0};
 506    fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
 507                                 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
 508                                 dimDp,nBinsDp,lowBinDp,hiBinDp);
 509    fOutputList->Add(fHDphiVsJetPtAll);
 510
 511
 512    //analyze MC generator level
 513    if(fIsMC){
 514       fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 100,0,200, 100,0,200);
 515       fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC jet pt spectrum versus reconstructed pt spectrum
 516
 517       fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",100,0,200); //MC generator jet pt spectrum
 518       fOutputList->Add(fhJetPtGen);
 519
 520       fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
 521       fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
 522       fOutputList->Add(fh2NtriggersGen);
 523
 524       //Centrality, A, pT, pTtrigg
 525       fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
 526       fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
 527       fOutputList->Add(fHJetSpecGen);
 528
 529       //track efficiency/contamination histograms  eta versus pT
 530       fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
 531       fOutputList->Add(fhPtTrkTruePrimRec);
 532
 533       fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
 534       fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
 535       fOutputList->Add(fhPtTrkTruePrimGen);
 536
 537       fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
 538       fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
 539       fOutputList->Add(fhPtTrkSecOrFakeRec);
 540    }
 541    //-------------------------------------
 542    //     pythia histograms
 543    const Int_t nBinPt = 150;
 544    Double_t binLimitsPt[nBinPt+1];
 545    for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 546       if(iPt == 0){
 547          binLimitsPt[iPt] = -50.0;
 548       }else{// 1.0
 549          binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
 550       }
 551    }
 552
 553    fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 554    fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 555    fOutputList->Add(fh1Xsec);
 556    fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 557    fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 558    fOutputList->Add(fh1Trials);
 559    fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
 560    fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
 561    fOutputList->Add(fh1AvgTrials);
 562    fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 563    fOutputList->Add(fh1PtHard);
 564    fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 565    fOutputList->Add(fh1PtHardNoW);
 566    fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 567    fOutputList->Add(fh1PtHardTrials);
 568
 569
 570    // =========== Switch on Sumw2 for all histos ===========
 571    for(Int_t i=0; i<fOutputList->GetEntries(); i++){
 572       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 573       if(h1){
 574          h1->Sumw2();
 575          continue;
 576       }
 577       THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 578       if(hn){
 579          hn->Sumw2();
 580       }
 581    }
 582    TH1::AddDirectory(oldStatus);
 583
 584    PostData(1, fOutputList);
 585 }
 586
 587 //--------------------------------------------------------------------
 588
 589 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
 590 {
 591    //Event loop
 592    Double_t eventW  = 1.0;
 593    Double_t ptHard  = 0.0;
 594    Double_t nTrials = 1.0; // Trials for MC trigger
 595    if(fIsMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
 596
 597    if(TMath::Abs((Float_t) fJetParamR)<0.00001){
 598       AliError("Cone radius is set to zero.");
 599       return;
 600    }
 601    if(!strlen(fJetBranchName.Data())){
 602       AliError("Jet branch name not set.");
 603       return;
 604    }
 605
 606    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 607    if(!fESD){
 608       if(fDebug>1) AliError("ESD not available");
 609       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 610    }
 611
 612    fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 613    AliAODEvent* aod = NULL;
 614    // take all other information from the aod we take the tracks from
 615    if(!aod){
 616       if(!fESD) aod = fAODIn;
 617       else      aod = fAODOut;
 618    }
 619
 620    if(fNonStdFile.Length()!=0){
 621       // case that we have an AOD extension we can fetch the jets from the extended output
 622       AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 623       fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
 624       if(!fAODExtension){
 625          if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
 626       }
 627    }
 628
 629    // ----------------- EVENT SELECTION  --------------------------
 630    fHistEvtSelection->Fill(1); // number of events before event selection
 631
 632    // physics selection
 633    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 634            ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 635
 636    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 637       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 638       fHistEvtSelection->Fill(2);
 639       PostData(1, fOutputList);
 640       return;
 641    }
 642
 643    //check AOD pointer
 644    if(!aod){
 645       if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
 646       fHistEvtSelection->Fill(3);
 647       PostData(1, fOutputList);
 648       return;
 649    }
 650
 651    // vertex selection
 652    AliAODVertex* primVtx = aod->GetPrimaryVertex();
 653
 654    if(!primVtx){
 655       if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
 656       fHistEvtSelection->Fill(3);
 657       PostData(1, fOutputList);
 658       return;
 659    }
 660
 661    Int_t nTracksPrim = primVtx->GetNContributors();
 662    Float_t vtxz = primVtx->GetZ();
 663    //Input events
 664    fhContribVtx->Fill(nTracksPrim);
 665    if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
 666
 667    if((nTracksPrim < fMinContribVtx) ||
 668       (vtxz < fVtxZMin) ||
 669       (vtxz > fVtxZMax)){
 670       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
 671                          (char*)__FILE__,__LINE__,vtxz);
 672       fHistEvtSelection->Fill(3);
 673       PostData(1, fOutputList);
 674       return;
 675    }else{
 676       //Accepted events
 677       fhContribVtxAccept->Fill(nTracksPrim);
 678       fhVertexZAccept->Fill(vtxz);
 679    }
 680    //FK// No event class selection imposed
 681    // event class selection (from jet helper task)
 682    //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 683    //if(fDebug) Printf("Event class %d", eventClass);
 684    //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 685    //   fHistEvtSelection->Fill(4);
 686    //   PostData(1, fOutputList);
 687    //   return;
 688    //}
 689
 690    //------------------ CENTRALITY SELECTION ---------------
 691    AliCentrality *cent = 0x0;
 692    Double_t centValue  = 0.0;
 693    if(fSystem){  //fSystem=0 for pp,   fSystem=1 for pPb
 694       if(fESD){
 695          cent = fESD->GetCentrality();
 696          if(cent) centValue = cent->GetCentralityPercentile("V0M");
 697       }else{
 698          centValue = aod->GetHeader()->GetCentrality();
 699       }
 700       if(fDebug) printf("centrality: %f\n", centValue);
 701       //Input events
 702       fhCentrality->Fill(centValue);
 703
 704       if(centValue < fCentMin || centValue > fCentMax){
 705          fHistEvtSelection->Fill(4);
 706          PostData(1, fOutputList);
 707          return;
 708       }else{
 709          //Accepted events
 710          fhCentralityAccept->Fill( centValue );
 711       }
 712    }
 713
 714    //-----------------select disjunct event subsamples ----------------
 715    Int_t eventnum   = aod->GetHeader()->GetEventNumberESDFile();
 716    Int_t lastdigit = eventnum % 10;
 717    if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
 718       fHistEvtSelection->Fill(5);
 719       PostData(1, fOutputList);
 720       return;
 721    }
 722
 723    if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
 724    fHistEvtSelection->Fill(0); // accepted events
 725    // ------------------- end event selection --------------------
 726
 727
 728    // fetch RECONSTRUCTED jets
 729    TClonesArray *aodJets = 0x0;
 730
 731    if(fAODOut && !aodJets){
 732       aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
 733    }
 734    if(fAODExtension && !aodJets){
 735       aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
 736    }
 737    if(fAODIn && !aodJets){
 738       aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
 739    }
 740
 741    //--------- Fill list of RECONSTRUCTED jets --------------
 742    fListJets->Clear();
 743    if(aodJets){
 744       if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(), aodJets->GetEntriesFast());
 745       for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
 746          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
 747          if (jet) fListJets->Add(jet);
 748       }
 749    }
 750
 751    //--------- Fill list of MC GENERATOR LEVEL jets --------------
 752    TList particleListGen; //list of tracks in MC
 753
 754    if(fIsMC){ //analyze generator level MC
 755       // fetch MC generator level jets
 756       TClonesArray *aodGenJets = NULL;
 757
 758       if(fAODOut&&!aodGenJets){
 759          aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMC.Data()));
 760       }
 761       if(fAODExtension&&!aodGenJets){
 762          aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMC.Data()));
 763       }
 764       if(fAODIn&&!aodGenJets){
 765          aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMC.Data()));
 766       }
 767
 768       if(!aodGenJets){
 769          Printf("%s:%d no generated Jet array with name %s in AOD",
 770                  (char*)__FILE__,__LINE__, fJetBranchNameMC.Data());
 771          PostData(1, fOutputList);
 772          return;
 773       }
 774
 775       fListJetsGen->Clear();
 776
 777       //serarch for charged trigger at the MC generator level
 778       Int_t    indexTriggGen = -1;
 779       Double_t ptTriggGen    = -1;
 780       Int_t    iCounterGen   =  0;
 781       Int_t    triggersMC[200];//list of trigger candidates
 782       Int_t    ntriggersMC   = 0; //index in triggers array
 783
 784       if(fESD){//ESD input
 785
 786          AliMCEvent* mcEvent = MCEvent();
 787          if(!mcEvent){
 788             PostData(1, fOutputList);
 789             return;
 790          }
 791
 792          AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 793          if(pythiaGenHeader){
 794             nTrials = pythiaGenHeader->Trials();
 795             ptHard  = pythiaGenHeader->GetPtHard();
 796             fh1PtHard->Fill(ptHard,eventW);
 797             fh1PtHardNoW->Fill(ptHard,1);
 798             fh1PtHardTrials->Fill(ptHard,nTrials);
 799          }
 800
 801          for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
 802             if(!mcEvent->IsPhysicalPrimary(it)) continue;
 803             AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
 804             if(!part) continue;
 805             if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
 806
 807                if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
 808                   if(indexTriggGen > -1){//trigger candidater was found
 809                      triggersMC[ntriggersMC] = indexTriggGen;
 810                      ntriggersMC++;
 811                   }
 812                }
 813
 814                iCounterGen++;
 815             }
 816          }
 817       }else{  //AOD input
 818          if(!fAODIn){
 819             PostData(1, fOutputList);
 820             return;
 821          }
 822          TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
 823          if(!tca){
 824             PostData(1, fOutputList);
 825             return;
 826          }
 827          for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
 828             AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
 829             if(!part) continue;
 830             if(!part->IsPhysicalPrimary()) continue;
 831             if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
 832
 833                if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
 834                   if(indexTriggGen > -1){ //trigger candidater was found
 835                      triggersMC[ntriggersMC] = indexTriggGen;
 836                      ntriggersMC++;
 837                   }
 838                }
 839
 840                iCounterGen++;
 841             }
 842          }
 843       }
 844
 845       if(fHardest==0){ //single inclusive trigger
 846          if(ntriggersMC>0){ //there is at least one trigger
 847             Int_t rnd     = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
 848             indexTriggGen = triggersMC[rnd];
 849          }else{
 850             indexTriggGen = -1; //trigger not found
 851          }
 852       }
 853
 854       //================== Fill jet list ===================
 855       if(aodGenJets){
 856          if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMC.Data(), aodGenJets->GetEntriesFast());
 857
 858          for(Int_t igJet = 0; igJet < aodGenJets->GetEntriesFast(); igJet++) {
 859             AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJets)[igJet]);
 860             if(jetGen) fListJetsGen->Add(jetGen);
 861          }
 862       }
 863
 864       //============  Generator trigger+jet ==================
 865       Int_t ilowGen  = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
 866       Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
 867
 868       for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
 869          indexTriggGen = igen; //trigger hadron
 870
 871          if(indexTriggGen == -1) continue;
 872          AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
 873          if(!triggerGen) continue;
 874
 875          if(fHardest >= 2){
 876             if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6
 877          }
 878          if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
 879
 880          ptTriggGen = triggerGen->Pt(); //count triggers
 881          fh2NtriggersGen->Fill(centValue, ptTriggGen);
 882
 883          //Count jets and trigger-jet pairs at MC  generator level
 884          if(!aodGenJets) continue;
 885          for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
 886             AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
 887             if(!jet) continue;
 888             Double_t etaJetGen = jet->Eta();
 889
 890             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
 891
 892                Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
 893                if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
 894
 895                //Centrality, A, pT, pTtrigg
 896                Double_t fillspecgen[] = { centValue,
 897                                           jet->EffectiveAreaCharged(),
 898                                           jet->Pt(),
 899                                           ptTriggGen,
 900                                           TMath::Abs((Double_t) dphi)
 901                                         };
 902               fHJetSpecGen->Fill(fillspecgen);
 903             }//back to back jet-trigger pair
 904          }//jet loop
 905       }//trigger loop
 906
 907
 908       //================ RESPONSE MATRIX ===============
 909       if(aodGenJets){
 910          //Count jets and trigger-jet pairs at MC  generator level
 911          for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
 912             AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
 913             if(!jet) continue;
 914             Double_t etaJetGen = jet->Eta();
 915             Double_t ptJetGen  = jet->Pt();
 916
 917             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
 918                fhJetPtGen->Fill(ptJetGen); // generator level pt spectum of jets response mx normalization
 919             }
 920          }
 921          if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
 922             Int_t ng = (Int_t) fListJetsGen->GetEntries();
 923             Int_t nr = (Int_t) fListJets->GetEntries();
 924
 925             //Find closest MC generator - reconstructed jets
 926             if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
 927             if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
 928
 929             if(fDebug){
 930                Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
 931                Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
 932             }
 933             //matching of MC genrator level and reconstructed jets
 934             AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
 935
 936             // Fill response matrix
 937             for(Int_t ir = 0; ir < nr; ir++){
 938                AliAODJet *recJet  = (AliAODJet*) fListJets->At(ir);
 939                Double_t etaJetRec = recJet->Eta();
 940                Double_t ptJetRec  = recJet->Pt();
 941                //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
 942
 943                if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
 944                   Int_t ig = faGenIndex[ir]; //associated generator level jet
 945                   if(ig >= 0 && ig < ng){
 946                      if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 947                      AliAODJet *genJet  = (AliAODJet*) fListJetsGen->At(ig);
 948                      Double_t ptJetGen  = genJet->Pt();
 949                      Double_t etaJetGen = genJet->Eta();
 950
 951                   //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
 952                      if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
 953                         fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
 954                      }
 955                   }//ig>=0
 956                }//rec jet in eta acceptance
 957             }//loop over reconstructed jets
 958          }// # of  rec jets >0
 959       }//pointer MC generator jets
 960    } //analyze generator level MC
 961
 962    //=============  RECONSTRUCTED INCLUSIVE JETS ===============
 963
 964    Double_t etaJet  = 0.0;
 965    Double_t pTJet   = 0.0;
 966    Double_t areaJet = 0.0;
 967    Double_t phiJet  = 0.0;
 968    Int_t indexLeadingJet     = -1;
 969    Double_t pTLeadingJet     = -10.0;
 970    Double_t areaLeadingJet   = -10.0;
 971
 972    for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
 973       AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
 974       if(!jet) continue;
 975       etaJet  = jet->Eta();
 976       phiJet  = jet->Phi();
 977       pTJet   = jet->Pt();
 978       if(pTJet==0) continue;
 979
 980       if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
 981       areaJet = jet->EffectiveAreaCharged();
 982
 983       //Jet Diagnostics---------------------------------
 984       fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
 985       fhJetEta->Fill(pTJet, etaJet);
 986       //search for leading jet
 987       if(pTJet > pTLeadingJet){
 988          indexLeadingJet  = ij;
 989          pTLeadingJet     = pTJet;
 990          areaLeadingJet   = areaJet;
 991       }
 992
 993       // raw spectra of INCLUSIVE jets
 994       //Centrality, pTjet, A
 995       Double_t fillraw[] = { centValue,
 996                              pTJet,
 997                              areaJet
 998                            };
 999       fHJetPtRaw->Fill(fillraw);
1000    }
1001
1002    if(indexLeadingJet > -1){
1003       // raw spectra of LEADING jets
1004       //Centrality, pTjet,  A
1005       Double_t fillleading[] = { centValue,
1006                                  pTLeadingJet,
1007                                  areaLeadingJet
1008                                };
1009       fHLeadingJetPtRaw->Fill(fillleading);
1010    }
1011
1012
1013    // ===============  RECONSTRUCTED TRIGGER-JET PAIRS ================
1014    //Find Hadron trigger
1015    TList particleList; //list of tracks
1016    Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1017    if(fIsMC) FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1018
1019    //set ranges of the trigger loop
1020    Int_t ilow  = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1021    Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
1022
1023    for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1024       indexTrigg = itrk; //trigger hadron with pT >10 GeV
1025
1026       if(indexTrigg < 0) continue;
1027
1028       AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1029       if(!triggerHadron){
1030          PostData(1, fOutputList);
1031          return;
1032       }
1033       if(fHardest >= 2){
1034          if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
1035       }
1036       if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
1037
1038       //Fill trigger histograms
1039       fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1040       fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1041       fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
1042
1043
1044       //---------- make trigger-jet pairs ---------
1045       Int_t injet4     = 0;
1046       Int_t injet      = 0;
1047
1048       for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1049          AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1050          if(!jet) continue;
1051          etaJet  = jet->Eta();
1052          phiJet  = jet->Phi();
1053          pTJet   = jet->Pt();
1054          if(pTJet==0) continue;
1055
1056          if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1057          areaJet = jet->EffectiveAreaCharged();
1058          if(areaJet >= 0.07) injet++;
1059          if(areaJet >= 0.4)  injet4++;
1060
1061          Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1062          fhDphiTriggerJet->Fill(dphi); //Input
1063
1064          //Dphi versus jet pT
1065          //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
1066          Double_t filldp[] = { centValue,
1067                                dphi,
1068                                pTJet,
1069                                triggerHadron->Pt()
1070                              };
1071          fHDphiVsJetPtAll->Fill(filldp);
1072
1073          // Select back to back trigger - jet pairs
1074          if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1075          fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1076
1077
1078          //Centrality, A, pTjet, pTtrigg
1079          Double_t fillspec[] = { centValue,
1080                                  areaJet,
1081                                  pTJet,
1082                                  triggerHadron->Pt(),
1083                                  TMath::Abs((Double_t) dphi)
1084                                };
1085          fHJetSpec->Fill(fillspec);
1086       }//jet loop
1087
1088       //Fill Jet Density In the Event A>0.07
1089       if(injet>0){
1090          Double_t filldens[]={ (Double_t) particleList.GetEntries(),
1091                             injet/fkAcceptance,
1092                             triggerHadron->Pt()
1093                           };
1094          fHJetDensity->Fill(filldens);
1095       }
1096
1097       //Fill Jet Density In the Event A>0.4
1098       if(injet4>0){
1099          Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
1100                                 injet4/fkAcceptance,
1101                                 triggerHadron->Pt()
1102                               };
1103          fHJetDensityA4->Fill(filldens4);
1104       }
1105    }
1106
1107
1108    PostData(1, fOutputList);
1109 }
1110
1111 //----------------------------------------------------------------------------
1112 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
1113 {
1114    // Draw result to the screen
1115    // Called once at the end of the query
1116
1117    if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
1118    if(!GetOutputData(1)) return;
1119 }
1120
1121
1122 //----------------------------------------------------------------------------
1123 Int_t  AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
1124    //Fill the list of accepted tracks (passed track cut)
1125    //return consecutive index of the hardest ch hadron in the list
1126    Int_t iCount        = 0;
1127    AliAODEvent *aodevt = NULL;
1128
1129    if(!fESD) aodevt = fAODIn;
1130    else      aodevt = fAODOut;
1131
1132    if(!aodevt) return -1;
1133
1134    Int_t    index = -1; //index of the highest particle in the list
1135    Double_t ptmax = -10;
1136    Int_t    triggers[200];
1137    Int_t    ntriggers = 0; //index in triggers array
1138
1139    for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){
1140       AliAODTrack *tr = aodevt->GetTrack(it);
1141
1142       //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
1143       if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
1144       if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
1145       if(tr->Pt() < fTrackLowPtCut) continue;
1146       list->Add(tr);
1147
1148       //Search trigger:
1149       if(fHardest>0){ //leading particle
1150          if(tr->Pt()>ptmax){
1151             ptmax = tr->Pt();
1152             index = iCount;
1153          }
1154       }
1155
1156       if(fHardest==0 && ntriggers<200){ //single inclusive
1157          if(fTriggerPtRangeLow <= tr->Pt() &&
1158             tr->Pt() < fTriggerPtRangeHigh){
1159             triggers[ntriggers] = iCount;
1160             ntriggers++;
1161          }
1162       }
1163
1164       iCount++;
1165    }
1166
1167    if(fHardest==0 && ntriggers>0){ //select random inclusive trigger
1168       Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1
1169       index     = triggers[rnd];
1170    }
1171
1172    return index;
1173 }
1174
1175 //----------------------------------------------------------------------------
1176 /*
1177 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
1178    //calculate sum of track pT in the cone perpendicular in phi to the jet
1179    //jetR = cone radius
1180    // jetPhi, jetEta = direction of the jet
1181    Int_t numberOfTrks = trkList->GetEntries();
1182    Double_t pTsum = 0.0;
1183    Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
1184    for(Int_t it=0; it<numberOfTrks; it++){
1185       AliVParticle *trk = (AliVParticle*) trkList->At(it);
1186       Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
1187       Double_t deta = trk->Eta()-jetEta;
1188
1189       if( (dphi*dphi + deta*deta)< (jetR*jetR)){
1190          pTsum += trk->Pt();
1191       }
1192    }
1193
1194    return pTsum;
1195 }
1196 */
1197 //----------------------------------------------------------------------------
1198
1199 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
1200    //Get relative azimuthal angle of two particles -pi to pi
1201    if      (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1202    else if (vphi > TMath::Pi())  vphi -= TMath::TwoPi();
1203
1204    if      (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1205    else if (mphi > TMath::Pi())  mphi -= TMath::TwoPi();
1206
1207    Double_t dphi = mphi - vphi;
1208    if      (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1209    else if (dphi > TMath::Pi())  dphi -= TMath::TwoPi();
1210
1211    return dphi;//dphi in [-Pi, Pi]
1212 }
1213
1214
1215 //----------------------------------------------------------------------------
1216 Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){
1217    //fill track efficiency denominator
1218    if(trk->Pt() < fTrackLowPtCut) return kFALSE;
1219    if(trk->Charge()==0)        return kFALSE;
1220    if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE;
1221    trkList->Add(trk);
1222    fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator
1223
1224    //Search trigger:
1225    if(fHardest>0){ //leading particle
1226       if(ptLeading < trk->Pt()){
1227          index      = counter;
1228          ptLeading  = trk->Pt();
1229       }
1230    }
1231
1232    if(fHardest==0){ //single inclusive
1233       index = -1;
1234       if(fTriggerPtRangeLow <= trk->Pt() &&
1235             trk->Pt() < fTriggerPtRangeHigh){
1236          index      = counter;
1237       }
1238    }
1239
1240    return kTRUE;
1241 }
1242
1243 //----------------------------------------------------------------------------
1244 void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){
1245    //fill track efficiency numerator
1246    if(!recList) return;
1247    if(!genList) return;
1248    Int_t nRec = recList->GetEntries();
1249    Int_t nGen = genList->GetEntries();
1250    for(Int_t ir=0; ir<nRec; ir++){
1251       AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir);
1252       if(!trkRec) continue;
1253       Int_t recLabel = TMath::Abs(trkRec->GetLabel());
1254       Bool_t matched = kFALSE;
1255
1256       for(Int_t ig=0; ig<nGen; ig++){
1257         AliVParticle *trkGen = (AliVParticle*) genList->At(ig);
1258         if(!trkGen) continue;
1259         Int_t genLabel = TMath::Abs(trkGen->GetLabel());
1260         if(recLabel==genLabel){
1261           fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta());
1262           matched = kTRUE;
1263           break;
1264         }
1265       }
1266       if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta());
1267    }
1268
1269    return;
1270 }
1271
1272