PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetSample.cxx

   1 // $Id$
   2 //
   3 // Jet sample analysis task.
   4 //
   5 // Author: S.Aiola
   6
   7 #include <TClonesArray.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TList.h>
  11 #include <TLorentzVector.h>
  12
  13 #include "AliVCluster.h"
  14 #include "AliVTrack.h"
  15 #include "AliEmcalJet.h"
  16 #include "AliRhoParameter.h"
  17 #include "AliLog.h"
  18
  19 #include "AliAnalysisTaskEmcalJetSample.h"
  20
  21 ClassImp(AliAnalysisTaskEmcalJetSample)
  22
  23 //________________________________________________________________________
  24 AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample() :
  25   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetSample", kTRUE)
  26
  27 {
  28   // Default constructor.
  29
  30   for (Int_t i = 0; i < 4; i++) {
  31     fHistTracksPt[i] = 0;
  32     fHistClustersPt[i] = 0;
  33     fHistLeadingJetPt[i] = 0;
  34     fHistJetsPhiEta[i] = 0;
  35     fHistJetsPtArea[i] = 0;
  36     fHistJetsPtLeadHad[i] = 0;
  37     fHistJetsCorrPtArea[i] = 0;
  38   }
  39
  40   SetMakeGeneralHistograms(kTRUE);
  41 }
  42
  43 //________________________________________________________________________
  44 AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample(const char *name) :
  45   AliAnalysisTaskEmcalJet(name, kTRUE)
  46 {
  47   // Standard constructor.
  48
  49   for (Int_t i = 0; i < 4; i++) {
  50     fHistTracksPt[i] = 0;
  51     fHistClustersPt[i] = 0;
  52     fHistLeadingJetPt[i] = 0;
  53     fHistJetsPhiEta[i] = 0;
  54     fHistJetsPtArea[i] = 0;
  55     fHistJetsPtLeadHad[i] = 0;
  56     fHistJetsCorrPtArea[i] = 0;
  57   }
  58
  59   SetMakeGeneralHistograms(kTRUE);
  60 }
  61
  62 //________________________________________________________________________
  63 AliAnalysisTaskEmcalJetSample::~AliAnalysisTaskEmcalJetSample()
  64 {
  65   // Destructor.
  66 }
  67
  68 //________________________________________________________________________
  69 void AliAnalysisTaskEmcalJetSample::UserCreateOutputObjects()
  70 {
  71   // Create user output.
  72
  73   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
  74
  75   TString histname;
  76
  77   for (Int_t i = 0; i < 4; i++) {
  78     if (!fTracksName.IsNull()) {
  79       histname = "fHistTracksPt_";
  80       histname += i;
  81       fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
  82       fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
  83       fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
  84       fOutput->Add(fHistTracksPt[i]);
  85     }
  86
  87     if (!fCaloName.IsNull()) {
  88       histname = "fHistClustersPt_";
  89       histname += i;
  90       fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
  91       fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
  92       fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
  93       fOutput->Add(fHistClustersPt[i]);
  94     }
  95
  96     if (!fJetsName.IsNull()) {
  97       histname = "fHistLeadingJetPt_";
  98       histname += i;
  99       fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 100       fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 101       fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 102       fOutput->Add(fHistLeadingJetPt[i]);
 103
 104       histname = "fHistJetsPhiEta_";
 105       histname += i;
 106       fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
 107       fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
 108       fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
 109       fOutput->Add(fHistJetsPhiEta[i]);
 110
 111       histname = "fHistJetsPtArea_";
 112       histname += i;
 113       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
 114       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 115       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 116       fOutput->Add(fHistJetsPtArea[i]);
 117
 118       histname = "fHistJetsPtLeadHad_";
 119       histname += i;
 120       fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 121       fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 122       fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
 123       fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
 124       fOutput->Add(fHistJetsPtLeadHad[i]);
 125
 126       if (!fRhoName.IsNull()) {
 127         histname = "fHistJetsCorrPtArea_";
 128         histname += i;
 129         fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
 130         fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
 131         fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
 132         fOutput->Add(fHistJetsCorrPtArea[i]);
 133       }
 134     }
 135   }
 136   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 137 }
 138
 139 //________________________________________________________________________
 140 Bool_t AliAnalysisTaskEmcalJetSample::FillHistograms()
 141 {
 142   // Fill histograms.
 143
 144   if (fTracks) {
 145     const Int_t ntracks = fTracks->GetEntriesFast();
 146
 147     for (Int_t it = 0; it < ntracks; it++) {
 148       AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(it));
 149
 150       if (!track) {
 151         AliError(Form("Could not receive track %d", it));
 152         continue;
 153       }
 154
 155       if (!AcceptTrack(track))
 156         continue;
 157
 158       fHistTracksPt[fCentBin]->Fill(track->Pt());
 159     }
 160   }
 161
 162   if (fCaloClusters) {
 163     const Int_t nclusters = fCaloClusters->GetEntriesFast();
 164
 165     for (Int_t ic = 0; ic < nclusters; ic++) {
 166       AliVCluster *cluster = static_cast<AliVCluster*>(fCaloClusters->At(ic));
 167
 168       if (!cluster) {
 169         AliError(Form("Could not receive cluster %d", ic));
 170         continue;
 171       }
 172
 173       TLorentzVector nPart;
 174       cluster->GetMomentum(nPart, fVertex);
 175       fHistClustersPt[fCentBin]->Fill(nPart.Pt());
 176     }
 177   }
 178
 179   if (fJets) {
 180     static Int_t sortedJets[9999] = {-1};
 181     Bool_t r = GetSortedArray(sortedJets, fJets);
 182
 183     if (r && sortedJets[0]>=0) {
 184       AliEmcalJet* leadJet = static_cast<AliEmcalJet*>(fJets->At(sortedJets[0]));
 185       if (leadJet)
 186         fHistLeadingJetPt[fCentBin]->Fill(leadJet->Pt());
 187       else
 188         AliError("Could not retrieve leading jet!");
 189     }
 190
 191     const Int_t njets = fJets->GetEntriesFast();
 192     for (Int_t ij = 0; ij < njets; ij++) {
 193
 194       AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
 195       if (!jet) {
 196         AliError(Form("Could not receive jet %d", ij));
 197         continue;
 198       }
 199
 200       if (!AcceptJet(jet))
 201         continue;
 202
 203       fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
 204       fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
 205
 206       Float_t ptLeading = GetLeadingHadronPt(jet);
 207       fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
 208
 209       if (fRho) {
 210         Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();
 211         fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
 212       }
 213     }
 214   }
 215
 216   return kTRUE;
 217 }
 218
 219 //________________________________________________________________________
 220 Bool_t AliAnalysisTaskEmcalJetSample::Run()
 221 {
 222   // Run analysis code here, if needed. It will be executed before FillHistograms().
 223
 224   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 225 }
 226
 227 //________________________________________________________________________
 228 void AliAnalysisTaskEmcalJetSample::Terminate(Option_t *)
 229 {
 230   // Called once at the end of the analysis.
 231 }