PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx

   1 // $Id$
   2 //
   3 // Jet deltaPt task.
   4 //
   5 // Author: S.Aiola
   6
   7 #include <TClonesArray.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13 #include <TRandom3.h>
  14
  15 #include "AliVCluster.h"
  16 #include "AliVParticle.h"
  17 #include "AliVTrack.h"
  18 #include "AliEmcalJet.h"
  19 #include "AliRhoParameter.h"
  20 #include "AliLog.h"
  21 #include "AliJetContainer.h"
  22 #include "AliParticleContainer.h"
  23 #include "AliClusterContainer.h"
  24
  25 #include "AliAnalysisTaskDeltaPt.h"
  26
  27 ClassImp(AliAnalysisTaskDeltaPt)
  28
  29 //________________________________________________________________________
  30 AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() :
  31   AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPt", kTRUE),
  32   fMCJetPtThreshold(1),
  33   fMinRC2LJ(-1),
  34   fRCperEvent(-1),
  35   fConeRadius(0.2),
  36   fConeMinEta(-0.9),
  37   fConeMaxEta(0.9),
  38   fConeMinPhi(0),
  39   fConeMaxPhi(TMath::Pi()*2),
  40   fJetsCont(0),
  41   fTracksCont(0),
  42   fCaloClustersCont(0),
  43   fEmbJetsCont(0),
  44   fEmbTracksCont(0),
  45   fEmbCaloClustersCont(0),
  46   fRandTracksCont(0),
  47   fRandCaloClustersCont(0),
  48   fHistRhovsCent(0),
  49   fHistRCPhiEta(0),
  50   fHistRCPt(0),
  51   fHistRCPtExLJ(0),
  52   fHistRCPtExPartialLJ(0),
  53   fHistRCPtRand(0),
  54   fHistRhoVSRCPt(0),
  55   fHistDeltaPtRCvsEP(0),
  56   fHistDeltaPtRCExLJ(0),
  57   fHistDeltaPtRCExPartialLJ(0),
  58   fHistDeltaPtRCRand(0),
  59   fHistEmbJetsPtArea(0),
  60   fHistEmbJetsCorrPtArea(0),
  61   fHistEmbPartPtvsJetPt(0),
  62   fHistEmbPartPtvsJetCorrPt(0),
  63   fHistJetPtvsJetCorrPt(0),
  64   fHistDistLeadPart2JetAxis(0),
  65   fHistEmbBkgArea(0),
  66   fHistRhoVSEmbBkg(0),
  67   fHistDeltaPtEmbArea(0),
  68   fHistDeltaPtEmbvsEP(0),
  69   fHistRCPtExLJVSDPhiLJ(0),
  70   fHistRCPtExPartialLJVSDPhiLJ(0),
  71   fHistEmbJetsPhiEta(0),
  72   fHistLeadPartPhiEta(0)
  73 {
  74   // Default constructor.
  75
  76   fHistRCPt = 0;
  77   fHistRCPtExLJ = 0;
  78   fHistRCPtExPartialLJ = 0;
  79   fHistRCPtRand = 0;
  80   fHistRhoVSRCPt = 0;
  81   fHistDeltaPtRCvsEP = 0;
  82   fHistDeltaPtRCExLJ = 0;
  83   fHistDeltaPtRCExPartialLJ = 0;
  84   fHistDeltaPtRCRand = 0;
  85   fHistEmbJetsPtArea = 0;
  86   fHistEmbJetsCorrPtArea = 0;
  87   fHistEmbPartPtvsJetPt = 0;
  88   fHistEmbPartPtvsJetCorrPt = 0;
  89   fHistJetPtvsJetCorrPt = 0;
  90   fHistDistLeadPart2JetAxis = 0;
  91   fHistEmbBkgArea = 0;
  92   fHistRhoVSEmbBkg = 0;
  93   fHistDeltaPtEmbArea = 0;
  94   fHistDeltaPtEmbvsEP = 0;
  95
  96   SetMakeGeneralHistograms(kTRUE);
  97 }
  98
  99 //________________________________________________________________________
 100 AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) :
 101   AliAnalysisTaskEmcalJet(name, kTRUE),
 102   fMCJetPtThreshold(1),
 103   fMinRC2LJ(-1),
 104   fRCperEvent(-1),
 105   fConeRadius(0.2),
 106   fConeMinEta(-0.9),
 107   fConeMaxEta(0.9),
 108   fConeMinPhi(0),
 109   fConeMaxPhi(TMath::Pi()*2),
 110   fJetsCont(0),
 111   fTracksCont(0),
 112   fCaloClustersCont(0),
 113   fEmbJetsCont(0),
 114   fEmbTracksCont(0),
 115   fEmbCaloClustersCont(0),
 116   fRandTracksCont(0),
 117   fRandCaloClustersCont(0),
 118   fHistRhovsCent(0),
 119   fHistRCPhiEta(0),
 120   fHistRCPt(0),
 121   fHistRCPtExLJ(0),
 122   fHistRCPtExPartialLJ(0),
 123   fHistRCPtRand(0),
 124   fHistRhoVSRCPt(0),
 125   fHistDeltaPtRCvsEP(0),
 126   fHistDeltaPtRCExLJ(0),
 127   fHistDeltaPtRCExPartialLJ(0),
 128   fHistDeltaPtRCRand(0),
 129   fHistEmbJetsPtArea(0),
 130   fHistEmbJetsCorrPtArea(0),
 131   fHistEmbPartPtvsJetPt(0),
 132   fHistEmbPartPtvsJetCorrPt(0),
 133   fHistJetPtvsJetCorrPt(0),
 134   fHistDistLeadPart2JetAxis(0),
 135   fHistEmbBkgArea(0),
 136   fHistRhoVSEmbBkg(0),
 137   fHistDeltaPtEmbArea(0),
 138   fHistDeltaPtEmbvsEP(0),
 139   fHistRCPtExLJVSDPhiLJ(0),
 140   fHistRCPtExPartialLJVSDPhiLJ(0),
 141   fHistEmbJetsPhiEta(0),
 142   fHistLeadPartPhiEta(0)
 143 {
 144   // Standard constructor.
 145
 146   fHistRCPt = 0;
 147   fHistRCPtExLJ = 0;
 148   fHistRCPtExPartialLJ = 0;
 149   fHistRCPtRand = 0;
 150   fHistRhoVSRCPt = 0;
 151   fHistDeltaPtRCvsEP = 0;
 152   fHistDeltaPtRCExLJ = 0;
 153   fHistDeltaPtRCExPartialLJ = 0;
 154   fHistDeltaPtRCRand = 0;
 155   fHistEmbJetsPtArea = 0;
 156   fHistEmbJetsCorrPtArea = 0;
 157   fHistEmbPartPtvsJetPt = 0;
 158   fHistEmbPartPtvsJetCorrPt = 0;
 159   fHistJetPtvsJetCorrPt = 0;
 160   fHistDistLeadPart2JetAxis = 0;
 161   fHistEmbBkgArea = 0;
 162   fHistRhoVSEmbBkg = 0;
 163   fHistDeltaPtEmbArea = 0;
 164   fHistDeltaPtEmbvsEP = 0;
 165
 166   SetMakeGeneralHistograms(kTRUE);
 167 }
 168
 169 //________________________________________________________________________
 170 void AliAnalysisTaskDeltaPt::AllocateHistogramArrays()
 171 {
 172   fHistRCPt = new TH1*[fNcentBins];
 173   fHistRCPtExLJ = new TH1*[fNcentBins];
 174   fHistRCPtExPartialLJ = new TH1*[fNcentBins];
 175   fHistRCPtRand = new TH1*[fNcentBins];
 176   fHistRhoVSRCPt = new TH2*[fNcentBins];
 177   fHistDeltaPtRCvsEP = new TH2*[fNcentBins];
 178   fHistDeltaPtRCExLJ = new TH1*[fNcentBins];
 179   fHistDeltaPtRCExPartialLJ = new TH1*[fNcentBins];
 180   fHistDeltaPtRCRand = new TH1*[fNcentBins];
 181   fHistEmbJetsPtArea = new TH3*[fNcentBins];
 182   fHistEmbJetsCorrPtArea = new TH3*[fNcentBins];
 183   fHistEmbPartPtvsJetPt = new TH2*[fNcentBins];
 184   fHistEmbPartPtvsJetCorrPt = new TH2*[fNcentBins];
 185   fHistJetPtvsJetCorrPt = new TH2*[fNcentBins];
 186   fHistDistLeadPart2JetAxis = new TH1*[fNcentBins];
 187   fHistEmbBkgArea = new TH2*[fNcentBins];
 188   fHistRhoVSEmbBkg = new TH2*[fNcentBins];
 189   fHistDeltaPtEmbArea = new TH2*[fNcentBins];
 190   fHistDeltaPtEmbvsEP = new TH2*[fNcentBins];
 191
 192   for (Int_t i = 0; i < fNcentBins; i++) {
 193     fHistRCPt[i] = 0;
 194     fHistRCPtExLJ[i] = 0;
 195     fHistRCPtExPartialLJ[i] = 0;
 196     fHistRCPtRand[i] = 0;
 197     fHistRhoVSRCPt[i] = 0;
 198     fHistDeltaPtRCvsEP[i] = 0;
 199     fHistDeltaPtRCExLJ[i] = 0;
 200     fHistDeltaPtRCExPartialLJ[i] = 0;
 201     fHistDeltaPtRCRand[i] = 0;
 202     fHistEmbJetsPtArea[i] = 0;
 203     fHistEmbJetsCorrPtArea[i] = 0;
 204     fHistEmbPartPtvsJetPt[i] = 0;
 205     fHistEmbPartPtvsJetCorrPt[i] = 0;
 206     fHistJetPtvsJetCorrPt[i] = 0;
 207     fHistDistLeadPart2JetAxis[i] = 0;
 208     fHistEmbBkgArea[i] = 0;
 209     fHistRhoVSEmbBkg[i] = 0;
 210     fHistDeltaPtEmbArea[i] = 0;
 211     fHistDeltaPtEmbvsEP[i] = 0;
 212   }
 213 }
 214
 215 //________________________________________________________________________
 216 void AliAnalysisTaskDeltaPt::UserCreateOutputObjects()
 217 {
 218   // Create user output.
 219
 220   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 221
 222   AllocateHistogramArrays();
 223
 224   fHistRhovsCent = new TH2F("fHistRhovsCent", "fHistRhovsCent", 101, -1,  100, fNbins, 0, fMaxBinPt*2);
 225   fHistRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
 226   fHistRhovsCent->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
 227   fOutput->Add(fHistRhovsCent);
 228
 229   fJetsCont = GetJetContainer("Jets");
 230   fTracksCont = GetParticleContainer("Tracks");
 231   fCaloClustersCont = GetClusterContainer("CaloClusters");
 232   fEmbJetsCont = GetJetContainer("EmbJets");
 233   fEmbTracksCont = GetParticleContainer("EmbTracks");
 234   fEmbCaloClustersCont = GetClusterContainer("EmbCaloClusters");
 235   fRandTracksCont = GetParticleContainer("RandTracks");
 236   fRandCaloClustersCont = GetClusterContainer("RandCaloClusters");
 237
 238   if (fTracksCont || fCaloClustersCont) {
 239     fHistRCPhiEta = new TH2F("fHistRCPhiEta","fHistRCPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 240     fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
 241     fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
 242     fOutput->Add(fHistRCPhiEta);
 243
 244     if (fJetsCont) {
 245       fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
 246       fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 247       fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
 248       fOutput->Add(fHistRCPtExLJVSDPhiLJ);
 249
 250       fHistRCPtExPartialLJVSDPhiLJ = new TH2F("fHistRCPtExPartialLJVSDPhiLJ","fHistRCPtExPartialLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
 251       fHistRCPtExPartialLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 252       fHistRCPtExPartialLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
 253       fOutput->Add(fHistRCPtExPartialLJVSDPhiLJ);
 254     }
 255   }
 256
 257   if (fEmbJetsCont) {
 258     fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 259     fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
 260     fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
 261     fOutput->Add(fHistEmbJetsPhiEta);
 262
 263     fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 264     fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
 265     fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
 266     fOutput->Add(fHistLeadPartPhiEta);
 267   }
 268
 269   TString histname;
 270
 271   const Int_t nbinsZ = 12;
 272   Double_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
 273
 274   Double_t *binsPt       = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
 275   Double_t *binsCorrPt   = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
 276   Double_t *binsArea     = GenerateFixedBinArray(50, 0, 2);
 277
 278   for (Int_t i = 0; i < fNcentBins; i++) {
 279     if (fTracksCont || fCaloClustersCont) {
 280       histname = "fHistRCPt_";
 281       histname += i;
 282       fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 283       fHistRCPt[i]->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 284       fHistRCPt[i]->GetYaxis()->SetTitle("counts");
 285       fOutput->Add(fHistRCPt[i]);
 286
 287       histname = "fHistRhoVSRCPt_";
 288       histname += i;
 289       fHistRhoVSRCPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 290       fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
 291       fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 292       fOutput->Add(fHistRhoVSRCPt[i]);
 293
 294       histname = "fHistDeltaPtRCvsEP_";
 295       histname += i;
 296       fHistDeltaPtRCvsEP[i] = new TH2F(histname.Data(), histname.Data(), 101, 0, TMath::Pi()*1.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
 297       fHistDeltaPtRCvsEP[i]->GetXaxis()->SetTitle("#phi_{RC} - #psi_{RP}");
 298       fHistDeltaPtRCvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
 299       fHistDeltaPtRCvsEP[i]->GetZaxis()->SetTitle("counts");
 300       fOutput->Add(fHistDeltaPtRCvsEP[i]);
 301
 302       if (fJetsCont) {
 303         histname = "fHistRCPtExLJ_";
 304         histname += i;
 305         fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 306         fHistRCPtExLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
 307         fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
 308         fOutput->Add(fHistRCPtExLJ[i]);
 309
 310         histname = "fHistDeltaPtRCExLJ_";
 311         histname += i;
 312         fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 313         fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
 314         fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
 315         fOutput->Add(fHistDeltaPtRCExLJ[i]);
 316
 317         histname = "fHistRCPtExPartialLJ_";
 318         histname += i;
 319         fHistRCPtExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 320         fHistRCPtExPartialLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
 321         fHistRCPtExPartialLJ[i]->GetYaxis()->SetTitle("counts");
 322         fOutput->Add(fHistRCPtExPartialLJ[i]);
 323
 324         histname = "fHistDeltaPtRCExPartialLJ_";
 325         histname += i;
 326         fHistDeltaPtRCExPartialLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 327         fHistDeltaPtRCExPartialLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
 328         fHistDeltaPtRCExPartialLJ[i]->GetYaxis()->SetTitle("counts");
 329         fOutput->Add(fHistDeltaPtRCExPartialLJ[i]);
 330       }
 331     }
 332
 333     if (fRandTracksCont || fRandCaloClustersCont) {
 334       histname = "fHistRCPtRand_";
 335       histname += i;
 336       fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 337       fHistRCPtRand[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
 338       fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
 339       fOutput->Add(fHistRCPtRand[i]);
 340
 341       histname = "fHistDeltaPtRCRand_";
 342       histname += i;
 343       fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 344       fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
 345       fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
 346       fOutput->Add(fHistDeltaPtRCRand[i]);
 347     }
 348
 349     if (fEmbJetsCont) {
 350       histname = "fHistEmbJetsPtArea_";
 351       histname += i;
 352       fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins, binsPt, nbinsZ, binsZ);
 353       fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
 354       fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
 355       fOutput->Add(fHistEmbJetsPtArea[i]);
 356
 357       histname = "fHistEmbJetsCorrPtArea_";
 358       histname += i;
 359       fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 50, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
 360       fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
 361       fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
 362       fOutput->Add(fHistEmbJetsCorrPtArea[i]);
 363
 364       histname = "fHistEmbPartPtvsJetPt_";
 365       histname += i;
 366       fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 367       fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 368       fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
 369       fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
 370       fOutput->Add(fHistEmbPartPtvsJetPt[i]);
 371
 372       histname = "fHistEmbPartPtvsJetCorrPt_";
 373       histname += i;
 374       fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
 375                                               fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
 376       fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 377       fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
 378       fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
 379       fOutput->Add(fHistEmbPartPtvsJetCorrPt[i]);
 380
 381       histname = "fHistJetPtvsJetCorrPt_";
 382       histname += i;
 383       fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(),
 384                                           fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
 385       fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
 386       fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
 387       fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
 388       fOutput->Add(fHistJetPtvsJetCorrPt[i]);
 389
 390       histname = "fHistDistLeadPart2JetAxis_";
 391       histname += i;
 392       fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
 393       fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
 394       fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
 395       fOutput->Add(fHistDistLeadPart2JetAxis[i]);
 396
 397       histname = "fHistEmbBkgArea_";
 398       histname += i;
 399       fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 50, 0, 2, fNbins, fMinBinPt, fMaxBinPt);
 400       fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
 401       fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 402       fOutput->Add(fHistEmbBkgArea[i]);
 403
 404       histname = "fHistRhoVSEmbBkg_";
 405       histname += i;
 406       fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 407       fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
 408       fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 409       fOutput->Add(fHistRhoVSEmbBkg[i]);
 410
 411       histname = "fHistDeltaPtEmbArea_";
 412       histname += i;
 413       fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(),
 414                                             50, 0, 2, fNbins * 2, -fMaxBinPt, fMaxBinPt);
 415       fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
 416       fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
 417       fHistDeltaPtEmbArea[i]->GetZaxis()->SetTitle("counts");
 418       fOutput->Add(fHistDeltaPtEmbArea[i]);
 419
 420       histname = "fHistDeltaPtEmbvsEP_";
 421       histname += i;
 422       fHistDeltaPtEmbvsEP[i] = new TH2F(histname.Data(), histname.Data(), 101, 0, TMath::Pi()*1.01, fNbins * 2, -fMaxBinPt, fMaxBinPt);
 423       fHistDeltaPtEmbvsEP[i]->GetXaxis()->SetTitle("#phi_{jet} - #Psi_{EP}");
 424       fHistDeltaPtEmbvsEP[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
 425       fHistDeltaPtEmbvsEP[i]->GetZaxis()->SetTitle("counts");
 426       fOutput->Add(fHistDeltaPtEmbvsEP[i]);
 427     }
 428   }
 429
 430   delete[] binsPt;
 431   delete[] binsCorrPt;
 432   delete[] binsArea;
 433
 434   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 435 }
 436
 437 //________________________________________________________________________
 438 Bool_t AliAnalysisTaskDeltaPt::FillHistograms()
 439 {
 440   // Fill histograms.
 441
 442   fHistRhovsCent->Fill(fCent, fRhoVal);
 443
 444   // ************
 445   // Random cones
 446   // _________________________________
 447
 448   const Float_t rcArea = fConeRadius * fConeRadius * TMath::Pi();
 449   Float_t RCpt = 0;
 450   Float_t RCeta = 0;
 451   Float_t RCphi = 0;
 452
 453   if (fTracksCont || fCaloClustersCont) {
 454
 455     for (Int_t i = 0; i < fRCperEvent; i++) {
 456       // Simple random cones
 457       RCpt = 0;
 458       RCeta = 0;
 459       RCphi = 0;
 460       GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0);
 461       if (RCpt > 0) {
 462         fHistRCPhiEta->Fill(RCeta, RCphi);
 463         fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt);
 464
 465         fHistRCPt[fCentBin]->Fill(RCpt);
 466
 467         Double_t ep = RCphi - fEPV0;
 468         while (ep < 0) ep += TMath::Pi();
 469         while (ep >= TMath::Pi()) ep -= TMath::Pi();
 470
 471         fHistDeltaPtRCvsEP[fCentBin]->Fill(ep, RCpt - rcArea * fRhoVal);
 472       }
 473
 474       if (fJetsCont) {
 475
 476         // Random cones far from leading jet
 477         AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
 478
 479         RCpt = 0;
 480         RCeta = 0;
 481         RCphi = 0;
 482         GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
 483         if (RCpt > 0) {
 484           if (jet) {
 485             Float_t dphi = RCphi - jet->Phi();
 486             if (dphi > 4.8) dphi -= TMath::Pi() * 2;
 487             if (dphi < -1.6) dphi += TMath::Pi() * 2;
 488             fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
 489           }
 490           fHistRCPtExLJ[fCentBin]->Fill(RCpt);
 491           fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 492         }
 493
 494         //partial exclusion
 495         if(fBeamType == kpA) {
 496
 497           RCpt = 0;
 498           RCeta = 0;
 499           RCphi = 0;
 500           GetRandomCone(RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet, kTRUE);
 501
 502           if (RCpt > 0) {
 503             if (jet) {
 504               Float_t dphi = RCphi - jet->Phi();
 505               if (dphi > 4.8) dphi -= TMath::Pi() * 2;
 506               if (dphi < -1.6) dphi += TMath::Pi() * 2;
 507               fHistRCPtExPartialLJVSDPhiLJ->Fill(RCpt, dphi);
 508             }
 509             fHistRCPtExPartialLJ[fCentBin]->Fill(RCpt);
 510             fHistDeltaPtRCExPartialLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 511           }
 512         }
 513       }
 514     }
 515   }
 516
 517   // Random cones with randomized particles
 518   if (fRandTracksCont || fRandCaloClustersCont) {
 519     RCpt = 0;
 520     RCeta = 0;
 521     RCphi = 0;
 522     GetRandomCone(RCpt, RCeta, RCphi, fRandTracksCont, fRandCaloClustersCont, 0);
 523     if (RCpt > 0) {
 524       fHistRCPtRand[fCentBin]->Fill(RCpt);
 525       fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 526     }
 527   }
 528
 529   // ************
 530   // Embedding
 531   // _________________________________
 532
 533   if (fEmbJetsCont) {
 534
 535     AliEmcalJet *embJet = NextEmbeddedJet(kTRUE);
 536
 537     while (embJet != 0) {
 538       TLorentzVector mom;
 539       fEmbJetsCont->GetLeadingHadronMomentum(mom,embJet);
 540
 541       Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - mom.Eta()) * (embJet->Eta() - mom.Eta()) + (embJet->Phi() - mom.Phi()) * (embJet->Phi() - mom.Phi()));
 542
 543       fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt());
 544       fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fRhoVal);
 545       fHistLeadPartPhiEta->Fill(mom.Eta(), mom.Phi());
 546       fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
 547
 548       fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), mom.Pt());
 549       fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area(), mom.Pt());
 550       fHistEmbJetsPhiEta->Fill(embJet->Eta(), embJet->Phi());
 551       fHistJetPtvsJetCorrPt[fCentBin]->Fill(embJet->Pt(), embJet->Pt() - fRhoVal * embJet->Area());
 552
 553       fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->MCPt());
 554       fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * embJet->Area(), embJet->Pt() - embJet->MCPt());
 555       fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
 556
 557       Double_t ep = embJet->Phi() - fEPV0;
 558       while (ep < 0) ep += TMath::Pi();
 559       while (ep >= TMath::Pi()) ep -= TMath::Pi();
 560
 561       fHistDeltaPtEmbvsEP[fCentBin]->Fill(ep, embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
 562
 563       embJet = NextEmbeddedJet();
 564     }
 565   }
 566
 567   return kTRUE;
 568 }
 569
 570 //________________________________________________________________________
 571 AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Bool_t reset)
 572 {
 573   // Get the next accepted embedded jet.
 574
 575   Int_t i = reset ? 0 : -1;
 576
 577   AliEmcalJet* jet = fEmbJetsCont->GetNextAcceptJet(i);
 578   while (jet && jet->MCPt() < fMCJetPtThreshold) jet = fEmbJetsCont->GetNextAcceptJet();
 579
 580   return jet;
 581 }
 582
 583 //________________________________________________________________________
 584 void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
 585                                            AliParticleContainer* tracks, AliClusterContainer* clusters,
 586                                            AliEmcalJet *jet, Bool_t bPartialExclusion) const
 587 {
 588   // Get rigid cone.
 589
 590   eta = -999;
 591   phi = -999;
 592   pt = 0;
 593
 594   if (!tracks && !clusters)
 595     return;
 596
 597   Float_t LJeta = 999;
 598   Float_t LJphi = 999;
 599
 600   if (jet) {
 601     LJeta = jet->Eta();
 602     LJphi = jet->Phi();
 603   }
 604
 605   Float_t maxEta = fConeMaxEta;
 606   Float_t minEta = fConeMinEta;
 607   Float_t maxPhi = fConeMaxPhi;
 608   Float_t minPhi = fConeMinPhi;
 609
 610   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 611   if (minPhi < 0) minPhi = 0;
 612
 613   Float_t dLJ = 0;
 614   Int_t repeats = 0;
 615   Bool_t reject = kTRUE;
 616   do {
 617     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
 618     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
 619     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 620
 621     if(bPartialExclusion) {
 622       reject = kFALSE;
 623
 624       TRandom3 rnd;
 625       rnd.SetSeed(0);
 626
 627       Double_t ncoll = GetNColl();
 628
 629       Double_t prob = 0.;
 630       if(ncoll>0)
 631         prob = 1./ncoll;
 632
 633       if(rnd.Rndm()<=prob) reject = kTRUE; //reject cone
 634     }
 635
 636     repeats++;
 637   } while (dLJ < fMinRC2LJ && repeats < 999 && reject);
 638
 639   if (repeats == 999) {
 640     AliWarning(Form("%s: Could not get random cone!", GetName()));
 641     return;
 642   }
 643
 644   if (clusters) {
 645     AliVCluster* cluster = clusters->GetNextAcceptCluster(0);
 646     while (cluster) {
 647       TLorentzVector nPart;
 648       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 649
 650       Float_t cluseta = nPart.Eta();
 651       Float_t clusphi = nPart.Phi();
 652
 653       if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi()))
 654         clusphi += 2 * TMath::Pi();
 655       if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi()))
 656         clusphi -= 2 * TMath::Pi();
 657
 658       Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi));
 659       if (d <= fConeRadius)
 660         pt += nPart.Pt();
 661
 662       cluster = clusters->GetNextAcceptCluster();
 663     }
 664   }
 665
 666   if (tracks) {
 667     AliVParticle* track = tracks->GetNextAcceptParticle(0);
 668     while(track) {
 669       Float_t tracketa = track->Eta();
 670       Float_t trackphi = track->Phi();
 671
 672       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
 673         trackphi += 2 * TMath::Pi();
 674       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
 675         trackphi -= 2 * TMath::Pi();
 676
 677       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
 678       if (d <= fConeRadius)
 679         pt += track->Pt();
 680
 681       track = tracks->GetNextAcceptParticle();
 682     }
 683   }
 684 }
 685
 686 //________________________________________________________________________
 687 void AliAnalysisTaskDeltaPt::SetConeEtaPhiEMCAL()
 688 {
 689   // Set default cuts for full cones
 690
 691   SetConeEtaLimits(-0.7+fConeRadius,0.7-fConeRadius);
 692   SetConePhiLimits(1.4+fConeRadius,TMath::Pi()-fConeRadius);
 693 }
 694
 695 //________________________________________________________________________
 696 void AliAnalysisTaskDeltaPt::SetConeEtaPhiTPC()
 697 {
 698   // Set default cuts for charged cones
 699
 700   SetConeEtaLimits(-0.9+fConeRadius, 0.9-fConeRadius);
 701   SetConePhiLimits(-10, 10);
 702 }
 703
 704 //________________________________________________________________________
 705 void AliAnalysisTaskDeltaPt::ExecOnce()
 706 {
 707   // Initialize the analysis.
 708
 709   AliAnalysisTaskEmcalJet::ExecOnce();
 710
 711   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 712   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 713   if (fEmbTracksCont && fEmbTracksCont->GetArray() == 0) fEmbTracksCont = 0;
 714   if (fEmbCaloClustersCont && fEmbCaloClustersCont->GetArray() == 0) fEmbCaloClustersCont = 0;
 715   if (fRandTracksCont && fRandTracksCont->GetArray() == 0) fRandTracksCont = 0;
 716   if (fRandCaloClustersCont && fRandCaloClustersCont->GetArray() == 0) fRandCaloClustersCont = 0;
 717   if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
 718   if (fEmbJetsCont && fEmbJetsCont->GetArray() == 0) fEmbJetsCont = 0;
 719
 720   if (fRCperEvent < 0) {
 721     Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi);
 722     Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
 723     fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
 724     if (fRCperEvent == 0)
 725       fRCperEvent = 1;
 726   }
 727
 728   if (fMinRC2LJ < 0)
 729     fMinRC2LJ = fConeRadius * 1.5;
 730
 731   const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
 732   if (fMinRC2LJ > maxDist) {
 733     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
 734                     "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
 735     fMinRC2LJ = maxDist;
 736   }
 737 }
 738
 739 //________________________________________________________________________
 740 Double_t AliAnalysisTaskDeltaPt::GetNColl() const {
 741   // Get NColl - returns value of corresponding bin
 742   // only works for pA
 743   // values taken from V0A slicing https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Tables_with_centrality_bins_from
 744
 745   if(fBeamType == kpA) {
 746
 747     const Int_t nNCollBins = 7;
 748
 749     Double_t centMin[nNCollBins] = {0.,5.,10.,20.,40.,60.,80.};
 750     Double_t centMax[nNCollBins] = {5.,10.,20.,40.,60.,80.,100.};
 751
 752     Double_t nColl[nNCollBins] = {14.7,13.,11.7,9.38,6.49,3.96,1.52};
 753
 754     for(Int_t i = 0; i<nNCollBins; i++) {
 755       if(fCent>=centMin[i] && fCent<centMax[i])
 756         return nColl[i];
 757     }
 758
 759     return -1.;
 760   }
 761   else {
 762     AliWarning(Form("%s: Only works for pA analysis. Returning -1",GetName()));
 763     return -1.;
 764   }
 765 }