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