2855b7e8521e2478e075d8f926c25403732b2fae
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliAnalysisTaskRhoBase.cxx
1 // $Id$
2 //
3 // Base class for rho calculation.
4 // Calculates parameterized rho for given centrality independent of input.
5 //
6 // Author: S.Aiola
7
8 #include <TF1.h>
9 #include <TH1F.h>
10 #include <TH2F.h>
11 #include <TClonesArray.h>
12
13 #include "AliLog.h"
14 #include "AliRhoParameter.h"
15 #include "AliEmcalJet.h"
16
17 #include "AliAnalysisTaskRhoBase.h"
18
19 ClassImp(AliAnalysisTaskRhoBase)
20
21 //________________________________________________________________________
22 AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase() : 
23   AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoBase", kFALSE),
24   fRhoScaledName(),
25   fCompareRhoName(),
26   fCompareRhoScaledName(),
27   fRhoFunction(0),
28   fScaleFunction(0),
29   fInEventSigmaRho(35.83),
30   fAttachToEvent(kTRUE),
31   fRhoScaled(0),
32   fCompareRho(0),
33   fCompareRhoScaled(0),
34   fHistJetPtvsCent(0),
35   fHistJetAreavsCent(0),
36   fHistNjetvsCent(0),
37   fHistJetPtvsNtrack(0),
38   fHistJetAreavsNtrack(0),
39   fHistNjetvsNtrack(0),
40   fHistRhovsCent(0),
41   fHistRhoScaledvsCent(0),
42   fHistDeltaRhovsCent(0),
43   fHistDeltaRhoScalevsCent(0),
44   fHistRhovsNtrack(0),
45   fHistRhoScaledvsNtrack(0),
46   fHistDeltaRhovsNtrack(0),
47   fHistDeltaRhoScalevsNtrack(0),
48   fHistRhovsNcluster(0),
49   fHistRhoScaledvsNcluster(0)
50 {
51   // Constructor.
52
53   for (Int_t i = 0; i < 4; i++) {
54     fHistJetNconstVsPt[i] = 0;
55   }
56   for (Int_t i = 0; i < 12; i++) {
57     fHistNjUEoverNjVsNj[i] = 0;
58   }
59 }
60
61 //________________________________________________________________________
62 AliAnalysisTaskRhoBase::AliAnalysisTaskRhoBase(const char *name, Bool_t histo) :
63   AliAnalysisTaskEmcalJet(name, histo),
64   fRhoScaledName(),
65   fCompareRhoName(),
66   fCompareRhoScaledName(),
67   fRhoFunction(0),
68   fScaleFunction(0),
69   fInEventSigmaRho(35.83),
70   fAttachToEvent(kTRUE),
71   fRhoScaled(0),
72   fCompareRho(0),
73   fCompareRhoScaled(0),
74   fHistJetPtvsCent(0),
75   fHistJetAreavsCent(0),
76   fHistNjetvsCent(0),
77   fHistJetPtvsNtrack(0),
78   fHistJetAreavsNtrack(0),
79   fHistNjetvsNtrack(0),
80   fHistRhovsCent(0),
81   fHistRhoScaledvsCent(0),
82   fHistDeltaRhovsCent(0),
83   fHistDeltaRhoScalevsCent(0),
84   fHistRhovsNtrack(0),
85   fHistRhoScaledvsNtrack(0),
86   fHistDeltaRhovsNtrack(0),
87   fHistDeltaRhoScalevsNtrack(0),
88   fHistRhovsNcluster(0),
89   fHistRhoScaledvsNcluster(0)
90 {
91   // Constructor.
92
93   for (Int_t i = 0; i < 4; i++) {
94     fHistJetNconstVsPt[i] = 0;
95   }
96   for (Int_t i = 0; i < 12; i++) {
97     fHistNjUEoverNjVsNj[i] = 0;
98   }
99
100   SetMakeGeneralHistograms(histo);
101 }
102
103 //________________________________________________________________________
104 void AliAnalysisTaskRhoBase::UserCreateOutputObjects()
105 {
106   // User create output objects, called at the beginning of the analysis.
107
108   if (!fCreateHisto)
109     return;
110
111   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
112
113   fHistRhovsCent = new TH2F("RhovsCent", "RhovsCent", 101, -1,  100, fNbins, fMinBinPt, fMaxBinPt*2);
114   fOutput->Add(fHistRhovsCent);
115
116   if (!fTracksName.IsNull()) {
117     fHistRhovsNtrack = new TH2F("RhovsNtrack", "RhovsNtrack", 125, 0, 4000, fNbins, fMinBinPt, fMaxBinPt*2);
118     fOutput->Add(fHistRhovsNtrack);
119   }
120
121   if (!fCaloName.IsNull()) {
122     fHistRhovsNcluster = new TH2F("RhovsNcluster", "RhovsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
123     fOutput->Add(fHistRhovsNcluster);
124   }
125
126   if (!fJetsName.IsNull()) {
127     fHistJetPtvsCent            = new TH2F("JetPtvsCent",           "JetPtvsCent",           101, -1,  100,   fNbins, fMinBinPt, fMaxBinPt);
128     fHistJetAreavsCent          = new TH2F("JetAreavsCent",         "JetAreavsCent",         101, -1,  100,   30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
129     fHistNjetvsCent             = new TH2F("NjetvsCent",            "NjetvsCent",            101, -1,  100,   150, -0.5, 149.5);
130
131     fOutput->Add(fHistJetPtvsCent);
132     fOutput->Add(fHistJetAreavsCent);
133     fOutput->Add(fHistNjetvsCent);
134
135     if (!fTracksName.IsNull()) {
136       fHistJetPtvsNtrack        = new TH2F("JetPtvsNtrack",         "JetPtvsNtrack",         125,  0,  4000,  fNbins, fMinBinPt, fMaxBinPt);
137       fHistJetAreavsNtrack      = new TH2F("JetAreavsNtrack",       "JetAreavsNtrack",       125,  0,  4000,  30, 0, fJetRadius * fJetRadius * TMath::Pi() * 3);
138       fHistNjetvsNtrack         = new TH2F("NjetvsNtrack",          "rNjetvsNtrack",         125,  0,  4000,  150, -0.5, 149.5);
139
140       fOutput->Add(fHistJetPtvsNtrack);
141       fOutput->Add(fHistJetAreavsNtrack);
142       fOutput->Add(fHistNjetvsNtrack);
143     }
144
145
146     TString name;
147     for (Int_t i = 0; i < 4; i++) {
148       name = Form("fHistJetNconstVsPt_%d",i);
149       fHistJetNconstVsPt[i] = new TH2F(name, name, 150, -0.5, 149.5, fNbins, fMinBinPt, fMaxBinPt);
150       fHistJetNconstVsPt[i]->GetXaxis()->SetTitle("# of constituents");
151       fHistJetNconstVsPt[i]->GetYaxis()->SetTitle("p_{T,jet} (GeV/c)");
152       fOutput->Add(fHistJetNconstVsPt[i]);
153
154       for (Int_t j = 0; j < 3; j++) {
155         name = Form("NjUEoverNjVsNj_%d_%d",i,j+1);
156         fHistNjUEoverNjVsNj[i*3+j] = new TH2F(name, name, 150, -0.5, 149.5, 120, 0.01, 1.21);
157         fHistNjUEoverNjVsNj[i*3+j]->GetXaxis()->SetTitle("N_{jet}");
158         fHistNjUEoverNjVsNj[i*3+j]->GetYaxis()->SetTitle("N_{jet_{UE}} / N_{jet}");
159         fOutput->Add(fHistNjUEoverNjVsNj[i*3+j]);
160       }
161     }
162   }
163   
164   if (!fCompareRhoName.IsNull()) {
165     fHistDeltaRhovsCent = new TH2F("DeltaRhovsCent", "DetlaRhovsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
166     fOutput->Add(fHistDeltaRhovsCent);
167     if (!fTracksName.IsNull()) {
168       fHistDeltaRhovsNtrack = new TH2F("DeltaRhovsNtrack", "DeltaRhovsNtrack", 125, 0, 4000, fNbins, -fMaxBinPt, fMaxBinPt);
169       fOutput->Add(fHistDeltaRhovsNtrack);
170     }
171   }
172
173   if (fScaleFunction) {
174     fHistRhoScaledvsCent = new TH2F("RhoScaledvsCent", "RhoScalevsCent", 101, -1, 100, fNbins, fMinBinPt , fMaxBinPt*2);
175     fOutput->Add(fHistRhoScaledvsCent);
176
177     if (!fTracksName.IsNull()) {
178       fHistRhoScaledvsNtrack = new TH2F("RhoScaledvsNtrack", "RhoScaledvsNtrack", 125, 0, 4000, fNbins, fMinBinPt, fMaxBinPt*2);
179       fOutput->Add(fHistRhoScaledvsNtrack);
180     }
181
182     if (!fCaloName.IsNull()) {
183       fHistRhoScaledvsNcluster = new TH2F("RhoScaledvsNcluster", "RhoScaledvsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
184       fOutput->Add(fHistRhoScaledvsNcluster);
185     }
186
187     if (!fCompareRhoScaledName.IsNull()) {
188       fHistDeltaRhoScalevsCent = new TH2F("DeltaRhoScalevsCent", "DeltaRhoScalevsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
189       fOutput->Add(fHistDeltaRhoScalevsCent);
190       
191       if (!fTracksName.IsNull()) {
192         fHistDeltaRhoScalevsNtrack = new TH2F("DeltaRhoScalevsNtrack", "DeltaRhoScalevsNtrack", 125, 0, 4000, fNbins, -fMaxBinPt, fMaxBinPt);
193         fOutput->Add(fHistDeltaRhoScalevsNtrack);
194       }
195     }
196   }
197 }
198
199 //________________________________________________________________________
200 Bool_t AliAnalysisTaskRhoBase::Run() 
201 {
202   // Run the analysis.
203
204   Double_t rho = GetRhoFactor(fCent);
205   fRho->SetVal(rho);
206
207   if (fScaleFunction) {
208     Double_t rhoScaled = rho * GetScaleFactor(fCent);
209     fRhoScaled->SetVal(rhoScaled);
210   }
211
212   return kTRUE;
213 }
214
215 //________________________________________________________________________
216 Bool_t AliAnalysisTaskRhoBase::FillHistograms() 
217 {
218   // Fill histograms.
219
220   Int_t Ntracks   = 0;
221   Int_t Nclusters = 0;
222
223   if (fTracks)
224     Ntracks   = fTracks->GetEntriesFast();
225   if (fCaloClusters)
226     Nclusters = fCaloClusters->GetEntriesFast();
227
228   if (fJets) {
229     Int_t    Njets         = fJets->GetEntries();
230     Int_t    NjetAcc       = 0;
231     Int_t    NjetUE1Sigma  = 0;
232     Int_t    NjetUE2Sigma  = 0;
233     Int_t    NjetUE3Sigma  = 0;
234     Double_t rhoPlus1Sigma = fRho->GetVal() + fInEventSigmaRho;
235     Double_t rhoPlus2Sigma = fRho->GetVal() + 2*fInEventSigmaRho;
236     Double_t rhoPlus3Sigma = fRho->GetVal() + 3*fInEventSigmaRho;
237
238     for (Int_t i = 0; i < Njets; ++i) {
239       
240       AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(i));
241       if (!jet) {
242         AliError(Form("%s: Could not receive jet %d", GetName(), i));
243         continue;
244       } 
245       
246       if (!AcceptJet(jet))
247         continue;
248       
249       fHistJetPtvsCent->Fill(fCent, jet->Pt());
250       fHistJetAreavsCent->Fill(fCent, jet->Area());
251       
252       if (fTracks) {
253         fHistJetPtvsNtrack->Fill(Ntracks, jet->Pt());
254         fHistJetAreavsNtrack->Fill(Ntracks, jet->Area());
255       }
256
257       fHistJetNconstVsPt[fCentBin]->Fill(jet->GetNumberOfConstituents(), jet->Pt());
258
259       if (jet->Pt() < rhoPlus1Sigma * jet->Area())
260         NjetUE1Sigma++;
261
262       if (jet->Pt() < rhoPlus2Sigma * jet->Area())
263         NjetUE2Sigma++;
264
265       if (jet->Pt() < rhoPlus3Sigma * jet->Area())
266         NjetUE3Sigma++;
267       
268       NjetAcc++;
269     }
270     
271     if (NjetAcc>0) {
272       fHistNjUEoverNjVsNj[fCentBin*3  ]->Fill(NjetAcc,1.*NjetUE1Sigma/NjetAcc);
273       fHistNjUEoverNjVsNj[fCentBin*3+1]->Fill(NjetAcc,1.*NjetUE2Sigma/NjetAcc);
274       fHistNjUEoverNjVsNj[fCentBin*3+2]->Fill(NjetAcc,1.*NjetUE3Sigma/NjetAcc);
275     }
276
277     fHistNjetvsCent->Fill(fCent, NjetAcc);
278     if (fTracks)
279       fHistNjetvsNtrack->Fill(Ntracks, NjetAcc);
280   }
281   
282   fHistRhovsCent->Fill(fCent, fRho->GetVal());
283
284   if (fTracks)
285     fHistRhovsNtrack->Fill(Ntracks, fRho->GetVal());
286   if (fCaloClusters)
287     fHistRhovsNcluster->Fill(Nclusters, fRho->GetVal());
288   if (fCompareRho) {
289     fHistDeltaRhovsCent->Fill(fCent, fRho->GetVal() - fCompareRho->GetVal());
290     if (fTracks)
291       fHistDeltaRhovsNtrack->Fill(Ntracks, fRho->GetVal() - fCompareRho->GetVal());
292   }
293
294   if (fRhoScaled) {
295     fHistRhoScaledvsCent->Fill(fCent, fRhoScaled->GetVal());
296     if (fTracks)
297       fHistRhoScaledvsNtrack->Fill(Ntracks, fRhoScaled->GetVal());
298     if (fCaloClusters)
299       fHistRhoScaledvsNcluster->Fill(Nclusters,  fRhoScaled->GetVal());
300     if (fCompareRhoScaled) {
301       fHistDeltaRhoScalevsCent->Fill(fCent, fRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
302       if (fTracks)
303         fHistDeltaRhoScalevsNtrack->Fill(Ntracks, fRhoScaled->GetVal() - fCompareRhoScaled->GetVal());
304     }
305   }
306
307   return kTRUE;
308 }      
309
310
311 //________________________________________________________________________
312 void AliAnalysisTaskRhoBase::ExecOnce() 
313 {
314   // Init the analysis.
315
316   if (!fRho) {
317     fRho = new AliRhoParameter(fRhoName, 0);
318
319     if (fAttachToEvent) {
320       if (!(InputEvent()->FindListObject(fRhoName))) {
321         InputEvent()->AddObject(fRho);
322       } else {
323         AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fRhoName.Data()));
324         return;
325       }
326     }
327     
328     if (fScaleFunction && !fRhoScaled) {
329       fRhoScaled = new AliRhoParameter(fRhoScaledName, 0);
330       if (fAttachToEvent) {
331         if (!(InputEvent()->FindListObject(fRhoScaledName))) {
332           InputEvent()->AddObject(fRhoScaled);
333         } else {
334           AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fRhoScaledName.Data()));
335           return;
336         }
337       }
338     }
339   }
340
341   if (!fCompareRhoName.IsNull() && !fCompareRho) {
342     fCompareRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoName));
343     if (!fCompareRho) {
344       AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoName.Data()));
345     }
346   }
347
348   if (!fCompareRhoScaledName.IsNull() && !fCompareRhoScaled) {
349     fCompareRhoScaled = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoScaledName));
350     if (!fCompareRhoScaled) {
351       AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoScaledName.Data()));
352     }
353   }
354
355   AliAnalysisTaskEmcalJet::ExecOnce();
356 }
357
358 //________________________________________________________________________
359 Double_t AliAnalysisTaskRhoBase::GetRhoFactor(Double_t cent)
360 {
361   // Return rho per centrality.
362
363   Double_t rho = 0;
364   if (fRhoFunction)
365     rho = fRhoFunction->Eval(cent);
366   return rho;
367 }
368
369 //________________________________________________________________________
370 Double_t AliAnalysisTaskRhoBase::GetScaleFactor(Double_t cent)
371 {
372   // Get scale factor.
373
374   Double_t scale = 1;
375   if (fScaleFunction)
376     scale = fScaleFunction->Eval(cent);
377   return scale;
378 }