PWGJE/EMCALJetTasks/AliAnalysisTaskDeltaPt.cxx

   1 // $Id$
   2 //
   3 // Jet deltaPt task.
   4 //
   5 // Author: S.Aiola
   6
   7 #include <TObject.h>
   8 #include <TChain.h>
   9 #include <TClonesArray.h>
  10 #include <TH1F.h>
  11 #include <TH2F.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14 #include <TRandom3.h>
  15 #include <TParameter.h>
  16
  17 #include "AliAnalysisManager.h"
  18 #include "AliCentrality.h"
  19 #include "AliVCluster.h"
  20 #include "AliVParticle.h"
  21 #include "AliVTrack.h"
  22 #include "AliEmcalJet.h"
  23 #include "AliVEventHandler.h"
  24 #include "AliRhoParameter.h"
  25 #include "AliLog.h"
  26
  27 #include "AliAnalysisTaskDeltaPt.h"
  28
  29 ClassImp(AliAnalysisTaskDeltaPt)
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt() :
  33   AliAnalysisTaskEmcalJet("AliAnalysisTaskDeltaPt", kTRUE),
  34   fMCAna(kFALSE),
  35   fMinRC2LJ(-1),
  36   fEmbJetsName(""),
  37   fEmbTracksName(""),
  38   fEmbCaloName(""),
  39   fRandTracksName("TracksRandomized"),
  40   fRandCaloName("CaloClustersRandomized"),
  41   fRCperEvent(-1),
  42   fEmbJets(0),
  43   fEmbTracks(0),
  44   fEmbCaloClusters(0),
  45   fRandTracks(0),
  46   fRandCaloClusters(0),
  47   fEmbeddedClusterId(-1),
  48   fEmbeddedTrackId(-1),
  49   fHistRCPhiEta(0),
  50   fHistRCPtExLJVSDPhiLJ(0),
  51   fHistRhoVSRCPt(0),
  52   fHistEmbJetPhiEta(0),
  53   fHistEmbPartPhiEta(0),
  54   fHistRhoVSEmbBkg(0)
  55 {
  56   // Default constructor.
  57
  58   for (Int_t i = 0; i < 4; i++) {
  59     fHistRCPt[i] = 0;
  60     fHistRCPtExLJ[i] = 0;
  61     fHistRCPtRand[i] = 0;
  62     fHistDeltaPtRC[i] = 0;
  63     fHistDeltaPtRCExLJ[i] = 0;
  64     fHistDeltaPtRCRand[i] = 0;
  65     fHistEmbNotFoundPhiEta[i] = 0;
  66     fHistEmbJetsPt[i] = 0;
  67     fHistEmbJetsCorrPt[i] = 0;
  68     fHistEmbJetsArea[i] = 0;
  69     fHistEmbPartPt[i] = 0;
  70     fHistDistEmbPartJetAxis[i] = 0;
  71     fHistDeltaPtEmb[i] = 0;
  72   }
  73
  74   SetMakeGeneralHistograms(kTRUE);
  75 }
  76
  77 //________________________________________________________________________
  78 AliAnalysisTaskDeltaPt::AliAnalysisTaskDeltaPt(const char *name) :
  79   AliAnalysisTaskEmcalJet(name, kTRUE),
  80   fMCAna(kFALSE),
  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   fEmbeddedClusterId(-1),
  94   fEmbeddedTrackId(-1),
  95   fHistRCPhiEta(0),
  96   fHistRCPtExLJVSDPhiLJ(0),
  97   fHistRhoVSRCPt(0),
  98   fHistEmbJetPhiEta(0),
  99   fHistEmbPartPhiEta(0),
 100   fHistRhoVSEmbBkg(0)
 101 {
 102   // Standard constructor.
 103
 104   for (Int_t i = 0; i < 4; i++) {
 105     fHistRCPt[i] = 0;
 106     fHistRCPtExLJ[i] = 0;
 107     fHistRCPtRand[i] = 0;
 108     fHistDeltaPtRC[i] = 0;
 109     fHistDeltaPtRCExLJ[i] = 0;
 110     fHistDeltaPtRCRand[i] = 0;
 111     fHistEmbNotFoundPhiEta[i] = 0;
 112     fHistEmbJetsPt[i] = 0;
 113     fHistEmbJetsCorrPt[i] = 0;
 114     fHistEmbJetsArea[i] = 0;
 115     fHistEmbPartPt[i] = 0;
 116     fHistDistEmbPartJetAxis[i] = 0;
 117     fHistDeltaPtEmb[i] = 0;
 118   }
 119
 120   SetMakeGeneralHistograms(kTRUE);
 121 }
 122
 123 //________________________________________________________________________
 124 AliAnalysisTaskDeltaPt::~AliAnalysisTaskDeltaPt()
 125 {
 126   // Destructor.
 127 }
 128
 129 //________________________________________________________________________
 130 void AliAnalysisTaskDeltaPt::UserCreateOutputObjects()
 131 {
 132   // Create user output.
 133
 134   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 135
 136   fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
 137   fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
 138   fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
 139   fOutput->Add(fHistRCPhiEta);
 140
 141   fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
 142   fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 143   fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
 144   fOutput->Add(fHistRCPtExLJVSDPhiLJ);
 145
 146   fHistRhoVSRCPt = new TH2F("fHistRhoVSRCPt","fHistRhoVSRCPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 147   fHistRhoVSRCPt->GetXaxis()->SetTitle("A#rho [GeV/c]");
 148   fHistRhoVSRCPt->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 149   fOutput->Add(fHistRhoVSRCPt);
 150
 151   if (!fJetsName.IsNull()) {
 152     fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
 153     fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
 154     fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
 155     fOutput->Add(fHistEmbJetPhiEta);
 156
 157     fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 50, -1, 1, 101, 0, TMath::Pi() * 2.02);
 158     fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
 159     fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
 160     fOutput->Add(fHistEmbPartPhiEta);
 161
 162     fHistRhoVSEmbBkg = new TH2F("fHistRhoVSEmbBkg","fHistRhoVSEmbBkg", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 163     fHistRhoVSEmbBkg->GetXaxis()->SetTitle("A#rho [GeV/c]");
 164     fHistRhoVSEmbBkg->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
 165     fOutput->Add(fHistRhoVSEmbBkg);
 166   }
 167
 168   TString histname;
 169
 170   for (Int_t i = 0; i < 4; i++) {
 171     histname = "fHistRCPt_";
 172     histname += i;
 173     fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 174     fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 175     fHistRCPt[i]->GetYaxis()->SetTitle("counts");
 176     fOutput->Add(fHistRCPt[i]);
 177
 178     histname = "fHistRCPtExLJ_";
 179     histname += i;
 180     fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 181     fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
 182     fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
 183     fOutput->Add(fHistRCPtExLJ[i]);
 184
 185     histname = "fHistRCPtRand_";
 186     histname += i;
 187     fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 188     fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
 189     fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
 190     fOutput->Add(fHistRCPtRand[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("#deltap_{T}^{RC} [GeV/c]");
 196     fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
 197     fOutput->Add(fHistDeltaPtRC[i]);
 198
 199     histname = "fHistDeltaPtRCExLJ_";
 200     histname += i;
 201     fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 202     fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
 203     fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
 204     fOutput->Add(fHistDeltaPtRCExLJ[i]);
 205
 206     histname = "fHistDeltaPtRCRand_";
 207     histname += i;
 208     fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 209     fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
 210     fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
 211     fOutput->Add(fHistDeltaPtRCRand[i]);
 212
 213     if (!fEmbJetsName.IsNull()) {
 214       histname = "fHistEmbJetsPt_";
 215       histname += i;
 216       fHistEmbJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 217       fHistEmbJetsPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{raw} [GeV/c]");
 218       fHistEmbJetsPt[i]->GetYaxis()->SetTitle("counts");
 219       fOutput->Add(fHistEmbJetsPt[i]);
 220
 221       histname = "fHistEmbJetsCorrPt_";
 222       histname += i;
 223       fHistEmbJetsCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 224       fHistEmbJetsCorrPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{corr} [GeV/c]");
 225       fHistEmbJetsCorrPt[i]->GetYaxis()->SetTitle("counts");
 226       fOutput->Add(fHistEmbJetsCorrPt[i]);
 227
 228       histname = "fHistEmbJetsArea_";
 229       histname += i;
 230       fHistEmbJetsArea[i] = new TH1F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
 231       fHistEmbJetsArea[i]->GetXaxis()->SetTitle("area");
 232       fHistEmbJetsArea[i]->GetYaxis()->SetTitle("counts");
 233       fOutput->Add(fHistEmbJetsArea[i]);
 234
 235       histname = "fHistEmbPartPt_";
 236       histname += i;
 237       fHistEmbPartPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 238       fHistEmbPartPt[i]->GetXaxis()->SetTitle("embedded particle p_{T}^{emb} [GeV/c]");
 239       fHistEmbPartPt[i]->GetYaxis()->SetTitle("counts");
 240       fOutput->Add(fHistEmbPartPt[i]);
 241
 242       histname = "fHistDistEmbPartJetAxis_";
 243       histname += i;
 244       fHistDistEmbPartJetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
 245       fHistDistEmbPartJetAxis[i]->GetXaxis()->SetTitle("distance");
 246       fHistDistEmbPartJetAxis[i]->GetYaxis()->SetTitle("counts");
 247       fOutput->Add(fHistDistEmbPartJetAxis[i]);
 248
 249       histname = "fHistEmbNotFoundPhiEta_";
 250       histname += i;
 251       fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
 252       fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
 253       fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
 254       fOutput->Add(fHistEmbNotFoundPhiEta[i]);
 255
 256       histname = "fHistDeltaPtEmb_";
 257       histname += i;
 258       fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 259       fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T}^{emb} [GeV/c]");
 260       fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
 261       fOutput->Add(fHistDeltaPtEmb[i]);
 262     }
 263   }
 264
 265   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 266 }
 267
 268 //________________________________________________________________________
 269 Bool_t AliAnalysisTaskDeltaPt::FillHistograms()
 270 {
 271   // Fill histograms.
 272
 273   Int_t *sortedJets = GetSortedArray(fJets);
 274
 275   AliEmcalJet* jet = 0;
 276
 277   if (sortedJets && sortedJets[0] > 0)
 278     jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
 279
 280   // ************
 281   // Random cones
 282   // _________________________________
 283
 284   const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
 285
 286   // Simple random cones
 287   for (Int_t i = 0; i < fRCperEvent; i++) {
 288     Float_t RCpt = 0;
 289     Float_t RCeta = 0;
 290     Float_t RCphi = 0;
 291     GetRandomCone(RCpt, RCeta, RCphi, 0);
 292     if (RCpt > 0) {
 293       fHistRCPhiEta->Fill(RCeta, RCphi);
 294       fHistRhoVSRCPt->Fill(fRhoVal * rcArea, RCpt);
 295
 296       fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
 297       fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 298     }
 299
 300     // Random cones far from leading jet
 301     RCpt = 0;
 302     RCeta = 0;
 303     RCphi = 0;
 304     GetRandomCone(RCpt, RCeta, RCphi, jet);
 305     if (RCpt > 0) {
 306       if (jet) {
 307         Float_t dphi = RCphi - jet->Phi();
 308         if (dphi > 4.8) dphi -= TMath::Pi() * 2;
 309         if (dphi < -1.6) dphi += TMath::Pi() * 2;
 310         fHistRCPtExLJVSDPhiLJ->Fill(RCpt, dphi);
 311       }
 312       fHistRCPtExLJ[fCentBin]->Fill(RCpt / rcArea);
 313       fHistDeltaPtRCExLJ[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 314     }
 315
 316     // Random cones with randomized particles
 317     RCpt = 0;
 318     RCeta = 0;
 319     RCphi = 0;
 320     GetRandomCone(RCpt, RCeta, RCphi, 0, fRandTracks, fRandCaloClusters);
 321     if (RCpt > 0) {
 322       fHistRCPtRand[fCentBin]->Fill(RCpt / rcArea);
 323       fHistDeltaPtRCRand[fCentBin]->Fill(RCpt - rcArea * fRhoVal);
 324     }
 325   }
 326
 327   // ************
 328   // Embedding
 329   // _________________________________
 330
 331   if (!fJets)
 332     return kTRUE;
 333
 334   AliEmcalJet *embJet  = 0;
 335   TObject     *embPart = 0;
 336
 337   DoEmbJetLoop(embJet, embPart);
 338
 339   if (embJet) {
 340     Double_t probePt = 0;
 341     Double_t probeEta = 0;
 342     Double_t probePhi = 0;
 343
 344     AliVCluster *cluster = dynamic_cast<AliVCluster*>(embPart);
 345     if (cluster) {
 346       TLorentzVector nPart;
 347       cluster->GetMomentum(nPart, fVertex);
 348
 349       probeEta = nPart.Eta();
 350       probePhi = nPart.Phi();
 351       probePt = nPart.Pt();
 352     }
 353     else {
 354       AliVTrack *track = dynamic_cast<AliVTrack*>(embPart);
 355       if (track) {
 356         probeEta = track->Eta();
 357         probePhi = track->Phi();
 358         probePt = track->Pt();
 359       }
 360       else {
 361         AliWarning(Form("%s - Embedded jet found but embedded particle not found (wrong type?)!", GetName()));
 362         return kTRUE;
 363       }
 364     }
 365     Double_t distProbeJet = TMath::Sqrt((embJet->Eta() - probeEta) * (embJet->Eta() - probeEta) + (embJet->Phi() - probePhi) * (embJet->Phi() - probePhi));
 366
 367     fHistEmbPartPt[fCentBin]->Fill(probePt);
 368     fHistEmbPartPhiEta->Fill(probeEta, probePhi);
 369
 370     fHistDistEmbPartJetAxis[fCentBin]->Fill(distProbeJet);
 371
 372     fHistEmbJetsPt[fCentBin]->Fill(embJet->Pt());
 373     fHistEmbJetsCorrPt[fCentBin]->Fill(embJet->Pt() - fRhoVal * embJet->Area());
 374     fHistEmbJetsArea[fCentBin]->Fill(embJet->Area());
 375     fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
 376
 377     fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRhoVal - probePt);
 378     fHistRhoVSEmbBkg->Fill(embJet->Area() * fRhoVal, embJet->Pt() - probePt);
 379   }
 380   else {
 381     if (fEmbTracks)
 382       DoEmbTrackLoop();
 383     if (fEmbCaloClusters)
 384       DoEmbClusterLoop();
 385     if (fEmbeddedTrackId >= 0) {
 386       AliVTrack *track2 = static_cast<AliVTrack*>(fEmbTracks->At(fEmbeddedTrackId));
 387       fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
 388     }
 389     else if (fEmbeddedClusterId >= 0) {
 390       AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterId));
 391       TLorentzVector nPart;
 392       cluster2->GetMomentum(nPart, fVertex);
 393       fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
 394     }
 395     else {
 396       AliWarning(Form("%s - Embedded particle not found!", GetName()));
 397     }
 398   }
 399
 400   return kTRUE;
 401 }
 402
 403 //________________________________________________________________________
 404 void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
 405 {
 406   // Do track loop.
 407
 408   if (!fEmbTracks)
 409     return;
 410
 411   fEmbeddedTrackId = -1;
 412
 413   Int_t ntracks = fEmbTracks->GetEntriesFast();
 414
 415   for (Int_t i = 0; i < ntracks; i++) {
 416
 417     AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track
 418
 419     if (!track) {
 420       AliError(Form("Could not retrieve track %d",i));
 421       continue;
 422     }
 423
 424     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
 425
 426     if (vtrack && !AcceptTrack(vtrack, kTRUE))
 427       continue;
 428
 429     if (track->GetLabel() == 100) {
 430       fEmbeddedTrackId = i;
 431       continue;
 432     }
 433   }
 434 }
 435
 436 //________________________________________________________________________
 437 void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
 438 {
 439   // Do cluster loop.
 440
 441   if (!fEmbCaloClusters)
 442     return;
 443
 444   fEmbeddedClusterId = -1;
 445
 446   Int_t nclusters =  fEmbCaloClusters->GetEntriesFast();
 447
 448   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 449     AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
 450     if (!cluster) {
 451       AliError(Form("Could not receive cluster %d", iClusters));
 452       continue;
 453     }
 454
 455     if (!AcceptCluster(cluster, kTRUE))
 456       continue;
 457
 458     if (cluster->Chi2() == 100) {
 459       fEmbeddedClusterId = iClusters;
 460       continue;
 461     }
 462   }
 463 }
 464
 465 //________________________________________________________________________
 466 void AliAnalysisTaskDeltaPt::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
 467 {
 468   // Do the embedded jet loop.
 469
 470   if (!fEmbJets)
 471     return;
 472
 473   TLorentzVector maxClusVect;
 474
 475   Int_t nembjets = fEmbJets->GetEntriesFast();
 476
 477   for (Int_t ij = 0; ij < nembjets; ij++) {
 478
 479     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
 480
 481     if (!jet) {
 482       AliError(Form("Could not receive jet %d", ij));
 483       continue;
 484     }
 485
 486     if (!AcceptJet(jet))
 487       continue;
 488
 489     if (!jet->IsMC())
 490       continue;
 491
 492     AliVParticle *maxTrack = 0;
 493
 494     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 495       AliVParticle *track = jet->TrackAt(it, fEmbTracks);
 496
 497       if (!track)
 498         continue;
 499
 500       if (track->GetLabel() != 100)
 501         continue;
 502
 503       if (!maxTrack || track->Pt() > maxTrack->Pt())
 504         maxTrack = track;
 505     }
 506
 507     AliVCluster *maxClus = 0;
 508
 509     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 510       AliVCluster *cluster = jet->ClusterAt(ic, fEmbCaloClusters);
 511
 512       if (!cluster)
 513         continue;
 514
 515       if (cluster->Chi2() != 100)
 516         continue;
 517
 518       TLorentzVector nPart;
 519       cluster->GetMomentum(nPart, fVertex);
 520
 521       if (!maxClus || nPart.Et() > maxClusVect.Et()) {
 522         new (&maxClusVect) TLorentzVector(nPart);
 523         maxClus = cluster;
 524       }
 525     }
 526
 527     if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
 528       embPart = maxClus;
 529       embJet = jet;
 530     }
 531     else if (maxTrack) {
 532       embPart = maxTrack;
 533       embJet = jet;
 534     }
 535
 536     return;  // MC jet found, exit
 537   }
 538 }
 539
 540 //________________________________________________________________________
 541 void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
 542                                         AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
 543 {
 544   // Get rigid cone.
 545
 546   if (!tracks)
 547     tracks = fTracks;
 548
 549   if (!clusters)
 550     clusters = fCaloClusters;
 551
 552   eta = 0;
 553   phi = 0;
 554   pt = 0;
 555
 556   if (!tracks && !clusters)
 557     return;
 558
 559   Float_t LJeta = 999;
 560   Float_t LJphi = 999;
 561
 562   if (jet) {
 563     LJeta = jet->Eta();
 564     LJphi = jet->Phi();
 565   }
 566
 567   Float_t maxEta = fJetMaxEta;
 568   Float_t minEta = fJetMinEta;
 569   Float_t maxPhi = fJetMaxPhi;
 570   Float_t minPhi = fJetMinPhi;
 571
 572   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 573   if (minPhi < 0) minPhi = 0;
 574
 575   Float_t dLJ = 0;
 576   Int_t repeats = 0;
 577   do {
 578     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
 579     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
 580     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 581     repeats++;
 582   } while (dLJ < fMinRC2LJ && repeats < 999);
 583
 584   if (repeats == 999) {
 585     AliWarning(Form("%s: Could not get random cone!", GetName()));
 586     return;
 587   }
 588
 589   if (clusters) {
 590     Int_t nclusters =  clusters->GetEntriesFast();
 591     for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 592       AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
 593       if (!cluster) {
 594         AliError(Form("Could not receive cluster %d", iClusters));
 595         continue;
 596       }
 597
 598       if (!AcceptCluster(cluster, fMCAna))
 599         continue;
 600
 601       TLorentzVector nPart;
 602       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 603
 604       Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
 605
 606       if (d <= fJetRadius)
 607         pt += nPart.Pt();
 608     }
 609   }
 610
 611   if (tracks) {
 612     Int_t ntracks = tracks->GetEntriesFast();
 613     for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 614       AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 615       if(!track) {
 616         AliError(Form("Could not retrieve track %d",iTracks));
 617         continue;
 618       }
 619
 620       if (!AcceptTrack(track, fMCAna))
 621         continue;
 622
 623       Float_t tracketa = track->Eta();
 624       Float_t trackphi = track->Phi();
 625
 626       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
 627         trackphi += 2 * TMath::Pi();
 628       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
 629         trackphi -= 2 * TMath::Pi();
 630
 631       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
 632       if (d <= fJetRadius)
 633         pt += track->Pt();
 634     }
 635   }
 636 }
 637
 638 //________________________________________________________________________
 639 void AliAnalysisTaskDeltaPt::ExecOnce()
 640 {
 641   // Initialize the analysis.
 642
 643   if (!fEmbJetsName.IsNull() && !fEmbJets) {
 644     fEmbJets =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
 645     if (!fEmbJets) {
 646       AliError(Form("%s: Could not retrieve embedded jets %s!", GetName(), fEmbJetsName.Data()));
 647       return;
 648     }
 649     else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
 650       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data()));
 651       fEmbJets = 0;
 652       return;
 653     }
 654   }
 655
 656   if (!fEmbCaloName.IsNull() && !fEmbCaloClusters) {
 657     fEmbCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
 658     if (!fEmbCaloClusters) {
 659       AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
 660       return;
 661     }
 662     else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
 663       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data()));
 664       fEmbCaloClusters = 0;
 665       return;
 666     }
 667   }
 668
 669   if (!fEmbTracksName.IsNull() && !fEmbTracks) {
 670     fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
 671     if (!fEmbTracks) {
 672       AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
 673       return;
 674     }
 675     else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
 676       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data()));
 677       fEmbTracks = 0;
 678       return;
 679     }
 680   }
 681
 682   if (!fRandCaloName.IsNull() && !fRandCaloClusters) {
 683     fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
 684     if (!fRandCaloClusters) {
 685       AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
 686       return;
 687     }
 688     else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
 689       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data()));
 690       fRandCaloClusters = 0;
 691       return;
 692     }
 693   }
 694
 695   if (!fRandTracksName.IsNull() && !fRandTracks) {
 696     fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
 697     if (!fRandTracks) {
 698       AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
 699       return;
 700     }
 701     else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
 702       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data()));
 703       fRandTracks = 0;
 704       return;
 705     }
 706   }
 707
 708   AliAnalysisTaskEmcalJet::ExecOnce();
 709
 710   if (fRCperEvent < 0) {
 711     Double_t area = (fJetMaxEta - fJetMinEta) * (fJetMaxPhi - fJetMinPhi);
 712     Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
 713     fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
 714     if (fRCperEvent == 0)
 715       fRCperEvent = 1;
 716   }
 717
 718   if (fMinRC2LJ < 0)
 719     fMinRC2LJ = fJetRadius * 1.5;
 720
 721   const Float_t maxDist = TMath::Max(fJetMaxPhi - fJetMinPhi, fJetMaxEta - fJetMinEta) / 2;
 722   if (fMinRC2LJ > maxDist) {
 723     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
 724                     "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
 725     fMinRC2LJ = maxDist;
 726   }
 727 }
 728
 729 //________________________________________________________________________
 730 void AliAnalysisTaskDeltaPt::Terminate(Option_t *)
 731 {
 732   // Called once at the end of the analysis.
 733 }