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   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 330 }
 331
 332 //________________________________________________________________________
 333 Bool_t AliAnalysisTaskDeltaPt::FillHistograms()
 334 {
 335   // Fill histograms.
 336
 337   // ************
 338   // Random cones
 339   // _________________________________
 340
 341   const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
 342   Float_t RCpt = 0;
 343   Float_t RCeta = 0;
 344   Float_t RCphi = 0;
 345
 346   if (fTracks || fCaloClusters) {
 347
 348     for (Int_t i = 0; i < fRCperEvent; i++) {
 349       // Simple random cones
 350       RCpt = 0;
 351       RCeta = 0;
 352       RCphi = 0;
 353       GetRandomCone(RCpt, RCeta, RCphi, 0);
 354       if (RCpt > 0) {
 355         fHistRCPhiEta->Fill(RCeta, RCphi);
 356         fHistRhoVSRCPt[fCentBin]->Fill(fRhoVal * rcArea, RCpt);
 357
 358         fHistRCPt[fCentBin]->Fill(RCpt);
 359         fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 360       }
 361
 362       if (fJets) {
 363
 364         // Random cones far from leading jet
 365         static Int_t sortedJets[9999] = {-1};
 366         GetSortedArray(sortedJets, fJets);
 367
 368         AliEmcalJet* jet = 0;
 369
 370         if (sortedJets[0] >= 0)
 371           jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
 372
 373         RCpt = 0;
 374         RCeta = 0;
 375         RCphi = 0;
 376         GetRandomCone(RCpt, RCeta, RCphi, jet);
 377         if (RCpt > 0) {
 378           if (jet) {
 379             Float_t dphi = RCphi - jet->Phi();
 380             if (dphi > 4.8) dphi -= TMath::Pi() * 2;
 381             if (dphi < -1.6) dphi += TMath::Pi() * 2;
 382             fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
 383           }
 384           fHistRCPtExLJ[fCentBin]->Fill(RCpt);
 385           fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 386         }
 387       }
 388     }
 389   }
 390
 391   // Random cones with randomized particles
 392   if (fRandTracks || fRandCaloClusters) {
 393     RCpt = 0;
 394     RCeta = 0;
 395     RCphi = 0;
 396     GetRandomCone(RCpt, RCeta, RCphi, 0, fRandTracks, fRandCaloClusters);
 397     if (RCpt > 0) {
 398       fHistRCPtRand[fCentBin]->Fill(RCpt);
 399       fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 400     }
 401   }
 402
 403   // ************
 404   // Embedding
 405   // _________________________________
 406
 407   if (fEmbJets) {
 408
 409     AliEmcalJet *embJet = NextEmbeddedJet(0);
 410
 411     Int_t countEmbJets = 0;
 412
 413     while (embJet != 0) {
 414       AliDebug(2,Form("Elaborating embedded jet n. %d", countEmbJets));
 415       countEmbJets++;
 416
 417       if (!AcceptJet(embJet)) {
 418         AliDebug(2,"Embedded jet not accepted, skipping...");
 419         fHistEmbRejectedJetsPhiEta[fCentBin]->Fill(embJet->Eta(), embJet->Phi());
 420         fHistEmbRejectedJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt());
 421
 422         embJet = NextEmbeddedJet();
 423         continue;
 424       }
 425
 426       Double_t maxClusterPt = 0;
 427       Double_t maxClusterEta = 0;
 428       Double_t maxClusterPhi = 0;
 429
 430       Double_t maxTrackPt = 0;
 431       Double_t maxTrackEta = 0;
 432       Double_t maxTrackPhi = 0;
 433
 434       Double_t maxPartPt = 0;
 435       Double_t maxPartEta = 0;
 436       Double_t maxPartPhi = 0;
 437
 438       if (fLeadingHadronType == 1 || fLeadingHadronType == 2) {
 439         AliVCluster *cluster = embJet->GetLeadingCluster(fEmbCaloClusters);
 440         if (cluster) {
 441           TLorentzVector nPart;
 442           cluster->GetMomentum(nPart, fVertex);
 443
 444           maxClusterEta = nPart.Eta();
 445           maxClusterPhi = nPart.Phi();
 446           maxClusterPt = nPart.Pt();
 447         }
 448       }
 449
 450       if (fLeadingHadronType == 0 || fLeadingHadronType == 2) {
 451         AliVParticle *track = embJet->GetLeadingTrack(fEmbTracks);
 452         if (track) {
 453           maxTrackEta = track->Eta();
 454           maxTrackPhi = track->Phi();
 455           maxTrackPt = track->Pt();
 456         }
 457       }
 458
 459       if (maxTrackPt > maxClusterPt) {
 460         maxPartPt = maxTrackPt;
 461         maxPartEta = maxTrackEta;
 462         maxPartPhi = maxTrackPhi;
 463       }
 464       else {
 465         maxPartPt = maxClusterPt;
 466         maxPartEta = maxClusterEta;
 467         maxPartPhi = maxClusterPhi;
 468       }
 469
 470       Double_t distLeading2Jet = TMath::Sqrt((embJet->Eta() - maxPartEta) * (embJet->Eta() - maxPartEta) + (embJet->Phi() - maxPartPhi) * (embJet->Phi() - maxPartPhi));
 471
 472       fHistEmbPartPtvsJetPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt());
 473       fHistEmbPartPtvsJetCorrPt[fCentBin]->Fill(embJet->MCPt(), embJet->Pt() - embJet->Area() * fRhoVal);
 474       fHistLeadPartPhiEta->Fill(maxPartEta, maxPartPhi);
 475       fHistDistLeadPart2JetAxis[fCentBin]->Fill(distLeading2Jet);
 476
 477       fHistEmbJetsPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt(), maxPartPt);
 478       fHistEmbJetsCorrPtArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - fRhoVal * embJet->Area(), maxPartPt);
 479       fHistEmbJetsPhiEta->Fill(embJet->Eta(), embJet->Phi());
 480       fHistJetPtvsJetCorrPt[fCentBin]->Fill(embJet->Pt(), embJet->Pt() - fRhoVal * embJet->Area());
 481
 482       fHistEmbBkgArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->MCPt());
 483       fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal * embJet->Area(), embJet->Pt() - embJet->MCPt());
 484       fHistDeltaPtEmbArea[fCentBin]->Fill(embJet->Area(), embJet->Pt() - embJet->Area() * fRhoVal - embJet->MCPt());
 485
 486       embJet = NextEmbeddedJet();
 487     }
 488
 489     if (countEmbJets==0) {
 490       AliDebug(1,"No embedded jets found!");
 491       if (fEmbTracks) {
 492         DoEmbTrackLoop();
 493         for (Int_t i = 0; i < fEmbeddedTrackNIds; i++) {
 494           AliDebug(2,Form("Embedded track %d found!",i));
 495           AliVParticle *track2 = static_cast<AliVParticle*>(fEmbTracks->At(fEmbeddedTrackIds[i]));
 496           if (!track2) continue;
 497           fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
 498           fHistEmbNotFoundPt[fCentBin]->Fill(track2->Pt());
 499         }
 500       }
 501
 502       if (fEmbCaloClusters) {
 503         DoEmbClusterLoop();
 504         for (Int_t i = 0; i < fEmbeddedClusterNIds; i++) {
 505           AliDebug(2,Form("Embedded cluster %d found!",i));
 506           AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterIds[i]));
 507           TLorentzVector nPart;
 508           cluster2->GetMomentum(nPart, fVertex);
 509           fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
 510           fHistEmbNotFoundPt[fCentBin]->Fill(nPart.Pt());
 511         }
 512       }
 513     }
 514   }
 515
 516   return kTRUE;
 517 }
 518
 519 //________________________________________________________________________
 520 void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
 521 {
 522   // Do track loop.
 523
 524   fEmbeddedTrackNIds = 0;
 525
 526   if (!fEmbTracks)
 527     return;
 528
 529   Int_t ntracks = fEmbTracks->GetEntriesFast();
 530
 531   for (Int_t i = 0; i < ntracks; i++) {
 532
 533     AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track
 534
 535     if (!track) {
 536       AliError(Form("Could not retrieve track %d",i));
 537       continue;
 538     }
 539
 540     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
 541
 542     if (vtrack && !AcceptTrack(vtrack))
 543       continue;
 544
 545     if (track->GetLabel() > 0) {
 546       fEmbeddedTrackIds[fEmbeddedTrackNIds] = i;
 547       fEmbeddedTrackNIds++;
 548     }
 549   }
 550 }
 551
 552 //________________________________________________________________________
 553 void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
 554 {
 555   // Do cluster loop.
 556
 557   fEmbeddedClusterNIds = 0;
 558
 559   if (!fEmbCaloClusters)
 560     return;
 561
 562   Int_t nclusters =  fEmbCaloClusters->GetEntriesFast();
 563
 564   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 565     AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
 566     if (!cluster) {
 567       AliError(Form("Could not receive cluster %d", iClusters));
 568       continue;
 569     }
 570
 571     if (!AcceptCluster(cluster))
 572       continue;
 573
 574     if (cluster->GetLabel() > 0) {
 575       fEmbeddedClusterIds[fEmbeddedClusterNIds] = iClusters;
 576       fEmbeddedClusterNIds++;
 577     }
 578   }
 579 }
 580
 581 //________________________________________________________________________
 582 AliEmcalJet* AliAnalysisTaskDeltaPt::NextEmbeddedJet(Int_t i)
 583 {
 584   // Do the embedded jet loop.
 585
 586   static Int_t iJet = 0;
 587
 588   if (i >= 0)
 589     iJet = i;
 590   else
 591     iJet++;
 592
 593   if (!fEmbJets)
 594     return 0;
 595
 596   TLorentzVector maxClusVect;
 597
 598   const Int_t nembjets = fEmbJets->GetEntriesFast();
 599
 600   for (; iJet < nembjets; iJet++) {
 601
 602     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(iJet));
 603
 604     if (!jet) {
 605       AliError(Form("Could not receive jet %d", iJet));
 606       continue;
 607     }
 608
 609     if (jet->MCPt() < fMCJetPtThreshold)
 610       continue;
 611
 612     return jet;
 613   }
 614
 615   return 0;
 616 }
 617
 618 //________________________________________________________________________
 619 void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
 620                                         AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
 621 {
 622   // Get rigid cone.
 623
 624   if (!tracks)
 625     tracks = fTracks;
 626
 627   if (!clusters)
 628     clusters = fCaloClusters;
 629
 630   eta = 0;
 631   phi = 0;
 632   pt = 0;
 633
 634   if (!tracks && !clusters)
 635     return;
 636
 637   Float_t LJeta = 999;
 638   Float_t LJphi = 999;
 639
 640   if (jet) {
 641     LJeta = jet->Eta();
 642     LJphi = jet->Phi();
 643   }
 644
 645   Float_t maxEta = fJetMaxEta;
 646   Float_t minEta = fJetMinEta;
 647   Float_t maxPhi = fJetMaxPhi;
 648   Float_t minPhi = fJetMinPhi;
 649
 650   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 651   if (minPhi < 0) minPhi = 0;
 652
 653   Float_t dLJ = 0;
 654   Int_t repeats = 0;
 655   do {
 656     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
 657     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
 658     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 659     repeats++;
 660   } while (dLJ < fMinRC2LJ && repeats < 999);
 661
 662   if (repeats == 999) {
 663     AliWarning(Form("%s: Could not get random cone!", GetName()));
 664     return;
 665   }
 666
 667   if (clusters) {
 668     Int_t nclusters =  clusters->GetEntriesFast();
 669     for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 670       AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
 671       if (!cluster) {
 672         AliError(Form("Could not receive cluster %d", iClusters));
 673         continue;
 674       }
 675
 676       if (!AcceptCluster(cluster))
 677         continue;
 678
 679       TLorentzVector nPart;
 680       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 681
 682       Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
 683
 684       if (d <= fJetRadius)
 685         pt += nPart.Pt();
 686     }
 687   }
 688
 689   if (tracks) {
 690     Int_t ntracks = tracks->GetEntriesFast();
 691     for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 692       AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 693       if(!track) {
 694         AliError(Form("Could not retrieve track %d",iTracks));
 695         continue;
 696       }
 697
 698       if (!AcceptTrack(track))
 699         continue;
 700
 701       Float_t tracketa = track->Eta();
 702       Float_t trackphi = track->Phi();
 703
 704       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
 705         trackphi += 2 * TMath::Pi();
 706       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
 707         trackphi -= 2 * TMath::Pi();
 708
 709       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
 710       if (d <= fJetRadius)
 711         pt += track->Pt();
 712     }
 713   }
 714 }
 715
 716 //________________________________________________________________________
 717 void AliAnalysisTaskDeltaPt::ExecOnce()
 718 {
 719   // Initialize the analysis.
 720
 721   if (!fEmbJetsName.IsNull() && !fEmbJets) {
 722     fEmbJets =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
 723     if (!fEmbJets) {
 724       AliError(Form("%s: Could not retrieve embedded jets %s!", GetName(), fEmbJetsName.Data()));
 725       return;
 726     }
 727     else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
 728       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data()));
 729       fEmbJets = 0;
 730       return;
 731     }
 732   }
 733
 734   if (!fEmbCaloName.IsNull() && !fEmbCaloClusters) {
 735     fEmbCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
 736     if (!fEmbCaloClusters) {
 737       AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
 738       return;
 739     }
 740     else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
 741       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data()));
 742       fEmbCaloClusters = 0;
 743       return;
 744     }
 745   }
 746
 747   if (!fEmbTracksName.IsNull() && !fEmbTracks) {
 748     fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
 749     if (!fEmbTracks) {
 750       AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
 751       return;
 752     }
 753     else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
 754       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data()));
 755       fEmbTracks = 0;
 756       return;
 757     }
 758   }
 759
 760   if (!fRandCaloName.IsNull() && !fRandCaloClusters) {
 761     fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
 762     if (!fRandCaloClusters) {
 763       AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
 764       return;
 765     }
 766     else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
 767       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data()));
 768       fRandCaloClusters = 0;
 769       return;
 770     }
 771   }
 772
 773   if (!fRandTracksName.IsNull() && !fRandTracks) {
 774     fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
 775     if (!fRandTracks) {
 776       AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
 777       return;
 778     }
 779     else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
 780       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data()));
 781       fRandTracks = 0;
 782       return;
 783     }
 784   }
 785
 786   AliAnalysisTaskEmcalJet::ExecOnce();
 787
 788   if (fRCperEvent < 0) {
 789     Double_t area = (fJetMaxEta - fJetMinEta) * (fJetMaxPhi - fJetMinPhi);
 790     Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
 791     fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
 792     if (fRCperEvent == 0)
 793       fRCperEvent = 1;
 794   }
 795
 796   if (fMinRC2LJ < 0)
 797     fMinRC2LJ = fJetRadius * 1.5;
 798
 799   const Float_t maxDist = TMath::Max(fJetMaxPhi - fJetMinPhi, fJetMaxEta - fJetMinEta) / 2;
 800   if (fMinRC2LJ > maxDist) {
 801     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
 802                     "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
 803     fMinRC2LJ = maxDist;
 804   }
 805 }
 806
 807 //________________________________________________________________________
 808 void AliAnalysisTaskDeltaPt::Terminate(Option_t *)
 809 {
 810   // Called once at the end of the analysis.
 811 }