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