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