3 // Base class for rho calculation.
4 // Calculates parameterized rho for given centrality independent of input.
11 #include <TClonesArray.h>
14 #include "AliRhoParameter.h"
15 #include "AliEmcalJet.h"
17 #include "AliAnalysisTaskRhoBase.h"
19 ClassImp(AliAnalysisTaskRhoBase)
21 //________________________________________________________________________
22 AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase() :
23 AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoBase", kFALSE),
27 fCompareRhoScaledName(),
30 fInEventSigmaRho(35.83),
31 fAttachToEvent(kTRUE),
37 fHistJetAreavsCent(0),
40 fHistJetPtvsNtrack(0),
41 fHistJetAreavsNtrack(0),
44 fHistRhoScaledvsCent(0),
45 fHistDeltaRhovsCent(0),
46 fHistDeltaRhoScalevsCent(0),
48 fHistRhoScaledvsNtrack(0),
49 fHistDeltaRhovsNtrack(0),
50 fHistDeltaRhoScalevsNtrack(0),
51 fHistRhovsNcluster(0),
52 fHistRhoScaledvsNcluster(0)
56 for (Int_t i = 0; i < 4; i++) {
57 fHistJetNconstVsPt[i] = 0;
58 fHistJetRhovsEta[i] = 0;
60 for (Int_t i = 0; i < 12; i++) {
61 fHistNjUEoverNjVsNj[i] = 0;
65 //________________________________________________________________________
66 AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase(const char *name, Bool_t histo) :
67 AliAnalysisTaskEmcalJet(name, histo),
71 fCompareRhoScaledName(),
74 fInEventSigmaRho(35.83),
75 fAttachToEvent(kTRUE),
81 fHistJetAreavsCent(0),
84 fHistJetPtvsNtrack(0),
85 fHistJetAreavsNtrack(0),
88 fHistRhoScaledvsCent(0),
89 fHistDeltaRhovsCent(0),
90 fHistDeltaRhoScalevsCent(0),
92 fHistRhoScaledvsNtrack(0),
93 fHistDeltaRhovsNtrack(0),
94 fHistDeltaRhoScalevsNtrack(0),
95 fHistRhovsNcluster(0),
96 fHistRhoScaledvsNcluster(0)
100 for (Int_t i = 0; i < 4; i++) {
101 fHistJetNconstVsPt[i] = 0;
102 fHistJetRhovsEta[i] = 0;
104 for (Int_t i = 0; i < 12; i++) {
105 fHistNjUEoverNjVsNj[i] = 0;
107 SetMakeGeneralHistograms(histo);
110 //________________________________________________________________________
111 void AliAnalysisTaskRhoBase::UserCreateOutputObjects()
113 // User create output objects, called at the beginning of the analysis.
118 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
120 fHistRhovsCent = new TH2F("fHistRhovsCent", "fHistRhovsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt*2);
121 fHistRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
122 fHistRhovsCent->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
123 fOutput->Add(fHistRhovsCent);
125 if (fParticleCollArray.GetEntriesFast()>0) {
126 fHistRhovsNtrack = new TH2F("fHistRhovsNtrack", "fHistRhovsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt*2);
127 fHistRhovsNtrack->GetXaxis()->SetTitle("No. of tracks");
128 fHistRhovsNtrack->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
129 fOutput->Add(fHistRhovsNtrack);
132 if (fClusterCollArray.GetEntriesFast()>0) {
133 fHistRhovsNcluster = new TH2F("fHistRhovsNcluster", "fHistRhovsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
134 fHistRhovsNcluster->GetXaxis()->SetTitle("No. of tracks");
135 fHistRhovsNcluster->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
136 fOutput->Add(fHistRhovsNcluster);
139 if (fJetCollArray.GetEntriesFast()>0) {
140 fHistJetPtvsCent = new TH2F("fHistJetPtvsCent", "fHistJetPtvsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt);
141 fHistJetPtvsCent->GetXaxis()->SetTitle("Centrality (%)");
142 fHistJetPtvsCent->GetYaxis()->SetTitle("#it{p}_{T,jet} (GeV/c)");
143 fOutput->Add(fHistJetPtvsCent);
145 fHistJetAreavsCent = new TH2F("fHistJetAreavsCent", "fHistJetAreavsCent", 101, -1, 100, 100, 0, 1);
146 fHistJetAreavsCent->GetXaxis()->SetTitle("Centrality (%)");
147 fHistJetAreavsCent->GetYaxis()->SetTitle("Jet area");
148 fOutput->Add(fHistJetAreavsCent);
150 fHistJetRhovsCent = new TH2F("fHistJetRhovsCent", "fHistJetRhovsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt*2);
151 fHistJetRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
152 fHistJetRhovsCent->GetYaxis()->SetTitle("Jet #rho (GeV/c)");
153 fOutput->Add(fHistJetRhovsCent);
155 fHistNjetvsCent = new TH2F("fHistNjetvsCent", "fHistNjetvsCent", 101, -1, 100, 150, -0.5, 149.5);
156 fHistNjetvsCent->GetXaxis()->SetTitle("Centrality (%)");
157 fHistNjetvsCent->GetYaxis()->SetTitle("No. of jets");
158 fOutput->Add(fHistNjetvsCent);
161 if (fParticleCollArray.GetEntriesFast()>0) {
162 fHistJetPtvsNtrack = new TH2F("fHistJetPtvsNtrack", "fHistJetPtvsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt);
163 fHistJetPtvsNtrack->GetXaxis()->SetTitle("No. of tracks");
164 fHistJetPtvsNtrack->GetYaxis()->SetTitle("#it{p}_{T,jet} (GeV/c)");
165 fOutput->Add(fHistJetPtvsNtrack);
167 fHistJetAreavsNtrack = new TH2F("fHistJetAreavsNtrack", "fHistJetAreavsNtrack", 150, 0, 6000, 100, 0, 1);
168 fHistJetAreavsNtrack->GetXaxis()->SetTitle("No. of tracks");
169 fHistJetAreavsNtrack->GetYaxis()->SetTitle("Jet area");
170 fOutput->Add(fHistJetAreavsNtrack);
172 fHistNjetvsNtrack = new TH2F("fHistNjetvsNtrack", "fHistNjetvsNtrack", 150, 0, 6000, 150, -0.5, 149.5);
173 fHistNjetvsNtrack->GetXaxis()->SetTitle("No. of jets");
174 fHistNjetvsNtrack->GetYaxis()->SetTitle("No. of tracks");
175 fOutput->Add(fHistNjetvsNtrack);
180 for (Int_t i = 0; i < 4; i++) {
181 name = Form("fHistJetNconstVsPt_%d",i);
182 fHistJetNconstVsPt[i] = new TH2F(name, name, 150, -0.5, 149.5, fNbins, fMinBinPt, fMaxBinPt);
183 fHistJetNconstVsPt[i]->GetXaxis()->SetTitle("No. of constituents");
184 fHistJetNconstVsPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
185 fOutput->Add(fHistJetNconstVsPt[i]);
187 name = Form("fHistJetRhovsEta_%d",i);
188 fHistJetRhovsEta[i] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt*2, 16, -0.8, 0.8);
189 fHistJetRhovsEta[i]->GetXaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
190 fHistJetRhovsEta[i]->GetYaxis()->SetTitle("#eta");
191 fOutput->Add(fHistJetRhovsEta[i]);
193 for (Int_t j = 0; j < 3; j++) {
194 name = Form("NjUEoverNjVsNj_%d_%d",i,j+1);
195 fHistNjUEoverNjVsNj[i*3+j] = new TH2F(name, name, 150, -0.5, 149.5, 120, 0.01, 1.21);
196 fHistNjUEoverNjVsNj[i*3+j]->GetXaxis()->SetTitle("N_{jet}");
197 fHistNjUEoverNjVsNj[i*3+j]->GetYaxis()->SetTitle("N_{jet_{UE}} / N_{jet}");
198 fOutput->Add(fHistNjUEoverNjVsNj[i*3+j]);
203 if (!fCompareRhoName.IsNull()) {
204 fHistDeltaRhovsCent = new TH2F("fHistDeltaRhovsCent", "fHistDeltaRhovsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
205 fHistDeltaRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
206 fHistDeltaRhovsCent->GetYaxis()->SetTitle("#Delta#rho (GeV/c * rad^{-1})");
207 fOutput->Add(fHistDeltaRhovsCent);
209 if (fParticleCollArray.GetEntriesFast()>0) {
210 fHistDeltaRhovsNtrack = new TH2F("fHistDeltaRhovsNtrack", "fHistDeltaRhovsNtrack", 150, 0, 6000, fNbins, -fMaxBinPt, fMaxBinPt);
211 fHistDeltaRhovsNtrack->GetXaxis()->SetTitle("No. of tracks");
212 fHistDeltaRhovsNtrack->GetYaxis()->SetTitle("#Delta#rho (GeV/c * rad^{-1})");
213 fOutput->Add(fHistDeltaRhovsNtrack);
217 if (fScaleFunction) {
218 fHistRhoScaledvsCent = new TH2F("fHistRhoScaledvsCent", "fHistRhoScaledvsCent", 101, -1, 100, fNbins, fMinBinPt , fMaxBinPt*2);
219 fHistRhoScaledvsCent->GetXaxis()->SetTitle("Centrality (%)");
220 fHistRhoScaledvsCent->GetYaxis()->SetTitle("#rho_{scaled} (GeV/c * rad^{-1})");
221 fOutput->Add(fHistRhoScaledvsCent);
223 if (fParticleCollArray.GetEntriesFast()>0) {
224 fHistRhoScaledvsNtrack = new TH2F("fHistRhoScaledvsNtrack", "fHistRhoScaledvsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt*2);
225 fHistRhoScaledvsNtrack->GetXaxis()->SetTitle("No. of tracks");
226 fHistRhoScaledvsNtrack->GetYaxis()->SetTitle("#rho_{scaled} (GeV/c * rad^{-1})");
227 fOutput->Add(fHistRhoScaledvsNtrack);
230 if (fClusterCollArray.GetEntriesFast()>0) {
231 fHistRhoScaledvsNcluster = new TH2F("fHistRhoScaledvsNcluster", "fHistRhoScaledvsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
232 fHistRhoScaledvsNcluster->GetXaxis()->SetTitle("No. of clusters");
233 fHistRhoScaledvsNcluster->GetYaxis()->SetTitle("#rho_{scaled} (GeV/c * rad^{-1})");
234 fOutput->Add(fHistRhoScaledvsNcluster);
237 if (!fCompareRhoScaledName.IsNull()) {
238 fHistDeltaRhoScalevsCent = new TH2F("fHistDeltaRhoScalevsCent", "fHistDeltaRhoScalevsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
239 fHistDeltaRhoScalevsCent->GetXaxis()->SetTitle("Centrality (%)");
240 fHistDeltaRhoScalevsCent->GetYaxis()->SetTitle("#Delta#rho_{scaled} (GeV/c * rad^{-1})");
241 fOutput->Add(fHistDeltaRhoScalevsCent);
243 if (fParticleCollArray.GetEntriesFast()>0) {
244 fHistDeltaRhoScalevsNtrack = new TH2F("fHistDeltaRhoScalevsNtrack", "fHistDeltaRhoScalevsNtrack", 150, 0, 6000, fNbins, -fMaxBinPt, fMaxBinPt);
245 fHistDeltaRhoScalevsNtrack->GetXaxis()->SetTitle("No. of tracks");
246 fHistDeltaRhoScalevsNtrack->GetYaxis()->SetTitle("#Delta#rho_{scaled} (GeV/c * rad^{-1})");
247 fOutput->Add(fHistDeltaRhoScalevsNtrack);
253 //________________________________________________________________________
254 Bool_t AliAnalysisTaskRhoBase::Run()
258 Double_t rho = GetRhoFactor(fCent);
259 fOutRho->SetVal(rho);
261 if (fScaleFunction) {
262 Double_t rhoScaled = rho * GetScaleFactor(fCent);
263 fOutRhoScaled->SetVal(rhoScaled);
269 //________________________________________________________________________
270 Bool_t AliAnalysisTaskRhoBase::FillHistograms()
278 Ntracks = fTracks->GetEntries();
280 Nclusters = fCaloClusters->GetEntries();
283 Int_t Njets = fJets->GetEntries();
285 Int_t NjetUE1Sigma = 0;
286 Int_t NjetUE2Sigma = 0;
287 Int_t NjetUE3Sigma = 0;
288 Double_t rhoPlus1Sigma = fOutRho->GetVal() + fInEventSigmaRho;
289 Double_t rhoPlus2Sigma = fOutRho->GetVal() + 2*fInEventSigmaRho;
290 Double_t rhoPlus3Sigma = fOutRho->GetVal() + 3*fInEventSigmaRho;
292 for (Int_t i = 0; i < Njets; ++i) {
294 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(i));
296 AliError(Form("%s: Could not receive jet %d", GetName(), i));
303 fHistJetPtvsCent->Fill(fCent, jet->Pt());
304 fHistJetAreavsCent->Fill(fCent, jet->Area());
305 fHistJetRhovsCent->Fill(fCent, jet->Pt() / jet->Area());
306 fHistJetRhovsEta[fCentBin]->Fill(jet->Pt() / jet->Area(), jet->Eta());
309 fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
310 fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
313 fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt());
315 if (jet->Pt() < rhoPlus1Sigma * jet->Area())
318 if (jet->Pt() < rhoPlus2Sigma * jet->Area())
321 if (jet->Pt() < rhoPlus3Sigma * jet->Area())
328 fHistNjUEoverNjVsNj[fCentBin*3 ]->Fill(NjetAcc,1.*NjetUE1Sigma/NjetAcc);
329 fHistNjUEoverNjVsNj[fCentBin*3+1]->Fill(NjetAcc,1.*NjetUE2Sigma/NjetAcc);
330 fHistNjUEoverNjVsNj[fCentBin*3+2]->Fill(NjetAcc,1.*NjetUE3Sigma/NjetAcc);
333 fHistNjetvsCent->Fill(fCent, NjetAcc);
335 fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
338 fHistRhovsCent->Fill(fCent, fOutRho->GetVal());
341 fHistRhovsNtrack->Fill(Ntracks, fOutRho->GetVal());
343 fHistRhovsNcluster->Fill(Nclusters, fOutRho->GetVal());
345 fHistDeltaRhovsCent->Fill(fCent, fOutRho->GetVal() - fCompareRho->GetVal());
347 fHistDeltaRhovsNtrack->Fill(Ntracks, fOutRho->GetVal() - fCompareRho->GetVal());
351 fHistRhoScaledvsCent->Fill(fCent, fOutRhoScaled->GetVal());
353 fHistRhoScaledvsNtrack->Fill(Ntracks, fOutRhoScaled->GetVal());
355 fHistRhoScaledvsNcluster->Fill(Nclusters, fOutRhoScaled->GetVal());
356 if (fCompareRhoScaled) {
357 fHistDeltaRhoScalevsCent->Fill(fCent, fOutRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
359 fHistDeltaRhoScalevsNtrack->Fill(Ntracks, fOutRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
367 //________________________________________________________________________
368 void AliAnalysisTaskRhoBase::ExecOnce()
370 // Init the analysis.
373 fOutRho = new AliRhoParameter(fOutRhoName, 0);
375 if (fAttachToEvent) {
376 if (!(InputEvent()->FindListObject(fOutRhoName))) {
377 InputEvent()->AddObject(fOutRho);
379 AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fOutRhoName.Data()));
385 if (fScaleFunction && !fOutRhoScaled) {
386 fOutRhoScaled = new AliRhoParameter(fOutRhoScaledName, 0);
388 if (fAttachToEvent) {
389 if (!(InputEvent()->FindListObject(fOutRhoScaledName))) {
390 InputEvent()->AddObject(fOutRhoScaled);
392 AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fOutRhoScaledName.Data()));
398 if (!fCompareRhoName.IsNull() && !fCompareRho) {
399 fCompareRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoName));
401 AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoName.Data()));
405 if (!fCompareRhoScaledName.IsNull() && !fCompareRhoScaled) {
406 fCompareRhoScaled = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoScaledName));
407 if (!fCompareRhoScaled) {
408 AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoScaledName.Data()));
412 AliAnalysisTaskEmcalJet::ExecOnce();
415 //________________________________________________________________________
416 Double_t AliAnalysisTaskRhoBase::GetRhoFactor(Double_t cent)
418 // Return rho per centrality.
422 rho = fRhoFunction->Eval(cent);
426 //________________________________________________________________________
427 Double_t AliAnalysisTaskRhoBase::GetScaleFactor(Double_t cent)
433 scale = fScaleFunction->Eval(cent);