Make sure that histograms are obtained from output list in
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMultiplicity.cxx
CommitLineData
7e4038b5 1#include "AliForwardMultiplicity.h"
2#include "AliTriggerAnalysis.h"
3#include "AliPhysicsSelection.h"
4#include "AliLog.h"
5#include "AliFMDAnaParameters.h"
6#include "AliESDEvent.h"
7#include "AliAODHandler.h"
8#include "AliMultiplicity.h"
9#include "AliInputEventHandler.h"
10#include <TH1.h>
11#include <TDirectory.h>
12#include <TTree.h>
13
14//====================================================================
15AliForwardMultiplicity::AliForwardMultiplicity()
16 : AliAnalysisTaskSE(),
17 fHEventsTr(0),
18 fHEventsTrVtx(0),
19 fHTriggers(0),
20 fHData(0),
21 fFirstEvent(true),
22 fLowFluxCut(1000),
23 fESDFMD(),
24 fHistos(),
25 fAODFMD(),
26 fSharingFilter(),
27 fDensityCalculator(),
28 fCorrections(),
29 fHistCollector(),
30 fList(0),
31 fTree(0)
32{
33}
34
35//____________________________________________________________________
36AliForwardMultiplicity::AliForwardMultiplicity(const char* name)
37 : AliAnalysisTaskSE(name),
38 fHEventsTr(0),
39 fHEventsTrVtx(0),
40 fHTriggers(0),
41 fHData(0),
42 fFirstEvent(true),
43 fLowFluxCut(1000),
44 fESDFMD(),
45 fHistos(),
46 fAODFMD(kTRUE),
47 fSharingFilter("sharing"),
48 fDensityCalculator("density"),
49 fCorrections("corrections"),
50 fHistCollector("collector"),
51 fList(0),
52 fTree(0)
53{
54 DefineOutput(1, TList::Class());
55 // DefineOutput(2, TTree::Class());
56}
57
58//____________________________________________________________________
59AliForwardMultiplicity::AliForwardMultiplicity(const AliForwardMultiplicity& o)
60 : AliAnalysisTaskSE(o),
61 fHEventsTr(o.fHEventsTr),
62 fHEventsTrVtx(o.fHEventsTrVtx),
63 fHTriggers(o.fHTriggers),
64 fHData(o.fHData),
65 fFirstEvent(true),
66 fLowFluxCut(1000),
67 fESDFMD(o.fESDFMD),
68 fHistos(o.fHistos),
69 fAODFMD(o.fAODFMD),
70 fSharingFilter(o.fSharingFilter),
71 fDensityCalculator(o.fDensityCalculator),
72 fCorrections(o.fCorrections),
73 fHistCollector(o.fHistCollector),
74 fList(o.fList),
75 fTree(o.fTree)
76{
77}
78
79//____________________________________________________________________
80AliForwardMultiplicity&
81AliForwardMultiplicity::operator=(const AliForwardMultiplicity& o)
82{
83 fHEventsTr = o.fHEventsTr;
84 fHEventsTrVtx = o.fHEventsTrVtx;
85 fHTriggers = o.fHTriggers;
86 fHData = o.fHData;
87 fFirstEvent = o.fFirstEvent;
88 fSharingFilter = o.fSharingFilter;
89 fDensityCalculator = o.fDensityCalculator;
90 fCorrections = o.fCorrections;
91 fHistCollector = o.fHistCollector;
92 fHistos = o.fHistos;
93 fAODFMD = o.fAODFMD;
94 fList = o.fList;
95 fTree = o.fTree;
96
97 return *this;
98}
99
100//____________________________________________________________________
101void
102AliForwardMultiplicity::Init()
103{
104 fFirstEvent = true;
105}
106
107//____________________________________________________________________
108void
109AliForwardMultiplicity::InitializeSubs()
110{
111 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
112 pars->Init(kTRUE);
113
9d99b0dd 114 fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger",
7e4038b5 115 pars->GetNvtxBins(),
116 -pars->GetVtxCutZ(),
117 pars->GetVtxCutZ());
118 fHEventsTr->SetXTitle("v_{z} [cm]");
119 fHEventsTr->SetYTitle("# of events");
120 fHEventsTr->SetFillColor(kRed+1);
121 fHEventsTr->SetFillStyle(3001);
122 fHEventsTr->SetDirectory(0);
9d99b0dd 123 fList->Add(fHEventsTr);
7e4038b5 124 // fHEventsTr->Sumw2();
125
126 fHEventsTrVtx = new TH1I("nEventsTrVtx",
127 "Number of events w/trigger and vertex",
128 pars->GetNvtxBins(),
129 -pars->GetVtxCutZ(),
130 pars->GetVtxCutZ());
131 fHEventsTrVtx->SetXTitle("v_{z} [cm]");
132 fHEventsTrVtx->SetYTitle("# of events");
133 fHEventsTrVtx->SetFillColor(kBlue+1);
134 fHEventsTrVtx->SetFillStyle(3001);
135 fHEventsTrVtx->SetDirectory(0);
9d99b0dd 136 fList->Add(fHEventsTrVtx);
7e4038b5 137 // fHEventsTrVtx->Sumw2();
138
139
140 fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
141 fHTriggers->SetFillColor(kRed+1);
142 fHTriggers->SetFillStyle(3001);
143 fHTriggers->SetStats(0);
144 fHTriggers->SetDirectory(0);
145 fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
146 fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
147 fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
148 fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
149 fHTriggers->GetXaxis()->SetBinLabel(5,"A");
150 fHTriggers->GetXaxis()->SetBinLabel(6,"B");
151 fHTriggers->GetXaxis()->SetBinLabel(7,"C");
152 fHTriggers->GetXaxis()->SetBinLabel(8,"E");
9d99b0dd 153 fList->Add(fHTriggers);
7e4038b5 154
155 TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
156 fHistos.Init(e);
157 fAODFMD.Init(e);
158
159 fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
160 fHData->SetStats(0);
161 fHData->SetDirectory(0);
9d99b0dd 162 fList->Add(fHData);
163
7e4038b5 164 fHistCollector.Init(*(fHEventsTr->GetXaxis()));
165}
166
167//____________________________________________________________________
168void
169AliForwardMultiplicity::UserCreateOutputObjects()
170{
171 fList = new TList;
172
173 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
174 AliAODHandler* ah =
175 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
176 if (!ah)
177 AliFatal("No AOD output handler set in analysis manager");
178
179
180 TObject* obj = &fAODFMD;
181 ah->AddBranch("AliAODForwardMult", &obj);
182
9d99b0dd 183 fSharingFilter.DefineOutput(fList);
184 fDensityCalculator.DefineOutput(fList);
185 fCorrections.DefineOutput(fList);
186
7e4038b5 187 // fTree = new TTree("T", "T");
188 // fTree->Branch("forward", &fAODFMD);
189
190 PostData(1, fList);
191 // PostData(2, fTree);
192}
193//____________________________________________________________________
194void
195AliForwardMultiplicity::UserExec(Option_t*)
196{
197 // Get the input data
198 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
199 if (!esd) {
200 AliWarning("No ESD event found for input event");
201 return;
202 }
203
204#if 0
205 static Int_t nEvents = 0;
206 nEvents++;
207 if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
208#endif
209
210 // On the first event, initialize the parameters
211 if (fFirstEvent) {
212 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
213 pars->SetParametersFromESD(esd);
214 pars->PrintStatus();
215 fFirstEvent = false;
216
217 InitializeSubs();
218 }
219 // Clear stuff
220 fHistos.Clear();
221 fESDFMD.Clear();
222 fAODFMD.Clear();
223
224 // Read trigger information from the ESD and store in AOD object
9d99b0dd 225 UInt_t triggers = 0;
226 if (!AliForwardUtil::ReadTriggers(esd, triggers, fHTriggers)) {
7e4038b5 227#ifdef VERBOSE
228 AliWarning("Failed to read triggers from ESD");
229#endif
230 return;
231 }
9d99b0dd 232 fAODFMD.SetTriggerBits(triggers);
7e4038b5 233
234 // Mark this event for storage
235 MarkEventForStore();
236
237 // Check if this is a high-flux event
238 const AliMultiplicity* testmult = esd->GetMultiplicity();
239 if (!testmult) {
240#ifdef VERBOSE
241 AliWarning("No central multiplicity object found");
242#endif
243 return;
244 }
245 Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
246
247 // Get the FMD ESD data
248 AliESDFMD* esdFMD = esd->GetFMDData();
249 if (!esdFMD) {
250#ifdef VERBOSE
251 AliWarning("No FMD data found in ESD");
252#endif
253 return;
254 }
255
256 // Get the vertex information
257 Double_t vz = 0;
9d99b0dd 258 Bool_t vzOk = AliForwardUtil::ReadVertex(esd, vz);
7e4038b5 259
260 fHEventsTr->Fill(vz);
261 if (!vzOk) {
262#ifdef VERBOSE
263 AliWarning("Failed to read vertex from ESD");
264#endif
265 return;
266 }
267 fHEventsTrVtx->Fill(vz);
268
269 // Get the vertex bin
270 Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
271 fAODFMD.SetIpZ(vz);
272 if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) {
273#if 0
274 AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
275 vz, fHEventsTr->GetXaxis()->GetXmin(),
276 fHEventsTr->GetXaxis()->GetXmax()));
277#endif
278 return;
279 }
280
281 // Apply the sharing filter (or hit merging or clustering if you like)
282 if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
283#ifdef VERBOSE
284 AliWarning("Sharing filter failed!");
285#endif
286 return;
287 }
288
289 // Calculate the inclusive charged particle density
290 if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
291 AliWarning("Density calculator failed!");
292 return;
293 }
294
295 // Do the secondary and other corrections.
296 if (!fCorrections.Correct(fHistos, ivz)) {
297 AliWarning("Corrections failed");
298 return;
299 }
300
301 if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
302 AliWarning("Histogram collector failed");
303 return;
304 }
305
306 if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
307 fHData->Add(&(fAODFMD.GetHistogram()));
308}
309
310//____________________________________________________________________
311void
312AliForwardMultiplicity::Terminate(Option_t*)
313{
314 TList* list = dynamic_cast<TList*>(GetOutputData(1));
315 if (!list) {
316 AliError("No output list defined");
317 return;
318 }
9d99b0dd 319
320 // Get our histograms from the container
321 TH1I* hEventsTr = static_cast<TH1I*>(list->FindObject("nEventsTr"));
322 TH1I* hEventsTrVtx = static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
323 TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
324
325
7e4038b5 326 // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
9d99b0dd 327 TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
328 TH1D* norm = hData->ProjectionX("dNdeta", 0, 1, "");
7e4038b5 329 dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
330 dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
331 dNdeta->Divide(norm);
332 dNdeta->SetStats(0);
9d99b0dd 333 dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
7e4038b5 334 "width");
7e4038b5 335 list->Add(dNdeta);
336
9d99b0dd 337 fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
338 fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
339 fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
7e4038b5 340}
341
342//____________________________________________________________________
343void
344AliForwardMultiplicity::MarkEventForStore() const
345{
346 // Make sure the AOD tree is filled
347 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
348 AliAODHandler* ah =
349 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
350 if (!ah)
351 AliFatal("No AOD output handler set in analysis manager");
352
353 ah->SetFillAOD(kTRUE);
354}
7e4038b5 355
356//____________________________________________________________________
357void
358AliForwardMultiplicity::Print(Option_t*) const
359{
360}
361
362//
363// EOF
364//