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"
16 #include "AliParticleContainer.h"
17 #include "AliClusterContainer.h"
19 #include "AliAnalysisTaskRhoBase.h"
21 ClassImp(AliAnalysisTaskRhoBase)
23 //________________________________________________________________________
24 AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase() :
25 AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoBase", kFALSE),
29 fCompareRhoScaledName(),
32 fInEventSigmaRho(35.83),
33 fAttachToEvent(kTRUE),
39 fHistJetAreavsCent(0),
42 fHistJetPtvsNtrack(0),
43 fHistJetAreavsNtrack(0),
46 fHistRhoScaledvsCent(0),
47 fHistDeltaRhovsCent(0),
48 fHistDeltaRhoScalevsCent(0),
50 fHistRhoScaledvsNtrack(0),
51 fHistDeltaRhovsNtrack(0),
52 fHistDeltaRhoScalevsNtrack(0),
53 fHistRhovsNcluster(0),
54 fHistRhoScaledvsNcluster(0)
58 for (Int_t i = 0; i < 4; i++) {
59 fHistJetNconstVsPt[i] = 0;
60 fHistJetRhovsEta[i] = 0;
62 for (Int_t i = 0; i < 12; i++) {
63 fHistNjUEoverNjVsNj[i] = 0;
67 //________________________________________________________________________
68 AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase(const char *name, Bool_t histo) :
69 AliAnalysisTaskEmcalJet(name, histo),
73 fCompareRhoScaledName(),
76 fInEventSigmaRho(35.83),
77 fAttachToEvent(kTRUE),
83 fHistJetAreavsCent(0),
86 fHistJetPtvsNtrack(0),
87 fHistJetAreavsNtrack(0),
90 fHistRhoScaledvsCent(0),
91 fHistDeltaRhovsCent(0),
92 fHistDeltaRhoScalevsCent(0),
94 fHistRhoScaledvsNtrack(0),
95 fHistDeltaRhovsNtrack(0),
96 fHistDeltaRhoScalevsNtrack(0),
97 fHistRhovsNcluster(0),
98 fHistRhoScaledvsNcluster(0)
102 for (Int_t i = 0; i < 4; i++) {
103 fHistJetNconstVsPt[i] = 0;
104 fHistJetRhovsEta[i] = 0;
106 for (Int_t i = 0; i < 12; i++) {
107 fHistNjUEoverNjVsNj[i] = 0;
109 SetMakeGeneralHistograms(histo);
112 //________________________________________________________________________
113 void AliAnalysisTaskRhoBase::UserCreateOutputObjects()
115 // User create output objects, called at the beginning of the analysis.
120 AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
122 fHistRhovsCent = new TH2F("fHistRhovsCent", "fHistRhovsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt*2);
123 fHistRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
124 fHistRhovsCent->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
125 fOutput->Add(fHistRhovsCent);
127 if (fParticleCollArray.GetEntriesFast()>0) {
128 fHistRhovsNtrack = new TH2F("fHistRhovsNtrack", "fHistRhovsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt*2);
129 fHistRhovsNtrack->GetXaxis()->SetTitle("No. of tracks");
130 fHistRhovsNtrack->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
131 fOutput->Add(fHistRhovsNtrack);
134 if (fClusterCollArray.GetEntriesFast()>0) {
135 fHistRhovsNcluster = new TH2F("fHistRhovsNcluster", "fHistRhovsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
136 fHistRhovsNcluster->GetXaxis()->SetTitle("No. of tracks");
137 fHistRhovsNcluster->GetYaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
138 fOutput->Add(fHistRhovsNcluster);
141 if (fJetCollArray.GetEntriesFast()>0) {
142 fHistJetPtvsCent = new TH2F("fHistJetPtvsCent", "fHistJetPtvsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt);
143 fHistJetPtvsCent->GetXaxis()->SetTitle("Centrality (%)");
144 fHistJetPtvsCent->GetYaxis()->SetTitle("#it{p}_{T,jet} (GeV/c)");
145 fOutput->Add(fHistJetPtvsCent);
147 fHistJetAreavsCent = new TH2F("fHistJetAreavsCent", "fHistJetAreavsCent", 101, -1, 100, 100, 0, 1);
148 fHistJetAreavsCent->GetXaxis()->SetTitle("Centrality (%)");
149 fHistJetAreavsCent->GetYaxis()->SetTitle("Jet area");
150 fOutput->Add(fHistJetAreavsCent);
152 fHistJetRhovsCent = new TH2F("fHistJetRhovsCent", "fHistJetRhovsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt*2);
153 fHistJetRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
154 fHistJetRhovsCent->GetYaxis()->SetTitle("Jet #rho (GeV/c)");
155 fOutput->Add(fHistJetRhovsCent);
157 fHistNjetvsCent = new TH2F("fHistNjetvsCent", "fHistNjetvsCent", 101, -1, 100, 150, -0.5, 149.5);
158 fHistNjetvsCent->GetXaxis()->SetTitle("Centrality (%)");
159 fHistNjetvsCent->GetYaxis()->SetTitle("No. of jets");
160 fOutput->Add(fHistNjetvsCent);
163 if (fParticleCollArray.GetEntriesFast()>0) {
164 fHistJetPtvsNtrack = new TH2F("fHistJetPtvsNtrack", "fHistJetPtvsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt);
165 fHistJetPtvsNtrack->GetXaxis()->SetTitle("No. of tracks");
166 fHistJetPtvsNtrack->GetYaxis()->SetTitle("#it{p}_{T,jet} (GeV/c)");
167 fOutput->Add(fHistJetPtvsNtrack);
169 fHistJetAreavsNtrack = new TH2F("fHistJetAreavsNtrack", "fHistJetAreavsNtrack", 150, 0, 6000, 100, 0, 1);
170 fHistJetAreavsNtrack->GetXaxis()->SetTitle("No. of tracks");
171 fHistJetAreavsNtrack->GetYaxis()->SetTitle("Jet area");
172 fOutput->Add(fHistJetAreavsNtrack);
174 fHistNjetvsNtrack = new TH2F("fHistNjetvsNtrack", "fHistNjetvsNtrack", 150, 0, 6000, 150, -0.5, 149.5);
175 fHistNjetvsNtrack->GetXaxis()->SetTitle("No. of jets");
176 fHistNjetvsNtrack->GetYaxis()->SetTitle("No. of tracks");
177 fOutput->Add(fHistNjetvsNtrack);
182 for (Int_t i = 0; i < 4; i++) {
183 name = Form("fHistJetNconstVsPt_%d",i);
184 fHistJetNconstVsPt[i] = new TH2F(name, name, 150, -0.5, 149.5, fNbins, fMinBinPt, fMaxBinPt);
185 fHistJetNconstVsPt[i]->GetXaxis()->SetTitle("No. of constituents");
186 fHistJetNconstVsPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
187 fOutput->Add(fHistJetNconstVsPt[i]);
189 name = Form("fHistJetRhovsEta_%d",i);
190 fHistJetRhovsEta[i] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt*2, 16, -0.8, 0.8);
191 fHistJetRhovsEta[i]->GetXaxis()->SetTitle("#rho (GeV/c * rad^{-1})");
192 fHistJetRhovsEta[i]->GetYaxis()->SetTitle("#eta");
193 fOutput->Add(fHistJetRhovsEta[i]);
195 for (Int_t j = 0; j < 3; j++) {
196 name = Form("NjUEoverNjVsNj_%d_%d",i,j+1);
197 fHistNjUEoverNjVsNj[i*3+j] = new TH2F(name, name, 150, -0.5, 149.5, 120, 0.01, 1.21);
198 fHistNjUEoverNjVsNj[i*3+j]->GetXaxis()->SetTitle("N_{jet}");
199 fHistNjUEoverNjVsNj[i*3+j]->GetYaxis()->SetTitle("N_{jet_{UE}} / N_{jet}");
200 fOutput->Add(fHistNjUEoverNjVsNj[i*3+j]);
205 if (!fCompareRhoName.IsNull()) {
206 fHistDeltaRhovsCent = new TH2F("fHistDeltaRhovsCent", "fHistDeltaRhovsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
207 fHistDeltaRhovsCent->GetXaxis()->SetTitle("Centrality (%)");
208 fHistDeltaRhovsCent->GetYaxis()->SetTitle("#Delta#rho (GeV/c * rad^{-1})");
209 fOutput->Add(fHistDeltaRhovsCent);
211 if (fParticleCollArray.GetEntriesFast()>0) {
212 fHistDeltaRhovsNtrack = new TH2F("fHistDeltaRhovsNtrack", "fHistDeltaRhovsNtrack", 150, 0, 6000, fNbins, -fMaxBinPt, fMaxBinPt);
213 fHistDeltaRhovsNtrack->GetXaxis()->SetTitle("No. of tracks");
214 fHistDeltaRhovsNtrack->GetYaxis()->SetTitle("#Delta#rho (GeV/c * rad^{-1})");
215 fOutput->Add(fHistDeltaRhovsNtrack);
219 if (fScaleFunction) {
220 fHistRhoScaledvsCent = new TH2F("fHistRhoScaledvsCent", "fHistRhoScaledvsCent", 101, -1, 100, fNbins, fMinBinPt , fMaxBinPt*2);
221 fHistRhoScaledvsCent->GetXaxis()->SetTitle("Centrality (%)");
222 fHistRhoScaledvsCent->GetYaxis()->SetTitle("#rho_{scaled} (GeV/c * rad^{-1})");
223 fOutput->Add(fHistRhoScaledvsCent);
225 if (fParticleCollArray.GetEntriesFast()>0) {
226 fHistRhoScaledvsNtrack = new TH2F("fHistRhoScaledvsNtrack", "fHistRhoScaledvsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt*2);
227 fHistRhoScaledvsNtrack->GetXaxis()->SetTitle("No. of tracks");
228 fHistRhoScaledvsNtrack->GetYaxis()->SetTitle("#rho_{scaled} (GeV/c * rad^{-1})");
229 fOutput->Add(fHistRhoScaledvsNtrack);
232 if (fClusterCollArray.GetEntriesFast()>0) {
233 fHistRhoScaledvsNcluster = new TH2F("fHistRhoScaledvsNcluster", "fHistRhoScaledvsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
234 fHistRhoScaledvsNcluster->GetXaxis()->SetTitle("No. of clusters");
235 fHistRhoScaledvsNcluster->GetYaxis()->SetTitle("#rho_{scaled} (GeV/c * rad^{-1})");
236 fOutput->Add(fHistRhoScaledvsNcluster);
239 if (!fCompareRhoScaledName.IsNull()) {
240 fHistDeltaRhoScalevsCent = new TH2F("fHistDeltaRhoScalevsCent", "fHistDeltaRhoScalevsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
241 fHistDeltaRhoScalevsCent->GetXaxis()->SetTitle("Centrality (%)");
242 fHistDeltaRhoScalevsCent->GetYaxis()->SetTitle("#Delta#rho_{scaled} (GeV/c * rad^{-1})");
243 fOutput->Add(fHistDeltaRhoScalevsCent);
245 if (fParticleCollArray.GetEntriesFast()>0) {
246 fHistDeltaRhoScalevsNtrack = new TH2F("fHistDeltaRhoScalevsNtrack", "fHistDeltaRhoScalevsNtrack", 150, 0, 6000, fNbins, -fMaxBinPt, fMaxBinPt);
247 fHistDeltaRhoScalevsNtrack->GetXaxis()->SetTitle("No. of tracks");
248 fHistDeltaRhoScalevsNtrack->GetYaxis()->SetTitle("#Delta#rho_{scaled} (GeV/c * rad^{-1})");
249 fOutput->Add(fHistDeltaRhoScalevsNtrack);
255 //________________________________________________________________________
256 Bool_t AliAnalysisTaskRhoBase::Run()
260 Double_t rho = GetRhoFactor(fCent);
261 fOutRho->SetVal(rho);
263 if (fScaleFunction) {
264 Double_t rhoScaled = rho * GetScaleFactor(fCent);
265 fOutRhoScaled->SetVal(rhoScaled);
271 //________________________________________________________________________
272 Bool_t AliAnalysisTaskRhoBase::FillHistograms()
279 if (GetParticleContainer(0))
280 Ntracks = GetParticleContainer(0)->GetNAcceptedParticles();
281 if (GetClusterContainer(0))
282 Nclusters = GetClusterContainer(0)->GetNAcceptedClusters();
285 Int_t Njets = fJets->GetEntries();
287 Int_t NjetUE1Sigma = 0;
288 Int_t NjetUE2Sigma = 0;
289 Int_t NjetUE3Sigma = 0;
290 Double_t rhoPlus1Sigma = fOutRho->GetVal() + fInEventSigmaRho;
291 Double_t rhoPlus2Sigma = fOutRho->GetVal() + 2*fInEventSigmaRho;
292 Double_t rhoPlus3Sigma = fOutRho->GetVal() + 3*fInEventSigmaRho;
294 for (Int_t i = 0; i < Njets; ++i) {
296 AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(i));
298 AliError(Form("%s: Could not receive jet %d", GetName(), i));
305 fHistJetPtvsCent->Fill(fCent, jet->Pt());
306 fHistJetAreavsCent->Fill(fCent, jet->Area());
307 fHistJetRhovsCent->Fill(fCent, jet->Pt() / jet->Area());
308 fHistJetRhovsEta[fCentBin]->Fill(jet->Pt() / jet->Area(), jet->Eta());
311 fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
312 fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
315 fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt());
317 if (jet->Pt() < rhoPlus1Sigma * jet->Area())
320 if (jet->Pt() < rhoPlus2Sigma * jet->Area())
323 if (jet->Pt() < rhoPlus3Sigma * jet->Area())
330 fHistNjUEoverNjVsNj[fCentBin*3 ]->Fill(NjetAcc,1.*NjetUE1Sigma/NjetAcc);
331 fHistNjUEoverNjVsNj[fCentBin*3+1]->Fill(NjetAcc,1.*NjetUE2Sigma/NjetAcc);
332 fHistNjUEoverNjVsNj[fCentBin*3+2]->Fill(NjetAcc,1.*NjetUE3Sigma/NjetAcc);
335 fHistNjetvsCent->Fill(fCent, NjetAcc);
337 fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
340 fHistRhovsCent->Fill(fCent, fOutRho->GetVal());
343 fHistRhovsNtrack->Fill(Ntracks, fOutRho->GetVal());
345 fHistRhovsNcluster->Fill(Nclusters, fOutRho->GetVal());
347 fHistDeltaRhovsCent->Fill(fCent, fOutRho->GetVal() - fCompareRho->GetVal());
349 fHistDeltaRhovsNtrack->Fill(Ntracks, fOutRho->GetVal() - fCompareRho->GetVal());
353 fHistRhoScaledvsCent->Fill(fCent, fOutRhoScaled->GetVal());
355 fHistRhoScaledvsNtrack->Fill(Ntracks, fOutRhoScaled->GetVal());
357 fHistRhoScaledvsNcluster->Fill(Nclusters, fOutRhoScaled->GetVal());
358 if (fCompareRhoScaled) {
359 fHistDeltaRhoScalevsCent->Fill(fCent, fOutRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
361 fHistDeltaRhoScalevsNtrack->Fill(Ntracks, fOutRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
369 //________________________________________________________________________
370 void AliAnalysisTaskRhoBase::ExecOnce()
372 // Init the analysis.
375 fOutRho = new AliRhoParameter(fOutRhoName, 0);
377 if (fAttachToEvent) {
378 if (!(InputEvent()->FindListObject(fOutRhoName))) {
379 InputEvent()->AddObject(fOutRho);
381 AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fOutRhoName.Data()));
387 if (fScaleFunction && !fOutRhoScaled) {
388 fOutRhoScaled = new AliRhoParameter(fOutRhoScaledName, 0);
390 if (fAttachToEvent) {
391 if (!(InputEvent()->FindListObject(fOutRhoScaledName))) {
392 InputEvent()->AddObject(fOutRhoScaled);
394 AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fOutRhoScaledName.Data()));
400 if (!fCompareRhoName.IsNull() && !fCompareRho) {
401 fCompareRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoName));
403 AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoName.Data()));
407 if (!fCompareRhoScaledName.IsNull() && !fCompareRhoScaled) {
408 fCompareRhoScaled = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoScaledName));
409 if (!fCompareRhoScaled) {
410 AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoScaledName.Data()));
414 AliAnalysisTaskEmcalJet::ExecOnce();
417 //________________________________________________________________________
418 Double_t AliAnalysisTaskRhoBase::GetRhoFactor(Double_t cent)
420 // Return rho per centrality.
424 rho = fRhoFunction->Eval(cent);
428 //________________________________________________________________________
429 Double_t AliAnalysisTaskRhoBase::GetScaleFactor(Double_t cent)
435 scale = fScaleFunction->Eval(cent);