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
  22 #include "AliAnalysisTaskDeltaPt.h"
  23
  24 ClassImp(AliAnalysisTaskDeltaPt)
  25
  26 //________________________________________________________________________
  27 AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() :
  28   AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPt", kTRUE),
  29   fMCJetPtThreshold(1),
  30   fMinRC2LJ(-1),
  31   fEmbJetsName(""),
  32   fEmbTracksName(""),
  33   fEmbCaloName(""),
  34   fRandTracksName("TracksRandomized"),
  35   fRandCaloName("CaloClustersRandomized"),
  36   fRCperEvent(-1),
  37   fEmbJets(0),
  38   fEmbTracks(0),
  39   fEmbCaloClusters(0),
  40   fRandTracks(0),
  41   fRandCaloClusters(0),
  42   fEmbeddedClusterNIds(0),
  43   fEmbeddedTrackNIds(0),
  44   fHistRCPhiEta(0),
  45   fHistRCPtExLJVSDPhiLJ(0),
  46   fHistEmbJetsPhiEta(0),
  47   fHistLeadPartPhiEta(0)
  48 {
  49   // Default constructor.
  50
  51   for (Int_t i = 0; i < 4; i++) {
  52     fHistRCPt[i] = 0;
  53     fHistRCPtExLJ[i] = 0;
  54     fHistRCPtRand[i] = 0;
  55     fHistRhoVSRCPt[i] = 0;
  56     fHistDeltaPtRC[i] = 0;
  57     fHistDeltaPtRCExLJ[i] = 0;
  58     fHistDeltaPtRCRand[i] = 0;
  59     fHistEmbRejectedJetsPhiEta[i] = 0;
  60     fHistEmbRejectedJetsPtArea[i] = 0;
  61     fHistEmbNotFoundPt[i] = 0;
  62     fHistEmbNotFoundPhiEta[i] = 0;
  63     fHistEmbJetsPtArea[i] = 0;
  64     fHistEmbJetsCorrPtArea[i] = 0;
  65     fHistEmbPartPtvsJetPt[i] = 0;
  66     fHistEmbPartPtvsJetCorrPt[i] = 0;
  67     fHistJetPtvsJetCorrPt[i] = 0;
  68     fHistDistLeadPart2JetAxis[i] = 0;
  69     fHistEmbBkgArea[i] = 0;
  70     fHistRhoVSEmbBkg[i] = 0;
  71     fHistDeltaPtEmbArea[i] = 0;
  72   }
  73
  74   SetMakeGeneralHistograms(kTRUE);
  75 }
  76
  77 //________________________________________________________________________
  78 AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) :
  79   AliAnalysisTaskEmcalJet(name, kTRUE),
  80   fMCJetPtThreshold(1),
  81   fMinRC2LJ(-1),
  82   fEmbJetsName(""),
  83   fEmbTracksName(""),
  84   fEmbCaloName(""),
  85   fRandTracksName("TracksRandomized"),
  86   fRandCaloName("CaloClustersRandomized"),
  87   fRCperEvent(-1),
  88   fEmbJets(0),
  89   fEmbTracks(0),
  90   fEmbCaloClusters(0),
  91   fRandTracks(0),
  92   fRandCaloClusters(0),
  93   fEmbeddedClusterNIds(0),
  94   fEmbeddedTrackNIds(0),
  95   fHistRCPhiEta(0),
  96   fHistRCPtExLJVSDPhiLJ(0),
  97   fHistEmbJetsPhiEta(0),
  98   fHistLeadPartPhiEta(0)
  99 {
 100   // Standard constructor.
 101
 102   for (Int_t i = 0; i < 4; i++) {
 103     fHistRCPt[i] = 0;
 104     fHistRCPtExLJ[i] = 0;
 105     fHistRCPtRand[i] = 0;
 106     fHistRhoVSRCPt[i] = 0;
 107     fHistDeltaPtRC[i] = 0;
 108     fHistDeltaPtRCExLJ[i] = 0;
 109     fHistDeltaPtRCRand[i] = 0;
 110     fHistEmbRejectedJetsPhiEta[i] = 0;
 111     fHistEmbRejectedJetsPtArea[i] = 0;
 112     fHistEmbNotFoundPt[i] = 0;
 113     fHistEmbNotFoundPhiEta[i] = 0;
 114     fHistEmbJetsPtArea[i] = 0;
 115     fHistEmbJetsCorrPtArea[i] = 0;
 116     fHistEmbPartPtvsJetPt[i] = 0;
 117     fHistEmbPartPtvsJetCorrPt[i] = 0;
 118     fHistJetPtvsJetCorrPt[i] = 0;
 119     fHistDistLeadPart2JetAxis[i] = 0;
 120     fHistEmbBkgArea[i] = 0;
 121     fHistRhoVSEmbBkg[i] = 0;
 122     fHistDeltaPtEmbArea[i] = 0;
 123   }
 124
 125   SetMakeGeneralHistograms(kTRUE);
 126 }
 127
 128 //________________________________________________________________________
 129 AliAnalysisTaskDeltaPt::~AliAnalysisTaskDeltaPt()
 130 {
 131   // Destructor.
 132 }
 133
 134 //________________________________________________________________________
 135 void AliAnalysisTaskDeltaPt::UserCreateOutputObjects()
 136 {
 137   // Create user output.
 138
 139   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 140
 141   if (!fTracksName.IsNull() || !fCaloName.IsNull()) {
 142     fHistRCPhiEta = new TH2F("fHistRCPhiEta","fHistRCPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 143     fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
 144     fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
 145     fOutput->Add(fHistRCPhiEta);
 146
 147     if (!fJetsName.IsNull()) {
 148       fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
 149       fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 150       fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
 151       fOutput->Add(fHistRCPtExLJVSDPhiLJ);
 152     }
 153   }
 154
 155   if (!fEmbJetsName.IsNull()) {
 156     fHistEmbJetsPhiEta = new TH2F("fHistEmbJetsPhiEta","fHistEmbJetsPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 157     fHistEmbJetsPhiEta->GetXaxis()->SetTitle("#eta");
 158     fHistEmbJetsPhiEta->GetYaxis()->SetTitle("#phi");
 159     fOutput->Add(fHistEmbJetsPhiEta);
 160
 161     fHistLeadPartPhiEta = new TH2F("fHistLeadPartPhiEta","fHistLeadPartPhiEta", 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 162     fHistLeadPartPhiEta->GetXaxis()->SetTitle("#eta");
 163     fHistLeadPartPhiEta->GetYaxis()->SetTitle("#phi");
 164     fOutput->Add(fHistLeadPartPhiEta);
 165   }
 166
 167   TString histname;
 168
 169   const Int_t nbinsZ = 12;
 170   Float_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
 171
 172   Float_t *binsPt       = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
 173   Float_t *binsCorrPt   = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
 174   Float_t *binsArea     = GenerateFixedBinArray(40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
 175
 176   for (Int_t i = 0; i < fNcentBins; i++) {
 177     if (!fTracksName.IsNull() || !fCaloName.IsNull()) {
 178       histname = "fHistRCPt_";
 179       histname += i;
 180       fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 181       fHistRCPt[i]->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 182       fHistRCPt[i]->GetYaxis()->SetTitle("counts");
 183       fOutput->Add(fHistRCPt[i]);
 184
 185       histname = "fHistRhoVSRCPt_";
 186       histname += i;
 187       fHistRhoVSRCPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 188       fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
 189       fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
 190       fOutput->Add(fHistRhoVSRCPt[i]);
 191
 192       histname = "fHistDeltaPtRC_";
 193       histname += i;
 194       fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 195       fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
 196       fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
 197       fOutput->Add(fHistDeltaPtRC[i]);
 198
 199       if (!fJetsName.IsNull()) {
 200         histname = "fHistRCPtExLJ_";
 201         histname += i;
 202         fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 203         fHistRCPtExLJ[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
 204         fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
 205         fOutput->Add(fHistRCPtExLJ[i]);
 206
 207         histname = "fHistDeltaPtRCExLJ_";
 208         histname += i;
 209         fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 210         fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
 211         fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
 212         fOutput->Add(fHistDeltaPtRCExLJ[i]);
 213       }
 214     }
 215
 216     if (!fRandTracksName.IsNull() || !fRandCaloName.IsNull()) {
 217       histname = "fHistRCPtRand_";
 218       histname += i;
 219       fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 220       fHistRCPtRand[i]->GetXaxis()->SetTitle("#it{p}_{T}^{RC} (GeV/#it{c})");
 221       fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
 222       fOutput->Add(fHistRCPtRand[i]);
 223
 224       histname = "fHistDeltaPtRCRand_";
 225       histname += i;
 226       fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 227       fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#delta#it{p}_{T}^{RC} (GeV/#it{c})");
 228       fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
 229       fOutput->Add(fHistDeltaPtRCRand[i]);
 230     }
 231
 232     if (!fEmbJetsName.IsNull()) {
 233       histname = "fHistEmbJetsPtArea_";
 234       histname += i;
 235       fHistEmbJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(), 40, binsArea, fNbins, binsPt, nbinsZ, binsZ);
 236       fHistEmbJetsPtArea[i]->GetXaxis()->SetTitle("area");
 237       fHistEmbJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
 238       fOutput->Add(fHistEmbJetsPtArea[i]);
 239
 240       histname = "fHistEmbJetsCorrPtArea_";
 241       histname += i;
 242       fHistEmbJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(), 40, binsArea, fNbins * 2, binsCorrPt, nbinsZ, binsZ);
 243       fHistEmbJetsCorrPtArea[i]->GetXaxis()->SetTitle("area");
 244       fHistEmbJetsCorrPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,corr} (GeV/#it{c})");
 245       fOutput->Add(fHistEmbJetsCorrPtArea[i]);
 246
 247       histname = "fHistEmbPartPtvsJetPt_";
 248       histname += i;
 249       fHistEmbPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 250       fHistEmbPartPtvsJetPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 251       fHistEmbPartPtvsJetPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
 252       fHistEmbPartPtvsJetPt[i]->GetZaxis()->SetTitle("counts");
 253       fOutput->Add(fHistEmbPartPtvsJetPt[i]);
 254
 255       histname = "fHistEmbPartPtvsJetCorrPt_";
 256       histname += i;
 257       fHistEmbPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
 258       fHistEmbPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 259       fHistEmbPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
 260       fHistEmbPartPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
 261       fOutput->Add(fHistEmbPartPtvsJetCorrPt[i]);
 262
 263       histname = "fHistJetPtvsJetCorrPt_";
 264       histname += i;
 265       fHistJetPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins*2, -fMaxBinPt, fMaxBinPt);
 266       fHistJetPtvsJetCorrPt[i]->GetXaxis()->SetTitle("#it{p}_{T,jet}^{emb} (GeV/#it{c})");
 267       fHistJetPtvsJetCorrPt[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - A#rho (GeV/#it{c})");
 268       fHistJetPtvsJetCorrPt[i]->GetZaxis()->SetTitle("counts");
 269       fOutput->Add(fHistJetPtvsJetCorrPt[i]);
 270
 271       histname = "fHistDistLeadPart2JetAxis_";
 272       histname += i;
 273       fHistDistLeadPart2JetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
 274       fHistDistLeadPart2JetAxis[i]->GetXaxis()->SetTitle("distance");
 275       fHistDistLeadPart2JetAxis[i]->GetYaxis()->SetTitle("counts");
 276       fOutput->Add(fHistDistLeadPart2JetAxis[i]);
 277
 278       histname = "fHistEmbNotFoundPhiEta_";
 279       histname += i;
 280       fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 281       fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
 282       fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
 283       fOutput->Add(fHistEmbNotFoundPhiEta[i]);
 284
 285       histname = "fHistEmbNotFoundPt_";
 286       histname += i;
 287       fHistEmbNotFoundPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 288       fHistEmbNotFoundPt[i]->GetXaxis()->SetTitle("#it{p}_{T,const}^{emb} (GeV/#it{c})");
 289       fHistEmbNotFoundPt[i]->GetYaxis()->SetTitle("counts");
 290       fOutput->Add(fHistEmbNotFoundPt[i]);
 291
 292       histname = "fHistEmbRejectedJetsPhiEta_";
 293       histname += i;
 294       fHistEmbRejectedJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 100, -1, 1, 201, 0, TMath::Pi() * 2.01);
 295       fHistEmbRejectedJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
 296       fHistEmbRejectedJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
 297       fOutput->Add(fHistEmbRejectedJetsPhiEta[i]);
 298
 299       histname = "fHistEmbRejectedJetsPtArea_";
 300       histname += i;
 301       fHistEmbRejectedJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 302       fHistEmbRejectedJetsPtArea[i]->GetXaxis()->SetTitle("area");
 303       fHistEmbRejectedJetsPtArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb,raw} (GeV/#it{c})");
 304       fOutput->Add(fHistEmbRejectedJetsPtArea[i]);
 305
 306       histname = "fHistEmbBkgArea_";
 307       histname += i;
 308       fHistEmbBkgArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt);
 309       fHistEmbBkgArea[i]->GetXaxis()->SetTitle("area");
 310       fHistEmbBkgArea[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 311       fOutput->Add(fHistEmbBkgArea[i]);
 312
 313       histname = "fHistRhoVSEmbBkg_";
 314       histname += i;
 315       fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 316       fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho (GeV/#it{c})");
 317       fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("#it{p}_{T,jet}^{emb} - #sum#it{p}_{T,const}^{emb} (GeV/#it{c})");
 318       fOutput->Add(fHistRhoVSEmbBkg[i]);
 319
 320       histname = "fHistDeltaPtEmbArea_";
 321       histname += i;
 322       fHistDeltaPtEmbArea[i] = new TH2F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins * 2, -fMaxBinPt, fMaxBinPt);
 323       fHistDeltaPtEmbArea[i]->GetXaxis()->SetTitle("area");
 324       fHistDeltaPtEmbArea[i]->GetYaxis()->SetTitle("#delta#it{p}_{T}^{emb} (GeV/#it{c})");
 325       fOutput->Add(fHistDeltaPtEmbArea[i]);
 326     }
 327   }
 328
 329   delete[] binsPt;
 330   delete[] binsCorrPt;
 331   delete[] binsArea;
 332
 333   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 334 }
 335
 336 //________________________________________________________________________
 337 Bool_t AliAnalysisTaskDeltaPt::FillHistograms()
 338 {
 339   // Fill histograms.
 340
 341   // ************
 342   // Random cones
 343   // _________________________________
 344
 345   const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
 346   Float_t RCpt = 0;
 347   Float_t RCeta = 0;
 348   Float_t RCphi = 0;
 349
 350   if (fTracks || fCaloClusters) {
 351
 352     for (Int_t i = 0; i < fRCperEvent; i++) {
 353       // Simple random cones
 354       RCpt = 0;
 355       RCeta = 0;
 356       RCphi = 0;
 357       GetRandomCone(RCpt, RCeta, RCphi, 0);
 358       if (RCpt > 0) {
 359         fHistRCPhiEta->Fill(RCeta, RCphi);
 360         fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt);
 361
 362         fHistRCPt[fCentBin]->Fill(RCpt);
 363         fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 364       }
 365
 366       if (fJets) {
 367
 368         // Random cones far from leading jet
 369         static Int_t sortedJets[9999] = {-1};
 370         GetSortedArray(sortedJets, fJets);
 371
 372         AliEmcalJet* jet = 0;
 373
 374         if (sortedJets[0] >= 0)
 375           jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
 376
 377         RCpt = 0;
 378         RCeta = 0;
 379         RCphi = 0;
 380         GetRandomCone(RCpt, RCeta, RCphi, jet);
 381         if (RCpt > 0) {
 382           if (jet) {
 383             Float_t dphi = RCphi - jet->Phi();
 384             if (dphi > 4.8) dphi -= TMath::Pi() * 2;
 385             if (dphi < -1.6) dphi += TMath::Pi() * 2;
 386             fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
 387           }
 388           fHistRCPtExLJ[fCentBin]->Fill(RCpt);
 389           fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 390         }
 391       }
 392     }
 393   }
 394
 395   // Random cones with randomized particles
 396   if (fRandTracks || fRandCaloClusters) {
 397     RCpt = 0;
 398     RCeta = 0;
 399     RCphi = 0;
 400     GetRandomCone(RCpt, RCeta, RCphi, 0, fRandTracks, fRandCaloClusters);
 401     if (RCpt > 0) {
 402       fHistRCPtRand[fCentBin]->Fill(RCpt);
 403       fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 404     }
 405   }
 406
 407   // ************
 408   // Embedding
 409   // _________________________________
 410
 411   if (fEmbJets) {
 412
 413     AliEmcalJet *embJet = NextEmbeddedJet(0);
 414
 415     Int_t countEmbJets = 0;
 416
 417     while (embJet != 0) {
 418       AliDebug(2,Form("Elaborating embedded jet n. %d", countEmbJets));
 419       countEmbJets++;
 420
 421       if (!AcceptJet(embJet)) {
 422         AliDebug(2,"Embedded jet not accepted, skipping...");
 423         fHistEmbRejectedJetsPhiEta[fCentBin]->Fill(embJet->Eta(), embJet->Phi());
 424         fHistEmbRejectedJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt());
 425
 426         embJet = NextEmbeddedJet();
 427         continue;
 428       }
 429
 430       Double_t maxClusterPt = 0;
 431       Double_t maxClusterEta = 0;
 432       Double_t maxClusterPhi = 0;
 433
 434       Double_t maxTrackPt = 0;
 435       Double_t maxTrackEta = 0;
 436       Double_t maxTrackPhi = 0;
 437
 438       Double_t maxPartPt = 0;
 439       Double_t maxPartEta = 0;
 440       Double_t maxPartPhi = 0;
 441
 442       if (fLeadingHadronType == 1 || fLeadingHadronType == 2) {
 443         AliVCluster *cluster = embJet->GetLeadingCluster(fEmbCaloClusters);
 444         if (cluster) {
 445           TLorentzVector nPart;
 446           cluster->GetMomentum(nPart, fVertex);
 447
 448           maxClusterEta = nPart.Eta();
 449           maxClusterPhi = nPart.Phi();
 450           maxClusterPt = nPart.Pt();
 451         }
 452       }
 453
 454       if (fLeadingHadronType == 0 || fLeadingHadronType == 2) {
 455         AliVParticle *track = embJet->GetLeadingTrack(fEmbTracks);
 456         if (track) {
 457           maxTrackEta = track->Eta();
 458           maxTrackPhi = track->Phi();
 459           maxTrackPt = track->Pt();
 460         }
 461       }
 462
 463       if (maxTrackPt > maxClusterPt) {
 464         maxPartPt = maxTrackPt;
 465         maxPartEta = maxTrackEta;
 466         maxPartPhi = maxTrackPhi;
 467       }
 468       else {
 469         maxPartPt = maxClusterPt;
 470         maxPartEta = maxClusterEta;
 471         maxPartPhi = maxClusterPhi;
 472       }
 473
 474       Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - maxPartEta) * (embJet->Eta() - maxPartEta) + (embJet->Phi() - maxPartPhi) * (embJet->Phi() - maxPartPhi));
 475
 476       fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt());
 477       fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fRhoVal);
 478       fHistLeadPartPhiEta->Fill(maxPartEta, maxPartPhi);
 479       fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
 480
 481       fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), maxPartPt);
 482       fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area(), maxPartPt);
 483       fHistEmbJetsPhiEta->Fill(embJet->Eta(), embJet->Phi());
 484       fHistJetPtvsJetCorrPt[fCentBin]->Fill(embJet->Pt(), embJet->Pt() - fRhoVal * embJet->Area());
 485
 486       fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->MCPt());
 487       fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * embJet->Area(), embJet->Pt() - embJet->MCPt());
 488       fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
 489
 490       embJet = NextEmbeddedJet();
 491     }
 492
 493     if (countEmbJets==0) {
 494       AliDebug(1,"No embedded jets found!");
 495       if (fEmbTracks) {
 496         DoEmbTrackLoop();
 497         for (Int_t i = 0; i < fEmbeddedTrackNIds; i++) {
 498           AliDebug(2,Form("Embedded track %d found!",i));
 499           AliVParticle *track2 = static_cast<AliVParticle*>(fEmbTracks->At(fEmbeddedTrackIds[i]));
 500           if (!track2) continue;
 501           fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
 502           fHistEmbNotFoundPt[fCentBin]->Fill(track2->Pt());
 503         }
 504       }
 505
 506       if (fEmbCaloClusters) {
 507         DoEmbClusterLoop();
 508         for (Int_t i = 0; i < fEmbeddedClusterNIds; i++) {
 509           AliDebug(2,Form("Embedded cluster %d found!",i));
 510           AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterIds[i]));
 511           TLorentzVector nPart;
 512           cluster2->GetMomentum(nPart, fVertex);
 513           fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
 514           fHistEmbNotFoundPt[fCentBin]->Fill(nPart.Pt());
 515         }
 516       }
 517     }
 518   }
 519
 520   return kTRUE;
 521 }
 522
 523 //________________________________________________________________________
 524 void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
 525 {
 526   // Do track loop.
 527
 528   fEmbeddedTrackNIds = 0;
 529
 530   if (!fEmbTracks)
 531     return;
 532
 533   Int_t ntracks = fEmbTracks->GetEntriesFast();
 534
 535   for (Int_t i = 0; i < ntracks; i++) {
 536
 537     AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track
 538
 539     if (!track) {
 540       AliError(Form("Could not retrieve track %d",i));
 541       continue;
 542     }
 543
 544     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
 545
 546     if (vtrack && !AcceptTrack(vtrack))
 547       continue;
 548
 549     if (track->GetLabel() > 0) {
 550       fEmbeddedTrackIds[fEmbeddedTrackNIds] = i;
 551       fEmbeddedTrackNIds++;
 552     }
 553   }
 554 }
 555
 556 //________________________________________________________________________
 557 void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
 558 {
 559   // Do cluster loop.
 560
 561   fEmbeddedClusterNIds = 0;
 562
 563   if (!fEmbCaloClusters)
 564     return;
 565
 566   Int_t nclusters =  fEmbCaloClusters->GetEntriesFast();
 567
 568   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 569     AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
 570     if (!cluster) {
 571       AliError(Form("Could not receive cluster %d", iClusters));
 572       continue;
 573     }
 574
 575     if (!AcceptCluster(cluster))
 576       continue;
 577
 578     if (cluster->GetLabel() > 0) {
 579       fEmbeddedClusterIds[fEmbeddedClusterNIds] = iClusters;
 580       fEmbeddedClusterNIds++;
 581     }
 582   }
 583 }
 584
 585 //________________________________________________________________________
 586 AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Int_t i)
 587 {
 588   // Do the embedded jet loop.
 589
 590   static Int_t iJet = 0;
 591
 592   if (i >= 0)
 593     iJet = i;
 594   else
 595     iJet++;
 596
 597   if (!fEmbJets)
 598     return 0;
 599
 600   TLorentzVector maxClusVect;
 601
 602   const Int_t nembjets = fEmbJets->GetEntriesFast();
 603
 604   for (; iJet < nembjets; iJet++) {
 605
 606     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(iJet));
 607
 608     if (!jet) {
 609       AliError(Form("Could not receive jet %d", iJet));
 610       continue;
 611     }
 612
 613     if (jet->MCPt() < fMCJetPtThreshold)
 614       continue;
 615
 616     return jet;
 617   }
 618
 619   return 0;
 620 }
 621
 622 //________________________________________________________________________
 623 void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
 624                                         AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
 625 {
 626   // Get rigid cone.
 627
 628   if (!tracks)
 629     tracks = fTracks;
 630
 631   if (!clusters)
 632     clusters = fCaloClusters;
 633
 634   eta = 0;
 635   phi = 0;
 636   pt = 0;
 637
 638   if (!tracks && !clusters)
 639     return;
 640
 641   Float_t LJeta = 999;
 642   Float_t LJphi = 999;
 643
 644   if (jet) {
 645     LJeta = jet->Eta();
 646     LJphi = jet->Phi();
 647   }
 648
 649   Float_t maxEta = fJetMaxEta;
 650   Float_t minEta = fJetMinEta;
 651   Float_t maxPhi = fJetMaxPhi;
 652   Float_t minPhi = fJetMinPhi;
 653
 654   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 655   if (minPhi < 0) minPhi = 0;
 656
 657   Float_t dLJ = 0;
 658   Int_t repeats = 0;
 659   do {
 660     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
 661     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
 662     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 663     repeats++;
 664   } while (dLJ < fMinRC2LJ && repeats < 999);
 665
 666   if (repeats == 999) {
 667     AliWarning(Form("%s: Could not get random cone!", GetName()));
 668     return;
 669   }
 670
 671   if (clusters) {
 672     Int_t nclusters =  clusters->GetEntriesFast();
 673     for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 674       AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
 675       if (!cluster) {
 676         AliError(Form("Could not receive cluster %d", iClusters));
 677         continue;
 678       }
 679
 680       if (!AcceptCluster(cluster))
 681         continue;
 682
 683       TLorentzVector nPart;
 684       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 685
 686       Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
 687
 688       if (d <= fJetRadius)
 689         pt += nPart.Pt();
 690     }
 691   }
 692
 693   if (tracks) {
 694     Int_t ntracks = tracks->GetEntriesFast();
 695     for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 696       AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 697       if(!track) {
 698         AliError(Form("Could not retrieve track %d",iTracks));
 699         continue;
 700       }
 701
 702       if (!AcceptTrack(track))
 703         continue;
 704
 705       Float_t tracketa = track->Eta();
 706       Float_t trackphi = track->Phi();
 707
 708       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
 709         trackphi += 2 * TMath::Pi();
 710       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
 711         trackphi -= 2 * TMath::Pi();
 712
 713       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
 714       if (d <= fJetRadius)
 715         pt += track->Pt();
 716     }
 717   }
 718 }
 719
 720 //________________________________________________________________________
 721 void AliAnalysisTaskDeltaPt::ExecOnce()
 722 {
 723   // Initialize the analysis.
 724
 725   if (!fEmbJetsName.IsNull() && !fEmbJets) {
 726     fEmbJets =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
 727     if (!fEmbJets) {
 728       AliError(Form("%s: Could not retrieve embedded jets %s!", GetName(), fEmbJetsName.Data()));
 729       return;
 730     }
 731     else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
 732       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data()));
 733       fEmbJets = 0;
 734       return;
 735     }
 736   }
 737
 738   if (!fEmbCaloName.IsNull() && !fEmbCaloClusters) {
 739     fEmbCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
 740     if (!fEmbCaloClusters) {
 741       AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
 742       return;
 743     }
 744     else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
 745       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data()));
 746       fEmbCaloClusters = 0;
 747       return;
 748     }
 749   }
 750
 751   if (!fEmbTracksName.IsNull() && !fEmbTracks) {
 752     fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
 753     if (!fEmbTracks) {
 754       AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
 755       return;
 756     }
 757     else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
 758       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data()));
 759       fEmbTracks = 0;
 760       return;
 761     }
 762   }
 763
 764   if (!fRandCaloName.IsNull() && !fRandCaloClusters) {
 765     fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
 766     if (!fRandCaloClusters) {
 767       AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
 768       return;
 769     }
 770     else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
 771       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data()));
 772       fRandCaloClusters = 0;
 773       return;
 774     }
 775   }
 776
 777   if (!fRandTracksName.IsNull() && !fRandTracks) {
 778     fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
 779     if (!fRandTracks) {
 780       AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
 781       return;
 782     }
 783     else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
 784       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data()));
 785       fRandTracks = 0;
 786       return;
 787     }
 788   }
 789
 790   AliAnalysisTaskEmcalJet::ExecOnce();
 791
 792   if (fRCperEvent < 0) {
 793     Double_t area = (fJetMaxEta - fJetMinEta) * (fJetMaxPhi - fJetMinPhi);
 794     Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
 795     fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
 796     if (fRCperEvent == 0)
 797       fRCperEvent = 1;
 798   }
 799
 800   if (fMinRC2LJ < 0)
 801     fMinRC2LJ = fJetRadius * 1.5;
 802
 803   const Float_t maxDist = TMath::Max(fJetMaxPhi - fJetMinPhi, fJetMaxEta - fJetMinEta) / 2;
 804   if (fMinRC2LJ > maxDist) {
 805     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
 806                     "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
 807     fMinRC2LJ = maxDist;
 808   }
 809 }
 810
 811 //________________________________________________________________________
 812 void AliAnalysisTaskDeltaPt::Terminate(Option_t *)
 813 {
 814   // Called once at the end of the analysis.
 815 }