PWGGA/EMCALJetTasks/AliAnalysisTaskSAJF.cxx

   1 // $Id$
   2 //
   3 // Jet analysis 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 "AliLog.h"
  25
  26 #include "AliAnalysisTaskSAJF.h"
  27
  28 ClassImp(AliAnalysisTaskSAJF)
  29
  30 //________________________________________________________________________
  31 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
  32   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAJF"),
  33   fMinRC2LJ(1.0),
  34   fEmbJetsName("EmbJets"),
  35   fRandTracksName("TracksRandomized"),
  36   fRandCaloName("CaloClustersRandomized"),
  37   fRhoName("Rho"),
  38   fSkipEventsWithNoJets(kTRUE),
  39   fEmbJets(0),
  40   fRandTracks(0),
  41   fRandCaloClusters(0),
  42   fRho(0),
  43   fHistCentrality(0),
  44   fHistRejectedEvents(0),
  45   fHistRhoVSleadJetPt(0),
  46   fHistRCPhiEta(0),
  47   fHistRCPtExLJVSDPhiLJ(0),
  48   fHistRhoVSRCPt(0),
  49   fHistEmbJetPhiEta(0),
  50   fHistEmbPartPhiEta(0),
  51   fHistRhoVSEmbBkg(0)
  52
  53 {
  54   // Default constructor.
  55
  56   for (Int_t i = 0; i < 4; i++) {
  57     fHistJetPhiEta[i] = 0;
  58     fHistJetsPt[i] = 0;
  59     fHistJetsPtArea[i] = 0;
  60     fHistLeadingJetPt[i] = 0;
  61     fHist2LeadingJetPt[i] = 0;
  62     fHistJetsNEFvsPt[i] = 0;
  63     fHistJetsZvsPt[i] = 0;
  64     fHistRho[i] = 0;
  65     fHistCorrJetsPt[i] = 0;
  66     fHistCorrLeadingJetPt[i] = 0;
  67     fHistRCPt[i] = 0;
  68     fHistRCPtExLJ[i] = 0;
  69     fHistRCPtRand[i] = 0;
  70     fHistDeltaPtRC[i] = 0;
  71     fHistDeltaPtRCExLJ[i] = 0;
  72     fHistDeltaPtRCRand[i] = 0;
  73     fHistEmbJets[i] = 0;
  74     fHistEmbPart[i] = 0;
  75     fHistDeltaPtEmb[i] = 0;
  76   }
  77 }
  78
  79 //________________________________________________________________________
  80 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
  81   AliAnalysisTaskEmcalJet(name),
  82   fMinRC2LJ(1.0),
  83   fEmbJetsName("EmbJets"),
  84   fRandTracksName("TracksRandomized"),
  85   fRandCaloName("CaloClustersRandomized"),
  86   fRhoName("Rho"),
  87   fSkipEventsWithNoJets(kTRUE),
  88   fEmbJets(0),
  89   fRandTracks(0),
  90   fRandCaloClusters(0),
  91   fRho(0),
  92   fHistCentrality(0),
  93   fHistRejectedEvents(0),
  94   fHistRhoVSleadJetPt(0),
  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     fHistJetPhiEta[i] = 0;
 106     fHistJetsPt[i] = 0;
 107     fHistJetsPtArea[i] = 0;
 108     fHistLeadingJetPt[i] = 0;
 109     fHist2LeadingJetPt[i] = 0;
 110     fHistJetsNEFvsPt[i] = 0;
 111     fHistJetsZvsPt[i] = 0;
 112     fHistRho[i] = 0;
 113     fHistCorrJetsPt[i] = 0;
 114     fHistCorrLeadingJetPt[i] = 0;
 115     fHistRCPt[i] = 0;
 116     fHistRCPtExLJ[i] = 0;
 117     fHistRCPtRand[i] = 0;
 118     fHistDeltaPtRC[i] = 0;
 119     fHistDeltaPtRCExLJ[i] = 0;
 120     fHistDeltaPtRCRand[i] = 0;
 121     fHistEmbJets[i] = 0;
 122     fHistEmbPart[i] = 0;
 123     fHistDeltaPtEmb[i] = 0;
 124   }
 125 }
 126
 127 //________________________________________________________________________
 128 AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
 129 {
 130   // Destructor.
 131 }
 132
 133 //________________________________________________________________________
 134 void AliAnalysisTaskSAJF::UserCreateOutputObjects()
 135 {
 136   // Create user output.
 137
 138   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 139
 140   const Float_t binWidth = (fMaxBinPt - fMinBinPt) / fNbins;
 141
 142   OpenFile(1);
 143   fOutput = new TList();
 144   fOutput->SetOwner();
 145
 146   fHistCentrality = new TH1F("fHistCentrality","fHistCentrality", fNbins, 0, 100);
 147   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
 148   fHistCentrality->GetYaxis()->SetTitle("counts");
 149   fOutput->Add(fHistCentrality);
 150
 151   fHistRejectedEvents = new TH1F("fHistRejectedEvents","fHistRejectedEvents", 6, 0, 6);
 152   fHistRejectedEvents->GetXaxis()->SetTitle("Reason");
 153   fHistRejectedEvents->GetYaxis()->SetTitle("counts");
 154   fHistRejectedEvents->GetXaxis()->SetBinLabel(1, "Rho <= 0");
 155   fHistRejectedEvents->GetXaxis()->SetBinLabel(2, "Max Jet <= 0");
 156   fHistRejectedEvents->GetXaxis()->SetBinLabel(3, "Max Jet not found");
 157   fHistRejectedEvents->GetXaxis()->SetBinLabel(4, "No jets");
 158   fOutput->Add(fHistRejectedEvents);
 159
 160   fHistRhoVSleadJetPt = new TH2F("fHistRhoVSleadJetPt","fHistRhoVSleadJetPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 161   fHistRhoVSleadJetPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
 162   fHistRhoVSleadJetPt->GetYaxis()->SetTitle("Leading jet p_{T} [GeV/c]");
 163   fOutput->Add(fHistRhoVSleadJetPt);
 164
 165   fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 40, -2, 2, 64, 0, 6.4);
 166   fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
 167   fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
 168   fOutput->Add(fHistRCPhiEta);
 169
 170   fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
 171   fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 172   fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
 173   fOutput->Add(fHistRCPtExLJVSDPhiLJ);
 174
 175   fHistRhoVSRCPt = new TH2F("fHistRhoVSRCPt","fHistRhoVSRCPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 176   fHistRhoVSRCPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
 177   fHistRhoVSRCPt->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 178   fOutput->Add(fHistRhoVSRCPt);
 179
 180   fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 40, -2, 2, 64, 0, 6.4);
 181   fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
 182   fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
 183   fOutput->Add(fHistEmbJetPhiEta);
 184
 185   fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 40, -2, 2, 64, 0, 6.4);
 186   fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
 187   fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
 188   fOutput->Add(fHistEmbPartPhiEta);
 189
 190   fHistRhoVSEmbBkg = new TH2F("fHistRhoVSEmbBkg","fHistRhoVSEmbBkg", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 191   fHistRhoVSEmbBkg->GetXaxis()->SetTitle("rho * area [GeV/c]");
 192   fHistRhoVSEmbBkg->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
 193   fOutput->Add(fHistRhoVSEmbBkg);
 194
 195   TString histname;
 196
 197   for (Int_t i = 0; i < 4; i++) {
 198     histname = "fHistJetPhiEta_";
 199     histname += i;
 200     fHistJetPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
 201     fHistJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
 202     fHistJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
 203     fOutput->Add(fHistJetPhiEta[i]);
 204
 205     histname = "fHistJetsPt_";
 206     histname += i;
 207     fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 208     fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 209     fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
 210     fOutput->Add(fHistJetsPt[i]);
 211
 212     histname = "fHistJetsPtArea_";
 213     histname += i;
 214     fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
 215     fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 216     fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 217     fOutput->Add(fHistJetsPtArea[i]);
 218
 219     histname = "fHistLeadingJetPt_";
 220     histname += i;
 221     fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 222     fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV]");
 223     fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 224     fOutput->Add(fHistLeadingJetPt[i]);
 225
 226     histname = "fHist2LeadingJetPt_";
 227     histname += i;
 228     fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 229     fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV]");
 230     fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 231     fOutput->Add(fHist2LeadingJetPt[i]);
 232
 233     histname = "fHistJetsZvsPt_";
 234     histname += i;
 235     fHistJetsZvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 236     fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
 237     fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 238     fOutput->Add(fHistJetsZvsPt[i]);
 239
 240     if (fAnaType == kEMCAL) {
 241       histname = "fHistJetsNEFvsPt_";
 242       histname += i;
 243       fHistJetsNEFvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 244       fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
 245       fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 246       fOutput->Add(fHistJetsNEFvsPt[i]);
 247     }
 248
 249     histname = "fHistRho_";
 250     histname += i;
 251     fHistRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 252     fHistRho[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 253     fHistRho[i]->GetYaxis()->SetTitle("counts");
 254     fOutput->Add(fHistRho[i]);
 255
 256     histname = "fHistCorrJetsPt_";
 257     histname += i;
 258     fHistCorrJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 259     fHistCorrJetsPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 260     fHistCorrJetsPt[i]->GetYaxis()->SetTitle("counts");
 261     fOutput->Add(fHistCorrJetsPt[i]);
 262
 263     histname = "fHistCorrLeadingJetPt_";
 264     histname += i;
 265     fHistCorrLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 266     fHistCorrLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 267     fHistCorrLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 268     fOutput->Add(fHistCorrLeadingJetPt[i]);
 269
 270     histname = "fHistRCPt_";
 271     histname += i;
 272     fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 273     fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 274     fHistRCPt[i]->GetYaxis()->SetTitle("counts");
 275     fOutput->Add(fHistRCPt[i]);
 276
 277     histname = "fHistRCPtExLJ_";
 278     histname += i;
 279     fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 280     fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 281     fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
 282     fOutput->Add(fHistRCPtExLJ[i]);
 283
 284     histname = "fHistRCPtRand_";
 285     histname += i;
 286     fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 287     fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 288     fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
 289     fOutput->Add(fHistRCPtRand[i]);
 290
 291     histname = "fHistDeltaPtRC_";
 292     histname += i;
 293     fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
 294     fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
 295     fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
 296     fOutput->Add(fHistDeltaPtRC[i]);
 297
 298     histname = "fHistDeltaPtRCExLJ_";
 299     histname += i;
 300     fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
 301     fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
 302     fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
 303     fOutput->Add(fHistDeltaPtRCExLJ[i]);
 304
 305     histname = "fHistDeltaPtRCRand_";
 306     histname += i;
 307     fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
 308     fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
 309     fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
 310     fOutput->Add(fHistDeltaPtRCRand[i]);
 311
 312     histname = "fHistEmbJets_";
 313     histname += i;
 314     fHistEmbJets[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 315     fHistEmbJets[i]->GetXaxis()->SetTitle("embedded jet p_{T} [GeV/c]");
 316     fHistEmbJets[i]->GetYaxis()->SetTitle("counts");
 317     fOutput->Add(fHistEmbJets[i]);
 318
 319     histname = "fHistEmbPart_";
 320     histname += i;
 321     fHistEmbPart[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 322     fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle p_{T} [GeV/c]");
 323     fHistEmbPart[i]->GetYaxis()->SetTitle("counts");
 324     fOutput->Add(fHistEmbPart[i]);
 325
 326     histname = "fHistDeltaPtEmb_";
 327     histname += i;
 328     fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
 329     fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T} [GeV/c]");
 330     fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
 331     fOutput->Add(fHistDeltaPtEmb[i]);
 332   }
 333
 334   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 335 }
 336
 337 //________________________________________________________________________
 338 Bool_t AliAnalysisTaskSAJF::RetrieveEventObjects()
 339 {
 340   // Retrieve event objects.
 341
 342   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 343     return kFALSE;
 344
 345   fRho = -1;
 346
 347   if (strcmp(fRhoName,"")) {
 348     TParameter<Double_t> *rho = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhoName));
 349
 350     if (rho) {
 351       fRho = rho->GetVal();
 352     }
 353     else {
 354       AliWarning(Form("Could not retrieve rho %s!", fRhoName.Data()));
 355     }
 356   }
 357
 358   if (strcmp(fEmbJetsName,"")) {
 359     fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
 360     if (!fEmbJets) {
 361       AliWarning(Form("Could not retrieve emb jets %s!", fEmbJetsName.Data()));
 362     }
 363   }
 364
 365   if (strcmp(fRandCaloName,"") && fAnaType == kEMCAL) {
 366     fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
 367     if (!fRandCaloClusters) {
 368       AliWarning(Form("Could not retrieve randomized clusters %s!", fRandCaloName.Data()));
 369     }
 370   }
 371
 372   if (strcmp(fRandTracksName,"")) {
 373     fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
 374     if (!fRandTracks) {
 375       AliWarning(Form("Could not retrieve randomized tracks %s!", fRandTracksName.Data()));
 376     }
 377   }
 378
 379   return kTRUE;
 380 }
 381
 382 //________________________________________________________________________
 383 Bool_t AliAnalysisTaskSAJF::FillHistograms()
 384 {
 385   // Fill histograms.
 386
 387   if (fRho < 0) {
 388     fHistRejectedEvents->Fill("Rho <= 0", 1);
 389     return kFALSE;
 390   }
 391
 392   Int_t maxJetIndex  = -1;
 393   Int_t max2JetIndex = -1;
 394
 395   // ************
 396   // General histograms
 397   // _________________________________
 398
 399   GetLeadingJets(maxJetIndex, max2JetIndex);
 400
 401   if (fSkipEventsWithNoJets && maxJetIndex < 0) {
 402     fHistRejectedEvents->Fill("No jets", 1);
 403     return kFALSE;
 404   }
 405
 406   AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(maxJetIndex));
 407   if (fSkipEventsWithNoJets && !jet) {
 408     fHistRejectedEvents->Fill("Max Jet not found", 1);
 409     return kFALSE;
 410   }
 411
 412   Float_t maxJetCorrPt = 0;
 413
 414   if (jet)
 415     maxJetCorrPt = jet->Pt() - fRho * jet->Area();
 416
 417   if (fSkipEventsWithNoJets && maxJetCorrPt <= 0) {
 418     fHistRejectedEvents->Fill("Max Jet <= 0", 1);
 419     return kFALSE;
 420   }
 421
 422   fHistCentrality->Fill(fCent);
 423   fHistRho[fCentBin]->Fill(fRho);
 424
 425   if (jet) {
 426     fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
 427     fHistRhoVSleadJetPt->Fill(fRho * jet->Area(), jet->Pt());
 428     fHistCorrLeadingJetPt[fCentBin]->Fill(maxJetCorrPt);
 429   }
 430
 431   AliEmcalJet* jet2 = 0;
 432   if (max2JetIndex >= 0)
 433     jet2 = dynamic_cast<AliEmcalJet*>(fJets->At(max2JetIndex));
 434
 435   if (jet2)
 436     fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());
 437
 438   DoJetLoop();
 439
 440   // ************
 441   // Random cones
 442   // _________________________________
 443
 444   const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
 445
 446   // Simple random cones
 447   Float_t RCpt = 0;
 448   Float_t RCeta = 0;
 449   Float_t RCphi = 0;
 450   GetRigidCone(RCpt, RCeta, RCphi, kFALSE, 0);
 451   if (RCpt > 0) {
 452     fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
 453     fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * fRho);
 454   }
 455
 456   // Random cones far from leading jet
 457   Float_t RCptExLJ = 0;
 458   Float_t RCetaExLJ = 0;
 459   Float_t RCphiExLJ = 0;
 460   GetRigidCone(RCptExLJ, RCetaExLJ, RCphiExLJ, kFALSE, jet);
 461   if (RCptExLJ > 0) {
 462     fHistRCPhiEta->Fill(RCetaExLJ, RCphiExLJ);
 463     fHistRhoVSRCPt->Fill(fRho, RCptExLJ / rcArea);
 464
 465     Float_t dphi = RCphiExLJ - jet->Phi();
 466     if (dphi > 4.8) dphi -= TMath::Pi() * 2;
 467     if (dphi < -1.6) dphi += TMath::Pi() * 2;
 468     fHistRCPtExLJVSDPhiLJ->Fill(RCptExLJ, dphi);
 469
 470     fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / rcArea);
 471     fHistDeltaPtRCExLJ[fCentBin]->Fill(RCptExLJ - rcArea * fRho);
 472   }
 473
 474   // Random cones with randomized particles
 475   Float_t RCptRand = 0;
 476   Float_t RCetaRand = 0;
 477   Float_t RCphiRand = 0;
 478   GetRigidCone(RCptRand, RCetaRand, RCphiRand, kTRUE, 0, fRandTracks, fRandCaloClusters);
 479   if (RCptRand > 0) {
 480     fHistRCPtRand[fCentBin]->Fill(RCptRand / rcArea);
 481     fHistDeltaPtRCRand[fCentBin]->Fill(RCptRand - rcArea * fRho);
 482   }
 483
 484   // ************
 485   // Embedding
 486   // _________________________________
 487
 488   if (!fEmbJets)
 489     return kTRUE;
 490
 491   AliEmcalJet *embJet  = 0;
 492   TObject     *maxPart = 0;
 493
 494   DoEmbJetLoop(embJet, maxPart);
 495
 496   if (embJet) {
 497
 498     fHistEmbJets[fCentBin]->Fill(embJet->Pt());
 499     fHistEmbPart[fCentBin]->Fill(embJet->MCPt());
 500     fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
 501     fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * fRho - embJet->MCPt());
 502     fHistRhoVSEmbBkg->Fill(embJet->Area() * fRho, embJet->Pt() - embJet->MCPt());
 503
 504     AliVCluster *cluster = dynamic_cast<AliVCluster*>(maxPart);
 505     if (cluster) {
 506       Float_t pos[3];
 507       cluster->GetPosition(pos);
 508       TVector3 clusVec(pos);
 509
 510       fHistEmbPartPhiEta->Fill(clusVec.Eta(), clusVec.Phi());
 511     }
 512     else {
 513       AliVTrack *track = dynamic_cast<AliVTrack*>(maxPart);
 514       if (track) {
 515         fHistEmbPartPhiEta->Fill(track->Eta(), track->Phi());
 516       }
 517       else {
 518         AliWarning(Form("%s - Embedded particle type not found or not recognized (neither AliVCluster nor AliVTrack)!", GetName()));
 519         return kTRUE;
 520       }
 521     }
 522
 523   }
 524   else {
 525     AliWarning(Form("%s - Embedded jet not found in the event!", GetName()));
 526   }
 527
 528   return kTRUE;
 529 }
 530
 531 //________________________________________________________________________
 532 void AliAnalysisTaskSAJF::GetLeadingJets(Int_t &maxJetIndex, Int_t &max2JetIndex)
 533 {
 534   // Get the leading jets.
 535
 536   if (!fJets)
 537     return;
 538
 539   const Int_t njets = fJets->GetEntriesFast();
 540
 541   Float_t maxJetPt = 0;
 542   Float_t max2JetPt = 0;
 543   for (Int_t ij = 0; ij < njets; ij++) {
 544
 545     AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(ij));
 546
 547     if (!jet) {
 548       AliError(Form("Could not receive jet %d", ij));
 549       continue;
 550     }
 551
 552     if (!AcceptJet(jet))
 553       continue;
 554
 555     Float_t corrPt = jet->Pt() - fRho * jet->Area();
 556
 557     if (maxJetIndex == -1 || maxJetPt < corrPt) {
 558       max2JetPt = maxJetPt;
 559       max2JetIndex = maxJetIndex;
 560       maxJetPt = corrPt;
 561       maxJetIndex = ij;
 562     }
 563     else if (max2JetIndex == -1 || max2JetPt < corrPt) {
 564       max2JetPt = corrPt;
 565       max2JetIndex = ij;
 566     }
 567   }
 568 }
 569
 570 //________________________________________________________________________
 571 void AliAnalysisTaskSAJF::DoJetLoop()
 572 {
 573   // Do the jet loop.
 574
 575   if (!fJets)
 576     return;
 577
 578   const Int_t njets = fJets->GetEntriesFast();
 579
 580   for (Int_t ij = 0; ij < njets; ij++) {
 581
 582     AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(ij));
 583
 584     if (!jet) {
 585       AliError(Form("Could not receive jet %d", ij));
 586       continue;
 587     }
 588
 589     if (!AcceptJet(jet))
 590       continue;
 591
 592     Float_t corrPt = jet->Pt() - fRho * jet->Area();
 593
 594     fHistJetsPt[fCentBin]->Fill(jet->Pt());
 595     fHistJetsPtArea[fCentBin]->Fill(corrPt, jet->Area());
 596     fHistCorrJetsPt[fCentBin]->Fill(corrPt);
 597
 598     fHistJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
 599
 600     if (fAnaType == kEMCAL)
 601       fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt());
 602
 603     if (fTracks) {
 604       for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 605         AliVParticle *track = jet->TrackAt(it, fTracks);
 606         if (track)
 607           fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt());
 608       }
 609     }
 610
 611     if (fCaloClusters) {
 612       for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 613         AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 614
 615         if (cluster) {
 616           TLorentzVector nPart;
 617           cluster->GetMomentum(nPart, fVertex);
 618         fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), jet->Pt());
 619         }
 620       }
 621     }
 622   } //jet loop
 623 }
 624
 625 //________________________________________________________________________
 626 void AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &maxPart)
 627 {
 628   // Do the embedded jet loop.
 629
 630   if (!fEmbJets)
 631     return;
 632
 633   TLorentzVector maxClusVect;
 634
 635   Int_t nembjets = fEmbJets->GetEntriesFast();
 636
 637   for (Int_t ij = 0; ij < nembjets; ij++) {
 638
 639     AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fEmbJets->At(ij));
 640
 641     if (!jet) {
 642       AliError(Form("Could not receive jet %d", ij));
 643       continue;
 644     }
 645
 646     if (!AcceptJet(jet))
 647       continue;
 648
 649     if (!jet->IsMC())
 650       continue;
 651
 652     AliVParticle *maxTrack = 0;
 653
 654     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 655       AliVParticle *track = jet->TrackAt(it, fTracks);
 656
 657       if (!track)
 658         continue;
 659
 660       if (track->GetLabel() != 100)
 661         continue;
 662
 663       if (!maxTrack || track->Pt() > maxTrack->Pt())
 664         maxTrack = track;
 665     }
 666
 667     AliVCluster *maxClus = 0;
 668
 669     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 670       AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 671
 672       if (!cluster)
 673         continue;
 674
 675       if (cluster->Chi2() != 100)
 676         continue;
 677
 678       TLorentzVector nPart;
 679       cluster->GetMomentum(nPart, fVertex);
 680
 681       if (!maxClus || nPart.Et() > maxClusVect.Et()) {
 682         new (&maxClusVect) TLorentzVector(nPart);
 683         maxClus = cluster;
 684       }
 685     }
 686
 687     if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
 688       maxPart = maxClus;
 689       embJet = jet;
 690     }
 691     else if (maxTrack) {
 692       maxPart = maxTrack;
 693       embJet = jet;
 694     }
 695
 696     return;  // MC jets found, exit
 697   }
 698 }
 699
 700 //________________________________________________________________________
 701 void AliAnalysisTaskSAJF::GetRigidCone(Float_t &pt, Float_t &eta, Float_t &phi, Bool_t acceptMC,
 702                                        AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
 703 {
 704   // Get rigid cone.
 705
 706   if (!tracks)
 707     tracks = fTracks;
 708
 709   if (!clusters)
 710     clusters = fCaloClusters;
 711
 712   if (!tracks && !clusters)
 713     return;
 714
 715   eta = 0;
 716   phi = 0;
 717   pt = 0;
 718
 719   Float_t LJeta = 999;
 720   Float_t LJphi = 999;
 721
 722   if (jet) {
 723     LJeta = jet->Eta();
 724     LJphi = jet->Phi();
 725   }
 726
 727   Float_t maxEta = fMaxEta;
 728   Float_t minEta = fMinEta;
 729   Float_t maxPhi = fMaxPhi;
 730   Float_t minPhi = fMinPhi;
 731
 732   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 733   if (minPhi < 0) minPhi = 0;
 734
 735   Float_t dLJ = 0;
 736   Int_t repeats = 0;
 737   do {
 738     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
 739     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
 740     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 741     repeats++;
 742   } while (dLJ < fMinRC2LJ && repeats < 999);
 743
 744   if (repeats == 999)
 745     return;
 746
 747   if (fAnaType == kEMCAL && clusters) {
 748     Int_t nclusters =  clusters->GetEntriesFast();
 749     for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 750       AliVCluster* cluster = dynamic_cast<AliVCluster*>(clusters->At(iClusters));
 751       if (!cluster) {
 752         AliError(Form("Could not receive cluster %d", iClusters));
 753         continue;
 754       }
 755
 756       if (!(cluster->IsEMCAL())) continue;
 757
 758       if (!AcceptCluster(cluster, acceptMC)) continue;
 759
 760       TLorentzVector nPart;
 761       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 762
 763       Float_t pos[3];
 764       cluster->GetPosition(pos);
 765       TVector3 clusVec(pos);
 766
 767       Float_t d = TMath::Sqrt((clusVec.Eta() - eta) * (clusVec.Eta() - eta) + (clusVec.Phi() - phi) * (clusVec.Phi() - phi));
 768
 769       if (d <= fJetRadius)
 770         pt += nPart.Pt();
 771     }
 772   }
 773
 774   if (tracks) {
 775     Int_t ntracks = tracks->GetEntriesFast();
 776     for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 777       AliVTrack* track = dynamic_cast<AliVTrack*>(tracks->At(iTracks));
 778       if(!track) {
 779         AliError(Form("Could not retrieve track %d",iTracks));
 780         continue;
 781       }
 782
 783       if (!AcceptTrack(track, acceptMC)) continue;
 784
 785       Float_t tracketa = track->Eta();
 786       Float_t trackphi = track->Phi();
 787
 788       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
 789         trackphi += 2 * TMath::Pi();
 790       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
 791         trackphi -= 2 * TMath::Pi();
 792
 793       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
 794
 795       if (d <= fJetRadius)
 796         pt += track->Pt();
 797     }
 798   }
 799 }
 800 //________________________________________________________________________
 801 void AliAnalysisTaskSAJF::Init()
 802 {
 803   // Initialize the analysis.
 804
 805   AliAnalysisTaskEmcalJet::Init();
 806
 807   const Float_t semiDiag = TMath::Sqrt((fMaxPhi - fMinPhi) * (fMaxPhi - fMinPhi) +
 808                                        (fMaxEta - fMinEta) * (fMaxEta - fMinEta)) / 2;
 809   if (fMinRC2LJ > semiDiag * 0.5) {
 810     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
 811                     "Will use fMinRC2LJ = %f", fMinRC2LJ, semiDiag * 0.5));
 812     fMinRC2LJ = semiDiag * 0.5;
 813   }
 814 }
 815
 816 //________________________________________________________________________
 817 void AliAnalysisTaskSAJF::Terminate(Option_t *)
 818 {
 819   // Called once at the end of the analysis.
 820 }