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.6),
 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.6),
 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
 390 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
 391 {
 392
 393   // Create histograms
 394    // Called once
 395    fListJets = new TList();  //reconstructed level
 396
 397    fIsMC = (fJetBranchNameMC.Length()>0) ? kTRUE : kFALSE;
 398
 399    fRandom = new TRandom3(0);
 400
 401    if(fIsMC) fListJetsGen = new TList(); //generator level
 402    OpenFile(1);
 403    if(!fOutputList) fOutputList = new TList();
 404    fOutputList->SetOwner(kTRUE);
 405
 406    Bool_t oldStatus = TH1::AddDirectoryStatus();
 407    TH1::AddDirectory(kFALSE);
 408
 409    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 410    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 411    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 412    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 413    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 414    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 415    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)");
 416
 417    fOutputList->Add(fHistEvtSelection);
 418
 419    Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
 420    //trigger pt spectrum (reconstructed)
 421    fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
 422                              nBinsCentrality,0.0,100.0,50,0.0,50.0);
 423    fOutputList->Add(fh2Ntriggers);
 424
 425    //Centrality, A, pTjet, pTtrigg
 426    const Int_t dimSpec   = 4;
 427    const Int_t nBinsSpec[dimSpec]     = {nBinsCentrality, 100,   100,  50};
 428    const Double_t lowBinSpec[dimSpec] = {0.0,             0.0,     0, 0.0};
 429    const Double_t hiBinSpec[dimSpec]  = {100.0,           1.0, 200.0,50.0};
 430    fHJetSpec = new THnSparseF("fHJetSpec",
 431                    "Recoil jet spectrum [cent,A,pTjet,pTtrig]",
 432                    dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
 433    fOutputList->Add(fHJetSpec);
 434
 435    //------------------- HISTOS FOR DIAGNOSTIC ----------------------
 436    //Jet number density histos [Trk Mult, jet density, pT trigger]
 437    const Int_t    dimJetDens   = 3;
 438    const Int_t    nBinsJetDens[dimJetDens]  = {100,   100,   10};
 439    const Double_t lowBinJetDens[dimJetDens] = {0.0,   0.0,  0.0};
 440    const Double_t hiBinJetDens[dimJetDens]  = {500.0, 5.0, 50.0};
 441
 442    fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
 443                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
 444
 445    fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
 446                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
 447
 448    fOutputList->Add(fHJetDensity);
 449    fOutputList->Add(fHJetDensityA4);
 450
 451
 452    //inclusive azimuthal and pseudorapidity histograms
 453    fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet",
 454                         50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
 455    fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg",
 456                         25, 0, 50, 50,-TMath::Pi(),TMath::Pi());
 457    fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet",
 458                         50,0, 100, 40,-0.9,0.9);
 459    fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg",
 460                         25, 0, 50, 40,-0.9,0.9);
 461
 462    fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
 463    fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
 464    fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
 465    fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
 466    fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
 467    fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
 468    fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
 469    fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
 470
 471    fOutputList->Add(fhJetPhi);
 472    fOutputList->Add(fhTriggerPhi);
 473    fOutputList->Add(fhJetEta);
 474    fOutputList->Add(fhTriggerEta);
 475    fOutputList->Add(fhVertexZ);
 476    fOutputList->Add(fhVertexZAccept);
 477    fOutputList->Add(fhContribVtx);
 478    fOutputList->Add(fhContribVtxAccept);
 479    fOutputList->Add(fhDphiTriggerJet);
 480    fOutputList->Add(fhDphiTriggerJetAccept);
 481    fOutputList->Add(fhCentrality);
 482    fOutputList->Add(fhCentralityAccept);
 483
 484    // raw spectra of INCLUSIVE jets
 485    //Centrality, pTjet, A
 486    const Int_t dimRaw   = 3;
 487    const Int_t nBinsRaw[dimRaw]     = {nBinsCentrality,  50,   100};
 488    const Double_t lowBinRaw[dimRaw] = {0.0,             0.0,   0.0};
 489    const Double_t hiBinRaw[dimRaw]  = {100.0,           100,   1.0};
 490    fHJetPtRaw = new THnSparseF("fHJetPtRaw",
 491                                 "Incl. jet spectrum [cent,pTjet,A]",
 492                                 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
 493    fOutputList->Add(fHJetPtRaw);
 494
 495    // raw spectra of LEADING jets
 496    //Centrality, pTjet, A
 497    fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw",
 498                                 "Leading jet spectrum [cent,pTjet,A]",
 499                                 dimRaw,nBinsRaw,lowBinRaw,hiBinRaw);
 500    fOutputList->Add(fHLeadingJetPtRaw);
 501
 502    // Dphi versus pT jet
 503    //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
 504    const Int_t dimDp   = 4;
 505    const Int_t nBinsDp[dimDp]     = {nBinsCentrality,  50,     50,    50};
 506    const Double_t lowBinDp[dimDp] = {0.0,       -TMath::Pi(),   0.0,   0.0};
 507    const Double_t hiBinDp[dimDp]  = {100.0,      TMath::Pi(), 100.0, 100.0};
 508    fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll",
 509                                 "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]",
 510                                 dimDp,nBinsDp,lowBinDp,hiBinDp);
 511    fOutputList->Add(fHDphiVsJetPtAll);
 512
 513
 514    //analyze MC generator level
 515    if(fIsMC){
 516       fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 100,0,200, 100,0,200);
 517       fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC jet pt spectrum versus reconstructed pt spectrum
 518
 519       fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",100,0,200); //MC generator jet pt spectrum
 520       fOutputList->Add(fhJetPtGen);
 521
 522       fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen");
 523       fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle()));
 524       fOutputList->Add(fh2NtriggersGen);
 525
 526       //Centrality, A, pT, pTtrigg
 527       fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen");
 528       fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle()));
 529       fOutputList->Add(fHJetSpecGen);
 530
 531       //track efficiency/contamination histograms  eta versus pT
 532       fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9);
 533       fOutputList->Add(fhPtTrkTruePrimRec);
 534
 535       fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen");
 536       fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen");
 537       fOutputList->Add(fhPtTrkTruePrimGen);
 538
 539       fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec");
 540       fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec");
 541       fOutputList->Add(fhPtTrkSecOrFakeRec);
 542    }
 543    //-------------------------------------
 544    //     pythia histograms
 545    const Int_t nBinPt = 150;
 546    Double_t binLimitsPt[nBinPt+1];
 547    for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
 548       if(iPt == 0){
 549          binLimitsPt[iPt] = -50.0;
 550       }else{// 1.0
 551          binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.;
 552       }
 553    }
 554
 555    fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
 556    fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
 557    fOutputList->Add(fh1Xsec);
 558    fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
 559    fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
 560    fOutputList->Add(fh1Trials);
 561    fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
 562    fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
 563    fOutputList->Add(fh1AvgTrials);
 564    fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 565    fOutputList->Add(fh1PtHard);
 566    fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 567    fOutputList->Add(fh1PtHardNoW);
 568    fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 569    fOutputList->Add(fh1PtHardTrials);
 570
 571
 572    // =========== Switch on Sumw2 for all histos ===========
 573    for(Int_t i=0; i<fOutputList->GetEntries(); i++){
 574       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 575       if(h1){
 576          h1->Sumw2();
 577          continue;
 578       }
 579       THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 580       if(hn){
 581          hn->Sumw2();
 582       }
 583    }
 584    TH1::AddDirectory(oldStatus);
 585
 586    PostData(1, fOutputList);
 587 }
 588
 589 //--------------------------------------------------------------------
 590
 591 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
 592 {
 593    //Event loop
 594    Double_t eventW  = 1.0;
 595    Double_t ptHard  = 0.0;
 596    Double_t nTrials = 1.0; // Trials for MC trigger
 597    if(fIsMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
 598
 599    if(TMath::Abs((Float_t) fJetParamR)<0.00001){
 600       AliError("Cone radius is set to zero.");
 601       return;
 602    }
 603    if(!strlen(fJetBranchName.Data())){
 604       AliError("Jet branch name not set.");
 605       return;
 606    }
 607
 608    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 609    if(!fESD){
 610       if(fDebug>1) AliError("ESD not available");
 611       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 612    }
 613
 614    fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 615    AliAODEvent* aod = NULL;
 616    // take all other information from the aod we take the tracks from
 617    if(!aod){
 618       if(!fESD) aod = fAODIn;
 619       else      aod = fAODOut;
 620    }
 621
 622    if(fNonStdFile.Length()!=0){
 623       // case that we have an AOD extension we can fetch the jets from the extended output
 624       AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 625       fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
 626       if(!fAODExtension){
 627          if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
 628       }
 629    }
 630
 631    // ----------------- EVENT SELECTION  --------------------------
 632    fHistEvtSelection->Fill(1); // number of events before event selection
 633
 634    // physics selection
 635    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 636            ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 637
 638    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 639       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 640       fHistEvtSelection->Fill(2);
 641       PostData(1, fOutputList);
 642       return;
 643    }
 644
 645    //check AOD pointer
 646    if(!aod){
 647       if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
 648       fHistEvtSelection->Fill(3);
 649       PostData(1, fOutputList);
 650       return;
 651    }
 652
 653    // vertex selection
 654    AliAODVertex* primVtx = aod->GetPrimaryVertex();
 655
 656    if(!primVtx){
 657       if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
 658       fHistEvtSelection->Fill(3);
 659       PostData(1, fOutputList);
 660       return;
 661    }
 662
 663    Int_t nTracksPrim = primVtx->GetNContributors();
 664    Float_t vtxz = primVtx->GetZ();
 665    //Input events
 666    fhContribVtx->Fill(nTracksPrim);
 667    if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
 668
 669    if((nTracksPrim < fMinContribVtx) ||
 670       (vtxz < fVtxZMin) ||
 671       (vtxz > fVtxZMax)){
 672       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
 673                          (char*)__FILE__,__LINE__,vtxz);
 674       fHistEvtSelection->Fill(3);
 675       PostData(1, fOutputList);
 676       return;
 677    }else{
 678       //Accepted events
 679       fhContribVtxAccept->Fill(nTracksPrim);
 680       fhVertexZAccept->Fill(vtxz);
 681    }
 682    //FK// No event class selection imposed
 683    // event class selection (from jet helper task)
 684    //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 685    //if(fDebug) Printf("Event class %d", eventClass);
 686    //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 687    //   fHistEvtSelection->Fill(4);
 688    //   PostData(1, fOutputList);
 689    //   return;
 690    //}
 691
 692    //------------------ CENTRALITY SELECTION ---------------
 693    AliCentrality *cent = 0x0;
 694    Double_t centValue  = 0.0;
 695    if(fSystem){  //fSystem=0 for pp,   fSystem=1 for pPb
 696       if(fESD){
 697          cent = fESD->GetCentrality();
 698          if(cent) centValue = cent->GetCentralityPercentile("V0M");
 699       }else{
 700          centValue = aod->GetHeader()->GetCentrality();
 701       }
 702       if(fDebug) printf("centrality: %f\n", centValue);
 703       //Input events
 704       fhCentrality->Fill(centValue);
 705
 706       if(centValue < fCentMin || centValue > fCentMax){
 707          fHistEvtSelection->Fill(4);
 708          PostData(1, fOutputList);
 709          return;
 710       }else{
 711          //Accepted events
 712          fhCentralityAccept->Fill( centValue );
 713       }
 714    }
 715
 716    //-----------------select disjunct event subsamples ----------------
 717    Int_t eventnum   = aod->GetHeader()->GetEventNumberESDFile();
 718    Int_t lastdigit = eventnum % 10;
 719    if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){
 720       fHistEvtSelection->Fill(5);
 721       PostData(1, fOutputList);
 722       return;
 723    }
 724
 725    if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
 726    fHistEvtSelection->Fill(0); // accepted events
 727    // ------------------- end event selection --------------------
 728
 729
 730    // fetch RECONSTRUCTED jets
 731    TClonesArray *aodJets = 0x0;
 732
 733    if(fAODOut && !aodJets){
 734       aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
 735    }
 736    if(fAODExtension && !aodJets){
 737       aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
 738    }
 739    if(fAODIn && !aodJets){
 740       aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
 741    }
 742
 743    //--------- Fill list of RECONSTRUCTED jets --------------
 744    fListJets->Clear();
 745    if(aodJets){
 746       if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(), aodJets->GetEntriesFast());
 747       for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
 748          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
 749          if (jet) fListJets->Add(jet);
 750       }
 751    }
 752
 753    //--------- Fill list of MC GENERATOR LEVEL jets --------------
 754    TList particleListGen; //list of tracks in MC
 755
 756    if(fIsMC){ //analyze generator level MC
 757       // fetch MC generator level jets
 758       TClonesArray *aodGenJets = NULL;
 759
 760       if(fAODOut&&!aodGenJets){
 761          aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMC.Data()));
 762       }
 763       if(fAODExtension&&!aodGenJets){
 764          aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMC.Data()));
 765       }
 766       if(fAODIn&&!aodGenJets){
 767          aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMC.Data()));
 768       }
 769
 770       if(!aodGenJets){
 771          Printf("%s:%d no generated Jet array with name %s in AOD",
 772                  (char*)__FILE__,__LINE__, fJetBranchNameMC.Data());
 773          PostData(1, fOutputList);
 774          return;
 775       }
 776
 777       fListJetsGen->Clear();
 778
 779       //serarch for charged trigger at the MC generator level
 780       Int_t    indexTriggGen = -1;
 781       Double_t ptTriggGen    = -1;
 782       Int_t    iCounterGen   =  0;
 783       Int_t    triggersMC[200];//list of trigger candidates
 784       Int_t    ntriggersMC   = 0; //index in triggers array
 785
 786       if(fESD){//ESD input
 787
 788          AliMCEvent* mcEvent = MCEvent();
 789          if(!mcEvent){
 790             PostData(1, fOutputList);
 791             return;
 792          }
 793
 794          AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
 795          if(pythiaGenHeader){
 796             nTrials = pythiaGenHeader->Trials();
 797             ptHard  = pythiaGenHeader->GetPtHard();
 798             fh1PtHard->Fill(ptHard,eventW);
 799             fh1PtHardNoW->Fill(ptHard,1);
 800             fh1PtHardTrials->Fill(ptHard,nTrials);
 801          }
 802
 803          for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){
 804             if(!mcEvent->IsPhysicalPrimary(it)) continue;
 805             AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it);
 806             if(!part) continue;
 807             if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
 808
 809                if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
 810                   if(indexTriggGen > -1){//trigger candidater was found
 811                      triggersMC[ntriggersMC] = indexTriggGen;
 812                      ntriggersMC++;
 813                   }
 814                }
 815
 816                iCounterGen++;
 817             }
 818          }
 819       }else{  //AOD input
 820          if(!fAODIn){
 821             PostData(1, fOutputList);
 822             return;
 823          }
 824          TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName()));
 825          if(!tca){
 826             PostData(1, fOutputList);
 827             return;
 828          }
 829          for(Int_t it = 0; it < tca->GetEntriesFast(); it++){
 830             AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
 831             if(!part) continue;
 832             if(!part->IsPhysicalPrimary()) continue;
 833             if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){
 834
 835                if(fHardest==0 && ntriggersMC<200){//single inclusive trigger
 836                   if(indexTriggGen > -1){ //trigger candidater was found
 837                      triggersMC[ntriggersMC] = indexTriggGen;
 838                      ntriggersMC++;
 839                   }
 840                }
 841
 842                iCounterGen++;
 843             }
 844          }
 845       }
 846
 847       if(fHardest==0){ //single inclusive trigger
 848          if(ntriggersMC>0){ //there is at least one trigger
 849             Int_t rnd     = fRandom->Integer(ntriggersMC); //0 to ntriggers-1
 850             indexTriggGen = triggersMC[rnd];
 851          }else{
 852             indexTriggGen = -1; //trigger not found
 853          }
 854       }
 855
 856       //================== Fill jet list ===================
 857       if(aodGenJets){
 858          if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMC.Data(), aodGenJets->GetEntriesFast());
 859
 860          for(Int_t igJet = 0; igJet < aodGenJets->GetEntriesFast(); igJet++) {
 861             AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJets)[igJet]);
 862             if(jetGen) fListJetsGen->Add(jetGen);
 863          }
 864       }
 865
 866       //============  Generator trigger+jet ==================
 867       Int_t ilowGen  = (fHardest==0 || fHardest==1) ? indexTriggGen : 0;
 868       Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries();
 869
 870       for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger
 871          indexTriggGen = igen; //trigger hadron
 872
 873          if(indexTriggGen == -1) continue;
 874          AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen);
 875          if(!triggerGen) continue;
 876
 877          if(fHardest >= 2){
 878             if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>10
 879          }
 880          if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue;
 881
 882          ptTriggGen = triggerGen->Pt(); //count triggers
 883          fh2NtriggersGen->Fill(centValue, ptTriggGen);
 884
 885          //Count jets and trigger-jet pairs at MC  generator level
 886          if(!aodGenJets) continue;
 887          for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
 888             AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
 889             if(!jet) continue;
 890             Double_t etaJetGen = jet->Eta();
 891
 892             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
 893
 894                Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi());
 895                if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
 896
 897                //Centrality, A, pT, pTtrigg
 898                Double_t fillspecgen[] = { centValue,
 899                                           jet->EffectiveAreaCharged(),
 900                                           jet->Pt(),
 901                                           ptTriggGen
 902                                         };
 903               fHJetSpecGen->Fill(fillspecgen);
 904             }//back to back jet-trigger pair
 905          }//jet loop
 906       }//trigger loop
 907
 908
 909       //================ RESPONSE MATRIX ===============
 910       if(aodGenJets){
 911          //Count jets and trigger-jet pairs at MC  generator level
 912          for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){
 913             AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij));
 914             if(!jet) continue;
 915             Double_t etaJetGen = jet->Eta();
 916             Double_t ptJetGen  = jet->Pt();
 917
 918             if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){
 919                fhJetPtGen->Fill(ptJetGen); // generator level pt spectum of jets response mx normalization
 920             }
 921          }
 922          if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets
 923             Int_t ng = (Int_t) fListJetsGen->GetEntries();
 924             Int_t nr = (Int_t) fListJets->GetEntries();
 925
 926             //Find closest MC generator - reconstructed jets
 927             if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet
 928             if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet
 929
 930             if(fDebug){
 931                Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize());
 932                Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize());
 933             }
 934             //matching of MC genrator level and reconstructed jets
 935             AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug);
 936
 937             // Fill response matrix
 938             for(Int_t ir = 0; ir < nr; ir++){
 939                AliAODJet *recJet  = (AliAODJet*) fListJets->At(ir);
 940                Double_t etaJetRec = recJet->Eta();
 941                Double_t ptJetRec  = recJet->Pt();
 942                //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
 943
 944                if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){
 945                   Int_t ig = faGenIndex[ir]; //associated generator level jet
 946                   if(ig >= 0 && ig < ng){
 947                      if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
 948                      AliAODJet *genJet  = (AliAODJet*) fListJetsGen->At(ig);
 949                      Double_t ptJetGen  = genJet->Pt();
 950                      Double_t etaJetGen = genJet->Eta();
 951
 952                   //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc
 953                      if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){
 954                         fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen);
 955                      }
 956                   }//ig>=0
 957                }//rec jet in eta acceptance
 958             }//loop over reconstructed jets
 959          }// # of  rec jets >0
 960       }//pointer MC generator jets
 961    } //analyze generator level MC
 962
 963    //=============  RECONSTRUCTED INCLUSIVE JETS ===============
 964
 965    Double_t etaJet  = 0.0;
 966    Double_t pTJet   = 0.0;
 967    Double_t areaJet = 0.0;
 968    Double_t phiJet  = 0.0;
 969    Int_t indexLeadingJet     = -1;
 970    Double_t pTLeadingJet     = -10.0;
 971    Double_t areaLeadingJet   = -10.0;
 972
 973    for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
 974       AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
 975       if(!jet) continue;
 976       etaJet  = jet->Eta();
 977       phiJet  = jet->Phi();
 978       pTJet   = jet->Pt();
 979       if(pTJet==0) continue;
 980
 981       if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
 982       areaJet = jet->EffectiveAreaCharged();
 983
 984       //Jet Diagnostics---------------------------------
 985       fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi
 986       fhJetEta->Fill(pTJet, etaJet);
 987       //search for leading jet
 988       if(pTJet > pTLeadingJet){
 989          indexLeadingJet  = ij;
 990          pTLeadingJet     = pTJet;
 991          areaLeadingJet   = areaJet;
 992       }
 993
 994       // raw spectra of INCLUSIVE jets
 995       //Centrality, pTjet, A
 996       Double_t fillraw[] = { centValue,
 997                              pTJet,
 998                              areaJet
 999                            };
