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'));
240 fEventInspector.Init(*pv);
241 fSharingFilter.Init(*pe);
242 fDensityCalculator.Init(*pe);
243 fCorrections.Init(*pe);
244 fHistCollector.Init(*pv,*pe);
245 fEventPlaneFinder.Init(*pe);
252 //____________________________________________________________________
254 AliForwardMCMultiplicityTask::UserCreateOutputObjects()
257 // Create output objects
263 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
265 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
266 if (!ah) AliFatal("No AOD output handler set in analysis manager");
269 TObject* obj = &fAODFMD;
270 ah->AddBranch("AliAODForwardMult", &obj);
272 TObject* mcobj = &fMCAODFMD;
273 ah->AddBranch("AliAODForwardMult", &mcobj);
275 TObject* epobj = &fAODEP;
276 ah->AddBranch("AliAODForwardEP", &epobj);
278 fPrimary = new TH2D("primary", "MC Primaries",
279 200, -4, 6, 20, 0, 2*TMath::Pi());
280 fPrimary->SetXTitle("#eta");
281 fPrimary->SetYTitle("#varphi [radians]");
282 fPrimary->SetZTitle("d^{2}N_{ch}/d#etad#phi");
284 fPrimary->SetStats(0);
285 fPrimary->SetDirectory(0);
286 ah->AddBranch("TH2D", &fPrimary);
288 fEventInspector.DefineOutput(fList);
289 fSharingFilter.DefineOutput(fList);
290 fDensityCalculator.DefineOutput(fList);
291 fCorrections.DefineOutput(fList);
292 fHistCollector.DefineOutput(fList);
293 fEventPlaneFinder.DefineOutput(fList);
297 //____________________________________________________________________
299 AliForwardMCMultiplicityTask::UserExec(Option_t*)
302 // Process each event
308 // Read production details
310 fEventInspector.ReadProductionDetails(MCEvent());
312 // Get the input data
313 AliESDEvent* esd = GetESDEvent();
314 AliMCEvent* mcEvent = MCEvent();
315 if (!esd || !mcEvent) return;
327 Bool_t lowFlux = kFALSE;
332 UShort_t nClusters = 0;
333 UInt_t found = fEventInspector.Process(esd, triggers, lowFlux,
334 ivz, vz, cent, nClusters);
342 fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b,
344 fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
349 Bool_t isAccepted = true;
350 if (found & AliFMDEventInspector::kNoEvent) isAccepted = false; // return;
351 if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
352 //MarkEventForStore();
353 // Always set the B trigger - each MC event _is_ a collision
354 triggers |= AliAODForwardMult::kB;
355 // Set trigger bits, and mark this event for storage
356 fAODFMD.SetTriggerBits(triggers);
357 fAODFMD.SetSNN(fEventInspector.GetEnergy());
358 fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
359 fAODFMD.SetCentrality(cent);
360 fAODFMD.SetNClusters(nClusters);
362 fMCAODFMD.SetTriggerBits(triggers);
363 fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
364 fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
365 fMCAODFMD.SetCentrality(cent);
366 fMCAODFMD.SetNClusters(nClusters);
368 //All events should be stored - HHD
369 //if (isAccepted) MarkEventForStore();
371 // Disable this check on SPD - will bias data
372 // if (found & AliFMDEventInspector::kNoSPD) isAccepted = false; // return;
373 if (found & AliFMDEventInspector::kNoFMD) isAccepted = false; // return;
374 if (found & AliFMDEventInspector::kNoVertex) isAccepted = false; // return;
378 fMCAODFMD.SetIpZ(vz);
380 if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
382 // We we do not want to use low flux specific code, we disable it here.
383 if (!fEnableLowFlux) lowFlux = false;
388 AliESDFMD* esdFMD = esd->GetFMDData();
390 // Apply the sharing filter (or hit merging or clustering if you like)
391 if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
392 AliWarning("Sharing filter failed!");
395 if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) {
396 AliWarning("MC Sharing filter failed!");
399 if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
400 // HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
402 //MarkEventForStore();
403 fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
405 // Calculate the inclusive charged particle density
406 if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent, vz)) {
407 AliWarning("Density calculator failed!");
410 if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
411 AliWarning("MC Density calculator failed!");
414 fDensityCalculator.CompareResults(fHistos, fMCHistos);
416 if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
417 if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos))
418 AliWarning("Eventplane finder failed!");
421 // Do the secondary and other corrections.
422 if (!fCorrections.Correct(fHistos, ivz)) {
423 AliWarning("Corrections failed");
426 if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
427 AliWarning("MC Corrections failed");
430 fCorrections.CompareResults(fHistos, fMCHistos);
432 if (!fHistCollector.Collect(fHistos, fRingSums,
433 ivz, fAODFMD.GetHistogram())) {
434 AliWarning("Histogram collector failed");
437 if (!fHistCollector.Collect(fMCHistos, fMCRingSums,
438 ivz, fMCAODFMD.GetHistogram())) {
439 AliWarning("MC Histogram collector failed");
443 if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
444 fHData->Add(&(fAODFMD.GetHistogram()));
449 //____________________________________________________________________
451 AliForwardMCMultiplicityTask::Terminate(Option_t*)
459 TList* list = dynamic_cast<TList*>(GetOutputData(1));
461 AliError(Form("No output list defined (%p)", GetOutputData(1)));
462 if (GetOutputData(1)) GetOutputData(1)->Print();
466 // Get our histograms from the container
468 TH1I* hEventsTrVtx = 0;
470 if (!fEventInspector.FetchHistograms(list, hEventsTr,
471 hEventsTrVtx, hTriggers)) {
472 AliError(Form("Didn't get histograms from event selector "
473 "(hEventsTr=%p,hEventsTrVtx=%p)",
474 hEventsTr, hEventsTrVtx));
479 TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
481 AliError(Form("Couldn't get our summed histogram from output "
482 "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
487 // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
488 TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
489 TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
490 dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
491 dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
492 dNdeta->Divide(norm);
494 dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
499 MakeRingdNdeta(list, "ringSums", list, "ringResults");
500 MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24);
502 fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
503 fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
504 fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));