PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAJF.cxx

   1 // $Id$
   2 //
   3 // Jet analysis task.
   4 //
   5 // Author: S.Aiola
   6
   7 #include <TClonesArray.h>
   8 #include <TH2F.h>
   9 #include <THnSparse.h>
  10 #include <TList.h>
  11 #include <TLorentzVector.h>
  12
  13 #include "AliVCluster.h"
  14 #include "AliVParticle.h"
  15 #include "AliEmcalJet.h"
  16 #include "AliRhoParameter.h"
  17 #include "AliLog.h"
  18
  19 #include "AliAnalysisTaskSAJF.h"
  20
  21 ClassImp(AliAnalysisTaskSAJF)
  22
  23 //________________________________________________________________________
  24 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF() :
  25   AliAnalysisTaskEmcalJet("AliAnalysisTaskSAJF", kTRUE),
  26   fHistEventObservables(0),
  27   fHistJetObservables(0)
  28
  29 {
  30   // Default constructor.
  31
  32   for (Int_t i = 0; i < 4; i++) {
  33     fHistTracksJetPt[i] = 0;
  34     fHistClustersJetPt[i] = 0;
  35     fHistTracksPtDist[i] = 0;
  36     fHistClustersPtDist[i] = 0;
  37   }
  38
  39   SetMakeGeneralHistograms(kTRUE);
  40 }
  41
  42 //________________________________________________________________________
  43 AliAnalysisTaskSAJF::AliAnalysisTaskSAJF(const char *name) :
  44   AliAnalysisTaskEmcalJet(name, kTRUE),
  45   fHistEventObservables(0),
  46   fHistJetObservables(0)
  47 {
  48   // Standard constructor.
  49
  50   for (Int_t i = 0; i < 4; i++) {
  51     fHistTracksJetPt[i] = 0;
  52     fHistClustersJetPt[i] = 0;
  53     fHistTracksPtDist[i] = 0;
  54     fHistClustersPtDist[i] = 0;
  55   }
  56
  57   SetMakeGeneralHistograms(kTRUE);
  58 }
  59
  60 //________________________________________________________________________
  61 void AliAnalysisTaskSAJF::UserCreateOutputObjects()
  62 {
  63   // Create user output.
  64
  65   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
  66
  67   TString title[20]= {""};
  68   Int_t nbins[20]  = {0};
  69   Double_t min[20] = {0.};
  70   Double_t max[20] = {0.};
  71   Int_t dim = 0;
  72
  73   if (fForceBeamType != kpp) {
  74     title[dim] = "Centrality (%)";
  75     nbins[dim] = 101;
  76     min[dim] = 0;
  77     max[dim] = 101;
  78     dim++;
  79
  80     title[dim] = "#psi_{RP} (rad)";
  81     nbins[dim] = 200;
  82     min[dim] = -TMath::Pi();
  83     max[dim] = TMath::Pi();
  84     dim++;
  85   }
  86
  87   title[dim] = "No. of jets";
  88   nbins[dim] = 150;
  89   min[dim] = -0.5;
  90   max[dim] = 149.5;
  91   dim++;
  92
  93   title[dim] = "p_{T,lead jet} (GeV/c)";
  94   nbins[dim] = fNbins;
  95   min[dim] = 0;
  96   max[dim] = 250;
  97   dim++;
  98
  99   title[dim] = "A_{lead jet}";
 100   nbins[dim] = 100;
 101   min[dim] = 0;
 102   max[dim] = 1;
 103   dim++;
 104
 105   if (!fRhoName.IsNull()) {
 106     title[dim] = "#rho (GeV/c)";
 107     nbins[dim] = fNbins;
 108     min[dim] = 0;
 109     max[dim] = 500;
 110     dim++;
 111
 112     title[dim] = "p_{T,lead jet}^{corr} (GeV/c)";
 113     nbins[dim] = fNbins*2;
 114     min[dim] = -250;
 115     max[dim] = 250;
 116     dim++;
 117   }
 118
 119   fHistEventObservables = new THnSparseD("fHistEventObservables","fHistEventObservables",dim,nbins,min,max);
 120   for (Int_t i = 0; i < dim; i++)
 121     fHistEventObservables->GetAxis(i)->SetTitle(title[i]);
 122
 123   dim = 0;
 124
 125   if (fForceBeamType != kpp) {
 126     title[dim] = "Centrality (%)";
 127     nbins[dim] = 101;
 128     min[dim] = 0;
 129     max[dim] = 101;
 130     dim++;
 131
 132     title[dim] = "#psi_{RP} (rad)";
 133     nbins[dim] = 200;
 134     min[dim] = -TMath::Pi();
 135     max[dim] = TMath::Pi();
 136     dim++;
 137   }
 138
 139   title[dim] = "#phi_{jet} (rad)";
 140   nbins[dim] = 201;
 141   min[dim] = 0;
 142   max[dim] = TMath::Pi()*201/100;
 143   dim++;
 144
 145   title[dim] = "#eta";
 146   nbins[dim] = 200;
 147   min[dim] = -1;
 148   max[dim] = 1;
 149   dim++;
 150
 151   title[dim] = "p_{T} (GeV/c)";
 152   nbins[dim] = fNbins;
 153   min[dim] = 0;
 154   max[dim] = 250;
 155   dim++;
 156
 157   title[dim] = "A_{jet}";
 158   nbins[dim] = 100;
 159   min[dim] = 0;
 160   max[dim] = 1;
 161   dim++;
 162
 163   if (fIsEmbedded) {
 164     title[dim] = "p_{T}^{MC} (GeV/c)";
 165     nbins[dim] = fNbins;
 166     min[dim] = 0;
 167     max[dim] = 250;
 168     dim++;
 169   }
 170
 171   if (!fRhoName.IsNull()) {
 172     title[dim] = "p_{T}^{corr} (GeV/c)";
 173     nbins[dim] = fNbins*2;
 174     min[dim] = -250;
 175     max[dim] = 250;
 176     dim++;
 177   }
 178
 179   title[dim] = "No. of constituents";
 180   nbins[dim] = 250;
 181   min[dim] = -0.5;
 182   max[dim] = 249.5;
 183   dim++;
 184
 185   title[dim] = "NEF";
 186   nbins[dim] = 120;
 187   min[dim] = 0;
 188   max[dim] = 1.2;
 189   dim++;
 190
 191   title[dim] = "Z";
 192   nbins[dim] = 120;
 193   min[dim] = 0;
 194   max[dim] = 1.2;
 195   dim++;
 196
 197   title[dim] = "p_{T,particle}^{leading} (GeV/c)";
 198   nbins[dim] = 120;
 199   min[dim] = 0;
 200   max[dim] = 120;
 201   dim++;
 202
 203   fHistJetObservables = new THnSparseD("fHistJetObservables","fHistJetObservables",dim,nbins,min,max);
 204   for (Int_t i = 0; i < dim; i++)
 205     fHistJetObservables->GetAxis(i)->SetTitle(title[i]);
 206
 207   for (Int_t i = 0; i < 4; i++) {
 208     TString histname;
 209
 210     if (!fTracksName.IsNull()) {
 211       histname = "fHistTracksJetPt_";
 212       histname += i;
 213       fHistTracksJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
 214       fHistTracksJetPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 215       fHistTracksJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
 216       fHistTracksJetPt[i]->GetZaxis()->SetTitle("counts");
 217       fOutput->Add(fHistTracksJetPt[i]);
 218
 219       histname = "fHistTracksPtDist_";
 220       histname += i;
 221       fHistTracksPtDist[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, 100, 0, 5);
 222       fHistTracksPtDist[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 223       fHistTracksPtDist[i]->GetYaxis()->SetTitle("d");
 224       fHistTracksPtDist[i]->GetZaxis()->SetTitle("counts");
 225       fOutput->Add(fHistTracksPtDist[i]);
 226     }
 227
 228     if (!fCaloName.IsNull()) {
 229       histname = "fHistClustersJetPt_";
 230       histname += i;
 231       fHistClustersJetPt[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, fNbins, fMinBinPt, fMaxBinPt);
 232       fHistClustersJetPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
 233       fHistClustersJetPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
 234       fHistClustersJetPt[i]->GetZaxis()->SetTitle("counts");
 235       fOutput->Add(fHistClustersJetPt[i]);
 236
 237       histname = "fHistClustersPtDist_";
 238       histname += i;
 239       fHistClustersPtDist[i] = new TH2F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2, 100, 0, 5);
 240       fHistClustersPtDist[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
 241       fHistClustersPtDist[i]->GetYaxis()->SetTitle("d");
 242       fHistClustersPtDist[i]->GetZaxis()->SetTitle("counts");
 243       fOutput->Add(fHistClustersPtDist[i]);
 244     }
 245   }
 246
 247   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 248
 249 }
 250
 251 //________________________________________________________________________
 252 Bool_t AliAnalysisTaskSAJF::FillHistograms()
 253 {
 254   // Fill histograms.
 255
 256   if (!fJets) {
 257     AliError(Form("%s - Jet array not provided, returning...", GetName()));
 258     return kFALSE;
 259   }
 260
 261   static Int_t sortedJets[9999] = {-1};
 262   Int_t nacc = GetSortedArray(sortedJets, fJets, fRhoVal);
 263
 264   Float_t corrPt = 0;
 265
 266   if (nacc > 0) {
 267     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
 268
 269     if (jet) {
 270       // Don't need to do AcceptJet since it was done already in GetSortedArray
 271       corrPt = jet->Pt() - fRhoVal * jet->Area();
 272
 273       DoJetLoop(nacc, sortedJets);
 274     }
 275     else {
 276       AliError(Form("Could not receive jet %d", sortedJets[0]));
 277     }
 278   }
 279
 280   // Fill THnSparse
 281
 282   return kTRUE;
 283 }
 284
 285 //________________________________________________________________________
 286 void AliAnalysisTaskSAJF::DoJetLoop(const Int_t nacc, const Int_t* sortedJets)
 287 {
 288   // Do the jet loop.
 289
 290   for (Int_t ij = 0; ij < nacc; ij++) {
 291
 292     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[ij]));
 293
 294     if (!jet) {
 295       AliError(Form("Could not receive jet %d", ij));
 296       continue;
 297     }
 298
 299     // Don't need to do AcceptJet since it was done already in GetSortedArray
 300
 301     Float_t ptLeading = GetLeadingHadronPt(jet);
 302     Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();
 303
 304     // Fill THnSparse
 305
 306     if (fTracks) {
 307       for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 308         AliVParticle *track = jet->TrackAt(it, fTracks);
 309         if (track) {
 310           fHistTracksJetPt[fCentBin]->Fill(track->Pt(), jet->Pt());
 311           Double_t dist = TMath::Sqrt((track->Eta() - jet->Eta()) * (track->Eta() - jet->Eta()) + (track->Phi() - jet->Phi()) * (track->Phi() - jet->Phi()));
 312           fHistTracksPtDist[fCentBin]->Fill(track->Pt(), dist);
 313         }
 314       }
 315     }
 316
 317     if (fCaloClusters) {
 318       for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 319         AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 320
 321         if (cluster) {
 322           TLorentzVector nPart;
 323           cluster->GetMomentum(nPart, fVertex);
 324
 325           fHistClustersJetPt[fCentBin]->Fill(nPart.Pt(), jet->Pt());
 326           Double_t dist = TMath::Sqrt((nPart.Eta() - jet->Eta()) * (nPart.Eta() - jet->Eta()) + (nPart.Phi() - jet->Phi()) * (nPart.Phi() - jet->Phi()));
 327           fHistClustersPtDist[fCentBin]->Fill(nPart.Pt(), dist);
 328         }
 329       }
 330     }
 331   } //jet loop
 332 }