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("nEvents", "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 // fHEventsTr->Sumw2();
125 fHEventsTrVtx = new TH1I("nEventsTrVtx",
126 "Number of events w/trigger and vertex",
130 fHEventsTrVtx->SetXTitle("v_{z} [cm]");
131 fHEventsTrVtx->SetYTitle("# of events");
132 fHEventsTrVtx->SetFillColor(kBlue+1);
133 fHEventsTrVtx->SetFillStyle(3001);
134 fHEventsTrVtx->SetDirectory(0);
135 // fHEventsTrVtx->Sumw2();
138 fHTriggers = new TH1I("triggers", "Triggers", 10, 0, 10);
139 fHTriggers->SetFillColor(kRed+1);
140 fHTriggers->SetFillStyle(3001);
141 fHTriggers->SetStats(0);
142 fHTriggers->SetDirectory(0);
143 fHTriggers->GetXaxis()->SetBinLabel(1,"INEL");
144 fHTriggers->GetXaxis()->SetBinLabel(2,"INEL>0");
145 fHTriggers->GetXaxis()->SetBinLabel(3,"NSD");
146 fHTriggers->GetXaxis()->SetBinLabel(4,"Empty");
147 fHTriggers->GetXaxis()->SetBinLabel(5,"A");
148 fHTriggers->GetXaxis()->SetBinLabel(6,"B");
149 fHTriggers->GetXaxis()->SetBinLabel(7,"C");
150 fHTriggers->GetXaxis()->SetBinLabel(8,"E");
152 TAxis e(pars->GetNetaBins(), pars->GetEtaMin(), pars->GetEtaMax());
156 fHData = static_cast<TH2D*>(fAODFMD.GetHistogram().Clone("d2Ndetadphi"));
158 fHData->SetDirectory(0);
159 fSharingFilter.Init();
160 fHistCollector.Init(*(fHEventsTr->GetXaxis()));
163 //____________________________________________________________________
165 AliForwardMultiplicity::UserCreateOutputObjects()
169 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
171 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
173 AliFatal("No AOD output handler set in analysis manager");
176 TObject* obj = &fAODFMD;
177 ah->AddBranch("AliAODForwardMult", &obj);
179 // fTree = new TTree("T", "T");
180 // fTree->Branch("forward", &fAODFMD);
183 // PostData(2, fTree);
185 //____________________________________________________________________
187 AliForwardMultiplicity::UserExec(Option_t*)
189 // Get the input data
190 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
192 AliWarning("No ESD event found for input event");
197 static Int_t nEvents = 0;
199 if (nEvents % 100 == 0) AliInfo(Form("Event # %6d", nEvents));
202 // On the first event, initialize the parameters
204 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
205 pars->SetParametersFromESD(esd);
216 // Read trigger information from the ESD and store in AOD object
217 if (!ReadTriggers(esd)) {
219 AliWarning("Failed to read triggers from ESD");
224 // Mark this event for storage
227 // Check if this is a high-flux event
228 const AliMultiplicity* testmult = esd->GetMultiplicity();
231 AliWarning("No central multiplicity object found");
235 Bool_t lowFlux = testmult->GetNumberOfTracklets() < fLowFluxCut;
237 // Get the FMD ESD data
238 AliESDFMD* esdFMD = esd->GetFMDData();
241 AliWarning("No FMD data found in ESD");
246 // Get the vertex information
248 Bool_t vzOk = ReadVertex(esd, vz);
250 fHEventsTr->Fill(vz);
253 AliWarning("Failed to read vertex from ESD");
257 fHEventsTrVtx->Fill(vz);
259 // Get the vertex bin
260 Int_t ivz = fHEventsTr->GetXaxis()->FindBin(vz)-1;
262 if (ivz < 0 || ivz >= fHEventsTr->GetXaxis()->GetNbins()) {
264 AliWarning(Form("Vertex @ %f outside of range [%f,%f]",
265 vz, fHEventsTr->GetXaxis()->GetXmin(),
266 fHEventsTr->GetXaxis()->GetXmax()));
271 // Apply the sharing filter (or hit merging or clustering if you like)
272 if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD, vz)) {
274 AliWarning("Sharing filter failed!");
279 // Calculate the inclusive charged particle density
280 if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) {
281 AliWarning("Density calculator failed!");
285 // Do the secondary and other corrections.
286 if (!fCorrections.Correct(fHistos, ivz)) {
287 AliWarning("Corrections failed");
291 if (!fHistCollector.Collect(fHistos, ivz, fAODFMD.GetHistogram())) {
292 AliWarning("Histogram collector failed");
296 if (fAODFMD.IsTriggerBits(AliAODForwardMult::kInel))
297 fHData->Add(&(fAODFMD.GetHistogram()));
300 //____________________________________________________________________
302 AliForwardMultiplicity::Terminate(Option_t*)
304 TList* list = dynamic_cast<TList*>(GetOutputData(1));
306 AliError("No output list defined");
309 // TH1D* dNdeta = fHData->ProjectionX("dNdeta", 0, -1, "e");
310 TH1D* dNdeta = fHData->ProjectionX("dNdeta", 1, -1, "e");
311 TH1D* norm = fHData->ProjectionX("dNdeta", 0, 1, "");
312 dNdeta->SetTitle("dN_{ch}/d#eta in the forward regions");
313 dNdeta->SetYTitle("#frac{1}{N}#frac{dN_{ch}}{d#eta}");
314 dNdeta->Divide(norm);
316 dNdeta->Scale(Double_t(fHEventsTrVtx->GetEntries())/fHEventsTr->GetEntries(),
319 list->Add(fHEventsTr);
320 list->Add(fHEventsTrVtx);
321 list->Add(fHTriggers);
325 TList* last = new TList;
326 last->SetName("LastEvent");
328 last->Add(&fAODFMD.GetHistogram());
329 last->Add(fHistos.fFMD1i);
330 last->Add(fHistos.fFMD2i);
331 last->Add(fHistos.fFMD2o);
332 last->Add(fHistos.fFMD3i);
333 last->Add(fHistos.fFMD3o);
336 fSharingFilter.ScaleHistograms(fHEventsTr->Integral());
337 fSharingFilter.Output(list);
339 fDensityCalculator.ScaleHistograms(fHEventsTrVtx->Integral());
340 fDensityCalculator.Output(list);
342 fCorrections.ScaleHistograms(fHEventsTrVtx->Integral());
343 fCorrections.Output(list);
346 //____________________________________________________________________
348 AliForwardMultiplicity::MarkEventForStore() const
350 // Make sure the AOD tree is filled
351 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
353 dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
355 AliFatal("No AOD output handler set in analysis manager");
357 ah->SetFillAOD(kTRUE);
359 //____________________________________________________________________
361 AliForwardMultiplicity::ReadTriggers(AliESDEvent* esd)
363 // Get the analysis manager - should always be there
364 AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
366 AliWarning("No analysis manager defined!");
370 // Get the input handler - should always be there
371 AliInputEventHandler* ih =
372 static_cast<AliInputEventHandler*>(am->GetInputEventHandler());
374 AliWarning("No input handler");
378 // Get the physics selection - add that by using the macro
379 // AddTaskPhysicsSelection.C
380 AliPhysicsSelection* ps =
381 static_cast<AliPhysicsSelection*>(ih->GetEventSelection());
383 AliWarning("No physics selection");
387 // Check if this is a collision candidate (INEL)
388 Bool_t inel = ps->IsCollisionCandidate(esd);
390 fAODFMD.SetTriggerBits(AliAODForwardMult::kInel);
391 fHTriggers->Fill(.5);
395 // IF this is inel, see if we have a tracklet
397 const AliMultiplicity* spdmult = esd->GetMultiplicity();
399 AliWarning("No SPD multiplicity");
402 Int_t n = spdmult->GetNumberOfTracklets();
403 for (Int_t j = 0; j < n; j++) {
404 if(TMath::Abs(spdmult->GetEta(j)) < 1) {
405 fAODFMD.SetTriggerBits(AliAODForwardMult::kInelGt0);
406 fHTriggers->Fill(1.5);
413 // Analyse some trigger stuff
414 AliTriggerAnalysis ta;
415 if (ta.IsOfflineTriggerFired(esd, AliTriggerAnalysis::kNSD1)) {
416 fAODFMD.SetTriggerBits(AliAODForwardMult::kNSD);
417 fHTriggers->Fill(2.5);
421 TString triggers = esd->GetFiredTriggerClasses();
422 if (triggers.Contains("CBEAMB-ABCE-NOPF-ALL")) {
423 fAODFMD.SetTriggerBits(AliAODForwardMult::kEmpty);
424 fHTriggers->Fill(3.5);
427 if (triggers.Contains("CINT1A-ABCE-NOPF-ALL")) {
428 fAODFMD.SetTriggerBits(AliAODForwardMult::kA);
429 fHTriggers->Fill(4.5);
432 if (triggers.Contains("CINT1B-ABCE-NOPF-ALL")) {
433 fAODFMD.SetTriggerBits(AliAODForwardMult::kB);
434 fHTriggers->Fill(5.5);
438 if (triggers.Contains("CINT1C-ABCE-NOPF-ALL")) {
439 fAODFMD.SetTriggerBits(AliAODForwardMult::kC);
440 fHTriggers->Fill(6.5);
443 if (triggers.Contains("CINT1-E-NOPF-ALL")) {
444 fAODFMD.SetTriggerBits(AliAODForwardMult::kE);
445 fHTriggers->Fill(7.5);
450 //____________________________________________________________________
452 AliForwardMultiplicity::ReadVertex(AliESDEvent* esd, Double_t& vz)
455 const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
458 AliWarning("No SPD vertex found in ESD");
463 // Check that enough tracklets contributed
464 if(vertex->GetNContributors() <= 0) {
466 AliWarning(Form("Number of contributors to vertex is %d<=0",
467 vertex->GetNContributors()));
472 // Check that the uncertainty isn't too large
473 if (vertex->GetZRes() > 0.1) {
475 AliWarning(Form("Uncertaintity in Z of vertex is too large %f > 0.1",
481 // Get the z coordiante
487 //____________________________________________________________________
489 AliForwardMultiplicity::Print(Option_t*) const