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