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 "AliRhoParameter.h"
  25 #include "AliLog.h"
  26
  27 #include "AliAnalysisTaskSAJF.h"
  28
  29 ClassImp(AliAnalysisTaskSAJF)
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
  33   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAJF", kTRUE),
  34   fMinRC2LJ(1.0),
  35   fEmbJetsName("EmbJets"),
  36   fEmbTracksName(""),
  37   fEmbCaloName(""),
  38   fRandTracksName("TracksRandomized"),
  39   fRandCaloName("CaloClustersRandomized"),
  40   fRhoName("Rho"),
  41   fEmbJets(0),
  42   fEmbTracks(0),
  43   fEmbCaloClusters(0),
  44   fRandTracks(0),
  45   fRandCaloClusters(0),
  46   fRho(0),
  47   fHistCentrality(0),
  48   fHistRhoVSleadJetPt(0),
  49   fHistRCPhiEta(0),
  50   fHistRCPtExLJVSDPhiLJ(0),
  51   fHistRhoVSRCPt(0),
  52   fHistEmbJetPhiEta(0),
  53   fHistEmbPartPhiEta(0),
  54   fHistRhoVSEmbBkg(0)
  55
  56 {
  57   // Default constructor.
  58
  59   for (Int_t i = 0; i < 4; i++) {
  60     fHistEvents[i] = 0;
  61     fHistTracksPt[i] = 0;
  62     fHistClustersPt[i] = 0;
  63     fHistJetPhiEta[i] = 0;
  64     fHistJetsPt[i] = 0;
  65     fHistJetsPtArea[i] = 0;
  66     fHistLeadingJetPt[i] = 0;
  67     fHist2LeadingJetPt[i] = 0;
  68     fHistJetsNEFvsPt[i] = 0;
  69     fHistJetsZvsPt[i] = 0;
  70     fHistMaxTrackPtvsJetPt[i] = 0;
  71     fHistMaxClusPtvsJetPt[i] = 0;
  72     fHistMaxPartPtvsJetPt[i] = 0;
  73     fHistMaxTrackPtvsJetCorrPt[i] = 0;
  74     fHistMaxClusPtvsJetCorrPt[i] = 0;
  75     fHistMaxPartPtvsJetCorrPt[i] = 0;
  76     fHistRho[i] = 0;
  77     fHistCorrJetsPt[i] = 0;
  78     fHistCorrJetsPtArea[i] = 0;
  79     fHistCorrLeadingJetPt[i] = 0;
  80     fHistRCPt[i] = 0;
  81     fHistRCPtExLJ[i] = 0;
  82     fHistRCPtRand[i] = 0;
  83     fHistDeltaPtRC[i] = 0;
  84     fHistDeltaPtRCExLJ[i] = 0;
  85     fHistDeltaPtRCRand[i] = 0;
  86     fHistEmbJetsPt[i] = 0;
  87     fHistEmbJetsCorrPt[i] = 0;
  88     fHistEmbPart[i] = 0;
  89     fHistDeltaPtEmb[i] = 0;
  90   }
  91 }
  92
  93 //________________________________________________________________________
  94 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
  95   AliAnalysisTaskEmcalJet(name, kTRUE),
  96   fMinRC2LJ(1.0),
  97   fEmbJetsName("EmbJets"),
  98   fEmbTracksName(""),
  99   fEmbCaloName(""),
 100   fRandTracksName("TracksRandomized"),
 101   fRandCaloName("CaloClustersRandomized"),
 102   fRhoName("Rho"),
 103   fEmbJets(0),
 104   fEmbTracks(0),
 105   fEmbCaloClusters(0),
 106   fRandTracks(0),
 107   fRandCaloClusters(0),
 108   fRho(0),
 109   fHistCentrality(0),
 110   fHistRhoVSleadJetPt(0),
 111   fHistRCPhiEta(0),
 112   fHistRCPtExLJVSDPhiLJ(0),
 113   fHistRhoVSRCPt(0),
 114   fHistEmbJetPhiEta(0),
 115   fHistEmbPartPhiEta(0),
 116   fHistRhoVSEmbBkg(0)
 117 {
 118   // Standard constructor.
 119
 120   for (Int_t i = 0; i < 4; i++) {
 121     fHistEvents[i] = 0;
 122     fHistTracksPt[i] = 0;
 123     fHistClustersPt[i] = 0;
 124     fHistJetPhiEta[i] = 0;
 125     fHistJetsPt[i] = 0;
 126     fHistJetsPtArea[i] = 0;
 127     fHistLeadingJetPt[i] = 0;
 128     fHist2LeadingJetPt[i] = 0;
 129     fHistJetsNEFvsPt[i] = 0;
 130     fHistJetsZvsPt[i] = 0;
 131     fHistMaxTrackPtvsJetPt[i] = 0;
 132     fHistMaxClusPtvsJetPt[i] = 0;
 133     fHistMaxPartPtvsJetPt[i] = 0;
 134     fHistMaxTrackPtvsJetCorrPt[i] = 0;
 135     fHistMaxClusPtvsJetCorrPt[i] = 0;
 136     fHistMaxPartPtvsJetCorrPt[i] = 0;
 137     fHistRho[i] = 0;
 138     fHistCorrJetsPt[i] = 0;
 139     fHistCorrJetsPtArea[i] = 0;
 140     fHistCorrLeadingJetPt[i] = 0;
 141     fHistRCPt[i] = 0;
 142     fHistRCPtExLJ[i] = 0;
 143     fHistRCPtRand[i] = 0;
 144     fHistDeltaPtRC[i] = 0;
 145     fHistDeltaPtRCExLJ[i] = 0;
 146     fHistDeltaPtRCRand[i] = 0;
 147     fHistEmbJetsPt[i] = 0;
 148     fHistEmbJetsCorrPt[i] = 0;
 149     fHistEmbPart[i] = 0;
 150     fHistDeltaPtEmb[i] = 0;
 151   }
 152 }
 153
 154 //________________________________________________________________________
 155 AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
 156 {
 157   // Destructor.
 158 }
 159
 160 //________________________________________________________________________
 161 void AliAnalysisTaskSAJF::UserCreateOutputObjects()
 162 {
 163   // Create user output.
 164
 165   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 166
 167   const Float_t binWidth = (fMaxBinPt - fMinBinPt) / fNbins;
 168
 169   OpenFile(1);
 170   fOutput = new TList();
 171   fOutput->SetOwner();
 172
 173   fHistCentrality = new TH1F("fHistCentrality","fHistCentrality", fNbins, 0, 100);
 174   fHistCentrality->GetXaxis()->SetTitle("Centrality (%)");
 175   fHistCentrality->GetYaxis()->SetTitle("counts");
 176   fOutput->Add(fHistCentrality);
 177
 178   fHistRhoVSleadJetPt = new TH2F("fHistRhoVSleadJetPt","fHistRhoVSleadJetPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 179   fHistRhoVSleadJetPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
 180   fHistRhoVSleadJetPt->GetYaxis()->SetTitle("Leading jet p_{T} [GeV/c]");
 181   fOutput->Add(fHistRhoVSleadJetPt);
 182
 183   fHistRCPhiEta = new TH2F("fHistRCPhiEta","Phi-Eta distribution of rigid cones", 40, -2, 2, 64, 0, 6.4);
 184   fHistRCPhiEta->GetXaxis()->SetTitle("#eta");
 185   fHistRCPhiEta->GetYaxis()->SetTitle("#phi");
 186   fOutput->Add(fHistRCPhiEta);
 187
 188   fHistRCPtExLJVSDPhiLJ = new TH2F("fHistRCPtExLJVSDPhiLJ","fHistRCPtExLJVSDPhiLJ", fNbins, fMinBinPt, fMaxBinPt, 128, -1.6, 4.8);
 189   fHistRCPtExLJVSDPhiLJ->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 190   fHistRCPtExLJVSDPhiLJ->GetYaxis()->SetTitle("#Delta#phi");
 191   fOutput->Add(fHistRCPtExLJVSDPhiLJ);
 192
 193   fHistRhoVSRCPt = new TH2F("fHistRhoVSRCPt","fHistRhoVSRCPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 194   fHistRhoVSRCPt->GetXaxis()->SetTitle("#rho * area [GeV/c]");
 195   fHistRhoVSRCPt->GetYaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 196   fOutput->Add(fHistRhoVSRCPt);
 197
 198   if (!fEmbJetsName.IsNull()) {
 199     fHistEmbJetPhiEta = new TH2F("fHistEmbJetPhiEta","Phi-Eta distribution of embedded jets", 40, -2, 2, 64, 0, 6.4);
 200     fHistEmbJetPhiEta->GetXaxis()->SetTitle("#eta");
 201     fHistEmbJetPhiEta->GetYaxis()->SetTitle("#phi");
 202     fOutput->Add(fHistEmbJetPhiEta);
 203
 204     fHistEmbPartPhiEta = new TH2F("fHistEmbPartPhiEta","Phi-Eta distribution of embedded particles", 40, -2, 2, 64, 0, 6.4);
 205     fHistEmbPartPhiEta->GetXaxis()->SetTitle("#eta");
 206     fHistEmbPartPhiEta->GetYaxis()->SetTitle("#phi");
 207     fOutput->Add(fHistEmbPartPhiEta);
 208
 209     fHistRhoVSEmbBkg = new TH2F("fHistRhoVSEmbBkg","fHistRhoVSEmbBkg", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 210     fHistRhoVSEmbBkg->GetXaxis()->SetTitle("rho * area [GeV/c]");
 211     fHistRhoVSEmbBkg->GetYaxis()->SetTitle("background of embedded track [GeV/c]");
 212     fOutput->Add(fHistRhoVSEmbBkg);
 213   }
 214
 215   TString histname;
 216
 217   for (Int_t i = 0; i < 4; i++) {
 218     histname = "fHistEvents_";
 219     histname += i;
 220     fHistEvents[i] = new TH1F(histname,histname, 6, 0, 6);
 221     fHistEvents[i]->GetXaxis()->SetTitle("Event state");
 222     fHistEvents[i]->GetYaxis()->SetTitle("counts");
 223     fHistEvents[i]->GetXaxis()->SetBinLabel(1, "Rho <= 0");
 224     fHistEvents[i]->GetXaxis()->SetBinLabel(2, "No jets");
 225     fHistEvents[i]->GetXaxis()->SetBinLabel(3, "Max Jet not found");
 226     fHistEvents[i]->GetXaxis()->SetBinLabel(4, "Max Jet <= 0");
 227     fHistEvents[i]->GetXaxis()->SetBinLabel(5, "Accepted");
 228     fOutput->Add(fHistEvents[i]);
 229
 230     histname = "fHistTracksPt_";
 231     histname += i;
 232     fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
 233     fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 234     fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
 235     fOutput->Add(fHistTracksPt[i]);
 236
 237     if (fAnaType == kEMCAL) {
 238       histname = "fHistClustersPt_";
 239       histname += i;
 240       fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
 241       fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 242       fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
 243       fOutput->Add(fHistClustersPt[i]);
 244     }
 245
 246     histname = "fHistJetPhiEta_";
 247     histname += i;
 248     fHistJetPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 40, -2, 2, 64, 0, 6.4);
 249     fHistJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
 250     fHistJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
 251     fOutput->Add(fHistJetPhiEta[i]);
 252
 253     histname = "fHistJetsPt_";
 254     histname += i;
 255     fHistJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 256     fHistJetsPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
 257     fHistJetsPt[i]->GetYaxis()->SetTitle("counts");
 258     fOutput->Add(fHistJetsPt[i]);
 259
 260     histname = "fHistJetsPtArea_";
 261     histname += i;
 262     fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
 263     fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
 264     fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 265     fOutput->Add(fHistJetsPtArea[i]);
 266
 267     histname = "fHistLeadingJetPt_";
 268     histname += i;
 269     fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 270     fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV]");
 271     fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 272     fOutput->Add(fHistLeadingJetPt[i]);
 273
 274     histname = "fHist2LeadingJetPt_";
 275     histname += i;
 276     fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 277     fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} [GeV]");
 278     fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 279     fOutput->Add(fHist2LeadingJetPt[i]);
 280
 281     histname = "fHistJetsZvsPt_";
 282     histname += i;
 283     fHistJetsZvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 284     fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
 285     fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
 286     fOutput->Add(fHistJetsZvsPt[i]);
 287
 288     if (fAnaType == kEMCAL) {
 289       histname = "fHistJetsNEFvsPt_";
 290       histname += i;
 291       fHistJetsNEFvsPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 292       fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
 293       fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} [GeV/c]");
 294       fOutput->Add(fHistJetsNEFvsPt[i]);
 295     }
 296
 297     histname = "fHistMaxTrackPtvsJetPt_";
 298     histname += i;
 299     fHistMaxTrackPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
 300     fHistMaxTrackPtvsJetPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
 301     fHistMaxTrackPtvsJetPt[i]->GetYaxis()->SetTitle("p_{T}^{track} [GeV/c]");
 302     fOutput->Add(fHistMaxTrackPtvsJetPt[i]);
 303
 304     if (fAnaType == kEMCAL) {
 305       histname = "fHistMaxClusPtvsJetPt_";
 306       histname += i;
 307       fHistMaxClusPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
 308       fHistMaxClusPtvsJetPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
 309       fHistMaxClusPtvsJetPt[i]->GetYaxis()->SetTitle("p_{T}^{clus} [GeV/c]");
 310       fOutput->Add(fHistMaxClusPtvsJetPt[i]);
 311     }
 312
 313     histname = "fHistMaxPartPtvsJetPt_";
 314     histname += i;
 315     fHistMaxPartPtvsJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
 316     fHistMaxPartPtvsJetPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
 317     fHistMaxPartPtvsJetPt[i]->GetYaxis()->SetTitle("p_{T}^{part} [GeV/c]");
 318     fOutput->Add(fHistMaxPartPtvsJetPt[i]);
 319
 320     histname = "fHistMaxTrackPtvsJetCorrPt_";
 321     histname += i;
 322     fHistMaxTrackPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
 323     fHistMaxTrackPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
 324     fHistMaxTrackPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{track} [GeV/c]");
 325     fOutput->Add(fHistMaxTrackPtvsJetCorrPt[i]);
 326
 327     if (fAnaType == kEMCAL) {
 328       histname = "fHistMaxClusPtvsJetCorrPt_";
 329       histname += i;
 330       fHistMaxClusPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
 331       fHistMaxClusPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
 332       fHistMaxClusPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{clus} [GeV/c]");
 333       fOutput->Add(fHistMaxClusPtvsJetCorrPt[i]);
 334     }
 335
 336     histname = "fHistMaxPartPtvsJetCorrPt_";
 337     histname += i;
 338     fHistMaxPartPtvsJetCorrPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt, fNbins / 2.5, fMinBinPt, fMaxBinPt / 2.5);
 339     fHistMaxPartPtvsJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{jet} [GeV/c]");
 340     fHistMaxPartPtvsJetCorrPt[i]->GetYaxis()->SetTitle("p_{T}^{part} [GeV/c]");
 341     fOutput->Add(fHistMaxPartPtvsJetCorrPt[i]);
 342
 343     histname = "fHistRho_";
 344     histname += i;
 345     fHistRho[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 346     fHistRho[i]->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 347     fHistRho[i]->GetYaxis()->SetTitle("counts");
 348     fOutput->Add(fHistRho[i]);
 349
 350     histname = "fHistCorrJetsPt_";
 351     histname += i;
 352     fHistCorrJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 353     fHistCorrJetsPt[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
 354     fHistCorrJetsPt[i]->GetYaxis()->SetTitle("counts");
 355     fOutput->Add(fHistCorrJetsPt[i]);
 356
 357     histname = "fHistCorrJetsPtArea_";
 358     histname += i;
 359     fHistCorrJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 20, 0, fJetRadius * fJetRadius * TMath::Pi() * 1.5);
 360     fHistCorrJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
 361     fHistCorrJetsPtArea[i]->GetYaxis()->SetTitle("area");
 362     fOutput->Add(fHistCorrJetsPtArea[i]);
 363
 364     histname = "fHistCorrLeadingJetPt_";
 365     histname += i;
 366     fHistCorrLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 367     fHistCorrLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{RC} [GeV/c]");
 368     fHistCorrLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 369     fOutput->Add(fHistCorrLeadingJetPt[i]);
 370
 371     histname = "fHistRCPt_";
 372     histname += i;
 373     fHistRCPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 374     fHistRCPt[i]->GetXaxis()->SetTitle("rigid cone p_{T} [GeV/c]");
 375     fHistRCPt[i]->GetYaxis()->SetTitle("counts");
 376     fOutput->Add(fHistRCPt[i]);
 377
 378     histname = "fHistRCPtExLJ_";
 379     histname += i;
 380     fHistRCPtExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 381     fHistRCPtExLJ[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
 382     fHistRCPtExLJ[i]->GetYaxis()->SetTitle("counts");
 383     fOutput->Add(fHistRCPtExLJ[i]);
 384
 385     histname = "fHistRCPtRand_";
 386     histname += i;
 387     fHistRCPtRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt * 2);
 388     fHistRCPtRand[i]->GetXaxis()->SetTitle("rigid cone p_{T}^{RC} [GeV/c]");
 389     fHistRCPtRand[i]->GetYaxis()->SetTitle("counts");
 390     fOutput->Add(fHistRCPtRand[i]);
 391
 392     histname = "fHistDeltaPtRC_";
 393     histname += i;
 394     fHistDeltaPtRC[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
 395     fHistDeltaPtRC[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
 396     fHistDeltaPtRC[i]->GetYaxis()->SetTitle("counts");
 397     fOutput->Add(fHistDeltaPtRC[i]);
 398
 399     histname = "fHistDeltaPtRCExLJ_";
 400     histname += i;
 401     fHistDeltaPtRCExLJ[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
 402     fHistDeltaPtRCExLJ[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
 403     fHistDeltaPtRCExLJ[i]->GetYaxis()->SetTitle("counts");
 404     fOutput->Add(fHistDeltaPtRCExLJ[i]);
 405
 406     histname = "fHistDeltaPtRCRand_";
 407     histname += i;
 408     fHistDeltaPtRCRand[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
 409     fHistDeltaPtRCRand[i]->GetXaxis()->SetTitle("#deltap_{T}^{RC} [GeV/c]");
 410     fHistDeltaPtRCRand[i]->GetYaxis()->SetTitle("counts");
 411     fOutput->Add(fHistDeltaPtRCRand[i]);
 412
 413     if (!fEmbJetsName.IsNull()) {
 414       histname = "fHistEmbJetsPt_";
 415       histname += i;
 416       fHistEmbJetsPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 417       fHistEmbJetsPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{raw} [GeV/c]");
 418       fHistEmbJetsPt[i]->GetYaxis()->SetTitle("counts");
 419       fOutput->Add(fHistEmbJetsPt[i]);
 420
 421       histname = "fHistEmbJetsCorrPt_";
 422       histname += i;
 423       fHistEmbJetsCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 424       fHistEmbJetsCorrPt[i]->GetXaxis()->SetTitle("embedded jet p_{T}^{corr} [GeV/c]");
 425       fHistEmbJetsCorrPt[i]->GetYaxis()->SetTitle("counts");
 426       fOutput->Add(fHistEmbJetsCorrPt[i]);
 427
 428       histname = "fHistEmbPart_";
 429       histname += i;
 430       fHistEmbPart[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 431       fHistEmbPart[i]->GetXaxis()->SetTitle("embedded particle p_{T}^{emb} [GeV/c]");
 432       fHistEmbPart[i]->GetYaxis()->SetTitle("counts");
 433       fOutput->Add(fHistEmbPart[i]);
 434
 435       histname = "fHistDeltaPtEmb_";
 436       histname += i;
 437       fHistDeltaPtEmb[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt - fMaxBinPt / 2 + binWidth / 2, fMinBinPt + fMaxBinPt / 2 + binWidth / 2);
 438       fHistDeltaPtEmb[i]->GetXaxis()->SetTitle("#deltap_{T}^{emb} [GeV/c]");
 439       fHistDeltaPtEmb[i]->GetYaxis()->SetTitle("counts");
 440       fOutput->Add(fHistDeltaPtEmb[i]);
 441     }
 442   }
 443
 444   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 445 }
 446
 447 //________________________________________________________________________
 448 Bool_t AliAnalysisTaskSAJF::RetrieveEventObjects()
 449 {
 450   // Retrieve event objects.
 451
 452   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 453     return kFALSE;
 454
 455   if (!fRhoName.IsNull() && !fRho) {
 456     fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
 457     if (!fRho) {
 458       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
 459       return kFALSE;
 460     }
 461   }
 462
 463   if (!fEmbJetsName.IsNull() && !fEmbJets) {
 464     fEmbJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbJetsName));
 465     if (!fEmbJets) {
 466       AliError(Form("%s: Could not retrieve emb jets %s!", GetName(), fEmbJetsName.Data()));
 467       return kFALSE;
 468     }
 469     else if (!fEmbJets->GetClass()->GetBaseClass("AliEmcalJet")) {
 470       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fEmbJetsName.Data()));
 471       fEmbJets = 0;
 472       return kFALSE;
 473     }
 474   }
 475
 476   if (!fEmbCaloName.IsNull() && fAnaType == kEMCAL && !fEmbCaloClusters) {
 477     fEmbCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbCaloName));
 478     if (!fEmbCaloClusters) {
 479       AliError(Form("%s: Could not retrieve embedded clusters %s!", GetName(), fEmbCaloName.Data()));
 480       return kFALSE;
 481     }
 482     else if (!fEmbCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fEmbCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
 483       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fEmbCaloName.Data()));
 484       fEmbCaloClusters = 0;
 485       return kFALSE;
 486     }
 487   }
 488
 489   if (!fEmbTracksName.IsNull() && !fEmbTracks) {
 490     fEmbTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fEmbTracksName));
 491     if (!fEmbTracks) {
 492       AliError(Form("%s: Could not retrieve embedded tracks %s!", GetName(), fEmbTracksName.Data()));
 493       return kFALSE;
 494     }
 495     else if (!fEmbTracks->GetClass()->GetBaseClass("AliVParticle") && !fEmbTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
 496       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fEmbTracksName.Data()));
 497       fEmbTracks = 0;
 498       return kFALSE;
 499     }
 500   }
 501
 502   if (!fRandCaloName.IsNull() && fAnaType == kEMCAL && !fRandCaloClusters) {
 503     fRandCaloClusters =  dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandCaloName));
 504     if (!fRandCaloClusters) {
 505       AliError(Form("%s: Could not retrieve randomized clusters %s!", GetName(), fRandCaloName.Data()));
 506       return kFALSE;
 507     }
 508     else if (!fRandCaloClusters->GetClass()->GetBaseClass("AliVCluster") && !fRandCaloClusters->GetClass()->GetBaseClass("AliEmcalParticle")) {
 509       AliError(Form("%s: Collection %s does not contain AliVCluster nor AliEmcalParticle objects!", GetName(), fRandCaloName.Data()));
 510       fRandCaloClusters = 0;
 511       return kFALSE;
 512     }
 513   }
 514
 515   if (!fRandTracksName.IsNull() && !fRandTracks) {
 516     fRandTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fRandTracksName));
 517     if (!fRandTracks) {
 518       AliError(Form("%s: Could not retrieve randomized tracks %s!", GetName(), fRandTracksName.Data()));
 519       return kFALSE;
 520     }
 521     else if (!fRandTracks->GetClass()->GetBaseClass("AliVParticle") && !fRandTracks->GetClass()->GetBaseClass("AliEmcalParticle")) {
 522       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fRandTracksName.Data()));
 523       fRandTracks = 0;
 524       return kFALSE;
 525     }
 526   }
 527
 528   return kTRUE;
 529 }
 530
 531 //________________________________________________________________________
 532 Bool_t AliAnalysisTaskSAJF::FillHistograms()
 533 {
 534   // Fill histograms.
 535
 536   // Before filling any histogram, check if the event is interesting
 537   Double_t rho = fRho->GetVal();
 538
 539   if (rho < 0) {  // rho < 0, probably a diffractive event, skipping
 540     fHistEvents[fCentBin]->Fill("Rho <= 0", 1);
 541     return kTRUE;
 542   }
 543
 544   Int_t maxJetIndex  = -1;
 545   Int_t max2JetIndex = -1;
 546
 547   GetLeadingJets(maxJetIndex, max2JetIndex);
 548
 549   if (maxJetIndex < 0) { // no accept jet, skipping
 550     fHistEvents[fCentBin]->Fill("No jets", 1);
 551     return kTRUE;
 552   }
 553
 554   AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(maxJetIndex));
 555   if (!jet) {  // error, I cannot get the lead jet from collection (should never happen), skipping
 556     fHistEvents[fCentBin]->Fill("Max Jet not found", 1);
 557     return kTRUE;
 558   }
 559
 560   // OK, event accepted
 561
 562   Float_t maxJetCorrPt = jet->Pt() - rho * jet->Area();
 563
 564   if (maxJetCorrPt <= 0)
 565     fHistEvents[fCentBin]->Fill("Max Jet <= 0", 1);
 566
 567   fHistEvents[fCentBin]->Fill("Accepted", 1);
 568
 569   // ************
 570   // General histograms
 571   // _________________________________
 572
 573   DoJetLoop();
 574   DoTrackLoop();
 575   DoClusterLoop();
 576
 577   fHistCentrality->Fill(fCent);
 578   fHistRho[fCentBin]->Fill(rho);
 579
 580   if (jet) {
 581     fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
 582     fHistRhoVSleadJetPt->Fill(rho * jet->Area(), jet->Pt());
 583     fHistCorrLeadingJetPt[fCentBin]->Fill(maxJetCorrPt);
 584   }
 585
 586   AliEmcalJet* jet2 = 0;
 587   if (max2JetIndex >= 0)
 588     jet2 = static_cast<AliEmcalJet*>(fJets->At(max2JetIndex));
 589
 590   if (jet2)
 591     fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());
 592
 593   // ************
 594   // Random cones
 595   // _________________________________
 596
 597   const Float_t rcArea = fJetRadius * fJetRadius * TMath::Pi();
 598   const Bool_t IsMCEvent = (Bool_t)(fTracksName.Contains("Randomized") || fTracksName.Contains("Embedded"));
 599
 600   // Simple random cones
 601   Float_t RCpt = 0;
 602   Float_t RCeta = 0;
 603   Float_t RCphi = 0;
 604   GetRigidCone(RCpt, RCeta, RCphi, IsMCEvent, 0);
 605   if (RCpt > 0) {
 606     fHistRCPt[fCentBin]->Fill(RCpt / rcArea);
 607     fHistDeltaPtRC[fCentBin]->Fill(RCpt - rcArea * rho);
 608   }
 609
 610   // Random cones far from leading jet
 611   Float_t RCptExLJ = 0;
 612   Float_t RCetaExLJ = 0;
 613   Float_t RCphiExLJ = 0;
 614   GetRigidCone(RCptExLJ, RCetaExLJ, RCphiExLJ, IsMCEvent, jet);
 615   if (RCptExLJ > 0) {
 616     fHistRCPhiEta->Fill(RCetaExLJ, RCphiExLJ);
 617     fHistRhoVSRCPt->Fill(rho, RCptExLJ / rcArea);
 618
 619     Float_t dphi = RCphiExLJ - jet->Phi();
 620     if (dphi > 4.8) dphi -= TMath::Pi() * 2;
 621     if (dphi < -1.6) dphi += TMath::Pi() * 2;
 622     fHistRCPtExLJVSDPhiLJ->Fill(RCptExLJ, dphi);
 623
 624     fHistRCPtExLJ[fCentBin]->Fill(RCptExLJ / rcArea);
 625     fHistDeltaPtRCExLJ[fCentBin]->Fill(RCptExLJ - rcArea * rho);
 626   }
 627
 628   // Random cones with randomized particles
 629   Float_t RCptRand = 0;
 630   Float_t RCetaRand = 0;
 631   Float_t RCphiRand = 0;
 632   GetRigidCone(RCptRand, RCetaRand, RCphiRand, kTRUE, 0, fRandTracks, fRandCaloClusters);
 633   if (RCptRand > 0) {
 634     fHistRCPtRand[fCentBin]->Fill(RCptRand / rcArea);
 635     fHistDeltaPtRCRand[fCentBin]->Fill(RCptRand - rcArea * rho);
 636   }
 637
 638   // ************
 639   // Embedding
 640   // _________________________________
 641
 642   if (!fEmbJets)
 643     return kTRUE;
 644
 645   AliEmcalJet *embJet  = 0;
 646   TObject     *embPart = 0;
 647
 648   DoEmbJetLoop(embJet, embPart);
 649
 650   if (embJet) {
 651
 652     fHistEmbJetsPt[fCentBin]->Fill(embJet->Pt());
 653     fHistEmbJetsCorrPt[fCentBin]->Fill(embJet->Pt() - rho * embJet->Area());
 654     fHistEmbPart[fCentBin]->Fill(embJet->MCPt());
 655     fHistEmbJetPhiEta->Fill(embJet->Eta(), embJet->Phi());
 656     fHistDeltaPtEmb[fCentBin]->Fill(embJet->Pt() - embJet->Area() * rho - embJet->MCPt());
 657     fHistRhoVSEmbBkg->Fill(embJet->Area() * rho, embJet->Pt() - embJet->MCPt());
 658
 659     AliVCluster *cluster = dynamic_cast<AliVCluster*>(embPart);
 660     if (cluster) {
 661       Float_t pos[3];
 662       cluster->GetPosition(pos);
 663       TVector3 clusVec(pos);
 664
 665       fHistEmbPartPhiEta->Fill(clusVec.Eta(), clusVec.Phi());
 666     }
 667     else {
 668       AliVTrack *track = dynamic_cast<AliVTrack*>(embPart);
 669       if (track) {
 670         fHistEmbPartPhiEta->Fill(track->Eta(), track->Phi());
 671       }
 672       else {
 673         AliWarning(Form("%s - Embedded particle type not found or not recognized (neither AliVCluster nor AliVTrack)!", GetName()));
 674         return kTRUE;
 675       }
 676     }
 677
 678   }
 679   else {
 680     AliWarning(Form("%s - Embedded jet not found in the event!", GetName()));
 681   }
 682
 683   return kTRUE;
 684 }
 685
 686 //________________________________________________________________________
 687 void AliAnalysisTaskSAJF::GetLeadingJets(Int_t &maxJetIndex, Int_t &max2JetIndex)
 688 {
 689   // Get the leading jets.
 690
 691   if (!fJets)
 692     return;
 693
 694   const Double_t rho = fRho->GetVal();
 695
 696   const Int_t njets = fJets->GetEntriesFast();
 697
 698   Float_t maxJetPt = -999;
 699   Float_t max2JetPt = -999;
 700   for (Int_t ij = 0; ij < njets; ij++) {
 701
 702     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
 703
 704     if (!jet) {
 705       AliError(Form("Could not receive jet %d", ij));
 706       continue;
 707     }
 708
 709     if (!AcceptJet(jet))
 710       continue;
 711
 712     Float_t corrPt = jet->Pt() - rho * jet->Area();
 713
 714     if (maxJetIndex == -1 || maxJetPt < corrPt) {
 715       max2JetPt = maxJetPt;
 716       max2JetIndex = maxJetIndex;
 717       maxJetPt = corrPt;
 718       maxJetIndex = ij;
 719     }
 720     else if (max2JetIndex == -1 || max2JetPt < corrPt) {
 721       max2JetPt = corrPt;
 722       max2JetIndex = ij;
 723     }
 724   }
 725 }
 726
 727 //________________________________________________________________________
 728 void AliAnalysisTaskSAJF::DoClusterLoop()
 729 {
 730   // Do cluster loop.
 731
 732   if (!fCaloClusters)
 733     return;
 734
 735   Int_t nclusters =  fCaloClusters->GetEntriesFast();
 736
 737   for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 738     AliVCluster* cluster = static_cast<AliVCluster*>(fCaloClusters->At(iClusters));
 739     if (!cluster) {
 740       AliError(Form("Could not receive cluster %d", iClusters));
 741       continue;
 742     }
 743
 744     if (!AcceptCluster(cluster, kTRUE))
 745       continue;
 746
 747     fHistClustersPt[fCentBin]->Fill(cluster->E());
 748   }
 749 }
 750
 751 //________________________________________________________________________
 752 void AliAnalysisTaskSAJF::DoTrackLoop()
 753 {
 754   // Do track loop.
 755
 756   if (!fTracks)
 757     return;
 758
 759   Int_t ntracks = fTracks->GetEntriesFast();
 760
 761   for (Int_t i = 0; i < ntracks; i++) {
 762
 763     AliVParticle* track = static_cast<AliVParticle*>(fTracks->At(i)); // pointer to reconstructed to track
 764
 765     if (!track) {
 766       AliError(Form("Could not retrieve track %d",i));
 767       continue;
 768     }
 769
 770     AliVTrack* vtrack = dynamic_cast<AliVTrack*>(track);
 771
 772     if (vtrack && !AcceptTrack(vtrack, kTRUE))
 773       continue;
 774
 775     fHistTracksPt[fCentBin]->Fill(track->Pt());
 776   }
 777 }
 778
 779
 780 //________________________________________________________________________
 781 void AliAnalysisTaskSAJF::DoJetLoop()
 782 {
 783   // Do the jet loop.
 784
 785   if (!fJets)
 786     return;
 787
 788   const Double_t rho = fRho->GetVal();
 789
 790   const Int_t njets = fJets->GetEntriesFast();
 791
 792   for (Int_t ij = 0; ij < njets; ij++) {
 793
 794     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
 795
 796     if (!jet) {
 797       AliError(Form("Could not receive jet %d", ij));
 798       continue;
 799     }
 800
 801     if (!AcceptJet(jet))
 802       continue;
 803
 804     Float_t corrPt = jet->Pt() - rho * jet->Area();
 805
 806     fHistJetsPt[fCentBin]->Fill(jet->Pt());
 807     fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
 808     fHistCorrJetsPt[fCentBin]->Fill(corrPt);
 809     fHistCorrJetsPtArea[fCentBin]->Fill(corrPt, jet->Area());
 810
 811     fHistMaxTrackPtvsJetPt[fCentBin]->Fill(jet->Pt(), jet->MaxTrackPt());
 812     fHistMaxPartPtvsJetPt[fCentBin]->Fill(jet->Pt(), jet->MaxPartPt());
 813
 814     fHistMaxTrackPtvsJetCorrPt[fCentBin]->Fill(corrPt, jet->MaxTrackPt());
 815     fHistMaxPartPtvsJetCorrPt[fCentBin]->Fill(corrPt, jet->MaxPartPt());
 816
 817     fHistJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
 818
 819     if (fAnaType == kEMCAL) {
 820       fHistMaxClusPtvsJetPt[fCentBin]->Fill(jet->Pt(), jet->MaxClusterPt());
 821       fHistMaxClusPtvsJetCorrPt[fCentBin]->Fill(corrPt, jet->MaxClusterPt());
 822       fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt());
 823     }
 824
 825     if (fTracks) {
 826       for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 827         AliVParticle *track = jet->TrackAt(it, fTracks);
 828         if (track)
 829           fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt());
 830       }
 831     }
 832
 833     if (fCaloClusters) {
 834       for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 835         AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 836
 837         if (cluster) {
 838           TLorentzVector nPart;
 839           cluster->GetMomentum(nPart, fVertex);
 840           fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), jet->Pt());
 841         }
 842       }
 843     }
 844   } //jet loop
 845 }
 846
 847 //________________________________________________________________________
 848 void AliAnalysisTaskSAJF::DoEmbJetLoop(AliEmcalJet* &embJet, TObject* &embPart)
 849 {
 850   // Do the embedded jet loop.
 851
 852   if (!fEmbJets)
 853     return;
 854
 855   TLorentzVector maxClusVect;
 856
 857   Int_t nembjets = fEmbJets->GetEntriesFast();
 858
 859   for (Int_t ij = 0; ij < nembjets; ij++) {
 860
 861     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fEmbJets->At(ij));
 862
 863     if (!jet) {
 864       AliError(Form("Could not receive jet %d", ij));
 865       continue;
 866     }
 867
 868     if (!AcceptJet(jet))
 869       continue;
 870
 871     if (!jet->IsMC())
 872       continue;
 873
 874     AliVParticle *maxTrack = 0;
 875
 876     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 877       AliVParticle *track = jet->TrackAt(it, fEmbTracks);
 878
 879       if (!track)
 880         continue;
 881
 882       if (track->GetLabel() != 100)
 883         continue;
 884
 885       if (!maxTrack || track->Pt() > maxTrack->Pt())
 886         maxTrack = track;
 887     }
 888
 889     AliVCluster *maxClus = 0;
 890
 891     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 892       AliVCluster *cluster = jet->ClusterAt(ic, fEmbCaloClusters);
 893
 894       if (!cluster)
 895         continue;
 896
 897       if (cluster->Chi2() != 100)
 898         continue;
 899
 900       TLorentzVector nPart;
 901       cluster->GetMomentum(nPart, fVertex);
 902
 903       if (!maxClus || nPart.Et() > maxClusVect.Et()) {
 904         new (&maxClusVect) TLorentzVector(nPart);
 905         maxClus = cluster;
 906       }
 907     }
 908
 909     if ((maxClus && maxTrack && maxClusVect.Et() > maxTrack->Pt()) || (maxClus && !maxTrack)) {
 910       embPart = maxClus;
 911       embJet = jet;
 912     }
 913     else if (maxTrack) {
 914       embPart = maxTrack;
 915       embJet = jet;
 916     }
 917
 918     return;  // MC jets found, exit
 919   }
 920 }
 921
 922 //________________________________________________________________________
 923 void AliAnalysisTaskSAJF::GetRigidCone(Float_t &pt, Float_t &eta, Float_t &phi, Bool_t acceptMC,
 924                                        AliEmcalJet *jet, TClonesArray* tracks, TClonesArray* clusters) const
 925 {
 926   // Get rigid cone.
 927
 928   if (!tracks)
 929     tracks = fTracks;
 930
 931   if (!clusters)
 932     clusters = fCaloClusters;
 933
 934   if (!tracks && !clusters)
 935     return;
 936
 937   eta = 0;
 938   phi = 0;
 939   pt = 0;
 940
 941   Float_t LJeta = 999;
 942   Float_t LJphi = 999;
 943
 944   if (jet) {
 945     LJeta = jet->Eta();
 946     LJphi = jet->Phi();
 947   }
 948
 949   Float_t maxEta = fMaxEta;
 950   Float_t minEta = fMinEta;
 951   Float_t maxPhi = fMaxPhi;
 952   Float_t minPhi = fMinPhi;
 953
 954   if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
 955   if (minPhi < 0) minPhi = 0;
 956
 957   Float_t dLJ = 0;
 958   Int_t repeats = 0;
 959   do {
 960     eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
 961     phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
 962     dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
 963     repeats++;
 964   } while (dLJ < fMinRC2LJ && repeats < 999);
 965
 966   if (repeats == 999) {
 967     AliWarning("Could not get random cone!");
 968     return;
 969   }
 970
 971   if (fAnaType == kEMCAL && clusters) {
 972     Int_t nclusters =  clusters->GetEntriesFast();
 973     for (Int_t iClusters = 0; iClusters < nclusters; iClusters++) {
 974       AliVCluster* cluster = static_cast<AliVCluster*>(clusters->At(iClusters));
 975       if (!cluster) {
 976         AliError(Form("Could not receive cluster %d", iClusters));
 977         continue;
 978       }
 979
 980       if (!AcceptCluster(cluster, acceptMC))
 981         continue;
 982
 983       TLorentzVector nPart;
 984       cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
 985
 986       Float_t pos[3];
 987       cluster->GetPosition(pos);
 988       TVector3 clusVec(pos);
 989
 990       Float_t d = TMath::Sqrt((clusVec.Eta() - eta) * (clusVec.Eta() - eta) + (clusVec.Phi() - phi) * (clusVec.Phi() - phi));
 991
 992       if (d <= fJetRadius)
 993         pt += nPart.Pt();
 994     }
 995   }
 996
 997   if (tracks) {
 998     Int_t ntracks = tracks->GetEntriesFast();
 999     for(Int_t iTracks = 0; iTracks < ntracks; iTracks++) {
1000       AliVTrack* track = static_cast<AliVTrack*>(tracks->At(iTracks));
1001       if(!track) {
1002         AliError(Form("Could not retrieve track %d",iTracks));
1003         continue;
1004       }
1005
1006       if (!AcceptTrack(track, acceptMC))
1007         continue;
1008
1009       Float_t tracketa = track->Eta();
1010       Float_t trackphi = track->Phi();
1011
1012       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
1013         trackphi += 2 * TMath::Pi();
1014       if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
1015         trackphi -= 2 * TMath::Pi();
1016
1017       Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
1018       if (d <= fJetRadius)
1019         pt += track->Pt();
1020     }
1021   }
1022 }
1023 //________________________________________________________________________
1024 void AliAnalysisTaskSAJF::Init()
1025 {
1026   // Initialize the analysis.
1027
1028   AliAnalysisTaskEmcalJet::Init();
1029
1030   if (!fEmbJetsName.IsNull()) {
1031     if (fEmbTracksName.IsNull()) {
1032       fEmbTracksName = fTracksName;
1033       fEmbTracksName += "Embedded";
1034     }
1035     if (fEmbCaloName.IsNull() && fAnaType == kEMCAL) {
1036       fEmbCaloName = fCaloName;
1037       fEmbCaloName += "Embedded";
1038     }
1039   }
1040
1041   const Float_t semiDiag = TMath::Sqrt((fMaxPhi - fMinPhi) * (fMaxPhi - fMinPhi) +
1042                                        (fMaxEta - fMinEta) * (fMaxEta - fMinEta)) / 2;
1043   if (fMinRC2LJ > semiDiag * 0.5) {
1044     AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
1045                     "Will use fMinRC2LJ = %f", fMinRC2LJ, semiDiag * 0.5));
1046     fMinRC2LJ = semiDiag * 0.5;
1047   }
1048 }
1049
1050 //________________________________________________________________________
1051 void AliAnalysisTaskSAJF::Terminate(Option_t *)
1052 {
1053   // Called once at the end of the analysis.
1054 }