PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskHJetSpectra.cxx

   1 #ifndef ALIANALYSISTASKSE_H
   2
   3 #include <Riostream.h>
   4 #include <TROOT.h>
   5 #include <TFile.h>
   6 #include <TChain.h>
   7 #include <TTree.h>
   8 #include <TKey.h>
   9 #include <TProfile.h>
  10 #include <TProfile2D.h>
  11 #include <TH1.h>
  12 #include <TH1F.h>
  13 #include <TH2F.h>
  14 #include <TH1D.h>
  15 #include <TH1I.h>
  16 #include <TArrayF.h>
  17 #include <THnSparse.h>
  18 #include <TCanvas.h>
  19 #include <TList.h>
  20 #include <TClonesArray.h>
  21 #include <TObject.h>
  22 #include <TMath.h>
  23 #include <TSystem.h>
  24 #include <TInterpreter.h>
  25 #include "AliAnalysisTask.h"
  26 #include "AliCentrality.h"
  27 #include "AliStack.h"
  28 #include "AliESDEvent.h"
  29 #include "AliESDInputHandler.h"
  30 #include "AliAODEvent.h"
  31 #include "AliAODHandler.h"
  32 #include "AliAnalysisManager.h"
  33 #include "AliAnalysisTaskSE.h"
  34 #include "AliAnalysisHelperJetTasks.h"
  35 #include "AliInputEventHandler.h"
  36 #endif
  37
  38 #include <time.h>
  39 #include <TRandom3.h>
  40 #include "AliGenEventHeader.h"
  41 #include "AliGenPythiaEventHeader.h"
  42 #include "AliGenHijingEventHeader.h"
  43 #include "AliAODMCHeader.h"
  44 #include "AliMCEvent.h"
  45 #include "AliLog.h"
  46 #include <AliEmcalJet.h>
  47 #include <AliPicoTrack.h>
  48 #include "AliVEventHandler.h"
  49 #include "AliVParticle.h"
  50 #include "AliAODMCParticle.h"
  51 #include "AliAnalysisUtils.h"
  52 #include "AliRhoParameter.h"
  53 #include "TVector3.h"
  54 #include "AliVVertex.h"
  55
  56 #include <stdio.h>
  57 #include <stdlib.h>
  58 #include <vector>
  59
  60 #include "AliAnalysisTaskHJetSpectra.h"
  61 using std::min;
  62
  63 // ANALYSIS OF HIGH PT HADRON TRIGGER ASSOCIATED SPECTRUM OF RECOIL JETS IN P+PB
  64 // Author Filip Krizek   (17.May. 2014)
  65
  66 //TODO: Not accessing the particles when using MC
  67 //TODO: FillHistogram can be done better with virtual TH1(?)
  68 ClassImp(AliAnalysisTaskHJetSpectra)
  69 //________________________________________________________________________________________
  70
  71 AliAnalysisTaskHJetSpectra::AliAnalysisTaskHJetSpectra():
  72 AliAnalysisTaskSE(), fOutputList(0), fAnalyzePythia(0), fAnalyzeHijing(0),  fIsKinematics(0), fUseDefaultVertexCut(1), fUsePileUpCut(1),
  73 fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0),  fRhoTaskName(),
  74 fRandConeRadius(0.4),fRandConeRadiusSquared(fRandConeRadius*fRandConeRadius), fSignalJetRadius(0.4), fBackgroundJetRadius(0.3), fBackgroundJetPtMin(15.0),
  75 fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.5), fNumberOfCentralityBins(20), fCentralityType("V0A"),
  76 fCrossSection(0.0), fTrials(0.0), fImpParam(-1.0), fRandom(0), fHelperClass(0), fInitialized(0),
  77 fTTlow(8.0), fTThigh(9.0), fTTtype(0), fDphiCut(TMath::Pi()-0.6), fUseDoubleBinPrecision(0),
  78 fHistEvtSelection(0x0), fh2Ntriggers(0x0), fHJetSpec(0x0), fHJetSpecSubUeMedian(0x0), fHJetSpecSubUeCone(0x0), fHJetSpecSubUeCMS(0x0),
  79 fhRhoCellMedian(0x0), fhRhoCone(0x0), fhRhoCMS(0x0),
  80 fhRhoCellMedianIncl(0x0), fhRhoConeIncl(0x0), fhRhoCMSIncl(0x0),
  81 fARhoCellMedian(0x0), fARhoCone(0x0), fARhoCMS(0x0),
  82 fhDeltaPtMedian(0x0), fhDeltaPtCone(0x0), fhDeltaPtCMS(0x0),
  83 fhDeltaPtMedianIncl(0x0), fhDeltaPtConeIncl(0x0), fhDeltaPtCMSIncl(0x0),
  84 fhDeltaPtMedianNearSide(0x0), fhDeltaPtMedianAwaySide(0x0), fhDeltaPtCMSNearSide(0x0), fhDeltaPtCMSAwaySide(0x0),
  85 fhDeltaPtMedianExclTrigCone(0x0),fhDeltaPtCMSExclTrigCone(0x0), fhDeltaPtMedianExclAwayJet(0x0), fhDeltaPtCMSExclAwayJet(0x0),
  86 fhJetPhi(0x0), fhTrackPhi(0x0), fhJetEta(0x0), fhTrackEta(0x0), fhTrackCentVsPt(0x0), fhVertexZ(0x0), fhVertexZAccept(0x0),
  87 fhDphiTriggerJetMinBias(0x0),fhDphiTriggerJetCent20(0x0), fhDphiTriggerJetAccept(0x0),
  88 fhCentrality(0x0), fhCentralityV0M(0x0), fhCentralityV0A(0x0), fhCentralityV0C(0x0), fhCentralityZNA(0x0),
  89 fNofRndTrials(2000), fJetFreeAreaFrac(0.8), fnEta(2), fnPhi(11), fEtaSize(0.9), fPhiSize(2*TMath::Pi()/fnPhi), fCellArea(fPhiSize*fEtaSize),
  90 fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0), fhImpactParameter(0x0), fhImpactParameterTT(0x0),
  91 fNofRandomCones(1),
  92 fRConesR(0.1),fRConesRSquared(fRConesR*fRConesR),fnRCones(16)
  93 {
  94     //default constructor
  95     for(Int_t k=0; k<50; k++){
  96        fRConePhi[k] = 0.0;
  97        fRConeEta[k] = 0.0;
  98     }
  99 }
 100
 101 //________________________________________________________________________
 102 AliAnalysisTaskHJetSpectra::AliAnalysisTaskHJetSpectra(const char *name, const char* trackArrayName, const char* jetArrayName, const char* backgroundJetArrayName) :
 103 AliAnalysisTaskSE(name), fOutputList(0), fAnalyzePythia(0), fAnalyzeHijing(0), fIsKinematics(0), fUseDefaultVertexCut(1), fUsePileUpCut(1),
 104  fJetArray(0), fTrackArray(0), fBackgroundJetArray(0), fJetArrayName(0), fTrackArrayName(0), fBackgroundJetArrayName(0), fRhoTaskName(),
 105 fRandConeRadius(0.4),fRandConeRadiusSquared(fRandConeRadius*fRandConeRadius), fSignalJetRadius(0.4), fBackgroundJetRadius(0.3), fBackgroundJetPtMin(15.0),
 106 fSignalJetEtaWindow(0.5), fBackgroundJetEtaWindow(0.5), fTrackEtaWindow(0.9), fMinTrackPt(0.150), fMinJetArea(0.5),  fNumberOfCentralityBins(20), fCentralityType("V0A"),
 107 fCrossSection(0.0), fTrials(0.0), fImpParam(-1.0), fRandom(0), fHelperClass(0), fInitialized(0),
 108 fTTlow(8.0), fTThigh(9.0), fTTtype(0), fDphiCut(TMath::Pi()-0.6), fUseDoubleBinPrecision(0),
 109 fHistEvtSelection(0x0), fh2Ntriggers(0x0), fHJetSpec(0x0), fHJetSpecSubUeMedian(0x0), fHJetSpecSubUeCone(0x0), fHJetSpecSubUeCMS(0x0),
 110 fhRhoCellMedian(0x0), fhRhoCone(0x0), fhRhoCMS(0x0),
 111 fhRhoCellMedianIncl(0x0), fhRhoConeIncl(0x0), fhRhoCMSIncl(0x0),
 112 fARhoCellMedian(0x0), fARhoCone(0x0), fARhoCMS(0x0),
 113 fhDeltaPtMedian(0x0), fhDeltaPtCone(0x0), fhDeltaPtCMS(0x0),
 114 fhDeltaPtMedianIncl(0x0), fhDeltaPtConeIncl(0x0), fhDeltaPtCMSIncl(0x0),
 115 fhDeltaPtMedianNearSide(0x0), fhDeltaPtMedianAwaySide(0x0), fhDeltaPtCMSNearSide(0x0), fhDeltaPtCMSAwaySide(0x0),
 116 fhDeltaPtMedianExclTrigCone(0x0),fhDeltaPtCMSExclTrigCone(0x0), fhDeltaPtMedianExclAwayJet(0x0), fhDeltaPtCMSExclAwayJet(0x0),
 117 fhJetPhi(0x0), fhTrackPhi(0x0), fhJetEta(0x0), fhTrackEta(0x0), fhTrackCentVsPt(0x0), fhVertexZ(0x0), fhVertexZAccept(0x0),
 118 fhDphiTriggerJetMinBias(0x0), fhDphiTriggerJetCent20(0x0), fhDphiTriggerJetAccept(0x0),
 119 fhCentrality(0x0), fhCentralityV0M(0x0), fhCentralityV0A(0x0), fhCentralityV0C(0x0), fhCentralityZNA(0x0),
 120 fNofRndTrials(2000), fJetFreeAreaFrac(0.8), fnEta(2), fnPhi(11), fEtaSize(0.9), fPhiSize(2*TMath::Pi()/fnPhi), fCellArea(fPhiSize*fEtaSize),
 121 fh1Xsec(0x0), fh1Trials(0x0), fh1PtHard(0x0), fhImpactParameter(0x0), fhImpactParameterTT(0x0),
 122 fNofRandomCones(1),
 123 fRConesR(0.1),fRConesRSquared(fRConesR*fRConesR), fnRCones(16)
 124 {
 125    //constructor that is called
 126    //LIST OF TRACKS
 127    fTrackArrayName = new TString(trackArrayName);
 128    if((fTrackArrayName->Contains("MC") && fTrackArrayName->Contains("Particles")) ||
 129       (fTrackArrayName->Contains("mc") && fTrackArrayName->Contains("particles"))){
 130       fIsKinematics = kTRUE;
 131    }
 132
 133    //LIST of JETS
 134    fJetArrayName = new TString(jetArrayName);
 135    if(strcmp(fJetArrayName->Data(),"") == 0){
 136       AliError(Form("%s: Jet branch missing !", GetName()));
 137    }
 138
 139    //LIST OF JETS TO BE IGNORED WHILE RHO ESTIMATE
 140    fBackgroundJetArrayName = new TString(backgroundJetArrayName); //jets to be removed from cell median rho estimate
 141    if(strcmp(fBackgroundJetArrayName->Data(),"") == 0){
 142       AliError(Form("%s: Bg Jet branch missing !", GetName()));
 143    }
 144
 145    for(Int_t k=0; k<50; k++){
 146       fRConePhi[k] = 0.0;
 147       fRConeEta[k] = 0.0;
 148    }
 149
 150    DefineOutput(1, TList::Class());
 151 }
 152
 153 //________________________________________________________________________
 154 void  AliAnalysisTaskHJetSpectra::SetAnalyzeMC(Int_t val){
 155    if(val==1){
 156       fAnalyzePythia = kTRUE;
 157       fAnalyzeHijing = kFALSE;
 158       return;
 159    }
 160    if(val==2){
 161       fAnalyzeHijing = kTRUE;
 162       fAnalyzePythia = kFALSE;
 163       return;
 164    }
 165
 166    fAnalyzeHijing = kFALSE;
 167    fAnalyzePythia = kFALSE;
 168    return;
 169 }
 170 //________________________________________________________________________
 171 Double_t AliAnalysisTaskHJetSpectra::GetConePt(Double_t eta, Double_t phi, Double_t radius){
 172    //sum up pt inside a cone
 173    Double_t tmpConePt = 0.0;
 174    Double_t dphi      = 0.0;
 175    Double_t deta      = 0.0;
 176    Double_t radiussquared = radius*radius;
 177
 178    for(Int_t i = 0; i < fTrackArray->GetEntries(); i++){
 179       AliVTrack* tmpTrack = static_cast<AliVTrack*>(fTrackArray->At(i));
 180       if(!tmpTrack) continue;
 181       if(IsTrackInAcceptance(tmpTrack)){
 182          dphi = RelativePhi(tmpTrack->Phi(),phi);
 183          deta = tmpTrack->Eta() - eta;
 184          if( dphi*dphi + deta*deta < radiussquared ){
 185             tmpConePt = tmpConePt + tmpTrack->Pt();
 186          }
 187       }
 188    }
 189    return tmpConePt;
 190 }
 191
 192
 193 //________________________________________________________________________
 194 Double_t AliAnalysisTaskHJetSpectra::GetPtHard(){
 195    //get pt hard from pythia
 196    AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 197    if(MCEvent()){
 198       if(!pythiaHeader){
 199          // Check if AOD
 200          AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 201
 202          if(aodMCH){
 203             for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++){
 204                pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
 205                if(pythiaHeader) break;
 206             }
 207          }
 208       }
 209    }
 210    if(pythiaHeader){
 211
 212       TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 213
 214       if(tree){
 215          TFile *curfile = tree->GetCurrentFile();
 216          if(!curfile) {
 217             Error("Notify","No current file");
 218             return kFALSE;
 219          }
 220          Float_t xsection = 0.0;
 221          Float_t trials  = 1.0;
 222
 223          AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,trials);
 224
 225          fCrossSection = (Double_t) xsection;//pythiaHeader->GetXsection();
 226
 227          if(fCrossSection>0.){ //save cross-section and the number of trials
 228             fTrials = (Double_t) trials; //pythiaHeader->Trials();
 229             fh1Xsec->Fill("<#sigma>", fCrossSection);
 230             fh1Trials->Fill("#sum{ntrials}",fTrials);
 231          }
 232       }
 233       return pythiaHeader->GetPtHard();
 234    }
 235    AliWarning(Form("In task %s: GetPtHard() failed!", GetName()));
 236    return -1.0;
 237 }
 238 //________________________________________________________________________
 239 Double_t AliAnalysisTaskHJetSpectra::GetImpactParameter(){
 240    //get impact parameter from hijing
 241    AliGenHijingEventHeader* hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader());
 242    if(MCEvent()){
 243       if(!hijingHeader){
 244          // Check if AOD
 245          AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 246
 247          if(aodMCH){
 248             for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++){
 249                hijingHeader = dynamic_cast<AliGenHijingEventHeader*>(aodMCH->GetCocktailHeader(i));
 250                if(hijingHeader) break;
 251             }
 252          }
 253       }
 254    }
 255    if(hijingHeader){
 256       return (Double_t) (hijingHeader->ImpactParameter());
 257    }
 258    AliWarning(Form("In task %s: GetImpactParameter() failed!", GetName()));
 259    return -1.0;
 260 }
 261 //________________________________________________________________________
 262 Double_t AliAnalysisTaskHJetSpectra::GetSimPrimaryVertex(){
 263    //get generator level primary vertex
 264    AliGenEventHeader* mcHeader = NULL;
 265    AliAODMCHeader* aodMCH = NULL;
 266    if(MCEvent()){
 267       if(fAnalyzePythia){
 268          mcHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 269          if(!mcHeader){
 270             // Check if AOD
 271              aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 272
 273              if(aodMCH){
 274                 for(UInt_t i = 0; i<aodMCH->GetNCocktailHeaders(); i++){
 275                   mcHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
 276                   if(mcHeader) break;
 277                }
 278             }
 279          }
 280       }
 281
 282       if(fAnalyzeHijing){
 283          mcHeader = dynamic_cast<AliGenHijingEventHeader*>(MCEvent()->GenEventHeader());
 284          if(!mcHeader){
 285             // Check if AOD
 286              aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 287
 288              if(aodMCH){
 289                 for(UInt_t i = 0; i<aodMCH->GetNCocktailHeaders(); i++){
 290                   mcHeader = dynamic_cast<AliGenHijingEventHeader*>(aodMCH->GetCocktailHeader(i));
 291                   if(mcHeader) break;
 292                }
 293             }
 294          }
 295       }
 296    }
 297    if(mcHeader){
 298
 299       TArrayF pyVtx(3);
 300       mcHeader->PrimaryVertex(pyVtx);
 301       return (Double_t) (pyVtx[2]);
 302    }
 303    AliWarning(Form("In task %s: Pythia Vertex failed!", GetName()));
 304    return 9999.0;
 305 }
 306
 307
 308
 309 //________________________________________________________________________
 310 /*Double_t AliAnalysisTaskHJetSpectra::GetPythiaTrials()
 311 {
 312   #ifdef DEBUGMODE
 313     AliInfo("Starting GetPythiaTrials.");
 314   #endif
 315   AliGenPythiaEventHeader* pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 316   if (MCEvent())
 317     if (!pythiaHeader)
 318     {
 319       // Check if AOD
 320       AliAODMCHeader* aodMCH = dynamic_cast<AliAODMCHeader*>(InputEvent()->FindListObject(AliAODMCHeader::StdBranchName()));
 321
 322       if (aodMCH)
 323       {
 324         for(UInt_t i = 0;i<aodMCH->GetNCocktailHeaders();i++)
 325         {
 326           pythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(aodMCH->GetCocktailHeader(i));
 327           if (pythiaHeader) break;
 328         }
 329       }
 330     }
 331
 332   #ifdef DEBUGMODE
 333     AliInfo("Ending GetPythiaTrials.");
 334   #endif
 335   if (pythiaHeader)
 336     return pythiaHeader->Trials();
 337
 338   AliWarning(Form("In task %s: GetPythiaTrials() failed!", GetName()));
 339   return -1.0;
 340 }
 341 */
 342
 343 //________________________________________________________________________
 344 Double_t AliAnalysisTaskHJetSpectra::GetExternalRho(){
 345    // Get rho from event using CMS approach
 346
 347    AliRhoParameter *rho = 0;
 348    if(!fRhoTaskName.IsNull()) {
 349       rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoTaskName.Data()));
 350       if(!rho){
 351          AliWarning(Form("%s: Could not retrieve rho with name %s!", GetName(), fRhoTaskName.Data()));
 352          return 0.0;
 353       }
 354    }else return 0.0;
 355
 356    return (rho->GetVal());
 357 }
 358
 359 //________________________________________________________________________
 360 Bool_t AliAnalysisTaskHJetSpectra::IsEventInAcceptance(AliVEvent* event){
 361    //EVENT SELECTION
 362
 363
 364    if(!event) return kFALSE;
 365
 366    //___________________________________________________
 367
 368    if(fAnalyzePythia || fAnalyzeHijing){ //PURE MC
 369       if(!MCEvent()) return kFALSE;
 370
 371        //BEFORE VERTEX CUT
 372       Double_t vtxMC = GetSimPrimaryVertex();
 373       fhVertexZ->Fill(vtxMC);
 374
 375       if(TMath::Abs(vtxMC) > 10.0){
 376          fHistEvtSelection->Fill(3); //count events rejected by vertex cut
 377          return kFALSE;
 378       }
 379       fhVertexZAccept->Fill(vtxMC);
 380
 381       return kTRUE;
 382    }
 383    //___________________________________________________
 384    //TEST PILE UP
 385    if(fUsePileUpCut){
 386       if(!fHelperClass || fHelperClass->IsPileUpEvent(event)){
 387          fHistEvtSelection->Fill(2); //count events rejected by pileup
 388          return kFALSE;
 389       }
 390    }
 391    //___________________________________________________
 392    //BEFORE VERTEX CUT
 393    fhVertexZ->Fill(event->GetPrimaryVertex()->GetZ());
 394
 395    if(fUseDefaultVertexCut){
 396       if(!fHelperClass || !fHelperClass->IsVertexSelected2013pA(event)){
 397          fHistEvtSelection->Fill(3); //count events rejected by vertex cut
 398          return kFALSE;
 399       }
 400    }else{
 401       if(TMath::Abs(event->GetPrimaryVertex()->GetZ()) > 10.0){
 402          fHistEvtSelection->Fill(3); //count events rejected by vertex cut
 403          return kFALSE;
 404       }
 405    }
 406    //___________________________________________________
 407    //AFTER VERTEX CUT
 408    fhVertexZAccept->Fill(event->GetPrimaryVertex()->GetZ());
 409
 410   return kTRUE;
 411 }
 412
 413 //________________________________________________________________________
 414 Bool_t AliAnalysisTaskHJetSpectra::IsTrackInAcceptance(AliVParticle* track){
 415    // Check if the track pt and eta range
 416    if(track != 0){
 417       if(fIsKinematics){
 418          // TODO: Only working for AOD MC
 419          if((!track->Charge()) || (!(static_cast<AliAODMCParticle*>(track))->IsPhysicalPrimary()) )
 420             return kFALSE;
 421       }
 422       if(TMath::Abs(track->Eta()) <= fTrackEtaWindow){ //APPLY TRACK ETA CUT
 423          if(track->Pt() >= fMinTrackPt){   //APPLY TRACK CUT
 424             return kTRUE;
 425          }
 426       }
 427    }
 428    return kFALSE;
 429 }
 430
 431 //________________________________________________________________________
 432 Bool_t AliAnalysisTaskHJetSpectra::IsBackgroundJetInAcceptance(AliEmcalJet *jet){
 433    //find jets to be removed from bg calculation
 434    if(jet != 0){
 435       if(TMath::Abs(jet->Eta()) <= fBackgroundJetEtaWindow){
 436          if(jet->Pt() >= fBackgroundJetPtMin){ //accept only hard jets
 437             return kTRUE;
 438          }
 439       }
 440    }
 441    return kFALSE;
 442 }
 443
 444 //________________________________________________________________________
 445 Bool_t AliAnalysisTaskHJetSpectra::IsSignalJetInAcceptance(AliEmcalJet *jet){
 446    //select jets in acceptance
 447    if(jet == 0) return kFALSE;
 448    if(TMath::Abs(jet->Eta()) <= fSignalJetEtaWindow){
 449       if(jet->Pt() >= fMinTrackPt){
 450          if(jet->Area() >= fMinJetArea){
 451             return kTRUE;
 452          }
 453       }
 454    }
 455    return kFALSE;
 456 }
 457
 458
 459 //________________________________________________________________________
 460 void AliAnalysisTaskHJetSpectra::ExecOnce(){
 461    //Read arrays of jets and tracks
 462
 463
 464    fInitialized = kTRUE; //change flag to skip this function next time when processing UserExec
 465
 466
 467    fnRCones = TMath::Nint(fRandConeRadiusSquared/fRConesRSquared); //the number of small R=0.1 random cones
 468
 469    // Check for track array
 470    if(strcmp(fTrackArrayName->Data(), "") != 0){
 471       fTrackArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrayName->Data()));
 472       if(!fTrackArray){
 473          AliWarning(Form("%s: Could not retrieve tracks %s!", GetName(), fTrackArrayName->Data()));
 474       }else{
 475          TClass *cl = fTrackArray->GetClass();
 476          if(!cl->GetBaseClass("AliVParticle")){
 477             AliError(Form("%s: Collection %s does not contain AliVParticle objects!", GetName(), fTrackArrayName->Data()));
 478             fTrackArray = 0;
 479          }
 480       }
 481    }
 482
 483    // Check for jet array
 484    if(strcmp(fJetArrayName->Data(), "") != 0){
 485       fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrayName->Data()));
 486
 487       if(!fJetArray){
 488          AliWarning(Form("%s: Could not retrieve jets %s!", GetName(), fJetArrayName->Data()));
 489       }else{
 490          if(!fJetArray->GetClass()->GetBaseClass("AliEmcalJet")){
 491             AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fJetArrayName->Data()));
 492             fJetArray = 0;
 493          }
 494       }
 495    }
 496
 497    // Check for list of jets to be removed from background
 498    if(strcmp(fBackgroundJetArrayName->Data(), "") != 0){
 499       fBackgroundJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fBackgroundJetArrayName->Data()));
 500       if(!fBackgroundJetArray){
 501          AliInfo(Form("%s: Could not retrieve background jets %s!", GetName(), fBackgroundJetArrayName->Data()));
 502       }else{
 503          if(!fBackgroundJetArray->GetClass()->GetBaseClass("AliEmcalJet")){
 504             AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fBackgroundJetArrayName->Data()));
 505             fBackgroundJetArray = 0;
 506          }
 507       }
 508    }
 509
 510    // Look, if initialization is OK
 511
 512    // Initialize helper class (for vertex selection & pile up correction)
 513    fHelperClass = new AliAnalysisUtils();
 514    fHelperClass->SetCutOnZVertexSPD(kFALSE); // kFALSE: no cut; kTRUE: |zvtx-SPD - zvtx-TPC|<0.5cm
 515
 516    return;
 517 }
 518
 519
 520 //________________________________________________________________________
 521 void  AliAnalysisTaskHJetSpectra::GetDeltaPt(Double_t rho1, Double_t &dpt1, Double_t rho2, Double_t &dpt2,
 522                                            Double_t rho3, Double_t &dpt3,
 523                                            Double_t &rcPhi, Double_t &rcEta,
 524                                            Double_t leadingJetExclusionProbability){
 525
 526    //delta pt = random cone - rho
 527
 528    // Define an invalid delta pt
 529    dpt1 = -10000.0;
 530    dpt2 = -10000.0;
 531    dpt3 = -10000.0;
 532
 533    // Define eta range
 534    Double_t etaMin, etaMax;
 535    etaMin = -(fTrackEtaWindow-fRandConeRadius);
 536    etaMax = +(fTrackEtaWindow-fRandConeRadius);
 537
 538    // Define random cone Eta+Phi
 539    Bool_t coneValid = kTRUE;
 540    Double_t tmpRandConeEta = etaMin + fRandom->Rndm()*(etaMax-etaMin);
 541    Double_t tmpRandConePhi = fRandom->Rndm()*TMath::TwoPi();
 542
 543    // if there is a jet, check for overlap if demanded
 544    if(leadingJetExclusionProbability){
 545       AliEmcalJet* tmpLeading = NULL;
 546       Double_t lpt = -1.0;
 547       // Get leading jet (regardless of pT)
 548       for(Int_t i = 0; i<fJetArray->GetEntries(); i++){
 549          AliEmcalJet* tmpJet = static_cast<AliEmcalJet*>(fJetArray->At(i));
 550          if(!tmpJet) continue;
 551          if((TMath::Abs(tmpJet->Eta()) <= fSignalJetEtaWindow) && (tmpJet->Area() >= fMinJetArea)){
 552             if(tmpJet->Pt() > lpt){
 553                tmpLeading = tmpJet;
 554                lpt =  tmpJet->Pt();
 555             }
 556          }
 557       }
 558       if(tmpLeading){
 559          Double_t excludedJetPhi = tmpLeading->Phi();
 560          Double_t tmpDeltaPhi    = RelativePhi(tmpRandConePhi, excludedJetPhi);
 561          Double_t excludedJetEta = tmpLeading->Eta()-tmpRandConeEta;
 562
 563          // Check, if cone has overlap with jet
 564          if(tmpDeltaPhi*tmpDeltaPhi + excludedJetEta*excludedJetEta <= fRandConeRadiusSquared){
 565             // Define probability to exclude the RC
 566             Double_t probability = leadingJetExclusionProbability;
 567
 568             // Only exclude cone with a given probability
 569             if(fRandom->Rndm()<=probability)  coneValid = kFALSE;
 570          }
 571       }
 572    }
 573
 574    rcPhi = 9999.0;
 575    rcEta = 9999.0;
 576
 577    // Get the cones' pt and calculate delta pt
 578    if(coneValid){
 579       rcPhi = tmpRandConePhi;
 580       rcEta = tmpRandConeEta;
 581       Double_t conePt = GetConePt(tmpRandConeEta,tmpRandConePhi,fRandConeRadius);
 582       dpt1 =  conePt - (rho1*fRandConeRadiusSquared*TMath::Pi());
 583       dpt2 =  conePt - (rho2*fRandConeRadiusSquared*TMath::Pi());
 584       dpt3 =  conePt - (rho3*fRandConeRadiusSquared*TMath::Pi());
 585    }
 586
 587 }
 588
 589
 590 //________________________________________________________________________
 591 void AliAnalysisTaskHJetSpectra::Calculate(AliVEvent* event){
 592    //Analyze the event and Fill histograms
 593
 594    if(fAnalyzePythia){
 595       fh1PtHard->Fill(GetPtHard());
 596    }
 597
 598    if(fAnalyzeHijing){
 599       fImpParam = GetImpactParameter();
 600       fhImpactParameter->Fill(fImpParam);
 601    }
 602    //_________________________________________________________________
 603    //  FILL EVENT STATISTICS
 604    fHistEvtSelection->Fill(1); //Count input event
 605
 606    if(!IsEventInAcceptance(event)) return; //post data is in UserExec
 607
 608
 609    // Get centrality
 610    AliCentrality* tmpCentrality = event->GetCentrality();
 611    if(!tmpCentrality){
 612       fHistEvtSelection->Fill(4);
 613       return; //post data is in UserExec
 614    }
 615    Double_t centralityPercentile    = -1.0;
 616    Double_t centralityPercentileV0A = -1.0;
 617    Double_t centralityPercentileV0C = -1.0;
 618    Double_t centralityPercentileV0M = -1.0;
 619    Double_t centralityPercentileZNA = -1.0;
 620    if(tmpCentrality != NULL){
 621       centralityPercentile    = tmpCentrality->GetCentralityPercentile(fCentralityType.Data());
 622       centralityPercentileV0A = tmpCentrality->GetCentralityPercentile("V0A");
 623       centralityPercentileV0C = tmpCentrality->GetCentralityPercentile("V0C");
 624       centralityPercentileV0M = tmpCentrality->GetCentralityPercentile("V0M");
 625       centralityPercentileZNA = tmpCentrality->GetCentralityPercentile("ZNA");
 626    }
 627
 628    if((centralityPercentile < 0.0) || (centralityPercentile > 100.0)){
 629       AliWarning(Form("Centrality value not valid (c=%E)",centralityPercentile));
 630       fHistEvtSelection->Fill(4);
 631       return;
 632    }
 633    fhCentrality->Fill(centralityPercentile);
 634    fhCentralityV0M->Fill(centralityPercentileV0M);
 635    fhCentralityV0A->Fill(centralityPercentileV0A);
 636    fhCentralityV0C->Fill(centralityPercentileV0C);
 637    fhCentralityZNA->Fill(centralityPercentileZNA);
 638
 639    fHistEvtSelection->Fill(0); //Count input event
 640
 641    // END EVENT SELECTION
 642    //___________________________________________________________
 643
 644    //LOOP OVER TRACKS  SEARCH FOR TRIGGER
 645    std::vector<Int_t> trigTracks; //list pf trigger particle indices
 646    //Bool_t bContainesHighPtTrack = kFALSE;
 647
 648    Int_t nTracks = fTrackArray->GetEntries();
 649
 650    for(Int_t i = 0; i < nTracks; i++){
 651       AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(i));
 652
 653       if(!track) continue;
 654
 655       if(IsTrackInAcceptance(track)){
 656          fhTrackPhi->Fill(track->Pt(), RelativePhi(track->Phi(),0.0)); // phi = -pi az pi
 657          fhTrackEta->Fill(track->Pt(), track->Eta());
 658          fhTrackCentVsPt->Fill(track->Pt(), centralityPercentile);
 659
 660          if(fTTlow <= track->Pt() && track->Pt() < fTThigh){
 661             trigTracks.push_back(i);  //trigger candidates
 662          }
 663
 664          //if(track->Pt()>=8.0) bContainesHighPtTrack = kTRUE;
 665       }
 666    }
 667
 668    Int_t ntriggers = (Int_t) trigTracks.size();
 669    Int_t indexSingleRndTrig = -1; //index of single random trigger
 670    Double_t areaJet,  pTJet;
 671    Double_t tmpArray[3];
 672    Double_t rhoFromCellMedian = 0.0; //UE density cell median
 673    Double_t rhoCone           = 0.0; //UE density perp cone
 674    Double_t rhoCMS            = 0.0; //UE density ala CMS
 675    Double_t deltaptCellMedian, deltaptCone, deltaptCMS, randConePhi, randConeEta;
 676    Double_t distanceFromTrigger;
 677
 678    if(ntriggers>0){
 679       if(fTTtype==0){ //select single inclusive trigger
 680          indexSingleRndTrig = fRandom->Integer(ntriggers); //Integer 0 ... ntriggers-1
 681       }
 682    }
 683
 684    rhoFromCellMedian = EstimateBgRhoMedian();
 685    rhoCone           = EstimateBgCone();
 686    rhoCMS            = GetExternalRho();
 687
 688    fhRhoCellMedianIncl->Fill((Float_t) rhoFromCellMedian,(Float_t) centralityPercentile);
 689    fhRhoConeIncl->Fill(      (Float_t) rhoCone,          (Float_t) centralityPercentile);
 690    fhRhoCMSIncl->Fill(       (Float_t) rhoCMS,           (Float_t) centralityPercentile);
 691
 692    for(Int_t irc=0; irc<fNofRandomCones; irc++){ //generate 4 random cones per event
 693       GetDeltaPt(rhoFromCellMedian, deltaptCellMedian,rhoCone, deltaptCone, rhoCMS, deltaptCMS, randConePhi, randConeEta, 0);
 694
 695       fhDeltaPtMedianIncl->Fill(deltaptCellMedian, (Double_t) centralityPercentile);
 696       fhDeltaPtConeIncl->Fill( deltaptCone,        (Double_t) centralityPercentile);
 697       fhDeltaPtCMSIncl->Fill( deltaptCMS,          (Double_t) centralityPercentile);
 698
 699       if(ntriggers>0){
 700          //fill delta pt histograms near side + away side
 701          fhDeltaPtMedian->Fill( deltaptCellMedian, (Double_t) centralityPercentile);
 702          fhDeltaPtCone->Fill( deltaptCone,         (Double_t) centralityPercentile);
 703          fhDeltaPtCMS->Fill(  deltaptCMS,          (Double_t) centralityPercentile);
 704
 705          if(indexSingleRndTrig>-1){
 706             AliVTrack* triggHad = static_cast<AliVTrack*>(fTrackArray->At(trigTracks[indexSingleRndTrig]));
 707             Double_t dphiTrigRC =  RelativePhi(triggHad->Phi(), randConePhi);
 708             Double_t detaTrigRC =  triggHad->Eta()- randConeEta;
 709             if(TMath::Abs(dphiTrigRC)< TMath::Pi()/2){ //near side
 710                fhDeltaPtMedianNearSide->Fill( deltaptCellMedian, (Double_t) centralityPercentile);
 711                fhDeltaPtCMSNearSide->Fill(  deltaptCMS,          (Double_t) centralityPercentile);
 712             }else{ //away side
 713                fhDeltaPtMedianAwaySide->Fill( deltaptCellMedian, (Double_t) centralityPercentile);
 714                fhDeltaPtCMSAwaySide->Fill(  deltaptCMS,          (Double_t) centralityPercentile);
 715             }
 716
 717             distanceFromTrigger = sqrt(dphiTrigRC*dphiTrigRC+detaTrigRC*detaTrigRC);
 718             while(distanceFromTrigger<0.5 + fRandConeRadius){
 719                GetDeltaPt(rhoFromCellMedian, deltaptCellMedian,rhoCone, deltaptCone, rhoCMS, deltaptCMS, randConePhi, randConeEta, 0);
 720                dphiTrigRC =  RelativePhi(triggHad->Phi(), randConePhi);
 721                detaTrigRC =  triggHad->Eta()- randConeEta;
 722                distanceFromTrigger = sqrt(dphiTrigRC*dphiTrigRC+detaTrigRC*detaTrigRC);
 723             }
 724             if(distanceFromTrigger>0.5 + fRandConeRadius){
 725                fhDeltaPtMedianExclTrigCone->Fill( deltaptCellMedian, (Double_t) centralityPercentile);
 726                fhDeltaPtCMSExclTrigCone->Fill(  deltaptCMS,          (Double_t) centralityPercentile);
 727             }
 728          }
 729       }
 730    }
 731
 732    //_______________________________________
 733    Int_t    idxLeadingJetAwaySide = -1;
 734    Double_t ptLeadingJetAwaySide  = -1.0;
 735    Double_t phiTrigger            = -1.0; //-pi,pi
 736
 737    if(ntriggers>0){
 738       //Estimate UE density
 739       //Fill once per event
 740       fhRhoCellMedian->Fill((Float_t) rhoFromCellMedian,(Float_t) centralityPercentile);
 741       fhRhoCone->Fill(      (Float_t) rhoCone,          (Float_t) centralityPercentile);
 742       fhRhoCMS->Fill(       (Float_t) rhoCMS,           (Float_t) centralityPercentile);
 743
 744
 745       //TRIGGER PARTICLE LOOP
 746       for(Int_t it=0; it<ntriggers; it++){ //loop over trigger configurations
 747
 748          if(fTTtype==0){
 749             if(it != indexSingleRndTrig) continue;
 750          }
 751
 752          AliVTrack* triggerHadron = static_cast<AliVTrack*>(fTrackArray->At(trigTracks[it]));
 753          if(!triggerHadron) continue;
 754
 755          fh2Ntriggers->Fill((Float_t) centralityPercentile, (Float_t) triggerHadron->Pt()); //trigger p
 756
 757          if(fAnalyzeHijing){ //impact parameter for triggered events
 758             fhImpactParameterTT->Fill(fImpParam);
 759          }
 760
 761          phiTrigger = RelativePhi(triggerHadron->Phi(),0.0); //-pi,pi
 762
 763          //JET LOOP
 764          for(Int_t ij = 0; ij < fJetArray->GetEntries(); ij++){
 765             AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(ij));
 766             if(!jet){
 767                AliError(Form("%s: Could not receive jet %d", GetName(), ij));
 768                continue;
 769             }
 770             if(!IsSignalJetInAcceptance(jet)) continue;
 771
 772             areaJet = jet->Area();
 773             pTJet   = jet->Pt();
 774
 775             if(it==0 || it == indexSingleRndTrig){
 776                fhJetPhi->Fill( pTJet, RelativePhi(jet->Phi(),0.0));
 777                fhJetEta->Fill( pTJet, jet->Eta());
 778             }
 779
 780             Double_t dphi = RelativePhi(triggerHadron->Phi(), jet->Phi());
 781
 782             Double_t dfi = dphi; //-0.5*pi to 1.5*Pi
 783             if(dfi<-0.5*TMath::Pi()) dfi += 2*TMath::Pi();
 784             if(dfi> 1.5*TMath::Pi()) dfi -= 2*TMath::Pi();
 785             fhDphiTriggerJetMinBias->Fill((Float_t) jet->Pt(),(Float_t) dfi);
 786             if(centralityPercentile<20.) fhDphiTriggerJetCent20->Fill((Float_t) jet->Pt(),(Float_t) dfi);
 787             //-------------------------
 788
 789             if(TMath::Abs(dphi) < fDphiCut) continue;  //Dphi cut between trigger and assoc
 790             fhDphiTriggerJetAccept->Fill(dfi); //Accepted
 791
 792             if(pTJet > ptLeadingJetAwaySide){ //search for the leading away side jet
 793                idxLeadingJetAwaySide = ij;
 794                ptLeadingJetAwaySide  = pTJet;
 795             }
 796
 797            //Centrality, A, pTjet
 798             tmpArray[0] =  centralityPercentile;
 799             tmpArray[1] =  areaJet;
 800             tmpArray[2] =  pTJet;
 801             fHJetSpec->Fill(tmpArray);
 802
 803             //Subtract cell median
 804             tmpArray[2] = pTJet - areaJet*rhoFromCellMedian;
 805             fHJetSpecSubUeMedian->Fill(tmpArray);
 806             fARhoCellMedian->Fill((Float_t) (areaJet*rhoFromCellMedian));
 807
 808             //Subtract perp cone
 809             tmpArray[2] = pTJet - areaJet*rhoCone;
 810             fHJetSpecSubUeCone->Fill(tmpArray);
 811             fARhoCone->Fill((Float_t) (areaJet*rhoCone));
 812
 813             //Subtract CMS bg
 814             tmpArray[2] = pTJet - areaJet*rhoCMS;
 815             fHJetSpecSubUeCMS->Fill(tmpArray);
 816             fARhoCMS->Fill((Float_t) (areaJet*rhoCMS));
 817
 818          }//JET LOOP
 819       }
 820    }
 821
 822    //_______________________________________
 823    // Get delta phi from small R=0.1 cones
 824    if(ntriggers>0 && indexSingleRndTrig>-1){
 825
 826       AliEmcalJet* jet = NULL;
 827       Double_t phiExclJet =0., etaExclJet = 9999., rExclJet = 0.0;
 828
 829       if(idxLeadingJetAwaySide>-1){
 830          jet = static_cast<AliEmcalJet*>(fJetArray->At(idxLeadingJetAwaySide));
 831          if(!jet){
 832             AliError(Form("%s: Could not receive leading jet %d", GetName(), idxLeadingJetAwaySide));
 833          }else{
 834              phiExclJet = jet->Phi();
 835              etaExclJet = jet->Eta();
 836              rExclJet   = TMath::Sqrt(jet->Area()/TMath::Pi());
 837          }
 838       }
 839
 840       Int_t countPlacedRcones = 0;
 841       Int_t inwhile=0;
 842       Double_t dphiMaxFromPiForRcones = TMath::Pi() - fDphiCut + fRandConeRadius - fRConesR;
 843       Double_t detaMaxForRcones = fTrackEtaWindow - fRConesR;
 844       Double_t rcphi,rceta;
 845       Bool_t goodrc;
 846
 847       while(countPlacedRcones < fnRCones){ //generate fnRCones cones of radius R=0.1
 848          inwhile++;
 849          if(inwhile>500){
 850             AliError(Form("%s: Small space where to put another random cone %d", GetName(), idxLeadingJetAwaySide));
 851             break;
 852          }
 853          rcphi = RelativePhi(phiTrigger + TMath::Pi() + dphiMaxFromPiForRcones*fRandom->Uniform(-1.0,1.0), 0.0); //-pi,pi
 854          rceta = detaMaxForRcones*fRandom->Uniform(-1.0,1.0);
 855          if(jet){//do not merge random cones with the leading jet
 856             if(!DistantCones(rcphi,rceta,fRConesR, phiExclJet, etaExclJet, rExclJet)) continue;
 857          }
 858
 859          goodrc = kTRUE; //generate disjoint random cones
 860          for(int k=0; k<countPlacedRcones; k++){
 861             if(!DistantCones(rcphi, rceta, fRConesR, (Double_t) fRConePhi[k], (Double_t) fRConeEta[k], fRConesR)){
 862                goodrc = kFALSE;
 863                break;
 864             }
 865          }//end loop over already placed cones
 866
 867          if(!goodrc) continue;
 868          fRConePhi[countPlacedRcones] = rcphi;
 869          fRConeEta[countPlacedRcones] = rceta;
 870          countPlacedRcones++;
 871       }//end of  loop generating small R=0.1 random cones
 872       //sum track pT in the random cones
 873       Double_t sumPtofRandomCones = 0.0;
 874
 875       for(Int_t itrk = 0; itrk < nTracks; itrk++){
 876          AliVTrack* track = static_cast<AliVTrack*>(fTrackArray->At(itrk));
 877          if(!track) continue;
 878
 879          if(IsTrackInAcceptance(track)){
 880             for(int k=0; k<countPlacedRcones; k++){
 881
 882                rcphi = RelativePhi(track->Phi(), fRConePhi[k]); // phi = -pi az pi
 883                rceta = track->Eta() - fRConeEta[k];
 884                if(rcphi*rcphi + rceta*rceta < fRConesRSquared){
 885                   sumPtofRandomCones += track->Pt(); //track is in the cone
 886                   break;
 887                }
 888             }//loop over cones
 889          }
 890       }//loop over tracks
 891       Double_t totarea = countPlacedRcones*TMath::Pi()*fRConesRSquared;
 892       fhDeltaPtMedianExclAwayJet->Fill( sumPtofRandomCones - rhoFromCellMedian*totarea, (Double_t) centralityPercentile );
 893       fhDeltaPtCMSExclAwayJet->Fill(    sumPtofRandomCones - rhoCMS*totarea , (Double_t) centralityPercentile );
 894
 895    }//end delta phi from small R=0.1 random cones
 896
 897    return;
 898 }
 899
 900 //________________________________________________________________________
 901 Bool_t AliAnalysisTaskHJetSpectra::UserNotify(){
 902    // Implemented Notify() to read the cross sections
 903    // and number of trials from pyxsec.root
 904    /*
 905    if(fAnalyzePythia){
 906
 907      TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 908      TFile *currFile = tree->GetCurrentFile();
 909
 910      TString file(currFile->GetName());
 911
 912      if(file.Contains("root_archive.zip#")){
 913        Ssiz_t pos1 = file.Index("root_archive",12,TString::kExact);
 914        Ssiz_t pos = file.Index("#",1,pos1,TString::kExact);
 915        file.Replace(pos+1,20,"");
 916      }
 917      else {
 918        // not an archive take the basename....
 919        file.ReplaceAll(gSystem->BaseName(file.Data()),"");
 920      }
 921
 922      TFile *fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec.root")); // problem that we cannot really test the existance of a file in a archive so we have to lvie with open error message from root
 923      if(!fxsec){
 924        // next trial fetch the histgram file
 925        fxsec = TFile::Open(Form("%s%s",file.Data(),"pyxsec_hists.root"));
 926        if(!fxsec){
 927            // not a severe condition but inciate that we have no information
 928          return kFALSE;
 929        }
 930        else{
 931          // find the tlist we want to be independtent of the name so use the Tkey
 932          TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0);
 933          if(!key){
 934            fxsec->Close();
 935            return kFALSE;
 936          }
 937          TList *list = dynamic_cast<TList*>(key->ReadObj());
 938          if(!list){
 939            fxsec->Close();
 940            return kFALSE;
 941          }
 942          fCrossSection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
 943          fTrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
 944          fxsec->Close();
 945        }
 946      } // no tree pyxsec.root
 947      else {
 948        TTree *xtree = (TTree*)fxsec->Get("Xsection");
 949        if(!xtree){
 950          fxsec->Close();
 951          return kFALSE;
 952        }
 953        UInt_t   ntrials  = 0;
 954        Double_t  xsection  = 0;
 955        xtree->SetBranchAddress("xsection",&xsection);
 956        xtree->SetBranchAddress("ntrials",&ntrials);
 957        xtree->GetEntry(0);
 958        fTrials = ntrials;
 959        fCrossSection = xsection;
 960        fxsec->Close();
 961      }
 962
 963
 964      fh1Xsec->Fill("<#sigma>", fCrossSection);
 965      fh1Trials->Fill("#sum{ntrials}",fTrials);
 966
 967    }
 968    */
 969    return kTRUE;
 970 }
 971
 972 //________________________________________________________________________
 973 void AliAnalysisTaskHJetSpectra::Terminate(Option_t *){
 974    //Treminate
 975    PostData(1, fOutputList);
 976
 977    // Mandatory
 978    fOutputList = dynamic_cast<TList*> (GetOutputData(1)); // '1' refers to the output slot
 979    if(!fOutputList) {
 980       printf("ERROR: Output list not available\n");
 981       return;
 982    }
 983 }
 984
 985 //________________________________________________________________________
 986 AliAnalysisTaskHJetSpectra::~AliAnalysisTaskHJetSpectra(){
 987    // Destructor. Clean-up the output list, but not the histograms that are put inside
 988    // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 989    if(fOutputList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 990       delete fOutputList;
 991    }
 992    delete fRandom;
 993    delete fTrackArrayName;
 994    delete fJetArrayName;
 995    delete fBackgroundJetArrayName;
 996    delete fHelperClass;
 997
 998 }
 999
