]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis2/AliForwardMCMultiplicityTask.cxx
Code update for handling centrality in the dN/deta analysis.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardMCMultiplicityTask.cxx
CommitLineData
7984e5f7 1//
2// Calculate the multiplicity in the forward regions event-by-event
3//
4// Inputs:
5// - AliESDEvent
6// - Kinematics
7// - Track references
8//
9// Outputs:
10// - AliAODForwardMult
11//
12// Histograms
13//
14// Corrections used
15//
285e7b27 16#include "AliForwardMCMultiplicityTask.h"
17#include "AliTriggerAnalysis.h"
18#include "AliPhysicsSelection.h"
19#include "AliLog.h"
285e7b27 20#include "AliESDEvent.h"
21#include "AliAODHandler.h"
22#include "AliMultiplicity.h"
23#include "AliInputEventHandler.h"
24#include "AliForwardCorrectionManager.h"
25#include "AliAnalysisManager.h"
26#include <TH1.h>
27#include <TDirectory.h>
28#include <TTree.h>
29#include <TROOT.h>
285e7b27 30
31//====================================================================
32AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
33 : AliForwardMultiplicityBase(),
34 fHData(0),
35 fESDFMD(),
36 fHistos(),
37 fAODFMD(),
38 fMCESDFMD(),
39 fMCHistos(),
40 fMCAODFMD(),
4cbdf467 41 fPrimary(0),
285e7b27 42 fEventInspector(),
43 fEnergyFitter(),
44 fSharingFilter(),
45 fDensityCalculator(),
46 fCorrections(),
47 fHistCollector(),
48 fList(0)
49{
7984e5f7 50 //
51 // Constructor
52 //
285e7b27 53}
54
55//____________________________________________________________________
56AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
57 : AliForwardMultiplicityBase(name),
58 fHData(0),
59 fESDFMD(),
60 fHistos(),
61 fAODFMD(kFALSE),
62 fMCESDFMD(),
63 fMCHistos(),
64 fMCAODFMD(kTRUE),
4cbdf467 65 fPrimary(0),
285e7b27 66 fEventInspector("event"),
67 fEnergyFitter("energy"),
68 fSharingFilter("sharing"),
69 fDensityCalculator("density"),
70 fCorrections("corrections"),
71 fHistCollector("collector"),
72 fList(0)
73{
7984e5f7 74 //
75 // Constructor
76 //
77 // Parameters:
78 // name Name of task
79 //
285e7b27 80 DefineOutput(1, TList::Class());
81}
82
83//____________________________________________________________________
84AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMultiplicityTask& o)
85 : AliForwardMultiplicityBase(o),
86 fHData(o.fHData),
87 fESDFMD(o.fESDFMD),
88 fHistos(o.fHistos),
89 fAODFMD(o.fAODFMD),
90 fMCESDFMD(o.fMCESDFMD),
91 fMCHistos(o.fMCHistos),
92 fMCAODFMD(o.fMCAODFMD),
4cbdf467 93 fPrimary(o.fPrimary),
285e7b27 94 fEventInspector(o.fEventInspector),
95 fEnergyFitter(o.fEnergyFitter),
4cbdf467 96 fSharingFilter(o.fSharingFilter),
285e7b27 97 fDensityCalculator(o.fDensityCalculator),
98 fCorrections(o.fCorrections),
99 fHistCollector(o.fHistCollector),
100 fList(o.fList)
101{
7984e5f7 102 //
103 // Copy constructor
104 //
105 // Parameters:
106 // o Object to copy from
107 //
285e7b27 108 DefineOutput(1, TList::Class());
109}
110
111//____________________________________________________________________
112AliForwardMCMultiplicityTask&
113AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
114{
7984e5f7 115 //
116 // Assignment operator
117 //
118 // Parameters:
119 // o Object to assign from
120 //
121 // Return:
122 // Reference to this object
123 //
285e7b27 124 AliForwardMultiplicityBase::operator=(o);
125
126 fHData = o.fHData;
127 fEventInspector = o.fEventInspector;
128 fEnergyFitter = o.fEnergyFitter;
129 fSharingFilter = o.fSharingFilter;
130 fDensityCalculator = o.fDensityCalculator;
131 fCorrections = o.fCorrections;
132 fHistCollector = o.fHistCollector;
133 fHistos = o.fHistos;
134 fAODFMD = o.fAODFMD;
135 fMCHistos = o.fMCHistos;
136 fMCAODFMD = o.fMCAODFMD;
4cbdf467 137 fPrimary = o.fPrimary;
285e7b27 138 fList = o.fList;
139
140 return *this;
141}
142
143//____________________________________________________________________
144void
145AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
146{
7984e5f7 147 //
148 // Set debug level
149 //
150 // Parameters:
151 // dbg debug level
152 //
285e7b27 153 fEventInspector.SetDebug(dbg);
154 fEnergyFitter.SetDebug(dbg);
155 fSharingFilter.SetDebug(dbg);
156 fDensityCalculator.SetDebug(dbg);
157 fCorrections.SetDebug(dbg);
158 fHistCollector.SetDebug(dbg);
159}
160//____________________________________________________________________
161void
162AliForwardMCMultiplicityTask::InitializeSubs()
163{
7984e5f7 164 //
165 // Initialise the sub objects and stuff. Called on first event
166 //
167 //
7ec4d843 168 const TAxis* pe = 0;
169 const TAxis* pv = 0;
1174780f 170
19abe41d 171 if (!ReadCorrections(pe,pv,true)) return;
285e7b27 172
173 fHistos.Init(*pe);
174 fAODFMD.Init(*pe);
175 fMCHistos.Init(*pe);
176 fMCAODFMD.Init(*pe);
177
178 fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
179 fHData->SetStats(0);
180 fHData->SetDirectory(0);
4cbdf467 181
285e7b27 182 fList->Add(fHData);
183
184 fEnergyFitter.Init(*pe);
185 fEventInspector.Init(*pv);
186 fDensityCalculator.Init(*pe);
187 fCorrections.Init(*pe);
12fffad7 188 fHistCollector.Init(*pv,*pe);
285e7b27 189
190 this->Print();
191}
192
193//____________________________________________________________________
194void
195AliForwardMCMultiplicityTask::UserCreateOutputObjects()
196{
7984e5f7 197 //
198 // Create output objects
199 //
200 //
285e7b27 201 fList = new TList;
e1f47419 202 fList->SetOwner();
285e7b27 203
204 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
205 AliAODHandler* ah =
206 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
207 if (!ah) AliFatal("No AOD output handler set in analysis manager");
208
209
210 TObject* obj = &fAODFMD;
211 ah->AddBranch("AliAODForwardMult", &obj);
212
213 TObject* mcobj = &fMCAODFMD;
214 ah->AddBranch("AliAODForwardMult", &mcobj);
215
4cbdf467 216 fPrimary = new TH2D("primary", "MC Primaries",
217 200, -4, 6, 20, 0, 2*TMath::Pi());
218 fPrimary->SetXTitle("#eta");
219 fPrimary->SetYTitle("#varphi [radians]");
220 fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
221 fPrimary->Sumw2();
222 fPrimary->SetStats(0);
223 fPrimary->SetDirectory(0);
224 ah->AddBranch("TH2D", &fPrimary);
225
285e7b27 226 fEventInspector.DefineOutput(fList);
227 fEnergyFitter.DefineOutput(fList);
228 fSharingFilter.DefineOutput(fList);
229 fDensityCalculator.DefineOutput(fList);
230 fCorrections.DefineOutput(fList);
12fffad7 231 fHistCollector.DefineOutput(fList);
285e7b27 232
233 PostData(1, fList);
234}
235//____________________________________________________________________
236void
237AliForwardMCMultiplicityTask::UserExec(Option_t*)
238{
7984e5f7 239 //
240 // Process each event
241 //
242 // Parameters:
243 // option Not used
244 //
245
285e7b27 246 // Get the input data
e1f47419 247 AliESDEvent* esd = GetESDEvent();
248 AliMCEvent* mcEvent = MCEvent();
285e7b27 249
285e7b27 250 // Clear stuff
251 fHistos.Clear();
252 fESDFMD.Clear();
253 fAODFMD.Clear();
254 fMCHistos.Clear();
255 fMCESDFMD.Clear();
256 fMCAODFMD.Clear();
4cbdf467 257 fPrimary->Reset();
285e7b27 258
259 Bool_t lowFlux = kFALSE;
260 UInt_t triggers = 0;
261 UShort_t ivz = 0;
262 Double_t vz = 0;
5e4d905e 263 Double_t cent = 0;
264 UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
265 ivz, vz, cent);
e1f47419 266 UShort_t ivzMC = 0;
267 Double_t vzMC = 0;
268 Double_t phiR = 0;
269 Double_t b = 0;
270 // UInt_t foundMC =
271 fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, phiR);
1dbfc345 272
273
e333578d 274 //Store all events
148c39de 275 MarkEventForStore();
e333578d 276
277 Bool_t isAccepted = true;
278 if (found & AliFMDEventInspector::kNoEvent) isAccepted = false; // return;
279 if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
148c39de 280 //MarkEventForStore();
285e7b27 281 // Set trigger bits, and mark this event for storage
282 fAODFMD.SetTriggerBits(triggers);
283 fMCAODFMD.SetTriggerBits(triggers);
e58000b7 284 fAODFMD.SetCentrality(cent);
285
e333578d 286 //All events should be stored - HHD
148c39de 287 //if (isAccepted) MarkEventForStore();
285e7b27 288
e1f47419 289 // Disable this check on SPD - will bias data
290 // if (found & AliFMDEventInspector::kNoSPD) isAccepted = false; // return;
e333578d 291 if (found & AliFMDEventInspector::kNoFMD) isAccepted = false; // return;
292 if (found & AliFMDEventInspector::kNoVertex) isAccepted = false; // return;
293
294 if (isAccepted) {
295 fAODFMD.SetIpZ(vz);
296 fMCAODFMD.SetIpZ(vz);
297 }
298 if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
285e7b27 299
300 // We we do not want to use low flux specific code, we disable it here.
301 if (!fEnableLowFlux) lowFlux = false;
148c39de 302
303
285e7b27 304
305 // Get FMD data
306 AliESDFMD* esdFMD = esd->GetFMDData();
e1f47419 307
285e7b27 308 // Apply the sharing filter (or hit merging or clustering if you like)
e333578d 309 if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) {
285e7b27 310 AliWarning("Sharing filter failed!");
311 return;
312 }
4cbdf467 313 if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) {
285e7b27 314 AliWarning("MC Sharing filter failed!");
315 return;
316 }
e333578d 317 if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
148c39de 318 // HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
1dbfc345 319
148c39de 320 //MarkEventForStore();
285e7b27 321 fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
322
323 // Do the energy stuff
5e4d905e 324 if (!fEnergyFitter.Accumulate(*esdFMD, cent,
325 triggers & AliAODForwardMult::kEmpty)){
285e7b27 326 AliWarning("Energy fitter failed");
327 return;
328 }
329
330 // Calculate the inclusive charged particle density
331 if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
332 AliWarning("Density calculator failed!");
333 return;
334 }
335 if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
336 AliWarning("MC Density calculator failed!");
337 return;
338 }
339 fDensityCalculator.CompareResults(fHistos, fMCHistos);
340
341 // Do the secondary and other corrections.
342 if (!fCorrections.Correct(fHistos, ivz)) {
343 AliWarning("Corrections failed");
344 return;
345 }
346 if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
347 AliWarning("MC Corrections failed");
348 return;
349 }
350 fCorrections.CompareResults(fHistos, fMCHistos);
351
352 if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
353 AliWarning("Histogram collector failed");
354 return;
355 }
356 if (!fHistCollector.Collect(fMCHistos, ivz, fMCAODFMD.GetHistogram())) {
357 AliWarning("MC Histogram collector failed");
358 return;
359 }
360
361 if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
362 fHData->Add(&(fAODFMD.GetHistogram()));
363
364 PostData(1, fList);
365}
366
367//____________________________________________________________________
368void
369AliForwardMCMultiplicityTask::Terminate(Option_t*)
370{
7984e5f7 371 //
372 // End of job
373 //
374 // Parameters:
375 // option Not used
376 //
285e7b27 377 TList* list = dynamic_cast<TList*>(GetOutputData(1));
378 if (!list) {
379 AliError(Form("No output list defined (%p)", GetOutputData(1)));
380 if (GetOutputData(1)) GetOutputData(1)->Print();
381 return;
382 }
383
384 // Get our histograms from the container
385 TH1I* hEventsTr = 0;
386 TH1I* hEventsTrVtx = 0;
387 TH1I* hTriggers = 0;
388 if (!fEventInspector.FetchHistograms(list, hEventsTr,
389 hEventsTrVtx, hTriggers)) {
390 AliError(Form("Didn't get histograms from event selector "
391 "(hEventsTr=%p,hEventsTrVtx=%p)",
392 hEventsTr, hEventsTrVtx));
393 list->ls();
394 return;
395 }
396
397 TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
398 if (!hData) {
399 AliError(Form("Couldn't get our summed histogram from output "
400 "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
401 list->ls();
402 return;
403 }
404
405 // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
406 TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
407 TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
408 dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
409 dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
410 dNdeta->Divide(norm);
411 dNdeta->SetStats(0);
412 dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
413 "width");
414 list->Add(dNdeta);
415 list->Add(norm);
416
417 fEnergyFitter.Fit(list);
6feacd76 418 fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
419 fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
420 fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
285e7b27 421}
422
285e7b27 423
424//
425// EOF
426//