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