1000 //________________________________________________________________________
1001 void AliAnalysisTaskHJetSpectra::UserCreateOutputObjects(){
1002   // called once to create user defined output objects like histograms, plots etc.
1003   // and to put it on the output list.
1004   // Note: Saving to file with e.g. OpenFile(0) is must be before creating other objects.
1005
1006    fRandom = new TRandom3(0);
1007
1008    fOutputList = new TList();
1009    fOutputList->SetOwner(); // otherwise it produces leaks in merging
1010    Bool_t oldStatus = TH1::AddDirectoryStatus();
1011    TH1::AddDirectory(kFALSE);
1012
1013    //__________________________________________________________
1014    // Event statistics
1015    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
1016    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
1017    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
1018    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"pile up (rejected)");
1019    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
1020    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
1021
1022    fOutputList->Add(fHistEvtSelection);
1023    //___________________________________________________________
1024    // Hard trigger counter
1025    fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers",
1026                             fNumberOfCentralityBins,0.0,100.0,50,0.0,50.0);
1027    fOutputList->Add(fh2Ntriggers);
1028    //___________________________________________________________
1029    // trigger associated jet spectra (jet pT not corrected for UE)
1030    Int_t bw = (fUseDoubleBinPrecision==0) ? 1 : 2; //make larger bin width
1031
1032    //jet associated to given TT
1033    //Centrality, A, pTjet
1034    const Int_t    dimSpec   = 3;
1035    const Int_t    nBinsSpec[dimSpec]  = {fNumberOfCentralityBins,  50, bw*110};
1036    const Double_t lowBinSpec[dimSpec] = {0.0,                     0.0,  -20.0};
1037    const Double_t hiBinSpec[dimSpec]  = {100.0,                   1.5,  200.0};
1038    fHJetSpec = new THnSparseF("fHJetSpec",
1039                    "Recoil jet spectrum [cent,A,pTjet]",
1040                    dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
1041    fOutputList->Add(fHJetSpec);
1042    //___________________________________________________________
1043    //jet associated to given TT (jet pT corrected with rho from cell median)
1044    fHJetSpecSubUeMedian = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeMedian");
1045    fHJetSpecSubUeMedian->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe]");
1046    fOutputList->Add(fHJetSpecSubUeMedian);
1047    //___________________________________________________________
1048    //jet associated to given TT (jet pT corrected with rho from perp cone)
1049    fHJetSpecSubUeCone = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCone");
1050    fHJetSpecSubUeCone->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe]");
1051    fOutputList->Add(fHJetSpecSubUeCone);
1052    //___________________________________________________________
1053    //jet associated to given TT (jet pT corrected with rho from CMS approach)
1054    fHJetSpecSubUeCMS = (THnSparseF*) fHJetSpec->Clone("fHJetSpecSubUeCMS");
1055    fHJetSpecSubUeCMS->SetTitle("Recoil jet spectrum [cent,A,pTjet-pTUe]");
1056    fOutputList->Add(fHJetSpecSubUeCMS);
1057
1058    //____________________________________________________________________
1059    //UE from cell median  [Centrality, rho, pTUe ]
1060
1061    fhRhoCellMedian = new TH2F("fhRhoCellMedian","Rho",40, 0.0, 20.0, fNumberOfCentralityBins, 0.0, 100.);
1062    fOutputList->Add(fhRhoCellMedian);
1063
1064    fhRhoCone = (TH2F*) fhRhoCellMedian->Clone("fhRhoCone");
1065    fOutputList->Add(fhRhoCone);
1066
1067    fhRhoCMS = (TH2F*) fhRhoCellMedian->Clone("fhRhoCMS");
1068    fOutputList->Add(fhRhoCMS);
1069
1070    fhRhoCellMedianIncl = (TH2F*) fhRhoCellMedian->Clone("fhRhoCellMedianIncl");
1071    fOutputList->Add(fhRhoCellMedianIncl);
1072
1073    fhRhoConeIncl = (TH2F*) fhRhoCellMedian->Clone("fhRhoConeIncl");
1074    fOutputList->Add(fhRhoConeIncl);
1075
1076    fhRhoCMSIncl = (TH2F*) fhRhoCellMedian->Clone("fhRhoCMSIncl");
1077    fOutputList->Add(fhRhoCMSIncl);
1078
1079    //_______________________________________________________________________
1080    // rho times area
1081    fARhoCellMedian = new TH1F("fARhoCellMedian","Area times rho",40, 0.0, 20.0);
1082    fOutputList->Add(fARhoCellMedian);
1083
1084    fARhoCone = (TH1F*) fARhoCellMedian->Clone("fARhoCone");
1085    fOutputList->Add(fARhoCone);
1086
1087    fARhoCMS  = (TH1F*) fARhoCellMedian->Clone("fARhoCMS");
1088    fOutputList->Add(fARhoCMS);
1089
1090    //_______________________________________________________________________
1091    // Delta pt distributions
1092    fhDeltaPtMedian = new TH2D("fhDeltaPtMedian","DeltaPt", bw*110, -20, 200, fNumberOfCentralityBins,0.0,100.0);
1093    fOutputList->Add(fhDeltaPtMedian);
1094
1095    fhDeltaPtCone = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCone");
1096    fOutputList->Add(fhDeltaPtCone);
1097
1098    fhDeltaPtCMS = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMS");
1099    fOutputList->Add(fhDeltaPtCMS);
1100
1101    fhDeltaPtMedianIncl = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianIncl");
1102    fOutputList->Add(fhDeltaPtMedianIncl);
1103
1104    fhDeltaPtConeIncl = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtConeIncl");
1105    fOutputList->Add(fhDeltaPtConeIncl);
1106
1107    fhDeltaPtCMSIncl = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSIncl");
1108    fOutputList->Add(fhDeltaPtCMSIncl);
1109
1110    fhDeltaPtMedianNearSide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianNearSide");
1111    fOutputList->Add(fhDeltaPtMedianNearSide);
1112
1113    fhDeltaPtMedianAwaySide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianAwaySide");
1114    fOutputList->Add(fhDeltaPtMedianAwaySide);
1115
1116    fhDeltaPtCMSNearSide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSNearSide");
1117    fOutputList->Add(fhDeltaPtCMSNearSide);
1118
1119    fhDeltaPtCMSAwaySide= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSAwaySide");
1120    fOutputList->Add(fhDeltaPtCMSAwaySide);
1121
1122    fhDeltaPtMedianExclTrigCone= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianExclTrigCone");
1123    fOutputList->Add(fhDeltaPtMedianExclTrigCone);
1124
1125    fhDeltaPtCMSExclTrigCone= (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSExclTrigCone");
1126    fOutputList->Add(fhDeltaPtCMSExclTrigCone);
1127
1128    fhDeltaPtMedianExclAwayJet = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtMedianExclAwayJet");
1129    fOutputList->Add(fhDeltaPtMedianExclAwayJet);
1130
1131    fhDeltaPtCMSExclAwayJet = (TH2D*) fhDeltaPtMedian->Clone("fhDeltaPtCMSExclAwayJet");
1132    fOutputList->Add(fhDeltaPtCMSExclAwayJet);
1133
1134    //_______________________________________________________________________
1135    //inclusive azimuthal and pseudorapidity histograms
1136    fhJetPhi   = new TH2F("fhJetPhi","Azim dist jets vs pTjet", 50, 0, 100, 50,-TMath::Pi(),TMath::Pi());
1137    fOutputList->Add(fhJetPhi);
1138    //-------------------------
1139    fhTrackPhi = new TH2F("fhTrackPhi","azim dist trig had vs pT,trk", 50, 0, 50, 50,-TMath::Pi(),TMath::Pi());
1140    fOutputList->Add(fhTrackPhi);
1141    //-------------------------
1142    fhJetEta   = new TH2F("fhJetEta","Eta dist jets vs pTjet", 50,0, 100, 40,-0.9,0.9);
1143    fOutputList->Add(fhJetEta);
1144    //-------------------------
1145    fhTrackEta = new TH2F("fhTrackEta","Eta dist trig had vs pT,trk", 50, 0, 50, 40,-0.9,0.9);
1146    fOutputList->Add(fhTrackEta);
1147    //-------------------------
1148    fhTrackCentVsPt = new TH2F("fhTrackCentVsPt","pT,trk vs centrality", 50, 0, 50, fNumberOfCentralityBins,0,100);
1149    fOutputList->Add(fhTrackCentVsPt);
1150    //-------------------------
1151    fhVertexZ = new TH1F("fhVertexZ","z vertex",40,-20,20);
1152    fOutputList->Add(fhVertexZ);
1153    //-------------------------
1154    fhVertexZAccept = new TH1F("fhVertexZAccept","z vertex after cut",40,-20,20);
1155    fOutputList->Add(fhVertexZAccept);
1156    //-------------------------
1157    //fhContribVtx = new TH1F("fhContribVtx","contrib to vtx",200,0,200);
1158    //fOutputList->Add(fhContribVtx);
1159    //-------------------------
1160    //fhContribVtxAccept = new TH1F("fhContribVtxAccept","contrib to vtx after cut",200,0,200);
1161    //fOutputList->Add(fhContribVtxAccept);
1162    //-------------------------
1163    fhDphiTriggerJetMinBias = new TH2F("fhDphiTriggerJetMinBias","Deltaphi trig-jet",50,0,100, 100, -0.5*TMath::Pi(),1.5*TMath::Pi());
1164    fOutputList->Add(fhDphiTriggerJetMinBias);
1165
1166    fhDphiTriggerJetCent20 = (TH2F*) fhDphiTriggerJetMinBias->Clone("fhDphiTriggerJetCent20");
1167    fOutputList->Add(fhDphiTriggerJetCent20);
1168    //-------------------------
1169
1170    fhDphiTriggerJetAccept = new TH1F("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -0.5*TMath::Pi(),1.5*TMath::Pi());
1171    fOutputList->Add(fhDphiTriggerJetAccept);
1172    //-------------------------
1173    fhCentrality = new TH1F("fhCentrality","Centrality",100,0,100);
1174    fOutputList->Add(fhCentrality);
1175    //-------------------------
1176    fhCentralityV0M = new TH1F("hCentralityV0M","hCentralityV0M",100,0,100);
1177    fOutputList->Add(fhCentralityV0M);
1178    //-------------------------
1179    fhCentralityV0A = new TH1F("hCentralityV0A","hCentralityV0A",100,0,100);
1180    fOutputList->Add(fhCentralityV0A);
1181    //-------------------------
1182    fhCentralityV0C = new TH1F("hCentralityV0C","hCentralityV0C",100,0,100);
1183    fOutputList->Add(fhCentralityV0C);
1184    //-------------------------
1185    fhCentralityZNA = new TH1F("hCentralityZNA","hCentralityZNA",100,0,100);
1186    fOutputList->Add(fhCentralityZNA);
1187
1188    fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
1189    fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1190    fOutputList->Add(fh1Xsec);
1191
1192    fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
1193    fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1194    fOutputList->Add(fh1Trials);
1195
1196    fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",300,0,300);
1197    fOutputList->Add(fh1PtHard);
1198
1199    fhImpactParameter = new TH1D("fhImpactParameter","impact parameter distribution from HIJING",50,0,10);
1200    fOutputList->Add(fhImpactParameter);
1201
1202    fhImpactParameterTT = new TH1D("fhImpactParameterTT","b versus TT",50,0,10);
1203    fOutputList->Add(fhImpactParameterTT);
1204    // =========== Switch on Sumw2 for all histos ===========
1205    for(Int_t i=0; i<fOutputList->GetEntries(); i++){
1206       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
1207       if(h1){
1208          h1->Sumw2();
1209          continue;
1210       }
1211       THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
1212       if(hn){
1213          hn->Sumw2();
1214       }
1215    }
1216    TH1::AddDirectory(oldStatus);
1217
1218
1219    PostData(1, fOutputList);
1220 }
1221
1222 //________________________________________________________________________
1223 void AliAnalysisTaskHJetSpectra::UserExec(Option_t *){
1224    //executed in each event
1225
1226    if(!InputEvent()){
1227       AliError("??? Event pointer == 0 ???");
1228       return;
1229    }
1230
1231    //Execute only once:  Get tracks, jets, background from arrays if not already given
1232    if(!fInitialized) ExecOnce();
1233    if(fJetArray && fTrackArray && fBackgroundJetArray){
1234       Calculate(InputEvent());
1235    }
1236    PostData(1, fOutputList);
1237 }
1238 //________________________________________________________________________
1239
1240 Double_t AliAnalysisTaskHJetSpectra::RelativePhi(Double_t mphi,Double_t vphi){
1241    //Get relative azimuthal angle of two particles -pi to pi
1242    if      (vphi < -TMath::Pi()) vphi += TMath::TwoPi();
1243    else if (vphi > TMath::Pi())  vphi -= TMath::TwoPi();
1244
1245    if      (mphi < -TMath::Pi()) mphi += TMath::TwoPi();
1246    else if (mphi > TMath::Pi())  mphi -= TMath::TwoPi();
1247
1248    Double_t dphi = mphi - vphi;
1249    if      (dphi < -TMath::Pi()) dphi += TMath::TwoPi();
1250    else if (dphi > TMath::Pi())  dphi -= TMath::TwoPi();
1251
1252    return dphi;//dphi in [-Pi, Pi]
1253 }
1254
1255 //________________________________________________________________________
1256 Double_t AliAnalysisTaskHJetSpectra::EstimateBgRhoMedian(){
1257    //Estimate background rho by means of integrating track pT outside identified jet cones
1258    Double_t rhoMedian = 0.0;
1259
1260    //phi,eta and R2 of jets to be removed
1261    std::vector<Double_t>  jphi;
1262    std::vector<Double_t>  jeta;
1263    std::vector<Double_t> jRsquared;
1264
1265    if(!fBackgroundJetArray) return 0.0;
1266    if(!fTrackArray)         return 0.0;
1267
1268    for(Int_t i = 0; i < fBackgroundJetArray->GetEntries(); i++){
1269       AliEmcalJet* backgroundJet = static_cast<AliEmcalJet*>(fBackgroundJetArray->At(i));
1270
1271       if(!backgroundJet){
1272          AliError(Form("%s: Could not receive jet %d", GetName(), i));
1273          continue;
1274       }
1275       if(!IsBackgroundJetInAcceptance(backgroundJet)) continue; //apply minimum pT cut on jet to be removed from bg
1276       jphi.push_back(RelativePhi(backgroundJet->Phi(),0.0)); //-pi,pi
1277       jeta.push_back(backgroundJet->Eta());
1278       jRsquared.push_back(1.78*backgroundJet->Area()/TMath::Pi()); //1.78 = JetArea_R04/JetArea_R03
1279    }
1280
1281
1282    static Double_t nOutCone[10][4];
1283    static Double_t sumPtOutOfCone[10][4];
1284    Double_t rndphi, rndeta;
1285    Double_t rndphishift, rndetashift;
1286    Double_t dphi, deta;
1287    Bool_t   bIsInCone;
1288
1289
1290    for(Int_t ie=0; ie < fnEta; ie++){
1291       for(Int_t ip=0; ip < fnPhi; ip++){
1292          nOutCone[ip][ie]       = 0.0;     //initialize counter
1293          sumPtOutOfCone[ip][ie] = 0.0;
1294       }
1295    }
1296
1297    //get area in cells out of identified jet cones
1298    if(jphi.size()==0){ //no jet to be removed from the bg => all areas have their nominal area
1299       for(Int_t ie=0; ie < fnEta; ie++){
1300          for(Int_t ip=0; ip < fnPhi; ip++){
1301             nOutCone[ip][ie] = fNofRndTrials;
1302          }
1303       }
1304    }else{
1305       for(Int_t it=0; it<fNofRndTrials; it++){
1306
1307          rndphi = fRandom->Uniform(0, fPhiSize);
1308          rndeta = fRandom->Uniform(0, fEtaSize);
1309
1310          for(Int_t ip=0; ip<fnPhi; ip++){  //move radom position to each cell
1311             rndphishift = rndphi + ip*fPhiSize - TMath::Pi();
1312             for(Int_t ie=0; ie<fnEta; ie++){
1313                rndetashift = rndeta + ie*fEtaSize - fEtaSize;
1314
1315                bIsInCone = 0; //tag if trial is in the jet cone
1316                for(Int_t ij=0; ij< (Int_t) jRsquared.size(); ij++){
1317                   deta = jeta[ij] - rndetashift;
1318                   dphi = RelativePhi(rndphishift,jphi[ij]);
1319                   if((dphi*dphi + deta*deta) < jRsquared[ij]){
1320                      bIsInCone = 1;
1321                      break;
1322                   }
1323                }
1324                if(!bIsInCone) nOutCone[ip][ie]++;
1325             }
1326          }
1327       }
1328    }
1329
1330    Int_t phicell,etacell;
1331    for(Int_t ip=0; ip < fTrackArray->GetEntries(); ip++){
1332       AliVTrack* part = static_cast<AliVTrack*>(fTrackArray->At(ip));
1333       if(!part) continue;
1334       if(!IsTrackInAcceptance((AliVParticle*) part)) continue;
1335
1336       bIsInCone = 0; //init
1337       for(Int_t ij=0; ij<(Int_t) jRsquared.size(); ij++){
1338          dphi = RelativePhi(jphi[ij], part->Phi());
1339          deta = jeta[ij] - part->Eta();
1340          if((dphi*dphi + deta*deta) < jRsquared[ij]){
1341             bIsInCone = 1;
1342             break;
1343          }
1344       }
1345       if(!bIsInCone){
1346          phicell = TMath::Nint(TMath::Floor((RelativePhi(part->Phi(),0.0) + TMath::Pi())/fPhiSize));
1347          etacell = TMath::Nint(TMath::Floor((part->Eta()+fEtaSize)/fEtaSize));
1348          sumPtOutOfCone[phicell][etacell]+= part->Pt();
1349       }
1350    }
1351    // Calculate rho
1352    static Double_t rhoInCells[20];
1353    Double_t  relativeArea;
1354    Int_t  nCells=0;
1355    Double_t bufferArea=0.0, bufferPt=0.0; //sum cells where A< fJetFreeAreaFrac
1356    for(Int_t ip=0; ip<fnPhi; ip++){
1357       for(Int_t ie=0; ie<fnEta; ie++){
1358          relativeArea = nOutCone[ip][ie]/fNofRndTrials;
1359
1360          bufferArea += relativeArea;
1361          bufferPt   += sumPtOutOfCone[ip][ie];
1362          if(bufferArea > fJetFreeAreaFrac){
1363             rhoInCells[nCells] = bufferPt/(bufferArea*fCellArea);
1364
1365             bufferArea = 0.0;
1366             bufferPt   = 0.0;
1367             nCells++;
1368          }
1369       }
1370    }
1371
1372    if(nCells>0){
1373       rhoMedian = TMath::Median(nCells, rhoInCells);
1374    }
1375    return rhoMedian;
1376 }
1377
1378 //________________________________________________________________________
1379 Double_t AliAnalysisTaskHJetSpectra::EstimateBgCone(){
1380    //Estimate background rho by means of integrating track pT outside identified jet cones
1381    Double_t rhoPerpCone = 0.0;
1382
1383    Double_t pTleading  = -1.0;
1384    Double_t phiLeading = 1000.;
1385    Double_t etaLeading = 1000.;
1386
1387    if(!fJetArray) return 0.0;
1388
1389    for(Int_t ij = 0; ij < fJetArray->GetEntries(); ij++){
1390       AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJetArray->At(ij));
1391       if(!jet){
1392          AliError(Form("%s: Could not receive jet %d", GetName(), ij));
1393          continue;
1394       }
1395       if(!IsSignalJetInAcceptance(jet)) continue;
1396
1397       if(pTleading < jet->Pt()){
1398          pTleading  = jet->Pt();
1399          phiLeading = jet->Phi();
1400          etaLeading = jet->Eta();
1401       }
1402    }
1403    if(pTleading < 0.0) return 0.0;
1404
1405    Double_t phileftcone  = phiLeading + TMath::Pi()/2;
1406    Double_t phirightcone = phiLeading - TMath::Pi()/2;
1407
1408    /* Double_t dp, de;
1409
1410    for(Int_t ip=0; ip < fTrackArray->GetEntries(); ip++){
1411
1412       AliVTrack* part = static_cast<AliVTrack*>(fTrackArray->At(ip));
1413       if(!part) continue;
1414       if(!IsTrackInAcceptance((AliVParticle*) part)) continue;
1415
1416
1417       dp = RelativePhi(phileftcone, part->Phi());
1418       de = etaLeading - part->Eta();
1419       if( dp*dp + de*de < fRandConeRadiusSquared ) rhoPerpCone += part->Pt();
1420
1421       dp = RelativePhi(phirightcone, part->Phi());
1422       if( dp*dp + de*de < fRandConeRadiusSquared) rhoPerpCone += part->Pt();
1423
1424    }*/
1425
1426    rhoPerpCone +=  GetConePt(etaLeading, phileftcone,  fRandConeRadius);
1427    rhoPerpCone +=  GetConePt(etaLeading, phirightcone, fRandConeRadius);
1428
1429    //normalize total pT by two times cone are
1430    rhoPerpCone = rhoPerpCone/(2*TMath::Pi()*fRandConeRadiusSquared);
1431
1432
1433
1434    return rhoPerpCone;
1435 }
1436 //________________________________________________________________________
1437 Bool_t AliAnalysisTaskHJetSpectra::DistantCones(Double_t phi1, Double_t eta1, Double_t r1, Double_t phi2, Double_t eta2, Double_t r2){
1438    //checks if the two cones are farther away than the sum of their radii
1439
1440    Double_t dphi = RelativePhi(phi1,phi2);
1441    Double_t deta = eta1-eta2;
1442    Double_t d = r1+r2;
1443    if( dphi*dphi + deta*deta < d*d ) return kFALSE;
1444
1445    return kTRUE;
1446 }