1000       fHJetPtRaw->Fill(fillraw);
1001    }
1002
1003    if(indexLeadingJet > -1){
1004       // raw spectra of LEADING jets
1005       //Centrality, pTjet,  A
1006       Double_t fillleading[] = { centValue,
1007                                  pTLeadingJet,
1008                                  areaLeadingJet
1009                                };
1010       fHLeadingJetPtRaw->Fill(fillleading);
1011    }
1012
1013
1014    // ===============  RECONSTRUCTED TRIGGER-JET PAIRS ================
1015    //Find Hadron trigger
1016    TList particleList; //list of tracks
1017    Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
1018    if(fIsMC) FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos
1019
1020    //set ranges of the trigger loop
1021    Int_t ilow  = (fHardest==0 || fHardest==1) ? indexTrigg : 0;
1022    Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries();
1023
1024    for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger
1025       indexTrigg = itrk; //trigger hadron with pT >10 GeV
1026
1027       if(indexTrigg < 0) continue;
1028
1029       AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
1030       if(!triggerHadron){
1031          PostData(1, fOutputList);
1032          return;
1033       }
1034       if(fHardest >= 2){
1035          if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6
1036       }
1037       if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue;
1038
1039       //Fill trigger histograms
1040       fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
1041       fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
1042       fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta());
1043
1044
1045       //---------- make trigger-jet pairs ---------
1046       Int_t injet4     = 0;
1047       Int_t injet      = 0;
1048
1049       for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
1050          AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
1051          if(!jet) continue;
1052          etaJet  = jet->Eta();
1053          phiJet  = jet->Phi();
1054          pTJet   = jet->Pt();
1055          if(pTJet==0) continue;
1056
1057          if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
1058          areaJet = jet->EffectiveAreaCharged();
1059          if(areaJet >= 0.07) injet++;
1060          if(areaJet >= 0.4)  injet4++;
1061
1062          Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
1063          fhDphiTriggerJet->Fill(dphi); //Input
1064
1065          //Dphi versus jet pT
1066          //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg
1067          Double_t filldp[] = { centValue,
1068                                dphi,
1069                                pTJet,
1070                                triggerHadron->Pt()
1071                              };
1072          fHDphiVsJetPtAll->Fill(filldp);
1073
1074          // Select back to back trigger - jet pairs
1075          if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue;
1076          fhDphiTriggerJetAccept->Fill(dphi); //Accepted
1077
1078
1079          //Centrality, A, pTjet, pTtrigg
1080          Double_t fillspec[] = { centValue,
1081                                  areaJet,
1082                                  pTJet,
1083                                  triggerHadron->Pt()
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