2 // Calculate the multiplicity in the forward regions event-by-event
10 // - AliAODForwardMult
16 #include "AliForwardMCMultiplicityTask.h"
17 #include "AliTriggerAnalysis.h"
18 #include "AliPhysicsSelection.h"
20 #include "AliESDEvent.h"
21 #include "AliAODHandler.h"
22 #include "AliMultiplicity.h"
23 #include "AliInputEventHandler.h"
24 #include "AliForwardCorrectionManager.h"
25 #include "AliAnalysisManager.h"
27 #include <TDirectory.h>
31 //====================================================================
32 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask()
33 : AliForwardMultiplicityBase(),
58 //____________________________________________________________________
59 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const char* name)
60 : AliForwardMultiplicityBase(name),
72 fEventInspector("event"),
73 fSharingFilter("sharing"),
74 fDensityCalculator("density"),
75 fCorrections("corrections"),
76 fHistCollector("collector"),
77 fEventPlaneFinder("eventplane"),
86 DefineOutput(1, TList::Class());
89 //____________________________________________________________________
90 AliForwardMCMultiplicityTask::AliForwardMCMultiplicityTask(const AliForwardMCMultiplicityTask& o)
91 : AliForwardMultiplicityBase(o),
97 fMCESDFMD(o.fMCESDFMD),
98 fMCHistos(o.fMCHistos),
99 fMCAODFMD(o.fMCAODFMD),
100 fRingSums(o.fRingSums),
101 fMCRingSums(o.fMCRingSums),
102 fPrimary(o.fPrimary),
103 fEventInspector(o.fEventInspector),
104 fSharingFilter(o.fSharingFilter),
105 fDensityCalculator(o.fDensityCalculator),
106 fCorrections(o.fCorrections),
107 fHistCollector(o.fHistCollector),
108 fEventPlaneFinder(o.fEventPlaneFinder),
115 // o Object to copy from
117 DefineOutput(1, TList::Class());
120 //____________________________________________________________________
121 AliForwardMCMultiplicityTask&
122 AliForwardMCMultiplicityTask::operator=(const AliForwardMCMultiplicityTask& o)
125 // Assignment operator
128 // o Object to assign from
131 // Reference to this object
133 if (&o == this) return *this;
134 AliForwardMultiplicityBase::operator=(o);
137 fEventInspector = o.fEventInspector;
138 fSharingFilter = o.fSharingFilter;
139 fDensityCalculator = o.fDensityCalculator;
140 fCorrections = o.fCorrections;
141 fHistCollector = o.fHistCollector;
142 fEventPlaneFinder = o.fEventPlaneFinder;
146 fMCHistos = o.fMCHistos;
147 fMCAODFMD = o.fMCAODFMD;
148 fRingSums = o.fRingSums;
149 fMCRingSums = o.fMCRingSums;
150 fPrimary = o.fPrimary;
156 //____________________________________________________________________
158 AliForwardMCMultiplicityTask::SetDebug(Int_t dbg)
166 fEventInspector.SetDebug(dbg);
167 fSharingFilter.SetDebug(dbg);
168 fDensityCalculator.SetDebug(dbg);
169 fCorrections.SetDebug(dbg);
170 fHistCollector.SetDebug(dbg);
171 fEventPlaneFinder.SetDebug(dbg);
173 //____________________________________________________________________
175 AliForwardMCMultiplicityTask::SetOnlyPrimary(Bool_t use)
177 fSharingFilter.GetTrackDensity().SetUseOnlyPrimary(use);
178 fCorrections.SetSecondaryForMC(!use);
181 //____________________________________________________________________
183 AliForwardMCMultiplicityTask::InitializeSubs()
186 // Initialise the sub objects and stuff. Called on first event
192 if (!ReadCorrections(pe,pv,true)) return false;
200 fMCRingSums.Init(*pe);
202 fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
204 fHData->SetDirectory(0);
208 TList* rings = new TList;
209 rings->SetName("ringSums");
213 rings->Add(fRingSums.Get(1, 'I'));
214 rings->Add(fRingSums.Get(2, 'I'));
215 rings->Add(fRingSums.Get(2, 'O'));
216 rings->Add(fRingSums.Get(3, 'I'));
217 rings->Add(fRingSums.Get(3, 'O'));
218 fRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
219 fRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
220 fRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
221 fRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
222 fRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
224 TList* mcRings = new TList;
225 mcRings->SetName("mcRingSums");
229 mcRings->Add(fMCRingSums.Get(1, 'I'));
230 mcRings->Add(fMCRingSums.Get(2, 'I'));
231 mcRings->Add(fMCRingSums.Get(2, 'O'));
232 mcRings->Add(fMCRingSums.Get(3, 'I'));
233 mcRings->Add(fMCRingSums.Get(3, 'O'));
234 fMCRingSums.Get(1, 'I')->SetMarkerColor(AliForwardUtil::RingColor(1, 'I'));
235 fMCRingSums.Get(2, 'I')->SetMarkerColor(AliForwardUtil::RingColor(2, 'I'));
236 fMCRingSums.Get(2, 'O')->SetMarkerColor(AliForwardUtil::RingColor(2, 'O'));
237 fMCRingSums.Get(3, 'I')->SetMarkerColor(AliForwardUtil::RingColor(3, 'I'));
238 fMCRingSums.Get(3, 'O')->SetMarkerColor(AliForwardUtil::RingColor(3, 'O'));
242 fEventInspector.Init(*pv);
243 fSharingFilter.Init(*pe);
244 fDensityCalculator.Init(*pe);
245 fCorrections.Init(*pe);
246 fHistCollector.Init(*pv,*pe);
247 fEventPlaneFinder.Init(*pe);
254 //____________________________________________________________________
256 AliForwardMCMultiplicityTask::UserCreateOutputObjects()
259 // Create output objects
265 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
267 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
268 if (!ah) AliFatal("No AOD output handler set in analysis manager");
271 TObject* obj = &fAODFMD;
272 ah->AddBranch("AliAODForwardMult", &obj);
274 TObject* mcobj = &fMCAODFMD;
275 ah->AddBranch("AliAODForwardMult", &mcobj);
277 TObject* epobj = &fAODEP;
278 ah->AddBranch("AliAODForwardEP", &epobj);
280 fPrimary = new TH2D("primary", "MC Primaries",
281 200, -4, 6, 20, 0, 2*TMath::Pi());
282 fPrimary->SetXTitle("#eta");
283 fPrimary->SetYTitle("#varphi [radians]");
284 fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
286 fPrimary->SetStats(0);
287 fPrimary->SetDirectory(0);
288 ah->AddBranch("TH2D", &fPrimary);
290 fEventInspector.DefineOutput(fList);
291 fSharingFilter.DefineOutput(fList);
292 fDensityCalculator.DefineOutput(fList);
293 fCorrections.DefineOutput(fList);
294 fHistCollector.DefineOutput(fList);
295 fEventPlaneFinder.DefineOutput(fList);
299 //____________________________________________________________________
301 AliForwardMCMultiplicityTask::UserExec(Option_t*)
304 // Process each event
310 // Read production details
312 fEventInspector.ReadProductionDetails(MCEvent());
314 // Get the input data
315 AliESDEvent* esd = GetESDEvent();
316 AliMCEvent* mcEvent = MCEvent();
317 if (!esd || !mcEvent) return;
329 Bool_t lowFlux = kFALSE;
334 UShort_t nClusters = 0;
335 UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
336 ivz, vz, cent, nClusters);
344 fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b,
346 fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
351 Bool_t isAccepted = true;
352 if (found & AliFMDEventInspector::kNoEvent) isAccepted = false; // return;
353 if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
354 //MarkEventForStore();
355 // Always set the B trigger - each MC event _is_ a collision
356 triggers |= AliAODForwardMult::kB;
357 // Set trigger bits, and mark this event for storage
358 fAODFMD.SetTriggerBits(triggers);
359 fAODFMD.SetSNN(fEventInspector.GetEnergy());
360 fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
361 fAODFMD.SetCentrality(cent);
362 fAODFMD.SetNClusters(nClusters);
364 fMCAODFMD.SetTriggerBits(triggers);
365 fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
366 fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
367 fMCAODFMD.SetCentrality(cent);
368 fMCAODFMD.SetNClusters(nClusters);
370 //All events should be stored - HHD
371 //if (isAccepted) MarkEventForStore();
373 // Disable this check on SPD - will bias data
374 // if (found & AliFMDEventInspector::kNoSPD) isAccepted = false; // return;
375 if (found & AliFMDEventInspector::kNoFMD) isAccepted = false; // return;
376 if (found & AliFMDEventInspector::kNoVertex) isAccepted = false; // return;
380 fMCAODFMD.SetIpZ(vz);
382 if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
384 // We we do not want to use low flux specific code, we disable it here.
385 if (!fEnableLowFlux) lowFlux = false;
390 AliESDFMD* esdFMD = esd->GetFMDData();
392 // Apply the sharing filter (or hit merging or clustering if you like)
393 if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
394 AliWarning("Sharing filter failed!");
397 if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) {
398 AliWarning("MC Sharing filter failed!");
401 if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
402 // HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
404 //MarkEventForStore();
405 fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
407 // Calculate the inclusive charged particle density
408 if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent, vz)) {
409 AliWarning("Density calculator failed!");
412 if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
413 AliWarning("MC Density calculator failed!");
416 fDensityCalculator.CompareResults(fHistos, fMCHistos);
418 if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
419 if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos))
420 AliWarning("Eventplane finder failed!");
423 // Do the secondary and other corrections.
424 if (!fCorrections.Correct(fHistos, ivz)) {
425 AliWarning("Corrections failed");
428 if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
429 AliWarning("MC Corrections failed");
432 fCorrections.CompareResults(fHistos, fMCHistos);
434 if (!fHistCollector.Collect(fHistos, fRingSums,
435 ivz, fAODFMD.GetHistogram())) {
436 AliWarning("Histogram collector failed");
439 if (!fHistCollector.Collect(fMCHistos, fMCRingSums,
440 ivz, fMCAODFMD.GetHistogram())) {
441 AliWarning("MC Histogram collector failed");
445 if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
446 fHData->Add(&(fAODFMD.GetHistogram()));
451 //____________________________________________________________________
453 AliForwardMCMultiplicityTask::Terminate(Option_t*)
461 TList* list = dynamic_cast<TList*>(GetOutputData(1));
463 AliError(Form("No output list defined (%p)", GetOutputData(1)));
464 if (GetOutputData(1)) GetOutputData(1)->Print();
468 // Get our histograms from the container
470 TH1I* hEventsTrVtx = 0;
472 if (!fEventInspector.FetchHistograms(list, hEventsTr,
473 hEventsTrVtx, hTriggers)) {
474 AliError(Form("Didn't get histograms from event selector "
475 "(hEventsTr=%p,hEventsTrVtx=%p)",
476 hEventsTr, hEventsTrVtx));
481 TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
483 AliError(Form("Couldn't get our summed histogram from output "
484 "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
489 // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
490 TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
491 TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
492 dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
493 dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
494 dNdeta->Divide(norm);
496 dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
501 MakeRingdNdeta(list, "ringSums", list, "ringResults");
502 MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24);
504 fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
505 fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
506 fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));