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