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