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