PWGGA/EMCALJetTasks/AliAnalysisTaskSAJF.cxx

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