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     fHistEmbJetsPt[i] = 0;
  60     fHistEmbJetsCorrPt[i] = 0;
  61     fHistEmbJetsArea[i] = 0;
  62     fHistEmbPartPt[i] = 0;
  63     fHistDistEmbPartJetAxis[i] = 0;
  64     fHistRhoVSEmbBkg[i] = 0;
  65     fHistDeltaPtEmb[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     fHistEmbJetsPt[i] = 0;
 106     fHistEmbJetsCorrPt[i] = 0;
 107     fHistEmbJetsArea[i] = 0;
 108     fHistEmbPartPt[i] = 0;
 109     fHistDistEmbPartJetAxis[i] = 0;
 110     fHistRhoVSEmbBkg[i] = 0;
 111     fHistDeltaPtEmb[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 p_{T} [GeV/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 p_{T} [GeV/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 p_{T}^{RC} [GeV/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 p_{T}^{RC} [GeV/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*2, fNbins, fMinBinPt, fMaxBinPt);
 179     fHistRhoVSRCPt[i]->GetXaxis()->SetTitle("A#rho [GeV/c]");
 180     fHistRhoVSRCPt[i]->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/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("#deltap_{T}^{RC} [GeV/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("#deltap_{T}^{RC} [GeV/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("#deltap_{T}^{RC} [GeV/c]");
 201     fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
 202     fOutput->Add(fHistDeltaPtRCRand[i]);
 203
 204     if (!fEmbJetsName.IsNull()) {
 205       histname = "fHistEmbJetsPt_";
 206       histname += i;
 207       fHistEmbJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 208       fHistEmbJetsPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{raw} [GeV/c]");
 209       fHistEmbJetsPt[i]->GetYaxis()->SetTitle("counts");
 210       fOutput->Add(fHistEmbJetsPt[i]);
 211
 212       histname = "fHistEmbJetsCorrPt_";
 213       histname += i;
 214       fHistEmbJetsCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 215       fHistEmbJetsCorrPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{corr} [GeV/c]");
 216       fHistEmbJetsCorrPt[i]->GetYaxis()->SetTitle("counts");
 217       fOutput->Add(fHistEmbJetsCorrPt[i]);
 218
 219       histname = "fHistEmbJetsArea_";
 220       histname += i;
 221       fHistEmbJetsArea[i] = new TH1F(histname.Data(), histname.Data(), 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
 222       fHistEmbJetsArea[i]->GetXaxis()->SetTitle("area");
 223       fHistEmbJetsArea[i]->GetYaxis()->SetTitle("counts");
 224       fOutput->Add(fHistEmbJetsArea[i]);
 225
 226       histname = "fHistEmbPartPt_";
 227       histname += i;
 228       fHistEmbPartPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 229       fHistEmbPartPt[i]->GetXaxis()->SetTitle("embedded particle p_{T}^{emb} [GeV/c]");
 230       fHistEmbPartPt[i]->GetYaxis()->SetTitle("counts");
 231       fOutput->Add(fHistEmbPartPt[i]);
 232
 233       histname = "fHistDistEmbPartJetAxis_";
 234       histname += i;
 235       fHistDistEmbPartJetAxis[i] = new TH1F(histname.Data(), histname.Data(), 50, 0, 0.5);
 236       fHistDistEmbPartJetAxis[i]->GetXaxis()->SetTitle("distance");
 237       fHistDistEmbPartJetAxis[i]->GetYaxis()->SetTitle("counts");
 238       fOutput->Add(fHistDistEmbPartJetAxis[i]);
 239
 240       histname = "fHistEmbNotFoundPhiEta_";
 241       histname += i;
 242       fHistEmbNotFoundPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
 243       fHistEmbNotFoundPhiEta[i]->GetXaxis()->SetTitle("#eta");
 244       fHistEmbNotFoundPhiEta[i]->GetYaxis()->SetTitle("#phi");
 245       fOutput->Add(fHistEmbNotFoundPhiEta[i]);
 246
 247       histname = "fHistRhoVSEmbBkg_";
 248       histname += i;
 249       fHistRhoVSEmbBkg[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt*2, fNbins, fMinBinPt, fMaxBinPt);
 250       fHistRhoVSEmbBkg[i]->GetXaxis()->SetTitle("A#rho [GeV/c]");
 251       fHistRhoVSEmbBkg[i]->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
 252       fOutput->Add(fHistRhoVSEmbBkg[i]);
 253
 254       histname = "fHistDeltaPtEmb_";
 255       histname += i;
 256       fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 257       fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T}^{emb} [GeV/c]");
 258       fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
 259       fOutput->Add(fHistDeltaPtEmb[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   // Simple random cones
 285   for (Int_t i = 0; i < fRCperEvent; i++) {
 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, RCpt);
 293
 294       fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
 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 / rcArea);
 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 / rcArea);
 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     fHistEmbJetsPt[fCentBin]->Fill(embJet->Pt());
 371     fHistEmbJetsCorrPt[fCentBin]->Fill(embJet->Pt() - fRhoVal * embJet->Area());
 372     fHistEmbJetsArea[fCentBin]->Fill(embJet->Area());
 373     fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
 374
 375     fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRhoVal - probePt);
 376     fHistRhoVSEmbBkg[fCentBin]->Fill(fRhoVal, embJet->Pt() - probePt);
 377   }
 378   else {
 379     if (fEmbTracks)
 380       DoEmbTrackLoop();
 381     if (fEmbCaloClusters)
 382       DoEmbClusterLoop();
 383     if (fEmbeddedTrackId >= 0) {
 384       AliVTrack *track2 = static_cast<AliVTrack*>(fEmbTracks->At(fEmbeddedTrackId));
 385       fHistEmbNotFoundPhiEta[fCentBin]->Fill(track2->Eta(), track2->Phi());
 386     }
 387     else if (fEmbeddedClusterId >= 0) {
 388       AliVCluster *cluster2 = static_cast<AliVCluster*>(fEmbCaloClusters->At(fEmbeddedClusterId));
 389       TLorentzVector nPart;
 390       cluster2->GetMomentum(nPart, fVertex);
 391       fHistEmbNotFoundPhiEta[fCentBin]->Fill(nPart.Eta(), nPart.Phi());
 392     }
 393     else {
 394       AliWarning(Form("%s - Embedded particle not found!", GetName()));
 395     }
 396   }
 397
 398   return kTRUE;
 399 }
 400
 401 //________________________________________________________________________
 402 void AliAnalysisTaskDeltaPt::DoEmbTrackLoop()
 403 {
 404   // Do track loop.
 405
 406   if (!fEmbTracks)
 407     return;
 408
 409   fEmbeddedTrackId = -1;
 410
 411   Int_t ntracks = fEmbTracks->GetEntriesFast();
 412
 413   for (Int_t i = 0; i < ntracks; i++) {
 414
 415     AliVParticle* track = static_cast<AliVParticle*>(fEmbTracks->At(i)); // pointer to reconstructed to track
 416
 417     if (!track) {
 418       AliError(Form("Could not retrieve track %d",i));
 419       continue;
 420     }
 421
 422     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
 423
 424     if (vtrack && !AcceptTrack(vtrack, kTRUE))
 425       continue;
 426
 427     if (track->GetLabel() == 100) {
 428       fEmbeddedTrackId = i;
 429       continue;
 430     }
 431   }
 432 }
 433
 434 //________________________________________________________________________
 435 void AliAnalysisTaskDeltaPt::DoEmbClusterLoop()
 436 {
 437   // Do cluster loop.
 438
 439   if (!fEmbCaloClusters)
 440     return;
 441
 442   fEmbeddedClusterId = -1;
 443
 444   Int_t nclusters =  fEmbCaloClusters->GetEntriesFast();
 445
 446   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 447     AliVCluster* cluster = static_cast<AliVCluster*>(fEmbCaloClusters->At(iClusters));
 448     if (!cluster) {
 449       AliError(Form("Could not receive cluster %d", iClusters));
 450       continue;
 451     }
 452
 453     if (!AcceptCluster(cluster, kTRUE))
 454       continue;
 455
 456     if (cluster->Chi2() == 100) {
 457       fEmbeddedClusterId = iClusters;
 458       continue;
 459     }
 460   }
 461 }
 462
 463 //________________________________________________________________________
 464 void AliAnalysisTaskDeltaPt::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
 465 {
 466   // Do the embedded jet loop.
 467
 468   if (!fEmbJets)
 469     return;
 470
 471   TLorentzVector maxClusVect;
 472
 473   Int_t nembjets = fEmbJets->GetEntriesFast();
 474
 475   for (Int_t ij = 0; ij < nembjets; ij++) {
 476
 477     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
 478
 479     if (!jet) {
 480       AliError(Form("Could not receive jet %d", ij));
 481       continue;
 482     }
 483
 484     if (!AcceptJet(jet))
 485       continue;
 486
 487     if (!jet->IsMC())
 488       continue;
 489
 490     AliVParticle *maxTrack = 0;
 491
 492     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 493       AliVParticle *track = jet->TrackAt(it, fEmbTracks);
 494
 495       if (!track)
 496         continue;
 497
 498       if (track->GetLabel() != 100)
 499         continue;
 500
 501       if (!maxTrack || track->Pt() > maxTrack->Pt())
 502         maxTrack = track;
 503     }
 504
 505     AliVCluster *maxClus = 0;
 506
 507     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 508       AliVCluster *cluster = jet->ClusterAt(ic, fEmbCaloClusters);
 509
 510       if (!cluster)
 511         continue;
 512
 513       if (cluster->Chi2() != 100)
 514         continue;
 515
 516       TLorentzVector nPart;
 517       cluster->GetMomentum(nPart, fVertex);
 518
 519       if (!maxClus || nPart.Et() > maxClusVect.Et()) {
 520         new (&maxClusVect) TLorentzVector(nPart);
 521         maxClus = cluster;
 522       }
 523     }
 524
 525     if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
 526       embPart = maxClus;
 527       embJet = jet;
 528     }
 529     else if (maxTrack) {
 530       embPart = maxTrack;
 531       embJet = jet;
 532     }
 533
 534     return;  // MC jet found, exit
 535   }
 536 }
 537
 538 //________________________________________________________________________
 539 void AliAnalysisTaskDeltaPt::GetRandomCone(Float_t &pt, Float_t &eta, Float_t &phi,
 540                                         AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
 541 {
 542   // Get rigid cone.
 543
 544   if (!tracks)
 545     tracks = fTracks;
 546
 547   if (!clusters)
 548     clusters = fCaloClusters;
 549
 550   eta = 0;
 551   phi = 0;
 552   pt = 0;
 553
 554   if (!tracks && !clusters)
 555     return;
 556
 557   Float_t LJeta = 999;
 558   Float_t LJphi = 999;
 559
 560   if (jet) {
 561     LJeta = jet->Eta();
 562     LJphi = jet->Phi();
 563   }
 564
 565   Float_t maxEta = fJetMaxEta;
 566   Float_t minEta = fJetMinEta;
 567   Float_t maxPhi = fJetMaxPhi;
 568   Float_t minPhi = fJetMinPhi;
 569
 570   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 571   if (minPhi < 0) minPhi = 0;
 572
 573   Float_t dLJ = 0;
 574   Int_t repeats = 0;
 575   do {
 576     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
 577     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
 578     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 579     repeats++;
 580   } while (dLJ < fMinRC2LJ && repeats < 999);
 581
 582   if (repeats == 999) {
 583     AliWarning(Form("%s: Could not get random cone!", GetName()));
 584     return;
 585   }
 586
 587   if (clusters) {
 588     Int_t nclusters =  clusters->GetEntriesFast();
 589     for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 590       AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
 591       if (!cluster) {
 592         AliError(Form("Could not receive cluster %d", iClusters));
 593         continue;
 594       }
 595
 596       if (!AcceptCluster(cluster, fMCAna))
 597         continue;
 598
 599       TLorentzVector nPart;
 600       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 601
 602       Float_t d = TMath::Sqrt((nPart.Eta() - eta) * (nPart.Eta() - eta) + (nPart.Phi() - phi) * (nPart.Phi() - phi));
 603
 604       if (d <= fJetRadius)
 605         pt += nPart.Pt();
 606     }
 607   }
 608
 609   if (tracks) {
 610     Int_t ntracks = tracks->GetEntriesFast();
 611     for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 612       AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
 613       if(!track) {
 614         AliError(Form("Could not retrieve track %d",iTracks));
 615         continue;
 616       }
 617
 618       if (!AcceptTrack(track, fMCAna))
 619         continue;
 620
 621       Float_t tracketa = track->Eta();
 622       Float_t trackphi = track->Phi();
 623
 624       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
 625         trackphi += 2 * TMath::Pi();
 626       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
 627         trackphi -= 2 * TMath::Pi();
 628
 629       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
 630       if (d <= fJetRadius)
 631         pt += track->Pt();
 632     }
 633   }
 634 }
 635
 636 //________________________________________________________________________
 637 void AliAnalysisTaskDeltaPt::ExecOnce()
 638 {
 639   // Initialize the analysis.
 640
 641   if (!fEmbJetsName.IsNull() && !fEmbJets) {
 642     fEmbJets =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
 643     if (!fEmbJets) {
 644       AliError(Form("%s: Could not retrieve embedded jets %s!", GetName(), fEmbJetsName.Data()));
 645       return;
 646     }
 647     else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
 648       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data()));
 649       fEmbJets = 0;
 650       return;
 651     }
 652   }
 653
 654   if (!fEmbCaloName.IsNull() && !fEmbCaloClusters) {
 655     fEmbCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
 656     if (!fEmbCaloClusters) {
 657       AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
 658       return;
 659     }
 660     else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
 661       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data()));
 662       fEmbCaloClusters = 0;
 663       return;
 664     }
 665   }
 666
 667   if (!fEmbTracksName.IsNull() && !fEmbTracks) {
 668     fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
 669     if (!fEmbTracks) {
 670       AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
 671       return;
 672     }
 673     else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
 674       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data()));
 675       fEmbTracks = 0;
 676       return;
 677     }
 678   }
 679
 680   if (!fRandCaloName.IsNull() && !fRandCaloClusters) {
 681     fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
 682     if (!fRandCaloClusters) {
 683       AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
 684       return;
 685     }
 686     else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
 687       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data()));
 688       fRandCaloClusters = 0;
 689       return;
 690     }
 691   }
 692
 693   if (!fRandTracksName.IsNull() && !fRandTracks) {
 694     fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
 695     if (!fRandTracks) {
 696       AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
 697       return;
 698     }
 699     else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
 700       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data()));
 701       fRandTracks = 0;
 702       return;
 703     }
 704   }
 705
 706   AliAnalysisTaskEmcalJet::ExecOnce();
 707
 708   if (fRCperEvent < 0) {
 709     Double_t area = (fJetMaxEta - fJetMinEta) * (fJetMaxPhi - fJetMinPhi);
 710     Double_t jetArea = TMath::Pi() * fJetRadius * fJetRadius;
 711     fRCperEvent = TMath::FloorNint(area / jetArea - 0.5);
 712     if (fRCperEvent == 0)
 713       fRCperEvent = 1;
 714   }
 715
 716   if (fMinRC2LJ < 0)
 717     fMinRC2LJ = fJetRadius * 1.5;
 718
 719   const Float_t maxDist = TMath::Max(fJetMaxPhi - fJetMinPhi, fJetMaxEta - fJetMinEta) / 2;
 720   if (fMinRC2LJ > maxDist) {
 721     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
 722                     "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
 723     fMinRC2LJ = maxDist;
 724   }
 725 }
 726
 727 //________________________________________________________________________
 728 void AliAnalysisTaskDeltaPt::Terminate(Option_t *)
 729 {
 730   // Called once at the end of the analysis.
 731 }