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 "TMath.h"
  30 #include "TH1F.h"
  31 #include "TH1D.h"
  32 #include "TH2F.h"
  33 #include "TH3F.h"
  34 #include "THnSparse.h"
  35 #include "TCanvas.h"
  36
  37 #include "AliLog.h"
  38
  39 #include "AliAnalysisTask.h"
  40 #include "AliAnalysisManager.h"
  41
  42 #include "AliVEvent.h"
  43 #include "AliESDEvent.h"
  44 #include "AliESDInputHandler.h"
  45 #include "AliCentrality.h"
  46 #include "AliAnalysisHelperJetTasks.h"
  47 #include "AliInputEventHandler.h"
  48 #include "AliAODJetEventBackground.h"
  49 #include "AliAODMCParticle.h"
  50 //#include "AliAnalysisTaskFastEmbedding.h"
  51 #include "AliAODEvent.h"
  52 #include "AliAODHandler.h"
  53 #include "AliAODJet.h"
  54
  55 #include "AliAnalysisTaskJetCorePP.h"
  56
  57 using std::cout;
  58 using std::endl;
  59
  60 ClassImp(AliAnalysisTaskJetCorePP)
  61
  62 //Filip Krizek 1st March 2013
  63
  64 //---------------------------------------------------------------------
  65 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() :
  66 AliAnalysisTaskSE(),
  67 fESD(0x0),
  68 fAODIn(0x0),
  69 fAODOut(0x0),
  70 fAODExtension(0x0),
  71 fJetBranchName(""),
  72 //fListJets(0x0),
  73 fNonStdFile(""),
  74 fSystem(0), //pp=0  pPb=1
  75 fJetParamR(0.4),
  76 fOfflineTrgMask(AliVEvent::kAny),
  77 fMinContribVtx(1),
  78 fVtxZMin(-10.0),
  79 fVtxZMax(10.0),
  80 fFilterMask(0),
  81 fCentMin(0.0),
  82 fCentMax(100.0),
  83 fJetEtaMin(-0.5),
  84 fJetEtaMax(0.5),
  85 fTriggerEtaCut(0.9),
  86 fTrackEtaCut(0.9),
  87 fTrackLowPtCut(0.15),
  88 fOutputList(0x0),
  89 fHistEvtSelection(0x0),
  90 fh2Ntriggers(0x0),
  91 fHJetSpec(0x0),
  92 fHJetDensity(0x0),
  93 fHJetDensityA4(0x0),
  94 fhJetPhi(0x0),
  95 fhTriggerPhi(0x0),
  96 fhJetEta(0x0),
  97 fhTriggerEta(0x0),
  98 fhVertexZ(0x0),
  99 fhVertexZAccept(0x0),
 100 fhContribVtx(0x0),
 101 fhContribVtxAccept(0x0),
 102 fhDphiTriggerJet(0x0),
 103 fhDphiTriggerJetAccept(0x0),
 104 fhCentrality(0x0),
 105 fhCentralityAccept(0x0),
 106 fkAcceptance(2.0*TMath::Pi()*1.8)
 107 {
 108    // default Constructor
 109    fListJets = new TList();
 110 }
 111
 112 //---------------------------------------------------------------------
 113
 114 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) :
 115 AliAnalysisTaskSE(name),
 116 fESD(0x0),
 117 fAODIn(0x0),
 118 fAODOut(0x0),
 119 fAODExtension(0x0),
 120 fJetBranchName(""),
 121 //fListJets(0x0),
 122 fNonStdFile(""),
 123 fSystem(0),  //pp=0   pPb=1
 124 fJetParamR(0.4),
 125 fOfflineTrgMask(AliVEvent::kAny),
 126 fMinContribVtx(1),
 127 fVtxZMin(-10.0),
 128 fVtxZMax(10.0),
 129 fFilterMask(0),
 130 fCentMin(0.0),
 131 fCentMax(100.0),
 132 fJetEtaMin(-0.5),
 133 fJetEtaMax(0.5),
 134 fTriggerEtaCut(0.9),
 135 fTrackEtaCut(0.9),
 136 fTrackLowPtCut(0.15),
 137 fOutputList(0x0),
 138 fHistEvtSelection(0x0),
 139 fh2Ntriggers(0x0),
 140 fHJetSpec(0x0),
 141 fHJetDensity(0x0),
 142 fHJetDensityA4(0x0),
 143 fhJetPhi(0x0),
 144 fhTriggerPhi(0x0),
 145 fhJetEta(0x0),
 146 fhTriggerEta(0x0),
 147 fhVertexZ(0x0),
 148 fhVertexZAccept(0x0),
 149 fhContribVtx(0x0),
 150 fhContribVtxAccept(0x0),
 151 fhDphiTriggerJet(0x0),
 152 fhDphiTriggerJetAccept(0x0),
 153 fhCentrality(0x0),
 154 fhCentralityAccept(0x0),
 155 fkAcceptance(2.0*TMath::Pi()*1.8)
 156 {
 157    // Constructor
 158    fListJets = new TList();
 159
 160    DefineOutput(1, TList::Class());
 161 }
 162
 163 //--------------------------------------------------------------
 164 AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a):
 165 AliAnalysisTaskSE(a.GetName()),
 166 fESD(a.fESD),
 167 fAODIn(a.fAODIn),
 168 fAODOut(a.fAODOut),
 169 fAODExtension(a.fAODExtension),
 170 fJetBranchName(a.fJetBranchName),
 171 fListJets(a.fListJets),
 172 fNonStdFile(a.fNonStdFile),
 173 fSystem(a.fSystem),
 174 fJetParamR(a.fJetParamR),
 175 fOfflineTrgMask(a.fOfflineTrgMask),
 176 fMinContribVtx(a.fMinContribVtx),
 177 fVtxZMin(a.fVtxZMin),
 178 fVtxZMax(a.fVtxZMax),
 179 fFilterMask(a.fFilterMask),
 180 fCentMin(a.fCentMin),
 181 fCentMax(a.fCentMax),
 182 fJetEtaMin(a.fJetEtaMin),
 183 fJetEtaMax(a.fJetEtaMax),
 184 fTriggerEtaCut(a.fTriggerEtaCut),
 185 fTrackEtaCut(a.fTrackEtaCut),
 186 fTrackLowPtCut(a.fTrackLowPtCut),
 187 fOutputList(a.fOutputList),
 188 fHistEvtSelection(a.fHistEvtSelection),
 189 fh2Ntriggers(a.fh2Ntriggers),
 190 fHJetSpec(a.fHJetSpec),
 191 fHJetDensity(a.fHJetDensity),
 192 fHJetDensityA4(a.fHJetDensityA4),
 193 fhJetPhi(a.fhJetPhi),
 194 fhTriggerPhi(a.fhTriggerPhi),
 195 fhJetEta(a.fhJetEta),
 196 fhTriggerEta(a.fhTriggerEta),
 197 fhVertexZ(a.fhVertexZ),
 198 fhVertexZAccept(a.fhVertexZAccept),
 199 fhContribVtx(a.fhContribVtx),
 200 fhContribVtxAccept(a.fhContribVtxAccept),
 201 fhDphiTriggerJet(a.fhDphiTriggerJet),
 202 fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept),
 203 fhCentrality(a.fhCentrality),
 204 fhCentralityAccept(a.fhCentralityAccept),
 205 fkAcceptance(a.fkAcceptance)
 206 {
 207    //Copy Constructor
 208 }
 209 //--------------------------------------------------------------
 210
 211 AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){
 212   // assignment operator
 213   this->~AliAnalysisTaskJetCorePP();
 214   new(this) AliAnalysisTaskJetCorePP(a);
 215   return *this;
 216 }
 217 //--------------------------------------------------------------
 218
 219 AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP()
 220 {
 221    //Destructor
 222    delete fListJets;
 223    delete fOutputList; // ????
 224 }
 225
 226 //--------------------------------------------------------------
 227
 228 void AliAnalysisTaskJetCorePP::Init()
 229 {
 230    // check for jet branches
 231    if(!strlen(fJetBranchName.Data())){
 232       AliError("Jet branch name not set.");
 233    }
 234
 235 }
 236
 237 //--------------------------------------------------------------
 238
 239 void AliAnalysisTaskJetCorePP::UserCreateOutputObjects()
 240 {
 241
 242
 243   // Create histograms
 244    // Called once
 245    OpenFile(1);
 246    if(!fOutputList) fOutputList = new TList();
 247    fOutputList->SetOwner(kTRUE);
 248
 249    Bool_t oldStatus = TH1::AddDirectoryStatus();
 250    TH1::AddDirectory(kFALSE);
 251
 252    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 253    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 254    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 255    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 256    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 257    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 258    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
 259
 260    fOutputList->Add(fHistEvtSelection);
 261
 262    Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10
 263
 264    fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
 265                              nBinsCentrality,0.0,100.0,50,0.0,50.0);
 266
 267    //Centrality, A, pT - rho*A, pTtrigg, pT, rho*A
 268    const Int_t dimSpec   = 6;
 269    const Int_t nBinsSpec[dimSpec]     = {nBinsCentrality, 100,  140,  50, 100,  50};
 270    const Double_t lowBinSpec[dimSpec] = {0.0,             0.0,-80.0, 0.0, 0.0, 0.0};
 271    const Double_t hiBinSpec[dimSpec]  = {100.0,           1.0,200.0,50.0, 200, 20.0};
 272    fHJetSpec = new THnSparseF("fHJetSpec","Recoil jet spectrum [cent,A,pT-rho*A,pTtrig]",
 273                                           dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
 274    fOutputList->Add(fHJetSpec);
 275
 276    //------------------- HISTOS FOR DIAGNOSTIC ----------------------
 277    //Jet number density histos [Trk Mult, jet density, pT trigger]
 278    const Int_t    dimJetDens   = 3;
 279    const Int_t    nBinsJetDens[dimJetDens]  = {100,   100, 10};
 280    const Double_t lowBinJetDens[dimJetDens] = {0.0,   0.0,  0.0};
 281    const Double_t hiBinJetDens[dimJetDens]  = {500.0, 5.0, 50.0 };
 282
 283    fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07",
 284                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
 285
 286    fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4",
 287                                    dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens);
 288
 289    fOutputList->Add(fh2Ntriggers);
 290    fOutputList->Add(fHJetDensity);
 291    fOutputList->Add(fHJetDensityA4);
 292
 293
 294    //inclusive azimuthal and pseudorapidity histograms
 295    fhJetPhi = new TH1D("fhJetPhi","Azim dist jets",50,-TMath::Pi(),TMath::Pi());
 296    fhTriggerPhi= new TH1D("fhTriggerPhi","azim dist trig had",50,-TMath::Pi(),TMath::Pi());
 297    fhJetEta = new TH1D("fhJetEta","Eta dist jets",40,-0.9,0.9);
 298    fhTriggerEta = new TH1D("fhTriggerEta","Eta dist trig had",40,-0.9,0.9);
 299    fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20);
 300    fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20);
 301    fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200);
 302    fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
 303    fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi());
 304    fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi());
 305    fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100);
 306    fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100);
 307
 308    fOutputList->Add(fhJetPhi);
 309    fOutputList->Add(fhTriggerPhi);
 310    fOutputList->Add(fhJetEta);
 311    fOutputList->Add(fhTriggerEta);
 312    fOutputList->Add(fhVertexZ);
 313    fOutputList->Add(fhVertexZAccept);
 314    fOutputList->Add(fhContribVtx);
 315    fOutputList->Add(fhContribVtxAccept);
 316    fOutputList->Add(fhDphiTriggerJet);
 317    fOutputList->Add(fhDphiTriggerJetAccept);
 318    fOutputList->Add(fhCentrality);
 319    fOutputList->Add(fhCentralityAccept);
 320
 321
 322    // =========== Switch on Sumw2 for all histos ===========
 323    for(Int_t i=0; i<fOutputList->GetEntries(); i++){
 324       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 325       if(h1){
 326          h1->Sumw2();
 327          continue;
 328       }
 329       THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 330       if(hn){
 331          hn->Sumw2();
 332       }
 333    }
 334    TH1::AddDirectory(oldStatus);
 335
 336    PostData(1, fOutputList);
 337 }
 338
 339 //--------------------------------------------------------------------
 340
 341 void AliAnalysisTaskJetCorePP::UserExec(Option_t *)
 342 {
 343
 344    //Event loop
 345
 346    if(TMath::Abs((Float_t) fJetParamR)<0.00001){
 347       AliError("Cone radius is set to zero.");
 348       return;
 349    }
 350    if(!strlen(fJetBranchName.Data())){
 351       AliError("Jet branch name not set.");
 352       return;
 353    }
 354
 355    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 356    if(!fESD){
 357       AliError("ESD not available");
 358       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 359    }
 360
 361    fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 362    AliAODEvent* aod = NULL;
 363    // take all other information from the aod we take the tracks from
 364    if(!aod){
 365       if(!fESD) aod = fAODIn;
 366       else      aod = fAODOut;
 367    }
 368
 369    if(fNonStdFile.Length()!=0){
 370       // case that we have an AOD extension we can fetch the jets from the extended output
 371       AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 372       fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0;
 373       if(!fAODExtension){
 374          if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data());
 375       }
 376    }
 377
 378    // ----------------- event selection --------------------------
 379    fHistEvtSelection->Fill(1); // number of events before event selection
 380
 381    // physics selection
 382    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 383            ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 384
 385    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 386       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 387       fHistEvtSelection->Fill(2);
 388       PostData(1, fOutputList);
 389       return;
 390    }
 391
 392    //check AOD pointer
 393    if(!aod){
 394       if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
 395       fHistEvtSelection->Fill(3);
 396       PostData(1, fOutputList);
 397       return;
 398    }
 399
 400    // vertex selection
 401    AliAODVertex* primVtx = aod->GetPrimaryVertex();
 402
 403    if(!primVtx){
 404       if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
 405       fHistEvtSelection->Fill(3);
 406       PostData(1, fOutputList);
 407       return;
 408    }
 409
 410    Int_t nTracksPrim = primVtx->GetNContributors();
 411    Float_t vtxz = primVtx->GetZ();
 412    //Input events
 413    fhContribVtx->Fill(nTracksPrim);
 414    if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz);
 415
 416    if((nTracksPrim < fMinContribVtx) ||
 417       (vtxz < fVtxZMin) ||
 418       (vtxz > fVtxZMax)){
 419       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",
 420                          (char*)__FILE__,__LINE__,vtxz);
 421       fHistEvtSelection->Fill(3);
 422       PostData(1, fOutputList);
 423       return;
 424    }else{
 425       //Accepted events
 426       fhContribVtxAccept->Fill(nTracksPrim);
 427       fhVertexZAccept->Fill(vtxz);
 428    }
 429    //FK// No event class selection imposed
 430    // event class selection (from jet helper task)
 431    //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 432    //if(fDebug) Printf("Event class %d", eventClass);
 433    //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 434    //   fHistEvtSelection->Fill(4);
 435    //   PostData(1, fOutputList);
 436    //   return;
 437    //}
 438
 439    // centrality selection
 440    AliCentrality *cent = 0x0;
 441    Double_t centValue  = 0.0;
 442    if(fSystem){  //fSystem=0 for pp,   fSystem=1 for pPb
 443       if(fESD){
 444          cent = fESD->GetCentrality();
 445          if(cent) centValue = cent->GetCentralityPercentile("V0M");
 446       }else{
 447          centValue = aod->GetHeader()->GetCentrality();
 448       }
 449       if(fDebug) printf("centrality: %f\n", centValue);
 450       //Input events
 451       fhCentrality->Fill(centValue);
 452
 453       if(centValue < fCentMin || centValue > fCentMax){
 454          fHistEvtSelection->Fill(4);
 455          PostData(1, fOutputList);
 456          return;
 457       }else{
 458          //Accepted events
 459          fhCentralityAccept->Fill( centValue );
 460       }
 461    }
 462
 463    if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl;
 464
 465    fHistEvtSelection->Fill(0); // accepted events
 466
 467    // ------------------- end event selection --------------------
 468
 469    // fetch jets
 470    TClonesArray *aodJets = 0x0;
 471
 472    if(fAODOut && !aodJets){
 473       aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data()));
 474    }
 475    if(fAODExtension && !aodJets){
 476       aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data()));
 477    }
 478    if(fAODIn && !aodJets){
 479       aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data()));
 480    }
 481
 482    // ------------- Hadron trigger --------------
 483    TList particleList; //list of tracks
 484    Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list
 485
 486    if(indexTrigg<0) return; // no trigger track found
 487
 488    AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg);
 489    if(!triggerHadron){
 490       PostData(1, fOutputList);
 491       return;
 492    }
 493
 494
 495    fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT
 496
 497  //  if(triggerHadron->Pt() > 10.0){}
 498    if(triggerHadron->Pt() > 3.0){
 499       //Trigger Diagnostics---------------------------------
 500       fhTriggerPhi->Fill(RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi
 501       fhTriggerEta->Fill(triggerHadron->Eta());
 502    }
 503
 504    //--------- Fill list of jets --------------
 505    fListJets->Clear();
 506    if(aodJets){
 507       if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(),
 508                                       aodJets->GetEntriesFast());
 509       for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) {
 510          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]);
 511          if (jet) fListJets->Add(jet);
 512       }
 513    }
 514
 515    Double_t etaJet  = 0.0;
 516    Double_t pTJet   = 0.0;
 517    Double_t areaJet = 0.0;
 518    Double_t phiJet  = 0.0;
 519    Int_t injet4 = 0;
 520    Int_t injet  = 0;
 521    //---------- jet loop ---------
 522    for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){
 523       AliAODJet* jet = (AliAODJet*)(fListJets->At(ij));
 524       if(!jet) continue;
 525       etaJet  = jet->Eta();
 526       phiJet  = jet->Phi();
 527       pTJet   = jet->Pt();
 528       if(pTJet==0) continue;
 529
 530       if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue;
 531       areaJet = jet->EffectiveAreaCharged();
 532
 533       //Jet Diagnostics---------------------------------
 534       fhJetPhi->Fill(RelativePhi(phiJet,0.0)); //phi -pi,pi
 535       fhJetEta->Fill(etaJet);
 536       if(areaJet >= 0.07) injet++;
 537       if(areaJet >= 0.4)  injet4++;
 538       //--------------------------------------------------
 539
 540       Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet);
 541
 542       fhDphiTriggerJet->Fill(dphi); //Input
 543       if(TMath::Abs((Double_t) dphi) < TMath::Pi()-0.6) continue;
 544       fhDphiTriggerJetAccept->Fill(dphi); //Accepted
 545
 546       //Background w.r.t. jet axis
 547       Double_t pTBckWrtJet =
 548          GetBackgroundInPerpCone(fJetParamR, phiJet, etaJet, &particleList);
 549
 550       Double_t ratioOfAreas = areaJet/(TMath::Pi()*fJetParamR*fJetParamR);
 551       Double_t rhoA = ratioOfAreas*pTBckWrtJet; //bg activity in a cone of similar size
 552       Double_t ptcorr = pTJet - rhoA; //Pt Jet UE subtr
 553
 554
 555       //Centrality, A, pT - rho*A, pTtrigg, pT, rho*A
 556       Double_t fillspec[] = { centValue,
 557                               areaJet,
 558                               ptcorr,
 559                               triggerHadron->Pt(),
 560                               pTJet,
 561                               rhoA
 562                             };
 563       fHJetSpec->Fill(fillspec);
 564
 565    }//jet loop
 566
 567    //Fill Jet Density In the Event A>0.07
 568    if(injet>0){
 569       Double_t filldens[]={ (Double_t) particleList.GetEntries(),
 570                             injet/fkAcceptance,
 571                             triggerHadron->Pt()
 572                           };
 573       fHJetDensity->Fill(filldens);
 574    }
 575
 576    //Fill Jet Density In the Event A>0.4
 577    if(injet4>0){
 578       Double_t filldens4[]={ (Double_t) particleList.GetEntries(),
 579                              injet4/fkAcceptance,
 580                              triggerHadron->Pt()
 581                            };
 582       fHJetDensityA4->Fill(filldens4);
 583    }
 584
 585    PostData(1, fOutputList);
 586 }
 587
 588 //----------------------------------------------------------------------------
 589 void AliAnalysisTaskJetCorePP::Terminate(const Option_t *)
 590 {
 591    // Draw result to the screen
 592    // Called once at the end of the query
 593
 594    if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n");
 595    if(!GetOutputData(1)) return;
 596 }
 597
 598
 599 //----------------------------------------------------------------------------
 600 Int_t  AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){
 601    //Fill the list of accepted tracks (passed track cut)
 602    //return consecutive index of the hardest ch hadron in the list
 603    Int_t iCount        = 0;
 604    AliAODEvent *aodevt = NULL;
 605
 606    if(!fESD) aodevt = fAODIn;
 607    else      aodevt = fAODOut;
 608
 609    if(!aodevt) return -1;
 610
 611    Int_t    index = -1; //index of the highest particle in the list
 612    Double_t ptmax = -10;
 613
 614    for(int it = 0; it < aodevt->GetNumberOfTracks(); it++){
 615       AliAODTrack *tr = aodevt->GetTrack(it);
 616
 617       //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue;
 618       if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue;
 619       if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue;
 620       if(tr->Pt() < fTrackLowPtCut) continue;
 621       list->Add(tr);
 622       if(tr->Pt()>ptmax){
 623          ptmax = tr->Pt();
 624          index = iCount;
 625       }
 626       iCount++;
 627    }
 628
 629    if(index>-1){ //check pseudorapidity cut on trigger
 630       AliAODTrack *trigger = (AliAODTrack*) list->At(index);
 631       if(trigger && TMath::Abs((Float_t) trigger->Eta())< fTriggerEtaCut){ return index;}
 632       return -1;
 633    }else{
 634      return -1;
 635    }
 636 }
 637
 638 //----------------------------------------------------------------------------
 639
 640 Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){
 641    //calculate sum of track pT in the cone perpendicular in phi to the jet
 642    //jetR = cone radius
 643    // jetPhi, jetEta = direction of the jet
 644    Int_t numberOfTrks = trkList->GetEntries();
 645    Double_t pTsum = 0.0;
 646    Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi
 647    for(Int_t it=0; it<numberOfTrks; it++){
 648       AliVParticle *trk = (AliVParticle*) trkList->At(it);
 649       Double_t dphi = RelativePhi(perpConePhi,trk->Phi());
 650       Double_t deta = trk->Eta()-jetEta;
 651
 652       if( (dphi*dphi + deta*deta)< (jetR*jetR)){
 653          pTsum += trk->Pt();
 654       }
 655    }
 656
 657    return pTsum;
 658 }
 659
 660 //----------------------------------------------------------------------------
 661
 662 Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){
 663    //Get relative azimuthal angle of two particles -pi to pi
 664    if      (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
 665    else if (vphi > TMath::Pi())  vphi -= TMath::TwoPi();
 666
 667    if      (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
 668    else if (mphi > TMath::Pi())  mphi -= TMath::TwoPi();
 669
 670    Double_t dphi = mphi - vphi;
 671    if      (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
 672    else if (dphi > TMath::Pi())  dphi -= TMath::TwoPi();
 673
 674    return dphi;//dphi in [-Pi, Pi]
 675 }
 676
 677