PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAJF.cxx

   1 // $Id$
   2 //
   3 // Jet analysis task.
   4 //
   5 // Author: S.Aiola
   6
   7 #include <TClonesArray.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13
  14 #include "AliVCluster.h"
  15 #include "AliVParticle.h"
  16 #include "AliVTrack.h"
  17 #include "AliEmcalJet.h"
  18 #include "AliRhoParameter.h"
  19 #include "AliLog.h"
  20
  21 #include "AliAnalysisTaskSAJF.h"
  22
  23 ClassImp(AliAnalysisTaskSAJF)
  24
  25 //________________________________________________________________________
  26 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
  27   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAJF", kTRUE),
  28   fNjetsVsCent(0)
  29
  30 {
  31   // Default constructor.
  32
  33   for (Int_t i = 0; i < 4; i++) {
  34     fHistEvents[i] = 0;
  35     fHistLeadingJetPt[i] = 0;
  36     fHist2LeadingJetPt[i] = 0;
  37     fHistLeadingJetCorrPt[i] = 0;
  38     fHistRhoVSleadJetPt[i] = 0;
  39     fHistJetPhiEta[i] = 0;
  40     fHistJetsPtArea[i] = 0;
  41     fHistJetsCorrPtArea[i] = 0;
  42     fHistJetsNEFvsPt[i] = 0;
  43     fHistJetsZvsPt[i] = 0;
  44     fHistConstituents[i] = 0;
  45     fHistTracksJetPt[i] = 0;
  46     fHistClustersJetPt[i] = 0;
  47     fHistTracksPtDist[i] = 0;
  48     fHistClustersPtDist[i] = 0;
  49     fHistJetNconstVsPt[i] = 0;
  50   }
  51
  52   SetMakeGeneralHistograms(kTRUE);
  53 }
  54
  55 //________________________________________________________________________
  56 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
  57   AliAnalysisTaskEmcalJet(name, kTRUE),
  58   fNjetsVsCent(0)
  59 {
  60   // Standard constructor.
  61
  62   for (Int_t i = 0; i < 4; i++) {
  63     fHistEvents[i] = 0;
  64     fHistLeadingJetPt[i] = 0;
  65     fHist2LeadingJetPt[i] = 0;
  66     fHistLeadingJetCorrPt[i] = 0;
  67     fHistRhoVSleadJetPt[i] = 0;
  68     fHistJetPhiEta[i] = 0;
  69     fHistJetsPtArea[i] = 0;
  70     fHistJetsCorrPtArea[i] = 0;
  71     fHistJetsNEFvsPt[i] = 0;
  72     fHistJetsZvsPt[i] = 0;
  73     fHistConstituents[i] = 0;
  74     fHistTracksJetPt[i] = 0;
  75     fHistClustersJetPt[i] = 0;
  76     fHistTracksPtDist[i] = 0;
  77     fHistClustersPtDist[i] = 0;
  78     fHistJetNconstVsPt[i] = 0;
  79   }
  80
  81   SetMakeGeneralHistograms(kTRUE);
  82 }
  83
  84 //________________________________________________________________________
  85 AliAnalysisTaskSAJF::~AliAnalysisTaskSAJF()
  86 {
  87   // Destructor.
  88 }
  89
  90 //________________________________________________________________________
  91 void AliAnalysisTaskSAJF::UserCreateOutputObjects()
  92 {
  93   // Create user output.
  94
  95   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
  96
  97   fNjetsVsCent = new TH2F("fNjetsVsCent","fNjetsVsCent", 100, 0, 100, 150, -0.5, 149.5);
  98   fNjetsVsCent->GetXaxis()->SetTitle("Centrality (%)");
  99   fNjetsVsCent->GetYaxis()->SetTitle("# of jets");
 100   fOutput->Add(fNjetsVsCent);
 101
 102   const Int_t nbinsZ = 12;
 103   Float_t binsZ[nbinsZ+1] = {0,1,2,3,4,5,6,7,8,9,10,20,1000};
 104
 105   Float_t *binsPt       = GenerateFixedBinArray(fNbins, fMinBinPt, fMaxBinPt);
 106   Float_t *binsCorrPt   = GenerateFixedBinArray(fNbins*2, -fMaxBinPt, fMaxBinPt);
 107   Float_t *binsArea     = GenerateFixedBinArray(30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
 108   Float_t *binsEta      = GenerateFixedBinArray(50,-1, 1);
 109   Float_t *binsPhi      = GenerateFixedBinArray(101, 0, TMath::Pi() * 2.02);
 110   Float_t *bins120      = GenerateFixedBinArray(120, 0, 1.2);
 111   Float_t *bins150      = GenerateFixedBinArray(150, -0.5, 149.5);
 112
 113   TString histname;
 114
 115   for (Int_t i = 0; i < 4; i++) {
 116     histname = "fHistEvents_";
 117     histname += i;
 118     fHistEvents[i] = new TH1F(histname,histname, 6, 0, 6);
 119     fHistEvents[i]->GetXaxis()->SetTitle("Event state");
 120     fHistEvents[i]->GetYaxis()->SetTitle("counts");
 121     fHistEvents[i]->GetXaxis()->SetBinLabel(1, "No jets");
 122     fHistEvents[i]->GetXaxis()->SetBinLabel(2, "Max jet not found");
 123     fHistEvents[i]->GetXaxis()->SetBinLabel(3, "Rho == 0");
 124     fHistEvents[i]->GetXaxis()->SetBinLabel(4, "Max jet <= 0");
 125     fHistEvents[i]->GetXaxis()->SetBinLabel(5, "OK");
 126     fOutput->Add(fHistEvents[i]);
 127
 128     histname = "fHistLeadingJetPt_";
 129     histname += i;
 130     fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 131     fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 132     fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 133     fOutput->Add(fHistLeadingJetPt[i]);
 134
 135     histname = "fHist2LeadingJetPt_";
 136     histname += i;
 137     fHist2LeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 138     fHist2LeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 139     fHist2LeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 140     fOutput->Add(fHist2LeadingJetPt[i]);
 141
 142     histname = "fHistLeadingJetCorrPt_";
 143     histname += i;
 144     fHistLeadingJetCorrPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins * 2, -fMaxBinPt, fMaxBinPt);
 145     fHistLeadingJetCorrPt[i]->GetXaxis()->SetTitle("p_{T}^{corr} (GeV/c)");
 146     fHistLeadingJetCorrPt[i]->GetYaxis()->SetTitle("counts");
 147     fOutput->Add(fHistLeadingJetCorrPt[i]);
 148
 149     histname = "fHistRhoVSleadJetPt_";
 150     histname += i;
 151     fHistRhoVSleadJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt*2, fNbins, fMinBinPt, fMaxBinPt);
 152     fHistRhoVSleadJetPt[i]->GetXaxis()->SetTitle("#rho * area (GeV/c)");
 153     fHistRhoVSleadJetPt[i]->GetYaxis()->SetTitle("Leading jet p_{T} (GeV/c)");
 154     fOutput->Add(fHistRhoVSleadJetPt[i]);
 155
 156     histname = "fHistJetPhiEta_";
 157     histname += i;
 158     fHistJetPhiEta[i] = new TH3F(histname.Data(), histname.Data(),
 159                                  50, binsEta,
 160                                  101, binsPhi,
 161                                  nbinsZ, binsZ);
 162     fHistJetPhiEta[i]->GetXaxis()->SetTitle("#eta");
 163     fHistJetPhiEta[i]->GetYaxis()->SetTitle("#phi");
 164     fHistJetPhiEta[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
 165     fOutput->Add(fHistJetPhiEta[i]);
 166
 167     histname = "fHistJetsPtArea_";
 168     histname += i;
 169     fHistJetsPtArea[i] = new TH3F(histname.Data(), histname.Data(),
 170                                   fNbins, binsPt,
 171                                   30, binsArea,
 172                                   nbinsZ, binsZ);
 173     fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 174     fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 175     fHistJetsPtArea[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
 176     fOutput->Add(fHistJetsPtArea[i]);
 177
 178     histname = "fHistJetsZvsPt_";
 179     histname += i;
 180     fHistJetsZvsPt[i] = new TH3F(histname.Data(), histname.Data(),
 181                                  120, bins120,
 182                                  fNbins, binsPt,
 183                                  nbinsZ, binsZ);
 184     fHistJetsZvsPt[i]->GetXaxis()->SetTitle("Z");
 185     fHistJetsZvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 186     fHistJetsZvsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
 187     fOutput->Add(fHistJetsZvsPt[i]);
 188
 189     if (!fCaloName.IsNull()) {
 190       histname = "fHistJetsNEFvsPt_";
 191       histname += i;
 192       fHistJetsNEFvsPt[i] = new TH3F(histname.Data(), histname.Data(),
 193                                      120, bins120,
 194                                      fNbins, binsPt,
 195                                      nbinsZ, binsZ);
 196       fHistJetsNEFvsPt[i]->GetXaxis()->SetTitle("NEF");
 197       fHistJetsNEFvsPt[i]->GetYaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 198       fHistJetsNEFvsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
 199       fOutput->Add(fHistJetsNEFvsPt[i]);
 200     }
 201
 202     histname = "fHistJetsCorrPtArea_";
 203     histname += i;
 204     fHistJetsCorrPtArea[i] = new TH3F(histname.Data(), histname.Data(),
 205                                       fNbins * 2, binsCorrPt,
 206                                       30, binsArea,
 207                                       nbinsZ, binsZ);
 208     fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
 209     fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
 210     fOutput->Add(fHistJetsCorrPtArea[i]);
 211
 212     histname = "fHistConstituents_";
 213     histname += i;
 214     fHistConstituents[i] = new TH2F(histname.Data(), histname.Data(), 100, 1, 101, 100, -0.5, 99.5);
 215     fHistConstituents[i]->GetXaxis()->SetTitle("p_{T,part} (GeV/c)");
 216     fHistConstituents[i]->GetYaxis()->SetTitle("no. of particles");
 217     fHistConstituents[i]->GetZaxis()->SetTitle("counts");
 218     fOutput->Add(fHistConstituents[i]);
 219
 220     histname = "fHistTracksJetPt_";
 221     histname += i;
 222     fHistTracksJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
 223     fHistTracksJetPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 224     fHistTracksJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
 225     fHistTracksJetPt[i]->GetZaxis()->SetTitle("counts");
 226     fOutput->Add(fHistTracksJetPt[i]);
 227
 228     histname = "fHistTracksPtDist_";
 229     histname += i;
 230     fHistTracksPtDist[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, 100, 0, 5);
 231     fHistTracksPtDist[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 232     fHistTracksPtDist[i]->GetYaxis()->SetTitle("d");
 233     fHistTracksPtDist[i]->GetZaxis()->SetTitle("counts");
 234     fOutput->Add(fHistTracksPtDist[i]);
 235
 236     if (!fCaloName.IsNull()) {
 237       histname = "fHistClustersJetPt_";
 238       histname += i;
 239       fHistClustersJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
 240       fHistClustersJetPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
 241       fHistClustersJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
 242       fHistClustersJetPt[i]->GetZaxis()->SetTitle("counts");
 243       fOutput->Add(fHistClustersJetPt[i]);
 244
 245       histname = "fHistClustersPtDist_";
 246       histname += i;
 247       fHistClustersPtDist[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, 100, 0, 5);
 248       fHistClustersPtDist[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
 249       fHistClustersPtDist[i]->GetYaxis()->SetTitle("d");
 250       fHistClustersPtDist[i]->GetZaxis()->SetTitle("counts");
 251       fOutput->Add(fHistClustersPtDist[i]);
 252     }
 253
 254     histname = "fHistJetNconstVsPt_";
 255     histname += i;
 256     fHistJetNconstVsPt[i] = new TH3F(histname.Data(), histname.Data(), 150, bins150, fNbins, binsPt, nbinsZ, binsZ);
 257     fHistJetNconstVsPt[i]->GetXaxis()->SetTitle("# of constituents");
 258     fHistJetNconstVsPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
 259     fHistJetNconstVsPt[i]->GetZaxis()->SetTitle("p_{T,lead} (GeV/c)");
 260     fOutput->Add(fHistJetNconstVsPt[i]);
 261   }
 262
 263   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 264
 265   delete[] binsPt;
 266   delete[] binsCorrPt;
 267   delete[] binsArea;
 268   delete[] binsEta;
 269   delete[] binsPhi;
 270   delete[] bins120;
 271 }
 272
 273 //________________________________________________________________________
 274 Bool_t AliAnalysisTaskSAJF::FillHistograms()
 275 {
 276   // Fill histograms.
 277
 278   if (!fJets) {
 279     AliError(Form("%s - Jet array not provided, returning...", GetName()));
 280     return kFALSE;
 281   }
 282
 283   if (fJets->GetEntriesFast() < 1) { // no jets in array, skipping
 284     fHistEvents[fCentBin]->Fill("No jets", 1);
 285     return kTRUE;
 286   }
 287
 288   Int_t *sortedJets = GetSortedArray(fJets);
 289
 290   if (!sortedJets || sortedJets[0] < 0) { // no accepted jets, skipping
 291     fHistEvents[fCentBin]->Fill("No jets", 1);
 292     return kTRUE;
 293   }
 294
 295   AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
 296   if (!jet) {  // error, I cannot get the leading jet from collection (should never happen), skipping
 297     fHistEvents[fCentBin]->Fill("Max jet not found", 1);
 298     return kTRUE;
 299   }
 300
 301   // OK, event accepted
 302
 303   Float_t maxJetCorrPt = jet->Pt() - fRhoVal * jet->Area();
 304
 305   if (fRhoVal == 0)
 306     fHistEvents[fCentBin]->Fill("Rho == 0", 1);
 307
 308   else if (maxJetCorrPt <= 0)
 309     fHistEvents[fCentBin]->Fill("Max jet <= 0", 1);
 310
 311   else
 312     fHistEvents[fCentBin]->Fill("OK", 1);
 313
 314   if (jet) {
 315     fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
 316     fHistRhoVSleadJetPt[fCentBin]->Fill(fRhoVal, jet->Pt());
 317     fHistLeadingJetCorrPt[fCentBin]->Fill(maxJetCorrPt);
 318   }
 319
 320   AliEmcalJet* jet2 = 0;
 321   if (sortedJets[1] >= 0)
 322     jet2 = static_cast<AliEmcalJet*>(fJets->At(sortedJets[1]));
 323
 324   if (jet2)
 325     fHist2LeadingJetPt[fCentBin]->Fill(jet2->Pt());
 326
 327   Int_t njets = DoJetLoop();
 328
 329   fNjetsVsCent->Fill(fCent, njets);
 330
 331   return kTRUE;
 332 }
 333
 334 //________________________________________________________________________
 335 Int_t AliAnalysisTaskSAJF::DoJetLoop()
 336 {
 337   // Do the jet loop.
 338
 339   if (!fJets)
 340     return 0;
 341
 342   const Int_t njets = fJets->GetEntriesFast();
 343   Int_t nAccJets = 0;
 344
 345   TH1F constituents("constituents", "constituents",
 346                     fHistConstituents[0]->GetNbinsX(), fHistConstituents[0]->GetXaxis()->GetXmin(), fHistConstituents[0]->GetXaxis()->GetXmax());
 347
 348   for (Int_t ij = 0; ij < njets; ij++) {
 349
 350     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
 351
 352     if (!jet) {
 353       AliError(Form("Could not receive jet %d", ij));
 354       continue;
 355     }
 356
 357     if (!AcceptJet(jet))
 358       continue;
 359
 360     Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();
 361
 362     Float_t ptLeading = GetLeadingHadronPt(jet);
 363
 364     fHistJetPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi(), ptLeading);
 365     fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area(), ptLeading);
 366     fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area(), ptLeading);
 367     fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt(), ptLeading);
 368
 369     if (fCaloClusters)
 370       fHistJetsNEFvsPt[fCentBin]->Fill(jet->NEF(), jet->Pt(), ptLeading);
 371
 372     if (fTracks) {
 373       for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 374         AliVParticle *track = jet->TrackAt(it, fTracks);
 375         if (track) {
 376           fHistJetsZvsPt[fCentBin]->Fill(track->Pt() / jet->Pt(), jet->Pt(), ptLeading);
 377           constituents.Fill(track->Pt());
 378           fHistTracksJetPt[fCentBin]->Fill(track->Pt(), jet->Pt());
 379           Double_t dist = TMath::Sqrt((track->Eta() - jet->Eta()) * (track->Eta() - jet->Eta()) + (track->Phi() - jet->Phi()) * (track->Phi() - jet->Phi()));
 380           fHistTracksPtDist[fCentBin]->Fill(track->Pt(), dist);
 381         }
 382       }
 383     }
 384
 385     if (fCaloClusters) {
 386       for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 387         AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 388
 389         if (cluster) {
 390           TLorentzVector nPart;
 391           cluster->GetMomentum(nPart, fVertex);
 392           fHistJetsZvsPt[fCentBin]->Fill(nPart.Et() / jet->Pt(), jet->Pt(), ptLeading);
 393           constituents.Fill(nPart.Pt());
 394           fHistClustersJetPt[fCentBin]->Fill(nPart.Pt(), jet->Pt());
 395           Double_t dist = TMath::Sqrt((nPart.Eta() - jet->Eta()) * (nPart.Eta() - jet->Eta()) + (nPart.Phi() - jet->Phi()) * (nPart.Phi() - jet->Phi()));
 396           fHistClustersPtDist[fCentBin]->Fill(nPart.Pt(), dist);
 397         }
 398       }
 399     }
 400
 401     for (Int_t i = 1; i <= constituents.GetNbinsX(); i++) {
 402       fHistConstituents[fCentBin]->Fill(constituents.GetBinCenter(i), constituents.GetBinContent(i));
 403     }
 404
 405     constituents.Reset();
 406     nAccJets++;
 407   } //jet loop
 408
 409   return nAccJets;
 410 }
 411
 412 //________________________________________________________________________
 413 void AliAnalysisTaskSAJF::Terminate(Option_t *)
 414 {
 415   // Called once at the end of the analysis.
 416 }