1 #include "AliForwardMultiplicity.h"
2 #include "AliTriggerAnalysis.h"
3 #include "AliPhysicsSelection.h"
5 #include "AliFMDAnaParameters.h"
6 #include "AliESDEvent.h"
7 #include "AliAODHandler.h"
8 #include "AliMultiplicity.h"
9 #include "AliInputEventHandler.h"
11 #include <TDirectory.h>
14 //====================================================================
15 AliForwardMultiplicity::AliForwardMultiplicity()
16 : AliAnalysisTaskSE(),
35 //____________________________________________________________________
36 AliForwardMultiplicity::AliForwardMultiplicity(const char* name)
37 : AliAnalysisTaskSE(name),
47 fSharingFilter("sharing"),
48 fDensityCalculator("density"),
49 fCorrections("corrections"),
50 fHistCollector("collector"),
54 DefineOutput(1, TList::Class());
55 // DefineOutput(2, TTree::Class());
58 //____________________________________________________________________
59 AliForwardMultiplicity::AliForwardMultiplicity(const AliForwardMultiplicity& o)
60 : AliAnalysisTaskSE(o),
61 fHEventsTr(o.fHEventsTr),
62 fHEventsTrVtx(o.fHEventsTrVtx),
63 fHTriggers(o.fHTriggers),
70 fSharingFilter(o.fSharingFilter),
71 fDensityCalculator(o.fDensityCalculator),
72 fCorrections(o.fCorrections),
73 fHistCollector(o.fHistCollector),
79 //____________________________________________________________________
80 AliForwardMultiplicity&
81 AliForwardMultiplicity::operator=(const AliForwardMultiplicity& o)
83 fHEventsTr = o.fHEventsTr;
84 fHEventsTrVtx = o.fHEventsTrVtx;
85 fHTriggers = o.fHTriggers;
87 fFirstEvent = o.fFirstEvent;
88 fSharingFilter = o.fSharingFilter;
89 fDensityCalculator = o.fDensityCalculator;
90 fCorrections = o.fCorrections;
91 fHistCollector = o.fHistCollector;
100 //____________________________________________________________________
102 AliForwardMultiplicity::Init()
107 //____________________________________________________________________
109 AliForwardMultiplicity::InitializeSubs()
111 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
114 fHEventsTr = new TH1I("nEventsTr", "Number of events w/trigger",
118 fHEventsTr->SetXTitle("v_{z} [cm]");
119 fHEventsTr->SetYTitle("# of events");
120 fHEventsTr->SetFillColor(kRed+1);
121 fHEventsTr->SetFillStyle(3001);
122 fHEventsTr->SetDirectory(0);
123 fList->Add(fHEventsTr);
124 // fHEventsTr->Sumw2();
126 fHEventsTrVtx = new TH1I("nEventsTrVtx",
127 "Number of events w/trigger and vertex",
131 fHEventsTrVtx->SetXTitle("v_{z} [cm]");
132 fHEventsTrVtx->SetYTitle("# of events");
133 fHEventsTrVtx->SetFillColor(kBlue+1);
134 fHEventsTrVtx->SetFillStyle(3001);
135 fHEventsTrVtx->SetDirectory(0);
136 fList->Add(fHEventsTrVtx);
137 // fHEventsTrVtx->Sumw2();
140 fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
141 fHTriggers->SetFillColor(kRed+1);
142 fHTriggers->SetFillStyle(3001);
143 fHTriggers->SetStats(0);
144 fHTriggers->SetDirectory(0);
145 fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
146 fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
147 fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
148 fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
149 fHTriggers->GetXaxis()->SetBinLabel(5,"A");
150 fHTriggers->GetXaxis()->SetBinLabel(6,"B");
151 fHTriggers->GetXaxis()->SetBinLabel(7,"C");
152 fHTriggers->GetXaxis()->SetBinLabel(8,"E");
153 fList->Add(fHTriggers);
155 TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
159 fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
161 fHData->SetDirectory(0);
164 fHistCollector.Init(*(fHEventsTr->GetXaxis()));
167 //____________________________________________________________________
169 AliForwardMultiplicity::UserCreateOutputObjects()
173 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
175 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
177 AliFatal("No AOD output handler set in analysis manager");
180 TObject* obj = &fAODFMD;
181 ah->AddBranch("AliAODForwardMult", &obj);
183 fSharingFilter.DefineOutput(fList);
184 fDensityCalculator.DefineOutput(fList);
185 fCorrections.DefineOutput(fList);
187 // fTree = new TTree("T", "T");
188 // fTree->Branch("forward", &fAODFMD);
190 // PostData(1, fList);
191 // PostData(2, fTree);
193 //____________________________________________________________________
195 AliForwardMultiplicity::UserExec(Option_t*)
197 // Get the input data
198 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
200 AliWarning("No ESD event found for input event");
205 static Int_t nEvents = 0;
207 if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
210 // On the first event, initialize the parameters
212 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
213 pars->SetParametersFromESD(esd);
224 // Read trigger information from the ESD and store in AOD object
226 if (!AliForwardUtil::ReadTriggers(esd, triggers, fHTriggers)) {
228 AliWarning("Failed to read triggers from ESD");
232 fAODFMD.SetTriggerBits(triggers);
234 // Mark this event for storage
237 // Check if this is a high-flux event
238 const AliMultiplicity* testmult = esd->GetMultiplicity();
241 AliWarning("No central multiplicity object found");
245 Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
247 // Get the FMD ESD data
248 AliESDFMD* esdFMD = esd->GetFMDData();
251 AliWarning("No FMD data found in ESD");
256 // Get the vertex information
258 Bool_t vzOk = AliForwardUtil::ReadVertex(esd, vz);
260 fHEventsTr->Fill(vz);
263 AliWarning("Failed to read vertex from ESD");
267 fHEventsTrVtx->Fill(vz);
269 // Get the vertex bin
270 Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
272 if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) {
274 AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
275 vz, fHEventsTr->GetXaxis()->GetXmin(),
276 fHEventsTr->GetXaxis()->GetXmax()));
281 // Apply the sharing filter (or hit merging or clustering if you like)
282 if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
284 AliWarning("Sharing filter failed!");
289 // Calculate the inclusive charged particle density
290 if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
291 AliWarning("Density calculator failed!");
295 // Do the secondary and other corrections.
296 if (!fCorrections.Correct(fHistos, ivz)) {
297 AliWarning("Corrections failed");
301 if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
302 AliWarning("Histogram collector failed");
306 if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
307 fHData->Add(&(fAODFMD.GetHistogram()));
312 //____________________________________________________________________
314 AliForwardMultiplicity::Terminate(Option_t*)
316 TList* list = dynamic_cast<TList*>(GetOutputData(1));
318 AliError(Form("No output list defined (%p)", GetOutputData(1)));
319 if (GetOutputData(1)) GetOutputData(1)->Print();
323 // Get our histograms from the container
324 TH1I* hEventsTr = static_cast<TH1I*>(list->FindObject("nEventsTr"));
325 TH1I* hEventsTrVtx = static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
326 TH2D* hData = static_cast<TH2D*>(list->FindObject("d2Ndetadphi"));
327 if (!hData || !hEventsTr || !hEventsTrVtx) {
328 AliError(Form("one or more histograms could not be found in output "
329 "list %s (hEventsTr=%p,hEventsTrVtx=%p,d2Ndetadphi=%p)",
330 list->GetName(), hEventsTr, hEventsTrVtx, hData));
335 // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
336 TH1D* dNdeta = hData->ProjectionX("dNdeta", 1, -1, "e");
337 TH1D* norm = hData->ProjectionX("dNdeta", 0, 1, "");
338 dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
339 dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
340 dNdeta->Divide(norm);
342 dNdeta->Scale(Double_t(hEventsTrVtx->GetEntries())/hEventsTr->GetEntries(),
346 fSharingFilter.ScaleHistograms(list,hEventsTr->Integral());
347 fDensityCalculator.ScaleHistograms(list,hEventsTrVtx->Integral());
348 fCorrections.ScaleHistograms(list,hEventsTrVtx->Integral());
351 //____________________________________________________________________
353 AliForwardMultiplicity::MarkEventForStore() const
355 // Make sure the AOD tree is filled
356 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
358 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
360 AliFatal("No AOD output handler set in analysis manager");
362 ah->SetFillAOD(kTRUE);
365 //____________________________________________________________________
367 AliForwardMultiplicity::Print(Option_t*) const