PWGGA/EMCALJetTasks/AliAnalysisTaskSAJF.cxx

   1 // $Id$
   2 //
   3 // Jet finder analysis task (S.Aiola).
   4 //
   5 //
   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 "AliVTrack.h"
  21 #include "AliEmcalJet.h"
  22 #include "AliVEventHandler.h"
  23 #include "AliLog.h"
  24
  25 #include "AliAnalysisTaskSAJF.h"
  26
  27 ClassImp(AliAnalysisTaskSAJF)
  28
  29 //________________________________________________________________________
  30 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
  31   AliAnalysisTaskSE("AliAnalysisTaskSAJF"),
  32   fAnaType(kTPC),
  33   fMinEta(-0.9),
  34   fMaxEta(0.9),
  35   fMinPhi(-10),
  36   fMaxPhi(10),
  37   fJetRadius(0.4),
  38   fTracksName("Tracks"),
  39   fCaloName("CaloClusters"),
  40   fJetsName("Jets"),
  41   fEmbJetsName("EmbJets"),
  42   fRhoName("Rho"),
  43   fNbins(500),
  44   fMinPt(0),
  45   fMaxPt(250),
  46   fPtCut(0.15),
  47   fPtCutJetPart(10),
  48   fTracks(0),
  49   fCaloClusters(0),
  50   fJets(0),
  51   fEmbJets(0),
  52   fCent(0),
  53   fCentBin(-1),
  54   fRho(-1),
  55   fOutput(0),
  56   fHistCentrality(0),
  57   fHistJetPhiEta(0),
  58   fHistEmbJetPhiEta(0),
  59   fHistEmbPartPhiEta(0),
  60   fHistRhoPartVSleadJetPt(0),
  61   fHistMedKtVSRhoPart(0),
  62   fHistRCPtVSRhoPart(0),
  63   fHistMedKtVSEmbBkg(0),
  64   fHistMedKtVSRCPt(0),
  65   fHistRCPtExLJVSDPhiLJ(0),
  66   fHistRCPhiEta(0)
  67
  68 {
  69   // Default constructor.
  70
  71   for (Int_t i = 0; i < 4; i++) {
  72     fHistJetsPt[i] = 0;
  73     fHistCorrJetsPt[i] = 0;
  74     fHistUnfoldedJetsPt[i] = 0;
  75     fHistJetsPtNonBias[i] = 0;
  76     fHistCorrJetsPtNonBias[i] = 0;
  77     fHistJetsNEFvsPt[i] = 0;
  78     fHistJetsZvsPt[i] = 0;
  79     fHistLeadingJetPt[i] = 0;
  80     fHistCorrLeadingJetPt[i] = 0;
  81     fHist2LeadingJetPt[i] = 0;
  82     fHistTracksPtLJ[i] = 0;
  83     fHistClusEtLJ[i] = 0;
  84     fHistTracksPtBkg[i] = 0;
  85     fHistClusEtBkg[i] = 0;
  86     fHistBkgClusMeanRho[i] = 0;
  87     fHistBkgTracksMeanRho[i] = 0;
  88     fHistMedianPtKtJet[i] = 0;
  89     fHistDeltaPtRC[i] = 0;
  90     fHistDeltaPtRCExLJ[i] = 0;
  91     fHistRCPt[i] = 0;
  92     fHistRCPtExLJ[i] = 0;
  93     fHistEmbJets[i] = 0;
  94     fHistEmbPart[i] = 0;
  95     fHistDeltaPtMedKtEmb[i] = 0;
  96   }
  97
  98   // Output slot #1 writes into a TH1 container
  99   DefineOutput(1, TList::Class());
 100 }
 101
 102 //________________________________________________________________________
 103 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
 104   AliAnalysisTaskSE(name),
 105   fAnaType(kTPC),
 106   fMinEta(-0.9),
 107   fMaxEta(0.9),
 108   fMinPhi(-10),
 109   fMaxPhi(10),
 110   fJetRadius(0.4),
 111   fTracksName("Tracks"),
 112   fCaloName("CaloClusters"),
 113   fJetsName("Jets"),
 114   fEmbJetsName("EmbJets"),
 115   fRhoName("Rho"),
 116   fNbins(500),
 117   fMinPt(0),
 118   fMaxPt(250),
 119   fPtCut(0.15),
 120   fPtCutJetPart(10),
 121   fTracks(0),
 122   fCaloClusters(0),
 123   fJets(0),
 124   fEmbJets(0),
 125   fCent(0),
 126   fCentBin(-1),
 127   fRho(-1),
 128   fOutput(0),
 129   fHistCentrality(0),
 130   fHistJetPhiEta(0),
 131   fHistEmbJetPhiEta(0),
 132   fHistEmbPartPhiEta(0),
 133   fHistRhoPartVSleadJetPt(0),
 134   fHistMedKtVSRhoPart(0),
 135   fHistRCPtVSRhoPart(0),
 136   fHistMedKtVSEmbBkg(0),
 137   fHistMedKtVSRCPt(0),
 138   fHistRCPtExLJVSDPhiLJ(0),
 139   fHistRCPhiEta(0)
 140 {
 141   // Standard constructor.
 142
 143   for (Int_t i = 0; i < 4; i++) {
 144     fHistJetsPt[i] = 0;
 145     fHistCorrJetsPt[i] = 0;
 146     fHistUnfoldedJetsPt[i] = 0;
 147     fHistJetsPtNonBias[i] = 0;
 148     fHistCorrJetsPtNonBias[i] = 0;
 149     fHistJetsNEFvsPt[i] = 0;
 150     fHistJetsZvsPt[i] = 0;
 151     fHistLeadingJetPt[i] = 0;
 152     fHistCorrLeadingJetPt[i] = 0;
 153     fHist2LeadingJetPt[i] = 0;
 154     fHistTracksPtLJ[i] = 0;
 155     fHistClusEtLJ[i] = 0;
 156     fHistTracksPtBkg[i] = 0;
 157     fHistClusEtBkg[i] = 0;
 158     fHistBkgClusMeanRho[i] = 0;
 159     fHistBkgTracksMeanRho[i] = 0;
 160     fHistMedianPtKtJet[i] = 0;
 161     fHistDeltaPtRC[i] = 0;
 162     fHistDeltaPtRCExLJ[i] = 0;
 163     fHistRCPt[i] = 0;
 164     fHistRCPtExLJ[i] = 0;
 165     fHistEmbJets[i] = 0;
 166     fHistEmbPart[i] = 0;
 167     fHistDeltaPtMedKtEmb[i] = 0;
 168   }
 169
 170   // Output slot #1 writes into a TH1 container
 171   DefineOutput(1, TList::Class());
 172 }
 173
 174 //________________________________________________________________________
 175 AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
 176 {
 177   // Destructor
 178 }
 179
 180 //________________________________________________________________________
 181 void AliAnalysisTaskSAJF::UserCreateOutputObjects()
 182 {
 183   // Create histograms
 184
 185   Float_t binWidth = (fMaxPt - fMinPt) / fNbins;
 186
 187   OpenFile(1);
 188   fOutput = new TList();
 189   fOutput->SetOwner();
 190
 191   fHistCentrality = new TH1F("fHistCentrality","Event centrality distribution", fNbins, 0, 100);
 192   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
 193   fHistCentrality->GetYaxis()->SetTitle("counts");
 194   fOutput->Add(fHistCentrality);
 195
 196   fHistJetPhiEta = new TH2F("fHistJetPhiEta","Phi-Eta distribution of jets", 20, -2, 2, 32, 0, 6.4);
 197   fHistJetPhiEta->GetXaxis()->SetTitle("Eta");
 198   fHistJetPhiEta->GetYaxis()->SetTitle("Phi");
 199   fOutput->Add(fHistJetPhiEta);
 200
 201   fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 20, -2, 2, 32, 0, 6.4);
 202   fHistEmbJetPhiEta->GetXaxis()->SetTitle("Eta");
 203   fHistEmbJetPhiEta->GetYaxis()->SetTitle("Phi");
 204   fOutput->Add(fHistEmbJetPhiEta);
 205
 206   fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 20, -2, 2, 32, 0, 6.4);
 207   fHistEmbPartPhiEta->GetXaxis()->SetTitle("Eta");
 208   fHistEmbPartPhiEta->GetYaxis()->SetTitle("Phi");
 209   fOutput->Add(fHistEmbPartPhiEta);
 210
 211   fHistRhoPartVSleadJetPt = new TH2F("fHistRhoPartVSleadJetPt","fHistRhoPartVSleadJetPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
 212   fHistRhoPartVSleadJetPt->GetXaxis()->SetTitle("#rho [GeV/c]");
 213   fHistRhoPartVSleadJetPt->GetYaxis()->SetTitle("Leading jet P_{T} [GeV/c]");
 214   fOutput->Add(fHistRhoPartVSleadJetPt);
 215
 216   fHistMedKtVSRhoPart = new TH2F("fHistMedKtVSRhoPart","fHistMedKtVSRhoPart", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
 217   fHistMedKtVSRhoPart->GetXaxis()->SetTitle("median kt P_{T} [GeV/c]");
 218   fHistMedKtVSRhoPart->GetYaxis()->SetTitle("#rho [GeV/c]");
 219   fOutput->Add(fHistMedKtVSRhoPart);
 220
 221   fHistRCPtVSRhoPart = new TH2F("fHistRCPtVSRhoPart","fHistRCPtVSRhoPart", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
 222   fHistRCPtVSRhoPart->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
 223   fHistRCPtVSRhoPart->GetYaxis()->SetTitle("#rho [GeV/c]");
 224   fOutput->Add(fHistRCPtVSRhoPart);
 225
 226   fHistMedKtVSEmbBkg = new TH2F("fHistMedKtVSEmbBkg","fHistMedKtVSEmbBkg", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
 227   fHistMedKtVSEmbBkg->GetXaxis()->SetTitle("median kt P_{T} [GeV/c]");
 228   fHistMedKtVSEmbBkg->GetYaxis()->SetTitle("background of embedded track P_{T} [GeV]");
 229   fOutput->Add(fHistMedKtVSEmbBkg);
 230
 231   fHistMedKtVSRCPt = new TH2F("fHistMedKtVSRCPt","fHistMedKtVSRCPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
 232   fHistMedKtVSRCPt->GetXaxis()->SetTitle("median kt P_{T} [GeV/c]");
 233   fHistMedKtVSRCPt->GetYaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
 234   fOutput->Add(fHistMedKtVSRCPt);
 235
 236   fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinPt, fMaxPt, 128, -1.6, 4.8);
 237   fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
 238   fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
 239   fOutput->Add(fHistRCPtExLJVSDPhiLJ);
 240
 241   fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 20, -2, 2, 32, 0, 6.4);
 242   fHistRCPhiEta->GetXaxis()->SetTitle("Eta");
 243   fHistRCPhiEta->GetYaxis()->SetTitle("Phi");
 244   fOutput->Add(fHistRCPhiEta);
 245
 246   TString histname;
 247
 248   for (Int_t i = 0; i < 4; i++) {
 249     histname = "fHistJetsPt_";
 250     histname += i;
 251     fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 252     fHistJetsPt[i]->GetXaxis()->SetTitle("Pt [GeV/c]");
 253     fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
 254     fOutput->Add(fHistJetsPt[i]);
 255
 256     histname = "fHistCorrJetsPt_";
 257     histname += i;
 258     fHistCorrJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
 259     fHistCorrJetsPt[i]->GetXaxis()->SetTitle("Pt [GeV/c]");
 260     fHistCorrJetsPt[i]->GetYaxis()->SetTitle("counts");
 261     fOutput->Add(fHistCorrJetsPt[i]);
 262
 263     histname = "fHistUnfoldedJetsPt_";
 264     histname += i;
 265     fHistUnfoldedJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
 266     fHistUnfoldedJetsPt[i]->GetXaxis()->SetTitle("Pt [GeV/c]");
 267     fHistUnfoldedJetsPt[i]->GetYaxis()->SetTitle("counts");
 268     fOutput->Add(fHistUnfoldedJetsPt[i]);
 269
 270     histname = "fHistJetsPtNonBias_";
 271     histname += i;
 272     fHistJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 273     fHistJetsPtNonBias[i]->GetXaxis()->SetTitle("Pt [GeV/c]");
 274     fHistJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
 275     fOutput->Add(fHistJetsPtNonBias[i]);
 276
 277     histname = "fHistCorrJetsPtNonBias_";
 278     histname += i;
 279     fHistCorrJetsPtNonBias[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
 280     fHistCorrJetsPtNonBias[i]->GetXaxis()->SetTitle("Pt [GeV/c]");
 281     fHistCorrJetsPtNonBias[i]->GetYaxis()->SetTitle("counts");
 282     fOutput->Add(fHistCorrJetsPtNonBias[i]);
 283
 284     histname = "fHistJetsNEFvsPt_";
 285     histname += i;
 286     fHistJetsNEFvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins * 1.5, -fMaxPt * 0.75, fMaxPt * 0.75);
 287     fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
 288     fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("Pt [GeV/c]");
 289     fOutput->Add(fHistJetsNEFvsPt[i]);
 290
 291     histname = "fHistJetsZvsPt_";
 292     histname += i;
 293     fHistJetsZvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins * 1.5, -fMaxPt * 0.75, fMaxPt * 0.75);
 294     fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
 295     fHistJetsZvsPt[i]->GetYaxis()->SetTitle("Pt [GeV/c]");
 296     fOutput->Add(fHistJetsZvsPt[i]);
 297
 298     histname = "fHistLeadingJetPt_";
 299     histname += i;
 300     fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 301     fHistLeadingJetPt[i]->GetXaxis()->SetTitle("P_{T} [GeV]");
 302     fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 303     fOutput->Add(fHistLeadingJetPt[i]);
 304
 305     histname = "fHist2LeadingJetPt_";
 306     histname += i;
 307     fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 308     fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("P_{T} [GeV]");
 309     fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 310     fOutput->Add(fHist2LeadingJetPt[i]);
 311
 312     histname = "fHistCorrLeadingJetPt_";
 313     histname += i;
 314     fHistCorrLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxPt, fMaxPt);
 315     fHistCorrLeadingJetPt[i]->GetXaxis()->SetTitle("P_{T} [GeV]");
 316     fHistCorrLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 317     fOutput->Add(fHistCorrLeadingJetPt[i]);
 318
 319     histname = "fHistTracksPtLJ_";
 320     histname += i;
 321     fHistTracksPtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
 322     fHistTracksPtLJ[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 323     fHistTracksPtLJ[i]->GetYaxis()->SetTitle("counts");
 324     fOutput->Add(fHistTracksPtLJ[i]);
 325
 326     histname = "fHistClusEtLJ_";
 327     histname += i;
 328     fHistClusEtLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
 329     fHistClusEtLJ[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
 330     fHistClusEtLJ[i]->GetYaxis()->SetTitle("counts");
 331     fOutput->Add(fHistClusEtLJ[i]);
 332
 333     histname = "fHistTracksPtBkg_";
 334     histname += i;
 335     fHistTracksPtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
 336     fHistTracksPtBkg[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 337     fHistTracksPtBkg[i]->GetYaxis()->SetTitle("counts");
 338     fOutput->Add(fHistTracksPtBkg[i]);
 339
 340     histname = "fHistClusEtBkg_";
 341     histname += i;
 342     fHistClusEtBkg[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt / 5);
 343     fHistClusEtBkg[i]->GetXaxis()->SetTitle("E_{T} [GeV]");
 344     fHistClusEtBkg[i]->GetYaxis()->SetTitle("counts");
 345     fOutput->Add(fHistClusEtBkg[i]);
 346
 347     histname = "fHistBkgClusMeanRho_";
 348     histname += i;
 349     fHistBkgClusMeanRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 350     fHistBkgClusMeanRho[i]->GetXaxis()->SetTitle("#rho [GeV]");
 351     fHistBkgClusMeanRho[i]->GetYaxis()->SetTitle("counts");
 352     fOutput->Add(fHistBkgClusMeanRho[i]);
 353
 354     histname = "fHistBkgTracksMeanRho_";
 355     histname += i;
 356     fHistBkgTracksMeanRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 357     fHistBkgTracksMeanRho[i]->GetXaxis()->SetTitle("#rho [GeV/c]");
 358     fHistBkgTracksMeanRho[i]->GetYaxis()->SetTitle("counts");
 359     fOutput->Add(fHistBkgTracksMeanRho[i]);
 360
 361     histname = "fHistMedianPtKtJet_";
 362     histname += i;
 363     fHistMedianPtKtJet[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 364     fHistMedianPtKtJet[i]->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 365     fHistMedianPtKtJet[i]->GetYaxis()->SetTitle("counts");
 366     fOutput->Add(fHistMedianPtKtJet[i]);
 367
 368     histname = "fHistDeltaPtRC_";
 369     histname += i;
 370     fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
 371     fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
 372     fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
 373     fOutput->Add(fHistDeltaPtRC[i]);
 374
 375     histname = "fHistDeltaPtRCExLJ_";
 376     histname += i;
 377     fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
 378     fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
 379     fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
 380     fOutput->Add(fHistDeltaPtRCExLJ[i]);
 381
 382     histname = "fHistRCPt_";
 383     histname += i;
 384     fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 385     fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
 386     fHistRCPt[i]->GetYaxis()->SetTitle("counts");
 387     fOutput->Add(fHistRCPt[i]);
 388
 389     histname = "fHistRCPtExLJ_";
 390     histname += i;
 391     fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 392     fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone P_{T} [GeV/c]");
 393     fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
 394     fOutput->Add(fHistRCPtExLJ[i]);
 395
 396     histname = "fHistEmbJets_";
 397     histname += i;
 398     fHistEmbJets[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 399     fHistEmbJets[i]->GetXaxis()->SetTitle("embedded jet P_{T} [GeV/c]");
 400     fHistEmbJets[i]->GetYaxis()->SetTitle("counts");
 401     fOutput->Add(fHistEmbJets[i]);
 402
 403     histname = "fHistEmbPart_";
 404     histname += i;
 405     fHistEmbPart[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt, fMaxPt);
 406     fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle P_{T} [GeV/c]");
 407     fHistEmbPart[i]->GetYaxis()->SetTitle("counts");
 408     fOutput->Add(fHistEmbPart[i]);
 409
 410     histname = "fHistDeltaPtMedKtEmb_";
 411     histname += i;
 412     fHistDeltaPtMedKtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinPt - fMaxPt / 2 + binWidth / 2, fMinPt + fMaxPt / 2 + binWidth / 2);
 413     fHistDeltaPtMedKtEmb[i]->GetXaxis()->SetTitle("#deltaP_{T} [GeV/c]");
 414     fHistDeltaPtMedKtEmb[i]->GetYaxis()->SetTitle("counts");
 415     fOutput->Add(fHistDeltaPtMedKtEmb[i]);
 416   }
 417
 418   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 419 }
 420
 421 //________________________________________________________________________
 422 void AliAnalysisTaskSAJF::RetrieveEventObjects()
 423 {
 424   if (strcmp(fCaloName,"") && fAnaType == kEMCAL) {
 425     fCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
 426     if (!fCaloClusters) {
 427       AliWarning(Form("Could not retrieve clusters!"));
 428     }
 429   }
 430
 431   fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
 432   if (!fTracks) {
 433     AliWarning(Form("Could not retrieve tracks!"));
 434   }
 435
 436   fJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetsName));
 437   if (!fJets) {
 438     AliWarning(Form("Could not retrieve jets!"));
 439   }
 440
 441   if (strcmp(fEmbJetsName,"")) {
 442     fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
 443     if (!fEmbJets) {
 444       AliWarning(Form("Could not retrieve emb jets!"));
 445     }
 446   }
 447
 448   AliCentrality *aliCent = InputEvent()->GetCentrality();
 449   if (aliCent) {
 450     fCent = aliCent->GetCentralityPercentile("V0M");
 451     if      (fCent >=  0 && fCent <   10) fCentBin = 0;
 452     else if (fCent >= 10 && fCent <   30) fCentBin = 1;
 453     else if (fCent >= 30 && fCent <   50) fCentBin = 2;
 454     else if (fCent >= 50 && fCent <= 100) fCentBin = 3;
 455     else {
 456       AliWarning(Form("Negative centrality: %f. Assuming 99", fCent));
 457       fCentBin = 3;
 458     }
 459   }
 460   else {
 461     AliWarning(Form("Could not retrieve centrality information! Assuming 99"));
 462     fCentBin = 3;
 463   }
 464
 465   fRho = -1;
 466
 467   if (strcmp(fRhoName,"")) {
 468     TParameter<Double_t> *rho = dynamic_cast<TParameter<Double_t>*>(InputEvent()->FindListObject(fRhoName));
 469
 470     if (rho) {
 471       fRho = rho->GetVal();
 472     }
 473     else {
 474       AliWarning("Could not retrieve rho information.");
 475     }
 476   }
 477 }
 478
 479 //________________________________________________________________________
 480 AliVTrack* AliAnalysisTaskSAJF::GetTrack(const Int_t i) const
 481 {
 482   if (fTracks)
 483     return dynamic_cast<AliVTrack*>(fTracks->At(i));
 484   else
 485     return 0;
 486 }
 487
 488 //________________________________________________________________________
 489 Int_t AliAnalysisTaskSAJF::GetNumberOfTracks() const
 490 {
 491   if (fTracks)
 492     return fTracks->GetEntriesFast();
 493   else
 494     return 0;
 495 }
 496
 497 //________________________________________________________________________
 498 AliVCluster* AliAnalysisTaskSAJF::GetCaloCluster(const Int_t i) const
 499 {
 500   if (fCaloClusters)
 501     return dynamic_cast<AliVCluster*>(fCaloClusters->At(i));
 502   else
 503     return 0;
 504 }
 505
 506 //________________________________________________________________________
 507 Int_t AliAnalysisTaskSAJF::GetNumberOfCaloClusters() const
 508 {
 509   if (fCaloClusters)
 510     return fCaloClusters->GetEntriesFast();
 511   else
 512     return 0;
 513 }
 514
 515 //________________________________________________________________________
 516 AliEmcalJet* AliAnalysisTaskSAJF::GetJet(const Int_t i) const
 517 {
 518   if (fJets)
 519     return dynamic_cast<AliEmcalJet*>(fJets->At(i));
 520   else
 521     return 0;
 522 }
 523
 524 //________________________________________________________________________
 525 Int_t AliAnalysisTaskSAJF::GetNumberOfJets() const
 526 {
 527   if (fJets)
 528     return fJets->GetEntriesFast();
 529   else
 530     return 0;
 531 }
 532
 533 //________________________________________________________________________
 534 AliEmcalJet* AliAnalysisTaskSAJF::GetEmbJet(const Int_t i) const
 535 {
 536   if (fEmbJets)
 537     return dynamic_cast<AliEmcalJet*>(fEmbJets->At(i));
 538   else
 539     return 0;
 540 }
 541
 542 //________________________________________________________________________
 543 Int_t AliAnalysisTaskSAJF::GetNumberOfEmbJets() const
 544 {
 545   if (fEmbJets)
 546     return fEmbJets->GetEntriesFast();
 547   else
 548     return 0;
 549 }
 550
 551 //________________________________________________________________________
 552 void AliAnalysisTaskSAJF::FillHistograms()
 553 {
 554   Float_t A = fJetRadius * fJetRadius * TMath::Pi();
 555
 556   Int_t maxJetIndex  = -1;
 557   Int_t max2JetIndex = -1;
 558
 559   DoJetLoop(maxJetIndex, max2JetIndex);
 560
 561   if (maxJetIndex < 0)
 562     return;
 563
 564   AliEmcalJet* jet = GetJet(maxJetIndex);
 565   if (!jet)
 566     return;
 567
 568   fHistCentrality->Fill(fCent);
 569
 570   fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
 571   jet->SortConstituents();
 572
 573   AliEmcalJet* jet2 = 0;
 574   if (max2JetIndex >= 0)
 575     jet2 = GetJet(max2JetIndex);
 576
 577   if (jet2) {
 578     fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());
 579     jet2->SortConstituents();
 580   }
 581
 582   fHistMedianPtKtJet[fCentBin]->Fill(fRho);
 583
 584   Float_t maxJetCorrPt = jet->Pt() - fRho * jet->Area();
 585   if (maxJetCorrPt > 0)
 586     fHistCorrLeadingJetPt[fCentBin]->Fill(maxJetCorrPt);
 587
 588   Float_t rhoTracksLJex = 0;
 589   Float_t rhoTracks = 0;
 590   DoTrackLoop(rhoTracksLJex, rhoTracks, maxJetIndex, max2JetIndex);
 591   fHistBkgTracksMeanRho[fCentBin]->Fill(rhoTracksLJex);
 592
 593   Float_t rhoClusLJex = 0;
 594   Float_t rhoClus = 0;
 595   if (fAnaType == kEMCAL)
 596     DoClusterLoop(rhoClusLJex, rhoClus, maxJetIndex, max2JetIndex);
 597
 598   Float_t rhoPartLJex = rhoTracksLJex + rhoClusLJex;
 599
 600   fHistBkgClusMeanRho[fCentBin]->Fill(rhoClusLJex);
 601
 602   fHistRhoPartVSleadJetPt->Fill(jet->Area() * rhoPartLJex, jet->Pt());
 603
 604   fHistMedKtVSRhoPart->Fill(fRho, rhoPartLJex);
 605
 606   Int_t nRCs = 1; //(Int_t)(GetArea() / A - 1);
 607
 608   for (Int_t i = 0; i < nRCs; i++) {
 609     Float_t RCpt = 0;
 610     Float_t RCeta = 0;
 611     Float_t RCphi = 0;
 612     GetRigidConePt(RCpt, RCeta, RCphi, 0);
 613
 614     Float_t RCptExLJ = 0;
 615     Float_t RCetaExLJ = 0;
 616     Float_t RCphiExLJ = 0;
 617     GetRigidConePt(RCptExLJ, RCetaExLJ, RCphiExLJ, jet);
 618
 619     fHistDeltaPtRC[fCentBin]->Fill(RCpt - A * fRho);
 620     fHistDeltaPtRCExLJ[fCentBin]->Fill(RCptExLJ - A * fRho);
 621
 622     fHistRCPt[fCentBin]->Fill(RCpt / A);
 623     fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / A);
 624     fHistRCPtVSRhoPart->Fill(RCptExLJ / A, rhoPartLJex);
 625
 626     fHistMedKtVSRCPt->Fill(A * fRho, RCptExLJ);
 627
 628     fHistRCPhiEta->Fill(RCeta, RCphiExLJ);
 629
 630     Float_t dphi = RCphiExLJ - jet->Phi();
 631     if (dphi > 4.8) dphi -= TMath::Pi() * 2;
 632     if (dphi < -1.6) dphi += TMath::Pi() * 2;
 633     fHistRCPtExLJVSDPhiLJ->Fill(RCptExLJ, dphi);
 634   }
 635
 636   AliEmcalJet *maxJet  = 0;
 637   TObject     *maxPart = 0;
 638
 639   Bool_t embOK = DoEmbJetLoop(maxJet, maxPart);
 640
 641   if (embOK) {
 642     AliVCluster *cluster = dynamic_cast<AliVCluster*>(maxPart);
 643     AliVTrack *track = dynamic_cast<AliVTrack*>(maxPart);
 644     Float_t maxEmbPartPt = 0;
 645     Float_t maxEmbPartEta = 0;
 646     Float_t maxEmbPartPhi = 0;
 647     if (cluster) {
 648       Double_t vertex[3] = {0, 0, 0};
 649       InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 650       TLorentzVector nPart;
 651       cluster->GetMomentum(nPart, vertex);
 652       Float_t pos[3];
 653       cluster->GetPosition(pos);
 654       TVector3 clusVec(pos);
 655
 656       maxEmbPartPt = nPart.Et();
 657       maxEmbPartEta = clusVec.Eta();
 658       maxEmbPartPhi = clusVec.Phi();
 659     }
 660     else if (track) {
 661       maxEmbPartPt = track->Pt();
 662       maxEmbPartEta = track->Eta();
 663       maxEmbPartPhi = track->Phi();
 664     }
 665     else {
 666       AliWarning(Form("%s - Embedded particle type not recognized (neither AliVCluster nor AliVTrack)!", GetName()));
 667       return;
 668     }
 669     fHistEmbJets[fCentBin]->Fill(maxJet->Pt());
 670     fHistEmbPart[fCentBin]->Fill(maxEmbPartPt);
 671     fHistDeltaPtMedKtEmb[fCentBin]->Fill(maxJet->Pt() - maxJet->Area() * fRho - maxEmbPartPt);
 672     fHistMedKtVSEmbBkg->Fill(maxJet->Area() * fRho, maxJet->Pt() - maxEmbPartPt);
 673     fHistEmbJetPhiEta->Fill(maxJet->Eta(), maxJet->Phi());
 674     fHistEmbPartPhiEta->Fill(maxEmbPartEta, maxEmbPartPhi);
 675   }
 676   else {
 677     if (maxPart != 0)
 678       AliWarning(Form("%s - Embedded particle is not the leading particle of the leading jet!", GetName()));
 679   }
 680 }
 681
 682 //________________________________________________________________________
 683 void AliAnalysisTaskSAJF::DoJetLoop(Int_t &maxJetIndex, Int_t &max2JetIndex)
 684 {
 685   Double_t vertex[3] = {0, 0, 0};
 686   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 687
 688   Int_t njets = GetNumberOfJets();
 689
 690   Float_t maxJetPt = 0;
 691   Float_t max2JetPt = 0;
 692   for (Int_t ij = 0; ij < njets; ij++) {
 693
 694     AliEmcalJet* jet = GetJet(ij);
 695
 696     if (!jet) {
 697       AliError(Form("Could not receive jet %d", ij));
 698       continue;
 699     }
 700
 701     if (jet->Pt() <= 0)
 702       continue;
 703
 704     if (!AcceptJet(jet))
 705       continue;
 706
 707     Float_t corrPt = jet->Pt() - fRho * jet->Area();
 708
 709     fHistJetsPtNonBias[fCentBin]->Fill(jet->Pt());
 710     fHistCorrJetsPtNonBias[fCentBin]->Fill(corrPt);
 711
 712     if (!AcceptJet(jet, kTRUE))
 713       continue;
 714
 715     fHistJetsPt[fCentBin]->Fill(jet->Pt());
 716     fHistCorrJetsPt[fCentBin]->Fill(corrPt);
 717
 718     fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
 719     fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), corrPt);
 720
 721     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 722       AliVTrack *track = jet->TrackAt(it, fTracks);
 723       if (track)
 724         fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), corrPt);
 725     }
 726
 727     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 728       AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 729
 730       if (cluster) {
 731         TLorentzVector nPart;
 732         cluster->GetMomentum(nPart, vertex);
 733         fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), corrPt);
 734       }
 735     }
 736
 737     if (maxJetPt < corrPt) {
 738       max2JetPt = maxJetPt;
 739       max2JetIndex = maxJetIndex;
 740       maxJetPt = corrPt;
 741       maxJetIndex = ij;
 742     }
 743     else if (max2JetPt < corrPt) {
 744       max2JetPt = corrPt;
 745       max2JetIndex = ij;
 746     }
 747   } //jet loop
 748 }
 749
 750 //________________________________________________________________________
 751 Bool_t AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &maxJet, TObject* &maxPart)
 752 {
 753   Double_t vertex[3] = {0, 0, 0};
 754   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 755
 756   Int_t nembjets = GetNumberOfEmbJets();
 757
 758   Int_t maxJetIdx = -1;
 759   Int_t maxTrackIdx = -1;
 760   Int_t maxClusIdx = -1;
 761
 762   Float_t maxTrackPt = 0;
 763   Float_t maxClusEt = 0;
 764   Float_t maxJetPt = 0;
 765
 766   for (Int_t ij = 0; ij < nembjets; ij++) {
 767
 768     AliEmcalJet* jet = GetEmbJet(ij);
 769
 770     if (!jet) {
 771       AliError(Form("Could not receive jet %d", ij));
 772       continue;
 773     }
 774
 775     if (jet->Pt() <= 0)
 776         continue;
 777
 778     if (!AcceptJet(jet, kTRUE))
 779       continue;
 780
 781     if (jet->Pt() > maxJetPt) {
 782       maxJetPt = jet->Pt();
 783       maxJetIdx = ij;
 784     }
 785   }
 786
 787   if (maxJetPt <= 0)
 788     return kFALSE;
 789
 790   maxJet = GetEmbJet(maxJetIdx);
 791
 792   for (Int_t it = 0; it < maxJet->GetNumberOfTracks(); it++) {
 793     AliVTrack *track = maxJet->TrackAt(it, fTracks);
 794
 795     if (!track) continue;
 796
 797     if (track->Pt() > maxTrackPt) {
 798       maxTrackPt = track->Pt();
 799       maxTrackIdx = it;
 800     }
 801   }
 802
 803   for (Int_t ic = 0; ic < maxJet->GetNumberOfClusters(); ic++) {
 804     AliVCluster *cluster = maxJet->ClusterAt(ic, fCaloClusters);
 805
 806     if (!cluster) continue;
 807
 808     TLorentzVector nPart;
 809     cluster->GetMomentum(nPart, vertex);
 810
 811     if (nPart.Et() > maxClusEt) {
 812       maxClusEt = nPart.Et();
 813       maxClusIdx = ic;
 814     }
 815   }
 816
 817   if (maxClusEt > maxTrackPt) {
 818     AliVCluster *cluster = maxJet->ClusterAt(maxClusIdx, fCaloClusters);
 819     maxPart = cluster;
 820     return (Bool_t)(cluster->Chi2() == 100);
 821   }
 822   else {
 823     AliVTrack *track = maxJet->TrackAt(maxTrackIdx, fTracks);
 824     maxPart = track;
 825     return (Bool_t)(track->GetLabel() == 100);
 826   }
 827 }
 828
 829 //________________________________________________________________________
 830 void AliAnalysisTaskSAJF::DoTrackLoop(Float_t &rhoTracksLJex, Float_t &rhoTracks, Int_t maxJetIndex, Int_t max2JetIndex)
 831 {
 832   AliEmcalJet* jet = 0;
 833   if (max2JetIndex >= 0)
 834     jet = GetJet(maxJetIndex);
 835
 836   AliEmcalJet* jet2 = 0;
 837   if (max2JetIndex >= 0)
 838     jet2 = GetJet(max2JetIndex);
 839
 840   Float_t area = GetArea();
 841   if (jet) area -= jet->Area();
 842   if (jet2) area -= jet2->Area();
 843
 844   Int_t ntracks = GetNumberOfTracks();
 845
 846   rhoTracksLJex = 0;
 847   rhoTracks = 0;
 848
 849   for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
 850     AliVTrack* track = GetTrack(iTracks);
 851     if(!track) {
 852       AliError(Form("Could not retrieve track %d",iTracks));
 853       continue;
 854     }
 855
 856     if (!AcceptTrack(track)) continue;
 857
 858     rhoTracks += track->Pt();
 859
 860     if (jet && IsJetTrack(jet, iTracks)) {
 861       fHistTracksPtLJ[fCentBin]->Fill(track->Pt());
 862     }
 863     else if (!jet2 || !IsJetTrack(jet2, iTracks)) {
 864       fHistTracksPtBkg[fCentBin]->Fill(track->Pt());
 865       rhoTracksLJex += track->Pt();
 866     }
 867   }
 868   rhoTracksLJex /= area;
 869   rhoTracks /= area;
 870 }
 871
 872 //________________________________________________________________________
 873 void AliAnalysisTaskSAJF::DoClusterLoop(Float_t &rhoClusLJex, Float_t &rhoClus, Int_t maxJetIndex, Int_t max2JetIndex)
 874 {
 875   Double_t vertex[3] = {0, 0, 0};
 876   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 877
 878   AliEmcalJet* jet = 0;
 879   if (max2JetIndex >= 0)
 880     jet = GetJet(maxJetIndex);
 881
 882   AliEmcalJet* jet2 = 0;
 883   if (max2JetIndex >= 0)
 884     jet2 = GetJet(max2JetIndex);
 885
 886   Float_t area = GetArea();
 887   if (jet) area -= jet->Area();
 888   if (jet2) area -= jet2->Area();
 889
 890   rhoClusLJex = 0;
 891   rhoClus = 0;
 892
 893   Int_t nclusters =  GetNumberOfCaloClusters();
 894   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 895     AliVCluster* cluster = GetCaloCluster(iClusters);
 896     if (!cluster) {
 897       AliError(Form("Could not receive cluster %d", iClusters));
 898       continue;
 899     }
 900
 901     if (!(cluster->IsEMCAL())) continue;
 902
 903     if (!AcceptCluster(cluster)) continue;
 904
 905     TLorentzVector nPart;
 906     cluster->GetMomentum(nPart, vertex);
 907
 908     rhoClus += nPart.Et();
 909
 910     if (jet && IsJetCluster(jet, iClusters)) {
 911       fHistClusEtLJ[fCentBin]->Fill(nPart.Et());
 912     }
 913     else if (!jet2 || !IsJetCluster(jet2, iClusters)) {
 914       fHistClusEtBkg[fCentBin]->Fill(nPart.Et());
 915       rhoClusLJex += nPart.Et();
 916     }
 917   } //cluster loop
 918   rhoClusLJex /= area;
 919   rhoClus /= area;
 920 }
 921
 922 //________________________________________________________________________
 923 void AliAnalysisTaskSAJF::Init()
 924 {
 925   if (fAnaType == kTPC) {
 926     fMinEta = -0.9;
 927     fMaxEta = 0.9;
 928     fMinPhi = -10;
 929     fMaxPhi = 10;
 930   }
 931   else if (fAnaType == kEMCAL) {
 932     fMinEta = -0.7;
 933     fMaxEta = 0.7;
 934     fMinPhi = 80 * TMath::DegToRad();
 935     fMaxPhi = 180 * TMath::DegToRad();
 936   }
 937   else {
 938     AliWarning("Analysis type not recognized! Assuming kTPC...");
 939     fAnaType = kTPC;
 940     Init();
 941   }
 942 }
 943
 944 //________________________________________________________________________
 945 Bool_t AliAnalysisTaskSAJF::GetRigidConePt(Float_t &pt, Float_t &eta, Float_t &phi, AliEmcalJet *jet, Float_t minD)
 946 {
 947   static TRandom3 random;
 948
 949   Double_t vertex[3] = {0, 0, 0};
 950   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 951
 952   eta = 0;
 953   phi = 0;
 954
 955   Float_t LJeta = 999;
 956   Float_t LJphi = 999;
 957
 958   if (jet) {
 959     LJeta = jet->Eta();
 960     LJphi = jet->Phi();
 961   }
 962
 963   Float_t maxEta = fMaxEta - fJetRadius;
 964   Float_t minEta = fMinEta + fJetRadius;
 965   Float_t maxPhi = fMaxPhi - fJetRadius;
 966   Float_t minPhi = fMinPhi + fJetRadius;
 967
 968   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 969   if (minPhi < 0) minPhi = 0;
 970
 971   Float_t dLJ = 0;
 972   do {
 973     eta = random.Rndm() * (maxEta - minEta) + minEta;
 974     phi = random.Rndm() * (maxPhi - minPhi) + minPhi;
 975     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 976   } while (dLJ < minD && !AcceptJet(eta, phi));
 977
 978   pt = 0;
 979
 980   if (fAnaType == kEMCAL) {
 981     Int_t nclusters =  GetNumberOfCaloClusters();
 982     for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 983       AliVCluster* cluster = GetCaloCluster(iClusters);
 984       if (!cluster) {
 985         AliError(Form("Could not receive cluster %d", iClusters));
 986         continue;
 987       }
 988
 989       if (!(cluster->IsEMCAL())) continue;
 990
 991       if (!AcceptCluster(cluster)) continue;
 992
 993       TLorentzVector nPart;
 994       cluster->GetMomentum(nPart, vertex);
 995
 996       Float_t pos[3];
 997       cluster->GetPosition(pos);
 998       TVector3 clusVec(pos);
 999
1000       Float_t d = TMath::Sqrt((clusVec.Eta() - eta) * (clusVec.Eta() - eta) + (clusVec.Phi() - phi) * (clusVec.Phi() - phi));
1001
1002       if (d <= fJetRadius)
1003         pt += nPart.Pt();
1004     }
1005   }
1006
1007   Int_t ntracks = GetNumberOfTracks();
1008   for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1009     AliVTrack* track = GetTrack(iTracks);
1010     if(!track) {
1011       AliError(Form("Could not retrieve track %d",iTracks));
1012       continue;
1013     }
1014
1015     if (!AcceptTrack(track)) continue;
1016
1017     Float_t tracketa = track->Eta();
1018     Float_t trackphi = track->Phi();
1019
1020     if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
1021       trackphi += 2 * TMath::Pi();
1022     if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
1023       trackphi -= 2 * TMath::Pi();
1024
1025     Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
1026
1027     if (d <= fJetRadius)
1028       pt += track->Pt();
1029   }
1030
1031   return kTRUE;
1032 }
1033
1034 //________________________________________________________________________
1035 Float_t AliAnalysisTaskSAJF::GetArea() const
1036 {
1037   Float_t dphi = fMaxPhi - fMinPhi;
1038   if (dphi > TMath::Pi() * 2) dphi = TMath::Pi() * 2;
1039   Float_t deta = fMaxEta - fMinEta;
1040   return deta * dphi;
1041 }
1042
1043 //________________________________________________________________________
1044 Bool_t AliAnalysisTaskSAJF::IsJetTrack(AliEmcalJet* jet, Int_t itrack, Bool_t sorted) const
1045 {
1046   for (Int_t i = 0; i < jet->GetNumberOfTracks(); i++) {
1047     Int_t ijettrack = jet->TrackAt(i);
1048     if (sorted && ijettrack > itrack)
1049       return kFALSE;
1050     if (ijettrack == itrack)
1051       return kTRUE;
1052   }
1053   return kFALSE;
1054 }
1055
1056 //________________________________________________________________________
1057 Bool_t AliAnalysisTaskSAJF::IsJetCluster(AliEmcalJet* jet, Int_t iclus, Bool_t sorted) const
1058 {
1059   for (Int_t i = 0; i < jet->GetNumberOfClusters(); i++) {
1060     Int_t ijetclus = jet->ClusterAt(i);
1061     if (sorted && ijetclus > iclus)
1062       return kFALSE;
1063     if (ijetclus == iclus)
1064       return kTRUE;
1065   }
1066   return kFALSE;
1067 }
1068
1069 //________________________________________________________________________
1070 Bool_t AliAnalysisTaskSAJF::AcceptJet(Float_t eta, Float_t phi) const
1071 {
1072   return (Bool_t)(eta > fMinEta + fJetRadius && eta < fMaxEta - fJetRadius && phi > fMinPhi + fJetRadius && phi < fMaxPhi - fJetRadius);
1073 }
1074
1075 //________________________________________________________________________
1076 Bool_t AliAnalysisTaskSAJF::AcceptJet(AliEmcalJet* jet, Bool_t checkPt) const
1077 {
1078   if (checkPt && jet->MaxTrackPt() < fPtCutJetPart && jet->MaxClusterPt() < fPtCutJetPart)
1079     return kFALSE;
1080
1081   return AcceptJet(jet->Eta(), jet->Phi());
1082 }
1083
1084 //________________________________________________________________________
1085 Bool_t AliAnalysisTaskSAJF::AcceptCluster(AliVCluster* clus, Bool_t acceptMC)
1086 {
1087   if (!acceptMC && clus->Chi2() == 100)
1088     return kFALSE;
1089
1090   Double_t vertex[3] = {0, 0, 0};
1091   InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1092   TLorentzVector nPart;
1093   clus->GetMomentum(nPart, vertex);
1094
1095   if (nPart.Et() < fPtCut)
1096     return kFALSE;
1097
1098   return kTRUE;
1099 }
1100
1101 //________________________________________________________________________
1102 Bool_t AliAnalysisTaskSAJF::AcceptTrack(AliVTrack* track, Bool_t acceptMC) const
1103 {
1104   if (!acceptMC && track->GetLabel() == 100)
1105     return kFALSE;
1106
1107   if (track->Pt() < fPtCut)
1108     return kFALSE;
1109
1110   return (Bool_t)(track->Eta() > fMinEta && track->Eta() < fMaxEta && track->Phi() > fMinPhi && track->Phi() < fMaxPhi);
1111 }
1112
1113 //________________________________________________________________________
1114 void AliAnalysisTaskSAJF::UserExec(Option_t *)
1115 {
1116   Init();
1117
1118   RetrieveEventObjects();
1119
1120   FillHistograms();
1121
1122   // information for this iteration of the UserExec in the container
1123   PostData(1, fOutput);
1124 }
1125
1126 //________________________________________________________________________
1127 void AliAnalysisTaskSAJF::Terminate(Option_t *)
1128 {
1129   // Called once at the end of the analysis.
1130 }