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);
343 fEventInspector.ProcessMC(mcEvent, triggers, ivzMC, vzMC, b, cMC,
345 fEventInspector.CompareResults(vz, vzMC, cent, b, npart, nbin);
350 Bool_t isAccepted = true;
351 if (found & AliFMDEventInspector::kNoEvent) isAccepted = false; // return;
352 if (found & AliFMDEventInspector::kNoTriggers) isAccepted = false; // return;
353 //MarkEventForStore();
354 // Always set the B trigger - each MC event _is_ a collision
355 triggers |= AliAODForwardMult::kB;
356 // Set trigger bits, and mark this event for storage
357 fAODFMD.SetTriggerBits(triggers);
358 fAODFMD.SetSNN(fEventInspector.GetEnergy());
359 fAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
360 fAODFMD.SetCentrality(cent);
361 fAODFMD.SetNClusters(nClusters);
363 fMCAODFMD.SetTriggerBits(triggers);
364 fMCAODFMD.SetSNN(fEventInspector.GetEnergy());
365 fMCAODFMD.SetSystem(fEventInspector.GetCollisionSystem());
366 fMCAODFMD.SetCentrality(cent);
367 fMCAODFMD.SetNClusters(nClusters);
369 //All events should be stored - HHD
370 //if (isAccepted) MarkEventForStore();
372 // Disable this check on SPD - will bias data
373 // if (found & AliFMDEventInspector::kNoSPD) isAccepted = false; // return;
374 if (found & AliFMDEventInspector::kNoFMD) isAccepted = false; // return;
375 if (found & AliFMDEventInspector::kNoVertex) isAccepted = false; // return;
379 fMCAODFMD.SetIpZ(vz);
381 if (found & AliFMDEventInspector::kBadVertex) isAccepted = false; // return;
383 // We we do not want to use low flux specific code, we disable it here.
384 if (!fEnableLowFlux) lowFlux = false;
389 AliESDFMD* esdFMD = esd->GetFMDData();
391 // Apply the sharing filter (or hit merging or clustering if you like)
392 if (isAccepted && !fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
393 AliWarning("Sharing filter failed!");
396 if (!fSharingFilter.FilterMC(*esdFMD, *mcEvent, vz, fMCESDFMD, fPrimary)) {
397 AliWarning("MC Sharing filter failed!");
400 if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
401 // HHD if (!isAccepted) return; // Exit on MC event w/o trigger, vertex, data
403 //MarkEventForStore();
404 fSharingFilter.CompareResults(fESDFMD, fMCESDFMD);
406 // Calculate the inclusive charged particle density
407 if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux, cent, vz)) {
408 AliWarning("Density calculator failed!");
411 if (!fDensityCalculator.CalculateMC(fMCESDFMD, fMCHistos)) {
412 AliWarning("MC Density calculator failed!");
415 fDensityCalculator.CompareResults(fHistos, fMCHistos);
417 if (fEventInspector.GetCollisionSystem() == AliFMDEventInspector::kPbPb) {
418 if (!fEventPlaneFinder.FindEventplane(esd, fAODEP, &(fAODFMD.GetHistogram()) , &fHistos))
419 AliWarning("Eventplane finder failed!");
422 // Do the secondary and other corrections.
423 if (!fCorrections.Correct(fHistos, ivz)) {
424 AliWarning("Corrections failed");
427 if (!fCorrections.CorrectMC(fMCHistos, ivz)) {
428 AliWarning("MC Corrections failed");
431 fCorrections.CompareResults(fHistos, fMCHistos);
433 if (!fHistCollector.Collect(fHistos, fRingSums,
434 ivz, fAODFMD.GetHistogram())) {
435 AliWarning("Histogram collector failed");
438 if (!fHistCollector.Collect(fMCHistos, fMCRingSums,
439 ivz, fMCAODFMD.GetHistogram())) {
440 AliWarning("MC Histogram collector failed");
444 if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
445 fHData->Add(&(fAODFMD.GetHistogram()));
450 //____________________________________________________________________
452 AliForwardMCMultiplicityTask::Terminate(Option_t*)
460 TList* list = dynamic_cast<TList*>(GetOutputData(1));
462 AliError(Form("No output list defined (%p)", GetOutputData(1)));
463 if (GetOutputData(1)) GetOutputData(1)->Print();
467 // Get our histograms from the container
469 TH1I* hEventsTrVtx = 0;
471 if (!fEventInspector.FetchHistograms(list, hEventsTr,
472 hEventsTrVtx, hTriggers)) {
473 AliError(Form("Didn't get histograms from event selector "
474 "(hEventsTr=%p,hEventsTrVtx=%p)",
475 hEventsTr, hEventsTrVtx));
480 TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
482 AliError(Form("Couldn't get our summed histogram from output "
483 "list %s (d2Ndetadphi=%p)", list->GetName(), hData));
488 // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
489 TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
490 TH1D* norm = hData->ProjectionX("norm", 0, 1, "");
491 dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
492 dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
493 dNdeta->Divide(norm);
495 dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
500 MakeRingdNdeta(list, "ringSums", list, "ringResults");
501 MakeRingdNdeta(list, "mcRingSums", list, "mcRingResults", 24);
503 fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
504 fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
505 fCorrections.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));