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